# =============================================================================
# Nama : Syamsiatul Ramadhani
# NIM : 2611018014
# DOSEN : DR.M.FATHURAHMAN,S.SI,M.Si
# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: Program manajemen stres dan skor Perceived Stress Scale (PSS-10)
# pada mahasiswa (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: studi intervensi berkelompok pada 90 mahasiswa dengan tingkat stres sedang-tinggi
# diacak ke tiga kelompok (n = 30 per kelompok):
# - Kontrol : edukasi kesehatan mental umum
# - Mindfulness : latihan mindfulness terstruktur
# - Mindfulness+AF: latihan mindfulness + aktivitas fisik terstruktur (3x/minggu)
# Skor stres PSS-10 (rentang 0-40; skor lebih tinggi = stres lebih berat) 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", "Mindfulness", "Mindfulness+AF")
minggu <- c(0, 4, 8, 12)
# Rerata populasi skor PSS-10 per kelompok x waktu
mu <- rbind(
"Kontrol" = c(26, 25.5, 25, 24.5),
"Mindfulness" = c(26, 23, 21, 19.5),
"Mindfulness+AF" = c(26, 22, 19, 16.5)
)
sd_int <- 4 # SD intersep acak (perbedaan skor stres dasar antar mahasiswa)
sd_slope <- 0.22 # SD kemiringan acak per minggu (perbedaan respons antar mahasiswa)
sd_eps <- 2.2 # 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("PSS_M", minggu)
data.frame(id = sprintf("P%03d", id),
kelompok = kel_lab[g],
usia = round(runif(n_per, 18, 25)),
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("PSS_M"), names_to = "waktu", values_to = "pss") |>
mutate(waktu = factor(waktu, levels = paste0("PSS_M", minggu),
labels = paste0("M", minggu)),
minggu = as.numeric(sub("M", "", waktu)))
head(dat_wide)
## id kelompok usia jk PSS_M0 PSS_M4 PSS_M8 PSS_M12
## 1 P001 Kontrol 20 L 27.8 22.8 30.3 25.2
## 2 P002 Kontrol 20 L 22.0 24.6 25.7 21.4
## 3 P003 Kontrol 24 P 26.1 24.8 25.6 26.6
## 4 P004 Kontrol 19 P 25.0 27.6 27.2 26.1
## 5 P005 Kontrol 22 P 19.5 19.0 18.8 23.6
## 6 P006 Kontrol 20 L 15.3 15.9 18.0 20.0
head(dat_long)
## # A tibble: 6 × 7
## id kelompok usia jk waktu pss minggu
## <fct> <fct> <dbl> <chr> <fct> <dbl> <dbl>
## 1 P001 Kontrol 20 L M0 27.8 0
## 2 P001 Kontrol 20 L M4 22.8 4
## 3 P001 Kontrol 20 L M8 30.3 8
## 4 P001 Kontrol 20 L M12 25.2 12
## 5 P002 Kontrol 20 L M0 22 0
## 6 P002 Kontrol 20 L M4 24.6 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","Mindfulness",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num [1:360] 20 20 20 20 20 20 20 20 24 24 ...
## $ 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 ...
## $ pss : num [1:360] 27.8 22.8 30.3 25.2 22 24.6 25.7 21.4 26.1 24.8 ...
## $ minggu : num [1:360] 0 4 8 12 0 4 8 12 0 4 ...
# 2. EKSPLORASI DATA
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(pss, type = "mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 pss 30 25.8 5.41
## 2 Kontrol M4 pss 30 24.6 4.89
## 3 Kontrol M8 pss 30 25.1 4.28
## 4 Kontrol M12 pss 30 23.9 4.87
## 5 Mindfulness M0 pss 30 26.6 4.24
## 6 Mindfulness M4 pss 30 23.6 4.77
## 7 Mindfulness M8 pss 30 21.5 4.37
## 8 Mindfulness M12 pss 30 20.0 4.52
## 9 Mindfulness+AF M0 pss 30 26.3 5.31
## 10 Mindfulness+AF M4 pss 30 22.4 5.41
## 11 Mindfulness+AF M8 pss 30 19.6 5.85
## 12 Mindfulness+AF M12 pss 30 17.7 6.71
# Matriks kovarians & korelasi antarwaktu (seluruh subjek, dalam kelompok)
# -> memberi gambaran awal apakah sfierisitas masuk akal
S <- cov(dat_wide[, paste0("PSS_M", minggu)])
R <- cor(dat_wide[, paste0("PSS_M", minggu)])
round(S, 1); round(R, 2)
## PSS_M0 PSS_M4 PSS_M8 PSS_M12
## PSS_M0 24.7 18.7 16.9 16.5
## PSS_M4 18.7 25.6 21.2 22.2
## PSS_M8 16.9 21.2 28.6 27.0
## PSS_M12 16.5 22.2 27.0 35.7
## PSS_M0 PSS_M4 PSS_M8 PSS_M12
## PSS_M0 1.00 0.74 0.64 0.56
## PSS_M4 0.74 1.00 0.78 0.73
## PSS_M8 0.64 0.78 1.00 0.85
## PSS_M12 0.56 0.73 0.85 1.00
# Varians selisih antarpasangan waktu (inti asumsi sfierisitas)
pasangan <- combn(paste0("PSS_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)
## PSS_M0 - PSS_M4 PSS_M0 - PSS_M8 PSS_M0 - PSS_M12 PSS_M4 - PSS_M8
## 12.9 19.4 27.3 11.7
## PSS_M4 - PSS_M12 PSS_M8 - PSS_M12
## 16.9 10.3
# Profile plot: rerata +/- 95% CI per kelompok
p_profil <- ggplot(dat_long, aes(minggu, pss, 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 = "Skor stres PSS-10", colour = "Kelompok",
title = "Profil rerata skor stres PSS-10 (± 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, pss, 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 = "Skor PSS-10", title = "Lintasan individu dan rerata kelompok")
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan: apakah skor stres berubah selama 12 minggu pada kelompok Mindfulness+AF?
d1 <- droplevels(filter(dat_long, kelompok == "Mindfulness+AF"))
d1w <- filter(dat_wide, kelompok == "Mindfulness+AF")
## 3a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per waktu (ekstrem = di luar Q1-3IQR / Q3+3IQR)
d1 |> group_by(waktu) |> identify_outliers(pss)
## [1] waktu id kelompok usia jk pss 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(pss)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 pss 0.958 0.279
## 2 M4 pss 0.950 0.171
## 3 M8 pss 0.985 0.930
## 4 M12 pss 0.964 0.398
ggpubr::ggqqplot(d1, "pss", facet.by = "waktu")
# (iii) Sfierisitas: Mauchly (dilaporkan otomatis oleh anova_test & afex)
aov1_rs <- anova_test(data = d1, dv = pss, 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 51.441 3.22e-19 * 0.639
##
## $`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]
## 1 waktu 0.708 2.12, 61.62 2.78e-14 * 0.766 2.3, 66.63 2.94e-15
## 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.12 61.62 51.441 2.78e-14 * 0.639
## 3b. ANOVA dengan afex (sumber utama laporan) -------------------------------
aov1 <- aov_ez(id = "id", dv = "pss", 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: pss
## Effect df MSE F ges pes p.value
## 1 waktu 2.12, 61.62 11.39 51.44 *** .239 .639 <.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) 55406 1 3264.4 29 492.204 < 2.2e-16 ***
## waktu 1245 3 701.9 87 51.441 < 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.5354 0.003957
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.70831 2.782e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.7658601 2.942126e-15
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.64 | [0.54, 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.23 | [0.10, 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.94436 492.20 1 29 < 2.2e-16 ***
## waktu 1 0.76422 29.17 3 27 1.272e-08 ***
## ---
## 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 26.3 0.970 29 24.3 28.2
## M4 22.4 0.987 29 20.4 24.4
## M8 19.6 1.070 29 17.4 21.8
## M12 17.7 1.230 29 15.2 20.2
##
## Confidence level used: 0.95
pairs(em1, adjust = "bonferroni") # semua pasangan waktu (6 perbandingan)
## contrast estimate SE df t.ratio p.value
## M0 - M4 3.88 0.704 29 5.508 <0.0001
## M0 - M8 6.64 0.739 29 8.990 <0.0001
## M0 - M12 8.57 1.000 29 8.571 <0.0001
## M4 - M8 2.77 0.491 29 5.637 <0.0001
## M4 - M12 4.69 0.745 29 6.301 <0.0001
## M8 - M12 1.93 0.625 29 3.084 0.0267
##
## 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 -3.88 0.704 29 -5.508 <0.0001
## M8 - M0 -6.64 0.739 29 -8.990 <0.0001
## M12 - M0 -8.57 1.000 29 -8.571 <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 -28.48 3.140 29 -9.058 <0.0001
## quadratic 1.95 0.864 29 2.258 0.0317
## cubic -0.27 1.590 29 -0.170 0.8662
## 3e. Alternatif nonparametrik ------------------------------------------------
friedman_test(d1, pss ~ waktu | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 pss 30 56.6 3 3.07e-12 Friedman test
friedman_effsize(d1, pss ~ waktu | id) # Kendall's W
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 pss 30 0.629 Kendall W large
d1 |> wilcox_test(pss ~ 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 pss M0 M4 30 30 428. 0.0000110 6.57e-5 ****
## 2 pss M0 M8 30 30 459 0.0000000261 1.56e-7 ****
## 3 pss M0 M12 30 30 458 0.0000000354 2.12e-7 ****
## 4 pss M4 M8 30 30 430. 0.00000824 4.94e-5 ****
## 5 pss M4 M12 30 30 438 0.00000226 1.35e-5 ****
## 6 pss M8 M12 30 30 371 0.00328 1.97e-2 *
# (Opsional) ANOVA robust berbasis trimmed mean -- paket WRS2
if (requireNamespace("WRS2", quietly = TRUE)) {
print(WRS2::rmanova(d1$pss, d1$waktu, d1$id, tr = 0.2))
}
## Call:
## WRS2::rmanova(y = d1$pss, groups = d1$waktu, blocks = d1$id,
## tr = 0.2)
##
## Test statistic: F = 36.5507
## Degrees of freedom 1: 2.17
## Degrees of freedom 2: 36.91
## p-value: 0
# 4. MIXED DESIGN ANOVA (Kelompok [between] x Waktu [within])
# Pertanyaan: apakah pola perubahan skor stres berbeda antarkelompok intervensi?
## 4a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per sel
dat_long |> group_by(kelompok, waktu) |> identify_outliers(pss)
## # A tibble: 3 × 9
## kelompok waktu id usia jk pss minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol M0 P030 22 P 38.2 0 TRUE FALSE
## 2 Kontrol M8 P015 23 P 14.3 8 TRUE FALSE
## 3 Kontrol M12 P015 23 P 11.2 12 TRUE FALSE
# (ii) Normalitas per sel (3 x 4 = 12 sel) dan residual model
dat_long |> group_by(kelompok, waktu) |> shapiro_test(pss)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 pss 0.988 0.977
## 2 Kontrol M4 pss 0.964 0.381
## 3 Kontrol M8 pss 0.970 0.552
## 4 Kontrol M12 pss 0.976 0.715
## 5 Mindfulness M0 pss 0.959 0.293
## 6 Mindfulness M4 pss 0.980 0.813
## 7 Mindfulness M8 pss 0.967 0.455
## 8 Mindfulness M12 pss 0.972 0.590
## 9 Mindfulness+AF M0 pss 0.958 0.279
## 10 Mindfulness+AF M4 pss 0.950 0.171
## 11 Mindfulness+AF M8 pss 0.985 0.930
## 12 Mindfulness+AF M12 pss 0.964 0.398
ggpubr::ggqqplot(dat_long, "pss", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antarkelompok pada TIAP waktu (Levene, median-centered)
dat_long |> group_by(waktu) |> levene_test(pss ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 0.960 0.387
## 2 M4 2 87 0.621 0.540
## 3 M8 2 87 2.20 0.117
## 4 M12 2 87 5.01 0.00870
# (iv) Homogenitas matriks kovarians antarkelompok (Box's M; uji pada alpha = .001)
box_m(dat_wide[, paste0("PSS_M", minggu)], dat_wide$kelompok)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 19.7 0.480 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 = "pss", data = dat_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: pss
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 84.24 4.02 * .070 .085 .021
## 2 waktu 2.60, 226.55 7.64 79.76 *** .149 .478 <.001
## 3 kelompok:waktu 5.21, 226.55 7.64 11.62 *** .049 .211 <.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) 191943 1 7328.8 87 2278.535 < 2.2e-16 ***
## kelompok 677 2 7328.8 87 4.021 0.02137 *
## waktu 1587 3 1730.8 261 79.755 < 2.2e-16 ***
## kelompok:waktu 462 6 1730.8 261 11.616 1.58e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.80905 0.0027491
## kelompok:waktu 0.80905 0.0027491
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.86799 < 2.2e-16 ***
## kelompok:waktu 0.86799 2.762e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8971592 3.753768e-33
## kelompok:waktu 0.8971592 1.466790e-10
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.96322 2278.53 1 87 < 2.2e-16 ***
## kelompok 2 0.08462 4.02 2 87 0.02137 *
## waktu 1 0.64568 51.63 3 85 < 2.2e-16 ***
## kelompok:waktu 2 0.36988 6.50 6 172 3.336e-06 ***
## ---
## 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 = pss, 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.021 2.10e-02 * 0.085
## 2 waktu 2.60 226.55 79.755 3.68e-32 * 0.478
## 3 kelompok:waktu 5.21 226.55 11.616 2.76e-10 * 0.211
# 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.48 | [0.41, 1.00]
## kelompok:waktu | 0.21 | [0.13, 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.15 | [0.08, 1.00]
## kelompok:waktu | 0.04 | [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 = "Skor PSS-10", 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.434 0.0703
##
## kelompok = Mindfulness:
## model term df1 df2 F.ratio p.value
## waktu 3 87 25.284 <0.0001
##
## kelompok = Mindfulness+AF:
## model term df1 df2 F.ratio p.value
## waktu 3 87 42.023 <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.223 0.8005
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 1.444 0.2416
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 9.792 0.0001
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 9.874 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 -1.197 0.630 87 -1.898 0.1220
## M8 - M0 -0.677 0.664 87 -1.019 0.3109
## M12 - M0 -1.893 0.813 87 -2.328 0.0667
##
## kelompok = Mindfulness:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -2.993 0.630 87 -4.748 <0.0001
## M8 - M0 -5.187 0.664 87 -7.813 <0.0001
## M12 - M0 -6.610 0.813 87 -8.126 <0.0001
##
## kelompok = Mindfulness+AF:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -3.877 0.630 87 -6.149 <0.0001
## M8 - M0 -6.643 0.664 87 -10.008 <0.0001
## M12 - M0 -8.570 0.813 87 -10.536 <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 - Mindfulness -0.863 1.30 87 -0.667 0.7835
## Kontrol - (Mindfulness+AF) -0.480 1.30 87 -0.371 0.9272
## Mindfulness - (Mindfulness+AF) 0.383 1.30 87 0.296 0.9529
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - Mindfulness 0.933 1.30 87 0.718 0.7534
## Kontrol - (Mindfulness+AF) 2.200 1.30 87 1.693 0.2137
## Mindfulness - (Mindfulness+AF) 1.267 1.30 87 0.975 0.5947
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - Mindfulness 3.647 1.26 87 2.889 0.0134
## Kontrol - (Mindfulness+AF) 5.487 1.26 87 4.347 0.0001
## Mindfulness - (Mindfulness+AF) 1.840 1.26 87 1.458 0.3162
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - Mindfulness 3.853 1.41 87 2.736 0.0204
## Kontrol - (Mindfulness+AF) 6.197 1.41 87 4.400 <0.0001
## Mindfulness - (Mindfulness+AF) 2.343 1.41 87 1.664 0.2249
##
## 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 - Mindfulness 4.72 1.15 87 4.100 0.0002
## M12-M0 Kontrol - (Mindfulness+AF) 6.68 1.15 87 5.804 <0.0001
## M12-M0 Mindfulness - (Mindfulness+AF) 1.96 1.15 87 1.704 0.0920
##
## 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 -5.16 2.61 87 -1.978 0.0511
## linear Mindfulness -22.02 2.61 87 -8.444 <0.0001
## linear Mindfulness+AF -28.48 2.61 87 -10.918 <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 - Mindfulness 16.863333 3.688574 87 4.571776
## 4 linear Kontrol - (Mindfulness+AF) 23.316667 3.688574 87 6.321323
## 7 linear Mindfulness - (Mindfulness+AF) 6.453333 3.688574 87 1.749547
## p.value p.holm
## 1 1.590100e-05 3.180201e-05
## 4 1.075498e-08 3.226493e-08
## 7 8.372277e-02 8.372277e-02
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Tidak mensyaratkan sfierisitas, menampung data hilang (MAR) dan waktu
# pengukuran yang tidak seragam.
lmm1 <- lmer(pss ~ kelompok * waktu + (1 | id), data = dat_long, REML = TRUE)
lmm2 <- lmer(pss ~ 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: pss ~ kelompok * waktu + (1 | id)
## lmm2: pss ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 1953 2007.4 -962.50 1925
## lmm2 16 1939 2001.2 -953.49 1907 18.012 2 0.0001226 ***
## ---
## 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 38.96 19.482 2 87.00 4.0210 0.02137 *
## waktu 773.13 257.711 3 185.37 52.9470 < 2.2e-16 ***
## kelompok:waktu 233.80 38.967 6 206.40 7.9968 9.075e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1) # korelasi intrakelas
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.745
## Unadjusted ICC: 0.577
# 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$pss[sample(which(dat_miss$waktu != "M0"), 30)] <- NA
lmm_miss <- lmer(pss ~ 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 37.92 18.960 2 86.976 3.8056 0.02604 *
## waktu 765.78 255.259 3 168.318 50.9923 < 2.2e-16 ***
## kelompok:waktu 249.07 41.512 6 185.986 8.2831 5.86e-08 ***
## ---
## 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$pss)])
## [1] 25
library(writexl)
write_xlsx(dat_long, "dat_long.xlsx")
write_xlsx(dat_wide, "dat_wide.xlsx")
getwd()
## [1] "D:/1. KULIAH/1. BIOSTAT/TUGAS 2"