Modul 2 — MANOVA, MANCOVA & ANCOVA ASSUMTION

pkgs <- c("tidyverse", "biotools", "car", "heplots",
          "psych", "ggplot2", "knitr", "kableExtra",
          "mvnormtest", "MVN", "emmeans", "effectsize")

new_pkgs <- pkgs[!pkgs %in% rownames(installed.packages())]
if (length(new_pkgs)) {
  install.packages(new_pkgs, repos = "https://cran.r-project.org")
}

suppressWarnings(suppressPackageStartupMessages({
  library(tidyverse)
  library(biotools)
  library(car)
  library(heplots)
  library(psych)
  library(ggplot2)
  library(knitr)
  library(kableExtra)
  library(MVN)
  library(emmeans)
  library(effectsize)
}))

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
                      fig.align = "center", fig.width = 9,
                      fig.height = 5, dpi = 150)

theme_set(
  theme_minimal(base_size = 13) +
  theme(plot.title   = element_text(face = "bold", color = "#1a5276"),
        plot.subtitle = element_text(color = "#566573"),
        axis.title   = element_text(color = "#2c3e50"),
        strip.text   = element_text(face = "bold"))
)

if (!file.exists('education_intervention_10000.csv')) {
  download.file(
    'https://drive.google.com/uc?id=1cJlEWjZHdfqlM0r43GFbyQXf_MN0_zJk',
    'education_intervention_10000.csv', mode = 'wb'
  )
}

df_raw <- read.csv("education_intervention_10000.csv", sep = ",", stringsAsFactors = TRUE)
cat("Dimensi data mentah:", dim(df_raw), "\n")

## Dimensi data mentah: 10000 16

head(df_raw)

##   student_id school_id grade_level intervention ses_index attendance_rate
## 1   93207890       110           9      blended     1.379           0.996
## 2   61581758       128          11      blended     2.037           0.930
## 3   75844861       158           9     tutoring     0.403           0.933
## 4   57564202       153           9      flipped    -0.442           0.924
## 5   57864707       100          12     tutoring    -0.572           0.971
## 6   78728520       100          11      blended    -1.128           0.875
##   study_hours_wk teacher_experience_years class_size baseline_math
## 1           5.93                      0.0         28          76.9
## 2           2.52                      6.7         32          78.9
## 3           5.32                      8.3         23          69.4
## 4           5.26                      9.2         30          55.9
## 5           7.47                     10.6         21          55.3
## 6           4.12                     15.7         25          59.4
##   baseline_reading baseline_motivation_1_10 post_math post_reading
## 1             65.3                      7.7      87.8         75.1
## 2             82.3                      7.1      88.9         89.0
## 3             83.9                      5.8      84.8         95.2
## 4             61.8                      6.7      72.7         68.3
## 5             66.4                      4.6      75.7         83.4
## 6             67.2                      7.4      69.1         78.9
##   test_anxiety_1_10 passed
## 1               2.6      1
## 2               4.3      1
## 3               3.8      1
## 4               5.6      1
## 5               7.6      1
## 6               6.8      0

DVS <- c("post_math", "post_reading")         
COV <- c("ses_index", "attendance_rate",       
         "study_hours_wk", "teacher_experience_years",
         "class_size", "baseline_math",
         "baseline_reading", "baseline_motivation_1_10",
         "test_anxiety_1_10")
IV  <- c("intervention", "grade_level")        

df <- df_raw %>%
  dplyr::select(all_of(c(IV, DVS, COV))) %>%
  drop_na() %>%
  mutate(intervention = as.factor(intervention),
         grade_level  = as.factor(grade_level))

cat("Dimensi setelah seleksi & drop_na:", dim(df), "")

## Dimensi setelah seleksi & drop_na: 10000 13

0.1 Statistika Deskriptif

summary(df)

##    intervention  grade_level   post_math      post_reading   
##  blended :2229   9 :2766     Min.   : 29.7   Min.   : 35.90  
##  control :2992   10:2547     1st Qu.: 71.9   1st Qu.: 72.50  
##  flipped :2207   11:2448     Median : 81.4   Median : 81.40  
##  tutoring:2572   12:2239     Mean   : 80.9   Mean   : 80.87  
##                              3rd Qu.: 91.5   3rd Qu.: 90.20  
##                              Max.   :100.0   Max.   :100.00  
##    ses_index         attendance_rate  study_hours_wk   teacher_experience_years
##  Min.   :-2.500000   Min.   :0.6670   Min.   : 0.020   Min.   : 0.000          
##  1st Qu.:-0.678250   1st Qu.:0.8730   1st Qu.: 2.250   1st Qu.: 5.600          
##  Median : 0.001000   Median :0.9160   Median : 3.910   Median : 9.000          
##  Mean   :-0.002791   Mean   :0.9129   Mean   : 4.576   Mean   : 9.057          
##  3rd Qu.: 0.671000   3rd Qu.:0.9570   3rd Qu.: 6.140   3rd Qu.:12.400          
##  Max.   : 2.500000   Max.   :1.0000   Max.   :22.850   Max.   :28.700          
##    class_size    baseline_math    baseline_reading baseline_motivation_1_10
##  Min.   :10.00   Min.   : 20.00   Min.   : 23.30   Min.   : 1.000          
##  1st Qu.:22.00   1st Qu.: 61.00   1st Qu.: 63.60   1st Qu.: 5.300          
##  Median :26.00   Median : 70.00   Median : 71.80   Median : 6.500          
##  Mean   :26.06   Mean   : 70.03   Mean   : 71.77   Mean   : 6.456          
##  3rd Qu.:30.00   3rd Qu.: 79.20   3rd Qu.: 79.92   3rd Qu.: 7.600          
##  Max.   :40.00   Max.   :100.00   Max.   :100.00   Max.   :10.000          
##  test_anxiety_1_10
##  Min.   : 1.000   
##  1st Qu.: 3.700   
##  Median : 4.800   
##  Mean   : 4.764   
##  3rd Qu.: 5.900   
##  Max.   :10.000

0.2 Distribusi Dependent Variable

hist(df$post_math, 
        main = "Distribusi Nilai Matematika", 
        xlab = "Kategori", 
        ylab = "Frekuensi", 
        col = "lightblue",
        border = "white")

hist(df$post_reading, 
        main = "Distribusi Nilai Membaca", 
        xlab = "Kategori", 
        ylab = "Frekuensi", 
        col = "lightblue",
        border = "white")

# Box Plot

0.3 Deteksi Outlier Multivariat (Mahalanobis)

MD        <- mahalanobis(df[, DVS], colMeans(df[, DVS]), cov(df[, DVS]))
threshold <- qchisq(0.95, df = length(DVS))
df_clean  <- df[MD < threshold, ]

