1 Penggunaan dan tujuan

Versi R dari notebook Python untuk mahasiswa S1 Kesehatan Masyarakat. Buka file ini di RStudio, lalu klik Knit → Knit to HTML. Jika belum tersedia, jalankan sekali pada Console:

install.packages(c("rmarkdown", "knitr"), repos = "https://cloud.r-project.org")

Seluruh analisis memakai fungsi dasar R sehingga tidak membutuhkan paket statistik tambahan. Data dummy dibuat langsung dalam dokumen: 120 peserta, 30 per kombinasi. Angka simulasi merupakan contoh pembelajaran, bukan hasil penelitian nyata. Seed yang sama pada R dan Python tidak menghasilkan sampel identik, karena generator dan urutan pengambilan angka acaknya berbeda. Rancangan dan mekanisme simulasi tetap sama.

Alur demo: rancangan, data, deskriptif, grafik, ANOVA, efek dan interaksi, pemeriksaan asumsi. ANCOVA merupakan pendalaman.

2 Rancangan dan mekanisme simulasi

\[D = 2 + 2A + B + 5AB + 0{,}15(SBP_0-150) + \epsilon, \quad \epsilon \sim N(0,3^2).\]

Pada baseline 150 mmHg, rata-rata populasi ditetapkan sebesar 2, 3, 4, dan 10 mmHg. Rata-rata sampel akan berbeda karena variasi acak. Jumlah 30 per sel dipilih untuk demonstrasi, bukan dari perhitungan power of test.

SEED <- 2026
N_PER_CELL <- 30
INTERACTION <- 5  # Ubah menjadi 0 untuk skenario tanpa interaksi
SD_ERROR <- 3
set.seed(SEED)

kombinasi <- expand.grid(A = 0:1, B = 0:1)
alokasi <- kombinasi[rep(seq_len(nrow(kombinasi)), each = N_PER_CELL), ]
alokasi <- alokasi[sample(seq_len(nrow(alokasi))), ]
rownames(alokasi) <- NULL
N <- nrow(alokasi)
dat <- data.frame(id = sprintf("P%03d", seq_len(N)), alokasi)
dat$usia <- sample(30:70, N, replace = TRUE)
dat$sbp_awal <- rnorm(N, 150, 10)
dat$penurunan <- with(dat, 2 + 2*A + B + INTERACTION*A*B +
                         0.15*(sbp_awal - 150) + rnorm(N, 0, SD_ERROR))
dat$sbp_minggu8 <- dat$sbp_awal - dat$penurunan
dat$edukasi <- factor(dat$A, levels = 0:1, labels = c("Standar", "Intensif"))
dat$whatsapp <- factor(dat$B, levels = 0:1, labels = c("Tanpa", "Dengan"))
dat$baseline_c <- dat$sbp_awal - 150
knitr::kable(head(dat, 10), digits = 2)
id A B usia sbp_awal penurunan sbp_minggu8 edukasi whatsapp baseline_c
P001 1 1 63 152.16 7.81 144.35 Intensif Dengan 2.16
P002 1 1 63 156.13 8.64 147.49 Intensif Dengan 6.13
P003 1 0 42 164.46 3.57 160.89 Intensif Tanpa 14.46
P004 1 0 39 156.58 9.96 146.62 Intensif Tanpa 6.58
P005 1 1 58 153.75 12.23 141.52 Intensif Dengan 3.75
P006 1 1 53 143.25 10.79 132.46 Intensif Dengan -6.75
P007 1 1 43 155.80 4.96 150.84 Intensif Dengan 5.80
P008 1 0 59 145.86 2.24 143.62 Intensif Tanpa -4.14
P009 1 0 67 171.81 10.74 161.06 Intensif Tanpa 21.81
P010 0 0 48 153.43 6.68 146.75 Standar Tanpa 3.43 
dim(dat)
## [1] 120  10

3 Audit data dan deskriptif

Satu baris mewakili satu peserta. Alokasi seimbang tidak menjamin karakteristik awal identik.

stopifnot(!anyDuplicated(dat$id), !anyNA(dat),
          isTRUE(all.equal(dat$penurunan, dat$sbp_awal - dat$sbp_minggu8)))
