# Nama : Desy Ulfayanti Kalauw
# NIM : 2611018046
# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: Program intervensi Diabetes Melitus dan Gula Darah Sewaktu
# pada pasien Diabetes Melitus Tipe 2 di 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 & session info
# =============================================================================
# 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
library(emmeans) # rerata marginal, post hoc, kontras
library(rstatix) # uji asumsi
library(car) # Levene/Anova
library(effectsize) # ukuran efek
library(lme4) # linear mixed model
library(lmerTest) # uji F/t LMM
})
options(contrasts = c("contr.sum", "contr.poly"))
afex_options(emmeans_model = "multivariate")
theme_set(theme_bw(base_size = 12))
# 1. SIMULASI DATA
# Skenario:
# 90 pasien Diabetes Melitus Tipe 2 dibagi menjadi tiga kelompok
# (n = 30 per kelompok):
# - Kontrol : edukasi standar pengelolaan diabetes
# - Diet : edukasi diet diabetes terstruktur
# - Diet+Aktivitas: edukasi diet + aktivitas fisik terstruktur
#
# Gula Darah Sewaktu (GDS, mg/dL) diukur pada minggu ke-0, 4, 8, dan 12.
# Nilai dibuat realistis sebagai data simulasi untuk latihan analisis statistik.
set.seed(2026)
n_per <- 30
kel_lab <- c("Kontrol", "Diet", "Diet+Aktivitas")
minggu <- c(0, 4, 8, 12)
# Rerata populasi GDS (mg/dL) per kelompok x waktu
# Semakin efektif intervensi, semakin besar penurunan GDS.
mu <- rbind(
"Kontrol" = c(245, 240, 237, 235),
"Diet" = c(245, 225, 210, 198),
"Diet+Aktivitas" = c(245, 215, 195, 178)
)
# Komponen variasi individual
sd_int <- 25 # variasi kadar GDS dasar antarpasien
sd_slope <- 1.20 # variasi respons perubahan per minggu
sd_eps <- 14 # galat pengukuran
# Membentuk data format wide
dat_wide <- lapply(seq_along(kel_lab), function(g) {
id <- (g - 1) * n_per + seq_len(n_per)
# efek acak individu
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("GDS_M", minggu)
data.frame(
id = sprintf("DM%03d", id),
kelompok = kel_lab[g],
usia = round(runif(n_per, 40, 70)),
jk = sample(c("L", "P"), n_per, replace = TRUE, prob = c(.45, .55)),
lama_dm = round(runif(n_per, 1, 12), 1),
round(y, 1)
)
}) |> bind_rows()
# Mengubah tipe variabel menjadi faktor
dat_wide$kelompok <- factor(dat_wide$kelompok, levels = kel_lab)
dat_wide$id <- factor(dat_wide$id)
# Mengubah data wide menjadi format long
# Satu baris = satu pengukuran pada satu pasien
dat_long <- dat_wide |>
pivot_longer(
starts_with("GDS_M"),
names_to = "waktu",
values_to = "gds"
) |>
mutate(
waktu = factor(
waktu,
levels = paste0("GDS_M", minggu),
labels = paste0("M", minggu)
),
minggu = as.numeric(sub("M", "", waktu))
)
# Melihat struktur data
head(dat_wide)
## id kelompok usia jk lama_dm GDS_M0 GDS_M4 GDS_M8 GDS_M12
## 1 DM001 Kontrol 47 L 3.1 256.0 223.1 271.6 240.7
## 2 DM002 Kontrol 48 L 1.6 220.4 234.4 241.5 215.1
## 3 DM003 Kontrol 67 P 9.6 245.7 235.6 241.3 248.7
## 4 DM004 Kontrol 45 P 5.7 238.8 253.2 250.6 244.7
## 5 DM005 Kontrol 58 P 4.0 204.1 198.7 197.7 229.0
## 6 DM006 Kontrol 47 L 4.2 178.3 179.0 191.8 204.9
head(dat_long)
## # A tibble: 6 × 8
## id kelompok usia jk lama_dm waktu gds minggu
## <fct> <fct> <dbl> <chr> <dbl> <fct> <dbl> <dbl>
## 1 DM001 Kontrol 47 L 3.1 M0 256 0
## 2 DM001 Kontrol 47 L 3.1 M4 223. 4
## 3 DM001 Kontrol 47 L 3.1 M8 272. 8
## 4 DM001 Kontrol 47 L 3.1 M12 241. 12
## 5 DM002 Kontrol 48 L 1.6 M0 220. 0
## 6 DM002 Kontrol 48 L 1.6 M4 234. 4
str(dat_long)
## tibble [360 × 8] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 90 levels "DM001","DM002",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok: Factor w/ 3 levels "Kontrol","Diet",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num [1:360] 47 47 47 47 48 48 48 48 67 67 ...
## $ jk : chr [1:360] "L" "L" "L" "L" ...
## $ lama_dm : num [1:360] 3.1 3.1 3.1 3.1 1.6 1.6 1.6 1.6 9.6 9.6 ...
## $ waktu : Factor w/ 4 levels "M0","M4","M8",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ gds : num [1:360] 256 223 272 241 220 ...
## $ minggu : num [1:360] 0 4 8 12 0 4 8 12 0 4 ...
# 2. EKSPLORASI DATA
# Statistik deskriptif GDS: mean dan SD per kelompok dan waktu
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(gds, type = "mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 gds 30 244. 34.0
## 2 Kontrol M4 gds 30 234. 30.7
## 3 Kontrol M8 gds 30 238. 26.5
## 4 Kontrol M12 gds 30 231. 29.8
## 5 Diet M0 gds 30 250. 24.2
## 6 Diet M4 gds 30 231. 24.6
## 7 Diet M8 gds 30 216. 26.5
## 8 Diet M12 gds 30 201. 28.2
## 9 Diet+Aktivitas M0 gds 30 254. 38.7
## 10 Diet+Aktivitas M4 gds 30 225. 33.5
## 11 Diet+Aktivitas M8 gds 30 207. 37.6
## 12 Diet+Aktivitas M12 gds 30 188. 36.4
# Matriks kovarians dan korelasi antarwaktu
S <- cov(dat_wide[, paste0("GDS_M", minggu)])
R <- cor(dat_wide[, paste0("GDS_M", minggu)])
round(S, 1)
## GDS_M0 GDS_M4 GDS_M8 GDS_M12
## GDS_M0 1072.5 762.6 643.3 620.6
## GDS_M4 762.6 885.5 731.0 814.6
## GDS_M8 643.3 731.0 1085.6 974.7
## GDS_M12 620.6 814.6 974.7 1317.2
round(R, 2)
## GDS_M0 GDS_M4 GDS_M8 GDS_M12
## GDS_M0 1.00 0.78 0.60 0.52
## GDS_M4 0.78 1.00 0.75 0.75
## GDS_M8 0.60 0.75 1.00 0.82
## GDS_M12 0.52 0.75 0.82 1.00
# Varians selisih antarpasangan waktu
# Digunakan sebagai gambaran awal asumsi sfierisitas
pasangan <- combn(paste0("GDS_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)
## GDS_M0 - GDS_M4 GDS_M0 - GDS_M8 GDS_M0 - GDS_M12 GDS_M4 - GDS_M8
## 432.8 871.6 1148.4 509.1
## GDS_M4 - GDS_M12 GDS_M8 - GDS_M12
## 573.5 453.5
# Profile plot: rerata +/- 95% CI per kelompok
p_profil <- ggplot(
dat_long,
aes(x = minggu, y = gds, 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 = "Gula Darah Sewaktu (mg/dL)",
colour = "Kelompok",
title = "Profil Rerata GDS (+/- 95% CI)"
) +
theme(legend.position = "bottom")
p_profil
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
# Spaghetti plot: lintasan GDS tiap pasien
p_spag <- ggplot(dat_long, aes(x = minggu, y = gds, 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 = "GDS (mg/dL)",
title = "Lintasan Individu dan Rerata GDS per Kelompok"
)
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan:
# Apakah GDS berubah selama 12 minggu pada kelompok Diet+Aktivitas?
d1 <- droplevels(filter(dat_long, kelompok == "Diet+Aktivitas"))
d1w <- filter(dat_wide, kelompok == "Diet+Aktivitas")
# 3a. UJI ASUMSI ---------------------------------------------------------------
# (i) Outlier per waktu
# Outlier ekstrem = nilai di luar Q1-3IQR atau Q3+3IQR
d1 |>
group_by(waktu) |>
identify_outliers(gds)
## # A tibble: 1 × 10
## waktu id kelompok usia jk lama_dm gds minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <chr> <dbl> <dbl> <dbl> <lgl> <lgl>
## 1 M0 DM082 Diet+Aktiv… 57 L 3.1 156. 0 TRUE FALSE
# (ii) Normalitas per waktu dengan Shapiro-Wilk
d1 |>
group_by(waktu) |>
shapiro_test(gds)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 gds 0.967 0.462
## 2 M4 gds 0.932 0.0545
## 3 M8 gds 0.953 0.200
## 4 M12 gds 0.931 0.0518
# Q-Q plot
ggpubr::ggqqplot(d1, "gds", facet.by = "waktu")
# (iii) Uji sfierisitas Mauchly
# anova_test akan menghasilkan ANOVA, Mauchly, serta koreksi GG/HF
aov1_rs <- anova_test(
data = d1,
dv = gds,
wid = id,
within = waktu,
effect.size = "pes"
)
aov1_rs
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 87 110.797 1.29e-29 * 0.793
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.752 0.162
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.852 2.55, 74.08 1.42e-25 * 0.941 2.82, 81.85 5.27e-28
## p[HF]<.05
## 1 *
# Secara otomatis memakai koreksi GG bila Mauchly p < 0.05
get_anova_table(aov1_rs, correction = "auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 87 110.797 1.29e-29 * 0.793
# 3b. REPEATED MEASURE ANOVA DENGAN AFEX --------------------------------------
aov1 <- aov_ez(
id = "id",
dv = "gds",
data = d1,
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
# Tabel utama repeated measure ANOVA
aov1
## Anova Table (Type 3 tests)
##
## Response: gds
## Effect df MSE F ges pes p.value
## 1 waktu 2.55, 74.08 249.91 110.80 *** .313 .793 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Output lengkap termasuk Mauchly dan epsilon GG/HF
summary(aov1)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 5726973 1 136735 29 1214.6 < 2.2e-16 ***
## waktu 70737 3 18515 87 110.8 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.75214 0.16229
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.85154 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.9407933 5.273886e-28
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.79 | [0.73, 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 | 0.30 | [0.16, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 3c. PENDEKATAN MULTIVARIAT (MANOVA) ----------------------------------------
# Tidak memerlukan asumsi sfierisitas
# Menampilkan Pillai, Wilks, Hotelling-Lawley, dan Roy
aov1$Anova
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.97668 1214.62 1 29 < 2.2e-16 ***
## waktu 1 0.88557 69.65 3 27 7.857e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3d. POST HOC DAN KONTRAS TREN ------------------------------------------------
# Estimated Marginal Means tiap waktu
em1 <- emmeans(aov1, ~ waktu)
em1
## waktu emmean SE df lower.CL upper.CL
## M0 254 7.06 29 240 268
## M4 225 6.11 29 212 237
## M8 207 6.86 29 193 221
## M12 188 6.65 29 174 201
##
## Confidence level used: 0.95
# Semua pasangan waktu
pairs(em1, adjust = "bonferroni")
## contrast estimate SE df t.ratio p.value
## M0 - M4 29.4 3.23 29 9.102 <0.0001
## M0 - M8 46.6 4.29 29 10.862 <0.0001
## M0 - M12 66.1 4.46 29 14.818 <0.0001
## M4 - M8 17.2 3.61 29 4.758 0.0003
## M4 - M12 36.7 3.25 29 11.298 <0.0001
## M8 - M12 19.5 3.57 29 5.461 <0.0001
##
## P value adjustment: bonferroni method for 6 tests
# Membandingkan setiap waktu dengan baseline (M0)
contrast(
em1,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## contrast estimate SE df t.ratio p.value
## M4 - M0 -29.4 3.23 29 -9.102 <0.0001
## M8 - M0 -46.6 4.29 29 -10.862 <0.0001
## M12 - M0 -66.1 4.46 29 -14.818 <0.0001
##
## P value adjustment: holm method for 3 tests
# Tren linear, kuadratik, dan kubik
contrast(em1, "poly")
## contrast estimate SE df t.ratio p.value
## linear -215.53 14.50 29 -14.890 <0.0001
## quadratic 9.92 4.38 29 2.264 0.0312
## cubic -14.53 11.00 29 -1.326 0.1952
# 3e. ALTERNATIF NONPARAMETRIK -------------------------------------------------
# Uji Friedman
friedman_test(d1, gds ~ waktu | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 gds 30 74.5 3 4.59e-16 Friedman test
# Ukuran efek Kendall's W
friedman_effsize(d1, gds ~ waktu | id)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 gds 30 0.828 Kendall W large
# Post hoc Wilcoxon berpasangan
d1 |>
wilcox_test(
gds ~ 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 gds M0 M4 30 30 458 0.0000000354 2.12e-7 ****
## 2 gds M0 M8 30 30 465 0.00000000186 1.12e-8 ****
## 3 gds M0 M12 30 30 465 0.00000000186 1.12e-8 ****
## 4 gds M4 M8 30 30 416 0.0000493 2.96e-4 ***
## 5 gds M4 M12 30 30 465 0.00000000186 1.12e-8 ****
## 6 gds M8 M12 30 30 430 0.00000799 4.80e-5 ****
# ANOVA robust berbasis trimmed mean (opsional)
if (requireNamespace("WRS2", quietly = TRUE)) {
print(WRS2::rmanova(d1$gds, d1$waktu, d1$id, tr = 0.2))
}
## Call:
## WRS2::rmanova(y = d1$gds, groups = d1$waktu, blocks = d1$id,
## tr = 0.2)
##
## Test statistic: F = 66.3848
## Degrees of freedom 1: 2.85
## Degrees of freedom 2: 48.52
## p-value: 0
# 4. MIXED DESIGN ANOVA
# Between-subject = Kelompok
# Within-subject = Waktu
#
# Pertanyaan:
# Apakah pola perubahan GDS selama 12 minggu berbeda antarkelompok?
# 4a. UJI ASUMSI ---------------------------------------------------------------
# (i) Outlier per kelompok x waktu
dat_long |>
group_by(kelompok, waktu) |>
identify_outliers(gds)
## # A tibble: 7 × 10
## kelompok waktu id usia jk lama_dm gds minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <chr> <dbl> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol M0 DM030 55 P 1.2 322. 0 TRUE FALSE
## 2 Kontrol M8 DM015 62 P 10.4 172. 8 TRUE FALSE
## 3 Kontrol M12 DM015 62 P 10.4 154. 12 TRUE FALSE
## 4 Diet M0 DM042 47 P 8.9 183. 0 TRUE FALSE
## 5 Diet M8 DM042 47 P 8.9 156. 8 TRUE FALSE
## 6 Diet M8 DM060 58 L 6.4 283. 8 TRUE FALSE
## 7 Diet+Aktiv… M0 DM082 57 L 3.1 156. 0 TRUE FALSE
# (ii) Normalitas setiap sel (3 kelompok x 4 waktu = 12 sel)
dat_long |>
group_by(kelompok, waktu) |>
shapiro_test(gds)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 gds 0.987 0.971
## 2 Kontrol M4 gds 0.963 0.361
## 3 Kontrol M8 gds 0.969 0.511
## 4 Kontrol M12 gds 0.974 0.643
## 5 Diet M0 gds 0.966 0.444
## 6 Diet M4 gds 0.972 0.594
## 7 Diet M8 gds 0.972 0.609
## 8 Diet M12 gds 0.950 0.165
## 9 Diet+Aktivitas M0 gds 0.967 0.462
## 10 Diet+Aktivitas M4 gds 0.932 0.0545
## 11 Diet+Aktivitas M8 gds 0.953 0.200
## 12 Diet+Aktivitas M12 gds 0.931 0.0518
# Q-Q plot per kelompok dan waktu
ggpubr::ggqqplot(dat_long, "gds", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antarkelompok pada setiap waktu
# Levene test
dat_long |>
group_by(waktu) |>
levene_test(gds ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 2.27 0.109
## 2 M4 2 87 0.939 0.395
## 3 M8 2 87 2.26 0.111
## 4 M12 2 87 1.47 0.237
# (iv) Homogenitas matriks kovarians antarkelompok
# Box's M biasanya dievaluasi secara konservatif, misalnya alpha 0.001
box_m(
dat_wide[, paste0("GDS_M", minggu)],
dat_wide$kelompok
)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 14.9 0.783 20 Box's M-test for Homogeneity of Covariance Matric…
# (v) Sfierisitas diperoleh dari summary model mixed ANOVA di bawah
# 4b. MIXED DESIGN ANOVA -------------------------------------------------------
aov2 <- aov_ez(
id = "id",
dv = "gds",
data = dat_long,
between = "kelompok",
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
# Tabel ANOVA campuran
aov2
## Anova Table (Type 3 tests)
##
## Response: gds
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 3199.10 3.29 * .058 .070 .042
## 2 waktu 2.65, 230.91 268.15 120.52 *** .201 .581 <.001
## 3 kelompok:waktu 5.31, 230.91 268.15 18.85 *** .073 .302 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Mauchly, epsilon GG/HF, dan output lengkap
summary(aov2)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 18489609 1 278322 87 5779.6305 < 2e-16 ***
## kelompok 21034 2 278322 87 3.2876 0.04203 *
## waktu 85776 3 61920 261 120.5198 < 2e-16 ***
## kelompok:waktu 26834 6 61920 261 18.8515 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.81372 0.0033904
## kelompok:waktu 0.81372 0.0033904
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.88473 < 2.2e-16 ***
## kelompok:waktu 0.88473 2.233e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.9151557 4.470672e-45
## kelompok:waktu 0.9151557 7.305154e-17
# Pendekatan multivariat
aov2$Anova
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.98517 5779.6 1 87 < 2.2e-16 ***
## kelompok 2 0.07027 3.3 2 87 0.04203 *
## waktu 1 0.74437 82.5 3 85 < 2.2e-16 ***
## kelompok:waktu 2 0.48960 9.3 6 172 8.03e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix
aov2_rs <- anova_test(
data = dat_long,
dv = gds,
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 3.288 4.20e-02 * 0.070
## 2 waktu 2.65 230.91 120.520 1.14e-43 * 0.581
## 3 kelompok:waktu 5.31 230.91 18.852 2.23e-16 * 0.302
# Ukuran efek
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.07 | [0.00, 1.00]
## waktu | 0.58 | [0.52, 1.00]
## kelompok:waktu | 0.30 | [0.22, 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.05 | [0.00, 1.00]
## waktu | 0.20 | [0.13, 1.00]
## kelompok:waktu | 0.07 | [0.01, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Plot interaksi kelompok x waktu
afex_plot(
aov2,
x = "waktu",
trace = "kelompok",
error = "within",
mapping = c("colour", "shape", "linetype")
) +
labs(
y = "GDS (mg/dL)",
x = "Waktu Pengukuran"
)
## 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 (SIMPLE EFFECTS) -----------------------------------------
# Digunakan terutama bila interaksi kelompok x waktu signifikan
# Estimated marginal means waktu dalam setiap kelompok
em2 <- emmeans(aov2, ~ waktu | kelompok)
em2
## kelompok = Kontrol:
## waktu emmean SE df lower.CL upper.CL
## M0 244 5.99 87 232 256
## M4 234 5.44 87 223 245
## M8 238 5.60 87 227 249
## M12 231 5.78 87 220 243
##
## kelompok = Diet:
## waktu emmean SE df lower.CL upper.CL
## M0 250 5.99 87 238 262
## M4 231 5.44 87 220 241
## M8 216 5.60 87 205 227
## M12 201 5.78 87 190 213
##
## kelompok = Diet+Aktivitas:
## waktu emmean SE df lower.CL upper.CL
## M0 254 5.99 87 242 266
## M4 225 5.44 87 214 235
## M8 207 5.60 87 196 219
## M12 188 5.78 87 176 199
##
## Confidence level used: 0.95
# Efek waktu di dalam masing-masing kelompok
joint_tests(aov2, by = "kelompok")
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## kelompok = Kontrol:
## model term df1 df2 F.ratio p.value
## waktu 3 87 3.025 0.0338
##
## kelompok = Diet:
## model term df1 df2 F.ratio p.value
## waktu 3 87 38.975 <0.0001
##
## kelompok = Diet+Aktivitas:
## model term df1 df2 F.ratio p.value
## waktu 3 87 69.360 <0.0001
# Efek kelompok pada masing-masing waktu
joint_tests(aov2, by = "waktu")
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.764 0.4687
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.808 0.4490
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 7.919 0.0007
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 14.895 <0.0001
# Setiap waktu dibandingkan baseline pada masing-masing kelompok
contrast(
em2,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -9.38 3.53 87 -2.659 0.0280
## M8 - M0 -5.72 4.43 87 -1.293 0.1995
## M12 - M0 -12.18 4.66 87 -2.616 0.0280
##
## kelompok = Diet:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -19.85 3.53 87 -5.625 <0.0001
## M8 - M0 -34.31 4.43 87 -7.750 <0.0001
## M12 - M0 -49.02 4.66 87 -10.524 <0.0001
##
## kelompok = Diet+Aktivitas:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -29.42 3.53 87 -8.338 <0.0001
## M8 - M0 -46.61 4.43 87 -10.528 <0.0001
## M12 - M0 -66.11 4.66 87 -14.193 <0.0001
##
## P value adjustment: holm method for 3 tests
# Perbandingan antarkelompok pada setiap waktu
em2b <- emmeans(aov2, ~ kelompok | waktu)
pairs(em2b, adjust = "tukey")
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet -6.71 8.48 87 -0.791 0.7095
## Kontrol - (Diet+Aktivitas) -10.33 8.48 87 -1.218 0.4456
## Diet - (Diet+Aktivitas) -3.62 8.48 87 -0.427 0.9043
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 3.76 7.70 87 0.488 0.8772
## Kontrol - (Diet+Aktivitas) 9.71 7.70 87 1.261 0.4212
## Diet - (Diet+Aktivitas) 5.95 7.70 87 0.773 0.7207
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 21.88 7.91 87 2.765 0.0189
## Kontrol - (Diet+Aktivitas) 30.56 7.91 87 3.861 0.0006
## Diet - (Diet+Aktivitas) 8.68 7.91 87 1.096 0.5189
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 30.13 8.18 87 3.683 0.0012
## Kontrol - (Diet+Aktivitas) 43.60 8.18 87 5.330 <0.0001
## Diet - (Diet+Aktivitas) 13.47 8.18 87 1.647 0.2318
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4d. KONTRAS INTERAKSI --------------------------------------------------------
# Membandingkan perubahan dari minggu ke-0 sampai minggu ke-12 antarkelompok
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 - Diet 36.8 6.59 87 5.592 <0.0001
## M12-M0 Kontrol - (Diet+Aktivitas) 53.9 6.59 87 8.187 <0.0001
## M12-M0 Diet - (Diet+Aktivitas) 17.1 6.59 87 2.595 0.0111
##
## P value adjustment: holm method for 3 tests
# Tren polinomial waktu di setiap kelompok
contrast(em2, "poly")
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## linear -32.89 15.00 87 -2.195 0.0308
## quadratic 2.92 4.76 87 0.614 0.5410
## cubic -23.16 11.70 87 -1.983 0.0505
##
## kelompok = Diet:
## contrast estimate SE df t.ratio p.value
## linear -161.53 15.00 87 -10.779 <0.0001
## quadratic 5.14 4.76 87 1.079 0.2835
## cubic -5.62 11.70 87 -0.481 0.6317
##
## kelompok = Diet+Aktivitas:
## contrast estimate SE df t.ratio p.value
## linear -215.53 15.00 87 -14.383 <0.0001
## quadratic 9.92 4.76 87 2.083 0.0402
## cubic -14.53 11.70 87 -1.244 0.2168
# Membandingkan tren linear antarkelompok
tren_int <- summary(
contrast(
em_full,
interaction = c(waktu = "poly", kelompok = "pairwise"),
adjust = "none"
)
)
tren_lin <- subset(tren_int, waktu_poly == "linear")
tren_lin$p.holm <- p.adjust(tren_lin$p.value, method = "holm")
tren_lin
## waktu_poly kelompok_pairwise estimate SE df t.ratio
## 1 linear Kontrol - Diet 128.63667 21.19312 87 6.069738
## 4 linear Kontrol - (Diet+Aktivitas) 182.64333 21.19312 87 8.618049
## 7 linear Diet - (Diet+Aktivitas) 54.00667 21.19312 87 2.548312
## p.value p.holm
## 1 3.256211e-08 6.512423e-08
## 4 2.717529e-13 8.152587e-13
## 7 1.257968e-02 1.257968e-02
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# LMM tidak mensyaratkan sfierisitas dan dapat menangani data hilang
# dengan lebih fleksibel dibanding repeated measure ANOVA klasik.
# Random intercept
lmm1 <- lmer(
gds ~ kelompok * waktu + (1 | id),
data = dat_long,
REML = TRUE
)
# Random intercept + random slope waktu numerik
lmm2 <- lmer(
gds ~ kelompok * waktu + (1 + minggu | id),
data = dat_long,
REML = TRUE
)
# Membandingkan struktur random effect
anova(lmm1, lmm2, refit = FALSE)
## Data: dat_long
## Models:
## lmm1: gds ~ kelompok * waktu + (1 | id)
## lmm2: gds ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 3203.1 3257.5 -1587.5 3175.1
## lmm2 16 3196.5 3258.7 -1582.2 3164.5 10.566 2 0.005077 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji efek tetap dengan derajat bebas Kenward-Roger
anova(lmm2, ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 1232 616.2 2 87.00 3.2875 0.04204 *
## waktu 48504 16167.9 3 185.37 85.8646 < 2.2e-16 ***
## kelompok:waktu 15119 2519.8 6 206.40 13.3674 8.722e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Intraclass Correlation Coefficient
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.757
## Unadjusted ICC: 0.549
# Diagnostik residual LMM
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))
# SIMULASI DATA HILANG UNTUK MENUNJUKKAN KEUNGGULAN LMM
set.seed(1)
dat_miss <- dat_long
# Membuat 30 nilai GDS pasca-baseline menjadi missing secara acak
# Kira-kira 11% dari pengukuran pasca-baseline
dat_miss$gds[
sample(which(dat_miss$waktu != "M0"), 30)
] <- NA
# LMM masih dapat memakai pasien dengan sebagian observasi tersedia
lmm_miss <- lmer(
gds ~ kelompok * waktu + (1 + minggu | id),
data = dat_miss,
control = lmerControl(optimizer = "bobyqa")
)
anova(lmm_miss, ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 1324 662.1 2 86.982 3.2981 0.04162 *
## waktu 58705 19568.4 3 168.607 97.0147 < 2.2e-16 ***
## kelompok:waktu 18936 3155.9 6 186.374 15.6282 1.681e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Jumlah pasien yang memiliki minimal satu data hilang
n_distinct(dat_miss$id[is.na(dat_miss$gds)])
## [1] 25
# 6. SIMPAN DATA & SESSION INFO
# Menyimpan data wide
write.csv(
dat_wide,
"data_diabetes_gds_wide.csv",
row.names = FALSE
)
# Menyimpan data long
write.csv(
dat_long,
"data_diabetes_gds_long.csv",
row.names = FALSE
)
# Menyimpan statistik deskriptif
write.csv(
desk,
"ringkasan_deskriptif_gds.csv",
row.names = FALSE
)
# Informasi versi R dan paket
sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_Indonesia.utf8 LC_CTYPE=English_Indonesia.utf8
## [3] LC_MONETARY=English_Indonesia.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_Indonesia.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] gtable_0.3.6 xfun_0.61 bslib_0.12.0
## [4] bayestestR_0.19.0 insight_1.5.4 lattice_0.22-9
## [7] numDeriv_2016.8-1.1 vctrs_0.7.3 tools_4.6.1
## [10] Rdpack_2.6.6 generics_0.1.4 pbkrtest_0.5.5
## [13] parallel_4.6.1 datawizard_1.4.0 tibble_3.3.1
## [16] pkgconfig_2.0.3 WRS2_1.1-7 RColorBrewer_1.1-3
## [19] S7_0.2.2 lifecycle_1.0.5 compiler_4.6.1
## [22] farver_2.1.2 stringr_1.6.0 htmltools_0.5.9
## [25] sass_0.4.10 yaml_2.3.12 Formula_1.2-6
## [28] ggpubr_1.0.0 pillar_1.11.1 nloptr_2.2.1
## [31] jquerylib_0.1.4 MASS_7.3-65 cachem_1.1.0
## [34] reformulas_0.4.4 boot_1.3-32 abind_1.4-8
## [37] nlme_3.1-169 tidyselect_1.2.1 digest_0.6.39
## [40] performance_0.18.2 mvtnorm_1.4-2 stringi_1.8.9
## [43] reshape2_1.4.5 purrr_1.2.2 labeling_0.4.3
## [46] splines_4.6.1 fastmap_1.2.0 grid_4.6.1
## [49] cli_3.6.6 magrittr_2.0.5 utf8_1.2.6
## [52] broom_1.0.13 withr_3.0.3 scales_1.4.0
## [55] backports_1.5.1 estimability_2.0.0 rmarkdown_2.32
## [58] ggsignif_0.6.4 evaluate_1.0.5 knitr_1.52
## [61] parameters_0.29.3 rbibutils_2.4.1 rlang_1.3.0
## [64] Rcpp_1.1.2 glue_1.8.1 reshape_0.8.10
## [67] minqa_1.2.8 jsonlite_2.0.0 R6_2.6.1
## [70] plyr_1.8.9
# =============================================================================
# CATATAN INTERPRETASI SINGKAT
# =============================================================================
# 1. Deskriptif:
# Perhatikan mean dan SD GDS tiap kelompok pada M0, M4, M8, dan M12.
#
# 2. Efek waktu:
# Jika p < 0.05, terdapat perubahan GDS yang signifikan sepanjang waktu.
#
# 3. Efek kelompok:
# Jika p < 0.05, secara keseluruhan terdapat perbedaan GDS antarkelompok.
#
# 4. Interaksi kelompok x waktu:
# Jika p < 0.05, pola perubahan GDS berbeda antarkelompok.
# Ini biasanya merupakan hasil paling penting dalam penelitian intervensi.
#
# 5. Post hoc:
# Digunakan untuk mengetahui waktu mana yang berbeda dan kelompok mana
# yang berbeda setelah hasil utama signifikan.
#
# 6. Effect size:
# Gunakan partial eta squared (pes), generalized eta squared (ges), atau
# omega squared untuk menjelaskan besar efek, bukan hanya signifikansi.
#
# 7. LMM:
# Digunakan sebagai analisis pembanding, terutama bila terdapat data hilang
# atau struktur korelasi antarwaktu yang lebih kompleks.
# =============================================================================