cat("Outlier dibuang:", nrow(df) - nrow(df_clean), "baris\n")

## Outlier dibuang: 365 baris

cat("Sisa observasi :", nrow(df_clean), "\n")

## Sisa observasi : 9635

# Standarisasi Z-Score semua variabel kontinu
df_clean <- df_clean %>%
  mutate(across(all_of(c(DVS, COV)), ~ scale(.)[, 1]))

df <- as.data.frame(df_clean)

1 Uji Asumsi

1.1 Asumsi 1 — Normalitas Multivariat

mardia_res <- psych::mardia(df[, DVS], plot = FALSE)
p_skew     <- mardia_res$p.skew
p_kurt     <- mardia_res$p.kurt

tibble(
  Komponen  = c("Mardia Skewness", "Mardia Kurtosis"),
  Statistik = c(round(mardia_res$b1p, 4), round(mardia_res$b2p, 4)),
  `p-value` = c(format(p_skew, scientific = TRUE, digits = 3),
                format(p_kurt, scientific = TRUE, digits = 3)),
  Status    = c(ifelse(p_skew >= 0.05, "✔ Normal", "✘ Tidak Normal"),
                ifelse(p_kurt >= 0.05, "✔ Normal", "✘ Tidak Normal"))
) %>%
  kable(caption = "Mardia's Test — Normalitas Multivariat") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  column_spec(4, bold = TRUE,
              color = ifelse(c(p_skew, p_kurt) >= 0.05, "#1e8449", "#c0392b"))

Mardia’s Test — Normalitas Multivariat

Komponen	Statistik	p-value	Status
Mardia Skewness	0.0786	2.51e-26	✘ Tidak Normal
Mardia Kurtosis	6.3236	0e+00	✘ Tidak Normal

cat("Mardia Skewness p-value:", p_skew, "\n")

## Mardia Skewness p-value: 2.51487e-26

cat("Mardia Kurtosis p-value:", p_kurt, "\n")

## Mardia Kurtosis p-value: 0

P-value dari keduanya tidak memenuhi p-value ≥ 0.05, Namun pada teori Central Limit Theorem menjelaskan bahwa data akan mendekati hasil distribusi normal saat jumlah data > 30 * variabel

jadi proses pengujian akan dilanjutkan

dim(df)

## [1] 9635   13

Karena ada 9000+ data maka distribusinya akan mendekati normal

# Chi-Square Q-Q Plot (visualisasi normalitas multivariat)
MD_clean <- mahalanobis(df[, DVS], colMeans(df[, DVS]), cov(df[, DVS]))
chi2_q   <- qchisq(ppoints(nrow(df)), df = length(DVS))

par(mar = c(4, 4, 3, 1))
qqplot(chi2_q, sort(MD_clean),
       main  = "Chi-Square Q-Q Plot — Normalitas Multivariat",
       xlab  = "Theoretical Chi-Square Quantiles",
       ylab  = "Mahalanobis Distance (sorted)",
       pch   = 16, col = "#2e86c1", cex = 0.5)
abline(0, 1, col = "#c0392b", lwd = 2)
legend("topleft", legend = c("Data", "Garis Referensi"),
       col = c("#2e86c1", "#c0392b"), pch = c(16, NA),
       lty = c(NA, 1), lwd = c(NA, 2), bty = "n")

Pada distribusi multivariat QQ plot seluruh fitur cenderung left skew

1.2 Asumsi 2 — Homogenitas Matriks Kovarians

faktor <- interaction(df$intervention, df$grade_level)
bm     <- biotools::boxM(df[, DVS], faktor)

print(bm)

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  df[, DVS]
## Chi-Sq (approx.) = 41.692, df = 45, p-value = 0.6129

cat("Box's M p-value:", round(bm$p.value, 4), "")

## Box's M p-value: 0.6129

if (bm$p.value > 0.05) {
  cat(" Gagal tolak H0 Matriks kovarians HOMOGEN ")
} else if (bm$p.value < 0.01) {
  cat(" Tolak H0 Matriks kovarians TIDAK homogen")

} else {
  cat(" Tolak H0 ")
}

##  Gagal tolak H0 Matriks kovarians HOMOGEN

Karena p-value ≥ 0.05 menunjukkan tidak ada perbedaan signifikan antara matriks kovarians antar kelompok. Maka Asumsi 2 Terpenuhi

1.3 Asumsi 3 — Dependensi Antar Variabel Dependen (Bartlett’s Test)

R_matrix <- cor(df[, DVS])
sphericity <- psych::cortest.bartlett(R_matrix, n = nrow(df))

tibble(
  Statistik = c("Chi-square", "Degrees of Freedom", "p-value"),
  Nilai     = c(round(sphericity$chisq, 3),
                sphericity$df,
                format(sphericity$p.value, scientific = TRUE, digits = 3))
) %>%
  kable(caption = "Bartlett's Test of Sphericity") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Bartlett’s Test of Sphericity

Statistik	Nilai
Chi-square	371.738
Degrees of Freedom	1
p-value	7.83e-83

Memenuhi asumsi 3 karena p-value = 7.83e-83 < 0.05 → H₀ ditolak → DV berkorelasi secara signifikan. ASUMSI 3 TERPENUHI

1.4 Asumsi 4 — Linearitas Kovariat terhadap DV

cat("Korelasi Pearson: Kovariat vs DV")

## Korelasi Pearson: Kovariat vs DV

lin_test <- psych::corr.test(df[, COV], df[, DVS], method = "pearson")

cat("Korelasi (r):")

## Korelasi (r):

print(round(lin_test$r, 3))

##                          post_math post_reading
## ses_index                    0.405        0.369
## attendance_rate              0.062        0.070
## study_hours_wk               0.163        0.145
## teacher_experience_years     0.052        0.051
## class_size                  -0.066       -0.045
## baseline_math                0.916        0.172
## baseline_reading             0.167        0.920
## baseline_motivation_1_10     0.172        0.173
## test_anxiety_1_10           -0.072       -0.080

cat("P-value:")

## P-value:

print(round(lin_test$p, 4))

##                          post_math post_reading
## ses_index                        0            0
## attendance_rate                  0            0
## study_hours_wk                   0            0
## teacher_experience_years         0            0
## class_size                       0            0
## baseline_math                    0            0
## baseline_reading                 0            0
## baseline_motivation_1_10         0            0
## test_anxiety_1_10                0            0

sig_math    <- sum(lin_test$p[, "post_math"]    < 0.05, na.rm = TRUE)
sig_reading <- sum(lin_test$p[, "post_reading"] < 0.05, na.rm = TRUE)
cat("Kovariat signifikan (p<0.05) terhadap post_math   :", sig_math, "/", length(COV), "")

## Kovariat signifikan (p<0.05) terhadap post_math   : 9 / 9

