Linear Mixed Model (LMM)

# Linear Mixed Model: Students in Schools
# Study Hours → Math Achievement (accounting for school clustering)
# Demonstrates random intercepts model

library(lme4)         # LMM fitting
library(lmerTest)     # p-values for lmer
library(performance)  # ICC, R²
library(sjPlot)       # Visualization
library(ggplot2)
library(dplyr)

# Simulate realistic hierarchical data
set.seed(2025)
n_schools <- 30
n_students_per_school <- 20
n_total <- n_schools * n_students_per_school  # 600 students

# Level-2 (school-level) data
school_effects <- data.frame(
  school_id = 1:n_schools,
  school_intercept = rnorm(n_schools, mean=0, sd=8)  # School random effects, τ=8
)

# Level-1 (student-level) data
data <- expand.grid(
  school_id = 1:n_schools,
  student_id = 1:n_students_per_school
) %>%
  left_join(school_effects, by="school_id") %>%
  mutate(
    study_hours = rnorm(n_total, mean=8, sd=4),
    study_hours = pmax(0, pmin(20, study_hours)),  # Constrain to 0-20
    # True model: Math = 50 + 2*study_hours + school_intercept + error
    # Fixed effect of study: β=2 points per hour
    # School random intercepts: SD=8
    # Residual error: SD=10
    math_score = 50 + 2*study_hours + school_intercept + rnorm(n_total, mean=0, sd=10),
    math_score = pmax(0, pmin(100, math_score))  # Constrain to 0-100
  )

# Check data structure
cat("=== Data Structure ===")
cat("\nTotal students:", nrow(data))
cat("\nNumber of schools:", length(unique(data$school_id)))
cat("\nStudents per school:", table(data$school_id)[1], "\n")

# Descriptive statistics
cat("\n=== Descriptive Statistics ===")
print(summary(data[, c("study_hours", "math_score")]))

cat("\n=== School-level means ===")
school_means <- data %>%
  group_by(school_id) %>%
  summarise(
    mean_math = mean(math_score),
    mean_study = mean(study_hours),
    n = n()
  )
print(head(school_means))
cat("\nMath score range across schools: [", 
    round(min(school_means$mean_math), 1), ", ", 
    round(max(school_means$mean_math), 1), "]\n", sep="")

# Visualize clustering
ggplot(data, aes(x=study_hours, y=math_score, group=school_id, color=factor(school_id))) +
  geom_point(alpha=0.4, size=1) +
  geom_smooth(method="lm", se=FALSE, size=0.5, alpha=0.6) +
  labs(title="Math Achievement by Study Hours(by School)",
       subtitle="Each line = one school; demonstrates clustering",
       x="Study Hours per Week", y="Math Test Score(0-100)") +
  theme_classic() +
  theme(legend.position="none")

# === STEP 1: Calculate ICC (Intraclass Correlation) ===
# ICC = proportion of total variance due to clustering
# Fit null model (intercept-only) to estimate variance components
null_model <- lmer(math_score ~ 1 + (1|school_id), data=data, REML=TRUE)
summary(null_model)

# Extract variance components
variance_components <- as.data.frame(VarCorr(null_model))
tau_squared <- variance_components$vcov[1]  # Between-school variance
sigma_squared <- sigma(null_model)^2         # Within-school variance

ICC <- tau_squared / (tau_squared + sigma_squared)
cat("\n=== Intraclass Correlation(ICC) ===")
cat("\nBetween-school variance(τ²):", round(tau_squared, 2))
cat("\nWithin-school variance(σ²):", round(sigma_squared, 2))
cat("\nICC:", round(ICC, 3))
cat("\n\nInterpretation: ICC =", round(ICC, 3), 
    "means", round(ICC*100, 1), "% of variance in math scores is due to school clustering.")
cat("\nThis indicates", 
    ifelse(ICC > 0.10, "substantial", "modest"), 
    "clustering; LMM is appropriate(ICC >", 
    ifelse(ICC > 0.05, "0.05", "0.02"), ").\n")

# Design effect
avg_cluster_size <- mean(table(data$school_id))
design_effect <- 1 + (avg_cluster_size - 1) * ICC
cat("\nDesign effect:", round(design_effect, 2))
cat("\nEffective sample size:", round(n_total / design_effect, 0), 
    "(clustering reduces effective n from", n_total, "to ~", round(n_total/design_effect, 0), ")\n")

