> rmvn.eigen <- function(n, mu, Sigma) {
+ d <- length(mu)
+ ev <- eigen(Sigma, symmetric = TRUE)
+ lambda <- ev$values
+ V <- ev$vectors
+ R <- V %*% diag(sqrt(lambda)) %*% t(V)
+ Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
+ X <- Z %*% R + matrix(mu, n, d, byrow = TRUE)
+ return(X)
+ }
> ref_FE <- function(VarRed, TargetYVar, XVar) {
+ std_VarRed <- sqrt(VarRed)
+ uns_FE <- std_VarRed * sqrt(TargetYVar)/sqrt(XVar)
+ return(uns_FE)
+ }
> data_gen <- function(nSchools, min_StuPerGrade, max_StuPerGrade) {
+ nGrades <- 4
+ NGrades <- nGrades * nSchools
+ GradeClusterSize <- replicate(n = NGrades, expr = round(runif(n = 1, 
+ 
+ min = min_StuPerGrade, max = max_StuPerGrade)))
+ sim_data0 <- cbind(school_seq = rep(1:nSchools, each = nGrades), 
+ 
+ grade_seq = 1:NGrades, GradeClusterSize)
+ sim_data1 <- NULL
+ for (r in 1:nrow(sim_data0)) {
+ 
+ sim_data1 <- rbind(sim_data1, data.frame(SchoolID = rep(sim_data0[r, 
+ 
+ 
+ 1], sim_data0[r, 3]), GradeID = rep(sim_data0[r, 
+ 
+ 
+ 2], sim_data0[r, 3])))
+ 
+ }
+ sim_data1$StudentID <- 1:nrow(sim_data1)
+ L3_data <- data.frame(SchoolID = unique(sim_data1$SchoolID))
+ L2_data <- data.frame(GradeID = unique(sim_data1$GradeID))
+ L1_data <- data.frame(StudentID = unique(sim_data1$StudentID))
+ TotalVarx1 <- 1
+ L3_PropVarx1 <- 0.06
+ L2_PropVarx1 <- 0.24
+ L3_Varx1 <- TotalVarx1 * L3_PropVarx1
+ L2_Varx1 <- TotalVarx1 * L2_PropVarx1
+ L1_Varx1 <- TotalVarx1 * (1 - L3_PropVarx1 - L2_PropVarx1)
+ L3_data$L3_x1_True <- rnorm(nSchools, 0, sqrt(L3_Varx1))
+ L2_data$L2_x1_True <- rnorm(NGrades, 0, sqrt(L2_Varx1))
+ L1_data$L1_x1_True <- rnorm(nrow(sim_data1), 0, sqrt(L1_Varx1))
+ TotalVarx2 <- 16
+ L3_PropVarx2 <- 0.12
+ L2_PropVarx2 <- 0.28
+ L3_Varx2 <- TotalVarx2 * L3_PropVarx2
+ L2_Varx2 <- TotalVarx2 * L2_PropVarx2
+ L1_Varx2 <- TotalVarx2 * (1 - L3_PropVarx2 - L2_PropVarx2)
+ L3_data$L3_x2_True <- rnorm(nSchools, 0, sqrt(L3_Varx2))
+ L2_data$L2_x2_True <- rnorm(NGrades, 0, sqrt(L2_Varx2))
+ L1_data$L1_x2_True <- rnorm(nrow(sim_data1), 0, sqrt(L1_Varx2))
+ L2_data$L2_x3_GradeLevel <- rep(9:12, times = nSchools)
+ L2_data$L2_x3_Cat2 <- ifelse(L2_data$L2_x3_GradeLevel == 
+ 
+ 10, 1, 0)
+ L2_data$L2_x3_Cat3 <- ifelse(L2_data$L2_x3_GradeLevel == 
+ 
+ 11, 1, 0)
+ L2_data$L2_x3_Cat4 <- ifelse(L2_data$L2_x3_GradeLevel == 
+ 
+ 12, 1, 0)
+ sim_data2 <- merge(sim_data1, L3_data, by = "SchoolID")
+ sim_data2 <- merge(sim_data2, L2_data, by = "GradeID")
+ sim_data2 <- merge(sim_data2, L1_data, by = "StudentID")
+ sim_data2$L1_x1_Comp <- sim_data2$L3_x1_True + sim_data2$L2_x1_True + 
+ 
+ sim_data2$L1_x1_True
+ l_bound_x1 <- 0.4
+ u_bound_x1 <- 0.75
+ sim_data2$L1_x1_Comp_Cat <- ifelse(sim_data2$L1_x1_Comp < 
+ 
+ qnorm(l_bound_x1), 1, ifelse(sim_data2$L1_x1_Comp > qnorm(u_bound_x1), 
+ 
+ 3, 2))
+ sim_data2$L1_x1_Cat2 <- ifelse(sim_data2$L1_x1_Comp_Cat == 
+ 
+ 2, 1, 0)
+ sim_data2$L1_x1_Cat3 <- ifelse(sim_data2$L1_x1_Comp_Cat == 
+ 
+ 3, 1, 0)
+ sim_data2$L2_x1_PropCat2 <- ave(sim_data2$L1_x1_Cat2, sim_data2$GradeID)
+ sim_data2$L2_x1_PropCat3 <- ave(sim_data2$L1_x1_Cat3, sim_data2$GradeID)
+ sim_data2$L3_x1_PropCat2 <- ave(sim_data2$L1_x1_Cat2, sim_data2$SchoolID)
+ sim_data2$L3_x1_PropCat3 <- ave(sim_data2$L1_x1_Cat3, sim_data2$SchoolID)
+ sim_data2$L1_x1_Cat2_CMC <- sim_data2$L1_x1_Cat2 - sim_data2$L2_x1_PropCat2
+ sim_data2$L1_x1_Cat3_CMC <- sim_data2$L1_x1_Cat3 - sim_data2$L2_x1_PropCat3
+ sim_data2$L2_x1_PropCat2_CMC <- sim_data2$L2_x1_PropCat2 - 
+ 
+ sim_data2$L3_x1_PropCat2
+ sim_data2$L2_x1_PropCat3_CMC <- sim_data2$L2_x1_PropCat3 - 
+ 
+ sim_data2$L3_x1_PropCat3
+ sim_data2$L1_x2_Comp <- 25 + sim_data2$L3_x2_True + sim_data2$L2_x2_True + 
+ 
+ sim_data2$L1_x2_True
+ sim_data2$L1_x2_Comp_Round.or <- round(sim_data2$L1_x2_Comp)
+ l_bound_x2 <- 10
+ u_bound_x2 <- 40
+ sim_data2$L1_x2_Comp_Round <- ifelse(sim_data2$L1_x2_Comp_Round.or < 
+ 
+ l_bound_x2, l_bound_x2, ifelse(sim_data2$L1_x2_Comp_Round.or > 
+ 
+ u_bound_x2, u_bound_x2, sim_data2$L1_x2_Comp_Round.or))
+ sim_data2$L2_x2_Mean <- ave(sim_data2$L1_x2_Comp_Round, sim_data2$GradeID)
+ sim_data2$L3_x2_Mean <- ave(sim_data2$L1_x2_Comp_Round, sim_data2$SchoolID)
+ sim_data2$L1_x2_Comp_Round_CMC <- sim_data2$L1_x2_Comp_Round - 
+ 
+ sim_data2$L2_x2_Mean
+ sim_data2$L2_x2_Mean_CMC <- sim_data2$L2_x2_Mean - sim_data2$L3_x2_Mean
+ sim_data2$L2_x2_Mean_CMCQuad <- sim_data2$L2_x2_Mean_CMC^2
+ TotalVary <- 5
+ L3_PropVary <- 0.15
+ L2_PropVary <- 0.35
+ L1_ResVar <- TotalVary * (1 - L3_PropVary - L2_PropVary)
+ L2_RanIntVar <- TotalVary * L2_PropVary
+ L3_RanIntVar <- TotalVary * L3_PropVary
+ U0s <- rnorm(nSchools, 0, sqrt(L3_RanIntVar))
+ L3_data_RE <- data.frame(SchoolID = unique(sim_data2$SchoolID), 
+ 
+ U0s)
+ sim_data2 <- merge(sim_data2, L3_data_RE, by = "SchoolID")
+ L2_x2_RSRel <- 0.5
+ L2_x2_RSVar <- (L2_x2_RSRel + (L1_ResVar/((nrow(sim_data2)/NGrades) * 
+ 
+ L1_Varx2)))/(1 - L2_x2_RSRel)
+ L2_mu <- c(0, 0)
+ L2_Sigma <- matrix(c(L2_RanIntVar, 0, 0, L2_x2_RSVar), nrow = 2, 
+ 
+ ncol = 2)
+ L2_data_RE <- rmvn.eigen(NGrades, L2_mu, L2_Sigma)
+ L2_data_RE <- data.frame(GradeID = unique(sim_data2$GradeID), 
+ 
+ U0g = L2_data_RE[, 1], U1g = L2_data_RE[, 2])
+ sim_data2 <- merge(sim_data2, L2_data_RE, by = "GradeID")
+ targetResVar <- L1_ResVar * 1.55
+ targetRandIntVar2 <- L2_RanIntVar * 1.77
+ targetRSVarx2 <- L2_x2_RSVar * 1.4
+ UFE_L1_x1_Cat2_CMC <- ref_FE(0.1, targetResVar, var(sim_data2$L1_x1_Cat2_CMC))
+ UFE_L1_x1_Cat3_CMC <- ref_FE(0.15, targetResVar, var(sim_data2$L1_x1_Cat3_CMC))
+ UFE_L1_x2_CMC <- ref_FE(0.3, targetResVar, var(sim_data2$L1_x2_Comp_Round_CMC))
+ UFE_crosslevel <- ref_FE(0.4, targetRSVarx2, var(sim_data2$L2_x2_Mean_CMC))
+ UFE_L2_x1_Cat2_CMC <- ref_FE(0.05, targetRandIntVar2, var(sim_data2$L2_x1_PropCat2_CMC))
+ UFE_L2_x1_Cat3_CMC <- ref_FE(0.08, targetRandIntVar2, var(sim_data2$L2_x1_PropCat3_CMC))
+ UFE_L2_x2_CMC <- ref_FE(0.2, targetRandIntVar2, var(sim_data2$L2_x2_Mean_CMC))
+ UFE_L2_x3_Cat2 <- ref_FE(0.02, targetRandIntVar2, var(sim_data2$L2_x3_Cat2))
+ UFE_L2_x3_Cat3 <- ref_FE(0.05, targetRandIntVar2, var(sim_data2$L2_x3_Cat3))
+ UFE_L2_x3_Cat4 <- ref_FE(0.07, targetRandIntVar2, var(sim_data2$L2_x3_Cat4))
+ UFE_quad <- ref_FE(0.3, targetRandIntVar2, var(sim_data2$L2_x2_Mean_CMCQuad))
+ sim_data2$y_CMC <- (25 + sim_data2$U0g + sim_data2$U0s + 
+ 
+ UFE_L1_x1_Cat2_CMC * sim_data2$L1_x1_Cat2_CMC + UFE_L1_x1_Cat3_CMC * 
+ 
+ sim_data2$L1_x1_Cat3_CMC + UFE_L1_x2_CMC * sim_data2$L1_x2_Comp_Round_CMC + 
+ 
+ UFE_crosslevel * sim_data2$L1_x2_Comp_Round_CMC * sim_data2$L2_x2_Mean_CMC + 
+ 
+ UFE_L2_x1_Cat2_CMC * sim_data2$L2_x1_PropCat2_CMC + UFE_L2_x1_Cat3_CMC * 
+ 
+ sim_data2$L2_x1_PropCat3_CMC + UFE_L2_x2_CMC * sim_data2$L2_x2_Mean_CMC + 
+ 
+ UFE_L2_x3_Cat2 * sim_data2$L2_x3_Cat2 + UFE_L2_x3_Cat3 * 
+ 
+ sim_data2$L2_x3_Cat3 + UFE_L2_x3_Cat4 * sim_data2$L2_x3_Cat4 + 
+ 
+ UFE_quad * sim_data2$L2_x2_Mean_CMCQuad + UFE_L2_x1_Cat2_CMC * 
+ 
+ 0.7 * sim_data2$L3_x1_PropCat2 + UFE_L2_x1_Cat3_CMC * 
+ 
+ 0.7 * sim_data2$L3_x1_PropCat3 + UFE_L2_x2_CMC * 0.7 * 
+ 
+ sim_data2$L3_x2_Mean + sim_data2$L1_x2_Comp_Round_CMC * 
+ 
+ sim_data2$U1g + rnorm(nrow(sim_data2), 0, sqrt(L1_ResVar)))
+ sim_data2$y_CMC_Round <- round(sim_data2$y_CMC)
+ l_bound_y <- 15
+ u_bound_y <- 45
+ sim_data2$y_CMC_Model <- ifelse(sim_data2$y_CMC_Round < l_bound_y, 
+ 
+ l_bound_y, ifelse(sim_data2$y_CMC_Round > u_bound_y, 
+ 
+ 
+ u_bound_y, sim_data2$y_CMC_Round))
+ return(sim_data2)
+ }
> sim_data <- data_gen(nSchools = 10, min_StuPerGrade = 25, max_StuPerGrade = 30)
> sim_names <- c("L1_x1_Comp_Cat", "L1_x2_Comp_Round", "y_CMC_Model", 
+ "L1_x1_Cat2_CMC", "L1_x1_Cat3_CMC", "L1_x2_Comp_Round_CMC", 
+ "L2_x3_Cat2", "L2_x3_Cat3", "L2_x3_Cat4", "L2_x1_PropCat2_CMC", 
+ "L2_x1_PropCat3_CMC", "L2_x2_Mean_CMC", "L2_x2_Mean_CMCQuad")
> model_names <- c("LunchBenefit", "Permissiveness", "CB", "CMCreducedprice", 
+ "CMCfreelunch", "CMCpermissiveness", "CMCgrade10th", "CMCgrade11th", 
+ "CMCgrade12th", "CMCpropreducedprice", "CMCpropfreelunch", 
+ "CMCavgpermissiveness", "CMCavgpermissivenessQuad")
> names_index <- match(sim_names, names(sim_data))
> names(sim_data)[names_index] <- model_names
> cat("\n========== DATA GENERATION COMPLETE ==========\n")