jumlah <- with(dat, table(edukasi, whatsapp))
stopifnot(all(jumlah == N_PER_CELL))
knitr::kable(jumlah)
Tanpa Dengan
Standar 30 30
Intensif 30 30 
grid <- expand.grid(A = 0:1, B = 0:1)
ringkasan <- do.call(rbind, lapply(seq_len(nrow(grid)), function(i) {
  d <- subset(dat, A == grid$A[i] & B == grid$B[i])
  n <- nrow(d); m <- mean(d$penurunan); sd_y <- sd(d$penurunan)
  se <- sd_y / sqrt(n); margin <- qt(.975, n - 1)*se
  data.frame(A = grid$A[i], B = grid$B[i],
             edukasi = as.character(d$edukasi[1]),
             whatsapp = as.character(d$whatsapp[1]), n = n,
             usia_mean = mean(d$usia), baseline_mean = mean(d$sbp_awal),
             akhir_mean = mean(d$sbp_minggu8), penurunan_mean = m,
             penurunan_sd = sd_y, se = se, ci_low = m-margin, ci_high = m+margin)
}))
knitr::kable(ringkasan, digits = 2)
A B edukasi whatsapp n usia_mean baseline_mean akhir_mean penurunan_mean penurunan_sd se ci_low ci_high
0 0 Standar Tanpa 30 47.97 150.20 148.01 2.19 3.21 0.59 0.99 3.39
1 0 Intensif Tanpa 30 47.07 152.21 148.02 4.19 4.26 0.78 2.60 5.79
0 1 Standar Dengan 30 46.90 148.37 145.23 3.14 2.74 0.50 2.12 4.16
1 1 Intensif Dengan 30 48.33 153.56 143.70 9.86 3.59 0.66 8.52 11.20

4 Grafik interaksi dan distribusi

CI pada grafik adalah CI rata-rata tiap sel. Tumpang tindih CI bukan uji formal interaksi. Garis tidak perlu berpotongan untuk menunjukkan pola interaksi.

op <- par(mfrow = c(1, 2), mar = c(5, 4, 3, 1))
warna <- c("#156082", "#E97132")
ylim <- range(c(ringkasan$ci_low, ringkasan$ci_high))
plot(NA, xlim = c(-.1, 1.1), ylim = ylim, xaxt = "n",
     xlab = "Bentuk edukasi", ylab = "Penurunan tekanan darah (mmHg)",
     main = "Rata-rata dan CI 95%")
axis(1, at = 0:1, labels = c("Standar", "Intensif"))
for (b in 0:1) {
  d <- ringkasan[ringkasan$B == b, ]; d <- d[order(d$A), ]
  lines(d$A, d$penurunan_mean, type = "b", pch = 19, col = warna[b+1], lwd = 2)
  arrows(d$A, d$ci_low, d$A, d$ci_high, angle = 90, code = 3,
         length = .07, col = warna[b+1])
}
legend("topleft", c("Tanpa WhatsApp", "Dengan WhatsApp"),
       col = warna, lty = 1, pch = 19, bty = "n", cex = .8)
kelompok <- factor(paste(dat$A, dat$B, sep = "_"),
                   levels = c("0_0", "0_1", "1_0", "1_1"),
                   labels = c("Standar\nTanpa", "Standar\nDengan",
                              "Intensif\nTanpa", "Intensif\nDengan"))
boxplot(dat$penurunan ~ kelompok, ylab = "Penurunan (mmHg)",
        xlab = "", main = "Variasi antarindividu", col = "#EDF4F7")

par(op)

5 Model dan kontras efek

Model dengan coding 0/1:

\[D = \beta_0 + \beta_A A + \beta_B B + \beta_{AB}AB + e.\]

model <- lm(penurunan ~ A*B, data = dat)

# Kolom mengikuti nama koefisien model, sehingga tidak bergantung pada urutan implisit.
L <- matrix(0, nrow = 7, ncol = length(coef(model)),
            dimnames = list(c("Edukasi saat tanpa WhatsApp",
                              "Edukasi saat dengan WhatsApp",
                              "WhatsApp saat edukasi standar",
                              "WhatsApp saat edukasi intensif",
                              "Interaksi: selisih efek edukasi",
                              "Edukasi marginal (bobot sama)",
                              "WhatsApp marginal (bobot sama)"), names(coef(model))))
L[1, "A"] <- 1
L[2, c("A", "A:B")] <- 1
L[3, "B"] <- 1
L[4, c("B", "A:B")] <- 1
L[5, "A:B"] <- 1
L[6, c("A", "A:B")] <- c(1, .5)
L[7, c("B", "A:B")] <- c(1, .5)