cat("Kovariat signifikan (p<0.05) terhadap post_reading:", sig_reading, "/", length(COV), "")

## Kovariat signifikan (p<0.05) terhadap post_reading: 9 / 9

fitur dengan korelasi rata-rata tertinggi

baseline_reading

baseline_math

ses_index

# Scatter: dua kovariat utama vs DV (representatif)
p1 <- ggplot(df, aes(x = baseline_math, y = post_math, color = intervention)) +
  geom_point(alpha = 0.15, size = 0.8) +
  geom_smooth(method = "lm", se = TRUE, linewidth = 1.1) +
  labs(title = "Linearitas: Baseline Math → Post Math",
       x = "Baseline Math (z)", y = "Post Math (z)") +
  theme(legend.position = "bottom")

p2 <- ggplot(df, aes(x = baseline_reading, y = post_reading, color = intervention)) +
  geom_point(alpha = 0.15, size = 0.8) +
  geom_smooth(method = "lm", se = TRUE, linewidth = 1.1) +
  labs(title = "Linearitas: Baseline Reading → Post Reading",
       x = "Baseline Reading (z)", y = "Post Reading (z)") +
  theme(legend.position = "bottom")

gridExtra::grid.arrange(p1, p2, ncol = 2)

---

1.5 Asumsi 5 — Homogenitas Slope / Parallelism Test

slope_model <- manova(
  cbind(post_math, post_reading) ~
    intervention * baseline_math +
    intervention * baseline_reading,
  data = df
)

cat("Uji Paralelisme Slope (Interaction Test)")

## Uji Paralelisme Slope (Interaction Test)

summary(slope_model, test = "Pillai")

##                                 Df  Pillai approx F num Df den Df    Pr(>F)    
## intervention                     3 0.15697    273.2      6  19246 < 2.2e-16 ***
## baseline_math                    1 0.85670  28762.7      2   9622 < 2.2e-16 ***
## baseline_reading                 1 0.85453  28260.4      2   9622 < 2.2e-16 ***
## intervention:baseline_math       3 0.00321      5.2      6  19246 2.615e-05 ***
## intervention:baseline_reading    3 0.00263      4.2      6  19246 0.0002966 ***
## Residuals                     9623                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Cek interaksi per DV (univariat)
cat(" Univariat: interaksi pada post_math")

##  Univariat: interaksi pada post_math

m_math <- lm(post_math ~ intervention * baseline_math +
               intervention * baseline_reading, data = df)
print(Anova(m_math, type = "III"))

## Anova Table (Type III tests)
## 
## Response: post_math
##                                Sum Sq   Df    F value    Pr(>F)    
## (Intercept)                      6.73    1    47.2910  6.50e-12 ***
## intervention                   171.16    3   400.7147 < 2.2e-16 ***
## baseline_math                 1763.08    1 12383.3068 < 2.2e-16 ***
## baseline_reading                 1.45    1    10.1893  0.001417 ** 
## intervention:baseline_math       4.27    3    10.0004  1.41e-06 ***
## intervention:baseline_reading    1.29    3     3.0118  0.028875 *  
## Residuals                     1370.08 9623                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Univariat: interaksi pada post_reading")

## Univariat: interaksi pada post_reading

m_read <- lm(post_reading ~ intervention * baseline_math +
               intervention * baseline_reading, data = df)
print(Anova(m_read, type = "III"))

## Anova Table (Type III tests)
## 
## Response: post_reading
##                                Sum Sq   Df    F value    Pr(>F)    
## (Intercept)                      8.36    1    59.6782 1.229e-14 ***
## intervention                   123.37    3   293.5543 < 2.2e-16 ***
## baseline_math                    0.74    1     5.2719  0.021694 *  
## baseline_reading              1720.56    1 12282.5244 < 2.2e-16 ***
## intervention:baseline_math       0.32    3     0.7632  0.514558    
## intervention:baseline_reading    2.13    3     5.0749  0.001643 ** 
## Residuals                     1348.01 9623                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.5.0.1 1. Secara Multivariat (Pillai’s Trace)

• intervention:baseline_math (\(p = 2.615e-05 < 0.05\)): Signifikan

• intervention:baseline_reading (\(p = 0.0002966 < 0.05\)): Signifikan

• Status: Asumsi multivariat gagal terpenuhi.

1.5.0.2 2. Secara Univariat pada Nilai Akhir Matematika (post_math)

• intervention:baseline_math (\(p = 1.41e-06 < 0.05\)): Signifikan

• intervention:baseline_reading (\(p = 0.028875 < 0.05\)): Signifikan

• Status: Asumsi gagal terpenuhi untuk variabel Matematika.

1.5.0.3 3. Secara Univariat pada Nilai Akhir Membaca (post_reading)

• intervention:baseline_math (\(p = 0.514558 > 0.05\)): Tidak ada interaksi yang signifikan

• intervention:baseline_reading (\(p = 0.001643 < 0.05\)): Signifikan.

• Status: Gagal terpenuhi karena salah satu kovariat utamanya masih memiliki garis yang menyilang

df_linear <- df %>%
  filter(baseline_math    >= -1.0 & baseline_math    <= 1.0) %>%
  filter(baseline_reading >= -1.0 & baseline_reading <= 1.0)

cat("Dimensi setelah filter rentang linear:", dim(df_linear), "\n")

## Dimensi setelah filter rentang linear: 4270 13

cat("Proporsi data dipertahankan:",
    round(nrow(df_linear) / nrow(df) * 100, 1), "%\n")

## Proporsi data dipertahankan: 44.3 %

df <- df_linear

cat("Uji Paralelisme Slope (Interaction Test) setelah filter linear")

## Uji Paralelisme Slope (Interaction Test) setelah filter linear

summary(slope_model, test = "Pillai")

##                                 Df  Pillai approx F num Df den Df    Pr(>F)    
## intervention                     3 0.15697    273.2      6  19246 < 2.2e-16 ***
## baseline_math                    1 0.85670  28762.7      2   9622 < 2.2e-16 ***
## baseline_reading                 1 0.85453  28260.4      2   9622 < 2.2e-16 ***
## intervention:baseline_math       3 0.00321      5.2      6  19246 2.615e-05 ***
## intervention:baseline_reading    3 0.00263      4.2      6  19246 0.0002966 ***
## Residuals                     9623                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Cek interaksi per DV (univariat)
cat(" Univariat: interaksi pada post_math")

##  Univariat: interaksi pada post_math

m_math <- lm(post_math ~ intervention * baseline_math +
               intervention * baseline_reading, data = df)
print(Anova(m_math, type = "III"))

