# ===========================================================================
# TUGAS BIOSTATISTIKA INTERMEDIATE
# Nama : Maria Sondang Hotmanginar Sinaga
# NIM : 2611018039
# REPEATED MEASURE ANALYSIS DENGAN R
# Data simulasi berdasarkan Yang, Park, & Jin (2026)
# Outcome: HbA1c (%) | 3 treatment x 4 waktu
# ============================================================================
# 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. DATA SIMULASI -----------------------------------------------------------
dat_wide <- read.csv("FINAL_CGM_HbA1c_Data.csv")
dat_wide$treatment <- factor(dat_wide$treatment,
levels = c("Control", "CGM only", "CGM + coaching"))
dat_wide$id <- factor(dat_wide$id)
dat_long <- dat_wide |>
pivot_longer(c(hba1c_baseline, hba1c_3m, hba1c_6m, hba1c_12m),
names_to = "time", values_to = "hba1c") |>
mutate(time = factor(time,
levels = c("hba1c_baseline","hba1c_3m","hba1c_6m","hba1c_12m"),
labels = c("Baseline","3 months","6 months","12 months")))
head(dat_wide)
## id treatment hba1c_baseline hba1c_3m hba1c_6m hba1c_12m
## 1 C01 Control 9.47 8.87 8.38 8.45
## 2 C02 Control 9.12 9.17 8.74 8.46
## 3 C03 Control 8.73 8.90 7.91 8.11
## 4 C04 Control 7.30 7.10 7.11 6.85
## 5 C05 Control 9.03 8.90 8.48 8.83
## 6 C06 Control 7.97 8.12 8.20 7.26
str(dat_long)
## tibble [360 × 4] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 90 levels "C01","C02","C03",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ treatment: Factor w/ 3 levels "Control","CGM only",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ time : Factor w/ 4 levels "Baseline","3 months",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ hba1c : num [1:360] 9.47 8.87 8.38 8.45 9.12 9.17 8.74 8.46 8.73 8.9 ...
# 2. EKSPLORASI --------------------------------------------------------------
desk <- dat_long |>
group_by(treatment, time) |>
get_summary_stats(hba1c, type = "mean_sd")
desk
## # A tibble: 12 × 6
## treatment time variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Control Baseline hba1c 30 8.39 1.11
## 2 Control 3 months hba1c 30 8.12 1.09
## 3 Control 6 months hba1c 30 8.01 1.04
## 4 Control 12 months hba1c 30 7.86 1.06
## 5 CGM only Baseline hba1c 30 8.14 1.19
## 6 CGM only 3 months hba1c 30 7.35 1.12
## 7 CGM only 6 months hba1c 30 7.76 1.07
## 8 CGM only 12 months hba1c 30 7.76 1.23
## 9 CGM + coaching Baseline hba1c 30 8.06 1.08
## 10 CGM + coaching 3 months hba1c 30 7.01 0.987
## 11 CGM + coaching 6 months hba1c 30 7.21 0.963
## 12 CGM + coaching 12 months hba1c 30 7.7 0.996
# Profile plot
ggplot(dat_long, aes(time, hba1c, colour = treatment, group = treatment)) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", size = 2.4) +
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = .12) +
labs(x = "Time", y = "HbA1c (%)", colour = "Treatment",
title = "Mean HbA1c profile (95% CI)")
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
# Spaghetti plot
ggplot(dat_long, aes(time, hba1c, group = id)) +
geom_line(alpha = .20) +
stat_summary(aes(group = treatment), fun = mean, geom = "line", linewidth = 1.1) +
stat_summary(aes(group = treatment), fun = mean, geom = "point", size = 2) +
facet_wrap(~ treatment) +
labs(x = "Time", y = "HbA1c (%)", title = "Individual trajectories and group means")
# 3. ONE-WAY REPEATED MEASURES: CGM + COACHING -------------------------------
d1 <- droplevels(filter(dat_long, treatment == "CGM + coaching"))
# Assumptions
d1 |> group_by(time) |> identify_outliers(hba1c)
## # A tibble: 3 × 6
## time id treatment hba1c is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <lgl> <lgl>
## 1 3 months I220 CGM + coaching 4.83 TRUE FALSE
## 2 3 months I230 CGM + coaching 4.79 TRUE FALSE
## 3 6 months I230 CGM + coaching 5.07 TRUE FALSE
d1 |> group_by(time) |> shapiro_test(hba1c)
## # A tibble: 4 × 4
## time variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline hba1c 0.981 0.839
## 2 3 months hba1c 0.964 0.385
## 3 6 months hba1c 0.988 0.978
## 4 12 months hba1c 0.985 0.934
ggpubr::ggqqplot(d1, "hba1c", facet.by = "time")
aov1_rs <- anova_test(data = d1, dv = hba1c, wid = id, within = time,
effect.size = "pes")
aov1_rs
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 time 3 87 65.72 2.74e-22 * 0.694
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 time 0.841 0.441
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 time 0.905 2.72, 78.76 2.12e-20 * 1.008 3.02, 87.68 2.74e-22
## p[HF]<.05
## 1 *
get_anova_table(aov1_rs, correction = "auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 time 3 87 65.72 2.74e-22 * 0.694
aov1 <- aov_ez(id = "id", dv = "hba1c", data = d1, within = "time",
anova_table = list(es = c("ges","pes"), correction = "GG"))
aov1
## Anova Table (Type 3 tests)
##
## Response: hba1c
## Effect df MSE F ges pes p.value
## 1 time 2.72, 78.76 0.11 65.72 *** .148 .694 <.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) 6741.3 1 108.455 29 1802.57 < 2.2e-16 ***
## time 20.3 3 8.966 87 65.72 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.8411 0.44137
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.90526 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 1.007872 2.736919e-22
# Post hoc
em1 <- emmeans(aov1, ~ time)
pairs(em1, adjust = "holm")
## contrast estimate SE df t.ratio p.value
## Baseline - X3.months 1.050 0.0851 29 12.345 <0.0001
## Baseline - X6.months 0.849 0.0686 29 12.384 <0.0001
## Baseline - X12.months 0.360 0.0808 29 4.454 0.0002
## X3.months - X6.months -0.201 0.0856 29 -2.345 0.0261
## X3.months - X12.months -0.690 0.0812 29 -8.498 <0.0001
## X6.months - X12.months -0.489 0.0940 29 -5.204 <0.0001
##
## P value adjustment: holm method for 6 tests
# Nonparametric sensitivity
friedman_test(d1, hba1c ~ time | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 hba1c 30 59.3 3 8.21e-13 Friedman test
d1 |> wilcox_test(hba1c ~ time, paired = TRUE, p.adjust.method = "holm")
## # 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 hba1c Baseline 3 months 30 30 464 3.73e-9 2.24e-8 ****
## 2 hba1c Baseline 6 months 30 30 464 3.73e-9 2.24e-8 ****
## 3 hba1c Baseline 12 months 30 30 415 5.59e-5 1.12e-4 ***
## 4 hba1c 3 months 6 months 30 30 124 2.44e-2 2.44e-2 *
## 5 hba1c 3 months 12 months 30 30 3 9.31e-9 3.73e-8 ****
## 6 hba1c 6 months 12 months 30 30 45 2.96e-5 8.89e-5 ****
# 4. MIXED DESIGN ANOVA: TREATMENT x TIME ------------------------------------
# Assumptions
dat_long |> group_by(treatment, time) |> shapiro_test(hba1c)
## # A tibble: 12 × 5
## treatment time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Control Baseline hba1c 0.984 0.918
## 2 Control 3 months hba1c 0.947 0.141
## 3 Control 6 months hba1c 0.968 0.473
## 4 Control 12 months hba1c 0.963 0.367
## 5 CGM only Baseline hba1c 0.976 0.726
## 6 CGM only 3 months hba1c 0.975 0.682
## 7 CGM only 6 months hba1c 0.978 0.769
## 8 CGM only 12 months hba1c 0.969 0.520
## 9 CGM + coaching Baseline hba1c 0.981 0.839
## 10 CGM + coaching 3 months hba1c 0.964 0.385
## 11 CGM + coaching 6 months hba1c 0.988 0.978
## 12 CGM + coaching 12 months hba1c 0.985 0.934
dat_long |> group_by(time) |> levene_test(hba1c ~ treatment)
## # A tibble: 4 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 87 0.0536 0.948
## 2 3 months 2 87 0.215 0.807
## 3 6 months 2 87 0.0702 0.932
## 4 12 months 2 87 0.786 0.459
box_m(dat_wide[, c("hba1c_baseline","hba1c_3m","hba1c_6m","hba1c_12m")],
dat_wide$treatment)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 13.6 0.851 20 Box's M-test for Homogeneity of Covariance Matric…
aov2 <- aov_ez(id = "id", dv = "hba1c", data = dat_long,
between = "treatment", within = "time",
anova_table = list(es = c("ges","pes"), correction = "GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: hba1c
## Effect df MSE F ges pes p.value
## 1 treatment 2, 87 4.37 2.49 + .051 .054 .089
## 2 time 2.82, 245.57 0.11 81.93 *** .056 .485 <.001
## 3 treatment:time 5.65, 245.57 0.11 16.63 *** .024 .277 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2) # includes Mauchly and GG/HF corrections
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 21794.4 1 380.20 87 4987.1088 < 2.2e-16 ***
## treatment 21.7 2 380.20 87 2.4861 0.08913 .
## time 24.3 3 25.82 261 81.9254 < 2.2e-16 ***
## treatment:time 9.9 6 25.82 261 16.6282 3.12e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.90092 0.11133
## treatment:time 0.90092 0.11133
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.9409 < 2.2e-16 ***
## treatment:time 0.9409 2.094e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.9757083 1.556878e-36
## treatment:time 0.9757083 6.821025e-16
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## treatment | 0.05 | [0.00, 1.00]
## time | 0.48 | [0.41, 1.00]
## treatment:time | 0.28 | [0.19, 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
## ------------------------------------------------
## treatment | 0.03 | [0.00, 1.00]
## time | 0.06 | [0.01, 1.00]
## treatment:time | 0.02 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Simple effects
joint_tests(aov2, by = "treatment") # effect of time within each treatment
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## treatment = Control:
## model term df1 df2 F.ratio p.value
## time 3 87 21.808 <0.0001
##
## treatment = CGM only:
## model term df1 df2 F.ratio p.value
## time 3 87 30.835 <0.0001
##
## treatment = CGM + coaching:
## model term df1 df2 F.ratio p.value
## time 3 87 65.799 <0.0001
joint_tests(aov2, by = "time") # effect of treatment at each time
## time = Baseline:
## model term df1 df2 F.ratio p.value
## treatment 2 87 0.703 0.4981
##
## time = X3.months:
## model term df1 df2 F.ratio p.value
## treatment 2 87 8.511 0.0004
##
## time = X6.months:
## model term df1 df2 F.ratio p.value
## treatment 2 87 4.758 0.0109
##
## time = X12.months:
## model term df1 df2 F.ratio p.value
## treatment 2 87 0.162 0.8503
# Pairwise treatments at each time
em_g <- emmeans(aov2, ~ treatment | time)
pairs(em_g, adjust = "tukey")
## time = Baseline:
## contrast estimate SE df t.ratio p.value
## Control - CGM only 0.2507 0.291 87 0.863 0.6652
## Control - (CGM + coaching) 0.3300 0.291 87 1.136 0.4950
## CGM only - (CGM + coaching) 0.0793 0.291 87 0.273 0.9598
##
## time = X3.months:
## contrast estimate SE df t.ratio p.value
## Control - CGM only 0.7700 0.276 87 2.794 0.0174
## Control - (CGM + coaching) 1.1097 0.276 87 4.026 0.0004
## CGM only - (CGM + coaching) 0.3397 0.276 87 1.232 0.4375
##
## time = X6.months:
## contrast estimate SE df t.ratio p.value
## Control - CGM only 0.2503 0.265 87 0.944 0.6139
## Control - (CGM + coaching) 0.7993 0.265 87 3.015 0.0093
## CGM only - (CGM + coaching) 0.5490 0.265 87 2.071 0.1019
##
## time = X12.months:
## contrast estimate SE df t.ratio p.value
## Control - CGM only 0.1000 0.284 87 0.353 0.9338
## Control - (CGM + coaching) 0.1600 0.284 87 0.564 0.8394
## CGM only - (CGM + coaching) 0.0600 0.284 87 0.212 0.9756
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# Pairwise time within each treatment
em_t <- emmeans(aov2, ~ time | treatment)
pairs(em_t, adjust = "holm")
## treatment = Control:
## contrast estimate SE df t.ratio p.value
## Baseline - X3.months 0.270333 0.0836 87 3.234 0.0069
## Baseline - X6.months 0.380000 0.0735 87 5.168 <0.0001
## Baseline - X12.months 0.530000 0.0726 87 7.296 <0.0001
## X3.months - X6.months 0.109667 0.0805 87 1.363 0.1866
## X3.months - X12.months 0.259667 0.0872 87 2.976 0.0113
## X6.months - X12.months 0.150000 0.0884 87 1.697 0.1866
##
## treatment = CGM only:
## contrast estimate SE df t.ratio p.value
## Baseline - X3.months 0.789667 0.0836 87 9.446 <0.0001
## Baseline - X6.months 0.379667 0.0735 87 5.164 <0.0001
## Baseline - X12.months 0.379333 0.0726 87 5.222 <0.0001
## X3.months - X6.months -0.410000 0.0805 87 -5.096 <0.0001
## X3.months - X12.months -0.410333 0.0872 87 -4.703 <0.0001
## X6.months - X12.months -0.000333 0.0884 87 -0.004 0.9970
##
## treatment = CGM + coaching:
## contrast estimate SE df t.ratio p.value
## Baseline - X3.months 1.050000 0.0836 87 12.560 <0.0001
## Baseline - X6.months 0.849333 0.0735 87 11.551 <0.0001
## Baseline - X12.months 0.360000 0.0726 87 4.956 <0.0001
## X3.months - X6.months -0.200667 0.0805 87 -2.494 0.0145
## X3.months - X12.months -0.690000 0.0872 87 -7.909 <0.0001
## X6.months - X12.months -0.489333 0.0884 87 -5.536 <0.0001
##
## P value adjustment: holm method for 6 tests
# Change contrasts relative to baseline
contrast(em_t, "trt.vs.ctrl", ref = 1, adjust = "holm")
## treatment = Control:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline -0.270 0.0836 87 -3.234 0.0017
## X6.months - Baseline -0.380 0.0735 87 -5.168 <0.0001
## X12.months - Baseline -0.530 0.0726 87 -7.296 <0.0001
##
## treatment = CGM only:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline -0.790 0.0836 87 -9.446 <0.0001
## X6.months - Baseline -0.380 0.0735 87 -5.164 <0.0001
## X12.months - Baseline -0.379 0.0726 87 -5.222 <0.0001
##
## treatment = CGM + coaching:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline -1.050 0.0836 87 -12.560 <0.0001
## X6.months - Baseline -0.849 0.0735 87 -11.551 <0.0001
## X12.months - Baseline -0.360 0.0726 87 -4.956 <0.0001
##
## P value adjustment: holm method for 3 tests
# 5. LINEAR MIXED MODEL -------------------------------------------------------
lmm <- lmer(hba1c ~ treatment * time + (1 | id), data = dat_long, REML = TRUE)
anova(lmm, 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.4919 0.2459 2 87 2.4861 0.08913 .
## time 24.3132 8.1044 3 261 81.9254 < 2.2e-16 ***
## treatment:time 9.8696 1.6449 6 261 16.6282 3.12e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: hba1c ~ treatment * time + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 570
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15367 -0.60310 -0.05678 0.59480 3.03537
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.06781 1.0333
## Residual 0.09892 0.3145
## Number of obs: 360, groups: id, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.78075 0.11018 87.00000 70.619 < 2e-16 ***
## treatment1 0.31417 0.15582 87.00000 2.016 0.04686 *
## treatment2 -0.02858 0.15582 87.00000 -0.183 0.85488
## time1 0.41569 0.02871 261.00000 14.478 < 2e-16 ***
## time2 -0.28764 0.02871 261.00000 -10.018 < 2e-16 ***
## time3 -0.12064 0.02871 261.00000 -4.202 3.64e-05 ***
## treatment1:time1 -0.12061 0.04060 261.00000 -2.970 0.00325 **
## treatment2:time1 -0.02853 0.04060 261.00000 -0.703 0.48295
## treatment1:time2 0.31239 0.04060 261.00000 7.693 2.94e-13 ***
## treatment2:time2 -0.11486 0.04060 261.00000 -2.829 0.00504 **
## treatment1:time3 0.03572 0.04060 261.00000 0.880 0.37980
## treatment2:time3 0.12814 0.04060 261.00000 3.156 0.00179 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn1 trtmn2 time1 time2 time3 trt1:1 trt2:1 trt1:2
## treatment1 0.000
## treatment2 0.000 -0.500
## time1 0.000 0.000 0.000
## time2 0.000 0.000 0.000 -0.333
## time3 0.000 0.000 0.000 -0.333 -0.333
## trtmnt1:tm1 0.000 0.000 0.000 0.000 0.000 0.000
## trtmnt2:tm1 0.000 0.000 0.000 0.000 0.000 0.000 -0.500
## trtmnt1:tm2 0.000 0.000 0.000 0.000 0.000 0.000 -0.333 0.167
## trtmnt2:tm2 0.000 0.000 0.000 0.000 0.000 0.000 0.167 -0.333 -0.500
## trtmnt1:tm3 0.000 0.000 0.000 0.000 0.000 0.000 -0.333 0.167 -0.333
## trtmnt2:tm3 0.000 0.000 0.000 0.000 0.000 0.000 0.167 -0.333 0.167
## trt2:2 trt1:3
## treatment1
## treatment2
## time1
## time2
## time3
## trtmnt1:tm1
## trtmnt2:tm1
## trtmnt1:tm2
## trtmnt2:tm2
## trtmnt1:tm3 0.167
## trtmnt2:tm3 -0.333 -0.500
performance::icc(lmm)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.915
## Unadjusted ICC: 0.807
# 6. SAVE --------------------------------------------------------------------
write.csv(desk, "CGM_HbA1c_descriptive_summary.csv", row.names = FALSE)