========== DATA GENERATION COMPLETE ==========
> cat("Sample size:", nrow(sim_data), "students\n")
Sample size: 1107 students
> cat("Number of schools:", length(unique(sim_data$SchoolID)), 
+ "\n")
Number of schools: 10 
> cat("Number of grades:", length(unique(sim_data$GradeID)), "\n")
Number of grades: 40 
> cat("Average students per grade:", round(nrow(sim_data)/length(unique(sim_data$GradeID))), 
+ "\n\n")
Average students per grade: 28 

> sim_data <- read.csv("Example2_SimulationData.csv")
> cat("========== STEP 1: UNCONDITIONAL MODELS ==========\n\n")
========== STEP 1: UNCONDITIONAL MODELS ==========

> cat("--- Three-Level Empty Model for Cyberbullying Risk ---\n")
--- Three-Level Empty Model for Cyberbullying Risk ---
> empty_y_3L <- lmer(formula = CB ~ 1 + (1 | GradeID) + (1 | SchoolID), 
+ control = lmerControl(optimizer = "Nelder_Mead"), REML = TRUE, 
+ data = sim_data)
> cat("\nVariance Components (for CB in Table 5):\n")

Variance Components (for CB in Table 5):
> print(data.frame(VarCorr(empty_y_3L))[, c(1, 4)])
       grp  vcov
