Part 1: Random intercepts

Q1

Create a table of summary statistics for the level 1 variables.

sumtable(data = hsb, vars = c("minority", "female", "ses", "mathach"))

Summary Statistics
Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 75	Max
minority	7185	0.275	0.446	0	0	1	1
female	7185	0.528	0.499	0	0	1	1
ses	7185	0	0.779	-3.758	-0.538	0.602	2.692
mathach	7185	12.748	6.878	-2.832	7.275	18.317	24.993

1-2 sentences about descriptive stats

Q2

Create a table of summary statistics for the level 2 variables.

hsb$schid <- as.factor(hsb$schid)
sumtable(data = hsb, vars = c("schid", "meanses", "size", "sector"))

Summary Statistics
Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 75	Max
schid	7185
… 1224	47	0.7%
… 1288	25	0.3%
… 1296	48	0.7%
… 1308	20	0.3%
… 1317	48	0.7%
… 1358	30	0.4%
… 1374	28	0.4%
… 1433	35	0.5%
… 1436	44	0.6%
… 1461	33	0.5%
… 1462	57	0.8%
… 1477	62	0.9%
… 1499	53	0.7%
… 1637	27	0.4%
… 1906	53	0.7%
… 1909	28	0.4%
… 1942	29	0.4%
… 1946	39	0.5%
… 2030	47	0.7%
… 2208	60	0.8%
… 2277	61	0.8%
… 2305	67	0.9%
… 2336	47	0.7%
… 2458	57	0.8%
… 2467	52	0.7%
… 2526	57	0.8%
… 2626	38	0.5%
… 2629	57	0.8%
… 2639	42	0.6%
… 2651	38	0.5%
… 2655	52	0.7%
… 2658	45	0.6%
… 2755	47	0.7%
… 2768	25	0.3%
… 2771	55	0.8%
… 2818	42	0.6%
… 2917	43	0.6%
… 2990	48	0.7%
… 2995	46	0.6%
… 3013	53	0.7%
… 3020	59	0.8%
… 3039	21	0.3%
… 3088	39	0.5%
… 3152	52	0.7%
… 3332	38	0.5%
… 3351	39	0.5%
… 3377	45	0.6%
… 3427	49	0.7%
… 3498	53	0.7%
… 3499	38	0.5%
… 3533	48	0.7%
… 3610	64	0.9%
… 3657	51	0.7%
… 3688	43	0.6%
… 3705	45	0.6%
… 3716	41	0.6%
… 3838	54	0.8%
… 3881	41	0.6%
… 3967	52	0.7%
… 3992	53	0.7%
… 3999	46	0.6%
… 4042	64	0.9%
… 4173	44	0.6%
… 4223	45	0.6%
… 4253	58	0.8%
… 4292	65	0.9%
… 4325	53	0.7%
… 4350	33	0.5%
… 4383	25	0.3%
… 4410	41	0.6%
… 4420	32	0.4%
… 4458	48	0.7%
… 4511	58	0.8%
… 4523	47	0.7%
… 4530	63	0.9%
… 4642	61	0.8%
… 4868	34	0.5%
… 4931	58	0.8%
… 5192	28	0.4%
… 5404	57	0.8%
… 5619	66	0.9%
… 5640	57	0.8%
… 5650	45	0.6%
… 5667	61	0.8%
… 5720	53	0.7%
… 5761	52	0.7%
… 5762	37	0.5%
… 5783	29	0.4%
… 5815	25	0.3%
… 5819	50	0.7%
… 5838	31	0.4%
… 5937	29	0.4%
… 6074	56	0.8%
… 6089	33	0.5%
… 6144	43	0.6%
… 6170	21	0.3%
… 6291	35	0.5%
… 6366	58	0.8%
… 6397	60	0.8%
… 6415	54	0.8%
… 6443	30	0.4%
… 6464	29	0.4%
… 6469	57	0.8%
… 6484	35	0.5%
… 6578	56	0.8%
… 6600	56	0.8%
… 6808	44	0.6%
… 6816	55	0.8%
… 6897	49	0.7%
… 6990	53	0.7%
… 7011	33	0.5%
… 7101	28	0.4%
… 7172	44	0.6%
… 7232	52	0.7%
… 7276	53	0.7%
… 7332	48	0.7%
… 7341	51	0.7%
… 7342	58	0.8%
… 7345	56	0.8%
… 7364	44	0.6%
… 7635	51	0.7%
… 7688	54	0.8%
… 7697	32	0.4%
… 7734	22	0.3%
… 7890	51	0.7%
… 7919	37	0.5%
… 8009	47	0.7%
… 8150	44	0.6%
… 8165	49	0.7%
… 8175	33	0.5%
… 8188	30	0.4%
… 8193	43	0.6%
… 8202	35	0.5%
… 8357	27	0.4%
… 8367	14	0.2%
… 8477	37	0.5%
… 8531	41	0.6%
… 8627	53	0.7%
… 8628	61	0.8%
… 8707	48	0.7%
… 8775	48	0.7%
… 8800	32	0.4%
… 8854	32	0.4%
… 8857	64	0.9%
… 8874	36	0.5%
… 8946	58	0.8%
… 8983	51	0.7%
… 9021	56	0.8%
… 9104	55	0.8%
… 9158	53	0.7%
… 9198	31	0.4%
… 9225	36	0.5%
… 9292	19	0.3%
… 9340	29	0.4%
… 9347	57	0.8%
… 9359	53	0.7%
… 9397	47	0.7%
… 9508	35	0.5%
… 9550	29	0.4%
… 9586	59	0.8%
meanses	7185	0.006	0.414	-1.188	-0.317	0.333	0.831
size	7185	1056.862	604.172	100	565	1436	2713
sector	7185	0.493	0.5	0	0	1	1