# === STEP 2: Fit Linear Mixed Model (Random Intercepts) ===
model <- lmer(math_score ~ study_hours + (1|school_id), data=data, REML=TRUE)
summary(model)

cat("\n=== Fixed Effects Interpretation ===")
fixed_effects <- fixef(model)
cat("\nIntercept:", round(fixed_effects[1], 2), 
    "(predicted math score when study_hours=0, for average school)")
cat("\nStudy hours:", round(fixed_effects[2], 2), 
    "(each additional hour → +", round(fixed_effects[2], 2), "points in math score, accounting for school clustering)\n")

# Confidence intervals
CI <- confint(model, method="Wald", level=0.95)
cat("\n=== 95% Confidence Intervals ===")
print(CI)

# === STEP 3: Check Assumptions ===

cat("\n=== ASSUMPTION CHECKS ===")

# 1. Normality of level-1 residuals
resid_level1 <- residuals(model)
cat("\n1. Normality of Level-1 Residuals")
shapiro_test <- shapiro.test(sample(resid_level1, min(5000, length(resid_level1))))
cat("\n   Shapiro-Wilk test: W =", round(shapiro_test$statistic, 4), ", p =", round(shapiro_test$p.value, 4))
cat("\n  ", ifelse(shapiro_test$p.value > 0.05, "✓ Normality assumption met", "⚠ Mild deviation(acceptable with large n)"))

# Q-Q plot
qqnorm(resid_level1, main="Q-Q Plot: Level-1 Residuals")
qqline(resid_level1, col="red", lwd=2)

# 2. Normality of random effects (level-2)
ranef_school <- ranef(model)$school_id[[1]]
cat("\n\n2. Normality of Random Effects(School Intercepts)")
shapiro_test_re <- shapiro.test(ranef_school)
cat("\n   Shapiro-Wilk test: W =", round(shapiro_test_re$statistic, 4), ", p =", round(shapiro_test_re$p.value, 4))
cat("\n  ", ifelse(shapiro_test_re$p.value > 0.05, "✓ Normality of random effects met", "⚠ Check for outlier schools"))

qqnorm(ranef_school, main="Q-Q Plot: Random Effects(School Intercepts)")
qqline(ranef_school, col="red", lwd=2)

# Caterpillar plot (outlier schools)
library(lattice)
dotplot(ranef(model), main="Caterpillar Plot: Random Intercepts by School")
cat("\n   Check caterpillar plot for schools with non-overlapping confidence intervals(outliers)\n")

# 3. Homoscedasticity: Residuals vs. Fitted
fitted_vals <- fitted(model)
plot(fitted_vals, resid_level1,
     xlab="Fitted Values", ylab="Residuals",
     main="Residuals vs. Fitted Values")
abline(h=0, col="red", lwd=2, lty=2)
cat("\n3. Homoscedasticity: Check for funnel shape(variance increasing/decreasing with fitted values)")
cat("\n   If funnel present → consider log transformation or variance weights\n")

# 4. Independence: Residuals by cluster
ggplot(data.frame(school_id=data$school_id, residuals=resid_level1), 
       aes(x=factor(school_id), y=residuals)) +
  geom_boxplot() +
  geom_hline(yintercept=0, color="red", linetype="dashed") +
  labs(title="Residuals by School(Check for Systematic Patterns)",
       x="School ID", y="Residuals") +
  theme_classic() +
  theme(axis.text.x = element_blank())
cat("\n4. Independence: Check for systematic patterns in residuals by cluster")
cat("\n   If patterns present → may need additional random effects or covariates\n")

# 5. Multicollinearity (for fixed effects)
library(car)
# VIF from auxiliary linear model (LMM doesn't have built-in VIF)
lm_auxiliary <- lm(math_score ~ study_hours, data=data)
if(length(coef(lm_auxiliary)) > 2) {
  vif_vals <- vif(lm_auxiliary)
  cat("\n5. Multicollinearity(VIF):")
  print(vif_vals)
  cat("   All VIF < 5 → No multicollinearity\n")
} else {
  cat("\n5. Multicollinearity: Only one predictor, no multicollinearity possible\n")
}