1  GradeID  2.95
2 SchoolID  0.48
3 Residual 16.94
> cat("\nICC Estimates:\n")

ICC Estimates:
> print(icc(empty_y_3L))
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.168
  Unadjusted ICC: 0.168
> cat("\nLikelihood Ratio Tests for Random Effects:\n")

Likelihood Ratio Tests for Random Effects:
> print(ranova(empty_y_3L))
ANOVA-like table for random-effects: Single term deletions

Model:
CB ~ 1 + (1 | GradeID) + (1 | SchoolID)
               npar logLik  AIC  LRT Df Pr(>Chisq)
<none>            4  -3174 6356                   
(1 | GradeID)     3  -3215 6435 81.2  1     <2e-16
(1 | SchoolID)    3  -3174 6355  0.7  1       0.41
> cat("\n\n--- Hybrid Two-Level Empty Model ---\n")


--- Hybrid Two-Level Empty Model ---
> cat("(Including school as fixed effect to absorb L3 variance in the intercept)\n")
(Including school as fixed effect to absorb L3 variance in the intercept)
> empty_y_2L <- lmer(formula = CB ~ 1 + as.factor(SchoolID) + (1 | 
+ GradeID), control = lmerControl(optimizer = "Nelder_Mead"), 
+ REML = TRUE, data = sim_data)
> cat("\nVariance Components (for CB in Table 5):\n")

