Advanced Two Level Hierarchical Models
Nirmal Ghimire, PhD
2022-05-02
This is a sample hierarchical linear model research. This model begins with the necessary nitty-gritty of data analyses in R and culminates in an advanced 2-level mixed linear model (MLM) which contains fixed slopes, random slopes and even moderators. It took quite a while for me to figure out ways to conduct these analyses, especially the MLMs. Thus, this report is for my own resource. Thus, it doesn’t have much description. I have some explanation especially of the one’s I think I may forget. If you happen to visit this page, you are on luck!!
Highlights
- Null Model Using
lmer()function - Random Intercept, Fixed Slope Models
- Random Intercept, Random Slope Models
- Best Fitting Model
- Model Comparison
- Plotting Models (Random effects and Fixed effects)
- APA Tables and Figures
- Calculating R-Squared, Omega-Squared values etc.
Select Working Directory
setwd("C:/Users/nirma/Desktop/FJER Copies")Loading Required Libraries
When you become an advance R user, you refrain loading all the
libraries. We only load the libraries that we use all the times like
{tidyverse} or {ggplot}, and we call other libraries on the go as we
need using {library_name::function()}. R stores all of
these objects in RAM (Random Access Memory), thus, loading all the
libraries may seriously compromise the speed of our machine.
library(tidyverse)
library(ggplot2)
library(sjlabelled)
library(sjPlot)
library(stringr)
library(r2mlm)
library(lme4)
library(readr)
library(here)
library(papaja)
library(geomtextpath)
library(janitor)
library(lmerTest)
library(sjstats)
# library(QuantPsyc)Uploading My Data
final_data <- read_csv("EDF7006_final_Study_TWS_Data_Spring_2018_1.csv")
# Checking the dimension of the data set
dim(final_data)[1] 5561 25# Checking the summary of the data
summary(final_data) T_ID STUDENT TPROGRAM Majors TGRADE
Min. : 1 Min. : 1 Min. :0.000 Min. :0.000 Min. : 0.00
1st Qu.: 68 1st Qu.:1391 1st Qu.:0.000 1st Qu.:0.000 1st Qu.: 2.00
Median :132 Median :2781 Median :0.000 Median :0.000 Median : 4.00
Mean :132 Mean :2781 Mean :0.308 Mean :0.028 Mean : 4.29
3rd Qu.:199 3rd Qu.:4171 3rd Qu.:0.000 3rd Qu.:0.000 3rd Qu.: 5.00
Max. :236 Max. :5561 Max. :2.000 Max. :1.000 Max. :12.00
SGrd_Lvl MALE ETHNICITY FRLUNCH EL
Min. :0.000 Min. :0.000 Min. :0.00 Min. :0.00 Min. :0.000
1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.00 1st Qu.:0.00 1st Qu.:0.000
Median :0.000 Median :0.000 Median :1.00 Median :1.00 Median :0.000
Mean :0.307 Mean :0.489 Mean :1.11 Mean :0.51 Mean :0.105
3rd Qu.:0.000 3rd Qu.:1.000 3rd Qu.:1.00 3rd Qu.:1.00 3rd Qu.:0.000
Max. :2.000 Max. :1.000 Max. :5.00 Max. :1.00 Max. :1.000
NA's :1 NA's :1 NA's :1
PRELG1 PRELG2 PRELG3 POSTLG1 POSTLG2
Min. : 0 Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.0
1st Qu.: 1 1st Qu.: 1.0 1st Qu.: 1.0 1st Qu.: 3.0 1st Qu.: 3.0
Median : 3 Median : 2.0 Median : 2.0 Median : 5.0 Median : 4.0
Mean : 7 Mean : 5.1 Mean : 5.1 Mean : 11.7 Mean : 9.4
3rd Qu.: 6 3rd Qu.: 5.0 3rd Qu.: 5.0 3rd Qu.: 10.0 3rd Qu.: 10.0
Max. :100 Max. :100.0 Max. :100.0 Max. :100.0 Max. :100.0
NA's :58 NA's :124 NA's :121 NA's :55 NA's :109
POSTLG3 PREPERCENT POSTPERCENT PPG
Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. :-45.8
1st Qu.: 2.0 1st Qu.: 26.7 1st Qu.: 73.8 1st Qu.: 20.0
Median : 5.0 Median : 43.8 Median : 86.7 Median : 36.8
Mean : 9.6 Mean : 45.2 Mean : 82.5 Mean : 37.3
3rd Qu.: 9.0 3rd Qu.: 62.5 3rd Qu.: 95.8 3rd Qu.: 53.4
Max. :100.0 Max. :100.0 Max. :100.0 Max. :100.0
NA's :96 NA's :90 NA's :90 NA's :90
White Elementary Math Taught_ELEM
Min. :0.000 Min. :-1.000 Min. :-1.000 Min. :-1.000
1st Qu.:0.000 1st Qu.: 1.000 1st Qu.: 0.000 1st Qu.: 0.000
Median :0.000 Median : 1.000 Median : 0.000 Median : 0.000
Mean :0.463 Mean : 0.692 Mean :-0.113 Mean :-0.062
3rd Qu.:1.000 3rd Qu.: 1.000 3rd Qu.: 0.000 3rd Qu.: 0.000
Max. :1.000 Max. : 1.000 Max. : 1.000 Max. : 1.000
Taught_Middle filter_$
Min. :-1.00 Min. :0.000
1st Qu.: 0.00 1st Qu.:0.000
Median : 0.00 Median :0.000
Mean :-0.04 Mean :0.406
3rd Qu.: 0.00 3rd Qu.:1.000
Max. : 1.00 Max. :1.000
Checking the Column Names
names(final_data) [1] "T_ID" "STUDENT" "TPROGRAM" "Majors"
[5] "TGRADE" "SGrd_Lvl" "MALE" "ETHNICITY"
[9] "FRLUNCH" "EL" "PRELG1" "PRELG2"
[13] "PRELG3" "POSTLG1" "POSTLG2" "POSTLG3"
[17] "PREPERCENT" "POSTPERCENT" "PPG" "White"
[21] "Elementary" "Math" "Taught_ELEM" "Taught_Middle"
[25] "filter_$" Subsetting the data by selecting the required columns
# Selecting only the Required Columns
final_data <- final_data |>
select(
T_ID, STUDENT, TPROGRAM,
SGrd_Lvl, MALE, ETHNICITY, FRLUNCH,
EL, PREPERCENT, POSTPERCENT, White,
Elementary, Math, Taught_ELEM, Taught_Middle
)
# Renaming the desired Column
final_data <- final_data |>
rename(TGRADE = SGrd_Lvl)
# Checking the Structure of the Entire Dataset
str(final_data)tibble [5,561 x 15] (S3: tbl_df/tbl/data.frame)
$ T_ID : num [1:5561] 1 1 1 1 1 1 1 1 1 1 ...
$ STUDENT : num [1:5561] 1 2 3 4 5 6 7 8 9 10 ...
$ TPROGRAM : num [1:5561] 0 0 0 0 0 0 0 0 0 0 ...
$ TGRADE : num [1:5561] 0 0 0 0 0 0 0 0 0 0 ...
$ MALE : num [1:5561] 0 1 0 1 1 1 0 0 1 0 ...
$ ETHNICITY : num [1:5561] 1 0 0 3 3 0 1 1 0 1 ...
$ FRLUNCH : num [1:5561] 1 1 0 1 1 1 1 1 1 1 ...
$ EL : num [1:5561] 1 0 1 0 1 0 1 1 0 1 ...
$ PREPERCENT : num [1:5561] 50 80 70 60 70 NA 50 40 70 60 ...
$ POSTPERCENT : num [1:5561] 90 80 70 80 80 NA 70 60 70 60 ...
$ White : num [1:5561] 0 1 1 0 0 1 0 0 1 0 ...
$ Elementary : num [1:5561] 1 1 1 1 1 1 1 1 1 1 ...
$ Math : num [1:5561] 0 0 0 0 0 0 0 0 0 0 ...
$ Taught_ELEM : num [1:5561] 0 0 0 0 0 0 0 0 0 0 ...
$ Taught_Middle: num [1:5561] 0 0 0 0 0 0 0 0 0 0 ...Cleaning the Column Names for Consistency
# Using clean columns names for consistency
final_data <- final_data |>
clean_names()
# Checking the column names
names(final_data) [1] "t_id" "student" "tprogram" "tgrade"
[5] "male" "ethnicity" "frlunch" "el"
[9] "prepercent" "postpercent" "white" "elementary"
[13] "math" "taught_elem" "taught_middle"## Effect Size
no_nas <- na.omit(final_data)
# effsize::cohen.d(no_nas)
# Calculating Pretest and Posttest Means
pretest_mean <- mean(final_data$prepercent, na.rm = TRUE)
posttest_mean <- mean(final_data$postpercent, na.rm = TRUE)
# Displaying the means together as a data frame
data.frame(pretest_mean, posttest_mean) pretest_mean posttest_mean
1 45.2 82.5Changing the labels and class of the columns
Sometimes, we have to change the dummy codes to the actual categories
for better understanding and interpretation of the output. The codes
below display the action where I am changing lavels into
the desired labels. They were dummy coded on the SPSS
files. The names did not get carried over to R but just the codes.
# Setting dummy codes to real names and changing the columns to factor
final_data$el <- factor(final_data$el, levels = c(0, 1), labels = c("non_ELs", "ELs"))
final_data$male <- factor(final_data$male, levels = c(0, 1), labels = c("Female", "Male"))
final_data$tgrade <- factor(final_data$tgrade, levels = c(0, 1, 2), labels = c("Elementary Grades", "Middle School", "High School"))
final_data$ethnicity <- factor(final_data$ethnicity, levels = c(0, 1, 2, 3, 4, 5), labels = c("White", "Black", "Hispanics", "Asian & Pacific Islanders", "Alaskan Natives or American Indians", "Other or Multiracial"))
final_data$frlunch <- factor(final_data$frlunch, levels = c(0, 1), labels = c("high_SES", "low_SES"))
final_data$tprogram <- factor(final_data$tprogram, levels = c(0, 1, 2), labels = c("Elementary Education", "Math Education", "English Language Arts"))
final_data$white <- factor(final_data$white, levels = c(0, 1), labels = c("Minority", "non_Minority"))Creating Long Data (Person Period Data) out of Wide Data (Person Level Data)
There are many instances when we need our data to be in a long format
(person period data - one row for each measurement/time having
multiple rows per person based on the number or times they were
measured) rather than in a wide one (person level data
- one row of data per participant, and one column per
variable), for example, while conducting ANOVA, or other
advanced regression. We can change the wide data to long data using the
pivot_longer() function in base R. However, its scope is
limited for me to use here. I want to be able to transfer some
columns as they are, while stacking some of the columns for combined
measurement. Use of `data.frame() & stack() function would allow me
to do so.
# Changing the wide data into the long data needed in some analyses
long_data <- data.frame(t_id = final_data$t_id, student = final_data$student, el = final_data$el, white = final_data$white, frlunch = final_data$frlunch, stack(final_data, select = prepercent:postpercent)) |>
mutate(
test_scores = as.numeric(values),
test_types = as.factor(ind)
) |>
select(t_id, student, el, white, frlunch, test_scores, test_types)
# changing the labels of the newly minted variable neamed test_types
long_data$test_types <- factor(long_data$test_types, labels = c("prepercent", "postpercent"))
# Checking if that worked
head(long_data) t_id student el white frlunch test_scores test_types
1 1 1 ELs Minority low_SES 50 prepercent
2 1 2 non_ELs non_Minority low_SES 80 prepercent
3 1 3 ELs non_Minority high_SES 70 prepercent
4 1 4 non_ELs Minority low_SES 60 prepercent
5 1 5 ELs Minority low_SES 70 prepercent
6 1 6 non_ELs non_Minority low_SES NA prepercentRunning some Exploratory Data Analyses using Visualization Techniques
Visualizing the distribution of pretest and posttest scores
pre_post_density <- long_data |>
na.omit() |>
ggplot(aes(test_scores, fill = test_types, label = test_types), alpha = 0.3) +
geom_textdensity(size = 4, fontface = 4, hjust = 0.2, vjust = 0.3) +
geom_textvline(xintercept = pretest_mean, label = "pretest mean score = 45.2", hjust = 0.8, vjust = 1.3) +
geom_textvline(xintercept = posttest_mean, label = "posttest mean score = 82.5", hjust = 0.8, vjust = 1.3) +
xlab("Test Scores") +
theme(legend.position = "none")
# APA Plot
pre_post_density + theme_apa()
Comparing the Pretest and Posttest Scores as Function of EL Status
# Creating EL pretest and posttest score plot
b_bar <- long_data |>
na.omit() |>
group_by(el, test_types) |>
summarize(
mean_value = mean(test_scores, na.rm = TRUE),
sd_value = sd(test_scores, na.rm = TRUE)
) |>
ggplot(aes(x = factor(test_types), mean_value, fill = test_types)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_errorbar(aes(ymin = mean_value - sd_value, ymax = mean_value + sd_value), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(values = c("darkgrey", "lightgray")) +
labs(y = "Test Scores", x = "") +
facet_wrap(~el)
# printing the plot in apa theme
b_bar + theme_apa()
Checking out various necessary statistics
# Calculating total students by teacher program
distinct_t_subject <- final_data |>
select(t_id, tprogram) |>
group_by(t_id, tprogram) |>
count() |>
print()# A tibble: 236 x 3
# Groups: t_id, tprogram [236]
t_id tprogram n
<dbl> <fct> <int>
1 1 Elementary Education 21
2 2 Elementary Education 18
3 3 Elementary Education 18
4 4 Elementary Education 22
5 5 Elementary Education 14
6 6 Elementary Education 16
7 7 Elementary Education 19
8 8 Elementary Education 15
9 9 Elementary Education 4
10 10 Elementary Education 17
# ... with 226 more rows# Calculating total teachers by teacher program
teacher_count_program <- distinct_t_subject |>
select(tprogram) |>
group_by(tprogram) |>
count() |>
mutate(pct = n / 236) |>
print()# A tibble: 3 x 3
# Groups: tprogram [3]
tprogram n pct
<fct> <int> <dbl>
1 Elementary Education 220 0.932
2 Math Education 5 0.0212
3 English Language Arts 11 0.0466# Calculating total students by teacher grades
distinct_grade <- final_data |>
select(tgrade, t_id) |>
group_by(t_id, tgrade) |>
count() |>
print()# A tibble: 236 x 3
# Groups: t_id, tgrade [236]
t_id tgrade n
<dbl> <fct> <int>
1 1 Elementary Grades 21
2 2 Elementary Grades 18
3 3 Elementary Grades 18
4 4 Elementary Grades 22
5 5 Elementary Grades 14
6 6 Elementary Grades 16
7 7 Elementary Grades 19
8 8 Elementary Grades 15
9 9 Elementary Grades 4
10 10 Elementary Grades 17
# ... with 226 more rows# Calculating total teacher by the grades they taught
teacher_count_grade <- distinct_grade |>
select(tgrade) |>
group_by(tgrade) |>
count() |>
mutate(pct = n / 236) |>
print()# A tibble: 3 x 3
# Groups: tgrade [3]
tgrade n pct
<fct> <int> <dbl>
1 Elementary Grades 218 0.924
2 Middle School 7 0.0297
3 High School 11 0.0466Regressions
Getting Rid of Missing Values for easy calculations
final_data <- na.omit(final_data)Prologue to Regression
pretest_fixed <- gls(prepercent ~ 1, data = final_data, method = "ML")
summary(pretest_fixed)Generalized least squares fit by maximum likelihood
Model: prepercent ~ 1
Data: final_data
AIC BIC logLik
50357 50370 -25177
Coefficients:
Value Std.Error t-value p-value
(Intercept) 45.2 0.327 138 0
Standardized residuals:
Min Q1 Med Q3 Max
-1.8703 -0.7650 -0.0571 0.7170 2.2694
Residual standard error: 24.2
Degrees of freedom: 5469 total; 5468 residualpretest_random_intercept <- nlme::lme(prepercent ~ 1, data = final_data, random = ~ 1 | t_id, method = "ML")
summary(pretest_random_intercept)Linear mixed-effects model fit by maximum likelihood
Data: final_data
AIC BIC logLik
47664 47683 -23829
Random effects:
Formula: ~1 | t_id
(Intercept) Residual
StdDev: 17 17.7
Fixed effects: prepercent ~ 1
Value Std.Error DF t-value p-value
(Intercept) 47.2 1.14 5233 41.5 0
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.0382 -0.6263 0.0173 0.6169 4.3799
Number of Observations: 5469
Number of Groups: 236 ## Comparing fixed intercept and random intercept null models
anova(pretest_fixed, pretest_random_intercept) Model df AIC BIC logLik Test L.Ratio p-value
pretest_fixed 1 2 50357 50370 -25177
pretest_random_intercept 2 3 47664 47683 -23829 1 vs 2 2696 <.0001## Calculating R-squared and ICC
summary(pretest_fixed)$r.squaredNULLsummary(pretest_random_intercept)$r.squaredNULLA. Pretest Null Model
# Null Model
pretest_null_model <- lmer(prepercent ~ 1 + (1 | t_id), data = final_data)
# Printing Summary Table
# summary(pretest_null_model)
# R-squared
# r2mlm(pretest_null_model)
# sjPlot::tab_model(pretest_null_model)A.1. Confidence Intervals
confint(pretest_null_model) 2.5 % 97.5 %
.sig01 15.5 18.7
.sigma 17.4 18.0
(Intercept) 44.9 49.4B. Posttest Null Model
# Null Model
posttest_null_model <- lmer(postpercent ~ 1 + (1 | t_id), data = final_data)
# Printing Summary Table
summary(posttest_null_model)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: postpercent ~ 1 + (1 | t_id)
Data: final_data
REML criterion at convergence: 45393
Scaled residuals:
Min 1Q Median 3Q Max
-5.975 -0.490 0.201 0.655 2.810
Random effects:
Groups Name Variance Std.Dev.
t_id (Intercept) 85.6 9.25
Residual 214.1 14.63
Number of obs: 5469, groups: t_id, 236
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 82.36 0.64 232.38 129 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# R-squared
r2mlm(posttest_null_model)
$Decompositions
total within between
fixed, within 0 0 NA
fixed, between 0 NA 0
slope variation 0 0 NA
mean variation 0.285498104191483 NA 1
sigma2 0.714501895808518 1 NA
$R2s
total within between
f1 0 0 NA
f2 0 NA 0
v 0 0 NA
m 0.285498104191483 NA 1
f 0 NA NA
fv 0 0 NA
fvm 0.285498104191483 NA NA # sjPlot::tab_model(posttest_null_model)B.1. Confidence Intervals Posttest Null Model
confint(posttest_null_model) 2.5 % 97.5 %
.sig01 8.34 10.2
.sigma 14.36 14.9
(Intercept) 81.10 83.6C. EL Pretest Fixed Slope Model
el_pretest_model <- lmer(prepercent ~ el + (1 | t_id), data = final_data)
# summary(el_pretest_model)
# r2mlm(el_pretest_model)
# sjPlot::tab_model(el_pretest_model)
effectsize::standardize_parameters(el_pretest_model, method = "pseudo")# Standardization method: pseudo
Parameter | Coefficient (std.) | 95% CI
-------------------------------------------------
(Intercept) | 0.00 | [ 0.00, 0.00]
elELs | -0.15 | [-0.18, -0.12]D. EL Pretest Random Slope Model
el_pretest_random <- lmer(prepercent ~ el + (el | t_id), data = final_data)
# sjPlot::tab_model(el_pretest_random)
effectsize::standardize_parameters(el_pretest_random, method = "pseudo")# Standardization method: pseudo
Parameter | Coefficient (std.) | 95% CI
-------------------------------------------------
(Intercept) | 0.00 | [ 0.00, 0.00]
elELs | -0.16 | [-0.19, -0.13]D.1. EL Pretest Random Slope Confidence Interval
confint(el_pretest_random) 2.5 % 97.5 %
.sig01 15.64 18.94
.sig02 -0.86 -0.13
.sig03 2.85 8.29
.sigma 17.12 17.79
(Intercept) 45.92 50.45
elELs -11.08 -7.24E. Putting them Together
sjPlot::tab_model(pretest_null_model, el_pretest_model, el_pretest_random)| prepercent | prepercent | prepercent | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 47.18 | 44.95 – 49.41 | <0.001 | 48.15 | 45.92 – 50.38 | <0.001 | 48.18 | 45.92 – 50.45 | <0.001 |
| el [ELs] | -8.59 | -10.23 – -6.94 | <0.001 | -9.15 | -11.05 – -7.26 | <0.001 | |||
| Random Effects | |||||||||
| σ2 | 313.29 | 307.37 | 304.48 | ||||||
| τ00 | 289.64 t_id | 286.11 t_id | 296.15 t_id | ||||||
| τ11 | 33.78 t_id.elELs | ||||||||
| ρ01 | -0.46 t_id | ||||||||
| ICC | 0.48 | 0.48 | 0.49 | ||||||
| N | 236 t_id | 236 t_id | 236 t_id | ||||||
| Observations | 5469 | 5469 | 5469 | ||||||
| Marginal R2 / Conditional R2 | 0.000 / 0.480 | 0.012 / 0.488 | 0.013 / 0.495 | ||||||
F. EL FRPL Models
el_frpl_model <- lmer(prepercent ~ el + frlunch + (el + frlunch | t_id), data = final_data)
el_frpl_interaction <- lmer(prepercent ~ el + frlunch + el * frlunch + (el + frlunch | t_id), data = final_data)
# sjPlot::tab_model(el_frpl_model, el_frpl_interaction)
plot_model(el_frpl_interaction, type = "pred", terms = c("el", "frlunch"))
G. Level 1 Models
level_1_model <- lmer(prepercent ~ el + frlunch + male + white + (el + frlunch | t_id), data = final_data)
# sjPlot::tab_model(el_frpl_model, el_frpl_interaction, level_1_model)H. Level 1 Level 2 Models
level1_level2_model <- lmer(prepercent ~ el + frlunch + male + ethnicity + (el + frlunch | t_id) + (1 + male | tprogram) + (1 + male | tgrade), data = final_data, control = lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5)))
summary(level1_level2_model)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: prepercent ~ el + frlunch + male + ethnicity + (el + frlunch |
t_id) + (1 + male | tprogram) + (1 + male | tgrade)
Data: final_data
Control: lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e+05))
REML criterion at convergence: 47346
Scaled residuals:
Min 1Q Median 3Q Max
-4.161 -0.614 0.009 0.616 4.776
Random effects:
Groups Name Variance Std.Dev. Corr
t_id (Intercept) 292.54 17.10
elELs 28.68 5.36 -0.54
frlunchlow_SES 31.92 5.65 -0.26 0.49
tprogram (Intercept) 0.00 0.00
maleMale 21.93 4.68 NaN
tgrade (Intercept) 55.59 7.46
maleMale 8.58 2.93 0.10
Residual 291.28 17.07
Number of obs: 5469, groups: t_id, 236; tprogram, 3; tgrade, 3
Fixed effects:
Estimate Std. Error df
(Intercept) 45.583 4.997 1.550
elELs -7.898 0.988 154.086
frlunchlow_SES -4.802 0.759 145.982
maleMale 0.837 3.342 0.907
ethnicityBlack -2.292 0.648 5267.155
ethnicityHispanics 13.729 4.503 5117.673
ethnicityAsian & Pacific Islanders -4.350 0.768 5207.432
ethnicityAlaskan Natives or American Indians 3.154 1.153 5169.965
ethnicityOther or Multiracial -0.182 1.357 5176.037
t value Pr(>|t|)
(Intercept) 9.12 0.02486 *
elELs -8.00 2.8e-13 ***
frlunchlow_SES -6.33 2.8e-09 ***
maleMale 0.25 0.84700
ethnicityBlack -3.54 0.00041 ***
ethnicityHispanics 3.05 0.00231 **
ethnicityAsian & Pacific Islanders -5.66 1.6e-08 ***
ethnicityAlaskan Natives or American Indians 2.74 0.00626 **
ethnicityOther or Multiracial -0.13 0.89329
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) elELs fr_SES maleMl ethncB ethncH etA&PI eANoAI
elELs -0.049
frlnchl_SES -0.080 0.023
maleMale 0.026 0.009 0.001
ethnctyBlck -0.034 -0.281 -0.126 -0.001
ethnctyHspn -0.006 -0.001 -0.001 -0.002 0.058
ethnctyA&PI -0.026 0.005 -0.142 0.001 0.336 0.043
ethnctANoAI -0.018 -0.098 0.025 -0.007 0.214 0.026 0.143
ethnctyOtoM -0.014 -0.043 -0.042 0.002 0.173 0.019 0.144 0.079
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')# Putting Results together
sjPlot::tab_model(pretest_null_model, level_1_model, level1_level2_model)| prepercent | prepercent | prepercent | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 47.18 | 44.95 – 49.41 | <0.001 | 49.02 | 46.48 – 51.56 | <0.001 | 45.58 | 35.79 – 55.38 | <0.001 |
| el [ELs] | -7.50 | -9.37 – -5.62 | <0.001 | -7.90 | -9.83 – -5.96 | <0.001 | |||
| frlunch [low SES] | -5.38 | -6.88 – -3.88 | <0.001 | -4.80 | -6.29 – -3.31 | <0.001 | |||
| male [Male] | 1.51 | 0.58 – 2.44 | 0.001 | 0.84 | -5.71 – 7.39 | 0.802 | |||
| white [non Minority] | 1.95 | 0.89 – 3.00 | <0.001 | ||||||
| ethnicity [Black] | -2.29 | -3.56 – -1.02 | <0.001 | ||||||
| ethnicity [Hispanics] | 13.73 | 4.90 – 22.56 | 0.002 | ||||||
|
ethnicity [Asian & Pacific Islanders] |
-4.35 | -5.86 – -2.84 | <0.001 | ||||||
|
ethnicity [Alaskan Natives or American Indians] |
3.15 | 0.89 – 5.42 | 0.006 | ||||||
|
ethnicity [Other or Multiracial] |
-0.18 | -2.84 – 2.48 | 0.893 | ||||||
| Random Effects | |||||||||
| σ2 | 313.29 | 294.13 | 291.28 | ||||||
| τ00 | 289.64 t_id | 307.66 t_id | 292.54 t_id | ||||||
| 0.00 tprogram | |||||||||
| 55.59 tgrade | |||||||||
| τ11 | 26.88 t_id.elELs | 28.68 t_id.elELs | |||||||
| 34.22 t_id.frlunchlow_SES | 31.92 t_id.frlunchlow_SES | ||||||||
| 21.93 tprogram.maleMale | |||||||||
| 8.58 tgrade.maleMale | |||||||||
| ρ01 | -0.53 | -0.54 t_id.elELs | |||||||
| -0.30 | -0.26 t_id.frlunchlow_SES | ||||||||
| 0.10 tgrade | |||||||||
| ICC | 0.48 | 0.50 | |||||||
| N | 236 t_id | 236 t_id | 236 t_id | ||||||
| 3 tprogram | |||||||||
| 3 tgrade | |||||||||
| Observations | 5469 | 5469 | 5469 | ||||||
| Marginal R2 / Conditional R2 | 0.000 / 0.480 | 0.032 / 0.511 | 0.069 / NA | ||||||
H.1. Confidence Interval Pretest Final with Ethnicity
# confint(level1_level2_model) This code gave an error. Thus, below method is used. The values, though, may differ from other method like method = "boot"
confint.merMod(level1_level2_model, method = "Wald") 2.5 % 97.5 %
.sig01 NA NA
.sig02 NA NA
.sig03 NA NA
.sig04 NA NA
.sig05 NA NA
.sig06 NA NA
.sig07 NA NA
.sig08 NA NA
.sig09 NA NA
.sig10 NA NA
.sig11 NA NA
.sig12 NA NA
.sigma NA NA
(Intercept) 35.788 55.38
elELs -9.833 -5.96
frlunchlow_SES -6.288 -3.31
maleMale -5.713 7.39
ethnicityBlack -3.563 -1.02
ethnicityHispanics 4.903 22.56
ethnicityAsian & Pacific Islanders -5.855 -2.84
ethnicityAlaskan Natives or American Indians 0.894 5.41
ethnicityOther or Multiracial -2.842 2.48I. Extra Model Comparison
options(scipen = 999) # To get rid of scientific notations
anova(pretest_null_model, level_1_model, level1_level2_model) # Comparing Pretest ModelsData: final_data
Models:
pretest_null_model: prepercent ~ 1 + (1 | t_id)
level_1_model: prepercent ~ el + frlunch + male + white + (el + frlunch | t_id)
level1_level2_model: prepercent ~ el + frlunch + male + ethnicity + (el + frlunch | t_id) + (1 + male | tprogram) + (1 + male | tgrade)
npar AIC BIC logLik deviance Chisq Df
pretest_null_model 3 47664 47683 -23829 47658
level_1_model 12 47442 47521 -23709 47418 240 9
level1_level2_model 22 47412 47557 -23684 47368 50 10
Pr(>Chisq)
pretest_null_model
level_1_model < 0.0000000000000002 ***
level1_level2_model 0.00000027 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(cowplot) # To put the plots togetherJ. Final
# plot(final_data)
# Checking the fixed and random effects of all predictors(Results not populated because of space issue)
# head(coef(level1_level2_model))
# Plotting the Random effects of the mentioned model (only Random effects)
pre_random_plot <- plot(ranef(level1_level2_model))
pre_random_plot$t_id # Just Printing random plots for t_id as grouping variable. It also had random effects for t_program and t_grade
# Plotting the whole model(fitted vs. residuals)
plot(level1_level2_model)
K. R-Squared
# Calculating R-Squared for the Mentioned Models
r2(pretest_null_model)# R2 for Mixed Models
Conditional R2: 0.480
Marginal R2: 0.000r2(level_1_model)# R2 for Mixed Models
Conditional R2: 0.511
Marginal R2: 0.032r2(level1_level2_model)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.069L. Level 1 Interaction Model
level_1_model_int <- lmer(prepercent ~ el + frlunch + male + white + el * frlunch + el * male + el * white + frlunch * male + frlunch * white + male * white + (el + frlunch | t_id), data = final_data)
summary(level_1_model_int)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: prepercent ~ el + frlunch + male + white + el * frlunch + el *
male + el * white + frlunch * male + frlunch * white + male *
white + (el + frlunch | t_id)
Data: final_data
REML criterion at convergence: 47389
Scaled residuals:
Min 1Q Median 3Q Max
-4.243 -0.602 0.009 0.610 4.613
Random effects:
Groups Name Variance Std.Dev. Corr
t_id (Intercept) 307.7 17.54
elELs 27.9 5.28 -0.52
frlunchlow_SES 34.0 5.83 -0.30 0.50
Residual 294.0 17.15
Number of obs: 5469, groups: t_id, 236
Fixed effects:
Estimate Std. Error df t value
(Intercept) 49.2810 1.4090 374.3784 34.98
elELs -8.1092 1.8916 436.5526 -4.29
frlunchlow_SES -6.2253 1.0909 502.9885 -5.71
maleMale 0.7262 0.9374 5156.7984 0.77
whitenon_Minority 2.5124 0.9090 5156.5537 2.76
elELs:frlunchlow_SES 0.8634 1.9826 531.6911 0.44
elELs:maleMale -0.0809 1.6659 3640.8564 -0.05
elELs:whitenon_Minority 0.6052 3.1068 2191.4970 0.19
frlunchlow_SES:maleMale 2.2185 0.9991 5161.7614 2.22
frlunchlow_SES:whitenon_Minority -0.6039 1.0727 4929.4219 -0.56
maleMale:whitenon_Minority -0.6818 1.0273 5165.9737 -0.66
Pr(>|t|)
(Intercept) < 0.0000000000000002 ***
elELs 0.00002230 ***
frlunchlow_SES 0.00000002 ***
maleMale 0.4386
whitenon_Minority 0.0057 **
elELs:frlunchlow_SES 0.6634
elELs:maleMale 0.9613
elELs:whitenon_Minority 0.8456
frlunchlow_SES:maleMale 0.0264 *
frlunchlow_SES:whitenon_Minority 0.5735
maleMale:whitenon_Minority 0.5069
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) elELs fr_SES maleMl whtn_M eEL:_S elEL:M eEL:_M f_SES:M
elELs -0.263
frlnchl_SES -0.501 0.205
maleMale -0.340 0.127 0.311
whtnn_Mnrty -0.422 0.248 0.438 0.394
elELs:f_SES 0.146 -0.753 -0.294 -0.027 -0.173
elELs:malMl 0.084 -0.385 0.042 -0.241 -0.133 -0.014
elELs:wht_M 0.028 -0.166 0.030 -0.040 -0.089 0.002 0.085
frlnc_SES:M 0.216 0.016 -0.451 -0.646 -0.145 0.039 -0.105 0.001
frln_SES:_M 0.253 -0.140 -0.525 -0.011 -0.567 0.221 -0.004 -0.086 -0.003
mlMl:whtn_M 0.234 -0.111 -0.128 -0.687 -0.564 0.009 0.238 0.069 0.245
f_SES:_
elELs
frlnchl_SES
maleMale
whtnn_Mnrty
elELs:f_SES
elELs:malMl
elELs:wht_M
frlnc_SES:M
frln_SES:_M
mlMl:whtn_M 0.005 r2(level_1_model_int)# R2 for Mixed Models
Conditional R2: 0.512
Marginal R2: 0.032# Plotting all the interactions in the mentioned models
# theme_set(theme_sjplot())
# Plotting all the interactions
pretest_interaction_plot <- plot_model(level_1_model_int, type = "int")L.1. Interaction of EL with Other Variables in Predicting Pretest Scores
# Using Combined Title
title1 <- ggdraw() +
draw_label("Interactional Effects of EL and Other L1\n Variables in Predicting Pretest Scores", x = 0, hjust = 0) + theme(plot.margin = margin(0, 5, 0, 5))
plot_grid(title1, pretest_interaction_plot[[1]], pretest_interaction_plot[[2]], pretest_interaction_plot[[3]], ncol = 2, rel_heights = c(0.5, 0.5))
L.2. Interaction of FRPL and Other Variables in Predicting Pretest Scores
plot_grid(pretest_interaction_plot[[4]], pretest_interaction_plot[[5]], pretest_interaction_plot[[6]], align = "h")
M. Posttest Model
posttest_model <- lmer(postpercent ~ prepercent + el + frlunch + male + ethnicity + (prepercent + el + frlunch | t_id) + (1 + male | tprogram) + (1 + male | tgrade), data = final_data, control = lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5)))
summary(posttest_model)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: postpercent ~ prepercent + el + frlunch + male + ethnicity +
(prepercent + el + frlunch | t_id) + (1 + male | tprogram) +
(1 + male | tgrade)
Data: final_data
Control: lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 200000))
REML criterion at convergence: 44050
Scaled residuals:
Min 1Q Median 3Q Max
-6.754 -0.464 0.124 0.610 3.948
Random effects:
Groups Name Variance Std.Dev. Corr
t_id (Intercept) 217.0624 14.733
prepercent 0.0188 0.137 -1.00
elELs 30.2508 5.500 0.28 -0.28
frlunchlow_SES 11.2576 3.355 -0.16 0.16 -0.30
tprogram (Intercept) 129.4065 11.376
maleMale 62.6082 7.913 0.14
tgrade (Intercept) 102.7192 10.135
maleMale 100.3981 10.020 -0.09
Residual 163.1701 12.774
Number of obs: 5469, groups: t_id, 236; tprogram, 3; tgrade, 3
Fixed effects:
Estimate Std. Error df
(Intercept) 65.3276 8.9612 751.9088
prepercent 0.3196 0.0134 145.3309
elELs -5.1100 0.8228 111.0696
frlunchlow_SES -1.4796 0.5117 140.9150
maleMale 1.0306 7.4138 36.0167
ethnicityBlack -0.4624 0.4821 5277.5468
ethnicityHispanics -2.2567 3.3180 4810.2639
ethnicityAsian & Pacific Islanders -1.2604 0.5735 5238.3554
ethnicityAlaskan Natives or American Indians -0.0848 0.8563 4898.7513
ethnicityOther or Multiracial 0.1756 1.0118 5243.3166
t value Pr(>|t|)
(Intercept) 7.29 0.00000000000079 ***
prepercent 23.90 < 0.0000000000000002 ***
elELs -6.21 0.00000000940898 ***
frlunchlow_SES -2.89 0.0044 **
maleMale 0.14 0.8902
ethnicityBlack -0.96 0.3375
ethnicityHispanics -0.68 0.4964
ethnicityAsian & Pacific Islanders -2.20 0.0280 *
ethnicityAlaskan Natives or American Indians -0.10 0.9211
ethnicityOther or Multiracial 0.17 0.8622
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) prprcn elELs fr_SES maleMl ethncB ethncH etA&PI eANoAI
prepercent -0.116
elELs 0.007 -0.014
frlnchl_SES -0.032 0.129 -0.108
maleMale 0.015 0.000 0.003 0.000
ethnctyBlck -0.016 0.034 -0.241 -0.148 -0.001
ethnctyHspn -0.002 -0.025 0.000 -0.010 0.000 0.052
ethnctyA&PI -0.015 0.056 0.020 -0.160 0.000 0.336 0.037
ethnctANoAI -0.009 -0.026 -0.087 0.022 -0.003 0.204 0.026 0.135
ethnctyOtoM -0.007 0.004 -0.035 -0.043 0.001 0.168 0.018 0.141 0.076
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')sjPlot::tab_model(posttest_model)| postpercent | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 65.33 | 47.76 – 82.90 | <0.001 |
| prepercent | 0.32 | 0.29 – 0.35 | <0.001 |
| el [ELs] | -5.11 | -6.72 – -3.50 | <0.001 |
| frlunch [low SES] | -1.48 | -2.48 – -0.48 | 0.004 |
| male [Male] | 1.03 | -13.50 – 15.56 | 0.889 |
| ethnicity [Black] | -0.46 | -1.41 – 0.48 | 0.338 |
| ethnicity [Hispanics] | -2.26 | -8.76 – 4.25 | 0.496 |
|
ethnicity [Asian & Pacific Islanders] |
-1.26 | -2.38 – -0.14 | 0.028 |
|
ethnicity [Alaskan Natives or American Indians] |
-0.08 | -1.76 – 1.59 | 0.921 |
|
ethnicity [Other or Multiracial] |
0.18 | -1.81 – 2.16 | 0.862 |
| Random Effects | |||
| σ2 | 163.17 | ||
| τ00 t_id | 217.06 | ||
| τ00 tprogram | 129.41 | ||
| τ00 tgrade | 102.72 | ||
| τ11 t_id.prepercent | 0.02 | ||
| τ11 t_id.elELs | 30.25 | ||
| τ11 t_id.frlunchlow_SES | 11.26 | ||
| τ11 tprogram.maleMale | 62.61 | ||
| τ11 tgrade.maleMale | 100.40 | ||
| ρ01 t_id.prepercent | -1.00 | ||
| ρ01 t_id.elELs | 0.28 | ||
| ρ01 t_id.frlunchlow_SES | -0.16 | ||
| ρ01 tprogram | 0.14 | ||
| ρ01 tgrade | -0.09 | ||
| N t_id | 236 | ||
| N tprogram | 3 | ||
| N tgrade | 3 | ||
| Observations | 5469 | ||
| Marginal R2 / Conditional R2 | 0.300 / NA | ||
r2(posttest_model)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.300# plot_model(posttest_model, type="pred", colors = "bw")
effectsize::omega_squared(posttest_model, ci.lvl = .95)# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI
--------------------------------------------
prepercent | 0.79 | [0.75, 1.00]
el | 0.25 | [0.14, 1.00]
frlunch | 0.05 | [0.01, 1.00]
male | -0.03 | [0.00, 1.00]
ethnicity | 1.11e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at (1).N. Posttest Final
posttest_model_final <- lmer(postpercent ~ prepercent + el + frlunch + male + white + (prepercent + el + frlunch | t_id) + (1 + male | tprogram) + (1 + male | tgrade), data = final_data, control = lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5)))
summary(posttest_model_final)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
postpercent ~ prepercent + el + frlunch + male + white + (prepercent +
el + frlunch | t_id) + (1 + male | tprogram) + (1 + male | tgrade)
Data: final_data
Control: lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 200000))
REML criterion at convergence: 44052
Scaled residuals:
Min 1Q Median 3Q Max
-6.691 -0.462 0.125 0.609 3.881
Random effects:
Groups Name Variance Std.Dev. Corr
t_id (Intercept) 191.6689 13.844
prepercent 0.0163 0.128 -1.00
elELs 31.8890 5.647 0.33 -0.33
frlunchlow_SES 11.0141 3.319 -0.08 0.08 -0.35
tprogram (Intercept) 56.1465 7.493
maleMale 2.9007 1.703 0.16
tgrade (Intercept) 40.3482 6.352
maleMale 16.1061 4.013 -0.03
Residual 163.6160 12.791
Number of obs: 5469, groups: t_id, 236; tprogram, 3; tgrade, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 64.6258 5.9049 4.2874 10.94 0.00027
prepercent 0.3201 0.0129 139.8924 24.82 < 0.0000000000000002
elELs -4.9208 0.8106 102.6222 -6.07 0.000000022
frlunchlow_SES -1.5760 0.5062 139.1503 -3.11 0.00224
maleMale 1.2611 2.6113 0.0543 0.48 0.92401
whitenon_Minority 0.6154 0.3985 5243.9647 1.54 0.12256
(Intercept) ***
prepercent ***
elELs ***
frlunchlow_SES **
maleMale
whitenon_Minority
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) prprcn elELs fr_SES maleMl
prepercent -0.164
elELs 0.005 -0.034
frlnchl_SES -0.051 0.111 -0.119
maleMale 0.012 -0.001 0.007 0.001
whtnn_Mnrty -0.037 -0.038 0.166 0.161 0.003
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')sjPlot::tab_model(posttest_model_final)| postpercent | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 64.63 | 53.05 – 76.20 | <0.001 |
| prepercent | 0.32 | 0.29 – 0.35 | <0.001 |
| el [ELs] | -4.92 | -6.51 – -3.33 | <0.001 |
| frlunch [low SES] | -1.58 | -2.57 – -0.58 | 0.002 |
| male [Male] | 1.26 | -3.86 – 6.38 | 0.629 |
| white [non Minority] | 0.62 | -0.17 – 1.40 | 0.123 |
| Random Effects | |||
| σ2 | 163.62 | ||
| τ00 t_id | 191.67 | ||
| τ00 tprogram | 56.15 | ||
| τ00 tgrade | 40.35 | ||
| τ11 t_id.prepercent | 0.02 | ||
| τ11 t_id.elELs | 31.89 | ||
| τ11 t_id.frlunchlow_SES | 11.01 | ||
| τ11 tprogram.maleMale | 2.90 | ||
| τ11 tgrade.maleMale | 16.11 | ||
| ρ01 t_id.prepercent | -1.00 | ||
| ρ01 t_id.elELs | 0.33 | ||
| ρ01 t_id.frlunchlow_SES | -0.08 | ||
| ρ01 tprogram | 0.16 | ||
| ρ01 tgrade | -0.03 | ||
| N t_id | 236 | ||
| N tprogram | 3 | ||
| N tgrade | 3 | ||
| Observations | 5469 | ||
| Marginal R2 / Conditional R2 | 0.300 / NA | ||
r2(posttest_model_final)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.300plot_model(posttest_model_final, type = "pred")$prepercent
$el
$frlunch
$male
$white
effectsize::omega_squared(posttest_model_final, ci.lvl = .95)# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI
--------------------------------------------
prepercent | 0.81 | [0.77, 1.00]
el | 0.26 | [0.14, 1.00]
frlunch | 0.06 | [0.01, 1.00]
male | -0.60 | [0.00, 1.00]
white | 2.64e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at (1).O. Posttest, pretest*el
posttest_pretest_el <- lmer(postpercent ~ prepercent + el + frlunch + male + white + prepercent * el + prepercent * frlunch + prepercent * male + prepercent * white + el * frlunch + el * white + el * male + frlunch * white + frlunch * male + white * male + (1 | t_id), data = final_data)
r2(posttest_pretest_el)# R2 for Mixed Models
Conditional R2: 0.454
Marginal R2: 0.231summary(posttest_pretest_el)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: postpercent ~ prepercent + el + frlunch + male + white + prepercent *
el + prepercent * frlunch + prepercent * male + prepercent *
white + el * frlunch + el * white + el * male + frlunch *
white + frlunch * male + white * male + (1 | t_id)
Data: final_data
REML criterion at convergence: 44194
Scaled residuals:
Min 1Q Median 3Q Max
-6.146 -0.458 0.089 0.626 3.480
Random effects:
Groups Name Variance Std.Dev.
t_id (Intercept) 70.1 8.37
Residual 171.7 13.10
Number of obs: 5469, groups: t_id, 236
Fixed effects:
Estimate Std. Error df t value
(Intercept) 67.7780 1.1912 2613.2557 56.90
prepercent 0.3190 0.0183 5428.0963 17.38
elELs -7.2096 1.7888 5296.5945 -4.03
frlunchlow_SES -3.1426 1.0606 5395.6418 -2.96
maleMale 1.4159 1.0212 5260.6555 1.39
whitenon_Minority 2.0176 1.0597 5294.6704 1.90
prepercent:elELs 0.0695 0.0279 5320.7234 2.49
prepercent:frlunchlow_SES 0.0461 0.0174 5428.6126 2.65
prepercent:maleMale -0.0237 0.0153 5261.1301 -1.55
prepercent:whitenon_Minority -0.0239 0.0167 5305.0136 -1.43
elELs:frlunchlow_SES -0.6092 1.4521 5298.8825 -0.42
elELs:whitenon_Minority 2.8501 2.2868 5271.9067 1.25
elELs:maleMale 0.6731 1.2550 5264.4600 0.54
frlunchlow_SES:whitenon_Minority -1.2322 0.8082 5311.6098 -1.52
frlunchlow_SES:maleMale 0.5713 0.7659 5270.0186 0.75
maleMale:whitenon_Minority 0.8050 0.7800 5265.9265 1.03
Pr(>|t|)
(Intercept) < 0.0000000000000002 ***
prepercent < 0.0000000000000002 ***
elELs 0.000056 ***
frlunchlow_SES 0.0031 **
maleMale 0.1656
whitenon_Minority 0.0570 .
prepercent:elELs 0.0127 *
prepercent:frlunchlow_SES 0.0081 **
prepercent:maleMale 0.1209
prepercent:whitenon_Minority 0.1515
elELs:frlunchlow_SES 0.6748
elELs:whitenon_Minority 0.2127
elELs:maleMale 0.5917
frlunchlow_SES:whitenon_Minority 0.1274
frlunchlow_SES:maleMale 0.4558
maleMale:whitenon_Minority 0.3021
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1a <- plot_model(posttest_pretest_el, type = "int")
a[[1]] + theme_apa()
a[[2]] + theme_apa()
a[[3]] + theme_apa()
a[[4]] + theme_apa()
a[[5]] + theme_apa()
a[[6]] + theme_apa()
a[[7]] + theme_apa()
a[[8]] + theme_apa()
a[[9]] + theme_apa()
a[[10]] + theme_apa()
N. Pretest Final with White
pretest_final_white <- lmer(prepercent ~ el + frlunch + male + white + (el + frlunch | t_id) + (1 + male | tprogram) + (1 + male | tgrade), data = final_data, control = lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5)))
summary(pretest_final_white)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: prepercent ~ el + frlunch + male + white + (el + frlunch | t_id) +
(1 + male | tprogram) + (1 + male | tgrade)
Data: final_data
Control: lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 200000))
REML criterion at convergence: 47406
Scaled residuals:
Min 1Q Median 3Q Max
-4.203 -0.603 0.009 0.609 4.732
Random effects:
Groups Name Variance Std.Dev. Corr
t_id (Intercept) 295.64 17.19
elELs 27.63 5.26 -0.52
frlunchlow_SES 34.06 5.84 -0.29 0.54
tprogram (Intercept) 0.00 0.00
maleMale 23.04 4.80 NaN
tgrade (Intercept) 55.63 7.46
maleMale 9.14 3.02 0.10
Residual 293.59 17.13
Number of obs: 5469, groups: t_id, 236; tprogram, 3; tgrade, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 43.657 5.011 1.571 8.71 0.02571 *
elELs -7.511 0.962 138.599 -7.81 0.0000000000013 ***
frlunchlow_SES -5.360 0.764 143.991 -7.02 0.0000000000821 ***
maleMale 0.969 3.426 0.967 0.28 0.82575
whitenon_Minority 1.967 0.539 5247.316 3.65 0.00027 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) elELs fr_SES maleMl
elELs -0.070
frlnchl_SES -0.100 0.041
maleMale 0.028 0.009 0.002
whtnn_Mnrty -0.069 0.203 0.138 0.002
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')sjPlot::tab_model(pretest_final_white)| prepercent | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 43.66 | 33.83 – 53.48 | <0.001 |
| el [ELs] | -7.51 | -9.40 – -5.63 | <0.001 |
| frlunch [low SES] | -5.36 | -6.86 – -3.86 | <0.001 |
| male [Male] | 0.97 | -5.75 – 7.68 | 0.777 |
| white [non Minority] | 1.97 | 0.91 – 3.02 | <0.001 |
| Random Effects | |||
| σ2 | 293.59 | ||
| τ00 t_id | 295.64 | ||
| τ00 tprogram | 0.00 | ||
| τ00 tgrade | 55.63 | ||
| τ11 t_id.elELs | 27.63 | ||
| τ11 t_id.frlunchlow_SES | 34.06 | ||
| τ11 tprogram.maleMale | 23.04 | ||
| τ11 tgrade.maleMale | 9.14 | ||
| ρ01 t_id.elELs | -0.52 | ||
| ρ01 t_id.frlunchlow_SES | -0.29 | ||
| ρ01 tprogram | |||
| ρ01 tgrade | 0.10 | ||
| N t_id | 236 | ||
| N tprogram | 3 | ||
| N tgrade | 3 | ||
| Observations | 5469 | ||
| Marginal R2 / Conditional R2 | 0.060 / NA | ||
r2(pretest_final_white)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.060# plot_model(posttest_model_final, type="pred")
effectsize::omega_squared(pretest_final_white, ci.lvl = .95)# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI
-------------------------------------------
el | 0.30 | [0.20, 1.00]
frlunch | 0.25 | [0.15, 1.00]
male | -0.45 | [0.00, 1.00]
white | 2.34e-03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at (1).effectsize::cohens_f(pretest_final_white, ci = .95, partial = TRUE, verbose = TRUE)# Effect Size for ANOVA (Type III)
Parameter | Cohen's f (partial) | 95% CI
---------------------------------------------
el | 0.66 | [0.51, Inf]
frlunch | 0.58 | [0.44, Inf]
male | 0.29 | [0.00, Inf]
white | 0.05 | [0.03, Inf]
- One-sided CIs: upper bound fixed at (Inf).el_es <- effectsize::cohens_d("prepercent", "el", pooled_sd = TRUE, alternative = "two.sided", data = final_data)
el_es_po <- effectsize::cohens_d("postpercent", "el", pooled_sd = TRUE, alternative = "two.sided", data = final_data)
rbind(el_es, el_es_po)Cohen's d | 95% CI
------------------------
0.42 | [0.34, 0.51]
0.48 | [0.40, 0.57]
- Estimated using pooled SD.O. Pretest Final with White Confidence Interval
# confint(level1_level2_model) This code gave an error. Thus, below method is used. The values, though, may differ from other method like method = "boot"
confint.merMod(pretest_final_white, method = "Wald") 2.5 % 97.5 %
.sig01 NA NA
.sig02 NA NA
.sig03 NA NA
.sig04 NA NA
.sig05 NA NA
.sig06 NA NA
.sig07 NA NA
.sig08 NA NA
.sig09 NA NA
.sig10 NA NA
.sig11 NA NA
.sig12 NA NA
.sigma NA NA
(Intercept) 33.84 53.48
elELs -9.40 -5.63
frlunchlow_SES -6.86 -3.86
maleMale -5.75 7.68
whitenon_Minority 0.91 3.02Checking Assumptions
I am going to test assumptions against the final posttest model because it includes all the variables in model.
A. Checking for the Outliers
It is important to note the number of grouping variables while
calculating outliers and/or the Influential
Cases in multilevel models. In general, assumption can be
tested using plot() function in R. But for the multilevel
model, having mixed effects, it is not as straightforward.
- Outliers can be assessed by looking into the standardized
residuals (using
rstandard()function), and the studentized residuals (usingrstudent()function).
- Residual is a difference between an observed value minus the
predicted value in a regression model [Residual = Observed Value
- Predicted Value]. We can access residuals by using the
resid(regression_model)function, we can call for all the residuals in any regression model. We can plot them in a histogram usinghist(resid(model))function and check for any anomaly.
hist(resid(posttest_model_final),
freq = FALSE,
main = "Density of Residuals from Final Posttest Model",
xlab = "Residuals"
)
lines(density(resid(posttest_model_final))) # adding density to get some sense of normality
# box() # drawing a box around the plot- Standardized residuals are the raw residuals divided by an estimate
of the standard deviation of the residuals. They are sometimes called
the internally studentized residuals, which is a way of
estimating error for a particular data point. Ordinary sample contains
95% of its standardized residuals within +2 and -2
limit. One of four plots as the outcomes of
plot(model)function tests standardized values against the fitted values, in which the values are supposed to fall closely within the dotted line. However, for the multilevel model, these things are not as straightforward. It is hard to getstudentized and standardized residualsbecause of the complexity of the models, e.g., levels, covariates, and even interactions. One of the quick way is to plot Binned Residuals which are calculated by ‘dividing the data into categories based on their fitted values, and then plotting the average residuals vs. the average fitted value for each bin’ (Gelman, 2007). We can create a data frame and use them to check a few other things.
# saving binned residuals in an object
b_residual <- data.frame(performance::binned_residuals(posttest_model_final))
head(b_residual) xbar ybar n x.lo x.hi se ci_range CI_low CI_high group
1 46.1 -2.7906 73 17.4 51.5 4.35 2.22 -7.14 1.556 yes
2 54.1 -4.3658 74 51.6 56.0 4.75 2.43 -9.12 0.389 yes
3 57.7 -2.1762 76 56.1 59.2 4.09 2.08 -6.26 1.910 yes
4 60.6 2.5160 72 59.2 61.9 3.81 1.94 -1.29 6.323 yes
5 63.2 0.0573 75 61.9 64.2 3.20 1.63 -3.14 3.253 yes
6 64.8 -5.8538 73 64.2 65.4 5.02 2.56 -10.87 -0.837 no# Binned Plot
arm::binnedplot(fitted(posttest_model_final), resid(posttest_model_final, type = "response"))
- In general, if a studentized residual is greater than the
absolute value of 3, it is considered an outlier. We
can definitely check the values themselves using
3 | model$studentized_residuals | -3syntax (need to be modified for exact models. In general, we check this assumption looking at some plots. We can generate a histogram of studentized residual using thehist(rstudent(model))function. The histogram should look like a normal distribution. We can double check the assumption by checking it against the normal qqplot using theqqPlot(model), in which the residuals are supposed to fall along the solid line. The overall assumption of a regression model can be populated in R using theplot(model)function. Which contains a plot featuring qqplot of studentized residuals against the theoretical quantiles.
- Influential cases, on the other hand, can be checked into
DFBETAS (using
dfbeta()function) or Leverage aka hatvalues (usinghatvalues()function), or Cook’s Distance (usingcooks.distance()function), covariance ratio (usingcovratio()function).