1-2 sentences about descriptive stats

Q3

Estimate an “empty” random intercept model that allows you to estimate the ICC of mathematics achievement across schools. Report and interpret the ICC.

m.empty <- lmer(mathach ~ 1 + (1 | schid), data = hsb, REML = F)
summary(m.empty)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mathach ~ 1 + (1 | schid)
##    Data: hsb
## 
##      AIC      BIC   logLik deviance df.resid 
##  47121.8  47142.4 -23557.9  47115.8     7182 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.06262 -0.75365  0.02676  0.76070  2.74184 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schid    (Intercept)  8.553   2.925   
##  Residual             39.148   6.257   
## Number of obs: 7185, groups:  schid, 160
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  12.6371     0.2436 157.6209   51.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

prop_between <- 8.553
prop_within <- 39.148

(ICC <- (prop_between / (prop_within + prop_between)))

## [1] 0.1793044

ANSWER: the icc describes how much of the variance in the outcome, math achievement, is between vs. within groups. In this case the ICC = .179, meaning some of our observations are from the same clusters and we should calculate an effective sample size (N_effective = 244).

n_j <- 160 #160 schools
(Deff <- 1 + (n_j - 1) * ICC) #29.87

## [1] 29.5094

N <- 7185
(N_effect <- N /Deff)

## [1] 243.4817

Q4

Estimate a regular OLS regression model predicting math achievement from students’ family SES and then estimate a random-intercept model predicting math achievement from students’ family SES. Create a regression table to compare the models and briefly comment on what you notice about similarities and differences between the two models.

m.ols <- lm(mathach ~ ses, data = hsb)
m.mlm <- lmer(mathach ~ ses + (1 | schid), data = hsb, REML = F)

tab_model(m.ols, m.mlm)

	mathach			mathach
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	12.75	12.60 – 12.90	<0.001	12.66	12.29 – 13.02	<0.001
ses	3.18	2.99 – 3.37	<0.001	2.39	2.18 – 2.60	<0.001
Random Effects
σ²				37.03
τ₀₀				4.73 _schid
ICC				0.11
N				160 _schid
Observations	7185			7185
R² / R² adjusted	0.130 / 0.130			0.077 / 0.181

ANSWER: Both models include the same number of observations and generally the same conclusion–that family SES significantly predicts math achievement. Yet, beyond these two points, most everything else is different, if only slightly. The intercept in the MLM model slightly decreases from the OLS while the slope increases. The MLM model also includes random effects which account for a portion of the error variance not accounted for in the OLS model. This includes the estimated variance of group means for schools (4.37), estimated variance of within group residuals (37.03).

Q5

Estimate a random intercept model predicting math achievement from SES that allows you to estimate both the within-school and between-school association between SES and math achievement. Interpret the within-school coefficient and the between-school coefficient. Indicate whether these associations are statistically significantly different. What can you conclude from this analysis? Your answer should also address the idea of the “contextual” effect of SES on math achievement.

## meanses: within school mean SES (to create ses.cwc)
## want level 1 SES, and school level as predictor
## BIG Q: how do you center level 1 SES?
hsb$ses.cwc <- hsb$ses - hsb$meanses # centering within cluster
hsb$ses.gmc <- hsb$ses - mean(hsb$ses) # grand mean centering