Variance Components (for CB in Table 5):
> print(data.frame(VarCorr(empty_y_2L))[, c(1, 4)])
       grp vcov
1  GradeID  2.9
2 Residual 16.9
> cat("\nICC at Level 2:\n")

ICC at Level 2:
> print(icc(empty_y_2L))
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.148
  Unadjusted ICC: 0.139
> cat("\n\n--- Unconditional Models for Predictors ---\n")


--- Unconditional Models for Predictors ---
> cat("\nSchool Lunch Benefit:\n")

School Lunch Benefit:
> empty_x1_RInt3 <- clmm(as.factor(LunchBenefit) ~ 1 + (1 | GradeID) + 
+ (1 | SchoolID), data = sim_data)
> cat("\nVariance Components (for Lunch Benefit in Table 5):\n")

Variance Components (for Lunch Benefit in Table 5):
> print(data.frame(VarCorr(empty_x1_RInt3)))
            X.Intercept. X.Intercept..1
(Intercept)          1.2           0.15
> cat("ICC L2:", round(icc(empty_x1_RInt3)[, 1], 3), "\n")
ICC L2: 0.3 
> cat("ICC L3:", data.frame(VarCorr(empty_x1_RInt3))[1, 2]/sum(data.frame(VarCorr(empty_x1_RInt3))[1, 
+ ]), "\n")
ICC L3: 0.11 
> cat("\nPermissiveness toward Cyberbullying:\n")