contrast_table <- function(fit, contrasts, covariance = vcov(fit)) {
  contrasts <- contrasts[, names(coef(fit)), drop = FALSE]
  effect <- drop(contrasts %*% coef(fit))
  se <- sqrt(diag(contrasts %*% covariance %*% t(contrasts)))
  t_stat <- effect/se
  df_res <- df.residual(fit)
  data.frame(Kontras = rownames(contrasts), Efek_mmHg = effect, SE = se,
             CI95_low = effect - qt(.975, df_res)*se,
             CI95_high = effect + qt(.975, df_res)*se,
             t = t_stat, p = 2*pt(-abs(t_stat), df_res), row.names = NULL)
}
effects <- contrast_table(model, L)
effects$p_Holm_4_efek_sederhana <- NA_real_
effects$p_Holm_4_efek_sederhana[1:4] <- p.adjust(effects$p[1:4], method = "holm")
knitr::kable(effects, digits = 4)
Kontras Efek_mmHg SE CI95_low CI95_high t p p_Holm_4_efek_sederhana
Edukasi saat tanpa WhatsApp 2.0016 0.9028 0.2134 3.7898 2.2170 0.0286 0.0571
Edukasi saat dengan WhatsApp 6.7241 0.9028 4.9359 8.5123 7.4476 0.0000 0.0000
WhatsApp saat edukasi standar 0.9461 0.9028 -0.8421 2.7343 1.0479 0.2969 0.2969
WhatsApp saat edukasi intensif 5.6686 0.9028 3.8804 7.4568 6.2785 0.0000 0.0000
Interaksi: selisih efek edukasi 4.7225 1.2768 2.1936 7.2513 3.6986 0.0003 NA
Edukasi marginal (bobot sama) 4.3628 0.6384 3.0984 5.6273 6.8339 0.0000 NA
WhatsApp marginal (bobot sama) 3.3073 0.6384 2.0429 4.5718 5.1806 0.0000 NA

CI adalah interval individual, belum dikoreksi untuk inferensi simultan. Koreksi Holm pada nilai p ditambahkan untuk keluarga empat efek sederhana.

6 ANOVA dua jalur dengan interaksi

Untuk mempertahankan pengujian seperti notebook Python, kita uji efek marginal dan interaksi sebagai kontras satu derajat bebas. Pada desain 2×2 seimbang ini, uji tersebut ekuivalen dengan pengujian efek A, B, dan AB pada Type III dengan sum coding.

F untuk satu kontras sama dengan kuadrat statistik t. Kita tidak memakai anova(model) sebagai pengganti Type III karena fungsi itu secara default menghasilkan jumlah kuadrat berurutan (Type I).

idx <- c(6, 7, 5)
F_values <- effects$t[idx]^2
mse <- deviance(model)/df.residual(model)
anova_table <- data.frame(
  Sumber = c("Edukasi (efek marginal)", "WhatsApp (efek marginal)",
             "Interaksi edukasi x WhatsApp", "Residual"),
  df = c(1, 1, 1, df.residual(model)),
  sum_sq = c(F_values*mse, deviance(model)),
  F = c(F_values, NA),
  p = c(pf(F_values, 1, df.residual(model), lower.tail = FALSE), NA))
knitr::kable(anova_table, digits = 5)
Sumber df sum_sq F p
Edukasi (efek marginal) 1 571.0306 46.70268 0.00000
WhatsApp (efek marginal) 1 328.1519 26.83844 0.00000
Interaksi edukasi x WhatsApp 1 167.2619 13.67979 0.00033
Residual 116 1418.3245 NA NA

H0 interaksi: perbedaan efek edukasi antar kondisi WhatsApp sama dengan nol. Nilai p bukan ukuran besar manfaat. Jika interaksi relevan, efek sederhana memberikan konteks penting untuk keputusan.

7 Contoh interpretasi otomatis

for (i in c(1, 2, 5)) {
  cat(sprintf("- **%s:** %.2f mmHg (CI 95%% %.2f sampai %.2f; p = %.4g).\n",
              effects$Kontras[i], effects$Efek_mmHg[i],
              effects$CI95_low[i], effects$CI95_high[i], effects$p[i]))
}
if (effects$p[5] < .05) {
  cat("\nPada simulasi ini, terdapat bukti statistik interaksi pada skala mmHg.\n")
} else {
  cat("\nPada simulasi ini, bukti interaksi belum kuat. Ini tidak membuktikan efek identik.\n")
}