m.a2 <- lmer(mathach ~ 1 + ses.cwc + (1 | schid), data = hsb, REML = F)
m.a3 <- lmer(mathach ~ 1 + ses.cwc + meanses + (1 | schid), data = hsb, REML = F)

tab_model(m.a2, m.a3,
          string.est = "Coef (s.e.)",
          show.ci = F,
          show.p = F, 
          linebreak = F,
          collapse.se = T, 
          show.dev = T)

	mathach	mathach
Predictors	Coef (s.e.)	Coef (s.e.)
(Intercept)	12.65 (0.24)	12.66 (0.15)
ses.cwc	2.19 (0.11)	2.19 (0.11)
meanses		5.87 (0.36)
Random Effects
σ²	37.01	37.01
τ₀₀	8.61 _schid	2.65 _schid
ICC	0.19	0.07
N	160 _schid	160 _schid
Observations	7185	7185
Marginal R² / Conditional R²	0.044 / 0.224	0.167 / 0.223
Deviance	46720.413	46563.805

i = student, s = school

\(ses.cwc\) = \(ses_{is}\) - \(mean(ses_{.s})\)

\(mathAch_{is}\) = \(\gamma_{00}\) + \(\gamma_{10}\)(\(ses.cwc\)) + \(\gamma_{01}\)(\(mean(ses_{.s})\)) + \(u_{0s}\) + \(u_{1s}\)(\(ses.cwc\)) + \(\epsilon_{is}\)

The within coefficient (\(\gamma_{10}\) = \(\beta^w\) = 2.19) represents the different in math achievement scores between two students in the same school who differ on 1 unit on SES.

Between coefficient (\(\gamma_{01}\) = \(\beta^b\) = 5.87) represents the expected different between the means of two schools which differ by 1 unit in average SES.

The “contextual” effect is the expected difference in the math achievement between two students who have the same individual SES, but who attend schools differing by one unit in the mean SES. This indicates the increment by learning that accrues by virtue of being educated in one schoool vs. another that are 1 unit separate from each other.

\(\beta^c\) = \(\beta^b\) -\(\beta^w\) = 3.68

to test if there is a significant difference between \(\beta^b\) and \(\beta^w\), I need to run a GMC model to test \(H_0\): \(\gamma_{01}\) = 0

m.a3 <- lmer(mathach ~ 1 + ses.gmc + meanses + (1 | schid), data = hsb, REML = F)
summary(m.a3)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mathach ~ 1 + ses.gmc + meanses + (1 | schid)
##    Data: hsb
## 
##      AIC      BIC   logLik deviance df.resid 
##  46573.8  46608.2 -23281.9  46563.8     7180 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1675 -0.7257  0.0176  0.7553  2.9445 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schid    (Intercept)  2.647   1.627   
##  Residual             37.014   6.084   
## Number of obs: 7185, groups:  schid, 160
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   12.6619     0.1484  155.6420  85.317   <2e-16 ***
## ses.gmc        2.1912     0.1087 7022.4974  20.165   <2e-16 ***
## meanses        3.6745     0.3754  184.9826   9.788   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) ss.gmc
## ses.gmc  0.004       
## meanses -0.005 -0.289

We find that the slope is significantly different from zero, so \(\beta^b\) is significantly larger than \(\beta^w\) by 3.68 units.

Part 2: Random coefficients

Question 6

Estimate conditional models (model with x variables) with the following variables:

#school centering for minority (because already have this for SES)
hsb <- hsb%>%
  group_by(schid)%>% #grouping by school ID
  mutate(
    meanminority = mean(minority), #get proportion minority in each school
    minority.cwc = minority - meanminority #center within school
  )%>%
  ungroup()

#Model 1: fit a random-int model (all slopes fixed)
m.6.1 <- lmer(mathach ~ 1 + ses.cwc + minority.cwc + meanses + sector +
                (1 | schid), data = hsb, REML = F)

#Model 2: fit a random-coef model allowing both level-1 slopes to vary randomly.
m.6.2 <- lmer(mathach ~ 1 + ses.cwc + minority.cwc + meanses + sector +
                (minority.cwc + ses.cwc | schid), data = hsb, REML = F)

## boundary (singular) fit: see ?isSingular

tab_model(m.6.1, m.6.2,
          string.est = "Coef (s.e.)",
          show.ci = F,
          show.p = T, 
          linebreak = F,
          collapse.se = T, 
          show.dev = T)

	mathach		mathach
