Efektivitas problem-based learning, pembelajaran kooperatif, dan ceramah terhadap skor tes fisika

# Contoh kode R: ANOVA, ANCOVA and Tukey Post Hoc Test
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(car)
  library(emmeans)
  library(scales)
})
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: tiga kelompok perlakuan dengan satu kovariat ----
n <- 200
lv <- c("A", "B", "C")
d <- data.frame(teaching_method = factor(sample(rep(lv, length.out = n)), levels = lv),
                pretest_score = round(rnorm(n, 30, 6), 2))
z2 <- as.numeric(scale(d$pretest_score))
efek <- c(0, runif(1, 0.25, 0.6), runif(1, 0.55, 1.1))[sample(3)]   # beda rata-rata antarkelompok (dalam SD)
b_kov <- runif(1, 0.35, 0.7)                                         # pengaruh kovariat
d$physics_test_score <- round(100 + 15 * (efek[as.integer(d$teaching_method)] - mean(efek) + b_kov * z2 +
                                    rnorm(n, 0, runif(1, 0.55, 0.8))), 2)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis 1: ANOVA satu arah ----
fit_aov <- aov(physics_test_score ~ teaching_method, data = d)
tab_aov <- summary(fit_aov)[[1]]
print(tab_aov)
print(aggregate(physics_test_score ~ teaching_method, data = d, FUN = function(v) c(rerata = mean(v), sd = sd(v), n = length(v))))

# ---- Uji asumsi ANOVA ----
lev <- car::leveneTest(physics_test_score ~ teaching_method, data = d, center = median)   # homogenitas varians (Levene/Brown-Forsythe)
sw <- shapiro.test(residuals(fit_aov))                              # normalitas residual
print(lev); print(sw)
eta2 <- tab_aov[1, "Sum Sq"] / sum(tab_aov[, "Sum Sq"])             # eta kuadrat (= partial eta2 pada ANOVA satu arah)

# ---- Uji lanjut Tukey HSD ----
tk <- TukeyHSD(fit_aov, conf.level = 0.95)
print(tk)
tkd <- as.data.frame(tk[[1]])
names(tkd) <- c("selisih", "bawah", "atas", "p_adj")

# ---- Analisis 2: ANCOVA (kovariat pretest_score) ----
fit_slope <- lm(physics_test_score ~ pretest_score * teaching_method, data = d)                  # uji homogenitas kemiringan regresi
a_slope <- car::Anova(fit_slope, type = 2)
print(a_slope)
p_slope <- a_slope[3, "Pr(>F)"]
fit_anc <- lm(physics_test_score ~ pretest_score + teaching_method, data = d)
a_anc <- car::Anova(fit_anc, type = 2)
print(a_anc)
p_anc <- a_anc[2, "Pr(>F)"]
peta2 <- a_anc[2, "Sum Sq"] / (a_anc[2, "Sum Sq"] + a_anc["Residuals", "Sum Sq"])   # partial eta kuadrat kelompok
em <- emmeans::emmeans(fit_anc, specs = "teaching_method")                   # rata-rata terkoreksi kovariat
print(em)
print(pairs(em, adjust = "tukey"))
emd <- as.data.frame(summary(em))
names(emd)[1:6] <- c("kelompok", "rerata", "se", "db", "bawah", "atas")


# ---- Grafik ----
bungkus <- scales::label_wrap(12)
pendek <- function(s, w = 12) ifelse(nchar(s) > w, paste0(trimws(substr(s, 1, w - 1)), "."), s)
rata <- aggregate(physics_test_score ~ teaching_method, data = d, FUN = mean)
p1 <- ggplot(d, aes(teaching_method, physics_test_score)) +
  geom_boxplot(aes(fill = teaching_method), alpha = 0.25, outlier.shape = NA, colour = warna$redup, width = 0.6) +
  geom_jitter(aes(colour = teaching_method), width = 0.15, alpha = 0.5, size = 1.3) +
  geom_point(data = rata, shape = 23, size = 3.2, fill = warna$panel, colour = warna$teks, stroke = 1.1) +
  scale_fill_manual(values = pal[1:3], guide = "none") +
  scale_colour_manual(values = pal[1:3], guide = "none") +
  scale_x_discrete(labels = bungkus) +
  labs(title = "Distribution by group (diamond = mean)", x = NULL, y = "physics_test_score")
pasangan <- combn(lv, 2)
tkd$pasangan <- factor(paste(pendek(pasangan[2, ]), "-", pendek(pasangan[1, ])),
                       levels = rev(paste(pendek(pasangan[2, ]), "-", pendek(pasangan[1, ]))))
tkd$nyata <- ifelse(tkd$p_adj < 0.05, "p < 0.05", "n.s.")
p2 <- ggplot(tkd, aes(selisih, pasangan, colour = nyata)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_errorbarh(aes(xmin = bawah, xmax = atas), height = 0.25, linewidth = 0.9) +
  geom_point(size = 2.8) +
  scale_colour_manual(values = c("p < 0.05" = pal[2], "n.s." = warna$redup), name = NULL, drop = FALSE) +
  labs(title = "Tukey HSD: mean differences (95% CI)", x = "Difference in physics_test_score", y = NULL) +
  theme(legend.position = "bottom")
garis <- expand.grid(pretest_score = seq(min(d$pretest_score), max(d$pretest_score), length.out = 30), teaching_method = factor(lv, levels = lv))
garis$prediksi <- predict(fit_anc, newdata = garis)
p3 <- ggplot(d, aes(pretest_score, physics_test_score, colour = teaching_method)) +
  geom_point(alpha = 0.45, size = 1.3) +
  geom_line(data = garis, aes(y = prediksi), linewidth = 1.1) +
  scale_colour_manual(values = pal[1:3], name = NULL, labels = pendek) +
  labs(title = "ANCOVA: parallel fitted lines", x = "pretest_score", y = "physics_test_score") +
  theme(legend.position = "bottom")
mentah <- aggregate(physics_test_score ~ teaching_method, data = d, FUN = function(v) c(m = mean(v), se = sd(v) / sqrt(length(v))))
dm <- rbind(data.frame(kelompok = factor(lv, levels = lv), rerata = mentah$physics_test_score[, "m"],
                       bawah = mentah$physics_test_score[, "m"] - 1.96 * mentah$physics_test_score[, "se"],
                       atas = mentah$physics_test_score[, "m"] + 1.96 * mentah$physics_test_score[, "se"], jenis = "Raw mean"),
            data.frame(kelompok = emd$kelompok, rerata = emd$rerata, bawah = emd$bawah, atas = emd$atas,
                       jenis = "Adjusted mean"))
dm$jenis <- factor(dm$jenis, levels = c("Raw mean", "Adjusted mean"))
p4 <- ggplot(dm, aes(kelompok, rerata, colour = jenis)) +
  geom_pointrange(aes(ymin = bawah, ymax = atas), position = position_dodge(width = 0.5), size = 0.45, linewidth = 0.9) +
  scale_colour_manual(values = c("Raw mean" = warna$redup, "Adjusted mean" = pal[4]), name = NULL) +
  scale_x_discrete(labels = bungkus) +
  labs(title = "Raw vs covariate-adjusted means", x = NULL, y = "physics_test_score") +
  theme(legend.position = "bottom")
p_hasil <- if (tata == 1) (p1 | p2) / (p3 | p4) else (p1 | p3) / (p4 | p2)

print(p_hasil)