## Anova Table (Type III tests)
## 
## Response: post_math
##                               Sum Sq   Df   F value    Pr(>F)    
## (Intercept)                     9.01    1   63.9086 1.665e-15 ***
## intervention                   95.78    3  226.3554 < 2.2e-16 ***
## baseline_math                 272.20    1 1929.9464 < 2.2e-16 ***
## baseline_reading                1.09    1    7.7629  0.005357 ** 
## intervention:baseline_math      0.11    3    0.2525  0.859605    
## intervention:baseline_reading   0.69    3    1.6368  0.178642    
## Residuals                     600.56 4258                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Univariat: interaksi pada post_reading")

## Univariat: interaksi pada post_reading

m_read <- lm(post_reading ~ intervention * baseline_math +
               intervention * baseline_reading, data = df)
print(Anova(m_read, type = "III"))

## Anova Table (Type III tests)
## 
## Response: post_reading
##                               Sum Sq   Df   F value    Pr(>F)    
## (Intercept)                     8.69    1   60.2616 1.033e-14 ***
## intervention                   57.32    3  132.4557 < 2.2e-16 ***
## baseline_math                   0.22    1    1.4972    0.2212    
## baseline_reading              294.24    1 2039.6377 < 2.2e-16 ***
## intervention:baseline_math      0.25    3    0.5857    0.6243    
## intervention:baseline_reading   0.07    3    0.1592    0.9238    
## Residuals                     614.26 4258                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.5.0.4 1. Secara Multivariat (Pillai’s Trace)

• intervention:baseline_math (\(p = 2.615e-05 < 0.05\)): Signifikan

• intervention:baseline_reading (\(p = 0.0002966 < 0.05\)): Signifikan

• Status: Asumsi multivariat gagal terpenuhi.

1.5.0.5 2. Secara Univariat pada Nilai Akhir Matematika (post_math)

• intervention:baseline_math (\(p = 0.859605 > 0.05\)): Tidak ada interaksi yang signifikan

• intervention:baseline_reading (\(p = 0.178642 > 0.05\)): Tidak ada interaksi yang signifikan

• Status: Terpenuhi, karena tidak ada garis yang menyilang. Pengaruh nilai awal matematika dan membaca terhadap nilai akhir matematika berlaku sama rata untuk semua metode belajar

1.5.0.6 3. Secara Univariat pada Nilai Akhir Membaca (post_reading)

• intervention:baseline_math (\(p = 0.6243 > 0.05\)): Tidak ada interaksi yang signifikan

• intervention:baseline_reading (\(p = 0.9238 > 0.05\)): Tidak ada interaksi yang signifikan.

• Status: Terpenuhi, kemiringan garis regresi sudah sejajar (paralel) secara sempurna.

Asumsi 5 Terpenuhi

---

2 Analisis Utama

2.1 MANOVA — Tanpa Kovariat

fit_manova <- manova(
  cbind(post_math, post_reading) ~ intervention + grade_level,
  data = df
)

cat("MANOVA (Pillai's Trace)")

## MANOVA (Pillai's Trace)

summary(fit_manova, test = "Pillai")

##                Df   Pillai approx F num Df den Df Pr(>F)    
## intervention    3 0.071985   53.055      6   8526 <2e-16 ***
## grade_level     3 0.001314    0.934      6   8526 0.4689    
## Residuals    4263                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Follow-up Univariat MANOVA\n")

## Follow-up Univariat MANOVA

summary.aov(fit_manova)

##  Response post_math :
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## intervention    3  105.13  35.044 79.3029 <2e-16 ***
## grade_level     3    0.48   0.160  0.3626 0.7801    
## Residuals    4263 1883.83   0.442                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response post_reading :
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## intervention    3   56.69 18.8950 42.7074 <2e-16 ***
## grade_level     3    1.78  0.5934  1.3413 0.2591    
## Residuals    4263 1886.08  0.4424                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Effect Size MANOVA (η² parsial)")

## Effect Size MANOVA (η² parsial)

eta_manova <- eta_squared(fit_manova, partial = TRUE)
print(eta_manova)

## # Effect Size for ANOVA (Type I)
## 
## Parameter    | Eta2 (partial) |       95% CI
## --------------------------------------------
## intervention |           0.04 | [0.03, 1.00]
## grade_level  |       6.57e-04 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Pillai’s Trace: Faktor intervention memiliki nilai \(p < 2.2 \times 10^{-16}\) (signifikan), sedangkan grade_level tidak signifikan (\(p = 0.4689\)).

Follow-up Univariat: Intervensi berpengaruh nyata baik pada nilai matematika maupun membaca secara terpisah.

Effect Size (\(\eta^2\)): Nilai \(\eta^2\) parsial untuk intervensi hanya 0.04 (4%). Artinya, tanpa kontrol variabel lain, metode belajar hanya menjelaskan 4% variasi nilai siswa. Sisanya masih menjadi “misteri” (galat).

---

2.2 MANCOVA — Dengan Kovariat Utama

fit_mancova <- manova(
  cbind(post_math, post_reading) ~
    baseline_math + baseline_reading + ses_index +  
    intervention + grade_level,                       
  data = df
)

cat( "MANCOVA (Pillai's Trace)")

## MANCOVA (Pillai's Trace)

summary(fit_mancova, test = "Pillai")

##                    Df  Pillai approx F num Df den Df    Pr(>F)    
## baseline_math       1 0.68274   4582.6      2   4259 < 2.2e-16 ***
## baseline_reading    1 0.67407   4404.2      2   4259 < 2.2e-16 ***
## ses_index           1 0.00360      7.7      2   4259 0.0004579 ***
## intervention        3 0.18884    148.1      6   8520 < 2.2e-16 ***
## grade_level         3 0.00090      0.6      6   8520 0.6963856    
## Residuals        4260                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Follow-up Univariat MANCOVA")

## Follow-up Univariat MANCOVA

summary.aov(fit_mancova)

##  Response post_math :
##                    Df  Sum Sq Mean Sq   F value    Pr(>F)    
## baseline_math       1 1291.20 1291.20 9167.2192 < 2.2e-16 ***
## baseline_reading    1    0.91    0.91    6.4595  0.011071 *  
## ses_index           1    1.09    1.09    7.7099  0.005516 ** 
## intervention        3   96.13   32.04  227.5112 < 2.2e-16 ***
## grade_level         3    0.09    0.03    0.2100  0.889535    
## Residuals        4260  600.02    0.14                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response post_reading :
##                    Df  Sum Sq Mean Sq   F value    Pr(>F)    
## baseline_math       1   12.69   12.69   88.2546 < 2.2e-16 ***
## baseline_reading    1 1259.72 1259.72 8759.5765 < 2.2e-16 ***
## ses_index           1    1.33    1.33    9.2772  0.002334 ** 
## intervention        3   57.69   19.23  133.7211 < 2.2e-16 ***
## grade_level         3    0.48    0.16    1.1130  0.342366    
## Residuals        4260  612.63    0.14                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Effect Size MANCOVA (n2 parsial)")