Permissiveness toward Cyberbullying:
> empty_x2_RInt <- lmer(formula = Permissiveness ~ 1 + (1 | GradeID) + 
+ (1 | SchoolID), control = lmerControl(optimizer = "Nelder_Mead"), 
+ REML = TRUE, data = sim_data)
> cat("\nVariance Components (for Permissiveness in Table 5):\n")

Variance Components (for Permissiveness in Table 5):
> print(data.frame(VarCorr(empty_x2_RInt))[, c(1, 4)])
       grp vcov
1  GradeID 4.35
2 SchoolID 0.91
3 Residual 9.06
> cat("ICC L2:", round(icc(empty_x2_RInt)[, 1], 3), "\n")
ICC L2: 0.37 
> cat("ICC L3:", data.frame(VarCorr(empty_x2_RInt))[2, 4]/sum(data.frame(VarCorr(empty_x2_RInt))[1:2, 
+ 4]), "\n")
ICC L3: 0.17 
> cat("\n\n========== STEP 2: BUILDING MULTILEVEL MODELS ==========\n\n")


========== STEP 2: BUILDING MULTILEVEL MODELS ==========

> cat("--- MODEL 1: Grade Level Effects ---\n")
--- MODEL 1: Grade Level Effects ---
> model1_grade <- lmer(formula = CB ~ 1 + CMCgrade10th + CMCgrade11th + 
+ CMCgrade12th + as.factor(SchoolID) + (1 | GradeID), control = lmerControl(optimizer = "Nelder_Mead"), 
+ REML = TRUE, data = sim_data)
> cat("\nFixed Effects:\n")

Fixed Effects:
> print(summary(model1_grade, ddf = "Satterthwaite")$coefficients[1:4, 
+ ])
             Estimate Std. Error df t value Pr(>|t|)
(Intercept)      36.5       0.98 27    37.3  1.3e-24
CMCgrade10th      1.3       0.77 27     1.7  1.0e-01
CMCgrade11th      1.3       0.77 27     1.6  1.1e-01
CMCgrade12th      2.3       0.77 27     3.0  5.6e-03
> cat("\nVariance Components:\n")

Variance Components:
> print(data.frame(VarCorr(model1_grade))[, c(1, 4)])
       grp vcov
1  GradeID  2.3
2 Residual 16.9
> VarComp_empty <- data.frame(VarCorr(empty_y_2L))[, c(1, 4)]
> VarComp_m1 <- data.frame(VarCorr(model1_grade))[, c(1, 4)]
> pseudoR2_L2_m1 <- (VarComp_empty[1, 2] - VarComp_m1[1, 2])/VarComp_empty[1, 
+ 2]
> cat("\nPseudo R² (Random Intercept):", round(pseudoR2_L2_m1 * 
+ 100, 2), "%\n")

Pseudo R² (Random Intercept): 20 %
> cat("\nMultivariate Wald F-test for all grade level effects:\n")

Multivariate Wald F-test for all grade level effects:
> print(contestMD(model1_grade, ddf = "Satterthwaite", L = rbind(c(0, 
+ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 1, 0, 0, 0, 
+ 0, 0, 0, 0, 0, 0, 0), c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
+ 0, 0))), digits = 4)
  Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
1  154.3   51.44     3 26.93   3.036 0.0463
> cat("\n\n--- MODEL 2: Add Level 1 Within-Grade Effects ---\n")


--- MODEL 2: Add Level 1 Within-Grade Effects ---
> model2_within <- lmer(formula = CB ~ 1 + CMCreducedprice + CMCfreelunch + 
+ CMCpermissiveness + CMCgrade10th + CMCgrade11th + CMCgrade12th + 
+ as.factor(SchoolID) + (1 | GradeID), control = lmerControl(optimizer = "Nelder_Mead"), 
+ REML = TRUE, data = sim_data)
> cat("\nFixed Effects:\n")

