Meta-analysis of project-based learning on student achievement

# Contoh kode R: Meta-Analysis (metafor)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(metafor)
})
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: ringkasan statistik dari k studi primer ----
k <- sample(15:25, 1)
log_or <- tolower("") == "log or"                              # ukuran efek: log odds ratio atau SMD
d <- data.frame(studi = sprintf("Studi %02d", 1:k),
                n_perlakuan = sample(20:120, k, replace = TRUE),
                n_kontrol = sample(20:120, k, replace = TRUE))
d$intervention_duration <- round(pmax(rnorm(k, 50, 10), 0.1 * abs(50)), 1)   # moderator utama
d$study_sample_size <- if (grepl("year", tolower(""))) {
  sample(2008:2024, k, replace = TRUE)
} else if (grepl("sample|participant", tolower("study_sample_size"))) {
  d$n_perlakuan + d$n_kontrol
} else round(pmax(rnorm(k, 30, 6), 0.1 * abs(30)), 1)
mu <- runif(1, 0.25, 0.55)                                            # efek rata-rata sebenarnya
slope <- runif(1, 0.08, 0.2)                                          # pengaruh moderator (per SD)
tau <- runif(1, 0.08, 0.2)                                            # heterogenitas antarstudi
theta <- mu + slope * as.numeric(scale(d$intervention_duration)) + rnorm(k, 0, tau)
if (log_or) {
  # kejadian pada kelompok kontrol dan perlakuan (odds kelompok perlakuan dinaikkan exp(2 * theta))
  pk <- runif(k, 0.15, 0.4)
  pp <- plogis(qlogis(pk) + 2 * theta)
  d$kejadian_perlakuan <- rbinom(k, d$n_perlakuan, pp)
  d$kejadian_kontrol <- rbinom(k, d$n_kontrol, pk)
  es <- metafor::escalc("OR", ai = kejadian_perlakuan, n1i = n_perlakuan, ci = kejadian_kontrol, n2i = n_kontrol,
                        data = d, add = 0.5, to = "only0")
} else {
  # rerata dan simpangan baku dua kelompok pada skala hasil studi
  d$sd_perlakuan <- round(runif(k, 8, 14), 2)
  d$sd_kontrol <- round(d$sd_perlakuan * runif(k, 0.85, 1.15), 2)
  d$rerata_kontrol <- round(runif(k, 40, 70), 2)
  d$rerata_perlakuan <- round(d$rerata_kontrol + theta * sqrt((d$sd_perlakuan^2 + d$sd_kontrol^2) / 2) +
                                rnorm(k, 0, sqrt(d$sd_perlakuan^2 / d$n_perlakuan + d$sd_kontrol^2 / d$n_kontrol)), 2)
  es <- metafor::escalc("SMD", m1i = rerata_perlakuan, sd1i = sd_perlakuan, n1i = n_perlakuan,
                        m2i = rerata_kontrol, sd2i = sd_kontrol, n2i = n_kontrol, data = d)
}
d$pjbl_effect_on_achievement <- round(as.numeric(es$yi), 4)
d$varians <- round(as.numeric(es$vi), 5)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Model efek acak (REML) ----
res <- metafor::rma(yi = pjbl_effect_on_achievement, vi = varians, data = d, method = "REML", slab = studi)
print(res)
pr <- predict(res)                                                    # efek gabungan + prediction interval
print(pr)

# ---- Bias publikasi: uji Egger dan trim-and-fill ----
eg <- metafor::regtest(res)
print(eg)
tf <- metafor::trimfill(res)
print(tf)

# ---- Meta-regresi dengan moderator intervention_duration ----
mr <- metafor::rma(yi = pjbl_effect_on_achievement, vi = varians, mods = ~ intervention_duration, data = d, method = "REML")
print(mr)


# ---- Grafik ----
bobot <- weights(res)
df <- data.frame(studi = sub("Studi", "Study", d$studi), yi = d$pjbl_effect_on_achievement, se = sqrt(d$varians), w = bobot)
df$bawah <- df$yi - 1.96 * df$se; df$atas <- df$yi + 1.96 * df$se
df$studi <- factor(df$studi, levels = c("Pooled (RE)", rev(df$studi)))
est <- as.numeric(res$b)
berlian <- data.frame(x = c(res$ci.lb, est, res$ci.ub, est), y = c(1, 1.35, 1, 0.65))
p1 <- ggplot(df, aes(yi, studi)) +
  geom_vline(xintercept = 0, colour = warna$redup) +
  geom_vline(xintercept = est, colour = pal[2], linetype = "dashed") +
  geom_errorbarh(aes(xmin = bawah, xmax = atas), height = 0, colour = warna$redup, linewidth = 0.6) +
  geom_point(aes(size = w), shape = 15, colour = pal[1]) +
  annotate("segment", x = pr$pi.lb, xend = pr$pi.ub, y = 1, yend = 1, colour = pal[2], linewidth = 0.6) +
  geom_polygon(data = berlian, aes(x, y), fill = pal[2], colour = pal[2]) +
  scale_size_continuous(range = c(1.2, 4), guide = "none") +
  scale_y_discrete(drop = FALSE) +
  labs(title = "Forest plot (random effects, REML)", x = "pjbl_effect_on_achievement", y = NULL) +
  theme(axis.text.y = element_text(size = rel(0.75)))
se_max <- max(df$se) * 1.05
corong <- data.frame(x = c(est - 1.96 * se_max, est, est + 1.96 * se_max), y = c(se_max, 0, se_max))
isi <- if (tf$k0 > 0) data.frame(yi = as.numeric(tf$yi[tf$fill]), se = sqrt(as.numeric(tf$vi[tf$fill]))) else NULL
p2 <- ggplot(df, aes(yi, se)) +
  geom_polygon(data = corong, aes(x, y), fill = warna$grid, colour = warna$redup, alpha = 0.6, linetype = "dashed") +
  geom_vline(xintercept = est, colour = pal[2]) +
  geom_point(colour = pal[1], size = 2, alpha = 0.85) +
  scale_y_reverse() +
  labs(title = sprintf("Funnel plot (Egger p = %.3f)", eg$pval), x = "pjbl_effect_on_achievement", y = "Standard error")
if (!is.null(isi)) p2 <- p2 + geom_point(data = isi, shape = 21, colour = pal[3], fill = NA, size = 2.2, stroke = 1)
gx <- seq(min(d$intervention_duration), max(d$intervention_duration), length.out = 50)
pg <- predict(mr, newmods = gx)
dg <- data.frame(x = gx, m = pg$pred, lo = pg$ci.lb, hi = pg$ci.ub)
p3 <- ggplot() +
  geom_ribbon(data = dg, aes(x, ymin = lo, ymax = hi), fill = pal[2], alpha = 0.18) +
  geom_line(data = dg, aes(x, m), colour = pal[2], linewidth = 1.1) +
  geom_point(data = data.frame(x = d$intervention_duration, y = d$pjbl_effect_on_achievement, w = bobot), aes(x, y, size = w),
             shape = 21, fill = pal[1], colour = warna$panel, alpha = 0.7) +
  scale_size_continuous(range = c(1.5, 7), guide = "none") +
  labs(title = "Meta-regression bubble plot", x = "intervention_duration", y = "pjbl_effect_on_achievement")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else (p2 / p3) | p1

print(p_hasil)