# =============================================================================
# Nama : Chiesa Safira Yoviana Hefny
# NIM : 2611018002
# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: Intervensi diet dan olahraga terhadap gula darah puasa (GDP)
# pada pasien diabetes (DATA SIMULASI, dibaca dari file Excel)
#
# Isi:
# 0. Paket & pengaturan
# 1. Impor data Excel (format lebar -> panjang)
# 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", "readxl", "afex", "emmeans", "rstatix", "car",
# "effectsize", "lme4", "lmerTest", "performance", "ggpubr"))
suppressPackageStartupMessages({
library(readxl) # membaca file Excel
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. IMPOR DATA EXCEL
# Skenario: 90 pasien diabetes diacak ke tiga kelompok (n = 30 per kelompok):
# - Kontrol : edukasi standar
# - Diet : diet terkontrol
# - Diet+Olahraga : diet terkontrol + olahraga terstruktur
# Gula darah puasa (GDP, mg/dL) diukur 4 kali: GDP_W0 (baseline) s.d. GDP_W3.
#
# Pastikan file "data_gdp_diabetes.xlsx" berada di working directory, atau
# ganti 'file_data' dengan path lengkap, mis. "C:/Users/nama/Documents/data_gdp_diabetes.xlsx".
# Cek working directory dengan getwd(); ubah dengan setwd("folder_anda").
file_data <- "data_gdp_diabetes.xlsx"
sheet_dat <- "Data"
kel_lab <- c("Kontrol", "Diet", "Diet+Olahraga")
minggu <- c(0, 4, 8, 12) # jarak waktu pengukuran W0-W3 (minggu); ubah sesuai desain
kol_gdp <- paste0("GDP_W", 0:3) # nama kolom GDP di Excel
dat_wide <- read_excel(file_data, sheet = sheet_dat) |>
as.data.frame()
# Cek struktur & data hilang sebelum lanjut
str(dat_wide)
## 'data.frame': 90 obs. of 8 variables:
## $ id : chr "P001" "P002" "P003" "P004" ...
## $ kelompok: chr "Kontrol" "Kontrol" "Kontrol" "Kontrol" ...
## $ usia : num 59 48 65 47 46 60 43 59 48 40 ...
## $ jk : chr "P" "P" "P" "L" ...
## $ GDP_W0 : num 167 191 143 216 202 178 176 190 188 180 ...
## $ GDP_W1 : num 165 195 156 225 202 185 180 186 183 174 ...
## $ GDP_W2 : num 163 191 145 220 200 185 187 192 182 185 ...
## $ GDP_W3 : num 159 197 152 227 192 175 176 195 185 182 ...
colSums(is.na(dat_wide))
## id kelompok usia jk GDP_W0 GDP_W1 GDP_W2 GDP_W3
## 0 0 0 0 0 0 0 0
dat_wide$kelompok <- factor(dat_wide$kelompok, levels = kel_lab)
dat_wide$jk <- factor(dat_wide$jk)
dat_wide$id <- factor(dat_wide$id)
# Format panjang (satu baris = satu pengukuran), dibutuhkan afex/rstatix/lme4
dat_long <- dat_wide |>
pivot_longer(all_of(kol_gdp), names_to = "waktu", values_to = "gdp") |>
mutate(waktu = factor(waktu, levels = kol_gdp,
labels = paste0("W", 0:3)),
minggu = minggu[as.integer(waktu)])
head(dat_wide)
## id kelompok usia jk GDP_W0 GDP_W1 GDP_W2 GDP_W3
## 1 P001 Kontrol 59 P 167 165 163 159
## 2 P002 Kontrol 48 P 191 195 191 197
## 3 P003 Kontrol 65 P 143 156 145 152
## 4 P004 Kontrol 47 L 216 225 220 227
## 5 P005 Kontrol 46 P 202 202 200 192
## 6 P006 Kontrol 60 P 178 185 185 175
head(dat_long)
## # A tibble: 6 × 7
## id kelompok usia jk waktu gdp minggu
## <fct> <fct> <dbl> <fct> <fct> <dbl> <dbl>
## 1 P001 Kontrol 59 P W0 167 0
## 2 P001 Kontrol 59 P W1 165 4
## 3 P001 Kontrol 59 P W2 163 8
## 4 P001 Kontrol 59 P W3 159 12
## 5 P002 Kontrol 48 P W0 191 0
## 6 P002 Kontrol 48 P W1 195 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",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num [1:360] 59 59 59 59 48 48 48 48 65 65 ...
## $ jk : Factor w/ 2 levels "L","P": 2 2 2 2 2 2 2 2 2 2 ...
## $ waktu : Factor w/ 4 levels "W0","W1","W2",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ gdp : num [1:360] 167 165 163 159 191 195 191 197 143 156 ...
## $ 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 W0 gdp 30 188. 19.2
## 2 Kontrol W1 gdp 30 187. 17.4
## 3 Kontrol W2 gdp 30 188. 17.6
## 4 Kontrol W3 gdp 30 189. 18.4
## 5 Diet W0 gdp 30 195 23.0
## 6 Diet W1 gdp 30 186. 22.3
## 7 Diet W2 gdp 30 183. 22.6
## 8 Diet W3 gdp 30 177. 24.5
## 9 Diet+Olahraga W0 gdp 30 185. 21.9
## 10 Diet+Olahraga W1 gdp 30 175. 22.2
## 11 Diet+Olahraga W2 gdp 30 168. 22.6
## 12 Diet+Olahraga W3 gdp 30 162. 22.4
# Matriks kovarians & korelasi antarwaktu (seluruh subjek)
# -> memberi gambaran awal apakah sfierisitas masuk akal
S <- cov(dat_wide[, kol_gdp])
R <- cor(dat_wide[, kol_gdp])
round(S, 1); round(R, 2)
## GDP_W0 GDP_W1 GDP_W2 GDP_W3
## GDP_W0 467.7 436.5 436.4 447.8
## GDP_W1 436.5 448.9 452.0 473.1
## GDP_W2 436.4 452.0 504.5 518.4
## GDP_W3 447.8 473.1 518.4 584.9
## GDP_W0 GDP_W1 GDP_W2 GDP_W3
## GDP_W0 1.00 0.95 0.90 0.86
## GDP_W1 0.95 1.00 0.95 0.92
## GDP_W2 0.90 0.95 1.00 0.95
## GDP_W3 0.86 0.92 0.95 1.00
# Varians selisih antarpasangan waktu (inti asumsi sfierisitas)
pasangan <- combn(kol_gdp, 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_W0 - GDP_W1 GDP_W0 - GDP_W2 GDP_W0 - GDP_W3 GDP_W1 - GDP_W2 GDP_W1 - GDP_W3
## 43.7 99.4 157.0 49.4 87.7
## GDP_W2 - GDP_W3
## 52.6
# 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 = "Gula 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 pengukuran pada kelompok Diet+Olahraga?
d1 <- droplevels(filter(dat_long, kelompok == "Diet+Olahraga"))
d1w <- filter(dat_wide, kelompok == "Diet+Olahraga")
## 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 W0 gdp 0.986 0.949
## 2 W1 gdp 0.985 0.945
## 3 W2 gdp 0.974 0.652
## 4 W3 gdp 0.982 0.870
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 124.727 2.1e-31 * 0.811
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.53 0.003 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.784 2.35, 68.18 3.95e-25 * 0.857 2.57, 74.57 2.9e-27
## 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 2.35 68.18 124.727 3.95e-25 * 0.811
## 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.35, 68.18 28.29 124.73 *** .126 .811 <.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) 3579380 1 55556 29 1868.42 < 2.2e-16 ***
## waktu 8296 3 1929 87 124.73 < 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.52956 0.0034797
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.78365 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8571779 2.902781e-27
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.81 | [0.75, 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.12 | [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.98472 1868.42 1 29 < 2.2e-16 ***
## waktu 1 0.91776 100.44 3 27 9.239e-15 ***
## ---
## 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
## W0 185 3.99 29 177 193
## W1 175 4.06 29 167 184
## W2 168 4.12 29 160 177
## W3 162 4.08 29 154 171
##
## Confidence level used: 0.95
pairs(em1, adjust = "bonferroni") # semua pasangan waktu (6 perbandingan)
## contrast estimate SE df t.ratio p.value
## W0 - W1 9.27 0.681 29 13.601 <0.0001
## W0 - W2 16.43 1.320 29 12.443 <0.0001
## W0 - W3 22.27 1.380 29 16.130 <0.0001
## W1 - W2 7.17 1.220 29 5.871 <0.0001
## W1 - W3 13.00 1.250 29 10.398 <0.0001
## W2 - W3 5.83 1.300 29 4.472 0.0007
##
## 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
## W1 - W0 -9.27 0.681 29 -13.601 <0.0001
## W2 - W0 -16.43 1.320 29 -12.443 <0.0001
## W3 - W0 -22.27 1.380 29 -16.130 <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 -73.967 4.70 29 -15.746 <0.0001
## quadratic 3.433 1.44 29 2.382 0.0240
## cubic -0.767 3.45 29 -0.222 0.8256
## 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 74.9 3 3.87e-16 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.832 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 W0 W1 30 30 465 0.00000000186 1.12e-8 ****
## 2 gdp W0 W2 30 30 465 0.00000000186 1.12e-8 ****
## 3 gdp W0 W3 30 30 465 0.00000000186 1.12e-8 ****
## 4 gdp W1 W2 30 30 434 0.00000407 2.44e-5 ****
## 5 gdp W1 W3 30 30 460. 0.0000000149 8.94e-8 ****
## 6 gdp W2 W3 30 30 404. 0.000163 9.77e-4 ***
# (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 = 75.6972
## Degrees of freedom 1: 2.75
## Degrees of freedom 2: 46.67
## 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: 6 × 9
## kelompok waktu id usia jk gdp minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol W0 P003 65 P 143 0 TRUE FALSE
## 2 Kontrol W1 P004 47 L 225 4 TRUE FALSE
## 3 Kontrol W1 P026 45 L 228 4 TRUE FALSE
## 4 Kontrol W2 P003 65 P 145 8 TRUE FALSE
## 5 Kontrol W3 P004 47 L 227 12 TRUE FALSE
## 6 Kontrol W3 P026 45 L 228 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 W0 gdp 0.985 0.933
## 2 Kontrol W1 gdp 0.962 0.338
## 3 Kontrol W2 gdp 0.981 0.864
## 4 Kontrol W3 gdp 0.970 0.528
## 5 Diet W0 gdp 0.969 0.521
## 6 Diet W1 gdp 0.968 0.474
## 7 Diet W2 gdp 0.970 0.536
## 8 Diet W3 gdp 0.969 0.505
## 9 Diet+Olahraga W0 gdp 0.986 0.949
## 10 Diet+Olahraga W1 gdp 0.985 0.945
## 11 Diet+Olahraga W2 gdp 0.974 0.652
## 12 Diet+Olahraga W3 gdp 0.982 0.870
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 W0 2 87 1.22 0.300
## 2 W1 2 87 1.37 0.259
## 3 W2 2 87 1.18 0.313
## 4 W3 2 87 1.55 0.219
# (iv) Homogenitas matriks kovarians antarkelompok (Box's M; uji pada alpha = .001)
box_m(dat_wide[, kol_gdp], dat_wide$kelompok)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 21.2 0.384 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 1744.52 4.55 * .091 .095 .013
## 2 waktu 2.68, 233.49 25.19 122.45 *** .050 .585 <.001
## 3 kelompok:waktu 5.37, 233.49 25.19 37.09 *** .031 .460 <.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) 11919908 1 151773 87 6832.7710 < 2e-16 ***
## kelompok 15874 2 151773 87 4.5495 0.01321 *
## waktu 8280 3 5883 261 122.4521 < 2e-16 ***
## kelompok:waktu 5016 6 5883 261 37.0885 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.82959 0.0068065
## kelompok:waktu 0.82959 0.0068065
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.89459 < 2.2e-16 ***
## kelompok:waktu 0.89459 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.9257711 4.721357e-46
## kelompok:waktu 0.9257711 3.546979e-30
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.98743 6832.8 1 87 < 2.2e-16 ***
## kelompok 2 0.09468 4.5 2 87 0.01321 *
## waktu 1 0.77152 95.7 3 85 < 2.2e-16 ***
## kelompok:waktu 2 0.68766 15.0 6 172 8.669e-14 ***
## ---
## 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 4.550 1.30e-02 * 0.095
## 2 waktu 2.68 233.49 122.452 1.36e-44 * 0.585
## 3 kelompok:waktu 5.37 233.49 37.088 3.04e-29 * 0.460
# Ukuran efek
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.09 | [0.01, 1.00]
## waktu | 0.58 | [0.52, 1.00]
## kelompok:waktu | 0.46 | [0.38, 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.07 | [0.00, 1.00]
## waktu | 0.05 | [0.01, 1.00]
## kelompok:waktu | 0.03 | [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 0.713 0.5467
##
## kelompok = Diet:
## model term df1 df2 F.ratio p.value
## waktu 3 87 60.387 <0.0001
##
## kelompok = Diet+Olahraga:
## model term df1 df2 F.ratio p.value
## waktu 3 87 92.279 <0.0001
# Efek KELOMPOK pada tiap waktu
joint_tests(aov2, by = "waktu")
## waktu = W0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 1.812 0.1694
##
## waktu = W1:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 2.839 0.0639
##
## waktu = W2:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 7.154 0.0013
##
## waktu = W3:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 10.905 <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
## W1 - W0 -0.900 0.98 87 -0.918 1.0000
## W2 - W0 0.233 1.29 87 0.181 1.0000
## W3 - W0 0.800 1.39 87 0.577 1.0000
##
## kelompok = Diet:
## contrast estimate SE df t.ratio p.value
## W1 - W0 -9.133 0.98 87 -9.317 <0.0001
## W2 - W0 -12.300 1.29 87 -9.555 <0.0001
## W3 - W0 -17.700 1.39 87 -12.764 <0.0001
##
## kelompok = Diet+Olahraga:
## contrast estimate SE df t.ratio p.value
## W1 - W0 -9.267 0.98 87 -9.453 <0.0001
## W2 - W0 -16.433 1.29 87 -12.766 <0.0001
## W3 - W0 -22.267 1.39 87 -16.057 <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 = W0:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet -7.07 5.53 87 -1.277 0.4119
## Kontrol - (Diet+Olahraga) 3.23 5.53 87 0.584 0.8289
## Diet - (Diet+Olahraga) 10.30 5.53 87 1.861 0.1562
##
## waktu = W1:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 1.17 5.36 87 0.218 0.9742
## Kontrol - (Diet+Olahraga) 11.60 5.36 87 2.164 0.0833
## Diet - (Diet+Olahraga) 10.43 5.36 87 1.946 0.1320
##
## waktu = W2:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 5.47 5.44 87 1.006 0.5753
## Kontrol - (Diet+Olahraga) 19.90 5.44 87 3.661 0.0012
## Diet - (Diet+Olahraga) 14.43 5.44 87 2.655 0.0253
##
## waktu = W3:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 11.43 5.65 87 2.024 0.1124
## Kontrol - (Diet+Olahraga) 26.30 5.65 87 4.657 <0.0001
## Diet - (Diet+Olahraga) 14.87 5.65 87 2.632 0.0268
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## 4d. Kontras interaksi -------------------------------------------------------
# Apakah penurunan (W3 - W0) berbeda antarkelompok? -- inti pertanyaan uji klinis
em_full <- emmeans(aov2, ~ waktu * kelompok)
contrast(em_full, interaction = list(waktu = list("W3-W0" = c(-1, 0, 0, 1)),
kelompok = "pairwise"),
adjust = "holm")
## waktu_custom kelompok_pairwise estimate SE df t.ratio p.value
## W3-W0 Kontrol - Diet 18.50 1.96 87 9.433 <0.0001
## W3-W0 Kontrol - (Diet+Olahraga) 23.07 1.96 87 11.762 <0.0001
## W3-W0 Diet - (Diet+Olahraga) 4.57 1.96 87 2.329 0.0222
##
## 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 3.53 4.61 87 0.767 0.4453
## linear Diet -56.27 4.61 87 -12.210 <0.0001
## linear Diet+Olahraga -73.97 4.61 87 -16.051 <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 59.8 6.516957 87 9.176062
## 4 linear Kontrol - (Diet+Olahraga) 77.5 6.516957 87 11.892053
## 7 linear Diet - (Diet+Olahraga) 17.7 6.516957 87 2.715992
## p.value p.holm
## 1 1.957867e-14 3.915734e-14
## 4 6.254122e-20 1.876237e-19
## 7 7.969973e-03 7.969973e-03
# 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)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00726823 (tol = 0.002, component 1)
## See ?lme4::convergence and ?lme4::troubleshooting.
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 2536.0 2590.4 -1254.0 2508.0
## lmm2 16 2529.6 2591.7 -1248.8 2497.6 10.401 2 0.005515 **
## ---
## 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 162.3 81.17 2 87.00 4.5428 0.01329 *
## waktu 4756.3 1585.45 3 185.37 88.3207 < 2e-16 ***
## kelompok:waktu 2843.6 473.94 6 206.40 26.3722 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1) # korelasi intrakelas
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.950
## Unadjusted ICC: 0.806
# 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 != "W0"), 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 144.4 72.21 2 87.00 4.6196 0.0124 *
## waktu 3813.8 1271.26 3 167.67 80.9348 <2e-16 ***
## kelompok:waktu 2240.8 373.47 6 185.23 23.7495 <2e-16 ***
## ---
## 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_wide.csv", row.names = FALSE)
write.csv(dat_long, "data_gdp_long.csv", row.names = FALSE)
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Ventura 13.3.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Asia/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
## [13] readxl_1.5.0
##
## 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] ggsignif_0.6.4 cellranger_1.1.0 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