Prestasi matematika siswa antarsekolah: model multilevel motivasi dan rasio guru-siswa

# Contoh kode R: Multilevel Regression (lme4)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(lme4)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Data sintetis: individu (level 1) tersarang dalam kelompok (level 2) ----
n <- 200
J <- sample(20:36, 1)                                         # banyak kelompok
kode <- sprintf("%s %02d", "Group", 1:J)
uk <- as.numeric(rmultinom(1, n - 3 * J, rep(1, J))) + 3      # ukuran kelompok tidak sama (minimal 3)
kelompok <- rep(kode, uk)
idx <- rep(1:J, uk)
z2g <- rnorm(J)                                               # prediktor level kelompok
u0 <- rnorm(J, 0, runif(1, 0.4, 0.7))                        # efek acak intersep
u1 <- rnorm(J, 0, runif(1, 0.12, 0.3))                        # efek acak slope X1
z1 <- rnorm(n, 0.3 * z2g[idx], 1)
b1 <- runif(1, 0.35, 0.6)
b2 <- runif(1, 0.45, 0.75) * sample(c(-1, 1), 1)
y <- b1 * z1 + b2 * z2g[idx] + u0[idx] + u1[idx] * z1 + rnorm(n, 0, runif(1, 0.55, 0.8))
d <- data.frame(kelompok = kelompok,
                learning_motivation = round(50 + 10 * z1, 3),
                students_per_teacher = round(30 + 6 * z2g[idx], 3),
                student_mathematics_score = round(100 + 15 * y, 3))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: regresi multilevel (lme4::lmer) ----
d$kelompok <- factor(d$kelompok)
d$x1_pusat <- d$learning_motivation - mean(d$learning_motivation)                       # X1 dipusatkan pada rata-rata keseluruhan
# 1) Model null untuk ICC
m0 <- lmer(student_mathematics_score ~ 1 + (1 | kelompok), data = d, REML = TRUE)
vc0 <- as.data.frame(VarCorr(m0))
icc <- vc0$vcov[1] / sum(vc0$vcov)
cat(sprintf("ICC model null = %.3f\n", icc))
# 2) Random intercept, lalu coba random slope X1 (cadangan bila singular/gagal)
m1 <- lmer(student_mathematics_score ~ x1_pusat + students_per_teacher + (1 | kelompok), data = d, REML = TRUE)
m2 <- tryCatch(lmer(student_mathematics_score ~ x1_pusat + students_per_teacher + (1 + x1_pusat | kelompok), data = d, REML = TRUE,
                    control = lmerControl(calc.derivs = FALSE)), error = function(e) NULL)
p_slope <- NA
if (!is.null(m2) && !isSingular(m2)) {
  a12 <- suppressMessages(anova(m1, m2, refit = FALSE))          # LR uji varians slope acak
  print(a12)
  p_slope <- a12[2, "Pr(>Chisq)"]
}
fit <- if (is.finite(p_slope) && p_slope < 0.05) m2 else m1
label_model <- if (identical(fit, m2)) "random intercept + random slope" else "random intercept"
print(summary(fit))
ct <- coef(summary(fit))
pval <- unname(2 * pnorm(abs(ct[, "t value"]), lower.tail = FALSE))  # p-value aproksimasi normal (Wald z)
print(cbind(ct, "p (Wald z)" = round(pval, 4)))
# 3) Uji LR efek acak: model terpilih vs OLS biasa (ML), p dibagi dua (uji di batas)
ols <- lm(student_mathematics_score ~ x1_pusat + students_per_teacher, data = d)
fit_ml <- refitML(fit)
lr_stat <- max(0, as.numeric(2 * (logLik(fit_ml) - logLik(ols))))
df_lr <- attr(logLik(fit_ml), "df") - attr(logLik(ols), "df")
p_lr <- 0.5 * pchisq(lr_stat, df_lr - 1, lower.tail = FALSE) + 0.5 * pchisq(lr_stat, df_lr, lower.tail = FALSE)
cat(sprintf("Model terpilih: %s; LR vs OLS = %.2f, p = %.4g\n", label_model, lr_stat, p_lr))


# ---- Grafik ----
# Garis regresi tiap kelompok (spaghetti) dan garis efek tetap
rentang <- do.call(rbind, lapply(split(d, d$kelompok), function(g)
  data.frame(kelompok = g$kelompok[1], x1_pusat = range(g$x1_pusat), students_per_teacher = g$students_per_teacher[1])))
rentang$prediksi <- predict(fit, newdata = rentang)
rentang$x <- rentang$x1_pusat + mean(d$learning_motivation)
umum <- data.frame(x1_pusat = range(d$x1_pusat), students_per_teacher = mean(d$students_per_teacher))
umum$prediksi <- predict(fit, newdata = umum, re.form = NA)
umum$x <- umum$x1_pusat + mean(d$learning_motivation)
p1 <- ggplot() +
  geom_point(data = d, aes(learning_motivation, student_mathematics_score), colour = warna$redup, alpha = 0.25, size = 1.1) +
  geom_line(data = rentang, aes(x, prediksi, group = kelompok), colour = pal[1], alpha = 0.55, linewidth = 0.6) +
  geom_line(data = umum, aes(x, prediksi), colour = pal[2], linewidth = 1.6) +
  labs(title = sprintf("%s-specific lines and fixed-effect line", substr("Group", 1, 18)), x = "learning_motivation", y = "student_mathematics_score")
# Caterpillar plot intersep acak dengan CI 95%
re <- ranef(fit, condVar = TRUE)$kelompok
se_re <- sqrt(attr(ranef(fit, condVar = TRUE)$kelompok, "postVar")[1, 1, ])
cat_df <- data.frame(kelompok = rownames(re), b = re[, 1], se = se_re)
cat_df <- cat_df[order(cat_df$b), ]
cat_df$urut <- seq_len(nrow(cat_df))
cat_df$sig <- ifelse(cat_df$b - 1.96 * cat_df$se > 0 | cat_df$b + 1.96 * cat_df$se < 0, "CI excludes 0", "CI includes 0")
p2 <- ggplot(cat_df, aes(b, urut, colour = sig)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_segment(aes(x = b - 1.96 * se, xend = b + 1.96 * se, yend = urut), linewidth = 0.6) +
  geom_point(size = 1.6) +
  scale_colour_manual(values = c("CI excludes 0" = pal[2], "CI includes 0" = pal[3]), name = NULL) +
  scale_y_continuous(breaks = NULL) +
  labs(title = "Random intercepts (caterpillar, 95% CI)", x = "Deviation from fixed intercept",
       y = sprintf("%s (sorted)", substr("Group", 1, 20)))
# Rata-rata kelompok terhadap prediktor level kelompok
rk <- aggregate(cbind(student_mathematics_score, students_per_teacher) ~ kelompok, data = d, FUN = mean)
rk$n <- as.numeric(table(d$kelompok)[as.character(rk$kelompok)])
p3 <- ggplot(rk, aes(students_per_teacher, student_mathematics_score)) +
  geom_hline(yintercept = mean(d$student_mathematics_score), linetype = "dotted", colour = warna$redup) +
  geom_smooth(method = "lm", formula = y ~ x, colour = pal[4], fill = pal[4], alpha = 0.15, linewidth = 0.9) +
  geom_point(aes(size = n), colour = pal[5], alpha = 0.8) +
  scale_size_continuous(range = c(1.2, 4), guide = "none") +
  labs(title = sprintf("%s means vs level-2 predictor", substr("Group", 1, 18)), x = "students_per_teacher", y = "Mean student_mathematics_score")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else p1 / (p2 | p3)

print(p_hasil)