# NAMA : Elita Rohutami
# NIM : 2611018013
# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: Program intervensi Anemia Pada Remaja Putri di Wilayah Kerja Puskesmas(DATA SIMULASI)
#
# Isi:
# 0. Paket & pengaturan
# 1. Simulasi data (format panjang & lebar)
# 2. Eksplorasi data: statistik deskriptif, profile plot, spaghetti plot
# 3. Repeated Measure ANOVA satu arah (within-subject: Waktu)
# 3a. Uji asumsi: outlier, normalitas, sfierisitas (Mauchly)
# 3b. ANOVA + koreksi Greenhouse-Geisser / Huynh-Feldt
# 3c. Pendekatan multivariat (MANOVA) sebagai pembanding
# 3d. Post hoc berpasangan & kontras polinomial (tren)
# 3e. Alternatif nonparametrik: uji Friedman
# 4. Mixed Design ANOVA (between: Kelompok x within: Waktu)
# 4a. Uji asumsi: outlier, normalitas, Levene, Box's M, Mauchly
# 4b. ANOVA campuran + ukuran efek
# 4c. Analisis efek sederhana (simple effects) & post hoc
# 4d. Kontras interaksi (perubahan dari baseline antarkelompok)
# 5. Pembanding: Linear Mixed Model (LMM)
# 6. Menyimpan data & ringkasan hasil
# =============================================================================
# 0. PAKET & PENGATURAN
# Jalankan sekali bila paket belum terpasang:
# install.packages(c("tidyverse", "afex", "emmeans", "rstatix", "car",
# "effectsize", "lme4", "lmerTest", "performance", "ggpubr"))
suppressPackageStartupMessages({
library(dplyr) # manipulasi data
library(tidyr) # format panjang <-> lebar
library(ggplot2) # grafik
library(afex) # ANOVA within/mixed (tipe III, koreksi GG/HF, MANOVA)
library(emmeans) # rerata marginal, efek sederhana, post hoc, kontras
library(rstatix) # uji asumsi yang ramah pipe (outlier, Shapiro, Levene, Box's M)
library(car) # leveneTest, Anova
library(effectsize) # eta kuadrat parsial, omega kuadrat
library(lme4) # linear mixed model
library(lmerTest) # uji F/t dengan derajat bebas Satterthwaite / Kenward-Roger
})
options(contrasts = c("contr.sum", "contr.poly")) # kontras jumlah-nol untuk SS tipe III
afex_options(emmeans_model = "multivariate") # post hoc memakai model multivariat (tahan terhadap non-sfierisitas)
theme_set(theme_bw(base_size = 12))
# 1. SIMULASI DATA
# Skenario: uji klinis berkelompok di puskesmas. 90 pasien Anemia derajat 1
# diacak ke tiga kelompok (n = 30 per kelompok):
# - Kontrol : edukasi standar puskesmas
# - TTD : edukasi TTD (rendah garam, tinggi sayur-buah)
# - TTD+KIE : edukasi TTD + Komunikasi Informasi dan Edukasi terstruktur
# Anemia (Hb, mg/dl) diukur pada minggu ke-0 (baseline), 4, 8, 12.
#
# Struktur korelasi dibangkitkan dari intersep acak + kemiringan acak per subjek,
# sehingga varians meningkat dan korelasi menurun seiring jarak waktu --
# kondisi realistis yang cenderung MELANGGAR asumsi sfierisitas.
set.seed(2026)
n_per <- 30
kel_lab <- c("Kontrol", "TTD", "TTD+KIE")
minggu <- c(0, 4, 8, 12)
# Rerata populasi (mg/dl) per kelompok x waktu
mu <- rbind(
"Kontrol" = c(11, 11.1, 11, 11.2),
"TTD" = c(10.8, 11.1, 11.3, 11.5),
"TTD+KIE" = c(10.9, 11.3, 11.7, 11.9)
)
sd_int <- 9 # SD intersep acak (perbedaan TTD dasar antarpasien)
sd_slope <- 0.45 # SD kemiringan acak per minggu (perbedaan respons antarpasien)
sd_eps <- 4.5 # SD galat pengukuran
dat_wide <- lapply(seq_along(kel_lab), function(g) {
id <- (g - 1) * n_per + seq_len(n_per)
b0 <- rnorm(n_per, 0, sd_int)
b1 <- rnorm(n_per, 0, sd_slope)
y <- sapply(seq_along(minggu), function(t)
mu[g, t] + b0 + b1 * minggu[t] + rnorm(n_per, 0, sd_eps))
colnames(y) <- paste0("Hb_M", minggu)
data.frame(id = sprintf("P%03d", id),
kelompok = kel_lab[g],
usia = round(runif(n_per, 35, 65)),
jk = sample(c("L", "P"), n_per, replace = TRUE, prob = c(0, 1)),
round(y, 1))
}) |> bind_rows()
dat_wide$kelompok <- factor(dat_wide$kelompok, levels = kel_lab)
dat_wide$id <- factor(dat_wide$id)
# Format panjang (satu baris = satu pengukuran), dibutuhkan afex/rstatix/lme4
dat_long <- dat_wide |>
pivot_longer(starts_with("Hb_M"), names_to = "waktu", values_to = "hb") |>
mutate(waktu = factor(waktu, levels = paste0("Hb_M", minggu),
labels = paste0("M", minggu)),
minggu = as.numeric(sub("M", "", waktu)))
head(dat_wide)
## id kelompok usia jk Hb_M0 Hb_M4 Hb_M8 Hb_M12
## 1 P001 Kontrol 42 P 15.0 6.0 22.4 13.1
## 2 P002 Kontrol 43 P 2.0 8.3 11.5 3.9
## 3 P003 Kontrol 62 P 11.4 9.8 12.4 15.5
## 4 P004 Kontrol 40 P 8.9 15.4 15.5 14.5
## 5 P005 Kontrol 53 P -2.8 -2.8 -2.2 8.8
## 6 P006 Kontrol 42 P -12.9 -10.6 -5.3 0.0
head(dat_long)
## # A tibble: 6 × 7
## id kelompok usia jk waktu hb minggu
## <fct> <fct> <dbl> <chr> <fct> <dbl> <dbl>
## 1 P001 Kontrol 42 P M0 15 0
## 2 P001 Kontrol 42 P M4 6 4
## 3 P001 Kontrol 42 P M8 22.4 8
## 4 P001 Kontrol 42 P M12 13.1 12
## 5 P002 Kontrol 43 P M0 2 0
## 6 P002 Kontrol 43 P M4 8.3 4
str(dat_long)
## tibble [360 × 7] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 90 levels "P001","P002",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok: Factor w/ 3 levels "Kontrol","TTD",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num [1:360] 42 42 42 42 43 43 43 43 62 62 ...
## $ jk : chr [1:360] "P" "P" "P" "P" ...
## $ waktu : Factor w/ 4 levels "M0","M4","M8",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ hb : num [1:360] 15 6 22.4 13.1 2 8.3 11.5 3.9 11.4 9.8 ...
## $ minggu : num [1:360] 0 4 8 12 0 4 8 12 0 4 ...
# 2. EKSPLORASI DATA
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(hb, type = "mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 hb 30 10.5 11.9
## 2 Kontrol M4 hb 30 9.15 10.8
## 3 Kontrol M8 hb 30 11.2 9.49
## 4 Kontrol M12 hb 30 9.89 10.7
## 5 TTD M0 hb 30 12.2 9.48
## 6 TTD M4 hb 30 12.5 10.5
## 7 TTD M8 hb 30 12.3 9.69
## 8 TTD M12 hb 30 12.7 9.91
## 9 TTD+KIE M0 hb 30 11.4 11.8
## 10 TTD+KIE M4 hb 30 12.1 11.9
## 11 TTD+KIE M8 hb 30 12.9 12.8
## 12 TTD+KIE M12 hb 30 14.3 14.5
# Matriks kovarians & korelasi antarwaktu (seluruh subjek, dalam kelompok)
# -> memberi gambaran awal apakah sfierisitas masuk akal
S <- cov(dat_wide[, paste0("Hb_M", minggu)])
R <- cor(dat_wide[, paste0("Hb_M", minggu)])
round(S, 1); round(R, 2)
## Hb_M0 Hb_M4 Hb_M8 Hb_M12
## Hb_M0 120.8 96.8 90.3 89.2
## Hb_M4 96.8 122.4 97.7 102.0
## Hb_M8 90.3 97.7 114.3 105.8
## Hb_M12 89.2 102.0 105.8 140.9
## Hb_M0 Hb_M4 Hb_M8 Hb_M12
## Hb_M0 1.00 0.80 0.77 0.68
## Hb_M4 0.80 1.00 0.83 0.78
## Hb_M8 0.77 0.83 1.00 0.83
## Hb_M12 0.68 0.78 0.83 1.00
# Varians selisih antarpasangan waktu (inti asumsi sfierisitas)
pasangan <- combn(paste0("Hb_M", minggu), 2)
var_selisih <- apply(pasangan, 2, function(p) var(dat_wide[[p[1]]] - dat_wide[[p[2]]]))
names(var_selisih) <- apply(pasangan, 2, paste, collapse = " - ")
round(var_selisih, 1)
## Hb_M0 - Hb_M4 Hb_M0 - Hb_M8 Hb_M0 - Hb_M12 Hb_M4 - Hb_M8 Hb_M4 - Hb_M12
## 49.6 54.4 83.3 41.2 59.2
## Hb_M8 - Hb_M12
## 43.6
# Profile plot: rerata +/- 95% CI per kelompok
p_profil <- ggplot(dat_long, aes(minggu, hb, colour = kelompok, group = kelompok)) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", size = 2.5) +
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = .6) +
scale_x_continuous(breaks = minggu) +
labs(x = "Minggu ke-", y = "Anemia (mg/dl)", colour = "Kelompok",
title = "Profil rerata TTD (± 95% CI)") +
theme(legend.position = "bottom")
p_profil
# Spaghetti plot: lintasan tiap pasien
p_spag <- ggplot(dat_long, aes(minggu, hb, group = id)) +
geom_line(alpha = .3) +
stat_summary(aes(group = kelompok), fun = mean, geom = "line",
colour = "firebrick", linewidth = 1.2) +
facet_wrap(~ kelompok) +
scale_x_continuous(breaks = minggu) +
labs(x = "Minggu ke-", y = "Hb (mg/dl)", title = "Lintasan individu dan rerata kelompok")
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan: apakah TTD berubah selama 12 minggu pada kelompok TTD+KIE?
d1 <- droplevels(filter(dat_long, kelompok == "TTD+KIE"))
d1w <- filter(dat_wide, kelompok == "TTD+KIE")
## 3a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per waktu (ekstrem = di luar Q1-3IQR / Q3+3IQR)
d1 |> group_by(waktu) |> identify_outliers(hb)
## [1] waktu id kelompok usia jk hb minggu
## [8] is.outlier is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas per waktu (Shapiro-Wilk) dan Q-Q plot
d1 |> group_by(waktu) |> shapiro_test(hb)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 hb 0.959 0.295
## 2 M4 hb 0.953 0.204
## 3 M8 hb 0.982 0.869
## 4 M12 hb 0.960 0.310
ggpubr::ggqqplot(d1, "hb", facet.by = "waktu")
# (iii) Sfierisitas: Mauchly (dilaporkan otomatis oleh anova_test & afex)
aov1_rs <- anova_test(data = d1, dv = hb, wid = id, within = waktu,
effect.size = "pes")
aov1_rs # berisi: ANOVA, Mauchly's Test, koreksi GG & HF
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 87 1.377 0.255 0.045
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.533 0.004 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05
## 1 waktu 0.707 2.12, 61.53 0.26 0.765 2.29, 66.52 0.26
get_anova_table(aov1_rs, correction = "auto") # auto: GG dipakai jika Mauchly p < .05
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2.12 61.53 1.377 0.26 0.045
## 3b. ANOVA dengan afex (sumber utama laporan) -------------------------------
aov1 <- aov_ez(id = "id", dv = "hb", data = d1, within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov1 # tabel ringkas (df sudah dikoreksi GG)
## Anova Table (Type 3 tests)
##
## Response: hb
## Effect df MSE F ges pes p.value
## 1 waktu 2.12, 61.53 47.78 1.38 .007 .045 .260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov1) # univariat tanpa koreksi + Mauchly + epsilon GG & HF
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 19304.0 1 16042.6 29 34.8957 2.055e-06 ***
## waktu 139.6 3 2939.9 87 1.3767 0.2553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.53324 0.0037739
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.70726 0.2603
##
## HF eps Pr(>F[HF])
## waktu 0.7645997 0.2597175
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.05 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
omega_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## -------------------------------------------
## waktu | 1.94e-03 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
## 3c. Pendekatan multivariat (tidak memerlukan sfierisitas) ------------------
aov1$Anova # Pillai, Wilks, Hotelling-Lawley, Roy
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.54614 34.896 1 29 2.055e-06 ***
## waktu 1 0.07878 0.770 3 27 0.5211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 3d. Post hoc & kontras tren ------------------------------------------------
em1 <- emmeans(aov1, ~ waktu)
em1
## waktu emmean SE df lower.CL upper.CL
## M0 11.4 2.15 29 7.04 15.8
## M4 12.1 2.17 29 7.62 16.5
## M8 12.9 2.35 29 8.13 17.7
## M12 14.3 2.64 29 8.91 19.7
##
## Confidence level used: 0.95
pairs(em1, adjust = "bonferroni") # semua pasangan waktu (6 perbandingan)
## contrast estimate SE df t.ratio p.value
## M0 - M4 -0.630 1.44 29 -0.437 1.0000
## M0 - M8 -1.493 1.52 29 -0.984 1.0000
## M0 - M12 -2.877 2.05 29 -1.405 1.0000
## M4 - M8 -0.863 1.00 29 -0.861 1.0000
## M4 - M12 -2.247 1.52 29 -1.477 0.9034
## M8 - M12 -1.383 1.28 29 -1.085 1.0000
##
## P value adjustment: bonferroni method for 6 tests
contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm") # tiap waktu vs baseline
## contrast estimate SE df t.ratio p.value
## M4 - M0 0.63 1.44 29 0.437 0.6665
## M8 - M0 1.49 1.52 29 0.984 0.6665
## M12 - M0 2.88 2.05 29 1.405 0.5115
##
## P value adjustment: holm method for 3 tests
contrast(em1, "poly") # tren linear, kuadratik, kubik
## contrast estimate SE df t.ratio p.value
## linear 9.493 6.44 29 1.475 0.1511
## quadratic 0.753 1.77 29 0.426 0.6733
## cubic 0.287 3.24 29 0.089 0.9301
## 3e. Alternatif nonparametrik ------------------------------------------------
friedman_test(d1, hb ~ waktu | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 hb 30 2.04 3 0.564 Friedman test
friedman_effsize(d1, hb ~ waktu | id) # Kendall's W
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 hb 30 0.0227 Kendall W small
d1 |> wilcox_test(hb ~ waktu, paired = TRUE, p.adjust.method = "bonferroni")
## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 hb M0 M4 30 30 208. 0.630 1 ns
## 2 hb M0 M8 30 30 196. 0.455 1 ns
## 3 hb M0 M12 30 30 179 0.280 1 ns
## 4 hb M4 M8 30 30 186. 0.341 1 ns
## 5 hb M4 M12 30 30 167 0.182 1 ns
## 6 hb M8 M12 30 30 192. 0.419 1 ns
# (Opsional) ANOVA robust berbasis trimmed mean -- paket WRS2
if (requireNamespace("WRS2", quietly = TRUE)) {
print(WRS2::rmanova(d1$hb, d1$waktu, d1$id, tr = 0.2))
}
## Call:
## WRS2::rmanova(y = d1$hb, groups = d1$waktu, blocks = d1$id, tr = 0.2)
##
## Test statistic: F = 0.6974
## Degrees of freedom 1: 2.21
## Degrees of freedom 2: 37.51
## p-value: 0.51766
# 4. MIXED DESIGN ANOVA (Kelompok [between] x Waktu [within])
# Pertanyaan: apakah pola perubahan TTD berbeda antar kelompok intervensi?
## 4a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per sel
dat_long |> group_by(kelompok, waktu) |> identify_outliers(hb)
## # A tibble: 2 × 9
## kelompok waktu id usia jk hb minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol M8 P015 57 P -12.9 8 TRUE FALSE
## 2 Kontrol M12 P015 57 P -18 12 TRUE FALSE
# (ii) Normalitas per sel (3 x 4 = 12 sel) dan residual model
dat_long |> group_by(kelompok, waktu) |> shapiro_test(hb)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 hb 0.990 0.991
## 2 Kontrol M4 hb 0.963 0.366
## 3 Kontrol M8 hb 0.972 0.599
## 4 Kontrol M12 hb 0.976 0.716
## 5 TTD M0 hb 0.961 0.328
## 6 TTD M4 hb 0.980 0.819
## 7 TTD M8 hb 0.968 0.488
## 8 TTD M12 hb 0.974 0.646
## 9 TTD+KIE M0 hb 0.959 0.295
## 10 TTD+KIE M4 hb 0.953 0.204
## 11 TTD+KIE M8 hb 0.982 0.869
## 12 TTD+KIE M12 hb 0.960 0.310
ggpubr::ggqqplot(dat_long, "hb", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antarkelompok pada TIAP waktu (Levene, median-centered)
dat_long |> group_by(waktu) |> levene_test(hb ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 0.844 0.434
## 2 M4 2 87 0.672 0.514
## 3 M8 2 87 2.04 0.137
## 4 M12 2 87 4.76 0.0109
# (iv) Homogenitas matriks kovarians antarkelompok (Box's M; uji pada alpha = .001)
box_m(dat_wide[, paste0("Hb_M", minggu)], dat_wide$kelompok)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 19.6 0.486 20 Box's M-test for Homogeneity of Covariance Matric…
# (v) Sfierisitas: Mauchly (dari summary model di bawah)
## 4b. ANOVA campuran ---------------------------------------------------------
aov2 <- aov_ez(id = "id", dv = "hb", data = dat_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: hb
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 419.79 0.54 .010 .012 .583
## 2 waktu 2.60, 226.52 31.95 0.90 .002 .010 .429
## 3 kelompok:waktu 5.21, 226.52 31.95 0.81 .003 .018 .545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2) # Mauchly, epsilon GG/HF, p terkoreksi
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 49815 1 36522 87 118.6661 <2e-16 ***
## kelompok 457 2 36522 87 0.5438 0.5825
## waktu 75 3 7238 261 0.9032 0.4401
## kelompok:waktu 135 6 7238 261 0.8135 0.5602
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.80865 0.0026998
## kelompok:waktu 0.80865 0.0026998
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.86788 0.4288
## kelompok:waktu 0.86788 0.5454
##
## HF eps Pr(>F[HF])
## waktu 0.8970366 0.4314913
## kelompok:waktu 0.8970366 0.5488237
aov2$Anova # uji multivariat untuk efek within & interaksi
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.57698 118.666 1 87 <2e-16 ***
## kelompok 2 0.01235 0.544 2 87 0.5825
## waktu 1 0.02567 0.746 3 85 0.5274
## kelompok:waktu 2 0.05205 0.766 6 172 0.5976
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix (hasil identik; format ringkas untuk laporan)
aov2_rs <- anova_test(data = dat_long, dv = hb, wid = id,
between = kelompok, within = waktu, effect.size = "pes",
type = 3)
get_anova_table(aov2_rs, correction = "GG")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 kelompok 2.00 87.00 0.544 0.583 0.012
## 2 waktu 2.60 226.52 0.903 0.429 0.010
## 3 kelompok:waktu 5.21 226.52 0.814 0.545 0.018
# Ukuran efek
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.01 | [0.00, 1.00]
## waktu | 0.01 | [0.00, 1.00]
## kelompok:waktu | 0.02 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
omega_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## ------------------------------------------------
## kelompok | 0 | [0.00, 1.00]
## waktu | 0 | [0.00, 1.00]
## kelompok:waktu | 0 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Plot interaksi dari model
afex_plot(aov2, x = "waktu", trace = "kelompok", error = "within",
mapping = c("colour", "shape", "linetype")) +
labs(y = "TTD (mg/dl)", x = "Waktu")
## Warning: Panel(s) show a mixed within-between-design.
## Error bars do not allow comparisons across all means.
## Suppress error bars with: error = "none"
## 4c. Efek sederhana (karena interaksi signifikan) ---------------------------
em2 <- emmeans(aov2, ~ waktu | kelompok)
# Efek WAKTU di dalam tiap kelompok (uji F gabungan per kelompok)
joint_tests(aov2, by = "kelompok")
## kelompok = Kontrol:
## model term df1 df2 F.ratio p.value
## waktu 3 87 1.160 0.3299
##
## kelompok = TTD:
## model term df1 df2 F.ratio p.value
## waktu 3 87 0.046 0.9870
##
## kelompok = TTD+KIE:
## model term df1 df2 F.ratio p.value
## waktu 3 87 1.119 0.3460
# Efek KELOMPOK pada tiap waktu
joint_tests(aov2, by = "waktu")
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.177 0.8384
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.816 0.4455
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.205 0.8149
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 1.067 0.3486
# Post hoc: tiap waktu vs baseline di dalam tiap kelompok
contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm")
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -1.353 1.29 87 -1.048 0.8920
## M8 - M0 0.667 1.36 87 0.491 1.0000
## M12 - M0 -0.617 1.66 87 -0.371 1.0000
##
## kelompok = TTD:
## contrast estimate SE df t.ratio p.value
## M4 - M0 0.317 1.29 87 0.245 1.0000
## M8 - M0 0.113 1.36 87 0.083 1.0000
## M12 - M0 0.473 1.66 87 0.285 1.0000
##
## kelompok = TTD+KIE:
## contrast estimate SE df t.ratio p.value
## M4 - M0 0.630 1.29 87 0.488 0.6267
## M8 - M0 1.493 1.36 87 1.100 0.5491
## M12 - M0 2.877 1.66 87 1.730 0.2618
##
## P value adjustment: holm method for 3 tests
# Post hoc: perbandingan antarkelompok pada tiap waktu
em2b <- emmeans(aov2, ~ kelompok | waktu)
pairs(em2b, adjust = "tukey") # Tukey per waktu
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Kontrol - TTD -1.700 2.86 87 -0.593 0.8240
## Kontrol - (TTD+KIE) -0.930 2.86 87 -0.325 0.9436
## TTD - (TTD+KIE) 0.770 2.86 87 0.269 0.9610
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - TTD -3.370 2.86 87 -1.177 0.4697
## Kontrol - (TTD+KIE) -2.913 2.86 87 -1.018 0.5676
## TTD - (TTD+KIE) 0.457 2.86 87 0.160 0.9861
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - TTD -1.147 2.78 87 -0.412 0.9109
## Kontrol - (TTD+KIE) -1.757 2.78 87 -0.631 0.8036
## TTD - (TTD+KIE) -0.610 2.78 87 -0.219 0.9739
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - TTD -2.790 3.06 87 -0.911 0.6349
## Kontrol - (TTD+KIE) -4.423 3.06 87 -1.444 0.3230
## TTD - (TTD+KIE) -1.633 3.06 87 -0.533 0.8552
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## 4d. Kontras interaksi -------------------------------------------------------
# Apakah penurunan (M12 - M0) berbeda antarkelompok? -- inti pertanyaan uji klinis
em_full <- emmeans(aov2, ~ waktu * kelompok)
contrast(em_full, interaction = list(waktu = list("M12-M0" = c(-1, 0, 0, 1)),
kelompok = "pairwise"),
adjust = "holm")
## waktu_custom kelompok_pairwise estimate SE df t.ratio p.value
## M12-M0 Kontrol - TTD -1.09 2.35 87 -0.463 0.6442
## M12-M0 Kontrol - (TTD+KIE) -3.49 2.35 87 -1.485 0.4234
## M12-M0 TTD - (TTD+KIE) -2.40 2.35 87 -1.022 0.6195
##
## P value adjustment: holm method for 3 tests
# Tren linear per kelompok dan perbandingannya
contrast(em2, "poly")[c(1, 4, 7)]
## contrast kelompok estimate SE df t.ratio p.value
## linear Kontrol 0.17 5.33 87 0.032 0.9746
## linear TTD 1.22 5.33 87 0.228 0.8201
## linear TTD+KIE 9.49 5.33 87 1.780 0.0786
tren_int <- summary(contrast(em_full, interaction = c(waktu = "poly", kelompok = "pairwise"),
adjust = "none"))
tren_lin <- subset(tren_int, waktu_poly == "linear") # apakah laju penurunan linear berbeda?
tren_lin$p.holm <- p.adjust(tren_lin$p.value, "holm") # koreksi Holm untuk 3 perbandingan
tren_lin
## waktu_poly kelompok_pairwise estimate SE df t.ratio p.value
## 1 linear Kontrol - TTD -1.046667 7.542958 87 -0.1387608 0.8899599
## 4 linear Kontrol - (TTD+KIE) -9.323333 7.542958 87 -1.2360314 0.2197737
## 7 linear TTD - (TTD+KIE) -8.276667 7.542958 87 -1.0972707 0.2755509
## p.holm
## 1 0.8899599
## 4 0.6593212
## 7 0.6593212
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Tidak mensyaratkan sfierisitas, menampung data hilang (MAR) dan waktu
# pengukuran yang tidak seragam.
lmm1 <- lmer(hb ~ kelompok * waktu + (1 | id), data = dat_long, REML = TRUE)
lmm2 <- lmer(hb ~ kelompok * waktu + (1 + minggu | id), data = dat_long, REML = TRUE)
anova(lmm1, lmm2, refit = FALSE) # uji rasio kemungkinan struktur acak
## Data: dat_long
## Models:
## lmm1: hb ~ kelompok * waktu + (1 | id)
## lmm2: hb ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 2466.2 2520.6 -1219.1 2438.2
## lmm2 16 2452.4 2514.6 -1210.2 2420.4 17.742 2 0.0001404 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm2, ddf = "Kenward-Roger") # uji F tipe III efek tetap
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 22.038 11.019 2 87.00 0.5438 0.5825
## waktu 44.064 14.688 3 185.37 0.7215 0.5403
## kelompok:waktu 94.263 15.710 6 206.40 0.7709 0.5936
performance::icc(lmm1) # korelasi intrakelas
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.779
## Unadjusted ICC: 0.768
# Diagnostik residual LMM (normalitas & homogenitas)
par(mfrow = c(1, 3))
qqnorm(resid(lmm2), main = "Q-Q residual"); qqline(resid(lmm2))
qqnorm(ranef(lmm2)$id[, 1], main = "Q-Q intersep acak"); qqline(ranef(lmm2)$id[, 1])
plot(fitted(lmm2), resid(lmm2), xlab = "Nilai prediksi", ylab = "Residual",
main = "Residual vs prediksi"); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# (paket 'see' + performance::check_model(lmm2) memberi panel diagnostik lengkap)
# Simulasi 30 nilai hilang (MCAR, +/- 11% pengukuran pasca-baseline)
# untuk menunjukkan keunggulan LMM
set.seed(1)
dat_miss <- dat_long
dat_miss$hb[sample(which(dat_miss$waktu != "M0"), 30)] <- NA
lmm_miss <- lmer(hb ~ kelompok * waktu + (1 + minggu | id), data = dat_miss,
control = lmerControl(optimizer = "bobyqa"))
anova(lmm_miss, ddf = "Kenward-Roger") # semua pasien tetap dianalisis
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 26.502 13.251 2 86.984 0.6371 0.5313
## waktu 43.991 14.664 3 168.209 0.7017 0.5523
## kelompok:waktu 152.739 25.456 6 185.860 1.2167 0.2995
# RM ANOVA akan membuang seluruh pasien yang punya >= 1 nilai hilang:
n_distinct(dat_miss$id[is.na(dat_miss$hb)])
## [1] 25
# 6. SIMPAN DATA & SESSION INFO
write.csv(dat_wide, "data_hb_anemia_wide.csv", row.names = FALSE)
write.csv(dat_long, "data_hb_anemia_long.csv", row.names = FALSE)
sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Asia/Makassar
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lmerTest_3.2-1 effectsize_1.0.3 car_3.1-5 carData_3.0-6
## [5] rstatix_1.1.0 emmeans_2.0.4 afex_1.5-1 lme4_2.0-6
## [9] Matrix_1.7-5 ggplot2_4.0.3 tidyr_1.3.2 dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 farver_2.1.2 S7_0.2.2
## [4] fastmap_1.2.0 reshape_0.8.10 bayestestR_0.19.0
## [7] digest_0.6.39 rpart_4.1.27 estimability_2.0.0
## [10] lifecycle_1.0.5 cluster_2.1.8.2 magrittr_2.0.5
## [13] compiler_4.6.1 rlang_1.3.0 Hmisc_5.3-0
## [16] sass_0.4.10 tools_4.6.1 utf8_1.2.6
## [19] yaml_2.3.12 data.table_1.18.6.1 ggsignif_0.6.4
## [22] knitr_1.51 labeling_0.4.3 htmlwidgets_1.6.4
## [25] plyr_1.8.9 RColorBrewer_1.1-3 abind_1.4-8
## [28] withr_3.0.3 foreign_0.8-91 purrr_1.2.2
## [31] numDeriv_2016.8-1.1 nnet_7.3-20 grid_4.6.1
## [34] datawizard_1.4.0 ggpubr_1.0.0 colorspace_2.1-3
## [37] scales_1.4.0 MASS_7.3-65 insight_1.5.4
## [40] cli_3.6.6 mvtnorm_1.4-2 rmarkdown_2.31
## [43] reformulas_0.4.4 generics_0.1.4 performance_0.18.2
## [46] rstudioapi_0.19.0 reshape2_1.4.5 parameters_0.29.3
## [49] minqa_1.2.8 cachem_1.1.0 stringr_1.6.0
## [52] splines_4.6.1 parallel_4.6.1 WRS2_1.1-7
## [55] base64enc_0.1-6 vctrs_0.7.3 boot_1.3-32
## [58] jsonlite_2.0.0 pbkrtest_0.5.5 Formula_1.2-6
## [61] htmlTable_2.5.0 jquerylib_0.1.4 glue_1.8.1
## [64] nloptr_2.2.1 stringi_1.8.9 gtable_0.3.6
## [67] tibble_3.3.1 pillar_1.11.1 htmltools_0.5.9
## [70] R6_2.6.1 Rdpack_2.6.6 evaluate_1.0.5
## [73] lattice_0.22-9 rbibutils_2.4.1 backports_1.5.1
## [76] broom_1.0.13 bslib_0.12.0 Rcpp_1.1.2
## [79] gridExtra_2.3.1 nlme_3.1-169 checkmate_2.3.4
## [82] xfun_0.60 pkgconfig_2.0.3