Module 9: Binary/Generalized Linear Mixed Models

1. Treating “readprof” as the DV, and “schcode” as the level-2 clustering variable, estimate a null model. Compute, report, and interpret the ICC. Is it worth conducting MLM on these data?

The ICC is 0.06489173. Since it is above 0.05, then we would move forward with MLM.

2. Interpret the estimated fixed effect of the intercept (_cons) in the null model. What is this number measuring- what does it mean? You can be more vague with this answer! Logits are tough.

The logit is 0.1316. This is the odds of a student being proficient in reading.

Part 2: Adding Level-1 Covariates

3. Conduct a multilevel model with “readprof” as the DV, “lowses”, “female”, and “minor” as level-1 IVs. Include all three IVs as fixed effects (use a random intercept only). Interpret your results in odds ratios using the or option in Stata or the formula I provided in R…

Odds Ratio = 1.477114. This means that for a student who does not have a low socioeconomic status, is male, and is a nonminority, the predicted probability associated with reading proficiency has the odds of 1.48 to 1, which equates to the P(non lowses, male, and non-minority) = 0.5967.

Part 2: Adding Level-2 Covariates

4. Building from the previous section, test the impact of “smallsch” as a level-2 IVs, keeping all IVs as fixed effects (random intercept only). Interpret the results for the smallsch coefficient in terms of odds ratios. If you were a supporter of creating smaller schools, how would you interpret this result? (0.13854+ 0.32449 = .46303) The odds ratio is exp(.46303) = 1.588881. This would suggest that a small school increases the chance of proficient readers.

5. Next, add a second level-2 IV, schcomp to the model along with all variables from #4. Interpret the smallsch and schcomp coefficients. Return to your “small schools” argument- has it changed? The coefficients of smallsch is 0.10372. The coefficient of schcomp is -.20725. The significance of smallsch has changed with introducing schcomp. It does not appear to be significant in the model.

Part 3: Calculating Predicted Probabilities (See separate Stata and R directions)

6. R folks: Using the model you ran for #3 above, calculate the predicted probability of passing for students with the following characteristics: 1) A male student, who qualifies as low SES, who identifies as a member of a minority racial/ethnic group; and 2) A female student, who does not qualify as low SES, who does not identify as a member of a minority racial/ethnic group.

P(male, low SES, minority) = exp(.39009+0.32357)/(1+exp(.39009+.32357)) = exp(.71366)/(1+exp(.71366)) = 0.6712094

P(female, no low SES, nonminority) = exp(.39009-.436850-.43906)/(1+exp(.39009-.436850-.43906)) = exp(-0.48582)/(1+exp(-0.48582)) = 0.3809

Load in Our MVP Packages

library(tidyverse)
library(lme4)
library(Hmisc)

nsch.2 <- haven::read_dta("readprof.dta")
glimpse(nsch.2)

Rows: 6,528
Columns: 8
$ schcode  <dbl> 100, 100, 100, 100, 100, 100, 100, ...
$ lowses   <dbl+lbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ female   <dbl+lbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ minor    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ schcomp  <dbl> -2.4, -2.4, -2.4, -2.4, -2.4, -2.4,...
$ smallsch <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
$ id       <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ...
$ readprof <dbl> 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0,...

nsch.2.clean <- nsch.2 %>%
   mutate(.,
         lowses.fac = as.numeric(lowses), 
      female.fac = as.numeric(female), 
      minor.fac = as.numeric(minor), 
      smalsch.fac = as.numeric(smallsch), 
      readprof.fac = as.numeric(readprof), 
         )

glimpse(nsch.2.clean)

Rows: 6,528
Columns: 13
$ schcode      <dbl> 100, 100, 100, 100, 100, 100, 1...
$ lowses       <dbl+lbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
$ female       <dbl+lbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
$ minor        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ schcomp      <dbl> -2.4, -2.4, -2.4, -2.4, -2.4, -...
$ smallsch     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ id           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ...
$ readprof     <dbl> 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1...
$ lowses.fac   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ female.fac   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ minor.fac    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ smalsch.fac  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ readprof.fac <dbl> 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1...

Part 2: Multilevel Logistic Regression Models

model.null <- glmer(readprof.fac ~ (1|schcode), family = binomial, data = nsch.2.clean)
summary(model.null)

