STA303 - Assignment 2

STA303 - Assignment 2

This assignment is worth 5% of your final grade. It is also intended as preparation for Test 2
(worth 20%) and your final exam, so making a good effort here can help you get up to 33%
of your final grade. You will get your feedback on Assignment 2 before Test 2.
Submission is via Crowdmark NOT Quercus. You will receive an email from Crowdmark.
Contact Head TA Crystal Chen (chy.chen@mail.utoronto.ca) if you do not receive an email.
You should be able to do Question 1 by the end of week 4, Question 2 by the end of week 5
and Question 3 by the end of week 7.
• Question 1 uses data about the IQ and language test scores for students in the
netherlands, (school.csv). You will need to download this from the Assignment 2
Quercus page.
• Question 2 uses smoking data, (smoking.RData) and instructions for obtaining the
data are at the beginning of the question.
• Question 3 uses Road accident data, (pedestrians.rds) and instructions for obtaining
the data are at the beginning of the question.
Note: You can use whatever packages are useful to you, i.e., tidyverse is not required if you
prefer base R or something else. Just make sure you show which packages you are loading
in a libraries chunk. Example code for this assignment is shown with tidyverse functions in
the sta303_Assignment2_example-code.Rmd file on the Assignment 2 Quercus page.
Libraries used:
library(tidyverse)
install.packages("Pmisc", repos = "http://R-Forge.R-project.org",
type = "source")
1
Question 1: Linear mixed models
The file school.csv (available on Quercus) contains data on 760 Grade 8 students (i.e., most
are 11 years old) in 32 primary schools in the Netherlands. The data are adapted from
Snijders and Boskers’ Multilevel Analysis, 2nd Edition (Sage, 2012).
Table 1: Variables in the school.csv data set
Variable Description
school an ID number indicating which school the student attends
test the student’s score on an end-of-year language test
iq the student’s verbal IQ score
ses the socioeconomic status of the student’s family
sex the student’s sex
minority_status 1 if the student is an ethnic minority, 0 otherwise
Question of interest: Which variables are associated with Grade 8 students’
scores on an end-of-year language test?
Question 1a
Briefly describe why, without even looking at these data, you would have a concern about
one of the assumptions of linear regression.
Question 1b
Create a scatter plot to examine the relationship between verbal IQ scores and end-of-year
language scores. Include a line of best fit. Briefly describe what you see in the plot in the
context of the question of interest.
Question 1c
Create two new variables in the data set, mean_ses that is the mean of ses for each school,
and mean_iq that is mean of iq for each school.
Question 1d
Fit a linear model with test as the response and use iq, sex, ses, minority_status,
mean_ses and mean_iq as the covariates. Show the code for the model you fit and the
results of running summary() and confint() on the model you fit and briefly interpret the
results. (A complete interpretation here should discuss what the intercept means, and for
which subgroup of students it applies, as well as the location of the confidence intervals for
each covariate, i.e. below 0, includes 0 or above zero. Address the question of interest.)
2
Question 1e
Fit a linear mixed model with the same fixed effects as 1c and with a random intercept for
school.
Show the code for the model you fit and the results of running summary() and confint()
on the model you fit and briefly interpret the results.
Hint 1: Consider the estimated standard deviations in the summary to make sure you
understand the first two rows of the confint output.
Hint 2: If you want to suppress the ‘Computing profile confidence intervals …’ message you
can use message=FALSE in the chunk.
Question 1f
Briefly describe similarities and differences between the coefficients of the fixed effects in
the results from 1d and 1e and what causes the differences. You may wish to use the use
summaries of the data to help you. See the example code document.
Question 1g
Plot the random effects for the different schools. Does it seem reasonable to have included
these random effects?
Question 1h
Write a short paragraph summarising, what you have learned from this analysis. Focus
on answering the question of interest. Remember that interpreting confidence intervals is
preferred to point estimates and make sure any discussion of p-values and confidence intervals
are statistically correct. Also mention what proportion of the residual variation, after fitting
the fixed effects, the differences between schools accounts for.
Question 2: Generalised linear mixed models
Data from the 2014 American National Youth Tobacco Survey is available on http://pbro
wn.ca/teaching/303/data, where there is an R version of the 2014 dataset smoke.RData, a
pdf documentation file 2014-Codebook.pdf, and the code used to create the R version of
the data smokingData.R.
You can obtain the data with:
smokeFile = "smokeDownload.RData"
if (!file.exists(smokeFile)) {
download.file("http://pbrown.ca/teaching/303/data/smoke.RData",
smokeFile)
3
}
(load(smokeFile))
## [1] "smoke" "smokeFormats"
The smoke object is a data.frame containing the data, the smokeFormats gives some explanation of the variables. The colName and label columns of smokeFormats contain variable
names in smoke and descriptions respectively.
smokeFormats[smokeFormats[, "colName"] == "chewing_tobacco_snuff_or",
c("colName", "label")]
## colName
## 151 chewing_tobacco_snuff_or
## label## 151 RECODE: Used chewing tobacco, snuff, or dip on 1 or more days in the past 30 daysConsider the following model and set of results
# get rid of 9, 10 year olds and missing age and race
smokeSub = smoke[which(smoke$Age > 10 & !is.na(smoke$Race)),
]
smokeSub$ageC = smokeSub$Age - 16
library("glmmTMB")
smokeModelT = glmmTMB(chewing_tobacco_snuff_or ~ ageC * Sex +
RuralUrban + Race + (1 | state/school), data = smokeSub,
family = binomial(link = "logit"))
knitr::kable(summary(smokeModelT)$coef$cond, digits = 2)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.08 0.17 -17.91 0.00
ageC 0.36 0.03 11.97 0.00
SexF -2.04 0.13 -16.21 0.00
RuralUrbanRural 1.00 0.19 5.28 0.00
Raceblack -1.53 0.19 -8.17 0.00
Racehispanic -0.51 0.12 -4.29 0.00
Raceasian -1.12 0.35 -3.16 0.00
Racenative 0.03 0.29 0.10 0.92
Racepacific 1.12 0.39 2.87 0.00
ageC:SexF -0.33 0.06 -5.91 0.00
4
Table 3: Output of Pmisc::coefTable(smokeModelT)
est 2.5 % 97.5 %
ref prob
M:Urban:white 0.04 0.03 0.06
ageC
1.43 1.35 1.52
Sex
F 0.13 0.10 0.17
RuralUrban
Rural 2.72 1.88 3.95
Race
black 0.22 0.15 0.31
hispanic 0.60 0.47 0.76
asian 0.33 0.16 0.65
native 1.03 0.58 1.82
pacific 3.07 1.43 6.60
ageC:Sex
F 0.72 0.65 0.80
sd
school:state 0.75 0.59 0.95
state 0.31 0.13 0.74
The results from this code are shown in fig. 1.
Pmisc::ranefPlot(smokeModelT, grpvar = "state", level = 0.5,
maxNames = 12)
Pmisc::ranefPlot(smokeModelT, grpvar = "school:state", level = 0.5,
maxNames = 12, xlim = c(-1, 2.2))
Question 2a
Write down a statistical model corresponding to smokeModelT. Briefly explain the difference
between this model and a generalized linear model.
Question 2b
Briefly explain why this generalized linear mixed model with a logit link is more appropriate
for this dataset than a linear mixed model.
Question 2c
Write a paragraph assessing the hypothesis that state-level differences in chewing tobacco
usage amongst high school students are much larger than differences between schools within
a state. If one was interested in identifying locations with many tobacco chewers (in order to
5
−0.4 −0.2 0.0 0.2 0.4
x
NJ
CA
UT
IL
NH
ND
SD
NV
WA
WV
MI
NC
MT
ID
FL
MO
AZ
NY
IA
WI
CT
MN
AR
MA
VA
TN
AL
SC
GA
PA
OR
TX
OH
LA
OK
(a) state
−1.0 −0.5 0.0 0.5 1.0 1.5 2.0
x
mdr_01067788:UT mdr_01556466:NJ mdr_00212875:GA
mdr_00047084:AZ mdr_01481653:VA mdr_00302066:IL
mdr_04290788:ND mdr_00925515:PA
mdr_00663949:NH mdr_11103609:NV mdr_01405556:IL
mdr_04871112:NC mdr_04916310:TX
mdr_00605571:MT mdr_00906662:PA mdr_00197196:FL
mdr_01016856:TX mdr_04457091:FL mdr_00237942:IA
mdr_00968517:SD mdr_01144544:WI mdr_00015500:AL mdr_01540455:FL mdr_04032263:IL
mdr_01102144:WA mdr_01426926:CA
mdr_00844973:OK mdr_00065672:CA mdr_11561544:MI
mdr_05101001:CA mdr_04924771:FL
mdr_00139417:CA mdr_04886131:TX
mdr_01156418:WV mdr_00891689:PA
mdr_01556040:GA mdr_04949472:AZ
mdr_01396292:MO mdr_01398264:IL
mdr_00081781:CA mdr_00684802:NJ
mdr_00660882:NH mdr_05345295:UT
mdr_00432013:MA pss_00008044:IA
mdr_10910368:NY mdr_04875390:NJ
mdr_11540875:NC mdr_03018757:CA mdr_00069018:CA mdr_00794029:OH mdr_05351086:CA mdr_01548574:VA
mdr_00527917:MN mdr_00286967:IL
mdr_02231633:AL mdr_03004536:MI mdr_05274505:FL mdr_01821479:FL mdr_00319629:IL
mdr_01079901:VA mdr_00271297:IL
mdr_00445321:MA
mdr_03252927:CA mdr_00768068:NY mdr_00128339:CA mdr_05348819:NY mdr_00115930:CA mdr_11708431:AZ mdr_11070806:ID mdr_00527321:MN mdr_01145419:WI mdr_11068322:IL
mdr_05279646:MO mdr_01556105:IL
mdr_00057273:CA mdr_01067178:UT
mdr_00073203:CA mdr_10019718:TX mdr_00446375:MA mdr_00689618:NJ mdr_00195150:FL
mdr_00963646:SD mdr_01016234:TX
mdr_01115517:WA mdr_00680533:NJ
mdr_01087312:VA mdr_00668511:NJ mdr_00687311:NJ mdr_03048180:MO mdr_00066688:CA
mdr_00767947:NY mdr_03404522:MI mdr_00269153:IL mdr_00832281:OH mdr_00873118:OR mdr_01398111:GA mdr_00040347:AZ mdr_00996681:TX mdr_02042664:LA mdr_00884959:PA
mdr_05272521:FL mdr_03400473:AR mdr_11454525:AZ pss_00016465:NJ mdr_00065517:CA pss_00007641:GA mdr_00308620:IL mdr_99999999:CA pss_00006042:FL
mdr_04755594:FL pss_00006046:FL
mdr_00872906:OR mdr_11235422:CA mdr_04451774:ID mdr_00014922:AL mdr_00184735:FL mdr_01837492:AZ mdr_00998275:TX mdr_01011492:TX mdr_05350800:WI mdr_11129483:OK mdr_00308474:IL mdr_00986545:TN mdr_00173267:CT mdr_04776653:IL mdr_10010102:GA mdr_00767351:NY mdr_01051686:TX mdr_00657639:ND mdr_11717523:TX mdr_04948076:FL mdr_10009464:FL mdr_00420852:MA mdr_05348091:NV mdr_00196960:FL mdr_00015378:AL mdr_00986636:TN mdr_02044739:MA mdr_02889711:MI mdr_11457802:CA mdr_00862822:OK mdr_05280877:MO mdr_00663858:NH mdr_00139431:CA mdr_04870297:IL mdr_00674089:NJ mdr_04286452:VA 420681007335:PA mdr_01548392:TX mdr_05348431:FL mdr_00924262:PA mdr_00184711:FL mdr_00206981:GA mdr_00767624:NY mdr_00690265:NJ mdr_00872968:OR mdr_00112134:CA mdr_00508868:MI mdr_11128518:TX mdr_00125818:CA mdr_01439014:GA mdr_00406545:LA mdr_00194364:FL mdr_11454965:AL mdr_00430766:MA mdr_01398290:IL mdr_00184072:FL mdr_00605363:MT mdr_00558306:MO mdr_00023739:AR mdr_04029515:CA mdr_00308931:IL mdr_00406478:LA mdr_00746864:NY mdr_01087324:VA mdr_00917374:PA mdr_00957659:SC mdr_05092216:GA mdr_00422549:MA mdr_00291601:IL mdr_00960864:SC mdr_00224402:GA mdr_00527333:MN mdr_11079620:TX mdr_00223185:GA mdr_00013057:AL mdr_01071002:VA mdr_11684990:OH 550960001221:WI mdr_01051636:TX mdr_00263549:IL mdr_00634053:NC mdr_00013045:AL mdr_00986765:TN mdr_01070917:VA mdr_00230047:IA mdr_00183626:FL mdr_00872944:OR mdr_00850075:OK mdr_00073758:CA mdr_10910643:LA mdr_00862535:OK mdr_01007647:TX mdr_00786503:OH 040080102941:AZ mdr_00885264:PA
(b) school
Figure 1: Conditional mean and 50pct prediction interval for random effects
6
sell chewing tobacco to children, or if you prefer to implement programs to reduce tobacco
chewing), would it be important to find individual schools with high chewing rates or would
targeting those states where chewing is most common be sufficient?
Question 3: Death on the roads
The dataset below is a subset of the data from www.gov.uk/government/statistical-datasets/ras30-reported-casualties-in-road-accidents, with all of the road traffic accidents in the
UK from 1979 to 2015. The data below consist of all pedestrians involved in motor vehicle
accidents with either fatal or slight injuries (pedestrians with moderate injuries have been
removed).
dim(pedestrians)
## [1] 1159371 7
pedestrians[1:3, ]
## time age sex Casualty_Severity Light_Conditions
## 54 1979-01-01 22:40:00 26 - 35 Male Slight Darkness - lights lit
## 65 1979-01-02 10:40:00 26 - 35 Male Slight Daylight
## 79 1979-01-02 14:25:00 46 - 55 Male Slight Daylight
## Weather_Conditions y
## 54 Snowing no high winds FALSE
## 65 Raining no high winds FALSE
## 79 Raining no high winds FALSE
table(pedestrians$Casualty_Severity, pedestrians$sex)
##
## Male Female
## Slight 637919 481811
## Fatal 24429 15212
range(pedestrians$time)
## [1] "1979-01-01 01:00:00 EST" "2015-12-31 23:35:00 EST"
Notice that men are involved in accidents more than women, and the proportion of accidents
which are fatal is higher for men than for women. This might be due in part to women being
more reluctant than men to walk outdoors late at night or in poor weather, and could also
reflect men being on average more likely to engage in risky behaviour than women.
A glm adjusting for weather and light conditions is below.
theGlm = glm(y ~ sex + age + Light_Conditions + Weather_Conditions,
data = pedestrians, family = binomial(link = "logit"))
knitr::kable(summary(theGlm)$coef, digits = 3)
7
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.177 0.020 -203.929 0.000
sexFemale -0.275 0.011 -24.665 0.000
age0 - 5 0.186 0.032 5.831 0.000
age6 - 10 -0.357 0.030 -12.030 0.000
age11 - 15 -0.504 0.029 -17.668 0.000
age16 - 20 -0.338 0.027 -12.298 0.000
age21 - 25 -0.159 0.029 -5.457 0.000
age36 - 45 0.324 0.027 12.213 0.000
age46 - 55 0.660 0.026 25.030 0.000
age56 - 65 1.138 0.025 45.355 0.000
age66 - 75 1.760 0.023 75.234 0.000
ageOver 75 2.328 0.022 104.302 0.000
Light_ConditionsDarkness - lights lit 0.995 0.012 81.220 0.000
Light_ConditionsDarkness - lights unlit 1.176 0.052 22.415 0.000
Light_ConditionsDarkness - no lighting 2.765 0.021 131.303 0.000
Light_ConditionsDarkness - lighting unknown 0.259 0.068 3.788 0.000
Weather_ConditionsRaining no high winds -0.214 0.017 -12.957 0.000
Weather_ConditionsSnowing no high winds -0.751 0.092 -8.136 0.000
Weather_ConditionsFine + high winds 0.175 0.037 4.774 0.000
Weather_ConditionsRaining + high winds -0.066 0.040 -1.648 0.099
Weather_ConditionsSnowing + high winds -0.550 0.172 -3.193 0.001
Weather_ConditionsFog or mist 0.069 0.069 0.989 0.323
Here’s another GLM with interactions.
theGlmInt = glm(y ~ sex * age + Light_Conditions + Weather_Conditions,
data = pedestrians, family = binomial(link = "logit"))
knitr::kable(summary(theGlmInt)$coef, digits = 3)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.103 0.023 -179.887 0.000
sexFemale -0.545 0.044 -12.425 0.000
age0 - 5 0.021 0.039 0.544 0.587
age6 - 10 -0.460 0.035 -13.105 0.000
age11 - 15 -0.582 0.035 -16.625 0.000
age16 - 20 -0.369 0.032 -11.461 0.000
age21 - 25 -0.149 0.033 -4.501 0.000
age36 - 45 0.322 0.031 10.508 0.000
age46 - 55 0.656 0.031 21.281 0.000
age56 - 65 1.075 0.030 35.727 0.000
age66 - 75 1.622 0.029 56.315 0.000
ageOver 75 2.180 0.027 79.597 0.000
Light_ConditionsDarkness - lights lit 0.990 0.012 80.676 0.000
8
0 10 20 30 40 50 60 70
0.01 0.02 0.05 0.10
age
prob
male
female
Figure 2: Predicted probability of being a case in baseline conditions (daylight, fine no wind)
with 99% CI using theGlmInt
Estimate Std. Error z value Pr(>|z|)
Light_ConditionsDarkness - lights unlit 1.174 0.052 22.399 0.000
Light_ConditionsDarkness - no lighting 2.746 0.021 130.165 0.000
Light_ConditionsDarkness - lighting unknown 0.257 0.068 3.759 0.000
Weather_ConditionsRaining no high winds -0.211 0.017 -12.764 0.000
Weather_ConditionsSnowing no high winds -0.746 0.092 -8.075 0.000
Weather_ConditionsFine + high winds 0.176 0.037 4.803 0.000
Weather_ConditionsRaining + high winds -0.062 0.040 -1.545 0.122
Weather_ConditionsSnowing + high winds -0.548 0.172 -3.189 0.001
Weather_ConditionsFog or mist 0.065 0.069 0.943 0.346
sexFemale:age0 - 5 0.546 0.068 7.970 0.000
sexFemale:age6 - 10 0.367 0.066 5.606 0.000
sexFemale:age11 - 15 0.285 0.062 4.603 0.000
sexFemale:age16 - 20 0.150 0.062 2.408 0.016
sexFemale:age21 - 25 -0.041 0.069 -0.596 0.551
sexFemale:age36 - 45 0.029 0.062 0.475 0.635
sexFemale:age46 - 55 0.059 0.060 0.976 0.329
sexFemale:age56 - 65 0.246 0.056 4.417 0.000
sexFemale:age66 - 75 0.406 0.052 7.877 0.000
sexFemale:ageOver 75 0.411 0.049 8.348 0.000
9
Table 6: Odds ratios for theGlm and theGlmInt.
model 1 model 2
est 2.5 97.5 est 2.5 97.5
ref prob
Male:26 - 35:Daylight:Fine no 0.02 0.01 0.02 0.02 0.02 0.02
sex
Female 0.76 0.74 0.78 0.58 0.53 0.63
age
0 - 5 1.20 1.13 1.28 1.02 0.95 1.10
11 - 15 0.60 0.57 0.64 0.56 0.52 0.60
16 - 20 0.71 0.68 0.75 0.69 0.65 0.74
21 - 25 0.85 0.81 0.90 0.86 0.81 0.92
36 - 45 1.38 1.31 1.46 1.38 1.30 1.47
46 - 55 1.93 1.84 2.04 1.93 1.81 2.05
56 - 65 3.12 2.97 3.28 2.93 2.76 3.11
6 - 10 0.70 0.66 0.74 0.63 0.59 0.68
66 - 75 5.81 5.55 6.08 5.06 4.78 5.36
Over 75 10.26 9.82 10.71 8.84 8.38 9.33
Light Conditions
Darkness - lighting unknown 1.30 1.13 1.48 1.29 1.13 1.48
Darkness - lights lit 2.70 2.64 2.77 2.69 2.63 2.76
Darkness - lights unlit 3.24 2.92 3.59 3.23 2.92 3.58
Darkness - no lighting 15.89 15.24 16.56 15.58 14.95 16.24
Weather Conditions
Fine + high winds 1.19 1.11 1.28 1.19 1.11 1.28
Fog or mist 1.07 0.93 1.23 1.07 0.93 1.22
Raining + high winds 0.94 0.87 1.01 0.94 0.87 1.02
Raining no high winds 0.81 0.78 0.83 0.81 0.78 0.84
Snowing + high winds 0.58 0.41 0.81 0.58 0.41 0.81
Snowing no high winds 0.47 0.39 0.57 0.47 0.40 0.57
sex:age
Female:0 - 5 1.73 1.51 1.97
Female:11 - 15 1.33 1.18 1.50
Female:16 - 20 1.16 1.03 1.31
Female:21 - 25 0.96 0.84 1.10
Female:36 - 45 1.03 0.91 1.16
Female:46 - 55 1.06 0.94 1.19
Female:56 - 65 1.28 1.15 1.43
Female:6 - 10 1.44 1.27 1.64
Female:66 - 75 1.50 1.36 1.66
Female:Over 75 1.51 1.37 1.66
10
Question 3a
Write a short paragraph describing a case/control model (not the results) corresponding the
theGlm and theGlmInt objects. Be sure to specify the case definition and the control group,
and what the covariates are.
Question 3b
Write a short report assessing whether the UK road accident data are consistent with the
hypothesis that women tend to be, on average, safer as pedestrians than men, particularly as
teenagers and in early adulthood. Explain which of the two models fit is more appropriate
for addressing this research question.
Question 3c
It is well established that women are generally more willing to seek medical attention for
health problems than men, and it is hypothesized that men are less likely than women to
report minor injuries caused by road accidents. Write a critical assessment of whether or not
the control group is a valid one for assessing whether women are on average better at road
safety than man.
Some code
download data
pedestrainFile = Pmisc::downloadIfOld(
'http://pbrown.ca/teaching/303/data/pedestrians.rds')
pedestrians = readRDS(pedestrainFile)
pedestrians = pedestrians[!is.na(pedestrians$time), ]
pedestrians$y = pedestrians$Casualty_Severity == 'Fatal'
Code for fig. 2
newData = expand.grid(
age = levels(pedestrians$age),
sex = c('Male', 'Female'),
Light_Conditions = levels(pedestrians$Light_Conditions)[1],
Weather_Conditions = levels(pedestrians$Weather_Conditions)[1])
thePred = as.matrix(as.data.frame(
predict(theGlmInt, newData, se.fit=TRUE)[1:2])) %*% Pmisc::ciMat(0.99)
thePred = as.data.frame(thePred)
thePred$sex =newData$sex
thePred$age = as.numeric(gsub("[[:punct:]].*|[[:alpha:]]", "", newData$age))
toPlot2 = reshape2::melt(thePred, id.vars = c('age','sex'))
toPlot3 = reshape2::dcast(toPlot2, age ~ sex + variable)
11
matplot(toPlot3$age, exp(toPlot3[,-1]),
type='l', log='y', col=rep(c('black','red'), each=3),
lty=rep(c(1,2,2),2),
ylim = c(0.007, 0.11), xaxs='i',
xlab= 'age', ylab='prob')
legend('topleft', lty=1, col=c('black','red'), legend = c('male','female'), bty='n')
12