EdPsy/Psych/Stat 587
Spring 2023
C.J. Anderson
R: Computer Lab Session 3
Wednesday March 22, 2023
In this lab you will
• Fit random intercept and slopes models
• More practice using lmer and lmerTest.
• Test parameters from the marginal model.
tau and sigma
• Robust estimation using the “robust” function.
• Compute R12 and R22 using “hlmRsq” function.
• Use of the “contrast” function.
• Computer alternate BIC and compare with BIC from lmer [Link]()
• Model refinement
For this computer lab, we’ll continue to use the USA TIMSS data where science scores are
the outcome measure.
You will need to use the following R packages:
lme4
lmerTest
texreg
optimx
1. Download various function from the course web-site, which are now in a single file.
The following command will make them all available to use (i.e., you don’t have to
run each through the console window).
sources(‘‘C:.../All [Link]’’)
This loads all functions that we will use: icc, hlmRsq, roubust, contrast, and
[Link].
1
2. To jumpstart this lab, use the lab3 template that is on the course web-site. This will
read in data, create variables and fit model of the models from labs 1 and 2 (at least
the ones you’ll need).
The lab template will read in data set “[Link]”, set up variables, and fit models
from labs 1 & 2, specificaly steps (a)–(d) below.
(a) Compute school mean math scores, school mean math.
(b) Add school center math scores, school center math.
(c) Re-scale chool mean math and school center math so you don’t get annoying
messages about scale.
• lab3$xSchoolCenterMath ← lab3$school center math /
sd(lab3$school center math)
• lab3$xSchoolMeanMath ← lab3$school mean math /
sd(lab3$school mean math)
• Check data:
sd(lab3$xSchoolCenterMath)
sd(lab3$xSchoolMeanMath)
head(lab3)
(d) Set the following variables to be factors:
lab3$idschool ← [Link](lab3$idschool)
lab3$grade ← [Link](lab3$grade)
lab3$gender ← [Link](lab3$gender)
lab3$typecommunity ← [Link](lab3$typecommunity)
lab3$genshortages ← [Link](lab3$genshortages)
3. Compute R12 , R22 , and various BIC statistics for all the models from labs 1 & 2
(except ANOVA and ANCOVA)
(a) hlmRsq([Link], [Link])
(b) For hlmRsq to run correctly, the first set of variables in the formula are the
random ones.
For example:
id ← lab3$idschool #only need to do this once
hlmRsq(lab3, model.c, id)
(c) Running [Link] is pretty easy,
[Link]([Link], id)
Note that in homework 5, you need to report R12 , R22 , AIC, and various BIC statistics.
4. Fit the model below, which we will call the “base model”. Note,
2
(a) Use hlmRsq function to get R12 and R22 .
(b) Use [Link] to get various information criteria.
(c) Use the “robust” function to compute emprirical/sandwiche/robust standard
errors. The syntax is
• sand ← robust([Link], [Link],
[Link], [Link])
• round(sand$table,digits=4)
• For example,
[Link] ← robust([Link], lab3$science,
lab3$idschool, "between/within")
round([Link]$table,digits=4)
• ***The first set of variables in the formula must be the random
ones, same order***
The “base model”:
Hierarchical model :
Level 1 :
(science)ij = β0j + β1j (xSchoolCenterMath) ij + β2j (gender)ij + β3j (grade)ij
+β4j (hoursTV)ij + β5j (hourscomputergames) ij + Rij
where Rij ∼ N (0, σ 2 ) and independent.
Level 2 :
β0j = γ00 + γ01 (xSchoolMeanMath) j + γ02 (rural)j + γ03 (suburb)j
+γ04 (urban)j + U0j
β1j = γ10 + γ11 (xSchoolMeanMath) j + γ12 (rural)j + γ13 (suburb)j
+γ14 (urban)j + U1j
β2j = γ20
β3j = γ30
β4j = γ40
β5j = γ5j
where ! ! !!
U0j 0 τ02 τ01
∼N , i.i.d.
U1j 0 τ01 τ12
and independent of Rij .
3
Linear mixed model :
(science)ij = γ00 + γ10 (xSchoolCenterMath) ij + γ20 genderij + γ30 (grade)ij
+γ40 (hoursTV)ij + γ5j (hourscomputergames) ij
+γ01 (xSchoolMeanMath) j + γ02 (rural)j + γ03 (suburb)j
+γ04 (urban)j + γ11 (xSchoolMeanMath) j (xSchoolCenterMath) ij
+γ12 (rural)j (xSchoolCenterMath) ij
+γ13 (suburb)j (xSchoolCenterMath) ij
+γ14 (urban)j (xSchoolCenterMath) ij
+U0j + U1j (xSchoolCenterMath) ij + Rij
5. Consider the above model and fit a model that you can use to check whether you
need a random slope for (xSchoolCenterMath) ij .
6. You will need to report a test of the random slope for (xSchoolCenterMath), to do
this you need to fit model with and without a random slope and then. . .
[Link] ← (anova([Link],[Link]))
p1 ← [Link][2,8]
df0 ← [Link][2,7] - 1
p0 ← pchisq([Link][2,6],df0,[Link]=FALSE)
round(pvalue ← .5*(p1+p0),digits=3)
#
# if you want a nice table
#
[Link] ←
matrix(c([Link][2,6],[Link][2,7],df0,p1,p0,pvalue),nrow=1)
[Link] ← [Link]([Link]) names([Link]) ← c("LR
statistic","df1","df0","p1","p0","pvalue") round([Link],digits=3)
7. Examine the parameter estimates for type of community (i.e., location). Using a
contrast statement, test whether the difference between the parameters for urban and
suburban are statistically different.
(a) To use the contrast function,
i. I assume that you have downloaded the contrast function from the course
web-site (file named “contrast [Link]”) and ran the contrast in the console
window (cut and paste).
ii. Create an L vector of weights to use in contrasts (or other linear
combination):
• L ← matrix(0,nrow=1,ncol=14)
4
• L[1, ]← 1
• L[1, ] ← -1
• Check it
• contrast([Link], L)
This is not sufficient to consider urban and suburban having the same effect. We also
need to test locations as main effects and in the interaction effects; that is, L would
have 2 rows.
8. You should try some other contrasts, e.g., urban & suburban fixed effects same as
rural? or Others?
If you want to enter more than one row of L where the null hypothesis is
Ho : LΓ = 0, then, for example,
L ← matrix(0, nrow=3, ncol=14)
L[1, ] ← 1
L[2, ] ← 1
L[3, ] ← 1
9. For fixed effects that are not statistically different, for example, sub-urban and urban.
The variable typecommunity can be recoded and the model re-fit to the data used
the recoded variable. In this case, to re-code, the following commands will to the
trick. Also, to check, do a likelihood ratio test to see if re-coding changed the fit of
the model to the data (i.e., γ03 = γ04 = γ13 = γ14 )
#recode 4 to 3
lab3$location3 ← ifelse(lab3$typecommunity==4,3,lab3$typecommunity)
#check distribution
table(lab3$location3)
10. For this next part use your re-coded variables for typecommunity (i.e. location3)
and
(a) Consider hoursTV as a numerical variable.
(b) Consider hoursTV as a factor. To change this to cateogorical use
lab3$hrTV <- [Link](lab3$hrTv)
(c) Do contrasts if and as warranted.
(d) If warranted, re-code hoursTV.
(e) Do a likelihood ratio test to see if re-coding leads to a significant change in the
model fit to data. To do this the easy way:
i. Fit the null model
5
ii. anova([Link],[Link])
Recall
TIMSS R data Values Description
asbgday1 hourstv Time spent watching TV or
videos
1= no time
2= < 1 hr.
3= 1–2 hrs.
4= 3–4 hrs.
5= >4 hrs.
11. Using any re-coded variables you have, for this next part,
(a) Consider hourcomputergames as a numerical variable and as a categorical
variable. Test whether it should be numeric or categorical.
(b) Do contrasts if warranted.
(c) If warranted, re-code hourcomputergames .
(d) Do a likelihood ratio test to see if re-coding leads to a significant change in the
model fit to data.
Recall
TIMSS R data Values Description
asbgday2 hourscomputergames Time spent playing computer
games
1= no time
2= less than 1 hour”
3= 1–2 hours
4= 3–4 hours
5= more than 4 hours
12. Fit any other models or do any analyses that you think may be useful or needed.