# Nama : P. Rupian Nur Amin
# NIM : 2611018019
#=============================================================================
#REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: Program edukasi diet rendah purin dan kadar asam urat
# pada pasien dengan kadar asam urat tinggi 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: simulasi uji klinis di puskesmas. 90 pasien dengan asam urat tinggi
# diacak ke tiga kelompok (n = 30 per kelompok):
# - Kontrol : edukasi standar puskesmas
# - Diet_RP : edukasi diet rendah purin
# - Diet_RP+Konseling : edukasi diet rendah purin + konseling gizi berkala
# Kadar asam urat serum (AU, 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", "Diet_RP", "Diet_RP+Konseling")
minggu <- c(0, 4, 8, 12)
# Rerata SIMULASI (mg/dL) per kelompok x waktu; bukan bukti efek klinis
mu <- rbind(
"Kontrol" = c(8.50, 8.45, 8.40, 8.375),
"Diet_RP" = c(8.50, 8.275, 8.125, 8.025),
"Diet_RP+Konseling" = c(8.50, 8.175, 7.925, 7.750)
)
sd_int <- 0.45 # SD intersep acak (perbedaan asam urat dasar antarpasien)
sd_slope <- 0.0225 # SD kemiringan acak per minggu (perbedaan respons antarpasien)
sd_eps <- 0.225 # 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("AU_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(.45, .55)),
round(y, 2))
}) |> 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("AU_M"), names_to = "waktu", values_to = "asam_urat") |>
mutate(waktu = factor(waktu, levels = paste0("AU_M", minggu),
labels = paste0("M", minggu)),
minggu = as.numeric(sub("M", "", waktu)))
head(dat_wide)
## id kelompok usia jk AU_M0 AU_M4 AU_M8 AU_M12
## 1 P001 Kontrol 42 L 8.70 8.20 8.97 8.47
## 2 P002 Kontrol 43 L 8.05 8.31 8.42 8.01
## 3 P003 Kontrol 62 P 8.52 8.38 8.47 8.59
## 4 P004 Kontrol 40 P 8.40 8.66 8.62 8.54
## 5 P005 Kontrol 53 P 7.81 7.76 7.74 8.25
## 6 P006 Kontrol 42 L 7.31 7.36 7.58 7.81
head(dat_long)
## # A tibble: 6 × 7
## id kelompok usia jk waktu asam_urat minggu
## <fct> <fct> <dbl> <chr> <fct> <dbl> <dbl>
## 1 P001 Kontrol 42 L M0 8.7 0
## 2 P001 Kontrol 42 L M4 8.2 4
## 3 P001 Kontrol 42 L M8 8.97 8
## 4 P001 Kontrol 42 L M12 8.47 12
## 5 P002 Kontrol 43 L M0 8.05 0
## 6 P002 Kontrol 43 L M4 8.31 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","Diet_RP",..: 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] "L" "L" "L" "L" ...
## $ waktu : Factor w/ 4 levels "M0","M4","M8",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ asam_urat: num [1:360] 8.7 8.2 8.97 8.47 8.05 8.31 8.42 8.01 8.52 8.38 ...
## $ 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(asam_urat, type = "mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 asam_urat 30 8.48 0.594
## 2 Kontrol M4 asam_urat 30 8.35 0.54
## 3 Kontrol M8 asam_urat 30 8.41 0.476
## 4 Kontrol M12 asam_urat 30 8.31 0.536
## 5 Diet_RP M0 asam_urat 30 8.57 0.472
## 6 Diet_RP M4 asam_urat 30 8.35 0.525
## 7 Diet_RP M8 asam_urat 30 8.18 0.485
## 8 Diet_RP M12 asam_urat 30 8.08 0.496
## 9 Diet_RP+Konseling M0 asam_urat 30 8.53 0.588
## 10 Diet_RP+Konseling M4 asam_urat 30 8.21 0.595
## 11 Diet_RP+Konseling M8 asam_urat 30 7.99 0.642
## 12 Diet_RP+Konseling M12 asam_urat 30 7.87 0.722
# Matriks kovarians & korelasi antarwaktu (seluruh subjek, dalam kelompok)
# -> memberi gambaran awal apakah sfierisitas masuk akal
S <- cov(dat_wide[, paste0("AU_M", minggu)])
R <- cor(dat_wide[, paste0("AU_M", minggu)])
round(S, 1); round(R, 2)
## AU_M0 AU_M4 AU_M8 AU_M12
## AU_M0 0.3 0.2 0.2 0.2
## AU_M4 0.2 0.3 0.3 0.3
## AU_M8 0.2 0.3 0.3 0.3
## AU_M12 0.2 0.3 0.3 0.4
## AU_M0 AU_M4 AU_M8 AU_M12
## AU_M0 1.00 0.79 0.72 0.64
## AU_M4 0.79 1.00 0.81 0.77
## AU_M8 0.72 0.81 1.00 0.85
## AU_M12 0.64 0.77 0.85 1.00
# Varians selisih antarpasangan waktu (inti asumsi sfierisitas)
pasangan <- combn(paste0("AU_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)
## AU_M0 - AU_M4 AU_M0 - AU_M8 AU_M0 - AU_M12 AU_M4 - AU_M8 AU_M4 - AU_M12
## 0.1 0.2 0.2 0.1 0.2
## AU_M8 - AU_M12
## 0.1
# Profile plot: rerata +/- 95% CI per kelompok
p_profil <- ggplot(dat_long, aes(minggu, asam_urat, 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 = "Kadar asam urat serum (mg/dL)", colour = "Kelompok",
title = "Profil rerata asam urat (± 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, asam_urat, 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 = "asam urat (mg/dL)", title = "Lintasan individu dan rerata kelompok")
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan: apakah asam urat berubah selama 12 minggu pada kelompok Diet_RP+Konseling?
d1 <- droplevels(filter(dat_long, kelompok == "Diet_RP+Konseling"))
d1w <- filter(dat_wide, kelompok == "Diet_RP+Konseling")
## 3a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per waktu (ekstrem = di luar Q1-3IQR / Q3+3IQR)
d1 |> group_by(waktu) |> identify_outliers(asam_urat)
## [1] waktu id kelompok usia jk asam_urat 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(asam_urat)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 asam_urat 0.959 0.298
## 2 M4 asam_urat 0.953 0.200
## 3 M8 asam_urat 0.982 0.868
## 4 M12 asam_urat 0.960 0.301
ggpubr::ggqqplot(d1, "asam_urat", facet.by = "waktu")
# (iii) Sfierisitas: Mauchly (dilaporkan otomatis oleh anova_test & afex)
aov1_rs <- anova_test(data = d1, dv = asam_urat, 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 29.791 2.41e-13 * 0.507
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.535 0.004 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05
## 1 waktu 0.709 2.13, 61.7 4.2e-10 * 0.767 2.3, 66.73 9.52e-11 *
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.13 61.7 29.791 4.2e-10 * 0.507
## 3b. ANOVA dengan afex (sumber utama laporan) -------------------------------
aov1 <- aov_ez(id = "id", dv = "asam_urat", 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: asam_urat
## Effect df MSE F ges pes p.value
## 1 waktu 2.13, 61.70 0.12 29.79 *** .137 .507 <.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) 7968.9 1 40.079 29 5766.138 < 2.2e-16 ***
## waktu 7.5 3 7.323 87 29.791 2.409e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.53539 0.0039558
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.70922 4.203e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.76696 9.515012e-11
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.51 | [0.38, 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.13 | [0.02, 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.9950 5766.1 1 29 < 2.2e-16 ***
## waktu 1 0.6635 17.7 3 27 1.452e-06 ***
## ---
## 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 8.53 0.107 29 8.31 8.75
## M4 8.21 0.109 29 7.99 8.44
## M8 7.99 0.117 29 7.75 8.23
## M12 7.87 0.132 29 7.60 8.14
##
## 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.314 0.0721 29 4.350 0.0009
## M0 - M8 0.541 0.0756 29 7.154 <0.0001
## M0 - M12 0.656 0.1020 29 6.430 <0.0001
## M4 - M8 0.227 0.0502 29 4.526 0.0006
## M4 - M12 0.342 0.0761 29 4.496 0.0006
## M8 - M12 0.115 0.0636 29 1.815 0.4796
##
## 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.314 0.0721 29 -4.350 0.0002
## M8 - M0 -0.541 0.0756 29 -7.154 <0.0001
## M12 - M0 -0.656 0.1020 29 -6.430 <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 -2.195 0.3210 29 -6.839 <0.0001
## quadratic 0.198 0.0885 29 2.242 0.0328
## cubic 0.025 0.1620 29 0.154 0.8784
## 3e. Alternatif nonparametrik ------------------------------------------------
friedman_test(d1, asam_urat ~ waktu | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 asam_urat 30 44.6 3 0.00000000112 Friedman test
friedman_effsize(d1, asam_urat ~ waktu | id) # Kendall's W
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 asam_urat 30 0.496 Kendall W moderate
d1 |> wilcox_test(asam_urat ~ 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 asam_urat M0 M4 30 30 406. 0.000141 8.45e-4 ***
## 2 asam_urat M0 M8 30 30 445 0.000000691 4.15e-6 ****
## 3 asam_urat M0 M12 30 30 440 0.00000168 1.01e-5 ****
## 4 asam_urat M4 M8 30 30 401 0.000239 1.43e-3 **
## 5 asam_urat M4 M12 30 30 410 0.0000961 5.77e-4 ***
## 6 asam_urat M8 M12 30 30 322 0.0658 3.95e-1 ns
# (Opsional) ANOVA robust berbasis trimmed mean -- paket WRS2
if (requireNamespace("WRS2", quietly = TRUE)) {
print(WRS2::rmanova(d1$asam_urat, d1$waktu, d1$id, tr = 0.2))
}
## Call:
## WRS2::rmanova(y = d1$asam_urat, groups = d1$waktu, blocks = d1$id,
## tr = 0.2)
##
## Test statistic: F = 22.2234
## Degrees of freedom 1: 2.22
## Degrees of freedom 2: 37.74
## p-value: 0
# 4. MIXED DESIGN ANOVA (Kelompok [between] x Waktu [within])
# Pertanyaan: apakah pola perubahan asam urat berbeda antarkelompok intervensi?
## 4a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per sel
dat_long |> group_by(kelompok, waktu) |> identify_outliers(asam_urat)
## # A tibble: 2 × 9
## kelompok waktu id usia jk asam_urat minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol M8 P015 57 P 7.2 8 TRUE FALSE
## 2 Kontrol M12 P015 57 P 6.91 12 TRUE FALSE
# (ii) Normalitas per sel (3 x 4 = 12 sel) dan residual model
dat_long |> group_by(kelompok, waktu) |> shapiro_test(asam_urat)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 asam_urat 0.990 0.989
## 2 Kontrol M4 asam_urat 0.963 0.376
## 3 Kontrol M8 asam_urat 0.972 0.587
## 4 Kontrol M12 asam_urat 0.976 0.717
## 5 Diet_RP M0 asam_urat 0.961 0.322
## 6 Diet_RP M4 asam_urat 0.980 0.831
## 7 Diet_RP M8 asam_urat 0.968 0.486
## 8 Diet_RP M12 asam_urat 0.974 0.660
## 9 Diet_RP+Konseling M0 asam_urat 0.959 0.298
## 10 Diet_RP+Konseling M4 asam_urat 0.953 0.200
## 11 Diet_RP+Konseling M8 asam_urat 0.982 0.868
## 12 Diet_RP+Konseling M12 asam_urat 0.960 0.301
ggpubr::ggqqplot(dat_long, "asam_urat", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antarkelompok pada TIAP waktu (Levene, median-centered)
dat_long |> group_by(waktu) |> levene_test(asam_urat ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 0.852 0.430
## 2 M4 2 87 0.682 0.508
## 3 M8 2 87 2.03 0.138
## 4 M12 2 87 4.73 0.0112
# (iv) Homogenitas matriks kovarians antarkelompok (Box's M; uji pada alpha = .001)
box_m(dat_wide[, paste0("AU_M", minggu)], dat_wide$kelompok)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 19.5 0.487 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 = "asam_urat", data = dat_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: asam_urat
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 1.05 1.64 .030 .036 .201
## 2 waktu 2.61, 226.78 0.08 45.61 *** .080 .344 <.001
## 3 kelompok:waktu 5.21, 226.78 0.08 6.38 *** .024 .128 <.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) 24661.0 1 91.325 87 23493.1790 < 2.2e-16 ***
## kelompok 3.4 2 91.325 87 1.6369 0.2005
## waktu 9.5 3 18.058 261 45.6057 < 2.2e-16 ***
## kelompok:waktu 2.6 6 18.058 261 6.3756 2.798e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.80987 0.0028537
## kelompok:waktu 0.80987 0.0028537
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.86889 < 2.2e-16 ***
## kelompok:waktu 0.86889 1.052e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8981199 1.435661e-21
## kelompok:waktu 0.8981199 7.826445e-06
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.99631 23493.2 1 87 < 2.2e-16 ***
## kelompok 2 0.03627 1.6 2 87 0.2005195
## waktu 1 0.51456 30.0 3 85 2.463e-13 ***
## kelompok:waktu 2 0.25355 4.2 6 172 0.0006228 ***
## ---
## 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 = asam_urat, 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 1.637 2.01e-01 0.036
## 2 waktu 2.61 226.78 45.606 5.99e-21 * 0.344
## 3 kelompok:waktu 5.21 226.78 6.376 1.05e-05 * 0.128
# Ukuran efek
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.04 | [0.00, 1.00]
## waktu | 0.34 | [0.27, 1.00]
## kelompok:waktu | 0.13 | [0.06, 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.01 | [0.00, 1.00]
## waktu | 0.08 | [0.03, 1.00]
## kelompok:waktu | 0.02 | [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 = "asam urat (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 (bila 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.935 0.1298
##
## kelompok = Diet_RP:
## model term df1 df2 F.ratio p.value
## waktu 3 87 13.607 <0.0001
##
## kelompok = Diet_RP+Konseling:
## model term df1 df2 F.ratio p.value
## waktu 3 87 24.986 <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.224 0.8001
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.609 0.5461
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 4.612 0.0125
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 4.106 0.0198
# 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 -0.123 0.0646 87 -1.900 0.1461
## M8 - M0 -0.067 0.0677 87 -0.989 0.3252
## M12 - M0 -0.166 0.0830 87 -1.999 0.1461
##
## kelompok = Diet_RP:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -0.224 0.0646 87 -3.475 0.0008
## M8 - M0 -0.395 0.0677 87 -5.833 <0.0001
## M12 - M0 -0.487 0.0830 87 -5.865 <0.0001
##
## kelompok = Diet_RP+Konseling:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -0.314 0.0646 87 -4.858 <0.0001
## M8 - M0 -0.541 0.0677 87 -7.984 <0.0001
## M12 - M0 -0.656 0.0830 87 -7.900 <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 - Diet_RP -0.0957 0.143 87 -0.668 0.7825
## Kontrol - (Diet_RP+Konseling) -0.0513 0.143 87 -0.359 0.9317
## Diet_RP - (Diet_RP+Konseling) 0.0443 0.143 87 0.310 0.9486
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet_RP 0.0060 0.143 87 0.042 0.9990
## Kontrol - (Diet_RP+Konseling) 0.1397 0.143 87 0.976 0.5939
## Diet_RP - (Diet_RP+Konseling) 0.1337 0.143 87 0.934 0.6203
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet_RP 0.2323 0.139 87 1.668 0.2233
## Kontrol - (Diet_RP+Konseling) 0.4223 0.139 87 3.032 0.0089
## Diet_RP - (Diet_RP+Konseling) 0.1900 0.139 87 1.364 0.3642
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet_RP 0.2253 0.153 87 1.472 0.3094
## Kontrol - (Diet_RP+Konseling) 0.4387 0.153 87 2.865 0.0143
## Diet_RP - (Diet_RP+Konseling) 0.2133 0.153 87 1.393 0.3488
##
## 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 - Diet_RP 0.321 0.117 87 2.734 0.0152
## M12-M0 Kontrol - (Diet_RP+Konseling) 0.490 0.117 87 4.173 0.0002
## M12-M0 Diet_RP - (Diet_RP+Konseling) 0.169 0.117 87 1.439 0.1537
##
## 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.442 0.266 87 -1.662 0.1001
## linear Diet_RP -1.632 0.266 87 -6.131 <0.0001
## linear Diet_RP+Konseling -2.195 0.266 87 -8.248 <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
## 1 linear Kontrol - Diet_RP 1.1893333 0.3763644 87 3.160058
## 4 linear Kontrol - (Diet_RP+Konseling) 1.7526667 0.3763644 87 4.656834
## 7 linear Diet_RP - (Diet_RP+Konseling) 0.5633333 0.3763644 87 1.496776
## p.value p.holm
## 1 2.170353e-03 0.0043407057
## 4 1.144627e-05 0.0000343388
## 7 1.380707e-01 0.1380707024
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Tidak mensyaratkan sfierisitas, menampung data hilang (MAR) dan waktu
# pengukuran yang tidak seragam.
lmm1 <- lmer(asam_urat ~ kelompok * waktu + (1 | id), data = dat_long, REML = TRUE)
lmm2 <- lmer(asam_urat ~ 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: asam_urat ~ kelompok * waktu + (1 | id)
## lmm2: asam_urat ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 380.61 435.02 -176.31 352.61
## lmm2 16 367.07 429.24 -167.53 335.07 17.544 2 0.000155 ***
## ---
## 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 0.1659 0.08293 2 87.00 1.6369 0.2005208
## waktu 4.6994 1.56648 3 185.37 30.7776 3.413e-16 ***
## kelompok:waktu 1.3903 0.23171 6 206.40 4.5475 0.0002353 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1) # korelasi intrakelas
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.780
## Unadjusted ICC: 0.685
# 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