2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R

Edpsy 587: Random Intercept and Slope Models using R

Carolyn J. Anderson

1/19/2023

Set up
Libraries
Read in Data
hsb1
hsb2
Merging hsb1 and hsb2
Some exploratory looks at the data
Graphs of the data
Create More Variables
Model Fitting
Null HLM
Computing ICC (and function)
Random intercept and slope for centered ses
Random intercept and slope & meanses for ses
Graphs around page 20 of Lecture Notes
Discrete Level 1 predictors (fixed effects)
Plot of population average model
Plot of variance function
Alternative Treatment of categorical predictors
Random Descrete Predictors
Alternative Model b1
Another parameterization
Bad Parameteriation
This Works.
And this works.
Cross-level interactions
Models Refinements
Drop non-signicant Interaction
p127 Big model
Other tricks to get models to converge

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

1. [Link] that contains level 1 (tudent level) data.

2. [Link] contains level 2 (school level) data.

These two files will be read in and then merged.

hsb1
hsb1 <- [Link](file="[Link]", header=TRUE)

# variable names
names(hsb1)

[Link] 1/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## [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.583

tail(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.733

hsb2
hsb2 <- [Link](file="[Link]", 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.014

tail(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.627

Merging 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 11

names(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$meanses

Some exploratory looks at the data

# Number of schools
(nschools <- nrow(hsb2))

## [1] 160

Get a list of all the variables,

names(hsb)

[Link] 2/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## [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 1974

table(hsb$female)

##
## 0 1
## 3390 3795

table(hsb$sector)

##
## 0 1
## 3642 3543

Graphs 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')

[Link] 3/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R

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, [Link]='black',

type=c('p'),
main='Varability in Math ~ SES relationship')

[Link] 4/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R

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, [Link]='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, [Link]='black', type=c('p','r'), main='Math Achievement by School Centered SES',

[Link] 5/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R

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 <- [Link](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),
[Link] = 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

[Link] 6/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
x <- seq(from=-2, to=2, by=.1)

total<- 8 + 2*(-0.5)*x + (1.5)*(x**2) + 30

plot(x,total,type = 'l',
[Link] = 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,

[Link](hsb, '[Link]', [Link]=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 [Link]
## 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
##

[Link] 7/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## Random effects:
## Groups Name Variance [Link].
## 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 ' ' 1

Computing ICC (and function)

Cut and paste
icc <- 8.553/(8.553+39.148)
icc

## [1] 0.1793044

Or extract the information to compute

vars <- [Link](VarCorr(model0))[4]
total <- sum(vars)
tau00 <- vars[1,1]
icc <- tau00/total
icc

## [1] 0.1793109

Best yet, function since we do this over and over, use the following function that I wrote:
icc <- function(modl){
vars <- [Link](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.1793109

Random 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 [Link]
## 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 [Link].
## 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.003

icc(model1)

## [1] 0.1887876

# for easier comparison

screenreg(list(model0,model1),
[Link]=TRUE,
[Link]=c("Null","RI ses"),
[Link]=list("Random Intercept Models"=1:2))

[Link] 8/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
##
## ===============================================================
## 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.05

Random 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 [Link]
## 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 [Link]. 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.012

[Link](VarCorr(model2))[4]

## vcov
## 1 8.62127569
## 2 0.67826037
## 3 0.05031215
## 4 36.70001845

deviance(model2)

## [1] 46710.98

# for easier comparison with what we have already fit

screenreg(list(model0,model1,model2),
[Link]=TRUE,
[Link]=c("Null","cses", "ces"),
[Link]=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

[Link] 9/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## ==========================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05

Graphs 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 [Link](model2, level = 0:2): additional arguments to

## [Link] ignored: level

U0j <- random$id[,1]

U1j <- random$id[,2]

plot(hsb$mathach ~ hsb$cses, type = 'n',

[Link] = 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")
}

[Link] 10/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R

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',

[Link] = 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)
[Link] <- [Link](VarCorr(model1))[4]
[Link] <- rep(sum([Link]), 61)

[Link] 11/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
[Link] <- [Link](VarCorr(model2))[4]
[Link]<- [Link][1,] + 2*[Link][3,]*ses2 + [Link][2,]*(ses2**2) + [Link][4,]

plot(ses2,[Link],type = 'l',
[Link] = 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([Link]), 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 [Link]
## 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 [Link]. 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.012

[Link] 12/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
[Link](VarCorr(model.d1))[4]

## vcov
## 1 6.2628806
## 2 0.4800510
## 3 -0.3001918
## 4 35.6631590

n2lnlike <- -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),
[Link]=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.05

Plot 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')
)

[Link] 13/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R

Plot of variance function

ses2 <- seq(from=-2, to= 3, by=.1)
vars.d1 <- [Link](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',
[Link] = 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 <- [Link](hsb$female)

hsb$race <- [Link](hsb$minority)
hsb$boy <- [Link](-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))

[Link] 14/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## 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 [Link]
## 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 [Link]. 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
[Link](VarCorr(model.00))[4]

## vcov
## 1 6.2628806
## 2 0.4800510
## 3 -0.3001918
## 4 35.6631590

This 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 [Link]
## 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 [Link]. 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.014

[Link](VarCorr(model.b1))[4]

## vcov
## 1 6.9181349
## 2 0.4569283
## 3 0.7857687
## 4 -0.2066691
## 5 -0.9444765

[Link] 15/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## 6 -0.1727498
## 7 35.5342275

Alternative 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 [Link]
## 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 [Link]. 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.014

[Link](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.5342275

Another 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 [Link]
## 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 [Link]. 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.014

[Link](VarCorr(model.b2))[4]

[Link] 16/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## vcov
## 1 5.8148277
## 2 0.4570253
## 3 0.7858424
## 4 -0.3792555
## 5 0.1588552
## 6 0.1727952
## 7 35.5341774

Bad 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 [Link]
## 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 [Link]. 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 ?isSingular

[Link](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.5341774

This Works.
Why?
summary([Link] <- 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 [Link]
## 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 [Link]. Corr
## id (Intercept) 6.1723 2.484
## cses 0.3856 0.621 0.03

[Link] 17/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## 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.032

And 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 [Link]
## 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 [Link]. 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.032

You'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 [Link]
## 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 [Link]. 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

[Link] 18/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
##
## 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.080

[Link](VarCorr(model.3))[4]

## vcov
## 1 2.6443444
## 2 0.6496455
## 3 -0.2574434
## 4 36.7155992

Models 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 [Link]
## 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 [Link]. 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.004

p127 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 [Link]
## 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

[Link] 19/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
##
## Random effects:
## Groups Name Variance [Link]. 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 rescaling

We 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 [Link]
## 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 [Link]. 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.048

[Link](VarCorr(model.5))[4]

## vcov
## 1 1.95419176
[Link] 20/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## 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.33742752

Other 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 + minor
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 [Link]
## 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 [Link]. 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.085

[Link](VarCorr(model.6))[4]

[Link] 21/22
2/1/23, 10:10 AM Edpsy 587: Random Intercept and Slope Models using R
## vcov
## 1 1.22579595
## 2 0.16458404
## 3 0.77736156
## 4 0.15183165
## 5 -0.06946504
## 6 -0.13709806
## 7 35.46009963

Round 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 [Link]
## 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 [Link].
## 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.144

[Link](VarCorr(model.7))[4]

## vcov
## 1 1.300117
## 2 35.631569

[Link] 22/22