Fixed Effects:
> print(summary(model2_within, ddf = "Satterthwaite")$coefficients[1:7, 
+ ])
                  Estimate Std. Error   df t value Pr(>|t|)
(Intercept)          36.50       0.98   27    37.2  1.3e-24
CMCreducedprice       2.12       0.29 1064     7.3  7.8e-13
CMCfreelunch          3.04       0.36 1064     8.4  9.5e-17
CMCpermissiveness     0.25       0.04 1064     6.3  3.2e-10
CMCgrade10th          1.29       0.77   27     1.7  1.0e-01
CMCgrade11th          1.26       0.77   27     1.6  1.1e-01
CMCgrade12th          2.32       0.77   27     3.0  5.6e-03
> cat("\nVariance Components:\n")

Variance Components:
> print(data.frame(VarCorr(model2_within))[, c(1, 4)])
       grp vcov
1  GradeID  2.4
2 Residual 15.2
> VarComp_m2 <- data.frame(VarCorr(model2_within))[, c(1, 4)]
> pseudoR2_L1_m2 <- (VarComp_m1[2, 2] - VarComp_m2[2, 2])/VarComp_m1[2, 
+ 2]
> cat("\nPseudo R² (Residual Variance):", round(pseudoR2_L1_m2 * 
+ 100, 2), "%\n")

Pseudo R² (Residual Variance): 10 %
> cat("\nMultivariate Wald F-test for all level-1 effects:\n")

Multivariate Wald F-test for all level-1 effects:
> print(contestMD(model2_within, ddf = "Satterthwaite", L = rbind(c(0, 
+ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 1, 
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 0, 1, 0, 
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))), digits = 4)
  Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)
1   1908     636     3  1064   41.85 1.441e-25
> cat("\n\n--- MODEL 3: Add Level 2 Between-Grade Effects ---\n")


--- MODEL 3: Add Level 2 Between-Grade Effects ---
> model3_between <- lmer(formula = CB ~ 1 + CMCreducedprice + CMCfreelunch + 
+ CMCpermissiveness + CMCgrade10th + CMCgrade11th + CMCgrade12th + 
+ CMCpropreducedprice + CMCpropfreelunch + CMCavgpermissiveness + 
+ as.factor(SchoolID) + (1 | GradeID), control = lmerControl(optimizer = "Nelder_Mead"), 
+ REML = TRUE, data = sim_data)
> cat("\nFixed Effects:\n")

Fixed Effects:
> print(summary(model3_between, ddf = "Satterthwaite")$coefficients[1:10, 
+ ])
                     Estimate Std. Error   df t value Pr(>|t|)
(Intercept)             36.43       0.90   24   40.37  1.7e-23
CMCreducedprice          2.12       0.29 1064    7.25  7.8e-13
CMCfreelunch             3.04       0.36 1064    8.45  9.5e-17
CMCpermissiveness        0.25       0.04 1064    6.35  3.2e-10
CMCgrade10th             1.34       0.71   24    1.89  7.0e-02
CMCgrade11th             1.34       0.77   24    1.74  9.5e-02
CMCgrade12th             2.36       0.75   24    3.15  4.4e-03
CMCpropreducedprice      3.97       2.16   24    1.84  7.9e-02
CMCpropfreelunch         1.18       2.00   24    0.59  5.6e-01
CMCavgpermissiveness     0.29       0.14   24    2.16  4.1e-02
> cat("\nVariance Components:\n")

Variance Components:
> print(data.frame(VarCorr(model3_between))[, c(1, 4)])
       grp vcov
1  GradeID  1.9
2 Residual 15.2
> VarComp_m3 <- data.frame(VarCorr(model3_between))[, c(1, 4)]
> pseudoR2_L2_m3 <- (VarComp_m2[1, 2] - VarComp_m3[1, 2])/VarComp_m2[1, 
+ 2]
> cat("\nPseudo R² (Random Intercept):", round(pseudoR2_L2_m3 * 
+ 100, 2), "%\n")

Pseudo R² (Random Intercept): 21 %
> cat("\nMultivariate Wald F-test for all level-2 CM effects:\n")

Multivariate Wald F-test for all level-2 CM effects:
> print(contestMD(model3_between, ddf = "Satterthwaite", L = rbind(c(0, 
+ 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 
+ 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 
+ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0))), 
+ digits = 4)
  Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)
1  131.2   43.72     3 24.08   2.877 0.05695
> cat("\n\n--- MODEL 4: Add Random Slope for Permissiveness ---\n")


