# =============================================================================
# Nama : Rofiani Vitrilia Rachman
# NIM : 2611018018
# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: Program intervensi gaya hidup dan kadar glukosa darah puasa
# 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 & 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 diabetes melitus tipe 2
# diacak ke tiga kelompok (n = 30 per kelompok):
# - Kontrol : edukasi standar puskesmas
# - DM : edukasi diet diabetes dan pengaturan pola makan
# - DM+AF : edukasi diet diabetes + aktivitas fisik terstruktur (senam 3x/minggu)
# Glukosa darah puasa (GDP, 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", "DM", "DM+AF")
minggu <- c(0, 4, 8, 12)
# Rerata populasi (mmHg) per kelompok x waktu
mu <- rbind(
"Kontrol" = c(185, 183, 181, 180),
"DM" = c(185, 174, 166, 158),
"DM+AF" = c(185, 170, 157, 146)
)
# SD intersep acak (perbedaan GDP dasar antar pasien)
sd_int <- 18
sd_slope <- 0.55 # SD kemiringan acak per minggu (perbedaan respons antarpasien)
sd_eps <- 8 # 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("GDP_M", minggu)
data.frame(id = sprintf("P%03d", id),
kelompok = kel_lab[g],
usia = round(runif(n_per, 35, 70)),
jk = sample(c("L", "P"), n_per, replace = TRUE, prob = c(.45, .55)),
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("GDP_M"), names_to = "waktu", values_to = "gdp") |>
mutate(waktu = factor(waktu, levels = paste0("GDP_M", minggu),
labels = paste0("M", minggu)),
minggu = as.numeric(sub("M", "", waktu)))
head(dat_wide)
## id kelompok usia jk GDP_M0 GDP_M4 GDP_M8 GDP_M12
## 1 P001 Kontrol 43 L 193.2 175.5 203.3 186.1
## 2 P002 Kontrol 44 L 166.9 175.7 179.4 164.3
## 3 P003 Kontrol 67 P 185.9 181.1 184.2 188.7
## 4 P004 Kontrol 41 P 181.1 190.1 188.1 184.7
## 5 P005 Kontrol 56 P 159.1 156.8 155.9 173.8
## 6 P006 Kontrol 43 L 137.5 138.2 144.6 151.5
head(dat_long)
## # A tibble: 6 × 7
## id kelompok usia jk waktu gdp minggu
## <fct> <fct> <dbl> <chr> <fct> <dbl> <dbl>
## 1 P001 Kontrol 43 L M0 193. 0
## 2 P001 Kontrol 43 L M4 176. 4
## 3 P001 Kontrol 43 L M8 203. 8
## 4 P001 Kontrol 43 L M12 186. 12
## 5 P002 Kontrol 44 L M0 167. 0
## 6 P002 Kontrol 44 L M4 176. 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","DM",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num [1:360] 43 43 43 43 44 44 44 44 67 67 ...
## $ jk : chr [1:360] "L" "L" "L" "L" ...
## $ waktu : Factor w/ 4 levels "M0","M4","M8",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ gdp : num [1:360] 193 176 203 186 167 ...
## $ 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(gdp, type = "mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 gdp 30 184. 23.1
## 2 Kontrol M4 gdp 30 180. 21.1
## 3 Kontrol M8 gdp 30 181. 18.5
## 4 Kontrol M12 gdp 30 178. 20.2
## 5 DM M0 gdp 30 188. 18.8
## 6 DM M4 gdp 30 177. 20.8
## 7 DM M8 gdp 30 168. 19.4
## 8 DM M12 gdp 30 160. 19.3
## 9 DM+AF M0 gdp 30 186. 23.1
## 10 DM+AF M4 gdp 30 171. 22.9
## 11 DM+AF M8 gdp 30 159. 24.2
## 12 DM+AF M12 gdp 30 149. 25.6
# Matriks kovarians & korelasi antarwaktu (seluruh subjek, dalam kelompok)
# -> memberi gambaran awal apakah sfierisitas masuk akal
S <- cov(dat_wide[, paste0("GDP_M", minggu)])
R <- cor(dat_wide[, paste0("GDP_M", minggu)])
round(S, 1); round(R, 2)
## GDP_M0 GDP_M4 GDP_M8 GDP_M12
## GDP_M0 465.8 385.3 362.0 356.7
## GDP_M4 385.3 470.5 411.7 423.3
## GDP_M8 362.0 411.7 513.9 496.1
## GDP_M12 356.7 423.3 496.1 610.5
## GDP_M0 GDP_M4 GDP_M8 GDP_M12
## GDP_M0 1.00 0.82 0.74 0.67
## GDP_M4 0.82 1.00 0.84 0.79
## GDP_M8 0.74 0.84 1.00 0.89
## GDP_M12 0.67 0.79 0.89 1.00
# Varians selisih antarpasangan waktu (inti asumsi sfierisitas)
pasangan <- combn(paste0("GDP_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)
## GDP_M0 - GDP_M4 GDP_M0 - GDP_M8 GDP_M0 - GDP_M12 GDP_M4 - GDP_M8
## 165.7 255.6 363.1 161.0
## GDP_M4 - GDP_M12 GDP_M8 - GDP_M12
## 234.5 132.2
# Profile plot: rerata +/- 95% CI per kelompok
p_profil <- ggplot(dat_long, aes(minggu, gdp, 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 = "Glukosa darah puasa (mg/dL)", colour = "Kelompok",
title = "Profil rerata GDP (± 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 tiap pasien
p_spag <- ggplot(dat_long, aes(minggu, gdp, 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 = "GDP (mg/dL)", title = "Lintasan individu dan rerata kelompok")
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan: apakah GDP berubah selama 12 minggu pada kelompok DM+AF?
d1 <- droplevels(filter(dat_long, kelompok == "DM+AF"))
d1w <- filter(dat_wide, kelompok == "DM+AF")
## 3a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per waktu (ekstrem = di luar Q1-3IQR / Q3+3IQR)
d1 |> group_by(waktu) |> identify_outliers(gdp)
## [1] waktu id kelompok usia jk gdp 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(gdp)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 gdp 0.960 0.303
## 2 M4 gdp 0.957 0.262
## 3 M8 gdp 0.974 0.640
## 4 M12 gdp 0.969 0.505
ggpubr::ggqqplot(d1, "gdp", facet.by = "waktu")
# (iii) Sfierisitas: Mauchly (dilaporkan otomatis oleh anova_test & afex)
aov1_rs <- anova_test(data = d1, dv = gdp, 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 86.185 5.72e-26 * 0.748
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.695 0.073
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.812 2.44, 70.68 1.53e-21 * 0.892 2.68, 77.64 1.97e-23
## p[HF]<.05
## 1 *
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 3 87 86.185 5.72e-26 * 0.748
## 3b. ANOVA dengan afex (sumber utama laporan) -------------------------------
aov1 <- aov_ez(id = "id", dv = "gdp", 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: gdp
## Effect df MSE F ges pes p.value
## 1 waktu 2.44, 70.68 107.73 86.18 *** .253 .748 <.001
## ---
## 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) 3313729 1 59141 29 1624.893 < 2.2e-16 ***
## waktu 22629 3 7614 87 86.185 < 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.69491 0.072886
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.8124 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8924254 1.973737e-23
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.75 | [0.67, 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.24 | [0.11, 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.98247 1624.89 1 29 < 2.2e-16 ***
## waktu 1 0.86074 55.63 3 27 1.099e-11 ***
## ---
## 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 186 4.22 29 177 195
## M4 171 4.19 29 162 180
## M8 159 4.41 29 150 168
## M12 149 4.68 29 140 159
##
## Confidence level used: 0.95
pairs(em1, adjust = "bonferroni") # semua pasangan waktu (6 perbandingan)
## contrast estimate SE df t.ratio p.value
## M0 - M4 14.89 2.44 29 6.109 <0.0001
## M0 - M8 27.36 2.39 29 11.423 <0.0001
## M0 - M12 36.57 3.08 29 11.863 <0.0001
## M4 - M8 12.46 1.75 29 7.126 <0.0001
## M4 - M12 21.67 2.47 29 8.771 <0.0001
## M8 - M12 9.21 2.16 29 4.265 0.0012
##
## 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 -14.9 2.44 29 -6.109 <0.0001
## M8 - M0 -27.4 2.39 29 -11.423 <0.0001
## M12 - M0 -36.6 3.08 29 -11.863 <0.0001
##
## P value adjustment: holm method for 3 tests
contrast(em1, "poly") # tren linear, kuadratik, kubik
## contrast estimate SE df t.ratio p.value
## linear -122.163 9.61 29 -12.718 <0.0001
## quadratic 5.683 3.14 29 1.807 0.0811
## cubic 0.823 5.77 29 0.143 0.8876
## 3e. Alternatif nonparametrik ------------------------------------------------
friedman_test(d1, gdp ~ waktu | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 gdp 30 67.9 3 1.21e-14 Friedman test
friedman_effsize(d1, gdp ~ waktu | id) # Kendall's W
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 gdp 30 0.754 Kendall W large
d1 |> wilcox_test(gdp ~ 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 gdp M0 M4 30 30 436. 0.00000281 1.69e-5 ****
## 2 gdp M0 M8 30 30 464 0.00000000373 2.24e-8 ****
## 3 gdp M0 M12 30 30 464 0.00000000373 2.24e-8 ****
## 4 gdp M4 M8 30 30 452 0.000000164 9.83e-7 ****
## 5 gdp M4 M12 30 30 459 0.0000000261 1.56e-7 ****
## 6 gdp M8 M12 30 30 402 0.000226 1.36e-3 **
# (Opsional) ANOVA robust berbasis trimmed mean -- paket WRS2
if (requireNamespace("WRS2", quietly = TRUE)) {
print(WRS2::rmanova(d1$gdp, d1$waktu, d1$id, tr = 0.2))
}
## Call:
## WRS2::rmanova(y = d1$gdp, groups = d1$waktu, blocks = d1$id,
## tr = 0.2)
##
## Test statistic: F = 60.0673
## Degrees of freedom 1: 2.63
## Degrees of freedom 2: 44.77
## p-value: 0
# 4. MIXED DESIGN ANOVA (Kelompok [between] x Waktu [within])
# Pertanyaan: apakah pola perubahan GDP berbeda antarkelompok intervensi?
## 4a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per sel
dat_long |> group_by(kelompok, waktu) |> identify_outliers(gdp)
## # A tibble: 2 × 9
## kelompok waktu id usia jk gdp minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol M8 P015 61 P 135 8 TRUE FALSE
## 2 Kontrol M12 P015 61 P 125. 12 TRUE FALSE
# (ii) Normalitas per sel (3 x 4 = 12 sel) dan residual model
dat_long |> group_by(kelompok, waktu) |> shapiro_test(gdp)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 gdp 0.991 0.996
## 2 Kontrol M4 gdp 0.961 0.334
## 3 Kontrol M8 gdp 0.965 0.411
## 4 Kontrol M12 gdp 0.973 0.636
## 5 DM M0 gdp 0.962 0.348
## 6 DM M4 gdp 0.980 0.830
## 7 DM M8 gdp 0.974 0.643
## 8 DM M12 gdp 0.987 0.970
## 9 DM+AF M0 gdp 0.960 0.303
## 10 DM+AF M4 gdp 0.957 0.262
## 11 DM+AF M8 gdp 0.974 0.640
## 12 DM+AF M12 gdp 0.969 0.505
ggpubr::ggqqplot(dat_long, "gdp", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antarkelompok pada TIAP waktu (Levene, median-centered)
dat_long |> group_by(waktu) |> levene_test(gdp ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 0.709 0.495
## 2 M4 2 87 0.551 0.578
## 3 M8 2 87 1.56 0.216
## 4 M12 2 87 2.68 0.0746
# (iv) Homogenitas matriks kovarians antarkelompok (Box's M; uji pada alpha = .001)
box_m(dat_wide[, paste0("GDP_M", minggu)], dat_wide$kelompok)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 17.4 0.626 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 = "gdp", data = dat_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: gdp
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 1626.36 3.91 * .073 .082 .024
## 2 waktu 2.83, 245.97 81.33 116.96 *** .143 .573 <.001
## 3 kelompok:waktu 5.65, 245.97 81.33 19.98 *** .054 .315 <.001
## ---
## 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) 10811494 1 141494 87 6647.6487 < 2e-16 ***
## kelompok 12716 2 141494 87 3.9093 0.02367 *
## waktu 26894 3 20005 261 116.9599 < 2e-16 ***
## kelompok:waktu 9190 6 20005 261 19.9833 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.91481 0.17768
## kelompok:waktu 0.91481 0.17768
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.9424 < 2.2e-16 ***
## kelompok:waktu 0.9424 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.9773238 5.375425e-47
## kelompok:waktu 0.9773238 8.117757e-19
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.98708 6647.6 1 87 < 2.2e-16 ***
## kelompok 2 0.08246 3.9 2 87 0.02367 *
## waktu 1 0.75327 86.5 3 85 < 2.2e-16 ***
## kelompok:waktu 2 0.52070 10.1 6 172 1.509e-09 ***
## ---
## 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 = gdp, 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.909 2.40e-02 * 0.082
## 2 waktu 2.83 245.97 116.960 2.05e-45 * 0.573
## 3 kelompok:waktu 5.65 245.97 19.983 3.16e-18 * 0.315
# Ukuran efek
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.08 | [0.01, 1.00]
## waktu | 0.57 | [0.51, 1.00]
## kelompok:waktu | 0.31 | [0.23, 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.06 | [0.00, 1.00]
## waktu | 0.14 | [0.08, 1.00]
## kelompok:waktu | 0.05 | [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 = "GDP (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 2.280 0.0849
##
## kelompok = DM:
## model term df1 df2 F.ratio p.value
## waktu 3 87 42.425 <0.0001
##
## kelompok = DM+AF:
## model term df1 df2 F.ratio p.value
## waktu 3 87 75.216 <0.0001
# Efek KELOMPOK pada tiap waktu
joint_tests(aov2, by = "waktu")
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.222 0.8016
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 1.207 0.3040
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 9.189 0.0002
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 13.117 <0.0001
# 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 -4.44 2.24 87 -1.983 0.1011
## M8 - M0 -2.53 2.23 87 -1.132 0.2607
## M12 - M0 -6.01 2.59 87 -2.323 0.0676
##
## kelompok = DM:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -11.11 2.24 87 -4.967 <0.0001
## M8 - M0 -19.95 2.23 87 -8.938 <0.0001
## M12 - M0 -27.81 2.59 87 -10.755 <0.0001
##
## kelompok = DM+AF:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -14.89 2.24 87 -6.656 <0.0001
## M8 - M0 -27.36 2.23 87 -12.257 <0.0001
## M12 - M0 -36.57 2.59 87 -14.141 <0.0001
##
## 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 - DM -3.74 5.62 87 -0.666 0.7839
## Kontrol - (DM+AF) -1.91 5.62 87 -0.340 0.9382
## DM - (DM+AF) 1.83 5.62 87 0.325 0.9433
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - DM 2.93 5.59 87 0.525 0.8593
## Kontrol - (DM+AF) 8.54 5.59 87 1.529 0.2824
## DM - (DM+AF) 5.61 5.59 87 1.004 0.5763
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - DM 13.68 5.38 87 2.543 0.0338
## Kontrol - (DM+AF) 22.92 5.38 87 4.260 0.0002
## DM - (DM+AF) 9.24 5.38 87 1.717 0.2046
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - DM 18.06 5.66 87 3.193 0.0055
## Kontrol - (DM+AF) 28.65 5.66 87 5.065 <0.0001
## DM - (DM+AF) 10.59 5.66 87 1.872 0.1530
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## 4d. Kontras interaksi -------------------------------------------------------
# Apakah perubahan (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 - DM 21.80 3.66 87 5.962 <0.0001
## M12-M0 Kontrol - (DM+AF) 30.56 3.66 87 8.357 <0.0001
## M12-M0 DM - (DM+AF) 8.76 3.66 87 2.394 0.0188
##
## 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 -16.1 8.24 87 -1.956 0.0537
## linear DM -92.3 8.24 87 -11.203 <0.0001
## linear DM+AF -122.2 8.24 87 -14.833 <0.0001
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 - DM 76.15667 11.64755 87 6.538426 4.086522e-09
## 4 linear Kontrol - (DM+AF) 106.05333 11.64755 87 9.105203 2.734732e-14
## 7 linear DM - (DM+AF) 29.89667 11.64755 87 2.566777 1.197458e-02
## p.holm
## 1 8.173045e-09
## 4 8.204197e-14
## 7 1.197458e-02
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Tidak mensyaratkan sfierisitas, menampung data hilang (MAR) dan waktu
# pengukuran yang tidak seragam.
lmm1 <- lmer(gdp ~ kelompok * waktu + (1 | id), data = dat_long, REML = TRUE)
lmm2 <- lmer(gdp ~ 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: gdp ~ kelompok * waktu + (1 | id)
## lmm2: gdp ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 2849.3 2903.7 -1410.7 2821.3
## lmm2 16 2846.8 2909.0 -1407.4 2814.8 6.47 2 0.03936 *
## ---
## 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 501.1 250.6 2 87.00 3.9093 0.02367 *
## waktu 17060.6 5686.9 3 185.37 88.3172 < 2.2e-16 ***
## kelompok:waktu 5870.6 978.4 6 206.40 15.1783 2.23e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1) # korelasi intrakelas
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.835
## Unadjusted ICC: 0.646
# 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$gdp[sample(which(dat_miss$waktu != "M0"), 30)] <- NA
lmm_miss <- lmer(gdp ~ 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 491.0 245.5 2 86.994 3.7196 0.02818 *
## waktu 16814.3 5604.8 3 168.197 84.5202 < 2.2e-16 ***
## kelompok:waktu 6084.4 1014.1 6 185.876 15.2744 3.369e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# RM ANOVA akan membuang seluruh pasien yang punya >= 1 nilai hilang:
n_distinct(dat_miss$id[is.na(dat_miss$gdp)])
## [1] 25
# 6. SIMPAN DATA & SESSION INFO
write.csv(dat_wide, "data_gdp_dm_wide.csv", row.names = FALSE)
write.csv(dat_long, "data_gdp_dm_long.csv", row.names = FALSE)
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_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/Singapore
## 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.60 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] datawizard_1.4.0 parallel_4.6.1 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.31
## [58] otel_0.2.0 ggsignif_0.6.4 evaluate_1.0.5
## [61] knitr_1.51 parameters_0.29.3 rbibutils_2.4.1
## [64] rlang_1.3.0 Rcpp_1.1.2 glue_1.8.1
## [67] reshape_0.8.10 rstudioapi_0.19.0 minqa_1.2.8
## [70] jsonlite_2.0.0 R6_2.6.1 plyr_1.8.9