Multi-level Modeling with ELS:2002 data

In this exercise, I’ll be using the ELS:2002 dataset to perform multi-level modeling with geographic regions.

For this project, my continuous outcome variable will be “Base year reading scores” as defined in the ELS:2002 dataset as follows:

#reading scores
els$reading_score<-ifelse(els$bytxrstd==-8,NA,els$bytxrstd)

Reading scores are both continuous and (somewhat) normally distributed:

hist(els$reading_score)

The geographic regions of the base-year school samples is coded in the data as follows:

#geographic regions
els<-els %>% mutate(regions=case_when(.$byregion==1~"Northeast",
                                      .$byregion==2~"Midwest",
                                      .$byregion==3~"South",
                                      .$byregion==4~"West"))
els$regions<-as.factor(els$regions)

Are there differences in reading scores between regions?

We can run an ANOVA to find out if there are statistically significant differences in reading scores between these regions.

fit.an<-lm(reading_score~regions,els)
anova(fit.an)
## Analysis of Variance Table
## 
## Response: reading_score
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## regions      3   6707 2235.79  22.764 1.124e-14 ***
## Residuals 8987 882688   98.22                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As we can see, it does appear that there is some significant differences between the regions in regards to reading scores.

PART 1: Random Intercept Model

fit1<-lmer(reading_score~(1|regions), data=els)
arm::display(fit1, detail=T)
## lme4::lmer(formula = reading_score ~ (1 | regions), data = els)
##             coef.est coef.se t value
## (Intercept) 51.85     0.56   92.69  
## 
## Error terms:
##  Groups   Name        Std.Dev.
##  regions  (Intercept) 1.10    
##  Residual             9.91    
## ---
## number of obs: 8991, groups: regions, 4
## AIC = 66776.4, DIC = 66771.4
## deviance = 66770.9
meanscore<-mean(els$reading_score, na.rm=T)
sdscore<-sd(els$reading_score, na.rm=T)
els$pred_score<-sdscore*(fitted(fit1))+meanscore
rate<-fixef(fit1)+ranef(fit1)$regions
est<-data.frame(rate =rate, id=rownames(ranef(fit1)$regions))
regionmeans<-aggregate(cbind(reading_score, pred_score)~regions,els, mean)

The random intercept model has a standard deviation of 1.10 (among regions) and a residual of 9.91.

We can see the relationship between reading scores and predicted scores here:

plot(reading_score~pred_score,regionmeans)

Now we can add the individual-level characteristics to the model:

fit2<-lmer(reading_score~non_trad_family+repeated_grade+male+race_eth+(1|regions),els)
anova(fit1,fit2)
## refitting model(s) with ML (instead of REML)
## Data: els
## Models:
## object: reading_score ~ (1 | regions)
## ..1: reading_score ~ non_trad_family + repeated_grade + male + race_eth + 
## ..1:     (1 | regions)
##        Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## object  3 66777 66798 -33385    66771                             
## ..1    11 65266 65344 -32622    65244 1526.7      8  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As we can see, the individual-level characteristics do indeed improve the model.

Do reading scores have a significant variation across regions?

rand(fit2)
## Analysis of Random effects Table:
##         Chi.sq Chi.DF p.value    
## regions   13.1      1   3e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

They sure do!

Random slope model

Here, I’ll be creating a random slope model with the “male” variable which will be varying across regions:

fit.slope<-lmer(reading_score~non_trad_family+repeated_grade+male+race_eth+(male|regions), data=els, REML=F)
display(fit.slope, detail=T)
## lme4::lmer(formula = reading_score ~ non_trad_family + repeated_grade + 
##     male + race_eth + (male | regions), data = els, REML = F)
##                         coef.est coef.se t value
## (Intercept)              55.72     0.32  174.91 
## non_trad_family          -2.48     0.20  -12.15 
## repeated_grade           -6.71     0.31  -21.81 
## male                     -0.96     0.20   -4.77 
## race_ethasian_pacific    -1.98     0.38   -5.20 
## race_ethblack            -5.97     0.33  -18.32 
## race_ethhispanic         -6.19     0.29  -21.03 
## race_ethmultirace        -1.23     0.48   -2.56 
## race_ethnative_american  -6.33     1.15   -5.51 
## 
## Error terms:
##  Groups   Name        Std.Dev. Corr  
##  regions  (Intercept) 0.53           
##           male        0.12     -1.00 
##  Residual             9.11           
## ---
## number of obs: 8991, groups: regions, 4
## AIC = 65269.9, DIC = 65243.9
## deviance = 65243.9

The results are largely the same as before: non-traditional family, repeated grades, maleness and membership in a minority group tend to negatively impact reading scores, even across disparate regions. Interestingly, maleness has a perfect inverse coorelation with reading scores across regions, leading me to wonder if something in the model is off, somehow.

anova(fit1,fit.slope)
## refitting model(s) with ML (instead of REML)
## Data: els
## Models:
## object: reading_score ~ (1 | regions)
## ..1: reading_score ~ non_trad_family + repeated_grade + male + race_eth + 
## ..1:     (male | regions)
##        Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)    
## object  3 66777 66798 -33385    66771                            
## ..1    13 65270 65362 -32622    65244  1527     10  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We do see, however, that the Random Slope model does provide a better fit than the random intercept model from part 1, although not significantly different from the random intercept model plus individual-level predictors from part 1.