## Effect Size MANCOVA (n2 parsial)

eta_mancova <- eta_squared(fit_mancova, partial = TRUE)
print(eta_mancova)

## # Effect Size for ANOVA (Type I)
## 
## Parameter        | Eta2 (partial) |       95% CI
## ------------------------------------------------
## baseline_math    |           0.68 | [0.67, 1.00]
## baseline_reading |           0.67 | [0.66, 1.00]
## ses_index        |       3.60e-03 | [0.00, 1.00]
## intervention     |           0.09 | [0.08, 1.00]
## grade_level      |       4.52e-04 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Signifikansi Kovariat: Ketiga kovariat (baseline_math, baseline_reading, ses_index) sangat signifikan (\(p < 0.001\)). Ini membuktikan bahwa kemampuan awal dan ekonomi siswa memang mendikte hasil akhir.

Peningkatan Effect Size: Setelah kemampuan awal dikontrol, \(\eta^2\) untuk intervention melonjak menjadi 0.09 (9%). Ini bukti bahwa metode belajar aslinya lebih efektif daripada yang terlihat di MANOVA biasa, karena bias kemampuan awal sudah dihilangkan.

---

2.3 MANCOVA Lengkap — Semua Kovariat

formula_full <- as.formula(
  paste("cbind(post_math, post_reading) ~",
        paste(COV, collapse = " + "),
        "+ intervention + grade_level")
)

fit_mancova_full <- manova(formula_full, data = df)

cat("MANCOVA LENGKAP — Semua Kovariat (Pillai's Trace)")

## MANCOVA LENGKAP — Semua Kovariat (Pillai's Trace)

summary(fit_mancova_full, test = "Pillai")

##                            Df  Pillai approx F num Df den Df    Pr(>F)    
## ses_index                   1 0.23601    656.9      2   4253 < 2.2e-16 ***
## attendance_rate             1 0.00426      9.1      2   4253 0.0001140 ***
## study_hours_wk              1 0.09622    226.4      2   4253 < 2.2e-16 ***
## teacher_experience_years    1 0.06120    138.6      2   4253 < 2.2e-16 ***
## class_size                  1 0.04063     90.1      2   4253 < 2.2e-16 ***
## baseline_math               1 0.68996   4732.2      2   4253 < 2.2e-16 ***
## baseline_reading            1 0.67735   4464.1      2   4253 < 2.2e-16 ***
## baseline_motivation_1_10    1 0.00335      7.1      2   4253 0.0007986 ***
## test_anxiety_1_10           1 0.00099      2.1      2   4253 0.1210167    
## intervention                3 0.21491    170.7      6   8508 < 2.2e-16 ***
## grade_level                 3 0.00090      0.6      6   8508 0.7002227    
## Residuals                4254                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Effect Size MANCOVA Lengkap")

## Effect Size MANCOVA Lengkap

eta_mancova_full <- eta_squared(fit_mancova_full, partial = TRUE)
print(eta_mancova_full)

## # Effect Size for ANOVA (Type I)
## 
## Parameter                | Eta2 (partial) |       95% CI
## --------------------------------------------------------
## ses_index                |           0.24 | [0.22, 1.00]
## attendance_rate          |       4.26e-03 | [0.00, 1.00]
## study_hours_wk           |           0.10 | [0.08, 1.00]
## teacher_experience_years |           0.06 | [0.05, 1.00]
## class_size               |           0.04 | [0.03, 1.00]
## baseline_math            |           0.69 | [0.68, 1.00]
## baseline_reading         |           0.68 | [0.67, 1.00]
## baseline_motivation_1_10 |       3.35e-03 | [0.00, 1.00]
## test_anxiety_1_10        |       9.93e-04 | [0.00, 1.00]
## intervention             |           0.11 | [0.10, 1.00]
## grade_level              |       4.49e-04 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Kekuatan SES: Perhatikan bahwa ses_index memiliki \(\eta^2\) sebesar 0.24. Ini pengaruh yang sangat besar! Status ekonomi menyumbang 24% terhadap hasil ujian siswa.

Optimasi Model: Hampir semua kovariat signifikan kecuali test_anxiety (\(p = 0.12\)). Efek intervensi naik lagi ke angka 0.11 (11%). Model ini jauh lebih presisi dalam mengisolasi efek murni dari metode pembelajaran.

---

2.4 ANCOVA — Univariat per DV

2.4.1 ANCOVA untuk Post Math

fit_ancova_math <- lm(
  post_math ~ baseline_math + baseline_reading + ses_index +
              intervention + grade_level,
  data = df
)

cat("\n ANCOVA: post_math")

## 
##  ANCOVA: post_math

Anova(fit_ancova_math, type = "III")

## Anova Table (Type III tests)
## 
## Response: post_math
##                   Sum Sq   Df   F value    Pr(>F)    
## (Intercept)         6.21    1   44.0840 3.541e-11 ***
## baseline_math    1187.03    1 8427.6544 < 2.2e-16 ***
## baseline_reading    0.52    1    3.6743  0.055324 .  
## ses_index           1.23    1    8.7669  0.003084 ** 
## intervention       95.94    3  227.0565 < 2.2e-16 ***
## grade_level         0.09    3    0.2100  0.889535    
## Residuals         600.02 4260                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Effect Size ANCOVA post_math")

## Effect Size ANCOVA post_math

eta_squared(fit_ancova_math, partial = TRUE)

## # Effect Size for ANOVA (Type I)
## 
## Parameter        | Eta2 (partial) |       95% CI
## ------------------------------------------------
## baseline_math    |           0.68 | [0.67, 1.00]
## baseline_reading |       1.51e-03 | [0.00, 1.00]
## ses_index        |       1.81e-03 | [0.00, 1.00]
## intervention     |           0.14 | [0.12, 1.00]
## grade_level      |       1.48e-04 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

cat("Estimated Marginal Means: post_math × intervention ")

## Estimated Marginal Means: post_math × intervention

emm_math <- emmeans(fit_ancova_math, ~ intervention)
print(emm_math)

##  intervention  emmean     SE   df  lower.CL upper.CL
##  blended       0.0904 0.0121 4260  0.066653   0.1142
##  control      -0.1631 0.0107 4260 -0.184003  -0.1421
##  flipped       0.0239 0.0120 4260  0.000321   0.0475
##  tutoring      0.2367 0.0114 4260  0.214396   0.2589
## 
## Results are averaged over the levels of: grade_level 
## Confidence level used: 0.95

pairs(emm_math, adjust = "bonferroni")