# 6. Sample size adequacy
cat("\n6. Sample Size Adequacy")
cat("\n   Number of clusters(schools):", n_schools)
cat("\n   Minimum recommended: 30 (Maas & Hox, 2005)")
cat("\n  ", ifelse(n_schools >= 30, "✓ Adequate number of clusters", "⚠ Fewer than 30 clusters; interpret level-2 effects cautiously"))
cat("\n   Average cluster size:", round(avg_cluster_size, 1))
cat("\n   Minimum recommended: 5-10 per cluster")
cat("\n  ", ifelse(avg_cluster_size >= 5, "✓ Adequate cluster size\n", "⚠ Small cluster size\n"))

# === STEP 4: Model Comparison (Nested Models) ===

cat("\n=== MODEL COMPARISON ===")

# Null model (intercept only, random intercepts)
null_model <- lmer(math_score ~ 1 + (1|school_id), data=data, REML=FALSE)
# Full model (with predictor, random intercepts) 
full_model <- lmer(math_score ~ study_hours + (1|school_id), data=data, REML=FALSE)

# Likelihood ratio test
lrt <- anova(null_model, full_model)
cat("\nLikelihood Ratio Test(LRT): Does adding study_hours improve fit?\n")
print(lrt)
cat("\nInterpretation: χ²(", lrt$Df[2], ") = ", round(lrt$Chisq[2], 2), 
    ", p ", ifelse(lrt$`Pr(>Chisq)`[2] < .001, "< .001", paste("=", round(lrt$`Pr(>Chisq)`[2], 3))),
    "\n", sep="")
cat(ifelse(lrt$`Pr(>Chisq)`[2] < .05, 
           "✓ Study hours significantly improves model fit.\n",
           "Model with study_hours does not improve fit significantly.\n"))

# AIC/BIC comparison
cat("\nAIC/BIC Comparison(lower is better):\n")
cat("Null model - AIC:", round(AIC(null_model), 1), ", BIC:", round(BIC(null_model), 1), "\n")
cat("Full model - AIC:", round(AIC(full_model), 1), ", BIC:", round(BIC(full_model), 1), "\n")
cat("Δ AIC:", round(AIC(null_model) - AIC(full_model), 1), 
    "(full model is", round(AIC(null_model) - AIC(full_model), 1), "points better)\n")

# === STEP 5: Effect Sizes (R²) ===

library(performance)
r2_vals <- r2(model)
cat("\n=== Effect Sizes(Variance Explained) ===")
cat("\nMarginal R² (fixed effects only):", round(r2_vals$R2_marginal, 3))
cat("\n  Interpretation:", round(r2_vals$R2_marginal*100, 1), "% of variance explained by study_hours alone\n")
cat("\nConditional R² (fixed + random effects):", round(r2_vals$R2_conditional, 3))
cat("\n  Interpretation:", round(r2_vals$R2_conditional*100, 1), "% of variance explained by study_hours + school clustering\n")

cat("\nDifference(R²_conditional - R²_marginal):", round(r2_vals$R2_conditional - r2_vals$R2_marginal, 3))
cat("\n  Interpretation: School clustering accounts for", 
    round((r2_vals$R2_conditional - r2_vals$R2_marginal)*100, 1), "% of variance\n")

# === STEP 6: Prediction Example ===

cat("\n=== PREDICTION EXAMPLES ===")

# Predict for new student in average school
new_student_avg_school <- data.frame(study_hours = 10, school_id = NA)
pred_avg <- predict(model, newdata=new_student_avg_school, re.form=NA, allow.new.levels=TRUE)
cat("\nStudent studying 10 hours/week in AVERAGE school(population-level):")
cat("\nPredicted math score:", round(pred_avg, 1), "\n")

# Predict for new student in specific existing school (e.g., school 5)
new_student_school5 <- data.frame(study_hours = 10, school_id = 5)
pred_school5 <- predict(model, newdata=new_student_school5, re.form=NULL)
cat("\nStudent studying 10 hours/week in School 5 (school-specific):")
cat("\nPredicted math score:", round(pred_school5, 1))
cat("\nSchool 5 random intercept:", round(ranef(model)$school_id[5, 1], 2))
cat("\n(School 5 is", 
    ifelse(ranef(model)$school_id[5,1] > 0, "above", "below"), 
    "average by", abs(round(ranef(model)$school_id[5,1], 1)), "points)\n")