--- MODEL 4: Add Random Slope for Permissiveness ---
> model4_RS <- lmer(formula = CB ~ 1 + CMCreducedprice + CMCfreelunch + 
+ CMCpermissiveness + CMCgrade10th + CMCgrade11th + CMCgrade12th + 
+ CMCpropreducedprice + CMCpropfreelunch + CMCavgpermissiveness + 
+ as.factor(SchoolID) + (1 + CMCpermissiveness | GradeID), 
+ control = lmerControl(optimizer = "Nelder_Mead"), REML = TRUE, 
+ data = sim_data)
> cat("\nFixed Effects:\n")

Fixed Effects:
> print(summary(model4_RS, ddf = "Satterthwaite")$coefficients[1:10, 
+ ])
                     Estimate Std. Error   df t value Pr(>|t|)
(Intercept)             36.29       0.90   24   40.29  1.4e-23
CMCreducedprice          1.53       0.13 1026   12.18  5.7e-32
CMCfreelunch             2.25       0.15 1026   14.51  1.6e-43
CMCpermissiveness        0.28       0.19   39    1.45  1.5e-01
CMCgrade10th             1.40       0.70   24    1.98  5.9e-02
CMCgrade11th             1.51       0.77   24    1.97  6.1e-02
CMCgrade12th             2.45       0.75   24    3.28  3.1e-03
CMCpropreducedprice      3.75       2.16   24    1.74  9.5e-02
CMCpropfreelunch         0.74       1.99   24    0.37  7.1e-01
CMCavgpermissiveness     0.21       0.13   24    1.54  1.4e-01
> cat("\nVariance Components (including random slope):\n")

Variance Components (including random slope):
> print(data.frame(VarCorr(model4_RS))[, c(1, 4)])
       grp vcov
1  GradeID 2.39
2  GradeID 1.45
3  GradeID 0.29
4 Residual 2.70
> cat("\nLikelihood Ratio Test for Random Slope:\n")

Likelihood Ratio Test for Random Slope:
> anova_RS <- anova(model3_between, model4_RS, test = "LRT")
> print(anova_RS)
Data: sim_data
Models:
model3_between: CB ~ 1 + CMCreducedprice + CMCfreelunch + CMCpermissiveness + CMCgrade10th + CMCgrade11th + CMCgrade12th + CMCpropreducedprice + CMCpropfreelunch + CMCavgpermissiveness + as.factor(SchoolID) + (1 | GradeID)
model4_RS: CB ~ 1 + CMCreducedprice + CMCfreelunch + CMCpermissiveness + CMCgrade10th + CMCgrade11th + CMCgrade12th + CMCpropreducedprice + CMCpropfreelunch + CMCavgpermissiveness + as.factor(SchoolID) + (1 + CMCpermissiveness | GradeID)
               npar  AIC  BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model3_between   21 6232 6337  -3095      6190                    
model4_RS        23 4586 4701  -2270      4540  1650  2     <2e-16
> cat("\n\n--- MODEL 5: Full Model with Cross-Level Interaction ---\n")


--- MODEL 5: Full Model with Cross-Level Interaction ---
> model5_full <- lmer(formula = CB ~ 1 + CMCreducedprice + CMCfreelunch + 
+ CMCpermissiveness + CMCpermissiveness:CMCavgpermissiveness + 
+ CMCgrade10th + CMCgrade11th + CMCgrade12th + CMCpropreducedprice + 
+ CMCpropfreelunch + CMCavgpermissiveness + CMCavgpermissivenessQuad + 
+ as.factor(SchoolID) + (1 + CMCpermissiveness | GradeID), 
+ control = lmerControl(optimizer = "Nelder_Mead"), REML = TRUE, 
+ data = sim_data)
> cat("\nFINAL MODEL SUMMARY:\n")

FINAL MODEL SUMMARY:
> print(summary(model5_full, ddf = "Satterthwaite")$coefficients[1:11, 
+ ])
                         Estimate Std. Error   df t value Pr(>|t|)
(Intercept)                 36.08      0.809   23   44.60  7.2e-24
CMCreducedprice              1.54      0.126 1026   12.19  5.0e-32
CMCfreelunch                 2.24      0.155 1026   14.50  1.9e-43
CMCpermissiveness            0.28      0.146   38    1.94  6.0e-02
CMCgrade10th                 1.44      0.629   23    2.29  3.2e-02
CMCgrade11th                 2.06      0.718   23    2.87  8.6e-03
CMCgrade12th                 1.94      0.695   23    2.79  1.0e-02
CMCpropreducedprice          3.00      1.950   23    1.54  1.4e-01
CMCpropfreelunch             0.63      1.779   23    0.36  7.2e-01
CMCavgpermissiveness         0.35      0.123   23    2.85  9.1e-03
CMCavgpermissivenessQuad     0.18      0.066   23    2.65  1.4e-02
> cat("\n\nVariance Components:\n")


