Soc 712 Homework #8b - Random Effects Model/Multilevel Model
Raven Shan
Introduction
For this assignment, I will be utilizing a dataset retrieved from the Bristol University website, which can be accessed through this link. The study was conducted on 3,435 Scottish students who attended 148 primary schools and 19 secondary schools. The outcome variable of interest is the students’ verbal reasoning scores from tests they took upon entering secondary school.
This assignment involves examining information on 2 conceptual levels; the first level is the student or individual, and the 2nd level I selected was the primary school. To do so, I will be conducting several different regression models starting with complete-pooling and no-pooling models, then random intercept and random slope models using the lme4 package.
This assignment has two objectives. First, it seeks to examine how students’ verbal reasoning scores fluctuate as a function of their gender. In other words, what is the effect of the student’s gender on their scores? Does one gender score higher than the other? If so, to what extent? The second objective is to examine the role that the student’s primary school plays in this analysis. This will provide meaningful context as students within the same school tend to be more similar than students that are not. It is necessary to reflect this added complexity. Further in my analysis, I will also be incorporating another independent variable, social class, to my regression model to examine its impact.
Data and Variables
The variables used in this analysis are as follows:
- vr_score: This variable serves as the outcome variable of interest. It contains each student’s verbal reasoning test score.
- primID: The primary school identifying code
- sex: This binary variable represents the students sex (0 = male, 1 = female).
- class: This variable measures the students’ social class using a continuous scale score (from low to high social class).
str(studentdataset)'data.frame': 3435 obs. of 4 variables:
$ vr_score: num 11 0 -14 -6 -30 -17 -17 -11 -9 -19 ...
$ primID : num 1 1 1 1 1 1 1 1 1 1 ...
$ sex : int 0 1 0 0 1 1 1 0 0 0 ...
$ class : num 0 0 0 20 0 0 0 0 0 0 ...Complete-Pooling Model
First, I conducted a complete-pooling model. This model is oversimplified as it treats all students and their circumstances the same, and disregards any school-level factors that may provide further insight.
This complete-pooling model estimates the effect of a student’s gender on their verbal reasoning test score. The results suggest that on average, females score higher than males by 2.55 points. Males receive an average verbal reading test score of -3.45. The relationship is statistically significant. However, as previously stated, this model is not ideal as it does not take into account important between-school differences.
cpooling <- lm(vr_score ~ sex, data = studentdataset)
summary(cpooling)
Call:
lm(formula = vr_score ~ sex, data = studentdataset)
Residuals:
Min 1Q Median 3Q Max
-29.097 -9.097 -0.097 8.903 43.456
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.4560 0.3173 -10.891 < 2e-16 ***
sex 2.5527 0.4516 5.652 1.71e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 13.23 on 3433 degrees of freedom
Multiple R-squared: 0.009221, Adjusted R-squared: 0.008932
F-statistic: 31.95 on 1 and 3433 DF, p-value: 1.712e-08No-Pooling Model
Next, I conducted a no-pooling model. The no-pooling model has its shortcomings as separate regression models must be conducted for each school. In this case, this would require conducting 148 different regression models (as determined by the output below), which is clearly too complicated and cumbersome. Instead, to make the process more practical, I will use some dyplr tools to plot 148 estimated intercepts and then 148 estimated slopes, which will ultimately allow me to examine the distribution of the parameters.
length(unique(studentdataset$primID))[1] 148Intercept
The intercept histogram below shows the average verbal reading test scores for boys in each of 148 different schools. Each school has its own unique intercept. The results range from -30 to approximately 28.
dcoef <- studentdataset %>%
group_by(primID) %>%
do(mod = lm(vr_score ~ sex, data = .))
coef <- dcoef %>% do(data.frame(intc = coef(.$mod)[1]))
ggplot(coef, aes(x = intc)) + geom_histogram()
Slope
Next, I used the slope to examine the estimated gender differences in verbal reasoning test scores among the 148 different primary schools. The results range from around -27 to 40.
dcoef <- studentdataset %>%
group_by(primID) %>%
do(mod = lm(vr_score ~ sex, data = .))
coef <- dcoef %>% do(data.frame(sexc = coef(.$mod)[2]))
ggplot(coef, aes(x = sexc)) + geom_histogram()
After examining the distribution of both the estimated intercepts and slopes, it appears to be some between-school variation that warrants further examination.
Partial-Pooling/Multilevel Modeling
It is clear that the previous approaches are not sufficient. The complete-pooling model ignores 2nd level context, while the no-pooling model requires extreme measures. Instead, the most efficient and effective method of simplifying this multilevel data for better understanding is the partial-pooling model. They key to this model is its ability to normalize a skewed distribution, making it that much easier to analyze.
Multilevel Model 1: Random Intercept Model
Here, the intercept represents the standard deviation of the (hypothetical) normal distribution. It can only be assumed that the regression intercept between the 148 schools in this dataset are normally distributed with a mean of 0 and a standard deviation of 4.14. Still, this method is insufficient as it does not allow the gender difference in verbal reasoning scores to vary between primary schools.
m1_lme <- lme4::lmer(vr_score ~ sex + (1|primID), data = studentdataset)
summary(m1_lme)Linear mixed model fit by REML ['lmerMod']
Formula: vr_score ~ sex + (1 | primID)
Data: studentdataset
REML criterion at convergence: 27318.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.7758 -0.6878 0.0160 0.6513 3.2733
Random effects:
Groups Name Variance Std.Dev.
primID (Intercept) 17.18 4.145
Residual 158.86 12.604
Number of obs: 3435, groups: primID, 148
Fixed effects:
Estimate Std. Error t value
(Intercept) -3.4151 0.4852 -7.039
sex 2.4311 0.4356 5.581
Correlation of Fixed Effects:
(Intr)
sex -0.445Multilevel Model 2: Random Slope Model
The limitation in the previous (random intercept) model was that it did not allow for gender differences to fluctuate among primary schools. To make this possible, instead, I utilized the random slope model and replaced the 1 with the variable sex. This allows me to estimate both the random intercept for each school as well as the random slope for sex differences for each school.
In the model below, I now have the standard deviation for the normal distribution for the sex slope. The first step is to examine the fixed effects in order to determine the common trend among schools. The intercept suggests that for those who came from a typical primary school, the average verbal reasoning test score for boys is -3.41. In this same typical school, girls tend to score higher by 2.42 points. Then, to determine the predicted verbal reasoning score for a particular primary school (instead of just the overall trend), I would need to evaluate not only the fixed effects, but also the random effects. The random effects produce one unique numerical value for each primary school, which can then be added to the fixed effects or common trend.
m2_lme <- lme4::lmer(vr_score ~ sex + (sex|primID), data = studentdataset)
summary(m2_lme)Linear mixed model fit by REML ['lmerMod']
Formula: vr_score ~ sex + (sex | primID)
Data: studentdataset
REML criterion at convergence: 27316.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.7209 -0.6960 0.0181 0.6564 3.2969
Random effects:
Groups Name Variance Std.Dev. Corr
primID (Intercept) 14.3734 3.7912
sex 0.5389 0.7341 1.00
Residual 158.7140 12.5982
Number of obs: 3435, groups: primID, 148
Fixed effects:
Estimate Std. Error t value
(Intercept) -3.4109 0.4600 -7.415
sex 2.4275 0.4404 5.512
Correlation of Fixed Effects:
(Intr)
sex -0.348Multilevel Model 3: Adding another covariate
The following model incorporates the student-level independent variable, social class (“class”). The results suggests that the student’s social class has a statistically significant effect on their verbal reasoning score. In a typical school, one unit increase in social class increases the student’s score by 0.18.
m3_lme <- lme4::lmer(vr_score ~ sex + class + (sex|primID), data = studentdataset)
summary(m3_lme)Linear mixed model fit by REML ['lmerMod']
Formula: vr_score ~ sex + class + (sex | primID)
Data: studentdataset
REML criterion at convergence: 27239
Scaled residuals:
Min 1Q Median 3Q Max
-2.6305 -0.7052 0.0000 0.6574 3.4535
Random effects:
Groups Name Variance Std.Dev. Corr
primID (Intercept) 11.8464 3.4419
sex 0.4674 0.6837 1.00
Residual 155.6067 12.4742
Number of obs: 3435, groups: primID, 148
Fixed effects:
Estimate Std. Error t value
(Intercept) -4.7455 0.4592 -10.333
sex 2.5692 0.4355 5.899
class 0.1879 0.0203 9.253
Correlation of Fixed Effects:
(Intr) sex
sex -0.369
class -0.314 0.033Multilevel Model 4: Incorporating an interaction term
The following model tests the interaction between sex and social class. In other words, I am testing whether sex differences in verbal reasoning scores depend on the students’ social class. The results show no statistical significance.
m4_lme <- lmer(vr_score ~ sex*class + (sex|primID), data = studentdataset)
summary(m4_lme)Linear mixed model fit by REML t-tests use Satterthwaite approximations to
degrees of freedom [lmerMod]
Formula: vr_score ~ sex * class + (sex | primID)
Data: studentdataset
REML criterion at convergence: 27243.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.6438 -0.7071 0.0014 0.6592 3.4485
Random effects:
Groups Name Variance Std.Dev. Corr
primID (Intercept) 11.9412 3.4556
sex 0.4267 0.6533 1.00
Residual 155.6528 12.4761
Number of obs: 3435, groups: primID, 148
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -4.69805 0.47891 174.00000 -9.810 < 2e-16 ***
sex 2.47238 0.51179 1891.00000 4.831 1.47e-06 ***
class 0.18124 0.02745 3166.00000 6.602 4.75e-11 ***
sex:class 0.01423 0.03959 3350.00000 0.360 0.719
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) sex class
sex -0.451
class -0.409 0.375
sex:class 0.277 -0.526 -0.673Model Selection
Based on the AIC values, the third multilevel model proves to be the best fitting.
AIC(cpooling, m1_lme, m2_lme, m3_lme, m4_lme)htmlreg(list(cpooling, m1_lme, m2_lme, m3_lme, m4_lme))| Model 1 | Model 2 | Model 3 | Model 4 | Model 5 | ||
|---|---|---|---|---|---|---|
| (Intercept) | -3.46*** | -3.42*** | -3.41*** | -4.75*** | -4.70*** | |
| (0.32) | (0.49) | (0.46) | (0.46) | (0.48) | ||
| sex | 2.55*** | 2.43*** | 2.43*** | 2.57*** | 2.47*** | |
| (0.45) | (0.44) | (0.44) | (0.44) | (0.51) | ||
| class | 0.19*** | 0.18*** | ||||
| (0.02) | (0.03) | |||||
| sex:class | 0.01 | |||||
| (0.04) | ||||||
| R2 | 0.01 | |||||
| Adj. R2 | 0.01 | |||||
| Num. obs. | 3435 | 3435 | 3435 | 3435 | 3435 | |
| RMSE | 13.23 | |||||
| AIC | 27326.83 | 27328.78 | 27252.96 | 27259.45 | ||
| BIC | 27351.39 | 27365.63 | 27295.95 | 27308.59 | ||
| Log Likelihood | -13659.41 | -13658.39 | -13619.48 | -13621.73 | ||
| Num. groups: primID | 148 | 148 | 148 | 148 | ||
| Var: primID (Intercept) | 17.18 | 14.37 | 11.85 | 11.94 | ||
| Var: Residual | 158.86 | 158.71 | 155.61 | 155.65 | ||
| Var: primID sex | 0.54 | 0.47 | 0.43 | |||
| Cov: primID (Intercept) sex | 2.78 | 2.35 | 2.26 | |||
| p < 0.001, p < 0.01, p < 0.05 | ||||||
Random Effects Simulation with merTools
##estimating the model
m5_mer <- lme4::lmer(vr_score ~ sex + class + (1|primID), data = studentdataset)
fastdisp(m5_mer)lme4::lmer(formula = vr_score ~ sex + class + (1 | primID), data = studentdataset)
coef.est coef.se
(Intercept) -4.75 0.48
sex 2.57 0.43
class 0.19 0.02
Error terms:
Groups Name Std.Dev.
primID (Intercept) 3.77
Residual 12.48
---
number of obs: 3435, groups: primID, 148
AIC = 27250.6The following code produces a dataframe with estimates of the values of each of the random effects.
reEx <- REsim(m5_mer)
head(reEx)p1 <- plotREsim(reEx)
p1
As the above graph depicts, there is nearly equal variation in both directions. The primary school effect ranges from about -7 to 8.
Next, I generated confidence intervals, which provides a range of values that is expected to contain the value of parameters of interest. If we compare the output below with the estimates, this is indeed the case.
confint.merMod(m5_mer,method="profile")Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 3.0749544 4.5249501
.sigma 12.1807868 12.7834904
(Intercept) -5.6947345 -3.8067922
sex 1.7248589 3.4159230
class 0.1485141 0.2286467Overall, there is not a considerably wide variation in the random effects.
Centering Covariates
In the code below, I am centering the independent variable, ‘class’, which turns the mean value to 0 (Ideally, centering variables is best for models with interaction terms). I then refit the model and create a table to compare coefficients between the first model (with the original ‘class’ variable) and the second model (with the centered ‘class’ variable).
studentdataset %<>% mutate(c_class = class - mean(class))m3_lme2 <- lmer(vr_score ~ sex + (sex|primID) + c_class, data = studentdataset)
htmlreg(list(m3_lme, m3_lme2))| Model 1 | Model 2 | ||
|---|---|---|---|
| (Intercept) | -4.75*** | -3.46*** | |
| (0.46) | (0.44) | ||
| sex | 2.57*** | 2.57*** | |
| (0.44) | (0.44) | ||
| class | 0.19*** | ||
| (0.02) | |||
| c_class | 0.19*** | ||
| (0.02) | |||
| AIC | 27252.96 | 27252.96 | |
| BIC | 27295.95 | 27295.95 | |
| Log Likelihood | -13619.48 | -13619.48 | |
| Num. obs. | 3435 | 3435 | |
| Num. groups: primID | 148 | 148 | |
| Var: primID (Intercept) | 11.85 | 11.85 | |
| Var: primID sex | 0.47 | 0.47 | |
| Cov: primID (Intercept) sex | 2.35 | 2.35 | |
| Var: Residual | 155.61 | 155.61 | |
| p < 0.001, p < 0.01, p < 0.05 | |||
Based on the above table, although some coefficients differ, the AIC and BIC values are the exact same, indicating both models fit the data equally well. Essentially, the interpretation and fitness remains the same despite the covariate being centered in the latter model.
Intra-class Correlation (ICC)
The intra-class correlation (ICC) indicates how strongly the 2nd level effect is. The random effect for the intercept, 17.61, represents the amount of variation that exists at the primary school level. The residual value (160.19) represents the amount of variation at the student level.
m0_lme <-lme4::lmer(vr_score ~ 1 + (1|primID), data = studentdataset)
summary(m0_lme)Linear mixed model fit by REML ['lmerMod']
Formula: vr_score ~ 1 + (1 | primID)
Data: studentdataset
REML criterion at convergence: 27350
Scaled residuals:
Min 1Q Median 3Q Max
-2.8374 -0.6762 -0.0034 0.6627 3.2920
Random effects:
Groups Name Variance Std.Dev.
primID (Intercept) 17.61 4.196
Residual 160.19 12.657
Number of obs: 3435, groups: primID, 148
Fixed effects:
Estimate Std. Error t value
(Intercept) -2.2113 0.4388 -5.039The ICC is obtained by first adding the school variance and the residual variance, then dividing the school variance by the number obtained in the previous computation: 17.61/(17.61+160.19) =0.099. An ICC value of .099 suggests that nearly 10% of the variation in students’ verbal reasoning test scores can be attributed to differences between primary schools. Ultimately, while the primary school does contribute to a degree, the majority of the variation is due to individual-level differences between students.