# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
#
# Nama: Aisyah
# NIM: 2611018038
#
# Data : dataset_biostatistika_3perlakuan_3pengukuran.xlsx (sheet Data_Wide)
# Sumber : Cano-Montoya et al. (2025), J Cardiovasc Dev Dis 12(1):30
# DOI 10.3390/jcdd12010030
# Desain : 3 kelompok (Control, HIIT, RT; n = 13 per kelompok) x 3 waktu
# (baseline, minggu 4, minggu 8); outcome = tekanan darah sistolik (TDS, mmHg)
#
# CATATAN: nilai baseline (M0) dan minggu 8 (M8) berasal dari data individual
# artikel (Table 5); nilai minggu 4 (M4) DISIMULASIKAN dari rerata/SD kelompok
# (Table 2). Tuliskan hal ini dengan jelas di laporan.
#
# 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. Analisis sensitivitas hanya data asli (M0 dan M8)
# 7. Menyimpan output
# =============================================================================
# 0. PAKET & PENGATURAN
# Paket yang belum terpasang akan dipasang otomatis (butuh internet saat pertama kali).
paket <- c("readxl", "dplyr", "tidyr", "ggplot2", "afex", "emmeans", "rstatix",
"car", "effectsize", "lme4", "lmerTest", "pbkrtest", "performance",
"ggpubr")
baru <- paket[!paket %in% rownames(installed.packages())]
if (length(baru) > 0) install.packages(baru)
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
library(car) # leveneTest, Anova
library(effectsize) # eta kuadrat parsial, omega kuadrat
library(lme4) # linear mixed model
library(lmerTest) # uji F/t dengan derajat bebas Kenward-Roger
})
options(contrasts = c("contr.sum", "contr.poly")) # kontras jumlah-nol untuk SS tipe III
afex_options(emmeans_model = "multivariate") # post hoc tahan terhadap non-sfierisitas
theme_set(theme_bw(base_size = 12))
# 1. IMPOR DATA EXCEL
# Letakkan file Excel di working directory (cek dengan getwd(); ubah dengan
# setwd("folder_anda")) atau isi file_data dengan path lengkap, mis.
# "C:/Users/nama/Documents/dataset_biostatistika_3perlakuan_3pengukuran.xlsx"
file_data <- "dataset_biostatistika_3perlakuan_3pengukuran.xlsx"
if (!file.exists(file_data)) {
stop("File '", file_data, "' tidak ditemukan di: ", getwd(),
"\nPindahkan file ke folder tersebut atau gunakan setwd().")
}
kel_lab <- c("Control", "HIIT", "RT")
minggu <- c(0, 4, 8) # jarak waktu pengukuran (minggu)
kol_tds <- c("TDS_M0", "TDS_M4", "TDS_M8") # nama kolom TDS di sheet Data_Wide
dat_wide <- read_excel(file_data, sheet = "Data_Wide") |> as.data.frame()
# Pemeriksaan struktur & data hilang
stopifnot(all(c("id", "kelompok", kol_tds) %in% names(dat_wide)))
print(colSums(is.na(dat_wide[, c("id", "kelompok", kol_tds)])))
## id kelompok TDS_M0 TDS_M4 TDS_M8
## 0 0 0 0 0
dat_wide$kelompok <- factor(dat_wide$kelompok, levels = kel_lab)
stopifnot(!any(is.na(dat_wide$kelompok))) # berhenti jika ada label kelompok tak cocok
dat_wide$id <- factor(dat_wide$id)
print(table(dat_wide$kelompok))
##
## Control HIIT RT
## 13 13 13
# Format panjang (satu baris = satu pengukuran), dibutuhkan afex/rstatix/lme4
dat_long <- dat_wide |>
pivot_longer(all_of(kol_tds), names_to = "waktu", values_to = "tds") |>
mutate(waktu = factor(waktu, levels = kol_tds, labels = paste0("M", minggu)),
minggu = minggu[as.integer(waktu)])
print(head(dat_wide))
## id kelompok TDS_M0 TDS_M4 TDS_M8
## 1 Control_01 Control 144 135.1 113
## 2 Control_02 Control 148 142.2 123
## 3 Control_03 Control 158 159.2 159
## 4 Control_04 Control 158 163.9 154
## 5 Control_05 Control 152 159.3 143
## 6 Control_06 Control 146 134.9 148
print(head(dat_long))
## # A tibble: 6 × 5
## id kelompok waktu tds minggu
## <fct> <fct> <fct> <dbl> <dbl>
## 1 Control_01 Control M0 144 0
## 2 Control_01 Control M4 135. 4
## 3 Control_01 Control M8 113 8
## 4 Control_02 Control M0 148 0
## 5 Control_02 Control M4 142. 4
## 6 Control_02 Control M8 123 8
str(dat_long)
## tibble [117 × 5] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 39 levels "Control_01","Control_02",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ kelompok: Factor w/ 3 levels "Control","HIIT",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ waktu : Factor w/ 3 levels "M0","M4","M8": 1 2 3 1 2 3 1 2 3 1 ...
## $ tds : num [1:117] 144 135 113 148 142 ...
## $ minggu : num [1:117] 0 4 8 0 4 8 0 4 8 0 ...
# 2. EKSPLORASI DATA
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(tds, type = "mean_sd")
print(desk)
## # A tibble: 9 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Control M0 tds 13 137. 17.1
## 2 Control M4 tds 13 138 16.0
## 3 Control M8 tds 13 137. 13.8
## 4 HIIT M0 tds 13 141. 14.4
## 5 HIIT M4 tds 13 134 17.0
## 6 HIIT M8 tds 13 129. 15.5
## 7 RT M0 tds 13 135. 13.0
## 8 RT M4 tds 13 128 11.0
## 9 RT M8 tds 13 122. 10.4
# Matriks kovarians & korelasi antarwaktu (seluruh subjek)
S <- cov(dat_wide[, kol_tds])
R <- cor(dat_wide[, kol_tds])
print(round(S, 1)); print(round(R, 2))
## TDS_M0 TDS_M4 TDS_M8
## TDS_M0 217.8 143.5 103.8
## TDS_M4 143.5 227.8 112.5
## TDS_M8 103.8 112.5 205.6
## TDS_M0 TDS_M4 TDS_M8
## TDS_M0 1.00 0.64 0.49
## TDS_M4 0.64 1.00 0.52
## TDS_M8 0.49 0.52 1.00
# Varians selisih antarpasangan waktu (inti asumsi sfierisitas)
pasangan <- combn(kol_tds, 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 = " - ")
print(round(var_selisih, 1))
## TDS_M0 - TDS_M4 TDS_M0 - TDS_M8 TDS_M4 - TDS_M8
## 158.6 215.7 208.5
# Profile plot: rerata +/- 95% CI per kelompok
p_profil <- ggplot(dat_long, aes(minggu, tds, 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 = .4) +
scale_x_continuous(breaks = minggu) +
labs(x = "Minggu ke-", y = "Tekanan darah sistolik (mmHg)", colour = "Kelompok",
title = "Profil rerata TDS (± 95% CI)") +
theme(legend.position = "bottom")
print(p_profil)
# Spaghetti plot: lintasan tiap peserta
p_spag <- ggplot(dat_long, aes(minggu, tds, 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 = "TDS (mmHg)", title = "Lintasan individu dan rerata kelompok")
print(p_spag)
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan: apakah TDS berubah selama 8 minggu pada kelompok RT?
# (ganti "RT" dengan "HIIT" atau "Control" untuk kelompok lain)
d1 <- droplevels(filter(dat_long, kelompok == "RT"))
d1w <- filter(dat_wide, kelompok == "RT")
## 3a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per waktu (ekstrem = di luar Q1-3IQR / Q3+3IQR)
print(d1 |> group_by(waktu) |> identify_outliers(tds))
## # A tibble: 2 × 7
## waktu id kelompok tds minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 M0 RT_01 RT 169 0 TRUE FALSE
## 2 M8 RT_11 RT 145 8 TRUE FALSE
# (ii) Normalitas per waktu (Shapiro-Wilk) dan Q-Q plot
print(d1 |> group_by(waktu) |> shapiro_test(tds))
## # A tibble: 3 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 tds 0.881 0.0731
## 2 M4 tds 0.956 0.687
## 3 M8 tds 0.965 0.825
print(ggpubr::ggqqplot(d1, "tds", facet.by = "waktu"))
# (iii) Sfierisitas: Mauchly (dilaporkan otomatis oleh anova_test & afex)
aov1_rs <- anova_test(data = d1, dv = tds, wid = id, within = waktu,
effect.size = "pes")
print(aov1_rs) # ANOVA, Mauchly, koreksi GG & HF
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2 24 5.706 0.009 * 0.322
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.821 0.339
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05
## 1 waktu 0.848 1.7, 20.36 0.014 * 0.974 1.95, 23.36 0.01 *
print(get_anova_table(aov1_rs, correction = "auto")) # GG dipakai jika Mauchly p < .05
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2 24 5.706 0.009 * 0.322
## 3b. ANOVA dengan afex (sumber utama laporan) -------------------------------
aov1 <- aov_ez(id = "id", dv = "tds", data = d1, within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
print(aov1)
## Anova Table (Type 3 tests)
##
## Response: tds
## Effect df MSE F ges pes p.value
## 1 waktu 1.70, 20.36 108.61 5.71 * .180 .322 .014
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
print(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) 643849 1 2568.9 12 3007.5831 8.901e-16 ***
## waktu 1052 2 2211.8 24 5.7062 0.00939 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.82144 0.33897
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.84849 0.01385 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.9735391 0.01004706
# Ukuran efek tambahan
print(eta_squared(aov1, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.32 | [0.06, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
print(omega_squared(aov1, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## -------------------------------------------
## waktu | 0.14 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
## 3c. Pendekatan multivariat (tidak memerlukan sfierisitas) ------------------
print(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.99603 3007.58 1 12 8.901e-16 ***
## waktu 1 0.41986 3.98 2 11 0.05005 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 3d. Post hoc & kontras tren ------------------------------------------------
em1 <- emmeans(aov1, ~ waktu)
print(em1)
## waktu emmean SE df lower.CL upper.CL
## M0 135 3.61 12 127 143
## M4 128 3.06 12 121 135
## M8 122 2.88 12 116 129
##
## Confidence level used: 0.95
print(pairs(em1, adjust = "bonferroni")) # semua pasangan waktu (3)
## contrast estimate SE df t.ratio p.value
## M0 - M4 7.08 3.10 12 2.284 0.1241
## M0 - M8 12.69 4.45 12 2.851 0.0437
## M4 - M8 5.62 3.62 12 1.550 0.4413
##
## P value adjustment: bonferroni method for 3 tests
print(contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm")) # tiap waktu vs baseline
## contrast estimate SE df t.ratio p.value
## M4 - M0 -7.08 3.10 12 -2.284 0.0414
## M8 - M0 -12.69 4.45 12 -2.851 0.0292
##
## P value adjustment: holm method for 2 tests
print(contrast(em1, "poly")) # tren linear & kuadratik
## contrast estimate SE df t.ratio p.value
## linear -12.69 4.45 12 -2.851 0.0146
## quadratic 1.46 5.06 12 0.289 0.7777
## 3e. Alternatif nonparametrik ------------------------------------------------
print(friedman_test(d1, tds ~ waktu | id))
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 tds 13 2.46 2 0.292 Friedman test
print(friedman_effsize(d1, tds ~ waktu | id)) # Kendall's W
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 tds 13 0.0947 Kendall W small
print(d1 |> wilcox_test(tds ~ waktu, paired = TRUE, p.adjust.method = "bonferroni"))
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 tds M0 M4 13 13 73 0.0574 0.172 ns
## 2 tds M0 M8 13 13 78 0.0215 0.0645 ns
## 3 tds M4 M8 13 13 66 0.168 0.503 ns
# (Opsional) ANOVA robust berbasis trimmed mean -- paket WRS2
if (requireNamespace("WRS2", quietly = TRUE)) {
print(WRS2::rmanova(d1$tds, d1$waktu, d1$id, tr = 0.2))
}
## Call:
## WRS2::rmanova(y = d1$tds, groups = d1$waktu, blocks = d1$id,
## tr = 0.2)
##
## Test statistic: F = 6.8734
## Degrees of freedom 1: 2
## Degrees of freedom 2: 16
## p-value: 0.00701
# 4. MIXED DESIGN ANOVA (Kelompok [between] x Waktu [within])
# Pertanyaan: apakah pola perubahan TDS berbeda antarkelompok?
## 4a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per sel
print(dat_long |> group_by(kelompok, waktu) |> identify_outliers(tds))
## # A tibble: 8 × 7
## kelompok waktu id tds minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 Control M4 Control_03 159. 4 TRUE FALSE
## 2 Control M4 Control_04 164. 4 TRUE FALSE
## 3 Control M4 Control_05 159. 4 TRUE FALSE
## 4 Control M4 Control_10 108. 4 TRUE FALSE
## 5 HIIT M8 HIIT_08 167 8 TRUE TRUE
## 6 HIIT M8 HIIT_13 102 8 TRUE FALSE
## 7 RT M0 RT_01 169 0 TRUE FALSE
## 8 RT M8 RT_11 145 8 TRUE FALSE
# (ii) Normalitas per sel (3 x 3 = 9 sel) dan Q-Q plot
print(dat_long |> group_by(kelompok, waktu) |> shapiro_test(tds))
## # A tibble: 9 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Control M0 tds 0.948 0.562
## 2 Control M4 tds 0.954 0.656
## 3 Control M8 tds 0.965 0.825
## 4 HIIT M0 tds 0.975 0.944
## 5 HIIT M4 tds 0.973 0.927
## 6 HIIT M8 tds 0.918 0.233
## 7 RT M0 tds 0.881 0.0731
## 8 RT M4 tds 0.956 0.687
## 9 RT M8 tds 0.965 0.825
print(ggpubr::ggqqplot(dat_long, "tds", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok))
# (iii) Homogenitas varians antarkelompok pada TIAP waktu (Levene)
print(dat_long |> group_by(waktu) |> levene_test(tds ~ kelompok))
## # A tibble: 3 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 36 0.772 0.469
## 2 M4 2 36 0.873 0.426
## 3 M8 2 36 0.759 0.476
# (iv) Homogenitas matriks kovarians antarkelompok (Box's M; uji pada alpha = .001)
print(box_m(dat_wide[, kol_tds], dat_wide$kelompok))
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 9.15 0.690 12 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 = "tds", data = dat_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
print(aov2)
## Anova Table (Type 3 tests)
##
## Response: tds
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 36 439.53 1.75 .064 .089 .188
## 2 waktu 1.92, 69.11 96.36 7.41 ** .057 .171 .001
## 3 kelompok:waktu 3.84, 69.11 96.36 1.95 .031 .098 .114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
print(summary(aov2)) # Mauchly, epsilon GG/HF, p terkoreksi
## Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): HF eps > 1
## treated as 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 2081867 1 15823.2 36 4736.5538 < 2.2e-16 ***
## kelompok 1541 2 15823.2 36 1.7535 0.187641
## waktu 1371 2 6659.6 72 7.4118 0.001183 **
## kelompok:waktu 722 4 6659.6 72 1.9520 0.111083
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.95821 0.47378
## kelompok:waktu 0.95821 0.47378
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.95989 0.001401 **
## kelompok:waktu 0.95989 0.114221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 1.012783 0.0011831
## kelompok:waktu 1.012783 0.1110832
print(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.99246 4736.6 1 36 < 2e-16 ***
## kelompok 2 0.08877 1.8 2 36 0.18764
## waktu 1 0.29714 7.4 2 35 0.00209 **
## kelompok:waktu 2 0.19248 1.9 4 72 0.11687
## ---
## 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 = tds, wid = id,
between = kelompok, within = waktu, effect.size = "pes",
type = 3)
print(get_anova_table(aov2_rs, correction = "GG"))
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 kelompok 2.00 36.00 1.753 0.188 0.089
## 2 waktu 1.92 69.11 7.412 0.001 * 0.171
## 3 kelompok:waktu 3.84 69.11 1.952 0.114 0.098
# Ukuran efek
print(eta_squared(aov2, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.09 | [0.00, 1.00]
## waktu | 0.17 | [0.05, 1.00]
## kelompok:waktu | 0.10 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
print(omega_squared(aov2, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## ------------------------------------------------
## kelompok | 0.04 | [0.00, 1.00]
## waktu | 0.05 | [0.00, 1.00]
## kelompok:waktu | 0.01 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Plot interaksi dari model
print(afex_plot(aov2, x = "waktu", trace = "kelompok", error = "within",
mapping = c("colour", "shape", "linetype")) +
labs(y = "TDS (mmHg)", 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)
print(joint_tests(aov2, by = "kelompok"))
## kelompok = Control:
## model term df1 df2 F.ratio p.value
## waktu 2 36 0.082 0.9211
##
## kelompok = HIIT:
## model term df1 df2 F.ratio p.value
## waktu 2 36 5.850 0.0063
##
## kelompok = RT:
## model term df1 df2 F.ratio p.value
## waktu 2 36 5.968 0.0058
# Efek KELOMPOK pada tiap waktu
print(joint_tests(aov2, by = "waktu"))
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 36 0.562 0.5749
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 36 1.482 0.2407
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 36 3.774 0.0325
# Post hoc: tiap waktu vs baseline di dalam tiap kelompok
print(contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm"))
## kelompok = Control:
## contrast estimate SE df t.ratio p.value
## M4 - M0 1.3077 3.40 36 0.384 1.0000
## M8 - M0 0.0769 3.81 36 0.020 1.0000
##
## kelompok = HIIT:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -7.0769 3.40 36 -2.080 0.0447
## M8 - M0 -12.5385 3.81 36 -3.289 0.0045
##
## kelompok = RT:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -7.0769 3.40 36 -2.080 0.0447
## M8 - M0 -12.6923 3.81 36 -3.330 0.0040
##
## P value adjustment: holm method for 2 tests
# Post hoc: perbandingan antarkelompok pada tiap waktu
em2b <- emmeans(aov2, ~ kelompok | waktu)
print(pairs(em2b, adjust = "tukey"))
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Control - HIIT -4.38 5.86 36 -0.749 0.7363
## Control - RT 1.62 5.86 36 0.276 0.9590
## HIIT - RT 6.00 5.86 36 1.025 0.5664
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Control - HIIT 4.00 5.85 36 0.684 0.7742
## Control - RT 10.00 5.85 36 1.710 0.2152
## HIIT - RT 6.00 5.85 36 1.026 0.5654
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Control - HIIT 8.23 5.25 36 1.567 0.2729
## Control - RT 14.38 5.25 36 2.738 0.0253
## HIIT - RT 6.15 5.25 36 1.171 0.4776
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## 4d. Kontras interaksi -------------------------------------------------------
# Apakah penurunan (M8 - M0) berbeda antarkelompok?
em_full <- emmeans(aov2, ~ waktu * kelompok)
print(contrast(em_full, interaction = list(waktu = list("M8-M0" = c(-1, 0, 1)),
kelompok = "pairwise"),
adjust = "holm"))
## waktu_custom kelompok_pairwise estimate SE df t.ratio p.value
## M8-M0 Control - HIIT 12.615 5.39 36 2.340 0.0700
## M8-M0 Control - RT 12.769 5.39 36 2.369 0.0700
## M8-M0 HIIT - RT 0.154 5.39 36 0.029 0.9774
##
## P value adjustment: holm method for 3 tests
# Tren linear per kelompok (dengan 3 waktu: baris linear = 1, 3, 5) dan perbandingannya
print(contrast(em2, "poly")[c(1, 3, 5)])
## contrast kelompok estimate SE df t.ratio p.value
## linear Control 0.0769 3.81 36 0.020 0.9840
## linear HIIT -12.5385 3.81 36 -3.289 0.0023
## linear RT -12.6923 3.81 36 -3.330 0.0020
tren_int <- summary(contrast(em_full, interaction = c(waktu = "poly", kelompok = "pairwise"),
adjust = "none"))
tren_lin <- subset(tren_int, waktu_poly == "linear") # apakah laju perubahan linear berbeda?
tren_lin$p.holm <- p.adjust(tren_lin$p.value, "holm") # koreksi Holm untuk 3 perbandingan
print(tren_lin)
## waktu_poly kelompok_pairwise estimate SE df t.ratio p.value
## 1 linear Control - HIIT 12.6153846 5.391083 36 2.34004653 0.02494305
## 3 linear Control - RT 12.7692308 5.391083 36 2.36858368 0.02334418
## 5 linear HIIT - RT 0.1538462 5.391083 36 0.02853715 0.97739135
## p.holm
## 1 0.07003253
## 3 0.07003253
## 5 0.97739135
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Tidak mensyaratkan sfierisitas dan menampung data hilang (MAR).
# Dengan hanya 3 pengukuran per peserta, model intersep acak (lmm1) lebih
# stabil; model slope acak (lmm2) dapat menghasilkan peringatan singular fit.
lmm1 <- lmer(tds ~ kelompok * waktu + (1 | id), data = dat_long, REML = TRUE)
print(anova(lmm1, 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 324.37 162.19 2 36 1.7535 0.187641
## waktu 1371.09 685.55 2 72 7.4118 0.001183 **
## kelompok:waktu 722.19 180.55 4 72 1.9520 0.111083
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(lmm1))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: tds ~ kelompok * waktu + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 887.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.63822 -0.49859 0.04162 0.47457 1.90342
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 115.68 10.755
## Residual 92.49 9.617
## Number of obs: 117, groups: id, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 133.39316 1.93822 36.00000 68.823 < 2e-16 ***
## kelompok1 3.76068 2.74105 36.00000 1.372 0.17856
## kelompok2 1.14530 2.74105 36.00000 0.418 0.67855
## waktu1 4.22222 1.25742 72.00000 3.358 0.00126 **
## waktu2 -0.05983 1.25742 72.00000 -0.048 0.96218
## kelompok1:waktu1 -4.68376 1.77826 72.00000 -2.634 0.01033 *
## kelompok2:waktu1 2.31624 1.77826 72.00000 1.303 0.19689
## kelompok1:waktu2 0.90598 1.77826 72.00000 0.509 0.61198
## kelompok2:waktu2 -0.47863 1.77826 72.00000 -0.269 0.78858
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) klmpk1 klmpk2 waktu1 waktu2 klm1:1 klm2:1 klm1:2
## kelompok1 0.000
## kelompok2 0.000 -0.500
## waktu1 0.000 0.000 0.000
## waktu2 0.000 0.000 0.000 -0.500
## klmpk1:wkt1 0.000 0.000 0.000 0.000 0.000
## klmpk2:wkt1 0.000 0.000 0.000 0.000 0.000 -0.500
## klmpk1:wkt2 0.000 0.000 0.000 0.000 0.000 -0.500 0.250
## klmpk2:wkt2 0.000 0.000 0.000 0.000 0.000 0.250 -0.500 -0.500
print(performance::icc(lmm1)) # korelasi intrakelas
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.556
## Unadjusted ICC: 0.483
lmm2 <- lmer(tds ~ kelompok * waktu + (1 + minggu | id), data = dat_long, REML = TRUE,
control = lmerControl(optimizer = "bobyqa"))
## boundary (singular) fit: see help('isSingular')
print(anova(lmm1, lmm2, refit = FALSE)) # uji rasio kemungkinan struktur acak
## Data: dat_long
## Models:
## lmm1: tds ~ kelompok * waktu + (1 | id)
## lmm2: tds ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 11 909.80 940.18 -443.90 887.80
## lmm2 13 912.37 948.28 -443.19 886.37 1.426 2 0.4902
# Diagnostik residual LMM (normalitas & homogenitas)
par(mfrow = c(1, 2))
qqnorm(resid(lmm1), main = "Q-Q residual"); qqline(resid(lmm1))
plot(fitted(lmm1), resid(lmm1), xlab = "Nilai prediksi", ylab = "Residual",
main = "Residual vs prediksi"); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# 6. ANALISIS SENSITIVITAS: HANYA DATA ASLI ARTIKEL (M0 DAN M8)
# Karena M4 disimulasikan, hasil utama sebaiknya didukung analisis tanpa M4.
dat_asli <- dat_long |> filter(waktu %in% c("M0", "M8")) |> droplevels()
aov3 <- aov_ez(id = "id", dv = "tds", data = dat_asli,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes")))
print(aov3) # mixed ANOVA 3 x 2
## Anova Table (Type 3 tests)
##
## Response: tds
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 36 307.85 1.47 .059 .076 .243
## 2 waktu 1, 36 94.46 14.51 *** .086 .287 <.001
## 3 kelompok:waktu 2, 36 94.46 3.70 * .046 .170 .035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
print(contrast(emmeans(aov3, ~ waktu | kelompok), "revpairwise", adjust = "holm"))
## kelompok = Control:
## contrast estimate SE df t.ratio p.value
## M8 - M0 0.0769 3.81 36 0.020 0.9840
##
## kelompok = HIIT:
## contrast estimate SE df t.ratio p.value
## M8 - M0 -12.5385 3.81 36 -3.289 0.0023
##
## kelompok = RT:
## contrast estimate SE df t.ratio p.value
## M8 - M0 -12.6923 3.81 36 -3.330 0.0020
# ANCOVA: TDS minggu 8 dengan baseline sebagai kovariat
m_ancova <- lm(TDS_M8 ~ TDS_M0 + kelompok, data = dat_wide)
print(car::Anova(m_ancova, type = 3))
## Anova Table (Type III tests)
##
## Response: TDS_M8
## Sum Sq Df F value Pr(>F)
## (Intercept) 1682.7 1 12.7477 0.0010597 **
## TDS_M0 1838.7 1 13.9293 0.0006723 ***
## kelompok 1311.8 2 4.9691 0.0126010 *
## Residuals 4620.0 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(pairs(emmeans(m_ancova, ~ kelompok), adjust = "tukey"))
## contrast estimate SE df t.ratio p.value
## Control - HIIT 10.33 4.54 35 2.275 0.0728
## Control - RT 13.61 4.51 35 3.017 0.0128
## HIIT - RT 3.28 4.57 35 0.718 0.7546
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 7. SIMPAN OUTPUT & SESSION INFO
write.csv(desk, "Ringkasan_Deskriptif_TDS.csv", row.names = FALSE)
write.csv(dat_long, "data_tds_long.csv", row.names = FALSE)
ggsave("Profil_TDS.png", p_profil, width = 7, height = 4.5, dpi = 300)
ggsave("Spaghetti_TDS.png", p_spag, width = 9, height = 4.5, dpi = 300)
message("Selesai. Output disimpan di: ", getwd())
## Selesai. Output disimpan di: C:/Users/dianw/Downloads
print(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
## [13] readxl_1.5.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 farver_2.1.2 S7_0.2.2
## [4] fastmap_1.2.0 reshape_0.8.10 bayestestR_0.19.0
## [7] digest_0.6.39 rpart_4.1.27 estimability_2.0.0
## [10] lifecycle_1.0.5 cluster_2.1.8.2 magrittr_2.0.5
## [13] compiler_4.6.1 rlang_1.3.0 Hmisc_5.3-0
## [16] sass_0.4.10 tools_4.6.1 utf8_1.2.6
## [19] yaml_2.3.12 data.table_1.18.6.1 ggsignif_0.6.4
## [22] knitr_1.51 labeling_0.4.3 htmlwidgets_1.6.4
## [25] plyr_1.8.9 RColorBrewer_1.1-3 abind_1.4-8
## [28] withr_3.0.3 foreign_0.8-91 purrr_1.2.2
## [31] numDeriv_2016.8-1.1 nnet_7.3-20 grid_4.6.1
## [34] datawizard_1.4.0 ggpubr_1.0.0 colorspace_2.1-3
## [37] scales_1.4.0 MASS_7.3-65 insight_1.5.4
## [40] cli_3.6.6 mvtnorm_1.4-2 rmarkdown_2.31
## [43] reformulas_0.4.4 generics_0.1.4 performance_0.18.2
## [46] rstudioapi_0.19.0 reshape2_1.4.5 parameters_0.29.3
## [49] minqa_1.2.8 cachem_1.1.0 stringr_1.6.0
## [52] splines_4.6.1 parallel_4.6.1 WRS2_1.1-7
## [55] cellranger_1.1.0 base64enc_0.1-6 vctrs_0.7.3
## [58] boot_1.3-32 jsonlite_2.0.0 pbkrtest_0.5.5
## [61] Formula_1.2-6 htmlTable_2.5.0 jquerylib_0.1.4
## [64] glue_1.8.1 nloptr_2.2.1 stringi_1.8.9
## [67] gtable_0.3.6 tibble_3.3.1 pillar_1.11.1
## [70] htmltools_0.5.9 R6_2.6.1 Rdpack_2.6.6
## [73] evaluate_1.0.5 lattice_0.22-9 rbibutils_2.4.1
## [76] backports_1.5.1 broom_1.0.13 bslib_0.12.0
## [79] Rcpp_1.1.2 gridExtra_2.3.1 nlme_3.1-169
## [82] checkmate_2.3.4 xfun_0.60 pkgconfig_2.0.3