Predictors	Coef (s.e.)	p	Coef (s.e.)	p
(Intercept)	12.11 (0.20)	<0.001	12.04 (0.20)	<0.001
ses.cwc	1.95 (0.11)	<0.001	1.94 (0.12)	<0.001
minority.cwc	-2.90 (0.22)	<0.001	-2.95 (0.26)	<0.001
meanses	5.34 (0.37)	<0.001	5.22 (0.36)	<0.001
sector	1.22 (0.30)	<0.001	1.39 (0.30)	<0.001
Random Effects
σ²	36.13		35.66
τ₀₀	2.33 _schid		2.36 _schid
τ₁₁			1.97 _{schid.minority.cwc}
			0.46 _{schid.ses.cwc}
ρ₀₁			-0.05
			0.27
ICC	0.06		0.07
N	160 _schid		160 _schid
Observations	7185		7185
Marginal R² / Conditional R²	0.193 / 0.242		0.193 / 0.252
Deviance	46377.412		46356.097

anova(m.6.1, m.6.2)

Focusing on Model 2, interpret the estimated regression coefficients. Is there evidence of variation in the following quantities across schools: 1) racial inequality as measured by the difference between scores for minority and non-minority students 2) socioeconomic inequality as measured by differences in scores across SES levels?

Yes, there is evidence for variation in math achievement scores for centered within cluster minority vs. non-minority, B = -2.95, SE = 0.12, p < .001, such that as the proportion of minority students moves from 0 to 1 there is a 2.95 units score decrease in math achievement.

Yes, there’s evidence for variation in math achievement scores for centered within cluster SES across schools, B = 1.94, SE = 0.12, p < .001. This means that as you move 1 unit up in school SES, there is a 1.94 unit increase for average math achievement scores.

Compared to the model that solely includes random intercepts, the model that includes random slopes significantly improves predictive power. Comparing the two models, there is a significant difference in the models, p < .001, AIC decreases by 11 units and deviance decreases by 21.32 units when moving from the intercept only model to the slope model.

Part 3: Interactions

Question 7

Now extend Model 2 from question 6 so that you can answer the question: are there differences in either the racial/ethnic inequality or SES differentiation across school sectors? Estimate an appropriate MLM. Interpret the fixed effects estimates and appropriate variance components. Summarize your conclusions from the results.

m.7.1 <- lmer(mathach ~ 1 + (ses.cwc + minority.cwc) * sector + meanses +
                (minority.cwc + ses.cwc | schid), data = hsb, REML = F)

tab_model(
  m.6.2, m.7.1,
  string.est = "Coef (s.e.)",
  show.p = T,
  show.ci = F,
  linebreak = F,
  collapse.se = T,
  show.aic = T,
  show.dev = T
  )

	mathach		mathach
Predictors	Coef (s.e.)	p	Coef (s.e.)	p
(Intercept)	12.04 (0.20)	<0.001	12.09 (0.20)	<0.001
ses.cwc	1.94 (0.12)	<0.001	2.37 (0.15)	<0.001
minority.cwc	-2.95 (0.26)	<0.001	-3.92 (0.34)	<0.001
meanses	5.22 (0.36)	<0.001	5.20 (0.36)	<0.001
sector	1.39 (0.30)	<0.001	1.26 (0.30)	<0.001
ses.cwc * sector			-1.03 (0.23)	<0.001
minority.cwc * sector			2.05 (0.48)	<0.001
Random Effects
σ²	35.66		35.67
τ₀₀	2.36 _schid		2.35 _schid
τ₁₁	1.97 _{schid.minority.cwc}		0.85 _{schid.minority.cwc}
	0.46 _{schid.ses.cwc}		0.20 _{schid.ses.cwc}
ρ₀₁	-0.05		-0.09
	0.27		0.40
ICC	0.07		0.07
N	160 _schid		160 _schid
Observations	7185		7185
Marginal R² / Conditional R²	0.193 / 0.252		0.195 / 0.248
Deviance	46356.097		46320.754
AIC	46380.097		46348.754

anova(m.6.2, m.7.1)

plot_model(m.7.1, type = "pred", terms = c("minority.cwc", "sector")) + theme_bw()

plot_model(m.7.1, type = "pred", terms = c("ses.cwc", "sector")) + theme_bw()

Yes, there is an interaction between minority vs. non-minority status by school sector, B = 2.05, SE = 0.48, p < .001. Specifically, as you move from a Catholic school to a public school, the strength of the relationship between minority status and math achievement increases by 2.05 units. There is also a significant interaction between SES and sector, B = -1.03, SE = 0.23, p < .001, such that SES has a 1.03 greater impact on match achievement in public schools compared with catholic schools.

