ACS Y10 NAS: Multilevel Model Unpacked
Complete pooling, no pooling and partial pooling (multilevel) Poisson models for Q1_male, boys rating boys:
# load the data and create necessary variables and subsets:
library(foreign)
acsY10nas <- read.spss("ACS_all_schools to 2013_merge_6.sav", use.value.labels = FALSE,
to.data.frame = TRUE)
library(car)
acsY10nas$SEXMAIN <- recode(acsY10nas$SEXMAIN, "1=0;2=1")
acsY10nas$znas8 <- scale(acsY10nas$nas8)
acsY10nas$zhope <- scale(acsY10nas$hpmean_10_2012)
acsY10nas$zswb <- scale(acsY10nas$swbmean_10_2012)
acsY10nas$nasT <- recode(acsY10nas$nas8, "0:4=0; 4.01:6=1")
acsY10nas$schoolMAIN <- as.factor(acsY10nas$schoolMAIN)
Boys <- subset(acsY10nas, acsY10nas$SEXMAIN == "0", select = c(1:3046))
Girls <- subset(acsY10nas, acsY10nas$SEXMAIN == "1", select = c(1:3046))
Boys$zboys <- scale(Boys$n_boys)
Boys$zgirls <- scale(Boys$n_girls)
Girls$zboys <- scale(Girls$n_boys)
Girls$zgirls <- scale(Girls$n_girls)
# confirm that schoolMAIN is a factor (not numeric or integer)
class(Boys$schoolMAIN)
## [1] "factor"
# Complete pooling model: group indicator (schoolMAIN) not included in the
# model
summary(glm(Q1_male ~ znas8 + zboys, family = poisson, data = Boys))
##
## Call:
## glm(formula = Q1_male ~ znas8 + zboys, family = poisson, data = Boys)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.439 -0.963 -0.128 0.687 4.127
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7677 0.0223 34.46 < 2e-16 ***
## znas8 0.0937 0.0224 4.19 2.8e-05 ***
## zboys 0.0754 0.0218 3.46 0.00055 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1508.3 on 946 degrees of freedom
## Residual deviance: 1476.9 on 944 degrees of freedom
## (448 observations deleted due to missingness)
## AIC: 3600
##
## Number of Fisher Scoring iterations: 5
# No pooling model: separate models are fitted with each school
# the -1 after schoolMAIN factor removes the constant term, so school 1 is
# included as well
summary(glm(Q1_male ~ znas8 + zboys + schoolMAIN - 1, family = poisson, data = Boys))
##
## Call:
## glm(formula = Q1_male ~ znas8 + zboys + schoolMAIN - 1, family = poisson,
## data = Boys)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.447 -0.961 -0.168 0.665 4.130
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## znas8 0.0885 0.0225 3.93 8.7e-05 ***
## zboys 0.1052 0.0528 1.99 0.046 *
## schoolMAIN1 0.7220 0.0868 8.32 < 2e-16 ***
## schoolMAIN2 0.9216 0.1056 8.73 < 2e-16 ***
## schoolMAIN3 0.7655 0.0922 8.30 < 2e-16 ***
## schoolMAIN5 0.7059 0.0884 7.99 1.4e-15 ***
## schoolMAIN6 0.7307 0.0763 9.58 < 2e-16 ***
## schoolMAIN7 0.7930 0.0922 8.60 < 2e-16 ***
## schoolMAIN8 0.7875 0.0962 8.19 2.7e-16 ***
## schoolMAIN10 0.8314 0.0761 10.93 < 2e-16 ***
## schoolMAIN11 0.5784 0.0951 6.08 1.2e-09 ***
## schoolMAIN20 0.7663 0.0892 8.59 < 2e-16 ***
## schoolMAIN22 0.8427 0.1105 7.63 2.4e-14 ***
## schoolMAIN28 0.7966 0.1734 4.59 4.4e-06 ***
## schoolMAIN36 0.8718 0.1334 6.54 6.3e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2490.4 on 947 degrees of freedom
## Residual deviance: 1468.0 on 932 degrees of freedom
## (448 observations deleted due to missingness)
## AIC: 3616
##
## Number of Fisher Scoring iterations: 5
# Partial pooling: a compromise between complete pooling and no-pooling
# also known as multilevel modeling
library(lme4)
# Varying intercept with no predictors
# the predictor '1' after ~ allows the constant to vary by school
# (1 | schoolMAIN) is R's way of specifying varying intercept model
M0 <- glmer(Q1_male ~ 1 + (1 | schoolMAIN), family = poisson, data = Boys)
summary(M0)
## Generalized linear mixed model fit by maximum likelihood ['glmerMod']
## Family: poisson ( log )
## Formula: Q1_male ~ 1 + (1 | schoolMAIN)
## Data: Boys
##
## AIC BIC logLik deviance
## 3789 3799 -1893 3785
##
## Random effects:
## Groups Name Variance Std.Dev.
## schoolMAIN (Intercept) 0.0048 0.0693
## Number of obs: 993, groups: schoolMAIN, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7499 0.0298 25.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# to look at the estimated model within each school:
coef(M0)
## $schoolMAIN
## (Intercept)
## 1 0.8012
## 2 0.7757
## 3 0.7200
## 5 0.6927
## 6 0.7985
## 7 0.7312
## 8 0.7116
## 10 0.7843
## 11 0.6749
## 20 0.8358
## 22 0.7385
## 28 0.7455
## 36 0.7465
##
## attr(,"class")
## [1] "coef.mer"
# to get 'fixed effects' or the estimated model averaging over schools
fixef(M0)
## (Intercept)
## 0.7499
# to look at 'random effects' or the school-level errors
ranef(M0)
## $schoolMAIN
## (Intercept)
## 1 0.051242
## 2 0.025749
## 3 -0.029968
## 5 -0.057268
## 6 0.048559
## 7 -0.018694
## 8 -0.038327
## 10 0.034362
## 11 -0.074995
## 20 0.085884
## 22 -0.011439
## 28 -0.004482
## 36 -0.003476
##
## attr(,"class")
## [1] "ranef.mer"
# Varying intercept with a single predictor znas8
# note the slope is held constant here because not specified to vary by
# school
M1 <- glmer(Q1_male ~ znas8 + (1 | schoolMAIN), family = poisson, data = Boys)
summary(M1)
## Generalized linear mixed model fit by maximum likelihood ['glmerMod']
## Family: poisson ( log )
## Formula: Q1_male ~ znas8 + (1 | schoolMAIN)
## Data: Boys
##
## AIC BIC logLik deviance
## 3611 3626 -1803 3605
##
## Random effects:
## Groups Name Variance Std.Dev.
## schoolMAIN (Intercept) 0.00259 0.0509
## Number of obs: 947, groups: schoolMAIN, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7629 0.0269 28.39 < 2e-16 ***
## znas8 0.0950 0.0224 4.24 2.2e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## znas8 -0.100
# to look at the estimated model within each school:
coef(M1)
## $schoolMAIN
## (Intercept) znas8
## 1 0.7924 0.09502
## 2 0.7742 0.09502
## 3 0.7441 0.09502
## 5 0.7383 0.09502
## 6 0.7944 0.09502
## 7 0.7578 0.09502
## 8 0.7429 0.09502
## 10 0.7726 0.09502
## 11 0.7135 0.09502
## 20 0.8206 0.09502
## 22 0.7534 0.09502
## 28 0.7559 0.09502
## 36 0.7609 0.09502
##
## attr(,"class")
## [1] "coef.mer"
# to get 'fixed effects' or the estimated model averaging over schools
fixef(M1)
## (Intercept) znas8
## 0.76291 0.09502
# to look at 'random effects' or the school-level errors
ranef(M1)
## $schoolMAIN
## (Intercept)
## 1 0.029525
## 2 0.011276
## 3 -0.018776
## 5 -0.024567
## 6 0.031488
## 7 -0.005085
## 8 -0.020018
## 10 0.009681
## 11 -0.049417
## 20 0.057660
## 22 -0.009512
## 28 -0.007001
## 36 -0.002025
##
## attr(,"class")
## [1] "ranef.mer"
# Varying intercept with a single predictor znas8
# and adding nboys allowing its slope to vary
# 1 + zboys in (1 + zboys | schoolMAIN) is R's way of varying slope of zboys
# note that znas8 slope is still not allowed to vary
M2 <- glmer(Q1_male ~ znas8 + (1 + zboys | schoolMAIN), family = poisson, data = Boys)
summary(M2)
## Generalized linear mixed model fit by maximum likelihood ['glmerMod']
## Family: poisson ( log )
## Formula: Q1_male ~ znas8 + (1 + zboys | schoolMAIN)
## Data: Boys
##
## AIC BIC logLik deviance
## 3609 3634 -1800 3599
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schoolMAIN (Intercept) 0.00534 0.0731
## zboys 0.00767 0.0876 1.00
## Number of obs: 947, groups: schoolMAIN, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7218 0.0289 25.00 < 2e-16 ***
## znas8 0.0929 0.0224 4.14 3.4e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## znas8 -0.065
coef(M2) # estimated model within each school
## $schoolMAIN
## zboys (Intercept) znas8
## 1 0.062990 0.7743 0.09285
## 2 0.003385 0.7246 0.09285
## 3 -0.007920 0.7151 0.09285
## 5 -0.041466 0.6872 0.09285
## 6 0.055753 0.7683 0.09285
## 7 0.002422 0.7238 0.09285
## 8 0.066214 0.7770 0.09285
## 10 0.075061 0.7844 0.09285
## 11 -0.098425 0.6397 0.09285
## 20 0.079021 0.7877 0.09285
## 22 -0.006255 0.7165 0.09285
## 28 0.002022 0.7234 0.09285
## 36 0.021332 0.7395 0.09285
##
## attr(,"class")
## [1] "coef.mer"
fixef(M2) # estimated model averaging over schools
## (Intercept) znas8
## 0.72175 0.09285
ranef(M2) # school-level errors
## $schoolMAIN
## (Intercept) zboys
## 1 0.052534 0.062990
## 2 0.002823 0.003385
## 3 -0.006605 -0.007920
## 5 -0.034583 -0.041466
## 6 0.046498 0.055753
## 7 0.002020 0.002422
## 8 0.055223 0.066214
## 10 0.062602 0.075061
## 11 -0.082087 -0.098425
## 20 0.065904 0.079021
## 22 -0.005217 -0.006255
## 28 0.001686 0.002022
## 36 0.017791 0.021332
##
## attr(,"class")
## [1] "ranef.mer"
# Varying intercept AND slope for znas8
# finally allowing znas8 slope to vary!
# note zboys not included here (will be next)
M3 <- glmer(Q1_male ~ znas8 + (1 + znas8 | schoolMAIN), family = poisson, data = Boys)
# 1 + znas8 in (1 + znas8 | schoolMAIN) is R's way of varying slope of znas8
summary(M3)
## Generalized linear mixed model fit by maximum likelihood ['glmerMod']
## Family: poisson ( log )
## Formula: Q1_male ~ znas8 + (1 + znas8 | schoolMAIN)
## Data: Boys
##
## AIC BIC logLik deviance
## 3615 3639 -1803 3605
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schoolMAIN (Intercept) 0.002722 0.0522
## znas8 0.000133 0.0115 -1.00
## Number of obs: 947, groups: schoolMAIN, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7630 0.0271 28.19 < 2e-16 ***
## znas8 0.0971 0.0226 4.29 1.8e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## znas8 -0.186
coef(M3) # estimated model within each school
## $schoolMAIN
## (Intercept) znas8
## 1 0.7919 0.09070
## 2 0.7849 0.09226
## 3 0.7352 0.10321
## 5 0.7371 0.10280
## 6 0.8020 0.08847
## 7 0.7590 0.09796
## 8 0.7424 0.10163
## 10 0.7717 0.09517
## 11 0.7172 0.10720
## 20 0.8183 0.08487
## 22 0.7520 0.09951
## 28 0.7545 0.09897
## 36 0.7566 0.09850
##
## attr(,"class")
## [1] "coef.mer"
fixef(M3) # estimated model averaging over schools
## (Intercept) znas8
## 0.76301 0.09708
ranef(M3) # school-level errors
## $schoolMAIN
## (Intercept) znas8
## 1 0.028898 -0.0063783
## 2 0.021859 -0.0048247
## 3 -0.027781 0.0061318
## 5 -0.025889 0.0057141
## 6 0.039015 -0.0086114
## 7 -0.003994 0.0008815
## 8 -0.020586 0.0045436
## 10 0.008647 -0.0019086
## 11 -0.045826 0.0101147
## 20 0.055341 -0.0122146
## 22 -0.011011 0.0024303
## 28 -0.008539 0.0018847
## 36 -0.006432 0.0014197
##
## attr(,"class")
## [1] "ranef.mer"
# Varying intercept AND slope for znas8
# adding nboys allowing its slope to vary as well
# this is probably the final model, unless I've missed something else?
M4 <- glmer(Q1_male ~ znas8 + (1 + znas8 + zboys | schoolMAIN), family = poisson,
data = Boys)
summary(M4)
## Generalized linear mixed model fit by maximum likelihood ['glmerMod']
## Family: poisson ( log )
## Formula: Q1_male ~ znas8 + (1 + znas8 + zboys | schoolMAIN)
## Data: Boys
##
## AIC BIC logLik deviance
## 3615 3654 -1800 3599
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schoolMAIN (Intercept) 0.005327 0.0730
## znas8 0.000103 0.0101 1.00
## zboys 0.008258 0.0909 1.00 1.00
## Number of obs: 947, groups: schoolMAIN, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7213 0.0288 25.02 < 2e-16 ***
## znas8 0.0895 0.0225 3.97 7.1e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## znas8 -0.027
coef(M4) # estimated model within each school
## $schoolMAIN
## zboys (Intercept) znas8
## 1 0.064019 0.7728 0.09663
## 2 -0.014232 0.7099 0.08791
## 3 0.011286 0.7304 0.09075
## 5 -0.042033 0.6876 0.08481
## 6 0.054821 0.7654 0.09561
## 7 0.001153 0.7223 0.08962
## 8 0.070700 0.7781 0.09738
## 10 0.081315 0.7866 0.09856
## 11 -0.105094 0.6369 0.07778
## 20 0.080771 0.7862 0.09850
## 22 -0.004256 0.7179 0.08902
## 28 0.004903 0.7253 0.09004
## 36 0.028952 0.7446 0.09272
##
## attr(,"class")
## [1] "coef.mer"
fixef(M4) # estimated model averaging over schools
## (Intercept) znas8
## 0.7213 0.0895
ranef(M4) # school-level errors
## $schoolMAIN
## (Intercept) znas8 zboys
## 1 0.0514185 0.0071370 0.064019
## 2 -0.0114309 -0.0015866 -0.014232
## 3 0.0090643 0.0012581 0.011286
## 5 -0.0337597 -0.0046859 -0.042033
## 6 0.0440310 0.0061116 0.054821
## 7 0.0009261 0.0001285 0.001153
## 8 0.0567848 0.0078818 0.070700
## 10 0.0653098 0.0090651 0.081315
## 11 -0.0844093 -0.0117162 -0.105094
## 20 0.0648735 0.0090046 0.080771
## 22 -0.0034184 -0.0004745 -0.004256
## 28 0.0039381 0.0005466 0.004903
## 36 0.0232533 0.0032276 0.028952
##
## attr(,"class")
## [1] "ranef.mer"