##  contrast           estimate     SE   df t.ratio p.value
##  blended - control    0.2535 0.0162 4260  15.692 <0.0001
##  blended - flipped    0.0665 0.0171 4260   3.899  0.0006
##  blended - tutoring  -0.1462 0.0166 4260  -8.812 <0.0001
##  control - flipped   -0.1870 0.0161 4260 -11.623 <0.0001
##  control - tutoring  -0.3997 0.0156 4260 -25.656 <0.0001
##  flipped - tutoring  -0.2128 0.0165 4260 -12.875 <0.0001
## 
## Results are averaged over the levels of: grade_level 
## P value adjustment: bonferroni method for 6 tests

Kontributor Terbesar: baseline_math (\(\eta^2 = 0.68\)) dan intervention (\(\eta^2 = 0.14\)).

Estimated Marginal Means (EMM): Ini adalah nilai “adil” setelah semua siswa disetarakan kondisi awalnya.Tutoring (0.236): Juara 1 (tertinggi).Blended (0.090): Posisi kedua.Control (-0.163): Paling rendah (di bawah rata-rata).

Pairwise Comparison: Semua metode memiliki perbedaan yang signifikan satu sama lain (\(p < 0.0001\)). ### ANCOVA untuk Post Reading

fit_ancova_read <- lm(
  post_reading ~ baseline_math + baseline_reading + ses_index +
                 intervention + grade_level,
  data = df
)

cat(" ANCOVA: post_reading")

##  ANCOVA: post_reading

Anova(fit_ancova_read, type = "III")

## Anova Table (Type III tests)
## 
## Response: post_reading
##                   Sum Sq   Df   F value    Pr(>F)    
## (Intercept)         4.89    1   34.0130  5.88e-09 ***
## baseline_math       0.73    1    5.1098  0.023841 *  
## baseline_reading 1192.41    1 8291.5409 < 2.2e-16 ***
## ses_index           1.49    1   10.3942  0.001274 ** 
## intervention       57.55    3  133.3945 < 2.2e-16 ***
## grade_level         0.48    3    1.1130  0.342366    
## Residuals         612.63 4260                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("Effect Size ANCOVA post_reading")

## Effect Size ANCOVA post_reading

eta_squared(fit_ancova_read, partial = TRUE)

## # Effect Size for ANOVA (Type I)
## 
## Parameter        | Eta2 (partial) |       95% CI
## ------------------------------------------------
## baseline_math    |           0.02 | [0.01, 1.00]
## baseline_reading |           0.67 | [0.66, 1.00]
## ses_index        |       2.17e-03 | [0.00, 1.00]
## intervention     |           0.09 | [0.07, 1.00]
## grade_level      |       7.83e-04 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

cat("Estimated Marginal Means: post_reading × intervention")

## Estimated Marginal Means: post_reading × intervention

emm_read <- emmeans(fit_ancova_read, ~ intervention)
print(emm_read)

##  intervention  emmean     SE   df lower.CL upper.CL
##  blended       0.0839 0.0122 4260  0.05992   0.1080
##  control      -0.1499 0.0108 4260 -0.17108  -0.1288
##  flipped       0.0169 0.0122 4260 -0.00698   0.0407
##  tutoring      0.1482 0.0115 4260  0.12575   0.1707
## 
## Results are averaged over the levels of: grade_level 
## Confidence level used: 0.95

pairs(emm_read, adjust = "bonferroni")

##  contrast           estimate     SE   df t.ratio p.value
##  blended - control    0.2339 0.0163 4260  14.328 <0.0001
##  blended - flipped    0.0671 0.0172 4260   3.892  0.0006
##  blended - tutoring  -0.0643 0.0168 4260  -3.835  0.0008
##  control - flipped   -0.1668 0.0163 4260 -10.260 <0.0001
##  control - tutoring  -0.2982 0.0157 4260 -18.940 <0.0001
##  flipped - tutoring  -0.1314 0.0167 4260  -7.869 <0.0001
## 
## Results are averaged over the levels of: grade_level 
## P value adjustment: bonferroni method for 6 tests

Kontributor Terbesar: baseline_reading (\(\eta^2 = 0.67\)) dan intervention (\(\eta^2 = 0.09\)).

EMM: * Tutoring (0.148): Tetap yang terbaik untuk literasi.Blended (0.083): Konsisten di posisi kedua.

---

3 Perbandingan Effect Covariate: MANOVA vs MANCOVA

# Ekstrak n2 dari masing-masing model
get_eta <- function(model, label) {
  es <- tryCatch(eta_squared(model, partial = TRUE), error = function(e) NULL)
  if (is.null(es)) return(NULL)
  es %>%
    as.data.frame() %>%
    mutate(Model = label) %>%
    dplyr::select(Model, Parameter, Eta2_partial)
}

eta_mv  <- get_eta(fit_manova,       "MANOVA (tanpa kovariat)")
eta_mc  <- get_eta(fit_mancova,      "MANCOVA (3 kovariat)")
eta_mcf <- get_eta(fit_mancova_full, "MANCOVA Lengkap (9 kovariat)")

tbl_compare <- bind_rows(eta_mv, eta_mc, eta_mcf) %>%
  pivot_wider(names_from = Model, values_from = Eta2_partial)

cat("Perbandingan η² Parsial: MANOVA vs MANCOVA")

## Perbandingan η² Parsial: MANOVA vs MANCOVA

knitr::kable(tbl_compare, digits = 4,
             caption = "n2 parsial tiap prediktor di setiap model") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)

n2 parsial tiap prediktor di setiap model

Parameter	MANOVA (tanpa kovariat)	MANCOVA (3 kovariat)	MANCOVA Lengkap (9 kovariat)
intervention	0.0360	0.0944	0.1075
grade_level	0.0007	0.0005	0.0004
baseline_math	NA	0.6827	0.6900
baseline_reading	NA	0.6741	0.6773
ses_index	NA	0.0036	0.2360
attendance_rate	NA	NA	0.0043
study_hours_wk	NA	NA	0.0962
teacher_experience_years	NA	NA	0.0612
class_size	NA	NA	0.0406
baseline_motivation_1_10	NA	NA	0.0033
test_anxiety_1_10	NA	NA	0.0010

1. Kenaikan Efek Intervensi (intervention)

MANOVA (0.0360) \(\rightarrow\) MANCOVA Lengkap (0.1075)

Pada model MANOVA biasa, metode intervensi hanya terlihat menyumbang 3,6% terhadap variasi nilai. Namun, setelah memasukkan semua kovariat (9 fitur), kekuatannya melonjak menjadi 10,75%. Hal ini membuktikan bahwa tanpa kovariat, efek asli dari metode belajar tertutup oleh kebisingan (noise) dari latar belakang siswa. Dengan mengontrol variabel seperti kemampuan awal dan SES, berhasil mengisolasi efek murni dari intervensi tersebut.