Variance Components:
> print(data.frame(VarCorr(model5_full))[, c(1, 4)])
       grp vcov
1  GradeID 1.88
2  GradeID 0.84
3  GradeID 0.16
4 Residual 2.70
> VarComp_m4 <- data.frame(VarCorr(model4_RS))[, c(1, 4)]
> VarComp_m5 <- data.frame(VarCorr(model5_full))[, c(1, 4)]
> RSVar_m4 <- VarComp_m4[2, 2]
> RSVar_m5 <- VarComp_m5[2, 2]
> pseudoR2_RS <- (RSVar_m4 - RSVar_m5)/RSVar_m4
> cat("\nPseudo R² (Random Slope Variance):", round(pseudoR2_RS * 
+ 100, 2), "%\n")

Pseudo R² (Random Slope Variance): 42 %
> cat("\nMultivariate Wald F-test for 2 new interactions:\n")

Multivariate Wald F-test for 2 new interactions:
> print(contestMD(model5_full, ddf = "Satterthwaite", L = rbind(c(0, 
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
+ 0), c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
+ 0, 0, 0, 1))), digits = 4)
  Sum Sq Mean Sq NumDF DenDF F value      Pr(>F)
1  96.79   48.39     2 28.49    17.9 0.000009212
> cat("\n\n========== FINAL MODEL INTERPRETATION ==========\n\n")


========== FINAL MODEL INTERPRETATION ==========

> cat("Key Findings:\n")
Key Findings:
> cat("1. Students with school lunch benefits have higher cyberbullying risk\n")
1. Students with school lunch benefits have higher cyberbullying risk
> cat("2. Individual permissiveness predicts higher risk\n")
2. Individual permissiveness predicts higher risk
> cat("3. Grade-level permissiveness moderates individual effect (cross-level interaction)\n")
3. Grade-level permissiveness moderates individual effect (cross-level interaction)
> cat("4. Risk increases across grade levels (9th to 12th)\n")
4. Risk increases across grade levels (9th to 12th)
> cat("5. Quadratic effect of grade permissiveness shows accelerating pattern\n\n")
5. Quadratic effect of grade permissiveness shows accelerating pattern

> cat("Variance Explained:\n")
Variance Explained:
> cat("- Level 1 (Residual):", round(pseudoR2_L1_m2 * 100, 2), 
+ "%\n")
- Level 1 (Residual): 10 %
> cat("- Level 2 (Random Intercept):", round((1 - VarComp_m3[1, 
+ 2]/VarComp_empty[1, 2]) * 100, 2), "%\n")
- Level 2 (Random Intercept): 35 %
> cat("- Random Slope:", round(pseudoR2_RS * 100, 2), "%\n\n")
- Random Slope: 42 %

> cat("\n========== SESSION INFO ==========\n\n")

========== SESSION INFO ==========

> print(sessionInfo())
R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                           LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] TeachingDemos_2.13 ordinal_2025.12-29 performance_0.16.0 lmerTest_3.2-1     lme4_2.0-1         Matrix_1.7-4      
[7] pacman_0.5.1      

loaded via a namespace (and not attached):
 [1] gtable_0.3.6        ucminf_1.2.2        dplyr_1.2.0         compiler_4.5.2      tidyselect_1.2.1   
 [6] Rcpp_1.1.1          splines_4.5.2       scales_1.4.0        boot_1.3-32         lattice_0.22-7     
[11] ggplot2_4.0.1       R6_2.6.1            generics_0.1.4      rbibutils_2.4.1     MASS_7.3-65        
[16] tibble_3.3.0        nloptr_2.2.1        insight_1.4.6       minqa_1.2.8         pillar_1.11.1      
[21] RColorBrewer_1.1-3  rlang_1.1.7         S7_0.2.1            cli_3.6.5           magrittr_2.0.4     
[26] Rdpack_2.6.6        grid_4.5.2          rstudioapi_0.18.0   lifecycle_1.0.5     nlme_3.1-168       
[31] reformulas_0.4.4    vctrs_0.7.1         glue_1.8.0          farver_2.1.2        numDeriv_2016.8-1.1
[36] tools_4.5.2         pkgconfig_2.0.3    