interpret variance components

\(\sigma^2\)

\(\tau_{00}\)

\(\tau_{11}\)

Question 8

Calculate and interpret a measure of R^2 that quantifies the proportion of variation in level-1 student math scores explained by the model.

tab_model(
  m.empty, m.7.1,
  string.est = "Coef (s.e.)",
  show.p = T,
  show.ci = F,
  linebreak = F,
  collapse.se = T
  )

	mathach		mathach
Predictors	Coef (s.e.)	p	Coef (s.e.)	p
(Intercept)	12.64 (0.24)	<0.001	12.09 (0.20)	<0.001
ses.cwc			2.37 (0.15)	<0.001
minority.cwc			-3.92 (0.34)	<0.001
sector			1.26 (0.30)	<0.001
meanses			5.20 (0.36)	<0.001
ses.cwc * sector			-1.03 (0.23)	<0.001
minority.cwc * sector			2.05 (0.48)	<0.001
Random Effects
σ²	39.15		35.67
τ₀₀	8.55 _schid		2.35 _schid
τ₁₁			0.85 _{schid.minority.cwc}
			0.20 _{schid.ses.cwc}
ρ₀₁			-0.09
			0.40
ICC	0.18		0.07
N	160 _schid		160 _schid
Observations	7185		7185
Marginal R² / Conditional R²	0.000 / 0.179		0.195 / 0.248

sigma_sq_8 <- 39.15
sigma_sq_7 <- 35.67
(R_sq <- (sigma_sq_8 - sigma_sq_7)/sigma_sq_8)

## [1] 0.08888889

This means that our level 1 variables in this model explain around 8.9% of the variation in student math scores.

Question 9

Calculate and interpret a measure of R^2 that quantifies the proportion of variation in school-level SES differentiation that is explained by school sector.

tab_model(
  m.6.2, m.7.1, 
  string.est = "Coef (s.e.)",
  show.p = T,
  show.ci = F,
  linebreak = F,
  collapse.se = T
  )

	mathach		mathach
Predictors	Coef (s.e.)	p	Coef (s.e.)	p
(Intercept)	12.04 (0.20)	<0.001	12.09 (0.20)	<0.001
ses.cwc	1.94 (0.12)	<0.001	2.37 (0.15)	<0.001
minority.cwc	-2.95 (0.26)	<0.001	-3.92 (0.34)	<0.001
meanses	5.22 (0.36)	<0.001	5.20 (0.36)	<0.001
sector	1.39 (0.30)	<0.001	1.26 (0.30)	<0.001
ses.cwc * sector			-1.03 (0.23)	<0.001
minority.cwc * sector			2.05 (0.48)	<0.001
Random Effects
σ²	35.66		35.67
τ₀₀	2.36 _schid		2.35 _schid
τ₁₁	1.97 _{schid.minority.cwc}		0.85 _{schid.minority.cwc}
	0.46 _{schid.ses.cwc}		0.20 _{schid.ses.cwc}
ρ₀₁	-0.05		-0.09
	0.27		0.40
ICC	0.07		0.07
N	160 _schid		160 _schid
Observations	7185		7185
Marginal R² / Conditional R²	0.193 / 0.252		0.195 / 0.248

t_11_6 <- .46
t_11_7 <- .20

((t_11_6 - t_11_7)/t_11_6)

## [1] 0.5652174

Around 56.5% of the variation in SES is explained by school sector.

Question 10

Create and interpret:

a normal QQ plot for the level-1 residuals
a normal QQ plot for one set of level-2 residuals.

Explain why these plots are relevant to the use of the MLM in this section.

QQ plots help us determine if our model violates assumptions about error variance being distributed normally. First, our level 1 residual error seems to be more or less distributed normally, but could be slightly light tailed. Second, the level 2 residuals look similar but there are fewer data points near the higher and lower limits of the data set which makes it more difficult to conclude there is a light tail for error variance–rather it could just be some outliers.

resid.L1 <- resid(m.7.1)
resid.L2 <- ranef(m.7.1)
resid.L2 <- (resid.L2$schid)

qqnorm(resid.L1)
qqline(resid.L1, col = "purple")

qqnorm(resid.L2$`(Intercept)`)
qqline(resid.L2$`(Intercept)`, col = "red")

qqnorm(resid.L2$minority.cwc)
qqline(resid.L2$minority.cwc, col = "green")

qqnorm(resid.L2$ses.cwc)
qqline(resid.L2$ses.cwc, col = "blue")

PSET B

Dani Grant

2022-10-14