Pada simulasi ini, terdapat bukti statistik interaksi pada skala mmHg.

cat("\nKeputusan program juga mempertimbangkan biaya, akses digital, dan kepatuhan.\n")

Keputusan program juga mempertimbangkan biaya, akses digital, dan kepatuhan.

8 Pemeriksaan asumsi dan sensitivitas

Independensi berasal dari rancangan dan proses pengumpulan data, bukan dari uji normalitas. Periksa normalitas residual, kesamaan ragam, dan kemungkinan observasi berpengaruh. p>0,05 pada uji asumsi tidak membuktikan asumsi benar.

op <- par(mfrow = c(1, 2))
plot(fitted(model), resid(model), pch = 19, col = "#156082",
     xlab = "Prediksi penurunan (mmHg)", ylab = "Residual",
     main = "Residual vs fitted")
abline(h = 0, lty = 2)
qqnorm(resid(model), pch = 19, main = "Q-Q plot residual")
qqline(resid(model), col = "#E97132", lwd = 2)

par(op)
shapiro.test(resid(model))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.9806, p-value = 0.08095
# Levene dengan pusat median (Brown–Forsythe), ekuivalen dengan scipy center='median'.
abs_dev <- abs(dat$penurunan - ave(dat$penurunan, kelompok, FUN = median))
anova(lm(abs_dev ~ kelompok))
## Analysis of Variance Table
## 
## Response: abs_dev
##            Df Sum Sq Mean Sq F value  Pr(>F)  
## kelompok    3  29.40  9.8009  2.1833 0.09374 .
## Residuals 116 520.74  4.4891                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Jika ragam tidak seragam, standard error HC3 dapat menjadi analisis sensitivitas. HC3 tidak memperbaiki dependensi klaster, kesalahan model, atau bias alokasi.

X <- model.matrix(model)
h <- hatvalues(model)
bread <- solve(crossprod(X))
w <- (resid(model)/(1-h))^2
cov_hc3 <- bread %*% crossprod(X, X*w) %*% bread
robust_effects <- contrast_table(model, L, covariance = cov_hc3)
knitr::kable(robust_effects, digits = 4)
Kontras Efek_mmHg SE CI95_low CI95_high t p
Edukasi saat tanpa WhatsApp 2.0016 0.9910 0.0387 3.9645 2.0197 0.0457
Edukasi saat dengan WhatsApp 6.7241 0.8392 5.0619 8.3863 8.0122 0.0000
WhatsApp saat edukasi standar 0.9461 0.7835 -0.6058 2.4980 1.2075 0.2297
WhatsApp saat edukasi intensif 5.6686 1.0357 3.6173 7.7198 5.4734 0.0000
Interaksi: selisih efek edukasi 4.7225 1.2986 2.1503 7.2946 3.6365 0.0004
Edukasi marginal (bobot sama) 4.3628 0.6493 3.0768 5.6489 6.7191 0.0000
WhatsApp marginal (bobot sama) 3.3073 0.6493 2.0213 4.5934 5.0935 0.0000

9 Pendalaman: ANCOVA dengan baseline

Model tambahan: sbp_minggu8 ~ baseline_c + A*B. Perbandingan dilakukan pada baseline yang sama, yaitu 150 mmHg. Model mengasumsikan hubungan linear baseline–outcome dan slope baseline sama antarkelompok.

Pada outcome akhir, efek negatif berarti tekanan darah lebih rendah. Untuk konsistensi, tabel manfaat membalik tanda kontras sehingga nilai positif berarti manfaat tambahan. Rencanakan ANCOVA sebelum melihat hasil, bukan karena nilai p lebih kecil.

ancova <- lm(sbp_minggu8 ~ baseline_c + A*B, data = dat)
newdat <- grid
newdat$baseline_c <- 0
pred <- predict(ancova, newdata = newdat, interval = "confidence", level = .95)
adjusted <- cbind(newdat, as.data.frame(pred))
names(adjusted)[names(adjusted) == "fit"] <- "SBP_akhir_adjusted"
knitr::kable(adjusted, digits = 3)
A B baseline_c SBP_akhir_adjusted lwr upr
0 0 0 147.840 146.675 149.005
1 0 0 146.157 144.982 147.331
0 1 0 146.603 145.433 147.773
1 1 0 140.702 139.513 141.891 
L_ancova <- matrix(0, nrow = nrow(L), ncol = length(coef(ancova)),
                   dimnames = list(rownames(L), names(coef(ancova))))
