Edpsy 587: Random Intercept and Slope Models using R
Carolyn J. Anderson
1/19/2023
The document contains R script used for the high school and beyond analyses given in the lecture notes.
Set up
Libraries
I am adding one package that we haven't used before, "optimx". This allows us to change the algorithm and solves as least one estimation problem that we run into with these models and data.
library(lme4)
library(lmerTest)
library(lattice)
library(texreg)
library(optimx)
library(stargazer)Read in Data
This is where the data live on my computer:
setwd('D:/Dropbox/edps587/lectures/4 randomslope')The reading in of data and merging files is the same as in the "Random Intercepts Models in R" Rmarkdown file.
In HSB data for this example are contained in two files
- HSB1data.txt that contains level 1 (tudent level) data.
- HSB2data.txt contains level 2 (school level) data.
These two files will be read in and then merged.
hsb1
hsb1 <- read.table(file="HSB1data.txt", header=TRUE)
# variable names
names(hsb1)## [1] "id" "minority" "female" "ses" "mathach"# dimensions of data frame
dim(hsb1)## [1] 7185 5# Top & bottom of the file
head(hsb1)## id minority female ses mathach
## 1 1224 0 1 -1.528 5.876
## 2 1224 0 1 -0.588 19.708
## 3 1224 0 0 -0.528 20.349
## 4 1224 0 0 -0.668 8.781
## 5 1224 0 0 -0.158 17.898
## 6 1224 0 0 0.022 4.583tail(hsb1)## id minority female ses mathach
## 7180 9586 1 1 1.612 20.967
## 7181 9586 0 1 1.512 20.402
## 7182 9586 0 1 -0.038 14.794
## 7183 9586 0 1 1.332 19.641
## 7184 9586 0 1 -0.008 16.241
## 7185 9586 0 1 0.792 22.733hsb2
hsb2 <- read.table(file="HSB2data.txt", header=TRUE)
names(hsb2)## [1] "id" "size" "sector" "pracad" "disclim" "himinty" "meanses"dim(hsb2)## [1] 160 7# Top & bottom of the file
head(hsb2)## id size sector pracad disclim himinty meanses
## 1 1224 842 0 0.35 1.597 0 -0.428
## 2 1288 1855 0 0.27 0.174 0 0.128
## 3 1296 1719 0 0.32 -0.137 1 -0.420
## 4 1308 716 1 0.96 -0.622 0 0.534
## 5 1317 455 1 0.95 -1.694 1 0.351
## 6 1358 1430 0 0.25 1.535 0 -0.014tail(hsb2)## id size sector pracad disclim himinty meanses
## 155 9347 1067 1 0.58 -0.905 0 0.218
## 156 9359 1184 1 0.69 -0.475 0 0.360
## 157 9397 1314 0 0.44 -0.231 0 0.140
## 158 9508 1119 1 0.52 -1.138 0 -0.132
## 159 9550 1532 0 0.45 0.791 0 0.059
## 160 9586 262 1 1.00 -2.416 0 0.627Merging hsb1 and hsb2
We will make a long file.
I know that the data sets are already sorted by the school id, which will be the matching variable.
hsb <- merge(hsb1,hsb2, by=c('id'))
dim(hsb)## [1] 7185 11names(hsb)## [1] "id" "minority" "female" "ses" "mathach" "size"
## [7] "sector" "pracad" "disclim" "himinty" "meanses"We will want to center SES (eventually) but let's do it now
hsb$cses <- hsb$ses - hsb$meansesSome exploratory looks at the data
# Number of schools
(nschools <- nrow(hsb2))## [1] 160Get a list of all the variables,
names(hsb)## [1] "id" "minority" "female" "ses" "mathach" "size"
## [7] "sector" "pracad" "disclim" "himinty" "meanses" "cses"We would like some basic descriptive statistics for numeric variables
x <- hsb[, c(4,5,6,8,9,10,11)]
stargazer(x,type="text")##
## ==================================================================
## Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
## ------------------------------------------------------------------
## ses 7,185 0.0001 0.779 -3.758 -0.538 0.602 2.692
## mathach 7,185 12.748 6.878 -2.832 7.275 18.317 24.993
## size 7,185 1,056.862 604.172 100 565 1,436 2,713
## pracad 7,185 0.534 0.251 0.000 0.320 0.700 1.000
## disclim 7,185 -0.132 0.944 -2.416 -0.817 0.460 2.756
## himinty 7,185 0.280 0.449 0 0 1 1
## meanses 7,185 0.006 0.414 -1.188 -0.317 0.333 0.831
## ------------------------------------------------------------------For my discrete variables
table(hsb$minority)##
## 0 1
## 5211 1974table(hsb$female)##
## 0 1
## 3390 3795table(hsb$sector)##
## 0 1
## 3642 3543Graphs of the data
Is there variability in group (school) means for SES?
plot(hsb$id, hsb$meanses, type='p', pch=21,
main='HSB: Group Mean SES by School ID',
xlab='School ID',
ylab='Group Mean SES')
How much variability within school in terms of raw ses?
plot(hsb$id, hsb$ses, type='p', pch=21,
main='HSB: Student SES by School ID',
xlab='School ID',
ylab='Group Mean SES')
What does data look like for just one school and include linear regression line?
school.1224 <- hsb[ which(hsb$id==1224), ]
plot(school.1224$ses, school.1224$mathach, type='p', pch=21,
main="HSB: Math by SES for School 1224",
xlab='Student SES',
ylab='Math Achievement')
reg.1224 <- lm(school.1224$mathach~school.1224$ses, data=school.1224)
abline(reg.1224)
We might want to look at a bunch of schools but 160 graphs is a lot, so we'll take a random sample of 20 and just look at these.
and to get the random sample
groups <- unique(hsb$id)[sample(1:160,20)]
subset <- hsb[hsb$id%in%groups,]We can now plot the data for these 20 schools.
xyplot(mathach ~ ses | id, data=subset, col.line='black',
type=c('p'),
main='Varability in Math ~ SES relationship')
However, we would like to also put in regression lines for each of these schools, which mean that we have to change what is in type,
xyplot(mathach ~ ses | id, data=subset, col.line='black',
type=c('p','r'),
main='Varability in Math ~ SES relationship')
What to graphs look like when we used school mean centered SES?
r xyplot(mathach ~ cses | id, data=subset, col.line='black', type=c('p','r'), main='Math Achievement by School Centered SES', xlab='School Centered SES', ylab='Math Achievement Test Scores')
It is also useful to look at all the regressions in a single plot; however, I need consequtive integers to index the school (i.e., index that goes from 1 - 160)
sch_id <- seq(1:nschools)
id <- unique(hsb$id)
sch_info <- as.data.frame(cbind(id,sch_id))
hsb <- merge(hsb,sch_info,by=c("id"))To make the graph, we first set up a data frame that has the specifies the title, axis labels, etc. Note that type='n' doesn't plot any data
To make the graph, we first set up a data frame that has the specifies the title, axis labels, etc. ane then add in the regression lines Note that type='n' doesn't plot any data
plot(hsb$ses,hsb$mathach,type = 'n',
ylim = c(-10, 30),
xlim = c(-3, 3),
cex.main = 1.5,
xlab = 'SES',
ylab = "Math Scores",
main = "Separate Regression for Each School "
)
for(i in 1:nschools){
sub <- hsb[which(hsb$sch_id==i),]
fitted <- fitted(lm(mathach~ses,data=sub))
lines(sub$ses,fitted,col=i)
}
The following code prodcues the graph on page 45 of the lecture note on random intercept and slope models. This is the variance function
x <- seq(from=-2, to=2, by=.1)
total<- 8 + 2*(-0.5)*x + (1.5)*(x**2) + 30
plot(x,total,type = 'l',
cex.main = 1.5,
col='blue',
xlim=c(-2,2),
xlab = 'Value of x',
ylab = "Estimated Variance of x",
main = "Variance Function: Illustration of Heteroscedasticity"
)
Create More Variables
The data file already has school centered SES (ie., cses). But we also with create an aggregate variables for minority (which will be proportion of mintority students in a given school)
meanminority <- aggregate(minority ~ id , data=hsb, mean)
names(meanminority) <- c('id', 'meanminority')
hsb <- merge(hsb,meanminority, by =c('id'))We do the same for females,
meanfemale <- aggregate(female ~ id , data=hsb, mean)
names(meanfemale) <- c('id', 'meanfemale')
hsb <- merge(hsb,meanfemale, by =c('id'))At this point you might want to consider saving teh data to a file, especially if you plan to use it again,
write.table(hsb, 'hsb.txt', row.names=F, na='.')Model Fitting
Null HLM
#
# Random inercept / empty HLM / Null model / RANDOM effects ANOVA
#
summary(model0 <- lmer(mathach ~ 1 + (1 | id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 47121.8 47142.4 -23557.9 47115.8 7182
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.06262 -0.75365 0.02676 0.76070 2.74184
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 8.553 2.925
## Residual 39.148 6.257
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6371 0.2436 157.6209 51.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Computing ICC (and function)
Cut and paste
icc <- 8.553/(8.553+39.148)
icc## [1] 0.1793044Or extract the information to compute
vars <- as.data.frame(VarCorr(model0))[4]
total <- sum(vars)
tau00 <- vars[1,1]
icc <- tau00/total
icc## [1] 0.1793109Best yet, function since we do this over and over, use the following function that I wrote:
icc <- function(modl){
vars <- as.data.frame(VarCorr(modl))[4]
total <- sum(vars)
tau00 <- vars[1,1]
icc <- tau00/total
return(icc)
}
# example of how to use this
icc(model0)## [1] 0.1793109Random intercept and slope for centered ses
By centering here, we will be looking at the impact of where a student is relavtiave to their peers in terms of SES.
summary(model1 <- lmer(mathach ~ 1 + cses + (1 | id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46728.4 46755.9 -23360.2 46720.4 7181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09692 -0.73196 0.01945 0.75739 2.91423
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 8.612 2.935
## Residual 37.005 6.083
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6494 0.2437 157.7181 51.90 <2e-16 ***
## cses 2.1912 0.1086 7023.0025 20.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cses 0.003icc(model1)## [1] 0.1887876# for easier comparison
screenreg(list(model0,model1),
single.row=TRUE,
custom.model.names=c("Null","RI ses"),
custom.header=list("Random Intercept Models"=1:2))##
## ===============================================================
## Random Intercept Models
## ------------------------------------------
## Null RI ses
## ---------------------------------------------------------------
## (Intercept) 12.64 (0.24) *** 12.65 (0.24) ***
## cses 2.19 (0.11) ***
## ---------------------------------------------------------------
## AIC 47121.81 46728.41
## BIC 47142.45 46755.93
## Log Likelihood -23557.91 -23360.21
## Num. obs. 7185 7185
## Num. groups: id 160 160
## Var: id (Intercept) 8.55 8.61
## Var: Residual 39.15 37.01
## ===============================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Random intercept and slope & meanses for ses
To add a random slope for a level one predictor varaible, say \(x_{1ij}\), we add it to the random part of the model specification; that is, "(1 + \(x_{1ij}\) | id)". In this case, we still have the random intercept but have added a random slope for \(x_{1ij}\).
Below we fit a model with both a random intercept and random slope for cses.
summary(model2 <- lmer(mathach ~ 1 + cses + (1 + cses| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + (1 + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46723.0 46764.3 -23355.5 46711.0 7179
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09688 -0.73200 0.01796 0.75445 2.89906
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 8.6213 2.9362
## cses 0.6783 0.8236 0.02
## Residual 36.7000 6.0581
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6494 0.2437 157.7206 51.90 <2e-16 ***
## cses 2.1931 0.1278 156.1499 17.15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cses 0.012as.data.frame(VarCorr(model2))[4]## vcov
## 1 8.62127569
## 2 0.67826037
## 3 0.05031215
## 4 36.70001845deviance(model2)## [1] 46710.98# for easier comparison with what we have already fit
screenreg(list(model0,model1,model2),
single.row=TRUE,
custom.model.names=c("Null","cses", "ces"),
custom.header=list("Random Intercept"=1:2, "Random Slope"=3))##
## ==========================================================================================
## Random Intercept Random Slope
## ------------------------------------------ --------------------
## Null cses ces
## ------------------------------------------------------------------------------------------
## (Intercept) 12.64 (0.24) *** 12.65 (0.24) *** 12.65 (0.24) ***
## cses 2.19 (0.11) *** 2.19 (0.13) ***
## ------------------------------------------------------------------------------------------
## AIC 47121.81 46728.41 46722.98
## BIC 47142.45 46755.93 46764.26
## Log Likelihood -23557.91 -23360.21 -23355.49
## Num. obs. 7185 7185 7185
## Num. groups: id 160 160 160
## Var: id (Intercept) 8.55 8.61 8.62
## Var: Residual 39.15 37.01 36.70
## Var: id cses 0.68
## Cov: id (Intercept) cses 0.05
## ==========================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Graphs around page 20 of Lecture Notes
par(mfrow=c(1,1))
fixed <- fixef(model2)
mu2 <- fixed[1] + fixed[2]*hsb$cses
plot(hsb$cses, mu2, type='l',
col='black',
xlim=c(-2, 2),
ylim=c(0, 25),
xlab='School Mean Centered SES',
ylab='Fitted Math Scores',
main='HSB: Population Average Regressions (Model 2)'
)
For the graph on or about page 21 of lecture notes, sample of fitted school specific regression lines
random <- ranef(model2, level=0:2)## Warning in ranef.merMod(model2, level = 0:2): additional arguments to
## ranef.merMod ignored: levelU0j <- random$id[,1]
U1j <- random$id[,2]
plot(hsb$mathach ~ hsb$cses, type = 'n',
cex.main = 1.5,
xlab = 'School Mean Centered SES',
ylab = "Fitted Math Scores",
main = "HSB: Fitted School Regression Lines (private,white)"
)
# draw in scho0l specific lines: Private, white
for (i in 1:40){
abline(a=(fixed[1]+ U0j[i]), b=(fixed[2]+U1j[i]), col="blue")
}
The graph on or around page 22 of Side by Side Comparison
par(mfrow=c(1,2))
plot(hsb$cses, mu2, type='l',
col='black',
xlim=c(-2, 2),
ylim=c(0, 25),
xlab='School Mean Centered SES',
ylab='Fitted Math Scores',
main='HSB: Population Average Regressions (Model 2)'
)
plot(hsb$mathach ~ hsb$cses, type = 'n',
cex.main = 1.5,
xlab = 'School Mean Centered SES',
ylab = "Fitted Math Scores",
main = "HSB: Fitted School Regression Lines"
)
# draw in scho0l specific lines: Private, white
for (i in 1:40){
abline(a=(fixed[1]+ U0j[i]), b=(fixed[2]+U1j[i]), col='blue')
}
The graph on or about Page 23 of estimated variance functions for 2 models
par(mfrow=c(1,1))
ses2 <- seq(from=-2, to= 3, by=.1)
vars.ri <- as.data.frame(VarCorr(model1))[4]
total.ri <- rep(sum(vars.ri), 61)
vars.ris <- as.data.frame(VarCorr(model2))[4]
total.ris<- vars.ris[1,] + 2*vars.ris[3,]*ses2 + vars.ris[2,]*(ses2**2) + vars.ris[4,]
plot(ses2,total.ris,type = 'l',
cex.main = 1.5,
col='blue',
xlim=c(-2,3),
ylim=c(45,52),
xlab = 'School Mean Centered SES',
ylab = "Estimated Variance of Math Scores",
main = "HSB: Variance Function"
)
# random inercept
abline(a=sum(vars.ri), b=0, col='red')
legend(-2,52,
c('Random Intercept & Slope','Random Intercept'),
lty=c(1,1),
lwd=c(2.5,2.5),
col=c('blue','red'))
Discrete Level 1 predictors (fixed effects)
Note female and minority are dummy (i.e., I have not declared them to be factors in this model).
model.d1 <- lmer(mathach ~ 1 + cses + female + minority+ (1 + cses| id), data=hsb, REML=FALSE)
summary(model.d1)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + (1 + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46467.0 46522.0 -23225.5 46451.0 7177
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1821 -0.7199 0.0323 0.7602 2.8373
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 6.2629 2.5026
## cses 0.4801 0.6929 -0.17
## Residual 35.6632 5.9719
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.1422 0.2351 232.8885 60.163 < 2e-16 ***
## cses 1.8946 0.1224 157.1366 15.472 < 2e-16 ***
## female -1.2199 0.1643 7057.9081 -7.424 1.27e-13 ***
## minority -3.1188 0.2103 5764.5187 -14.834 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female
## cses -0.114
## female -0.367 0.045
## minority -0.250 0.139 0.012as.data.frame(VarCorr(model.d1))[4]## vcov
## 1 6.2628806
## 2 0.4800510
## 3 -0.3001918
## 4 35.6631590n2lnlike <- -2* logLik(model.d1)
n2lnlike # this is deviance## 'log Lik.' 46450.95 (df=8)Since we have fit a number of models, let's see where we are at.
screenreg(list(model1,model2,model.d1),
custom.model.names=c("RI","RI and RS","Discrete fixed"))##
## ======================================================================
## RI RI and RS Discrete fixed
## ----------------------------------------------------------------------
## (Intercept) 12.65 *** 12.65 *** 14.14 ***
## (0.24) (0.24) (0.24)
## cses 2.19 *** 2.19 *** 1.89 ***
## (0.11) (0.13) (0.12)
## female -1.22 ***
## (0.16)
## minority -3.12 ***
## (0.21)
## ----------------------------------------------------------------------
## AIC 46728.41 46722.98 46466.95
## BIC 46755.93 46764.26 46521.99
## Log Likelihood -23360.21 -23355.49 -23225.48
## Num. obs. 7185 7185 7185
## Num. groups: id 160 160 160
## Var: id (Intercept) 8.61 8.62 6.26
## Var: Residual 37.01 36.70 35.66
## Var: id cses 0.68 0.48
## Cov: id (Intercept) cses 0.05 -0.30
## ======================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Plot of population average model
fixed <- fixef(model.d1)
fem_mi <- fixed[1] + fixed[2]*hsb$cses + fixed[3] + fixed[4]
plot(hsb$cses, fem_mi, type='l',
col='black',
ylim = c(0,25),
xlim = c(-2,2),
xlab='School Centered SES',
ylab='Fitted Math Scores',
main='Population Average Regressions (discrete fixed)'
)
abline(a=(fixed[1] + fixed[3]), b=fixed[2], col='red')
abline(a=(fixed[1] + fixed[4]), b=fixed[2], col='blue')
abline(a=(fixed[1]), b=fixed[2], col='green')
legend(-2,25,
c('Female, Minority', 'Female, non-Minority', 'Male, Minority', 'Male, non-Minority'),
lty=c(1,1,1,1),
col=c('black', 'red', 'blue', 'green')
)
Plot of variance function
ses2 <- seq(from=-2, to= 3, by=.1)
vars.d1 <- as.data.frame(VarCorr(model.d1))[4]
total.d1<- vars.d1[1,] + 2*vars.d1[3,]*ses2 + vars.d1[2,]*(ses2**2) + vars.d1[4,]
plot(ses2,total.d1,type = 'l',
cex.main = 1.5,
col='blue',
xlim=c(-2,3),
xlab = 'School Mean Centered SES',
ylab = "Estimated Variance of Math Scores",
main = "HSB: Variance Function (discrete fixed)"
)
Alternative Treatment of categorical predictors
Showing effect of random categorical predictors when they are declared to be factors (rather than putting in say dummy codes)
hsb$gender <- as.factor(hsb$female)
hsb$race <- as.factor(hsb$minority)
hsb$boy <- as.factor(-1*(hsb$female-1))
hsb$male <- -1*(hsb$female-1)
summary(model.00 <- lmer(mathach ~ 1 + cses + gender + race + (1 + cses| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + gender + race + (1 + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46467.0 46522.0 -23225.5 46451.0 7177
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1821 -0.7199 0.0323 0.7602 2.8373
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 6.2629 2.5026
## cses 0.4801 0.6929 -0.17
## Residual 35.6632 5.9719
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.1422 0.2351 232.8885 60.163 < 2e-16 ***
## cses 1.8946 0.1224 157.1366 15.472 < 2e-16 ***
## gender1 -1.2199 0.1643 7057.9081 -7.424 1.27e-13 ***
## race1 -3.1188 0.2103 5764.5187 -14.834 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses gendr1
## cses -0.114
## gender1 -0.367 0.045
## race1 -0.250 0.139 0.012# To get variance
as.data.frame(VarCorr(model.00))[4]## vcov
## 1 6.2628806
## 2 0.4800510
## 3 -0.3001918
## 4 35.6631590This should yield that same fit statistics, interpretation of effects, significance of effects, etc.
Random Descrete Predictors
There are a number of ways to do this that lead to essentially the same result. Below we look at some.
summary(model.b1 <- lmer(mathach ~ 1 + cses + female + race + (1 + cses + female| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + race + (1 + cses + female | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46468.5 46544.2 -23223.3 46446.5 7174
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1767 -0.7197 0.0370 0.7562 2.8629
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 6.9181 2.6302
## cses 0.4569 0.6760 -0.12
## female 0.7858 0.8864 -0.41 -0.29
## Residual 35.5342 5.9611
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.1395 0.2453 170.9865 57.642 < 2e-16 ***
## cses 1.8977 0.1217 155.4180 15.589 < 2e-16 ***
## female -1.2454 0.1821 118.5827 -6.839 3.7e-10 ***
## race1 -3.1161 0.2098 5608.1500 -14.850 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female
## cses -0.089
## female -0.469 -0.011
## race1 -0.240 0.139 0.014as.data.frame(VarCorr(model.b1))[4]## vcov
## 1 6.9181349
## 2 0.4569283
## 3 0.7857687
## 4 -0.2066691
## 5 -0.9444765
## 6 -0.1727498
## 7 35.5342275Alternative Model b1
summary(model.ab1 <- lmer(mathach ~ 1 + cses + gender + race + (1 + cses + gender| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + gender + race + (1 + cses + gender | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46468.5 46544.2 -23223.3 46446.5 7174
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1767 -0.7197 0.0370 0.7562 2.8629
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 6.9181 2.6302
## cses 0.4569 0.6760 -0.12
## gender1 0.7858 0.8864 -0.41 -0.29
## Residual 35.5342 5.9611
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.1395 0.2453 170.9865 57.642 < 2e-16 ***
## cses 1.8977 0.1217 155.4180 15.589 < 2e-16 ***
## gender1 -1.2454 0.1821 118.5827 -6.839 3.7e-10 ***
## race1 -3.1161 0.2098 5608.1500 -14.850 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses gendr1
## cses -0.089
## gender1 -0.469 -0.011
## race1 -0.240 0.139 0.014as.data.frame(VarCorr(model.ab1))[4]## vcov
## 1 6.9181349
## 2 0.4569283
## 3 0.7857687
## 4 -0.2066691
## 5 -0.9444765
## 6 -0.1727498
## 7 35.5342275Another parameterization
summary(model.b2 <- lmer(mathach ~ 1 + cses + hsb$male + race + (1 + cses + hsb$male| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + hsb$male + race + (1 + cses + hsb$male | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46468.5 46544.2 -23223.3 46446.5 7174
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1767 -0.7197 0.0370 0.7562 2.8629
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 5.8148 2.4114
## cses 0.4570 0.6760 -0.23
## hsb$male 0.7858 0.8865 0.07 0.29
## Residual 35.5342 5.9611
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.8941 0.2268 160.0145 56.847 < 2e-16 ***
## cses 1.8977 0.1217 155.4232 15.588 < 2e-16 ***
## hsb$male 1.2454 0.1821 118.5627 6.839 3.7e-10 ***
## race1 -3.1161 0.2098 5608.1323 -14.850 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses hsb$ml
## cses -0.105
## hsb$male -0.296 0.011
## race1 -0.248 0.139 -0.014as.data.frame(VarCorr(model.b2))[4]## vcov
## 1 5.8148277
## 2 0.4570253
## 3 0.7858424
## 4 -0.3792555
## 5 0.1588552
## 6 0.1727952
## 7 35.5341774Bad Parameteriation
Note that this one does not work --- does not converge. Why?
# Model c.
summary(model.c <- lmer(mathach ~ 1 + cses + hsb$male + race + (0 + cses + male + female + minority| id), data=hsb, REML=FALSE))## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -1.3e+03## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + hsb$male + race + (0 + cses + male + female +
## minority | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46466.6 46569.8 -23218.3 46436.6 7170
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2390 -0.7144 0.0428 0.7584 2.8593
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id cses 0.000 0.000
## male 6.100 2.470 NaN
## female 4.921 2.218 NaN 0.92
## minority 2.093 1.447 NaN 0.25 0.36
## Residual 35.498 5.958
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.9020 0.2145 140.5681 60.135 < 2e-16 ***
## cses 1.8699 0.1086 7004.6639 17.225 < 2e-16 ***
## hsb$male 1.2824 0.1844 119.4346 6.954 2.03e-10 ***
## race1 -3.1701 0.2526 119.1022 -12.552 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses hsb$ml
## cses -0.018
## hsb$male -0.313 -0.045
## race1 -0.079 0.134 -0.047
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularas.data.frame(VarCorr(model.b2))[4]## vcov
## 1 5.8148277
## 2 0.4570253
## 3 0.7858424
## 4 -0.3792555
## 5 0.1588552
## 6 0.1727952
## 7 35.5341774This Works.
Why?
summary(model.all <- lmer(mathach ~ 1 + cses + female + minority + (1 + cses + female + minority | id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + (1 + cses + female +
## minority | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46459.8 46563.0 -23214.9 46429.8 7170
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2311 -0.7142 0.0401 0.7557 2.8849
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 6.1723 2.484
## cses 0.3856 0.621 0.03
## female 0.8501 0.922 -0.43 -0.23
## minority 1.9127 1.383 0.29 -0.66 0.12
## Residual 35.3086 5.942
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.1501 0.2358 149.3294 60.016 < 2e-16 ***
## cses 1.8726 0.1198 151.3902 15.635 < 2e-16 ***
## female -1.2581 0.1832 118.4607 -6.868 3.21e-10 ***
## minority -3.2011 0.2503 116.9638 -12.790 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female
## cses -0.037
## female -0.492 0.002
## minority -0.104 -0.007 0.032And this works.
summary(model.c <- lmer(mathach ~ 1 + cses + female + minority + (0 + cses + female + race| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + (0 + cses + female +
## race | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46459.8 46563.0 -23214.9 46429.8 7170
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2311 -0.7141 0.0401 0.7557 2.8849
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id cses 0.3857 0.6210
## female 0.8503 0.9221 -0.23
## race0 6.1724 2.4844 0.03 -0.43
## race1 10.0457 3.1695 -0.27 -0.29 0.91
## Residual 35.3085 5.9421
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.1501 0.2358 149.3246 60.015 < 2e-16 ***
## cses 1.8726 0.1198 151.4537 15.635 < 2e-16 ***
## female -1.2581 0.1832 118.4660 -6.868 3.21e-10 ***
## minority -3.2011 0.2503 116.9790 -12.789 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female
## cses -0.037
## female -0.492 0.002
## minority -0.103 -0.007 0.032You're on your own for the figures --- cut, paste & edit code for similar figures
Cross-level interactions
When a level 2 variable, say \(z_{j}\) is a predictor of a random effect (or fixed level 1 variable), we add the variable to the fixed/structural part of the formula. If \(z_{j}\) is a predictor of the slope for \(x_{ij}\), the formula would be \[y_{ij} \leftarrow 1 + x_{ij} + z_{j} + z_{j}*x_{ij} + (1 + x_{ij} | \mbox{id}) \] In other words, we add an interaction between \(z_{j}\) and \(x_{ij}\).
In the example below, we add meanses as a predictor of the effect of cses:
summary(model.3 <- lmer(mathach ~ 1 + cses + meanses + cses*meanses + (1 + cses| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + meanses + cses * meanses + (1 + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46568.4 46623.4 -23276.2 46552.4 7177
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1913 -0.7269 0.0157 0.7557 2.9203
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.6443 1.626
## cses 0.6496 0.806 -0.20
## Residual 36.7156 6.059
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6589 0.1482 155.6452 85.414 <2e-16 ***
## cses 2.1960 0.1272 156.2691 17.264 <2e-16 ***
## meanses 5.8700 0.3589 155.3156 16.357 <2e-16 ***
## cses:meanses 0.2863 0.3169 171.1168 0.903 0.368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses meanss
## cses -0.083
## meanses -0.004 0.000
## cses:meanss 0.000 0.042 -0.080as.data.frame(VarCorr(model.3))[4]## vcov
## 1 2.6443444
## 2 0.6496455
## 3 -0.2574434
## 4 36.7155992Models Refinements
Drop non-signicant Interaction
summary(model.4 <- lmer(mathach ~ 1 + cses + meanses + (1 + cses| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + meanses + (1 + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46567.2 46615.3 -23276.6 46553.2 7178
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1778 -0.7260 0.0158 0.7565 2.9403
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.6446 1.6262
## cses 0.6701 0.8186 -0.19
## Residual 36.7132 6.0591
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6588 0.1482 155.6389 85.41 <2e-16 ***
## cses 2.1912 0.1276 156.3840 17.17 <2e-16 ***
## meanses 5.8959 0.3577 154.8624 16.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses
## cses -0.083
## meanses -0.004 0.004p127 Big model
R starts to get slow, can switch to microsoft's OpenR (free), which is about 4 times faster. The speed up is due to the use of a faster BLAS (basic linear algebra system), in particular an Intel MKL libraries. Other alternative BLAS's are OpenBIAS and ATLAS.
This model fails
summary(model.5 <- lmer(mathach ~ 1 + cses + female + minority + meanses + size + sector
+ pracad + disclim + himinty +(1 + cses + female + minority| id), data=hsb, REML=FALSE))## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + size + sector +
## pracad + disclim + himinty + (1 + cses + female + minority | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46283.3 46427.8 -23120.7 46241.3 7164
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2820 -0.7157 0.0335 0.7531 2.9369
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 1.9538 1.3978
## cses 0.4039 0.6355 0.39
## female 0.7469 0.8642 -0.74 -0.14
## minority 1.6017 1.2656 -0.18 -0.84 -0.03
## Residual 35.3375 5.9445
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.148e+01 5.008e-01 1.743e+02 22.931 < 2e-16 ***
## cses 1.887e+00 1.203e-01 1.519e+02 15.688 < 2e-16 ***
## female -1.281e+00 1.737e-01 1.385e+02 -7.376 1.36e-11 ***
## minority -3.075e+00 2.429e-01 1.595e+02 -12.659 < 2e-16 ***
## meanses 3.168e+00 4.236e-01 1.428e+02 7.479 6.96e-12 ***
## size 6.937e-04 2.093e-04 1.511e+02 3.314 0.001150 **
## sector 8.130e-01 3.776e-01 1.410e+02 2.153 0.033003 *
## pracad 2.722e+00 7.952e-01 1.532e+02 3.423 0.000796 ***
## disclim -4.086e-01 1.753e-01 1.517e+02 -2.331 0.021073 *
## himinty 1.856e-01 3.142e-01 1.596e+02 0.591 0.555590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty meanss size sector pracad disclm
## cses 0.027
## female -0.287 0.018
## minority -0.020 -0.032 0.013
## meanses 0.490 0.005 0.002 0.089
## size -0.600 -0.005 0.025 -0.036 -0.181
## sector -0.147 -0.004 0.023 -0.044 0.025 0.251
## pracad -0.728 -0.001 0.060 -0.019 -0.609 0.112 -0.338
## disclim -0.306 0.000 0.093 -0.029 -0.022 -0.044 0.505 0.205
## himinty 0.157 -0.051 -0.027 -0.308 0.400 -0.178 -0.106 -0.198 -0.048
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingWe can do a bit better by rescaling school size which has a very different scale
hsb$xsize <- scale(hsb$size,center=FALSE,scale=TRUE)
summary(model.5 <- lmer(mathach ~ 1 + cses + female + minority + meanses + xsize + sector
+ pracad + disclim + himinty +(1 + cses + female + minority| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + xsize + sector +
## pracad + disclim + himinty + (1 + cses + female + minority | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46283.3 46427.8 -23120.7 46241.3 7164
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2820 -0.7157 0.0335 0.7531 2.9369
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 1.9542 1.3979
## cses 0.4039 0.6356 0.39
## female 0.7469 0.8643 -0.74 -0.14
## minority 1.6015 1.2655 -0.18 -0.84 -0.03
## Residual 35.3374 5.9445
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.4831 0.5008 174.3097 22.930 < 2e-16 ***
## cses 1.8867 0.1203 151.7666 15.687 < 2e-16 ***
## female -1.2814 0.1737 138.5124 -7.376 1.36e-11 ***
## minority -3.0752 0.2429 159.4606 -12.660 < 2e-16 ***
## meanses 3.1679 0.4236 142.7889 7.479 6.96e-12 ***
## xsize 0.8444 0.2548 151.0926 3.314 0.001151 **
## sector 0.8131 0.3776 141.0200 2.153 0.032998 *
## pracad 2.7217 0.7953 153.1523 3.422 0.000796 ***
## disclim -0.4086 0.1753 151.7272 -2.331 0.021069 *
## himinty 0.1856 0.3142 159.5661 0.591 0.555531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty meanss xsize sector pracad disclm
## cses 0.027
## female -0.287 0.018
## minority -0.020 -0.032 0.013
## meanses 0.490 0.005 0.002 0.089
## xsize -0.600 -0.005 0.025 -0.036 -0.181
## sector -0.147 -0.004 0.023 -0.044 0.025 0.251
## pracad -0.728 -0.001 0.060 -0.019 -0.609 0.112 -0.338
## disclim -0.306 0.000 0.093 -0.029 -0.022 -0.044 0.505 0.205
## himinty 0.157 -0.051 -0.027 -0.308 0.400 -0.178 -0.106 -0.198 -0.048as.data.frame(VarCorr(model.5))[4]## vcov
## 1 1.95419176
## 2 0.40392581
## 3 0.74693929
## 4 1.60152870
## 5 0.34521147
## 6 -0.88947968
## 7 -0.31798090
## 8 -0.07933781
## 9 -0.67179567
## 10 -0.02783238
## 11 35.33742752Other tricks to get models to converge
ALthough this is a topic for another lecture, but I did need an additional one to get the complex model above to converge.
In this model (see aroun) p127 of lecture notes), we do make some changes to the model by adding sector as cross-level (i.e. predictor at level 2) and dropping random gender.
(This does lead to some differences in estimates of tau relative to SAS.)
The thing that really help is changing the optimizer. This is done below using "control"
hsb$xhiminty <- scale(hsb$himinty,center=FALSE,scale=TRUE)
summary(model.6 <- lmer(mathach ~ 1 + cses + female + minority + meanses + xsize + sector + pracad + disclim + xhiminty + cses*sector + minority*sector + (1 + cses + minority| id), data=hsb, REML=FALSE,
control = lmerControl(optimizer ="Nelder_Mead") ))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + xsize + sector +
## pracad + disclim + xhiminty + cses * sector + minority *
## sector + (1 + cses + minority | id)
## Data: hsb
## Control: lmerControl(optimizer = "Nelder_Mead")
##
## AIC BIC logLik deviance df.resid
## 46250.7 46381.4 -23106.3 46212.7 7166
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2006 -0.7180 0.0336 0.7625 3.0629
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 1.2258 1.1072
## cses 0.1646 0.4057 0.34
## minority 0.7774 0.8817 -0.07 -0.38
## Residual 35.4601 5.9548
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.22677 0.50053 173.56520 22.430 < 2e-16 ***
## cses 2.31523 0.15283 153.47267 15.149 < 2e-16 ***
## female -1.24580 0.15764 5062.44507 -7.903 3.33e-15 ***
## minority -3.87559 0.31485 189.02045 -12.309 < 2e-16 ***
## meanses 2.90549 0.42846 148.36827 6.781 2.64e-10 ***
## xsize 1.07753 0.26341 163.64223 4.091 6.74e-05 ***
## sector 0.44953 0.39298 149.25275 1.144 0.2545
## pracad 3.21134 0.80350 156.84560 3.997 9.86e-05 ***
## disclim -0.36220 0.18026 165.44571 -2.009 0.0461 *
## xhiminty 0.07191 0.16945 174.64187 0.424 0.6718
## cses:sector -1.01031 0.22774 155.04826 -4.436 1.73e-05 ***
## minority:sector 1.78189 0.42240 124.78056 4.218 4.69e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty meanss xsize sector pracad disclm
## cses 0.013
## female -0.228 0.050
## minority 0.000 0.152 0.013
## meanses 0.481 0.020 0.004 0.118
## xsize -0.602 -0.024 0.023 -0.126 -0.179
## sector -0.138 -0.010 0.022 0.102 0.034 0.216
## pracad -0.732 -0.003 0.059 -0.036 -0.594 0.107 -0.348
## disclim -0.310 -0.003 0.090 -0.027 -0.025 -0.042 0.481 0.211
## xhiminty 0.145 -0.044 -0.030 -0.202 0.403 -0.184 -0.115 -0.176 -0.054
## cses:sector -0.013 -0.669 -0.021 -0.097 -0.016 0.023 0.024 0.003 0.002
## minrty:sctr -0.012 -0.102 0.002 -0.690 -0.078 0.149 -0.194 0.028 0.009
## xhmnty css:sc
## cses
## female
## minority
## meanses
## xsize
## sector
## pracad
## disclim
## xhiminty
## cses:sector 0.010
## minrty:sctr -0.031 0.085as.data.frame(VarCorr(model.6))[4]## vcov
## 1 1.22579595
## 2 0.16458404
## 3 0.77736156
## 4 0.15183165
## 5 -0.06946504
## 6 -0.13709806
## 7 35.46009963Round p127 A bit simpler model --- only random intercept
summary(model.7 <- lmer(mathach ~ 1 + cses + female + minority + meanses + xsize + sector
+ pracad + disclim + cses*sector + minority*sector
+ (1| id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + xsize + sector +
## pracad + disclim + cses * sector + minority * sector + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46242.8 46332.2 -23108.4 46216.8 7172
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1971 -0.7183 0.0326 0.7700 3.0843
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.30 1.140
## Residual 35.63 5.969
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.2320 0.4969 173.9880 22.606 < 2e-16 ***
## cses 2.3167 0.1457 7123.7634 15.902 < 2e-16 ***
## female -1.2371 0.1575 4973.0705 -7.854 4.90e-15 ***
## minority -3.8231 0.2863 3176.6041 -13.351 < 2e-16 ***
## meanses 2.8394 0.3856 167.7993 7.363 7.68e-12 ***
## xsize 1.1209 0.2593 166.6731 4.323 2.64e-05 ***
## sector 0.4775 0.3923 163.4302 1.217 0.2253
## pracad 3.1853 0.7880 157.9038 4.042 8.26e-05 ***
## disclim -0.3749 0.1793 165.8044 -2.091 0.0381 *
## cses:sector -1.0075 0.2175 7085.7694 -4.633 3.67e-06 ***
## minority:sector 1.7673 0.3862 2359.7732 4.576 4.98e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty meanss xsize sector pracad disclm
## cses -0.002
## female -0.224 0.051
## minority 0.033 0.210 0.007
## meanses 0.459 0.052 0.020 0.240
## xsize -0.592 -0.037 0.018 -0.180 -0.103
## sector -0.140 0.018 0.022 0.089 0.084 0.214
## pracad -0.720 -0.014 0.051 -0.080 -0.576 0.071 -0.370
## disclim -0.318 -0.003 0.089 -0.036 -0.012 -0.040 0.472 0.216
## cses:sector 0.001 -0.669 -0.022 -0.139 -0.025 0.027 -0.023 0.007 0.001
## minrty:sctr -0.008 -0.151 0.000 -0.719 -0.078 0.157 -0.214 0.025 0.006
## css:sc
## cses
## female
## minority
## meanses
## xsize
## sector
## pracad
## disclim
## cses:sector
## minrty:sctr 0.144as.data.frame(VarCorr(model.7))[4]## vcov
## 1 1.300117
## 2 35.631569