Generalized linear mixed model fit by maximum
  likelihood (Laplace Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: readprof.fac ~ (1 | schcode)
   Data: nsch.2.clean

     AIC      BIC   logLik deviance df.resid 
  8868.5   8882.1  -4432.3   8864.5     6526 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7663 -0.9815  0.6879  0.9145  1.3931 

Random effects:
 Groups  Name        Variance Std.Dev.
 schcode (Intercept) 0.2283   0.4778  
Number of obs: 6528, groups:  schcode, 122

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  0.13155    0.05198   2.531   0.0114 *
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Calculate ICC for logistic regression (L1 random variance is fixed)

null.icc <- 0.2283/(0.2283 + (pi^2/3))
null.icc

[1] 0.06489173

L1 predictors

model.2 <- glmer(readprof.fac ~ lowses.fac + female.fac + minor.fac +(1|schcode), family = binomial, data = nsch.2.clean)
summary(model.2)

Generalized linear mixed
  model fit by maximum
  likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: 
readprof.fac ~ lowses.fac + female.fac + minor.fac + (1 | schcode)
   Data: nsch.2.clean

     AIC      BIC   logLik 
  8696.8   8730.7  -4343.4 
deviance df.resid 
  8686.8     6523 

Scaled residuals: 
    Min      1Q  Median      3Q 
-1.9203 -0.9482  0.6194  0.8937 
    Max 
 1.6114 

Random effects:
 Groups  Name        Variance
 schcode (Intercept) 0.1033  
 Std.Dev.
 0.3214  
Number of obs: 6528, groups:  
schcode, 122

Fixed effects:
            Estimate Std. Error
(Intercept)  0.39009    0.05876
lowses.fac  -0.43685    0.05717
female.fac   0.32357    0.05168
minor.fac   -0.43906    0.05605
            z value Pr(>|z|)
(Intercept)   6.639 3.15e-11
lowses.fac   -7.642 2.14e-14
female.fac    6.260 3.84e-10
minor.fac    -7.833 4.77e-15
               
(Intercept) ***
lowses.fac  ***
female.fac  ***
minor.fac   ***
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’
  0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) lwss.f fml.fc
lowses.fac -0.359              
female.fac -0.415 -0.031       
minor.fac  -0.393 -0.171 -0.004

odds.ratio <- exp(.39009)
odds.ratio 

[1] 1.477114

exp(0.39009)/(1+exp(0.39009))

[1] 0.5963044

model.3 <- glmer(readprof.fac ~ lowses.fac + female.fac + minor.fac + smalsch.fac +(1|schcode), family = binomial, data = nsch.2.clean)
summary(model.3)

Generalized linear mixed model fit by
  maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: 
readprof.fac ~ lowses.fac + female.fac + minor.fac + smalsch.fac +  
    (1 | schcode)
   Data: nsch.2.clean

     AIC      BIC   logLik deviance 
  8695.8   8736.5  -4341.9   8683.8 
df.resid 
    6522 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9364 -0.9556  0.6210  0.8947  1.6222 

Random effects:
 Groups  Name        Variance Std.Dev.
 schcode (Intercept) 0.09668  0.3109  
Number of obs: 6528, groups:  
schcode, 122

Fixed effects:
            Estimate Std. Error z value
(Intercept)  0.32449    0.07001   4.635
lowses.fac  -0.43920    0.05713  -7.688
female.fac   0.32292    0.05167   6.249
minor.fac   -0.43983    0.05598  -7.857
smalsch.fac  0.13854    0.07941   1.745
            Pr(>|z|)    
(Intercept) 3.57e-06 ***
lowses.fac  1.50e-14 ***
female.fac  4.12e-10 ***
minor.fac   3.94e-15 ***
smalsch.fac   0.0811 .  
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’
  0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) lwss.f fml.fc mnr.fc
lowses.fac  -0.303                     
female.fac  -0.346 -0.031              
minor.fac   -0.336 -0.172 -0.003       
smalsch.fac -0.552 -0.002 -0.003  0.009

0.13854+  0.32449

[1] 0.46303

exp(.46303)

[1] 1.588881

Schcomp.

model.4 <- glmer(readprof.fac ~ lowses.fac + female.fac + minor.fac + smalsch.fac + schcomp +(1|schcode), family = binomial, data = nsch.2.clean)
summary(model.4)

Generalized linear mixed model fit by
  maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: 
readprof.fac ~ lowses.fac + female.fac + minor.fac + smalsch.fac +  
    schcomp + (1 | schcode)
   Data: nsch.2.clean

     AIC      BIC   logLik deviance 
  8665.5   8713.0  -4325.8   8651.5 
df.resid 
    6521 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0503 -0.9646  0.6064  0.8824  1.6455 

Random effects:
 Groups  Name        Variance Std.Dev.
 schcode (Intercept) 0.05658  0.2379  
Number of obs: 6528, groups:  
schcode, 122

Fixed effects:
            Estimate Std. Error z value
(Intercept)  0.28348    0.06401   4.429
lowses.fac  -0.36153    0.05778  -6.257
female.fac   0.32237    0.05161   6.246
minor.fac   -0.39549    0.05580  -7.087
smalsch.fac  0.10372    0.07008   1.480
schcomp     -0.20725    0.03511  -5.902
            Pr(>|z|)    
(Intercept) 9.48e-06 ***
lowses.fac  3.92e-10 ***
female.fac  4.21e-10 ***
minor.fac   1.37e-12 ***
smalsch.fac    0.139    
schcomp     3.59e-09 ***
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’
  0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) lwss.f fml.fc mnr.fc
lowses.fac  -0.329                     
female.fac  -0.377 -0.031              
minor.fac   -0.366 -0.159 -0.004       
smalsch.fac -0.500 -0.024 -0.005 -0.002
schcomp      0.154 -0.278 -0.011 -0.175
            smlsc.
lowses.fac        
female.fac        
minor.fac         
smalsch.fac       
schcomp      0.083

.39009+0.32357

[1] 0.71366

exp(.71366)/(1+exp(.71366))

[1] 0.6712094

exp(-0.8724)/(1+exp(-0.8724))