In-class Exercises: Multilevel models
Chang-Chun Chen
2020-10-10
1 Inclass exercise 1
1.1 Info
- Problem: Replicate the results of analysis reported the high school and beyond - null 2 3 4 examples.
- Data info: The High Schools and Beyond data set consists of measures on 7,185 students in 160 schools. The individual level outcome is a mathematical achievement score. The individual level covariate is a measure of socio-economic status. There are two school-level predictors: one is the school sector; and the other is the mean socio-economic status of students in a school.
Column 1: school: School ID
Column 2: minrty: a factor with levels(Minority)
Column 3: sx: a factor with levels Male and Female
Column 4: ses: Socio-economic status of the student (centered at grand mean)/a numeric vector of socio-economic scores
Column 5: mAch: a numeric vector of Mathematical achievement score of the student
Column 6: meanses: Mean socio-economic status (centered at grand mean) /a numeric vector of mean ses for the school
Column 7: sector: School sectors, 0 = public schools, 1 = catholic schools
Column 8: cses: Centered soco-economic status about the school mean (a numeric vector of centered ses values where the centering is with respect to the meanses for the school.)
1.2 Data input
# school-effects model version 0
# load package
pacman::p_load(mlmRev, tidyverse, lme4, nlme)
# data from mlmRev package
data(Hsb82, package= "mlmRev")
# data structure
str(Hsb82)## 'data.frame': 7185 obs. of 8 variables:
## $ school : Ord.factor w/ 160 levels "8367"<"8854"<..: 59 59 59 59 59 59 59 59 59 59 ...
## $ minrty : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ sx : Factor w/ 2 levels "Male","Female": 2 2 1 1 1 1 2 1 2 1 ...
## $ ses : num -1.528 -0.588 -0.528 -0.668 -0.158 ...
## $ mAch : num 5.88 19.71 20.35 8.78 17.9 ...
## $ meanses: num -0.434 -0.434 -0.434 -0.434 -0.434 ...
## $ sector : Factor w/ 2 levels "Public","Catholic": 1 1 1 1 1 1 1 1 1 1 ...
## $ cses : num -1.0936 -0.1536 -0.0936 -0.2336 0.2764 ...# make a copy of data to the dta1 object
dta1 <- Hsb82
# rename: use sensible variable names
names(dta1) <- c("School", "Minority", "Gender", "SES", "Math",
"Ave_SES", "Sector", "C_SES")
# head of 6 list
head(dta1)## School Minority Gender SES Math Ave_SES Sector C_SES
## 1 1224 No Female -1.528 5.876 -0.434383 Public -1.09361702
## 2 1224 No Female -0.588 19.708 -0.434383 Public -0.15361702
## 3 1224 No Male -0.528 20.349 -0.434383 Public -0.09361702
## 4 1224 No Male -0.668 8.781 -0.434383 Public -0.23361702
## 5 1224 No Male -0.158 17.898 -0.434383 Public 0.27638298
## 6 1224 No Male 0.022 4.583 -0.434383 Public 0.45638298# set theme
ot <- theme_set(theme_bw())1.3 Descriptive statistics
各校平均數學分數的平均值=12.62, 所有學生數學分數的總平均值=12.748
各校平均數學分數的sd=3.117, 所有學生數學分數的sd=6.878
# math mean by school
# ave_math = 各校數學的平均分數, dta1_mmath =school id and mean math of the school
dta1_mmath <- dta1 %>%
group_by(School) %>%
mutate( ave_math = mean(Math)) %>%
select( School, ave_math ) %>%
unique %>%
as.data.frame
summary(dta1_mmath$ave_math)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.24 10.47 12.90 12.62 14.65 19.72sd(dta1_mmath$ave_math)## [1] 3.117651# grand mean
summary(dta1$Math)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.832 7.275 13.131 12.748 18.317 24.993# sd for overall Math scores
sd(dta1$Math)## [1] 6.8782461.4 Visualization of raw
1.4.1 Histograms 1
所有學生數學分數分布的狀況
# freedman-diaconis rule for binwidth (用freedman-diaconis rule來決定histogram 中box的寬度)
# (dim(dta1)[1]^(1/3)): dim(x)= x data的列數 欄數, dim(x)[1], 取x data的列數,
# bin width =2*IQR(x)/3√n, IQR(x) is the interquartile range of the data, n = the number of observations in the sample x
bw <- 2*IQR(dta1$Math)/(dim(dta1)[1]^(1/3))
# Histograms
ggplot(dta1, aes(Math)) +
stat_bin(aes(fill = ..count..), binwidth=bw) +
labs(x="Math achievement scores")
1.4.2 Box plots
取public and catholic schools 共60間，並畫出各校Math score的boxplot
# seed the random number to reproduce the sample(取樣的種子碼)
set.seed(20160923)
# pick 60 out of 160 (pick one third sample )
n <- sample(160, 60)
# box plots for a sample of 60 schools
ggplot(filter(dta1, School %in% levels(School)[n]),
aes(reorder(School, Math, mean), Math, color=Sector)) +
geom_boxplot() +
geom_point(size= rel(0.6)) +
labs(x="School ID", y="Math achievement score") +
theme(axis.text.x=element_text(size=6)) +
ggtitle("A sample of 60 schools")
1.4.3 Histograms 2
所有學校平均數學分數分布的狀況
# freedman-diaconis rule for binwidth
bw <- 2*IQR(dta1_mmath$ave_math)/(160^(1/3))
# Histograms
ggplot(dta1_mmath, aes(ave_math)) +
geom_histogram(aes(y= ..density.., alpha= .2), binwidth=bw) +
geom_density(fill="NA") +
geom_vline(xintercept=mean(dta1_mmath$ave_math), linetype="longdash") +
geom_text(aes(7, 0.17, label= "ave = 12.62, sd = 3.12")) +
labs(x="Math achievement school means", y="Density") +
theme(legend.position="NONE")
1.5 LMER
1.5.1 m0
School-effects on Math score
# baseline model
summary(m0 <- lmer(Math ~ (1 | School), data = dta1))## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ (1 | School)
## Data: dta1
##
## REML criterion at convergence: 47116.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0631 -0.7539 0.0267 0.7606 2.7426
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 8.614 2.935
## Residual 39.148 6.257
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.6370 0.2444 51.71# CIs for model parameters by boot
confint(m0, oldNames=F, method="boot")## Computing bootstrap confidence intervals ...## 2.5 % 97.5 %
## sd_(Intercept)|School 2.596281 3.333500
## sigma 6.151525 6.360499
## (Intercept) 12.128943 13.114645The estimated overall mean math score is 12.637, CI=(12.128, 13.114) ≈ 12.62 (the mean of all school).
The estimated between school variance in mean math score by school is 8.614.
The variation between student within the same school is estimated to be 39.148.
The between school variance is about 8.614/(8.614 + 39.148) = 0.18 (18%)of the total variance. This is the intraclass correlation.
The intraclass correlation would be one if all the student in a school has the same score and zero if there were no differences in mean math scores between schools.
# default m0 residual plot
plot(m0, xlab="Fitted values", ylab="Pearson residuals",
pch=20, cex=.5, type = c("p","g"))
1.5.2 m1
School level: School-effects on Math score
# random effects anova - intercept only model (same as m0)
summary(m1 <- lmer(Math ~ (1 | School), data = dta1))## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ (1 | School)
## Data: dta1
##
## REML criterion at convergence: 47116.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0631 -0.7539 0.0267 0.7606 2.7426
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 8.614 2.935
## Residual 39.148 6.257
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.6370 0.2444 51.711.5.3 m2
School level: fixed-effect is mean SES by school and School as random effects on Math score
# random intercept; fixed-effect is MeanSES
summary(m2 <- lmer(Math ~ Ave_SES + (1 | School), data = dta1))## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ Ave_SES + (1 | School)
## Data: dta1
##
## REML criterion at convergence: 46961.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.13493 -0.75254 0.02413 0.76766 2.78515
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 2.639 1.624
## Residual 39.157 6.258
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.6846 0.1493 84.97
## Ave_SES 5.8635 0.3615 16.22
##
## Correlation of Fixed Effects:
## (Intr)
## Ave_SES 0.010The estimated overall mean school math score, controlling for school’s mean SES, is about 12.6846.
For an unit increase in a school’s mean SES, his or her math score, on average, by about 5.8635 point.
The variation in math score between schools (given mean SES by school) is estimated to be about 2.639 < 8.614 (estimated by the null model). “Why?”
The variation in the math score among pupils within schools (given mean SES by school) is estimated to be about 39.157 > 39.148 (The estimated variation between student within the same school).
The intraclass correlation is 2.639/(2.639 + 39.157) = 0.06314001, e.g.: about 6% of total variance is accounted for by the between school variance (controlling for mean SES by school).
(8.614-2.639)/8.614=0.6936383, About 69.36% of between school variation in math score is explained by the school’s mean SES.
(39.148-39.157)/39.148 = -0.0002298968, About 0.02% of between pupil variation (in the same school)in math score is explained by school’s mean SES.
# residual plot
plot(m2, xlab = "Fitted values", ylab = "Pearson residuals",
main = "Model 2", pch = 20, cex = .5, type = c("p", "g"))
1.5.4 m0a
School level: Mean SES by school effects on Math
# one regression for all
summary(m0a <- lm(Math ~ Ave_SES, data = dta1))##
## Call:
## lm(formula = Math ~ Ave_SES, data = dta1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.4520 -4.8332 0.1735 5.1117 16.2787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.74703 0.07621 167.27 <2e-16 ***
## Ave_SES 5.71688 0.18429 31.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.46 on 7183 degrees of freedom
## Multiple R-squared: 0.1181, Adjusted R-squared: 0.118
## F-statistic: 962.3 on 1 and 7183 DF, p-value: < 2.2e-16The estimated overall mean math score is 12.74703 ≈ 12.748 (the mean of all pupils).
For an unit increase in a school mean SES, his or her math score, on average, by about 5.71688 point < 5.8635 point (m2) 差不多?.
#
theme_set(theme_bw())
# individual math scores vs mean school ses
ggplot(dta1, aes(Ave_SES, Math)) +
geom_smooth(method = "lm") +
geom_point(alpha = .2) +
labs(x = "School mean SES", y = "Math achievement score")## `geom_smooth()` using formula 'y ~ x'
The higher school mean SES, the higher math achievement score.
# math mean by school
dta1_mmath <- dta1 %>%
group_by(School) %>%
mutate( ave_math = mean(Math)) %>%
select( School, Ave_SES, ave_math ) %>%
unique %>%
as.data.frameave_math = the mean math score by school Ave_SES = the mean ses by school
# regression of mean math over mean ses by school
ggplot(dta1_mmath, aes(Ave_SES, ave_math)) +
geom_smooth(method = "lm") +
geom_point(alpha = .2)+
labs(x = "Mean School SES", y = "Mean School Math Scores")## `geom_smooth()` using formula 'y ~ x'
The higher school mean SES, the higher mean math achievement score by school (the CI of the regression line is wider than the previous one).
# see the random number generator
set.seed(20160923)
# take a sample of size 31 from total number of schools
n <- sample(length(levels(dta1$School)), 31)
#
theme_set(theme_bw())
#
ggplot(data = filter(dta1, School %in% levels(School)[n]),
aes(C_SES, Math, color = Sector)) +
geom_smooth(method = "lm", se = F, lwd = .6) +
geom_point(size = rel(0.5)) +
facet_wrap(~ School, nrow = 2) +
labs(x = "Centered SES", y="Math achievement score")+
theme(axis.text.x = element_text(size = 6))## `geom_smooth()` using formula 'y ~ x'
In 31 school, most of them showed the higher the centered ses the higher math score. Only 3 schools (3377, 3498, 4292) showed the opposite results.
1.5.5 m1b
Two-level model with one pupil’s-level covariate
?? the interactions of school and C_SES effects on the math score without intercept ?? (the result of “School / C_SES” same as “C_SES : School”)
##
# one regression line per school
summary(m1b <- lm(Math ~ School / C_SES - 1, data = dta1))##
## Call:
## lm(formula = Math ~ School/C_SES - 1, data = dta1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.2662 -4.3241 0.1164 4.4434 18.3617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## School8367 4.55279 1.61953 2.811 0.004950 **
## School8854 4.23978 1.07122 3.958 7.64e-05 ***
## School4458 5.81140 0.87465 6.644 3.28e-11 ***
## School5762 4.32486 0.99621 4.341 1.44e-05 ***
## School6990 5.97679 0.83237 7.180 7.68e-13 ***
## School5815 7.27136 1.21194 6.000 2.08e-09 ***
## School7172 8.06682 0.91354 8.830 < 2e-16 ***
## School4868 12.31018 1.03923 11.845 < 2e-16 ***
## School7341 9.79418 0.84853 11.543 < 2e-16 ***
## School1358 11.20623 1.10635 10.129 < 2e-16 ***
## School4383 11.46568 1.21194 9.461 < 2e-16 ***
## School2305 11.13776 0.74031 15.045 < 2e-16 ***
## School8800 7.33594 1.07122 6.848 8.13e-12 ***
## School3088 9.14585 0.97033 9.425 < 2e-16 ***
## School8775 9.46700 0.87465 10.824 < 2e-16 ***
## School7890 8.34110 0.84853 9.830 < 2e-16 ***
## School6144 8.54507 0.92410 9.247 < 2e-16 ***
## School6443 9.47553 1.10635 8.565 < 2e-16 ***
## School5192 10.40950 1.14518 9.090 < 2e-16 ***
## School6808 9.28611 0.91354 10.165 < 2e-16 ***
## School2818 13.87283 0.93504 14.837 < 2e-16 ***
## School9340 11.17855 1.12526 9.934 < 2e-16 ***
## School4523 8.35174 0.88390 9.449 < 2e-16 ***
## School6816 14.53824 0.81709 17.793 < 2e-16 ***
## School2277 9.29761 0.77587 11.983 < 2e-16 ***
## School8009 14.08472 0.88390 15.935 < 2e-16 ***
## School5783 13.15034 1.12526 11.686 < 2e-16 ***
## School3013 12.61083 0.83237 15.151 < 2e-16 ***
## School7101 11.84993 1.14518 10.348 < 2e-16 ***
## School4530 9.05570 0.76345 11.861 < 2e-16 ***
## School9021 14.69666 0.80976 18.149 < 2e-16 ***
## School4511 13.40903 0.79568 16.852 < 2e-16 ***
## School2639 6.61548 0.93504 7.075 1.64e-12 ***
## School3377 9.18671 0.90333 10.170 < 2e-16 ***
## School6578 11.99400 0.80976 14.812 < 2e-16 ***
## School9347 13.53875 0.80263 16.868 < 2e-16 ***
## School3705 10.33169 0.90333 11.437 < 2e-16 ***
## School3533 10.40904 0.87465 11.901 < 2e-16 ***
## School1296 7.63596 0.87465 8.730 < 2e-16 ***
## School4350 11.85542 1.05486 11.239 < 2e-16 ***
## School9397 10.35547 0.88390 11.716 < 2e-16 ***
## School4253 9.41286 0.79568 11.830 < 2e-16 ***
## School2655 12.34517 0.84033 14.691 < 2e-16 ***
## School7342 11.16641 0.79568 14.034 < 2e-16 ***
## School9292 10.27963 1.39020 7.394 1.59e-13 ***
## School3499 13.27653 0.98302 13.506 < 2e-16 ***
## School7364 14.17214 0.91354 15.513 < 2e-16 ***
## School8983 10.99200 0.84853 12.954 < 2e-16 ***
## School5650 14.27353 0.90333 15.801 < 2e-16 ***
## School2658 13.39616 0.90333 14.830 < 2e-16 ***
## School8188 12.74097 1.10635 11.516 < 2e-16 ***
## School4410 13.47298 0.94637 14.236 < 2e-16 ***
## School9508 13.57466 1.02428 13.253 < 2e-16 ***
## School8707 12.88394 0.87465 14.730 < 2e-16 ***
## School1499 7.66036 0.83237 9.203 < 2e-16 ***
## School8477 12.52224 0.99621 12.570 < 2e-16 ***
## School1288 13.51080 1.21194 11.148 < 2e-16 ***
## School6291 10.10731 1.02428 9.868 < 2e-16 ***
## School1224 9.71545 0.88390 10.992 < 2e-16 ***
## School4292 12.86435 0.75162 17.116 < 2e-16 ***
## School8857 15.29694 0.75747 20.195 < 2e-16 ***
## School3967 12.03508 0.84033 14.322 < 2e-16 ***
## School6415 11.86020 0.82462 14.383 < 2e-16 ***
## School1317 13.17769 0.87465 15.066 < 2e-16 ***
## School2629 14.90777 0.80263 18.574 < 2e-16 ***
## School4223 14.62262 0.90333 16.187 < 2e-16 ***
## School1462 10.49556 0.80263 13.076 < 2e-16 ***
## School9550 11.08914 1.12526 9.855 < 2e-16 ***
## School6464 7.09162 1.12526 6.302 3.12e-10 ***
## School4931 13.79081 0.79568 17.332 < 2e-16 ***
## School5937 16.77597 1.12526 14.908 < 2e-16 ***
## School7919 14.84997 0.99621 14.906 < 2e-16 ***
## School3716 10.36866 0.94637 10.956 < 2e-16 ***
## School1909 14.42332 1.14518 12.595 < 2e-16 ***
## School2651 11.08432 0.98302 11.276 < 2e-16 ***
## School2467 10.14752 0.84033 12.076 < 2e-16 ***
## School1374 9.72846 1.14518 8.495 < 2e-16 ***
## School6600 11.70389 0.80976 14.453 < 2e-16 ***
## School5667 13.77823 0.77587 17.758 < 2e-16 ***
## School5720 14.28230 0.83237 17.159 < 2e-16 ***
## School3498 16.39045 0.83237 19.691 < 2e-16 ***
## School3881 11.94922 0.94637 12.626 < 2e-16 ***
## School2995 9.54611 0.89346 10.684 < 2e-16 ***
## School5838 13.68961 1.08836 12.578 < 2e-16 ***
## School3688 14.65626 0.92410 15.860 < 2e-16 ***
## School9158 8.54517 0.83237 10.266 < 2e-16 ***
## School8946 10.37509 0.79568 13.039 < 2e-16 ***
## School7232 12.54263 0.84033 14.926 < 2e-16 ***
## School2917 7.97895 0.92410 8.634 < 2e-16 ***
## School6170 14.18105 1.32234 10.724 < 2e-16 ***
## School8165 16.45122 0.86567 19.004 < 2e-16 ***
## School9104 16.83211 0.81709 20.600 < 2e-16 ***
## School2030 12.07819 0.88390 13.665 < 2e-16 ***
## School8150 14.85236 0.91354 16.258 < 2e-16 ***
## School4042 14.31542 0.75747 18.899 < 2e-16 ***
## School8357 14.38185 1.16619 12.332 < 2e-16 ***
## School8531 13.52868 0.94637 14.295 < 2e-16 ***
## School6074 13.77909 0.80976 17.016 < 2e-16 ***
## School4420 13.87416 1.07122 12.952 < 2e-16 ***
## School1906 15.98317 0.83237 19.202 < 2e-16 ***
## School3992 14.64521 0.83237 17.595 < 2e-16 ***
## School3999 10.94404 0.89346 12.249 < 2e-16 ***
## School4173 12.72466 0.91354 13.929 < 2e-16 ***
## School4325 13.24000 0.83237 15.906 < 2e-16 ***
## School5761 11.13806 0.84033 13.254 < 2e-16 ***
## School6484 12.91240 1.02428 12.606 < 2e-16 ***
## School6897 15.09763 0.86567 17.440 < 2e-16 ***
## School7635 15.06553 0.84853 17.755 < 2e-16 ***
## School7734 10.55964 1.29194 8.173 3.54e-16 ***
## School8175 11.69809 1.05486 11.090 < 2e-16 ***
## School8874 12.05503 1.00995 11.936 < 2e-16 ***
## School9225 14.66733 1.00995 14.523 < 2e-16 ***
## School2458 13.98568 0.80263 17.425 < 2e-16 ***
## School3610 15.35495 0.75747 20.271 < 2e-16 ***
## School5640 13.16011 0.80263 16.396 < 2e-16 ***
## School3838 16.06281 0.82462 19.479 < 2e-16 ***
## School9359 15.27062 0.83237 18.346 < 2e-16 ***
## School2208 15.40467 0.78231 19.691 < 2e-16 ***
## School6089 15.56958 1.05486 14.760 < 2e-16 ***
## School1477 14.22847 0.76959 18.488 < 2e-16 ***
## School2768 10.88692 1.21194 8.983 < 2e-16 ***
## School3039 16.96386 1.32234 12.829 < 2e-16 ***
## School5819 12.13890 0.85697 14.165 < 2e-16 ***
## School6397 12.79610 0.78231 16.357 < 2e-16 ***
## School1308 16.25550 1.35500 11.997 < 2e-16 ***
## School1433 19.71914 1.02428 19.252 < 2e-16 ***
## School1436 18.11161 0.91354 19.826 < 2e-16 ***
## School1461 16.84264 1.05486 15.967 < 2e-16 ***
## School1637 7.02411 1.16619 6.023 1.80e-09 ***
## School1942 18.11090 1.12526 16.095 < 2e-16 ***
## School1946 12.90844 0.97033 13.303 < 2e-16 ***
## School2336 16.51770 0.88390 18.687 < 2e-16 ***
## School2526 17.05300 0.80263 21.246 < 2e-16 ***
## School2626 13.39661 0.98302 13.628 < 2e-16 ***
## School2755 16.47651 0.88390 18.641 < 2e-16 ***
## School2771 11.84411 0.81709 14.495 < 2e-16 ***
## School2990 18.44792 0.87465 21.092 < 2e-16 ***
## School3020 14.39527 0.78891 18.247 < 2e-16 ***
## School3152 13.20904 0.84033 15.719 < 2e-16 ***
## School3332 14.27816 0.98302 14.525 < 2e-16 ***
## School3351 11.46518 0.97033 11.816 < 2e-16 ***
## School3427 19.71559 0.86567 22.775 < 2e-16 ***
## School3657 9.52118 0.84853 11.221 < 2e-16 ***
## School4642 14.59903 0.77587 18.816 < 2e-16 ***
## School5404 15.41498 0.80263 19.206 < 2e-16 ***
## School5619 15.41624 0.74590 20.668 < 2e-16 ***
## School6366 15.65640 0.79568 19.677 < 2e-16 ***
## School6469 18.45572 0.80263 22.994 < 2e-16 ***
## School7011 13.81358 1.05486 13.095 < 2e-16 ***
## School7276 12.67940 0.83237 15.233 < 2e-16 ***
## School7332 14.63610 0.87465 16.734 < 2e-16 ***
## School7345 11.33855 0.80976 14.002 < 2e-16 ***
## School7688 18.42231 0.82462 22.340 < 2e-16 ***
## School7697 15.72178 1.07122 14.677 < 2e-16 ***
## School8193 16.23226 0.92410 17.565 < 2e-16 ***
## School8202 11.71243 1.02428 11.435 < 2e-16 ***
## School8627 10.88372 0.83237 13.076 < 2e-16 ***
## School8628 16.52838 0.77587 21.303 < 2e-16 ***
## School9198 19.09229 1.08836 17.542 < 2e-16 ***
## School9586 14.86369 0.78891 18.841 < 2e-16 ***
## School8367:C_SES 0.25037 2.22488 0.113 0.910403
## School8854:C_SES 1.93884 1.35428 1.432 0.152292
## School4458:C_SES 1.13184 1.36793 0.827 0.408034
## School5762:C_SES -1.01410 1.95950 -0.518 0.604803
## School6990:C_SES 0.94769 0.98169 0.965 0.334393
## School5815:C_SES 3.01800 2.16485 1.394 0.163335
## School7172:C_SES 0.99448 1.36612 0.728 0.466662
## School4868:C_SES 1.28647 1.48571 0.866 0.386576
## School7341:C_SES 1.70370 1.00528 1.695 0.090169 .
## School1358:C_SES 5.06801 1.71168 2.961 0.003078 **
## School4383:C_SES 6.18019 2.11601 2.921 0.003504 **
## School2305:C_SES -0.78211 1.13763 -0.687 0.491795
## School8800:C_SES 2.56813 1.42668 1.800 0.071893 .
## School3088:C_SES 1.79134 1.23144 1.455 0.145807
## School8775:C_SES 1.00146 1.32013 0.759 0.448114
## School7890:C_SES -0.65597 1.44460 -0.454 0.649784
## School6144:C_SES 2.77027 1.30158 2.128 0.033340 *
## School6443:C_SES -0.74335 1.46027 -0.509 0.610734
## School5192:C_SES 1.60349 1.46102 1.098 0.272452
## School6808:C_SES 2.27610 1.15285 1.974 0.048384 *
## School2818:C_SES 2.80241 1.52348 1.839 0.065888 .
## School9340:C_SES 3.30950 2.00186 1.653 0.098334 .
## School4523:C_SES 2.38079 1.24545 1.912 0.055971 .
## School6816:C_SES 1.35272 1.16781 1.158 0.246767
## School2277:C_SES -2.01503 1.17745 -1.711 0.087062 .
## School8009:C_SES 1.55687 1.58675 0.981 0.326542
## School5783:C_SES 3.11406 1.49402 2.084 0.037165 *
## School3013:C_SES 3.83980 1.75094 2.193 0.028341 *
## School7101:C_SES 1.29545 1.42842 0.907 0.364484
## School4530:C_SES 1.64743 1.23926 1.329 0.183771
## School9021:C_SES 2.52416 1.39478 1.810 0.070384 .
## School4511:C_SES 0.04251 1.38066 0.031 0.975438
## School2639:C_SES -0.63010 1.52971 -0.412 0.680417
## School3377:C_SES -0.74685 1.31955 -0.566 0.571419
## School6578:C_SES 2.39054 1.30725 1.829 0.067492 .
## School9347:C_SES 2.68599 1.18051 2.275 0.022920 *
## School3705:C_SES 1.15850 1.33049 0.871 0.383933
## School3533:C_SES -0.31177 1.64469 -0.190 0.849658
## School1296:C_SES 1.07596 1.36610 0.788 0.430948
## School4350:C_SES 4.37192 1.59458 2.742 0.006127 **
## School9397:C_SES 2.44644 1.47583 1.658 0.097430 .
## School4253:C_SES -0.39954 1.15329 -0.346 0.729023
## School2655:C_SES 5.22704 1.35364 3.861 0.000114 ***
## School7342:C_SES 1.01246 1.42097 0.713 0.476173
## School9292:C_SES 0.75868 2.74811 0.276 0.782501
## School3499:C_SES 0.99238 1.58184 0.627 0.530445
## School7364:C_SES 0.25950 1.82833 0.142 0.887139
## School8983:C_SES 1.38594 1.34299 1.032 0.302119
## School5650:C_SES 0.68062 1.17460 0.579 0.562308
## School2658:C_SES 2.62990 1.42677 1.843 0.065335 .
## School8188:C_SES 4.39759 1.65556 2.656 0.007920 **
## School4410:C_SES 2.76020 1.55700 1.773 0.076312 .
## School9508:C_SES 3.95379 1.82178 2.170 0.030019 *
## School8707:C_SES 3.39153 1.09903 3.086 0.002037 **
## School1499:C_SES 3.63473 1.16428 3.122 0.001805 **
## School8477:C_SES 3.81216 1.25548 3.036 0.002403 **
## School1288:C_SES 3.25545 1.84816 1.761 0.078205 .
## School6291:C_SES 3.98087 1.69353 2.351 0.018770 *
## School1224:C_SES 2.50858 1.42433 1.761 0.078243 .
## School4292:C_SES -0.16061 1.16329 -0.138 0.890195
## School8857:C_SES 0.80222 1.43452 0.559 0.576025
## School3967:C_SES 3.31107 1.35069 2.451 0.014255 *
## School6415:C_SES 3.53008 1.21524 2.905 0.003686 **
## School1317:C_SES 1.27391 1.58930 0.802 0.422837
## School2629:C_SES 0.22235 1.14645 0.194 0.846225
## School4223:C_SES 2.48659 1.50770 1.649 0.099140 .
## School1462:C_SES -0.82881 1.39839 -0.593 0.553411
## School9550:C_SES 3.89194 1.45938 2.667 0.007675 **
## School6464:C_SES 1.00349 1.82459 0.550 0.582348
## School4931:C_SES 0.91185 1.17597 0.775 0.438133
## School5937:C_SES 1.03962 2.00531 0.518 0.604175
## School7919:C_SES 3.98937 1.88178 2.120 0.034042 *
## School3716:C_SES 5.86379 1.18088 4.966 7.01e-07 ***
## School1909:C_SES 2.85479 1.63722 1.744 0.081260 .
## School2651:C_SES 4.89906 1.55631 3.148 0.001652 **
## School2467:C_SES 3.13713 1.65830 1.892 0.058563 .
## School1374:C_SES 3.85432 1.66836 2.310 0.020904 *
## School6600:C_SES 4.70429 1.16489 4.038 5.44e-05 ***
## School5667:C_SES 3.52297 1.24071 2.839 0.004532 **
## School5720:C_SES 2.46631 1.26524 1.949 0.051302 .
## School3498:C_SES -0.13109 1.45712 -0.090 0.928320
## School3881:C_SES 2.39071 1.79271 1.334 0.182388
## School2995:C_SES 1.43231 1.01804 1.407 0.159492
## School5838:C_SES 1.85305 1.77992 1.041 0.297872
## School3688:C_SES 1.53672 1.72035 0.893 0.371749
## School9158:C_SES 3.86121 1.12532 3.431 0.000604 ***
## School8946:C_SES 1.69048 1.06749 1.584 0.113331
## School7232:C_SES 5.00160 1.47738 3.385 0.000715 ***
## School2917:C_SES 1.13585 1.02982 1.103 0.270083
## School6170:C_SES 4.81178 1.64962 2.917 0.003547 **
## School8165:C_SES 1.80224 1.38480 1.301 0.193148
## School9104:C_SES 1.49398 1.83440 0.814 0.415430
## School2030:C_SES 1.41198 1.31420 1.074 0.282681
## School8150:C_SES -0.18571 1.51064 -0.123 0.902162
## School4042:C_SES 1.69362 1.20501 1.405 0.159922
## School8357:C_SES 2.67578 1.15717 2.312 0.020788 *
## School8531:C_SES 3.31823 1.40287 2.365 0.018043 *
## School6074:C_SES 1.52909 1.30292 1.174 0.240603
## School4420:C_SES 2.95866 1.58878 1.862 0.062614 .
## School1906:C_SES 2.14551 1.36955 1.567 0.117259
## School3992:C_SES 0.53788 1.38802 0.388 0.698389
## School3999:C_SES 3.56698 0.98836 3.609 0.000310 ***
## School4173:C_SES 3.36567 1.45956 2.306 0.021143 *
## School4325:C_SES 2.75605 1.03273 2.669 0.007632 **
## School5761:C_SES 3.10801 1.19136 2.609 0.009106 **
## School6484:C_SES 0.60568 1.49351 0.406 0.685093
## School6897:C_SES 3.58049 1.17477 3.048 0.002314 **
## School7635:C_SES 2.44847 1.26927 1.929 0.053767 .
## School7734:C_SES 6.03523 1.48469 4.065 4.86e-05 ***
## School8175:C_SES 1.61237 1.54334 1.045 0.296186
## School8874:C_SES 4.09630 1.43512 2.854 0.004326 **
## School9225:C_SES 2.88589 1.39642 2.067 0.038806 *
## School2458:C_SES 2.95669 1.22988 2.404 0.016241 *
## School3610:C_SES 2.95585 1.20868 2.446 0.014489 *
## School5640:C_SES 3.82774 1.38890 2.756 0.005868 **
## School3838:C_SES 0.59790 1.23391 0.485 0.628005
## School9359:C_SES -0.83348 1.46679 -0.568 0.569895
## School2208:C_SES 2.63664 1.31898 1.999 0.045648 *
## School6089:C_SES 1.69245 1.51926 1.114 0.265318
## School1477:C_SES 1.23061 1.14871 1.071 0.284075
## School2768:C_SES 3.51228 1.37517 2.554 0.010668 *
## School3039:C_SES 2.95567 1.92045 1.539 0.123839
## School5819:C_SES 1.97252 1.43955 1.370 0.170658
## School6397:C_SES 2.75901 1.03849 2.657 0.007908 **
## School1308:C_SES 0.12602 2.89741 0.043 0.965308
## School1433:C_SES 1.85429 1.74798 1.061 0.288809
## School1436:C_SES 1.60056 1.61653 0.990 0.322150
## School1461:C_SES 6.26650 1.54684 4.051 5.15e-05 ***
## School1637:C_SES 3.11681 1.56230 1.995 0.046081 *
## School1942:C_SES 0.08938 2.17762 0.041 0.967260
## School1946:C_SES 3.58583 1.40776 2.547 0.010881 *
## School2336:C_SES 1.90497 1.39610 1.364 0.172457
## School2526:C_SES 0.15950 1.32609 0.120 0.904264
## School2626:C_SES 4.09968 1.77861 2.305 0.021197 *
## School2755:C_SES 0.56050 1.41945 0.395 0.692950
## School2771:C_SES 4.26819 1.60528 2.659 0.007859 **
## School2990:C_SES 1.32454 1.32327 1.001 0.316883
## School3020:C_SES 1.65368 1.21797 1.358 0.174595
## School3152:C_SES 2.76825 1.00158 2.764 0.005727 **
## School3332:C_SES 2.03095 1.71772 1.182 0.237106
## School3351:C_SES 2.45504 1.16971 2.099 0.035867 *
## School3427:C_SES -0.48817 1.31793 -0.370 0.711094
## School3657:C_SES 3.73591 1.10047 3.395 0.000691 ***
## School4642:C_SES 3.27239 1.17344 2.789 0.005306 **
## School5404:C_SES 1.21423 1.60712 0.756 0.449955
## School5619:C_SES 5.25753 1.25841 4.178 2.98e-05 ***
## School6366:C_SES 1.51752 1.21527 1.249 0.211812
## School6469:C_SES 1.75529 1.29272 1.358 0.174562
## School7011:C_SES 5.07465 1.72585 2.940 0.003289 **
## School7276:C_SES 3.77336 1.06358 3.548 0.000391 ***
## School7332:C_SES 2.46320 1.32282 1.862 0.062634 .
## School7345:C_SES 4.21192 0.98954 4.256 2.10e-05 ***
## School7688:C_SES 0.11634 1.47469 0.079 0.937119
## School7697:C_SES 3.13622 1.77439 1.767 0.077190 .
## School8193:C_SES 2.33521 1.62482 1.437 0.150704
## School8202:C_SES 3.70590 1.47606 2.511 0.012073 *
## School8627:C_SES 1.86956 1.18737 1.575 0.115408
## School8628:C_SES 1.23139 1.34892 0.913 0.361339
## School9198:C_SES 2.61055 1.69902 1.537 0.124462
## School9586:C_SES 1.67208 1.33730 1.250 0.211217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.06 on 6865 degrees of freedom
## Multiple R-squared: 0.8328, Adjusted R-squared: 0.825
## F-statistic: 106.8 on 320 and 6865 DF, p-value: < 2.2e-16# collect regression estimates from individual regression by school
m1b_betas <- m1b %>%
coef %>%
matrix(ncol = 2) %>%
as.data.frame %>%
rename( b0 = V1, b1 = V2) b0 = coefficients of m1b
b1 = residuals of m1b
#
colMeans(m1b_betas)## b0 b1
## 12.620755 2.201641the mean coefficients of m1b = 12.620755
the mean residuals of m1b = 2.201641
#
apply(m1b_betas, 2, sd)## b0 b1
## 3.117651 1.630998the sd of coefficients of m1b = 3.117651
the sd of residuals of m1b = 1.630998
# the correlations between residual and coefficients
cor(m1b_betas)## b0 b1
## b0 1.00000000 0.01346686
## b1 0.01346686 1.00000000The correlations between residual and coefficients is 0.0134 (1.34%) intraclass correlation?
# scatter plot of regression coefficient estimates
ggplot(m1b_betas, aes(b0, b1)) +
geom_point(pch = 20) +
stat_smooth(method = "lm", se = F) +
stat_ellipse() +
annotate("text", x = 5, y = 5, label = "r = 0.013") +
labs(x = "Intercept estimates", y = "Slope estimates")## `geom_smooth()` using formula 'y ~ x'
1.5.6 m3a
Two-level model with one pupil’s-level covariate fixed-effect is centered SES and School as random effects on Math score
# random intercept only
summary(m3a <- lmer(Math ~ C_SES + (1 | School), data = dta1))## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ C_SES + (1 | School)
## Data: dta1
##
## REML criterion at convergence: 46724
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0969 -0.7322 0.0194 0.7572 2.9147
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 8.672 2.945
## Residual 37.010 6.084
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.6361 0.2445 51.68
## C_SES 2.1912 0.1087 20.17
##
## Correlation of Fixed Effects:
## (Intr)
## C_SES 0.000The estimated overall mean school math score, controlling for centered SES, is about 12.6361.
For an unit increase in centered SES, his or her math score, on average, by about 2.1912 point.
The variation in math score between schools (given centered SES) is estimated to be about 8.672 > 8.614 (estimated by the null model).
The variation in the math score among pupils within schools (given centered SES) is estimated to be about 37.010 < 39.148 (The estimated variation between student within the same school).
The intraclass correlation is 8.672/(8.672 + 37.010) = 0.1898341, e.g.: about 18.9% of total variance is accounted for by the between school variance (controlling for centered SES).
(12.6361-2.1912)/12.6361=0.8265921, About 82.65% of between school variation in math score is explained by the centered SES.
(39.148-37.010)/39.148 = 0.05461326, About 5.5% of between pupil variation (in the same school)in math score is explained by centered SES.
1.5.7 m3b
fixed-effect is centered SES and the random effect of centered SES in School level on Math score
# random intercept and slope - correlated
summary(m3b <- lmer(Math ~ C_SES + (C_SES | School), data = dta1))## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ C_SES + (C_SES | School)
## Data: dta1
##
## REML criterion at convergence: 46714.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09680 -0.73193 0.01855 0.75386 2.89924
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## School (Intercept) 8.681 2.9464
## C_SES 0.694 0.8331 0.02
## Residual 36.700 6.0581
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.6362 0.2445 51.68
## C_SES 2.1932 0.1283 17.10
##
## Correlation of Fixed Effects:
## (Intr)
## C_SES 0.0091.5.8 m3c
# random intercept and slope - uncorrelated
summary(m3c <- lme (Math ~ C_SES, random = list(School = pdDiag(~ C_SES)), data = dta1))## Linear mixed-effects model fit by REML
## Data: dta1
## AIC BIC logLik
## 46724.25 46758.64 -23357.12
##
## Random effects:
## Formula: ~C_SES | School
## Structure: Diagonal
## (Intercept) C_SES Residual
## StdDev: 2.946303 0.8331144 6.058068
##
## Fixed effects: Math ~ C_SES
## Value Std.Error DF t-value p-value
## (Intercept) 12.636017 0.2445004 7024 51.68098 0
## C_SES 2.192649 0.1282614 7024 17.09516 0
## Correlation:
## (Intr)
## C_SES 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.09607430 -0.73095310 0.01694466 0.75390887 2.90167264
##
## Number of Observations: 7185
## Number of Groups: 160summary(m3c <- lmer (Math ~ C_SES + (1 | School) + (C_SES - 1 | School), data = dta1))## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ C_SES + (1 | School) + (C_SES - 1 | School)
## Data: dta1
##
## REML criterion at convergence: 46714.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09607 -0.73095 0.01694 0.75391 2.90167
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 8.6807 2.9463
## School.1 C_SES 0.6941 0.8331
## Residual 36.7002 6.0581
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.6360 0.2445 51.68
## C_SES 2.1926 0.1283 17.09
##
## Correlation of Fixed Effects:
## (Intr)
## C_SES 0.000# default residual plot
plot(m3c, xlab = "Fitted values", ylab = "Pearson residuals",
pch = 20, cex = .5, type = c("p", "g"))
不確定(C_SES - 1 | School)的意思是否等同 (C_SES | School)，output 結果好像一樣，以下解釋可能有誤
The similarity of pupils in the same C_SES of the same schools is estimated to be about 19% ≈ 8.6807/46.075= 0.1884037
The similarity of C_SES in the same schools is estimated to be about 0.3% ≈ 0.6941/46.075= 0.002771568
The similarity of pupils in the same schools is estimated to be about 80% ≈ 36.7002/46.075= 0.7965317
1.5.9 model comparisons
# model comparisons
anova(m3a, m3b, m3c)## refitting model(s) with ML (instead of REML)## Data: dta1
## Models:
## m3a: Math ~ C_SES + (1 | School)
## m3c: Math ~ C_SES + (1 | School) + (C_SES - 1 | School)
## m3b: Math ~ C_SES + (C_SES | School)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m3a 4 46728 46756 -23360 46720
## m3c 5 46721 46755 -23356 46711 9.4196 1 0.002147 **
## m3b 6 46723 46764 -23356 46711 0.0134 1 0.907954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sjPlot::tab_model(m3a, m3b, m3c, show.p=FALSE, show.r2=FALSE, show.obs=FALSE, show.ngroups=FALSE, show.se=TRUE, show.ci=FALSE)| Math | Math | Math | ||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Estimates | std. Error | Estimates | std. Error |
| (Intercept) | 12.64 | 0.24 | 12.64 | 0.24 | 12.64 | 0.24 |
| C_SES | 2.19 | 0.11 | 2.19 | 0.13 | 2.19 | 0.13 |
| Random Effects | ||||||
| σ2 | 37.01 | 36.70 | 36.70 | |||
| τ00 | 8.67 School | 8.68 School | 8.68 School | |||
| 0.69 School.1 | ||||||
| τ11 | 0.69 School.C_SES | |||||
| ρ01 | 0.02 School | |||||
| ICC | 0.19 | 0.20 | 0.19 | |||
fixed effect m3a same as m3b, ICC almost the same m3a m3b are similar and differ with m3c.
# CIs for model parameter
confint(m3c, oldNames=F, method="boot")## Computing bootstrap confidence intervals ...##
## 3 message(s): boundary (singular) fit: see ?isSingular
## 1 warning(s): Model failed to converge with max|grad| = 0.00255542 (tol = 0.002, component 1)## 2.5 % 97.5 %
## sd_(Intercept)|School 2.6100481 3.269883
## sd_C_SES|School 0.3848003 1.141918
## sigma 5.9425454 6.150370
## (Intercept) 12.1304752 13.088948
## C_SES 1.9525498 2.439487m3a m3b are similar and differ with m3c.
# Math vs SES by sector
ggplot(dta1, aes(C_SES, Math, color = Sector)) +
geom_smooth(method = "lm") +
geom_jitter(alpha = 0.2) +
labs(x = "Centered SES", y = "Math achievement score") +
theme(legend.position = c(.1, .8))## `geom_smooth()` using formula 'y ~ x'
# to facilitate comparison with sas output
# options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly"))The effect of C-SES on Math score for Pupil in public school is larger than the pupil in Catholic school.
1.5.10 m4a
Two level model with two school-level covariate(Sector and Ave_SES) and one pupil-level covariate (C_SES)
# full
summary(m4a <- lmer(Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES +
Sector:C_SES + (C_SES | School), data = dta1))## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +
## (C_SES | School)
## Data: dta1
##
## REML criterion at convergence: 46503.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.15926 -0.72319 0.01704 0.75444 2.95822
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## School (Intercept) 2.380 1.5426
## C_SES 0.101 0.3179 0.39
## Residual 36.721 6.0598
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.1279 0.1993 60.856
## Ave_SES 5.3329 0.3692 14.446
## SectorCatholic 1.2266 0.3063 4.005
## C_SES 2.9450 0.1556 18.928
## Ave_SES:C_SES 1.0393 0.2989 3.477
## SectorCatholic:C_SES -1.6427 0.2398 -6.851
##
## Correlation of Fixed Effects:
## (Intr) Av_SES SctrCt C_SES A_SES:
## Ave_SES 0.256
## SectorCthlc -0.699 -0.356
## C_SES 0.075 0.019 -0.053
## A_SES:C_SES 0.019 0.074 -0.026 0.293
## SctrC:C_SES -0.052 -0.027 0.077 -0.696 -0.351\(Public = 12.1279 + 5.3329*Ave_SES + 2.9450*"C_SES" + 1.0393*"Ave_SES*C_SES"\)
\(Catholic = (12.11+1.2266) + 5.3329*Ave_SES + (2.9450-1.6427)*"C_SES" + 1.0393*"Ave_SES*C_SES"\)
The between-group regression coefficient for meanses is 5.3329.
A pupil attending a higher mean SES school, on average, has a higher math score: an unit increas in school-level meanses leads, on average, to a gain of 5.3329 points in a pupil’s math score.
1.5.11 m4b
# reduced (C_SES | School) into (1 | School)
# 總model去除C_SES 在學校層次的隨機斜率(C_SES | School)，但保留學校層次的隨機截距(1 | School)
summary(m4b <- lmer(Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES +
Sector:C_SES + (1 | School), data = dta1))## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +
## (1 | School)
## Data: dta1
##
## REML criterion at convergence: 46504.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1701 -0.7249 0.0148 0.7542 2.9655
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 2.375 1.541
## Residual 36.766 6.064
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.1282 0.1992 60.885
## Ave_SES 5.3367 0.3690 14.463
## SectorCatholic 1.2245 0.3061 4.000
## C_SES 2.9421 0.1512 19.457
## Ave_SES:C_SES 1.0444 0.2910 3.589
## SectorCatholic:C_SES -1.6422 0.2331 -7.045
##
## Correlation of Fixed Effects:
## (Intr) Av_SES SctrCt C_SES A_SES:
## Ave_SES 0.256
## SectorCthlc -0.699 -0.356
## C_SES 0.000 0.000 0.000
## A_SES:C_SES 0.000 0.000 0.000 0.295
## SctrC:C_SES 0.000 0.000 0.000 -0.696 -0.351# will refit be ML estimation
sjPlot::tab_model(m4a, m4b, show.p=FALSE, show.r2=FALSE, show.obs=FALSE, show.ngroups=FALSE, show.se=TRUE, show.ci=FALSE)| Math | Math | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | Estimates | std. Error |
| (Intercept) | 12.13 | 0.20 | 12.13 | 0.20 |
| Ave_SES | 5.33 | 0.37 | 5.34 | 0.37 |
| Sector [Catholic] | 1.23 | 0.31 | 1.22 | 0.31 |
| C_SES | 2.95 | 0.16 | 2.94 | 0.15 |
| Ave_SES * C_SES | 1.04 | 0.30 | 1.04 | 0.29 |
| Sector [Catholic] * C_SES | -1.64 | 0.24 | -1.64 | 0.23 |
| Random Effects | ||||
| σ2 | 36.72 | 36.77 | ||
| τ00 | 2.38 School | 2.38 School | ||
| τ11 | 0.10 School.C_SES | |||
| ρ01 | 0.39 School | |||
| ICC | 0.06 | 0.06 | ||
anova(m4a, m4b)## refitting model(s) with ML (instead of REML)## Data: dta1
## Models:
## m4b: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +
## m4b: (1 | School)
## m4a: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +
## m4a: (C_SES | School)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m4b 8 46513 46568 -23249 46497
## m4a 10 46516 46585 -23248 46496 1.0016 2 0.6061C_SES 在學校層次的隨機斜率(C_SES | School) 對整體的model影響不大?
# residualual plot
plot(m4b, xlab = "Fitted values", ylab = "Pearson residuals",
pch = 20, cex = .5, type = c("p", "g"))
# quick summary of model parameter estimates "merTools"
library(merTools)## Loading required package: arm## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## arm (Version 1.11-2, built: 2020-7-27)## Working directory is C:/Users/88696/Desktop/R course/20201005## Registered S3 method overwritten by 'broom.mixed':
## method from
## tidy.gamlss broommerTools::fastdisp(m4b)## lmer(formula = Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES +
## Sector:C_SES + (1 | School), data = dta1)
## coef.est coef.se
## (Intercept) 12.13 0.20
## Ave_SES 5.34 0.37
## SectorCatholic 1.22 0.31
## C_SES 2.94 0.15
## Ave_SES:C_SES 1.04 0.29
## SectorCatholic:C_SES -1.64 0.23
##
## Error terms:
## Groups Name Std.Dev.
## School (Intercept) 1.54
## Residual 6.06
## ---
## number of obs: 7185, groups: School, 160
## AIC = 46520.8# estimate fixed-effects parameters by simulation
m4b_fe <- FEsim(m4b, 1000)
# plot for fixed effects
plotFEsim(m4b_fe) +
labs(title = "Coefficient Plot", x = "Median Effect Estimate")
# estimate random-effects parameters by simulation
m4b_re <- REsim(m4b)
# normality plot for random effects
plotREsim(m4b_re)
The largest effect size on the math score is the Ave_SES?
2 Inclass exercise 2
2.1 Info
- Problem: The following R script examines the effect of including predictors at different levels in the school-teacher-pupil example. Compare and explain the variance accounted for in each one of the fitted models.
- Data info: The data came from a study of school differences. The reading attainment score at the start and at the end of year one were recorded for a sample of 777 pupils from 20 London primary schools taught by 46 teachers in UK.
- Understand the meaning of predictors at different levels as well as the effect of the centering the teacher level predictor against the respective school-level means in a school-teacher-pupil three-level example.
Column 1: School ID (1-35)
Column 2: Teacher ID (11-352)
Column 3: Pupil ID (1-1948)
Column 4: Reading attainment at the end of reception (mean = 0, SD = 1)
Column 5: Reading attainment at the end of year one
2.2 Data input
# load package
pacman::p_load(mlmRev, tidyverse, lme4, merTools)
# input data
dta2 <- read.csv("pts.csv")
head(dta2)## School Teacher Pupil Read_0 Read_1
## 1 1 11 18 -0.16407000 4.2426
## 2 1 11 19 -0.98126000 1.0000
## 3 1 11 20 -1.25370000 2.2361
## 4 1 11 15 -0.87230000 3.1623
## 5 1 11 17 -0.00063104 3.4641
## 6 1 11 16 -0.92678000 4.0000# coerce variables to factor type and compute mean Read_0 by school
dta2 <- dta2 %>%
mutate(School = factor(School),
Teacher = factor(Teacher),
Pupil = factor(Pupil)) %>%
group_by( School ) %>%
mutate(msRead_0 = mean(Read_0))
# compute mean Read_0 by teacher and centered teacher mean of Read_0 (from respective school means)
dta2 <- dta2 %>%
group_by( Teacher ) %>%
mutate(mtRead_0 = mean(Read_0),
ctRead_0 = mean(Read_0) - msRead_0 )msRead_0 = mean Read_0 by school.
mtRead_0 = mean Read_0 by teacher.
ctRead_0 = centered teacher mean of Read_0 (from respective school means).
2.3 LMER
2.3.1 m0
Only Two level: The random effects are School and Teacher
# no predictor
summary(m0 <- lmer(Read_1 ~ (1 | School) + (1 | School:Teacher), data = dta2))## Linear mixed model fit by REML ['lmerMod']
## Formula: Read_1 ~ (1 | School) + (1 | School:Teacher)
## Data: dta2
##
## REML criterion at convergence: 2486.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.05438 -0.65686 0.06497 0.62523 2.47233
##
## Random effects:
## Groups Name Variance Std.Dev.
## School:Teacher (Intercept) 0.1370 0.3701
## School (Intercept) 0.1277 0.3574
## Residual 1.3250 1.1511
## Number of obs: 777, groups: School:Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.3755 0.1065 31.69The estimate reading attainment score at the end of year one is 3.3755.
The estimated between teacher in the same school(interaction of teacher and school) variance in reading attainment score is 0.1370.(School-level 的隨機斜率?)
The estimated between school variance in reading attainment score is 0.1277.
The variation between student within the same teacher and in the same school is estimated to be 1.3250.
The between teacher (in school level) variance is about 0.1370/(0.1370 + 0.1277 + 1.3250) = 0.086 (8.6%)of the total variance.
The between school variance is about 0.1277/(0.1370 + 0.1277 + 1.3250) = 0.08 (8%) of the total variance.
The between pupil variance is about 1.3250/(0.1370 + 0.1277 + 1.3250) = 0.83 (83%) of the total variance.
ICC= (0.1370+0.1277)/(0.1370 + 0.1277 + 1.3250) = 0.1665094 (0.17).
2.3.2 m0_s
add a school-level predictor: mean Read_0 by school
# add a school-level predictor
summary(m0_s <- update(m0, . ~ . + msRead_0))## boundary (singular) fit: see ?isSingular## Linear mixed model fit by REML ['lmerMod']
## Formula: Read_1 ~ (1 | School) + (1 | School:Teacher) + msRead_0
## Data: dta2
##
## REML criterion at convergence: 2467.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.06110 -0.65597 0.06822 0.62772 2.45391
##
## Random effects:
## Groups Name Variance Std.Dev.
## School:Teacher (Intercept) 0.1234 0.3513
## School (Intercept) 0.0000 0.0000
## Residual 1.3224 1.1500
## Number of obs: 777, groups: School:Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.39822 0.06808 49.918
## msRead_0 1.00159 0.17807 5.625
##
## Correlation of Fixed Effects:
## (Intr)
## msRead_0 0.086
## convergence code: 0
## boundary (singular) fit: see ?isSingularThe estimate reading attainment score at the end of year one is 3.39822 ≈ 3.3755 (m0).
The estimated between teacher in the same school(interaction of teacher and school) variance in reading attainment score is 0.1234 ≈ 0.1370 (m0).
The estimated between school variance in reading attainment score is 0 < 0.1277(m0). 加了msRead後，不同間School的閱讀分數的變異就不見了，?學校間的變異由每間學校平均的閱讀分數(msRead)決定? (有點廢話) The variation between student within the same teacher and in the same school is estimated to be 1.3224 ≈ 1.3250.
The between teacher (in school level) variance is about 0.1234/(0.1234 + 0 + 1.3224) = 0.08535067 (8.5%) ≈ 0.086 of the total variance.
The between school variance is about 0/(0.1234 + 0 + 1.3224) 0 < 0.08 of the total variance.
The between pupil variance is about 1.3224/(0.1234 + 0 + 1.3224) = 0.9146493 (91%) > 0.83 of the total variance.
2.3.3 m0_t
add a teacher-level predictor: mean Read_0 by teacher
# add a teacher-level predictor
summary(m0_t <- update(m0, . ~ . + mtRead_0))## boundary (singular) fit: see ?isSingular## Linear mixed model fit by REML ['lmerMod']
## Formula: Read_1 ~ (1 | School) + (1 | School:Teacher) + mtRead_0
## Data: dta2
##
## REML criterion at convergence: 2434.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.07539 -0.65979 0.06707 0.63259 2.49973
##
## Random effects:
## Groups Name Variance Std.Dev.
## School:Teacher (Intercept) 3.296e-15 5.741e-08
## School (Intercept) 4.471e-02 2.114e-01
## Residual 1.309e+00 1.144e+00
## Number of obs: 777, groups: School:Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.40649 0.06307 54.010
## mtRead_0 1.07319 0.11495 9.336
##
## Correlation of Fixed Effects:
## (Intr)
## mtRead_0 0.024
## convergence code: 0
## boundary (singular) fit: see ?isSingularThe estimate reading attainment score at the end of year one is 3.40649 ≈ 3.3755 (m0).
The estimated between teacher in the same school(interaction of teacher and school) variance in reading attainment score is 3.296e-15( ≈ 0) < 0.1370 (m0). 不同老師的學生的閱讀分數的變異幾乎為0。
The estimated between school variance in reading attainment score is 0.04471 < 0.1277(m0). 加了mtRead後，不同間School的學生的閱讀分數的變異減少。
The variation between student within the same teacher and in the same school is estimated to be 1.309 ≈ 1.3250.
The between teacher (in school level) variance is about 3.296e-15/1.35371 = 2.43479e-15 < 0.086 of the total variance.
The between school variance is about 4.471e-02/1.35371 = 0.033 < 0.08 of the total variance.
The between pupil variance is about 1.309/1.35371 = 0.9669722 (96%) > 0.83 of the total variance.
2.3.4 m0_ct
add a teacher-level predictor: ctRead_0 = centered teacher mean of Read_0 (from respective school means)
# add a teacher-level predictor away from respective school means
summary(m0_ct <- update(m0, . ~ . + ctRead_0))## boundary (singular) fit: see ?isSingular## Linear mixed model fit by REML ['lmerMod']
## Formula: Read_1 ~ (1 | School) + (1 | School:Teacher) + ctRead_0
## Data: dta2
##
## REML criterion at convergence: 2454.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1804 -0.6461 0.0791 0.6420 2.5272
##
## Random effects:
## Groups Name Variance Std.Dev.
## School:Teacher (Intercept) 2.462e-10 1.569e-05
## School (Intercept) 1.764e-01 4.200e-01
## Residual 1.311e+00 1.145e+00
## Number of obs: 777, groups: School:Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.3896 0.1029 32.932
## ctRead_0 1.1742 0.1581 7.427
##
## Correlation of Fixed Effects:
## (Intr)
## ctRead_0 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingularThe estimate reading attainment score at the end of year one is 3.3896 ≈ 3.3755 (m0).
The estimated between teacher in the same school(interaction of teacher and school) variance in reading attainment score is 2.462e-10( ≈ 0) < 0.1370 (m0).
The estimated between school variance in reading attainment score is 0.1764 > 0.1277(m0).
The variation between student within the same teacher and in the same school is estimated to be 1.311 ≈ 1.3250.
The between teacher (in school level) variance is about 2.462e-10/1.4874 = 1.655237e-10 < 0.086 of the total variance.
The between school variance is about 0.1764/1.4874 = 0.1185962 > 0.08 of the total variance.
The between pupil variance is about 1.311/1.4874 = 0.8814038 (88%) > 0.83 of the total variance.
ICC = (0.1764+2.462e-10)/1.4874 = 0.1185962 (0.12).
2.3.5 Comparison
sjPlot::tab_model(m0, m0_s, m0_t, m0_ct, show.p=FALSE, show.r2=FALSE, show.obs=FALSE, show.ngroups=FALSE, show.se=TRUE, show.ci=FALSE)| Read_1 | Read_1 | Read_1 | Read_1 | |||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Estimates | std. Error | Estimates | std. Error | Estimates | std. Error |
| (Intercept) | 3.38 | 0.11 | 3.40 | 0.07 | 3.41 | 0.06 | 3.39 | 0.10 |
| msRead_0 | 1.00 | 0.18 | ||||||
| mtRead_0 | 1.07 | 0.11 | ||||||
| ctRead_0 | 1.17 | 0.16 | ||||||
| Random Effects | ||||||||
| σ2 | 1.32 | 1.32 | 1.31 | 1.31 | ||||
| τ00 | 0.14 School:Teacher | 0.12 School:Teacher | 0.00 School:Teacher | 0.00 School:Teacher | ||||
| 0.13 School | 0.00 School | 0.04 School | 0.18 School | |||||
| ICC | 0.17 | 0.12 | ||||||
σ2 在這四個model都差不多。
msRead_0 (School level predictor): School level effect turn into 0.
mtRead_0 (Teacher level predictor): Teacher level effect turn into 0, school level effect is reduced.
ctRead_0 (Teacher level predictor): Teacher level effect turn into 0, school level effect is slightly increased, and the ICC slightly reduced.
3 Inclass exercise 3
3.1 Info
- Problem: Replicate the results of analysis reported the classroom data example.
- Data info: The Study of Instructional Improvement was carried out to examine the mathematics achievement scores of students in a random sample of classrooms from the U.S. elementary schools. Data presented here are for 1,190 first-grade students sampled from 312 classrooms in 107 schools. The response measure is the change in student’s mathematics achievement scores.
Column 1: Gender ID (B = Boy, G = Girl)
Column 2: Minority ID (M = Minority, N = Non-minority)
Column 3: Student’s math score in the spring of kindergarten
Column 4: Student’s gain in math score from the spring of kindergarten to the spring of first grade
Column 5: Student’s socio-economic status
Column 6: First-grade teacher’s teaching experience
Column 7: First-grade teacher’s mathematical content knowledge
Column 8: Percentage of households in the neighborhood of the school below poverty level
Column 9: First-grade teacher’s mathematics preparation
Column 10: Class ID
Column 11: School ID
Column 12: Student ID
3.2 Data input
# load package
pacman::p_load(mlmRev, tidyverse, lme4, merTools, sjPlot, sjmisc, sjlabelled)
data(classroom, package="WWGbook")
dta3 <- classroom[,c(11,10,12,4,3,1,2,5,6)]
names(dta3) <- c("School", "Class", "Pupil", "math1", "math0",
"Sex", "Minority", "SES", "Teacher_experience")
dta3 <- dta3 %>%
mutate(School = factor(School),
Class = factor(Class),
Pupil = factor(Pupil),
Sex = factor(Sex)) %>%
mutate(math0_a = mean(math0)) %>%
group_by(School) %>%
mutate(math0_sa = mean(math0),
math0_c = math0 - math0_sa)
head(dta3)## # A tibble: 6 x 12
## # Groups: School [1]
## School Class Pupil math1 math0 Sex Minority SES Teacher_experie~ math0_a
## <fct> <fct> <fct> <int> <int> <fct> <int> <dbl> <dbl> <dbl>
## 1 1 160 1 32 448 1 1 0.46 1 467.
## 2 1 160 2 109 460 0 1 -0.27 1 467.
## 3 1 160 3 56 511 1 1 -0.03 1 467.
## 4 1 217 4 83 449 0 1 -0.38 2 467.
## 5 1 217 5 53 425 0 1 -0.03 2 467.
## 6 1 217 6 65 450 1 1 0.76 2 467.
## # ... with 2 more variables: math0_sa <dbl>, math0_c <dbl>3.2.1 m0
only consider the School and Class level
summary(m0 <- lmer(math1 ~ (1 | School) + (1 | Class), data=dta3)) ## Linear mixed model fit by REML ['lmerMod']
## Formula: math1 ~ (1 | School) + (1 | Class)
## Data: dta3
##
## REML criterion at convergence: 11768.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6441 -0.5984 -0.0336 0.5334 5.6335
##
## Random effects:
## Groups Name Variance Std.Dev.
## Class (Intercept) 99.22 9.961
## School (Intercept) 77.50 8.804
## Residual 1028.23 32.066
## Number of obs: 1190, groups: Class, 312; School, 107
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 57.427 1.443 39.79tab_model(m0,
show.intercept = TRUE,
show.est = TRUE,
show.se = TRUE,
show.df = TRUE,
show.stat = TRUE,
show.p = TRUE,
show.aic = TRUE,
show.aicc = TRUE)| math1 | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 57.43 | 1.44 | 54.60 – 60.26 | 39.79 | <0.001 | 1186.00 |
| Random Effects | ||||||
| σ2 | 1028.23 | |||||
| τ00 Class | 99.22 | |||||
| τ00 School | 77.50 | |||||
| ICC | 0.15 | |||||
| N School | 107 | |||||
| N Class | 312 | |||||
| Observations | 1190 | |||||
| Marginal R2 / Conditional R2 | 0.000 / 0.147 | |||||
| AIC | 11776.799 | |||||
3.2.2 m1
fixed effect: math0, Sex, Minority, SES random effect: School, Class math1ijk(j) = b0j + b1k(j) + β2 × math0ijk(j) + β3 × sexijk(j) + β4 × minorityijk(j) + β5 × sesijk(j) + εijk(j),
where i is student index, j is school index and k is class index.
summary(m1 <- lmer(math1 ~ math0 + Sex + Minority + SES + (1 | School) + (1 |School:Class), data=dta3)) ## Linear mixed model fit by REML ['lmerMod']
## Formula:
## math1 ~ math0 + Sex + Minority + SES + (1 | School) + (1 | School:Class)
## Data: dta3
##
## REML criterion at convergence: 11385.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8257 -0.6110 -0.0337 0.5538 4.2678
##
## Random effects:
## Groups Name Variance Std.Dev.
## School:Class (Intercept) 83.28 9.126
## School (Intercept) 75.20 8.672
## Residual 734.57 27.103
## Number of obs: 1190, groups: School:Class, 312; School, 107
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 282.79033 10.85323 26.056
## math0 -0.46980 0.02227 -21.100
## Sex1 -1.25119 1.65773 -0.755
## Minority -8.26213 2.34011 -3.531
## SES 5.34638 1.24109 4.308
##
## Correlation of Fixed Effects:
## (Intr) math0 Sex1 Minrty
## math0 -0.978
## Sex1 -0.044 -0.031
## Minority -0.307 0.163 -0.018
## SES 0.140 -0.168 0.019 0.163tab_model(m1,
show.intercept = TRUE,
show.est = TRUE,
show.se = TRUE,
show.df = TRUE,
show.stat = TRUE,
show.p = TRUE,
show.aic = TRUE,
show.aicc = TRUE)| math1 | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 282.79 | 10.85 | 261.52 – 304.06 | 26.06 | <0.001 | 1182.00 |
| math0 | -0.47 | 0.02 | -0.51 – -0.43 | -21.10 | <0.001 | 1182.00 |
| Sex [1] | -1.25 | 1.66 | -4.50 – 2.00 | -0.75 | 0.450 | 1182.00 |
| Minority | -8.26 | 2.34 | -12.85 – -3.68 | -3.53 | <0.001 | 1182.00 |
| SES | 5.35 | 1.24 | 2.91 – 7.78 | 4.31 | <0.001 | 1182.00 |
| Random Effects | ||||||
| σ2 | 734.57 | |||||
| τ00 School:Class | 83.28 | |||||
| τ00 School | 75.20 | |||||
| ICC | 0.18 | |||||
| N School | 107 | |||||
| N Class | 312 | |||||
| Observations | 1190 | |||||
| Marginal R2 / Conditional R2 | 0.274 / 0.403 | |||||
| AIC | 11401.925 | |||||
sjPlot::tab_model(m0, m1, show.p=FALSE, show.r2=FALSE, show.obs=FALSE, show.ngroups=FALSE, show.se=TRUE, show.ci=FALSE)| math1 | math1 | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | Estimates | std. Error |
| (Intercept) | 57.43 | 1.44 | 282.79 | 10.85 |
| math0 | -0.47 | 0.02 | ||
| Sex [1] | -1.25 | 1.66 | ||
| Minority | -8.26 | 2.34 | ||
| SES | 5.35 | 1.24 | ||
| Random Effects | ||||
| σ2 | 1028.23 | 734.57 | ||
| τ00 | 99.22 Class | 83.28 School:Class | ||
| 77.50 School | 75.20 School | |||
| ICC | 0.15 | 0.18 | ||