# Compare predictions across schools
pred_range <- data.frame(
  study_hours = 10,
  school_id = 1:n_schools
)
pred_range$predicted <- predict(model, newdata=pred_range, re.form=NULL)
cat("\nPredicted math scores for 10-hour student across all schools:")
cat("\nRange: [", round(min(pred_range$predicted), 1), ", ", round(max(pred_range$predicted), 1), "]")
cat("\nSD of predictions:", round(sd(pred_range$predicted), 1), "(reflects school-level variability)\n")

# === STEP 7: Contrast with Naive OLS (Ignoring Clustering) ===

cat("\n=== COMPARISON: LMM vs. Naive OLS(Ignoring Clustering) ===")

ols_model <- lm(math_score ~ study_hours, data=data)
cat("\nNaive OLS(ignoring school clustering):")
cat("\nStudy hours coefficient:", round(coef(ols_model)[2], 3))
cat("\nStandard error:", round(summary(ols_model)$coefficients[2, 2], 3))
cat("\np-value:", format.pval(summary(ols_model)$coefficients[2, 4], digits=3))

cat("\n\nLMM(accounting for school clustering):")
lmm_summary <- summary(model)
cat("\nStudy hours coefficient:", round(fixef(model)[2], 3))
cat("\nStandard error:", round(lmm_summary$coefficients[2, 2], 3))
cat("\np-value:", format.pval(lmm_summary$coefficients[2, 5], digits=3))

se_ratio <- summary(ols_model)$coefficients[2, 2] / lmm_summary$coefficients[2, 2]
cat("\n\nSE ratio(OLS/LMM):", round(se_ratio, 3))
cat("\nInterpretation: OLS standard error is", 
    ifelse(se_ratio < 1, "SMALLER", "comparable to"), 
    "LMM SE.")
cat("\nOLS underestimates uncertainty by ignoring clustering(inflated Type I error).")
cat("\nLMM correctly accounts for clustering → more accurate inference.\n")

# === APA-Style Reporting ===

cat("\n\n=== APA-STYLE REPORT ===")
cat("\nA linear mixed model(LMM) was conducted to examine the effect of study hours\n")
cat("on math achievement, accounting for student clustering within schools. The sample\n")
cat("included", n_total, "students nested in", n_schools, "schools(average", round(avg_cluster_size, 0), "students per school).\n")
cat("The intraclass correlation(ICC) was", round(ICC, 3), ", indicating", round(ICC*100, 1), "% of variance\n")
cat("in math scores was due to school-level clustering, supporting the use of LMM.\n\n")

cat("Assumptions were examined: Level-1 residuals were approximately normally distributed\n")
cat("(Shapiro-Wilk p", ifelse(shapiro_test$p.value > .05, ">", "<"), ".05); ")
cat("random effects(school intercepts) were normally distributed(Shapiro-Wilk p", 
    ifelse(shapiro_test_re$p.value > .05, ">", "="), ".05); ")
cat("residual plots indicated homoscedasticity; no outlier schools were identified in\n")
cat("caterpillar plots; sample size was adequate(", n_schools, "≥30 clusters).\n\n")

cat("A random intercepts model was fitted with study hours as a fixed effect and random\n")
cat("intercepts for schools. Study hours significantly predicted math achievement(β=",
    round(fixef(model)[2], 2), ", SE=", round(lmm_summary$coefficients[2,2], 2),
    ", 95% CI [", round(CI[3,1], 2), ",", round(CI[3,2], 2), "], ")
cat("t(", round(lmm_summary$coefficients[2,3], 0), ")=", round(lmm_summary$coefficients[2,4], 2),
    ", p<.001): each additional study hour per week was associated with a", 
    round(fixef(model)[2], 2), "-point\n")
cat("increase in math scores, accounting for school clustering. The model explained\n")
cat(round(r2_vals$R2_marginal*100, 1), "% of variance via study hours(marginal R²) and", 
    round(r2_vals$R2_conditional*100, 1), "% including\n")
cat("school clustering(conditional R²). Likelihood ratio test confirmed study hours\n")
cat("significantly improved model fit(χ²(1)=", round(lrt$Chisq[2], 2), ", p<.001).\n\n")

cat("Compared to naive OLS regression(which ignores clustering), LMM provided more\n")
cat("accurate standard errors(SE_LMM=", round(lmm_summary$coefficients[2,2], 3), 
    "vs. SE_OLS=", round(summary(ols_model)$coefficients[2,2], 3), "),")
cat("\navoiding inflated Type I error rates. These findings support the importance of study\n")
cat("time for academic achievement while appropriately accounting for school-level clustering.\n")