# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
#TUGAS BIOSTATISTIKA INTERMEDIATE
#NAMA :AHMAD FADHLIL AZHIM
#NIM :2611018029
# Data simulasi/modifikasi berdasarkan Ory et al. (2025), Frontiers in Public Health
# Desain: 3 treatment x 3 waktu pengukuran HbA1c
# =============================================================================
# 0. PAKET & PENGATURAN
# install.packages(c("tidyverse", "afex", "emmeans", "rstatix", "car",
# "effectsize", "lme4", "lmerTest", "performance", "ggpubr"))
suppressPackageStartupMessages({
library(dplyr)
library(tidyr)
library(ggplot2)
library(afex)
library(emmeans)
library(rstatix)
library(car)
library(effectsize)
library(lme4)
library(lmerTest)
library(performance)
library(ggpubr)
})
options(contrasts = c("contr.sum", "contr.poly"))
afex_options(emmeans_model = "multivariate")
theme_set(theme_bw(base_size = 12))
# 1. IMPORT DATA SIMULASI
# Data individual pada file ini BUKAN data mentah artikel.
# Data disimulasikan untuk latihan dengan acuan desain, titik waktu, dan pola A1c
# yang dilaporkan Ory et al. (2025).
dat_wide <- read.csv("Data_Simulasi_A1c_DSMES_3x3.csv")
dat_wide$Treatment <- factor(dat_wide$Treatment,
levels = c("TBES", "vMMWD", "Combined"))
dat_wide$ID <- factor(dat_wide$ID)
dat_long <- dat_wide |>
pivot_longer(cols = starts_with("A1c_"), names_to = "waktu", values_to = "A1c") |>
mutate(waktu = factor(waktu,
levels = c("A1c_Baseline", "A1c_M3", "A1c_M6"),
labels = c("Baseline", "M3", "M6")),
bulan = case_when(waktu == "Baseline" ~ 0,
waktu == "M3" ~ 3,
waktu == "M6" ~ 6))
head(dat_wide)
## ID Treatment A1c_Baseline A1c_M3 A1c_M6
## 1 TB01 TBES 10.44 8.61 8.36
## 2 TB02 TBES 11.16 9.70 9.46
## 3 TB03 TBES 6.92 5.50 6.32
## 4 TB04 TBES 8.99 6.94 6.61
## 5 TB05 TBES 10.75 10.56 10.32
## 6 TB06 TBES 10.38 8.19 9.38
head(dat_long)
## # A tibble: 6 × 5
## ID Treatment waktu A1c bulan
## <fct> <fct> <fct> <dbl> <dbl>
## 1 TB01 TBES Baseline 10.4 0
## 2 TB01 TBES M3 8.61 3
## 3 TB01 TBES M6 8.36 6
## 4 TB02 TBES Baseline 11.2 0
## 5 TB02 TBES M3 9.7 3
## 6 TB02 TBES M6 9.46 6
str(dat_long)
## tibble [270 × 5] (S3: tbl_df/tbl/data.frame)
## $ ID : Factor w/ 90 levels "CO01","CO02",..: 31 31 31 32 32 32 33 33 33 34 ...
## $ Treatment: Factor w/ 3 levels "TBES","vMMWD",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ waktu : Factor w/ 3 levels "Baseline","M3",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ A1c : num [1:270] 10.44 8.61 8.36 11.16 9.7 ...
## $ bulan : num [1:270] 0 3 6 0 3 6 0 3 6 0 ...
# 2. EKSPLORASI DATA
(dat_long |>
group_by(Treatment, waktu) |>
get_summary_stats(A1c, type = "mean_sd"))
## # A tibble: 9 × 6
## Treatment waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 TBES Baseline A1c 30 9.43 1.22
## 2 TBES M3 A1c 30 8.26 1.32
## 3 TBES M6 A1c 30 8.28 1.26
## 4 vMMWD Baseline A1c 30 9.14 1.24
## 5 vMMWD M3 A1c 30 8.06 1.35
## 6 vMMWD M6 A1c 30 8.08 1.32
## 7 Combined Baseline A1c 30 8.86 1.26
## 8 Combined M3 A1c 30 7.99 1.47
## 9 Combined M6 A1c 30 7.69 1.22
# Profile plot
p_profil <- ggplot(dat_long,
aes(x = bulan, y = A1c, colour = Treatment, group = Treatment)) +
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 = .35) +
scale_x_continuous(breaks = c(0, 3, 6)) +
labs(x = "Bulan", y = "HbA1c (%)", title = "Profil rerata HbA1c (±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
p_spag <- ggplot(dat_long, aes(x = bulan, y = A1c, group = ID)) +
geom_line(alpha = .25) +
stat_summary(aes(group = Treatment), fun = mean, geom = "line",
linewidth = 1.2) +
facet_wrap(~ Treatment) +
scale_x_continuous(breaks = c(0, 3, 6)) +
labs(x = "Bulan", y = "HbA1c (%)", title = "Lintasan individual dan rerata")
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Contoh fokus: kelompok TBES
d1 <- droplevels(filter(dat_long, Treatment == "TBES"))
# 3a. Uji asumsi
d1 |> group_by(waktu) |> identify_outliers(A1c)
## # A tibble: 2 × 7
## waktu ID Treatment A1c bulan is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 M3 TB03 TBES 5.5 3 TRUE FALSE
## 2 M3 TB26 TBES 5.06 3 TRUE FALSE
d1 |> group_by(waktu) |> shapiro_test(A1c)
## # A tibble: 3 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline A1c 0.947 0.142
## 2 M3 A1c 0.964 0.393
## 3 M6 A1c 0.962 0.356
ggpubr::ggqqplot(d1, "A1c", facet.by = "waktu")
aov1_rs <- anova_test(data = d1, dv = A1c, wid = ID, within = waktu,
effect.size = "pes")
aov1_rs
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2 58 80.194 2e-17 * 0.734
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.967 0.628
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05
## 1 waktu 0.968 1.94, 56.17 6.02e-17 * 1.037 2.07, 60.12 2e-17 *
get_anova_table(aov1_rs, correction = "auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2 58 80.194 2e-17 * 0.734
# 3b. ANOVA dengan afex
aov1 <- aov_ez(id = "ID", dv = "A1c", data = d1, within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov1
## Anova Table (Type 3 tests)
##
## Response: A1c
## Effect df MSE F ges pes p.value
## 1 waktu 1.94, 56.17 0.17 80.19 *** .163 .734 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov1)
## 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) 6744.4 1 129.384 29 1511.682 < 2.2e-16 ***
## waktu 27.0 2 9.776 58 80.194 < 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.96736 0.62839
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.96839 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 1.036527 2.003621e-17
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.73 | [0.63, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 3c. Pendekatan multivariat
aov1$Anova
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.98118 1511.68 1 29 < 2.2e-16 ***
## waktu 1 0.82410 65.59 2 28 2.714e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3d. Post hoc berpasangan
em1 <- emmeans(aov1, ~ waktu)
pairs(em1, adjust = "bonferroni")
## contrast estimate SE df t.ratio p.value
## Baseline - M3 1.1713 0.111 29 10.593 <0.0001
## Baseline - M6 1.1537 0.111 29 10.413 <0.0001
## M3 - M6 -0.0177 0.096 29 -0.184 1.0000
##
## P value adjustment: bonferroni method for 3 tests
contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm")
## contrast estimate SE df t.ratio p.value
## M3 - Baseline -1.17 0.111 29 -10.593 <0.0001
## M6 - Baseline -1.15 0.111 29 -10.413 <0.0001
##
## P value adjustment: holm method for 2 tests
# 3e. Alternatif nonparametrik
friedman_test(d1, A1c ~ waktu | ID)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 A1c 30 42.1 2 7.33e-10 Friedman test
friedman_effsize(d1, A1c ~ waktu | ID)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 A1c 30 0.701 Kendall W large
d1 |> wilcox_test(A1c ~ 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 A1c Baseline M3 30 30 465 0.00000000186 5.59e-9 ****
## 2 A1c Baseline M6 30 30 463 0.00000000559 1.68e-8 ****
## 3 A1c M3 M6 30 30 235 0.964 1 e+0 ns
# 4. MIXED DESIGN ANOVA (Treatment [between] x Waktu [within])
# Pertanyaan utama: apakah pola perubahan HbA1c berbeda menurut treatment?
# 4a. Uji asumsi
dat_long |> group_by(Treatment, waktu) |> identify_outliers(A1c)
## # A tibble: 2 × 7
## Treatment waktu ID A1c bulan is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 TBES M3 TB03 5.5 3 TRUE FALSE
## 2 TBES M3 TB26 5.06 3 TRUE FALSE
dat_long |> group_by(Treatment, waktu) |> shapiro_test(A1c)
## # A tibble: 9 × 5
## Treatment waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 TBES Baseline A1c 0.947 0.142
## 2 TBES M3 A1c 0.964 0.393
## 3 TBES M6 A1c 0.962 0.356
## 4 vMMWD Baseline A1c 0.969 0.503
## 5 vMMWD M3 A1c 0.956 0.248
## 6 vMMWD M6 A1c 0.962 0.357
## 7 Combined Baseline A1c 0.985 0.938
## 8 Combined M3 A1c 0.976 0.705
## 9 Combined M6 A1c 0.966 0.431
ggpubr::ggqqplot(dat_long, "A1c", ggtheme = theme_bw()) +
facet_grid(waktu ~ Treatment)
dat_long |> group_by(waktu) |> levene_test(A1c ~ Treatment)
## # A tibble: 3 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 87 0.109 0.897
## 2 M3 2 87 0.497 0.610
## 3 M6 2 87 0.118 0.889
box_m(dat_wide[, c("A1c_Baseline", "A1c_M3", "A1c_M6")],
dat_wide$Treatment)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 16.0 0.192 12 Box's M-test for Homogeneity of Covariance Matric…
# 4b. Mixed Repeated Measures ANOVA
aov2 <- aov_ez(id = "ID", dv = "A1c", data = dat_long,
between = "Treatment", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: A1c
## Effect df MSE F ges pes p.value
## 1 Treatment 2, 87 4.71 1.09 .023 .024 .340
## 2 waktu 1.79, 156.13 0.18 215.25 *** .139 .712 <.001
## 3 Treatment:waktu 3.59, 156.13 0.18 1.85 .003 .041 .129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 19145.4 1 409.65 87 4066.0191 <2e-16 ***
## Treatment 10.3 2 409.65 87 1.0923 0.340
## waktu 71.0 2 28.70 174 215.2504 <2e-16 ***
## Treatment:waktu 1.2 4 28.70 174 1.8525 0.121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.88555 0.0053716
## Treatment:waktu 0.88555 0.0053716
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.8973 <2e-16 ***
## Treatment:waktu 0.8973 0.1288
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.9149949 5.923888e-44
## Treatment:waktu 0.9149949 1.274077e-01
aov2$Anova
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.97905 4066.0 1 87 <2e-16 ***
## Treatment 2 0.02450 1.1 2 87 0.3400
## waktu 1 0.86347 272.0 2 86 <2e-16 ***
## Treatment:waktu 2 0.06794 1.5 4 174 0.1956
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix
aov2_rs <- anova_test(data = dat_long, dv = A1c, wid = ID,
between = Treatment, 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 Treatment 2.00 87.00 1.092 3.40e-01 0.024
## 2 waktu 1.79 156.13 215.250 3.71e-43 * 0.712
## 3 Treatment:waktu 3.59 156.13 1.852 1.29e-01 0.041
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------------
## Treatment | 0.02 | [0.00, 1.00]
## waktu | 0.71 | [0.66, 1.00]
## Treatment:waktu | 0.04 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 4c. Efek sederhana / post hoc
# Karena artikel sumber tidak menemukan superioritas antarmodalitas,
# bagian ini dibaca sebagai analisis pendukung, terutama bila interaksi tidak signifikan.
em2 <- emmeans(aov2, ~ waktu | Treatment)
joint_tests(aov2, by = "Treatment")
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## Treatment = TBES:
## model term df1 df2 F.ratio p.value
## waktu 2 87 99.691 <0.0001
##
## Treatment = vMMWD:
## model term df1 df2 F.ratio p.value
## waktu 2 87 85.038 <0.0001
##
## Treatment = Combined:
## model term df1 df2 F.ratio p.value
## waktu 2 87 93.525 <0.0001
joint_tests(aov2, by = "waktu")
## waktu = Baseline:
## model term df1 df2 F.ratio p.value
## Treatment 2 87 1.610 0.2059
##
## waktu = M3:
## model term df1 df2 F.ratio p.value
## Treatment 2 87 0.309 0.7346
##
## waktu = M6:
## model term df1 df2 F.ratio p.value
## Treatment 2 87 1.681 0.1922
contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm")
## Treatment = TBES:
## contrast estimate SE df t.ratio p.value
## M3 - Baseline -1.171 0.1170 87 -10.030 <0.0001
## M6 - Baseline -1.154 0.0861 87 -13.404 <0.0001
##
## Treatment = vMMWD:
## contrast estimate SE df t.ratio p.value
## M3 - Baseline -1.080 0.1170 87 -9.250 <0.0001
## M6 - Baseline -1.066 0.0861 87 -12.386 <0.0001
##
## Treatment = Combined:
## contrast estimate SE df t.ratio p.value
## M3 - Baseline -0.868 0.1170 87 -7.435 <0.0001
## M6 - Baseline -1.171 0.0861 87 -13.602 <0.0001
##
## P value adjustment: holm method for 2 tests
em2b <- emmeans(aov2, ~ Treatment | waktu)
pairs(em2b, adjust = "tukey")
## waktu = Baseline:
## contrast estimate SE df t.ratio p.value
## TBES - vMMWD 0.289 0.320 87 0.905 0.6384
## TBES - Combined 0.573 0.320 87 1.794 0.1775
## vMMWD - Combined 0.284 0.320 87 0.889 0.6488
##
## waktu = M3:
## contrast estimate SE df t.ratio p.value
## TBES - vMMWD 0.198 0.356 87 0.557 0.8430
## TBES - Combined 0.270 0.356 87 0.760 0.7287
## vMMWD - Combined 0.072 0.356 87 0.202 0.9777
##
## waktu = M6:
## contrast estimate SE df t.ratio p.value
## TBES - vMMWD 0.202 0.327 87 0.616 0.8117
## TBES - Combined 0.590 0.327 87 1.804 0.1744
## vMMWD - Combined 0.389 0.327 87 1.188 0.4638
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4d. Kontras perubahan baseline -> 6 bulan
em_full <- emmeans(aov2, ~ waktu * Treatment)
contrast(em_full,
interaction = list(waktu = list("M6-Baseline" = c(-1, 0, 1)),
Treatment = "pairwise"),
adjust = "holm")
## waktu_custom Treatment_pairwise estimate SE df t.ratio p.value
## M6-Baseline TBES - vMMWD -0.0877 0.122 87 -0.720 1.0000
## M6-Baseline TBES - Combined 0.0170 0.122 87 0.140 1.0000
## M6-Baseline vMMWD - Combined 0.1047 0.122 87 0.860 1.0000
##
## P value adjustment: holm method for 3 tests
# 5. PEMBANDING: LINEAR MIXED MODEL
lmm1 <- lmer(A1c ~ Treatment * waktu + (1 | ID), data = dat_long, REML = TRUE)
lmm2 <- lmer(A1c ~ Treatment * waktu + (1 + bulan | ID), data = dat_long,
REML = TRUE, control = lmerControl(optimizer = "bobyqa"))
## boundary (singular) fit: see help('isSingular')
anova(lmm1, lmm2, refit = FALSE)
## Data: dat_long
## Models:
## lmm1: A1c ~ Treatment * waktu + (1 | ID)
## lmm2: A1c ~ Treatment * waktu + (1 + bulan | ID)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 11 627.70 667.28 -302.85 605.70
## lmm2 13 631.35 678.13 -302.67 605.35 0.3521 2 0.8386
anova(lmm2, ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 0.360 0.180 2 87.00 1.0923 0.3400
## waktu 70.769 35.384 2 115.30 213.7465 <2e-16 ***
## Treatment:waktu 1.222 0.305 4 129.03 1.8416 0.1248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.902
## Unadjusted ICC: 0.763
performance::check_model(lmm2)
# 6. SIMPAN RINGKASAN
write.csv(dat_long, "Data_Simulasi_A1c_DSMES_3x3_LONG.csv", row.names = FALSE)