L_ancova[, colnames(L)] <- -L
ancova_effects <- contrast_table(ancova, L_ancova)
knitr::kable(ancova_effects, digits = 4)
Kontras Efek_mmHg SE CI95_low CI95_high t p
Edukasi saat tanpa WhatsApp 1.6832 0.8345 0.0302 3.3363 2.0170 0.0460
Edukasi saat dengan WhatsApp 5.9011 0.8503 4.2169 7.5854 6.9402 0.0000
WhatsApp saat edukasi standar 1.2368 0.8341 -0.4153 2.8889 1.4829 0.1408
WhatsApp saat edukasi intensif 5.4547 0.8330 3.8047 7.1047 6.5484 0.0000
Interaksi: selisih efek edukasi 4.2179 1.1812 1.8782 6.5576 3.5708 0.0005
Edukasi marginal (bobot sama) 3.7922 0.6007 2.6022 4.9821 6.3125 0.0000
WhatsApp marginal (bobot sama) 3.3458 0.5882 2.1807 4.5108 5.6884 0.0000

10 Untuk diskusi dan review

  1. Mengapa ada 120 unit eksperimen? Apa yang berubah jika alokasi dilakukan pada empat desa?
  2. Ubah INTERACTION <- 0 dan Knit ulang. Apakah estimasi interaksi harus tepat nol?
  3. Ubah N_PER_CELL menjadi 10 atau 60. Apa yang terjadi pada lebar CI?
  4. Mengapa koefisien A berbeda dari efek edukasi marginal?
  5. Mengapa “signifikan di satu kelompok, tidak signifikan di kelompok lain” bukan bukti interaksi?
  6. Informasi apa lagi yang diperlukan jika akses WhatsApp rendah?

Petunjuk: unit mendapat alokasi independen; alokasi desa membutuhkan replikasi desa dan analisis klaster. Estimasi sampel mengandung noise meski interaksi populasi nol. CI biasanya lebih sempit saat n lebih besar, tetapi satu simulasi tidak menilai power. Koefisien A mengondisikan B=0, efek marginal merata-ratakan atas B. Perbedaan dua keputusan signifikansi bukan uji perbedaan efek.

11 Ekspor data dan hasil

CSV disimpan di folder hasil_faktorial_dummy_R di direktori kerja saat Knit. Di RStudio, buka tab Files untuk mengambil hasil. Data disimpan dengan presisi penuh.

out <- "hasil_faktorial_dummy_R"
dir.create(out, showWarnings = FALSE)
write.csv(dat, file.path(out, "data_dummy_faktorial_2x2.csv"), row.names = FALSE)
write.csv(ringkasan, file.path(out, "deskriptif.csv"), row.names = FALSE)
write.csv(anova_table, file.path(out, "anova_dua_jalur.csv"), row.names = FALSE)
write.csv(effects, file.path(out, "efek_dan_interaksi.csv"), row.names = FALSE)
write.csv(ancova_effects, file.path(out, "ancova_efek.csv"), row.names = FALSE)
writeLines(c("DATA SIMULASI untuk pembelajaran, bukan hasil penelitian nyata.",
             sprintf("Seed=%s; n per sel=%s; interaksi populasi=%s mmHg.",
                     SEED, N_PER_CELL, INTERACTION),
             capture.output(print(effects))), file.path(out, "ringkasan.txt"))
cat("Folder hasil:", normalizePath(out))
## Folder hasil: /Users/shofiandari/Library/CloudStorage/OneDrive-InstitutTeknologiSepuluhNopember/Dokumen/2024_work/narasumber/hasil_faktorial_dummy_R

12 Rujukan dan informasi sesi

sessionInfo()
## R version 4.6.0 (2026-04-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS 27.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Asia/Jakarta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39     R6_2.6.1          fastmap_1.2.0     xfun_0.59         cachem_1.1.0     
##  [6] knitr_1.52        htmltools_0.5.9   rmarkdown_2.32    lifecycle_1.0.5   cli_3.6.6        
## [11] sass_0.4.10       jquerylib_0.1.4   compiler_4.6.0    rstudioapi_0.19.0 tools_4.6.0      
## [16] evaluate_1.0.5    bslib_0.11.0      yaml_2.3.12       otel_0.2.0        jsonlite_2.0.0   
## [21] rlang_1.2.0