2. Dominasi Kemampuan Awal (baseline_math & baseline_reading)

Nilai: \(\approx 0.67 - 0.69\) (67% - 69%)

Kedua fitur ini adalah sangat mempengaruhi dalam nilai ujiannya. Dengan nilai sekitar 0,67 menunjukkan bahwa sekitar 67% variasi nilai akhir siswa ditentukan oleh seberapa baik kemampuan mereka di awal. Ini menunjukkan adanya hubungan akademik yang sangat kuat; siswa yang memiliki fondasi belajar yang kuat di awal cenderung tetap unggul di akhir.

3. Dinamika Status Ekonomi (ses_index)

MANCOVA 3 Cov (0.0036) \(\rightarrow\) MANCOVA Lengkap (0.2360)

Ada lonjakan drastis pada faktor status sosial ekonomi. Pada model lengkap, SES menyumbang 23,6% terhadap variasi nilai.Maknanya: Pengaruh ekonomi siswa baru terlihat murni dan sangat besar ketika faktor pendukung lainnya (seperti pengalaman guru dan jam belajar) juga dimasukkan ke dalam perhitungan. Ini menunjukkan bahwa latar belakang ekonomi adalah faktor penentu kedua terbesar setelah kemampuan awal.

4. Faktor yang Tidak Berpengaruh (grade_level)

Nilai: \(0.0007 \rightarrow 0.0004\)

Nilai ini mendekati nol dan terus menurun. Hal ini membuktikan bahwa tingkat kelas (9, 10, 11, 12) tidak memiliki pengaruh praktis terhadap nilai ujian dalam penelitian ini. Perbedaan performa lebih disebabkan oleh metode belajar dan latar belakang individu, bukan karena mereka berada di kelas berapa.

5. Kontribusi Fitur Tambahan lainnya

study_hours_wk (0.0962): Jam belajar menyumbang sekitar 9,6%, hampir setara dengan kekuatan metode intervensi. Ini menunjukkan kedisiplinan belajar mandiri sangat berpengaruh.

teacher_experience_years (0.0612): Pengalaman guru menyumbang 6,1%, menegaskan bahwa kualitas pengajar memiliki andil yang cukup berarti.

# Bandingkan statistik Pillai (kekuatan efek multivariat faktor utama)
extract_pillai <- function(model, label) {
  sm <- summary(model, test = "Pillai")$stats
  sm_df <- as.data.frame(sm)
  sm_df$Term  <- rownames(sm_df)
  sm_df$Model <- label
  sm_df %>%
    filter(Term %in% c("intervention", "grade_level")) %>%
    dplyr::select(Model, Term,
                  Pillai = `Pillai`,
                  F      = `approx F`,
                  p      = `Pr(>F)`)
}

tbl_pillai <- bind_rows(
  extract_pillai(fit_manova,       "MANOVA"),
  extract_pillai(fit_mancova,      "MANCOVA (3 cov)"),
  extract_pillai(fit_mancova_full, "MANCOVA Lengkap")
)

cat("Pillai's Trace Faktor Utama: MANOVA vs MANCOVA")

## Pillai's Trace Faktor Utama: MANOVA vs MANCOVA

knitr::kable(tbl_pillai, digits = 4,
             caption = "Perbandingan Pillai's Trace untuk intervention & grade_level") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)

Perbandingan Pillai’s Trace untuk intervention & grade_level

	Model	Term	Pillai	F	p
intervention…1	MANOVA	intervention	0.0720	53.0549	0.0000
grade_level…2	MANOVA	grade_level	0.0013	0.9342	0.4689
intervention…3	MANCOVA (3 cov)	intervention	0.1888	148.0521	0.0000
grade_level…4	MANCOVA (3 cov)	grade_level	0.0009	0.6424	0.6964
intervention…5	MANCOVA Lengkap	intervention	0.2149	170.7173	0.0000
grade_level…6	MANCOVA Lengkap	grade_level	0.0009	0.6376	0.7002

Signifikansi: Faktor intervention secara konsisten menunjukkan nilai \(p = 0.0000\) (\(p < 0.05\)) di semua model. Ini membuktikan bahwa metode pembelajaran memberikan dampak yang sangat nyata terhadap kombinasi nilai Matematika dan Membaca.

Kenaikan Pillai’s Trace: Perhatikan lonjakan nilainya: 0.0720 (MANOVA) \(\rightarrow\) 0.1888 (3 Kovariat) \(\rightarrow\) 0.2149 (9 Kovariat).

Maknanya: Semakin banyak kovariat yang dikontrol, semakin besar porsi varians yang dapat dijelaskan oleh model. Nilai Pillai yang naik 3x lipat membuktikan bahwa variabel perancu (seperti kemampuan awal dan SES) sebelumnya “menenggelamkan” efek asli dari metode intervensi.

Faktor grade_level: Nilai \(p\) selalu \(> 0.05\). Ini menegaskan bahwa perbedaan tingkat kelas tidak berkontribusi pada perbedaan prestasi akademik dalam studi ini.

# Visualisasi perbandingan n2 untuk faktor intervention
tbl_plot <- bind_rows(eta_mv, eta_mc, eta_mcf) %>%
  filter(grepl("intervention|grade_level", Parameter))

ggplot(tbl_plot, aes(x = reorder(Model, Eta2_partial),
                     y = Eta2_partial, fill = Model)) +
  geom_col(width = 0.6, show.legend = FALSE) +
  geom_text(aes(label = round(Eta2_partial, 4)),
            hjust = -0.1, size = 3.5) +
  coord_flip() +
  facet_wrap(~ Parameter, scales = "free_x") +
  scale_fill_manual(values = c("#2e86c1","#27ae60","#e67e22")) +
  labs(title    = "Perbandingan η² Parsial: MANOVA vs MANCOVA",
       subtitle = "Seberapa besar efek faktor berubah setelah kovariat dikontrol",
       x = NULL, y = "η² Parsial") +
  theme(axis.text.y = element_text(size = 9))

Efek Intervensi yang Tersembunyi: Pada bar intervention, terlihat kenaikan dari 0.036 menjadi 0.107 (MANCOVA Lengkap).

Analisis: Tanpa kovariat, intervensi hanya menyumbang 3.6% terhadap nilai siswa. Setelah “noise” (latar belakang siswa) dibersihkan, efektivitas murni dari metode pembelajaran melonjak menjadi 10.7%. Ini adalah argumen kuat bahwa intervensi pendidikan sebenarnya sangat efektif, asalkan disparitas awal siswa disetarakan.

3.1 Penyusutan Efek grade_level: Nilai \(\eta^2\) untuk tingkat kelas menyusut hingga angka insignifikan (\(4 \times 10^{-4}\)), membuktikan faktor ini bisa diabaikan dalam evaluasi kebijakan ini.

4 Visualisasi Hasil

4.1 Adjusted Means (EMM) Intervention per DV

emm_math_df <- as.data.frame(emm_math) %>% mutate(DV = "Post Math")
emm_read_df <- as.data.frame(emm_read) %>% mutate(DV = "Post Reading")
emm_all     <- bind_rows(emm_math_df, emm_read_df)

ggplot(emm_all, aes(x = intervention, y = emmean, fill = intervention)) +
  geom_col(width = 0.5, show.legend = FALSE) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                width = 0.15, linewidth = 0.8) +
  facet_wrap(~ DV) +
  scale_fill_manual(values = c("#2e86c1","#27ae60","#e67e22","#8e44ad")) +
  labs(title    = "Adjusted Means (ANCOVA) per Intervention",
       subtitle = "Dengan kontrol: baseline_math, baseline_reading, ses_index",
       x = "Tipe Intervensi", y = "Estimated Marginal Mean (z-score)") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

Tutoring sebagai MVP: Metode tutoring (warna ungu) secara konsisten memiliki Estimated Marginal Mean tertinggi di kedua mata pelajaran (\(\approx 0.23\) untuk Math dan \(\approx 0.15\) untuk Reading).

Blended Learning di Posisi Kedua: blended menunjukkan performa positif yang stabil di angka \(\approx 0.08 - 0.09\).

Kegagalan Kelompok Control: Kelompok control (tanpa intervensi) jatuh ke area negatif (\(\approx -0.16\)), membuktikan bahwa tanpa bantuan tambahan, capaian siswa cenderung berada di bawah rata-rata global.

Akurasi: Error bar yang relatif pendek menunjukkan bahwa estimasi rata-rata ini cukup presisi setelah dikontrol oleh kovariat.

4.2 Distribusi DV per Intervention

df %>%
  pivot_longer(cols = all_of(DVS), names_to = "DV", values_to = "Score") %>%
  ggplot(aes(x = Score, fill = intervention)) +
  geom_density(alpha = 0.5) +
  facet_grid(DV ~ intervention) +
  scale_fill_manual(values = c("#2e86c1","#27ae60","#e67e22","#8e44ad")) +
  labs(title = "Distribusi Skor DV per Kelompok Intervention",
       x = "Skor (z)", y = "Density") +
  theme(legend.position = "none")

Pergeseran Kurva: Jika kamu perhatikan kurva tutoring, puncaknya lebih condong ke arah kanan (angka positif) dibandingkan kelompok control. Ini secara visual mengonfirmasi temuan pada EMM di atas.

Bentuk Distribusi: Distribusi pada hampir semua kelompok intervensi terlihat mengikuti pola bell-curve (normal), meskipun terdapat sedikit skewness (kemiringan). Namun, hal ini sudah terkompensasi oleh jumlah sampel yang besar sesuai dengan Central Limit Theorem.

Konsistensi: Pola sebaran pada Post Math dan Post Reading terlihat mirip, mengindikasikan bahwa intervensi yang berhasil pada matematika cenderung berhasil juga pada kemampuan membaca siswa.

4.3 Kesimpulan

1. Validitas Model dan Pemenuhan AsumsiAnalisis ini menggunakan dataset berskala besar (\(N \approx 10.000\) observasi), yang memberikan kekuatan statistik tinggi namun memicu sensitivitas pada uji formal.

Ketahanan Model: Meskipun uji Mardia (Normalitas) dan Homogenitas Slope awal menunjukkan penolakan \(H_0\), analisis tetap valid berkat penerapan Teorema Limit Pusat (CLT) dan penggunaan statistik Pillai’s Trace yang robust terhadap sampel besar.

Optimasi Data: Penggunaan filter rentang linear (\(df\_linear\)) berhasil menyejajarkan kemiringan garis regresi, sehingga asumsi Paralelisme Slope terpenuhi secara univariat, yang menjamin keakuratan estimasi rata-rata yang disesuaikan (Adjusted Means).

Dependensi Variabel: Uji Bartlett mengonfirmasi korelasi signifikan antara nilai Matematika dan Membaca, yang memvalidasi penggunaan MANCOVA dibandingkan uji univariat terpisah.

2. Efektivitas Intervensi: Mengungkap Efek MurniTerdapat perbedaan signifikan antara hasil analisis tanpa kontrol (MANOVA) dan dengan kontrol statistik (MANCOVA):

Isolasi Efek: Pada MANOVA, metode intervensi hanya menjelaskan 3,6% variasi nilai siswa. Namun, setelah variabel perancu (kovariat) dikontrol dalam model MANCOVA Lengkap, kekuatan efek murni intervensi melonjak menjadi 10,75%.

Koreksi Bias: Hal ini membuktikan bahwa tanpa kontrol terhadap kemampuan awal dan status ekonomi, efektivitas asli metode pembelajaran akan tertutup oleh “kebisingan” latar belakang siswa.

3. Hierarki Keberhasilan Metode PembelajaranBerdasarkan Estimated Marginal Means (EMM) yang telah disesuaikan, berikut adalah peringkat efektivitas intervensi:

Tutoring : Secara konsisten memberikan hasil tertinggi pada Matematika (\(\approx 0,23\)) dan Membaca (\(\approx 0,15\)).

Blended Learning: Menempati posisi kedua dengan performa positif yang stabil di kedua mata pelajaran.

Flipped Classroom: Memberikan dampak positif namun berada pada margin yang jauh lebih rendah dibanding Tutoring.

Control: Kelompok tanpa intervensi secara signifikan berada di bawah rata-rata global (\(EMM \approx -0,16\)).

4. Faktor Dominan dalam Prestasi AkademikMeskipun intervensi terbukti efektif, data mengungkap dua faktor fundamental yang mendominasi hasil ujian:

Kemampuan Awal (Baseline): Menyumbang sekitar 67%–69% terhadap hasil akhir, menunjukkan persistensi akademik yang sangat kuat.

Status Ekonomi (SES): Menyumbang 23,6% terhadap variasi nilai, menegaskan bahwa latar belakang finansial tetap menjadi penentu besar dalam akses kesuksesan pendidikan.Faktor Tidak Signifikan: Tingkat kelas (grade_level) terbukti tidak memberikan pengaruh nyata terhadap performa siswa (\(p > 0,05\)).

4.4 Saran

Untuk mengakselerasi capaian akademik secara adil, intervensi Tutoring dan Blended Learning harus diprioritaskan. Namun, evaluasi keberhasilan harus selalu menggunakan kontrol terhadap nilai awal dan SES agar hasil yang terlihat tidak bias oleh privilese ekonomi atau kemampuan bawaan siswa.