# ================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Nama : Daivy Putri Anzelina Marbun
# Nim : 2611018040
# Data simulasi berdasarkan Bailey et al. (2026; online 2025)
# MODEL Study: 3 treatment x 4 repeated measurements
# Outcome: healthy eating (days/week)
# ================================================================
# 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
dat_wide <- read.csv("FINAL_MODEL_Data.csv")
dat_wide$treatment <- factor(dat_wide$treatment,
levels=c("EM only","HC + EM","TM + EM"))
dat_wide$id <- factor(dat_wide$id)
dat_long <- dat_wide |>
pivot_longer(c(healthy_baseline,healthy_3m,healthy_6m,healthy_12m),
names_to="time",values_to="healthy") |>
mutate(time=factor(time,
levels=c("healthy_baseline","healthy_3m","healthy_6m","healthy_12m"),
labels=c("Baseline","3 months","6 months","12 months")))
# 2. EKSPLORASI
desk <- dat_long |> group_by(treatment,time) |>
get_summary_stats(healthy,type="mean_sd")
print(desk)
## # A tibble: 12 × 6
## treatment time variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 EM only Baseline healthy 30 3.2 0.703
## 2 EM only 3 months healthy 30 3.40 0.855
## 3 EM only 6 months healthy 30 3.6 0.768
## 4 EM only 12 months healthy 30 3.87 0.882
## 5 HC + EM Baseline healthy 30 3.2 0.639
## 6 HC + EM 3 months healthy 30 3.5 0.654
## 7 HC + EM 6 months healthy 30 3.85 0.907
## 8 HC + EM 12 months healthy 30 4.19 0.956
## 9 TM + EM Baseline healthy 30 3.20 0.704
## 10 TM + EM 3 months healthy 30 3.6 0.847
## 11 TM + EM 6 months healthy 30 4.05 0.792
## 12 TM + EM 12 months healthy 30 4.56 0.742
ggplot(dat_long,aes(time,healthy,colour=treatment,group=treatment))+
stat_summary(fun=mean,geom="line",linewidth=1)+
stat_summary(fun=mean,geom="point",size=2.5)+
labs(y="Healthy eating (days/week)",x="Time")
ggplot(dat_long,aes(time,healthy,group=id,colour=treatment))+
geom_line(alpha=.18)+facet_wrap(~treatment)+
stat_summary(aes(group=treatment),fun=mean,geom="line",linewidth=1.2)
# 3. ONE-WAY REPEATED MEASURE: TM + EM
d1 <- droplevels(filter(dat_long,treatment=="TM + EM"))
d1 |> group_by(time) |> shapiro_test(healthy)
## # A tibble: 4 × 4
## time variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline healthy 0.975 0.678
## 2 3 months healthy 0.944 0.117
## 3 6 months healthy 0.978 0.773
## 4 12 months healthy 0.979 0.799
aov1_rs <- anova_test(data=d1,dv=healthy,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 34.422 9.2e-15 * 0.543
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 time 0.812 0.329
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 time 0.88 2.64, 76.57 2.92e-13 * 0.976 2.93, 84.94 1.82e-14
## 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 34.422 9.2e-15 * 0.543
aov1 <- aov_ez(id="id",dv="healthy",data=d1,within="time",
anova_table=list(es=c("ges","pes"),correction="GG"))
summary(aov1)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1781.55 1 43.301 29 1193.149 < 2.2e-16 ***
## time 30.87 3 26.005 87 34.422 9.197e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.81207 0.32942
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.88009 2.918e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.9763346 1.818781e-14
em1 <- emmeans(aov1,~time)
pairs(em1,adjust="holm")
## contrast estimate SE df t.ratio p.value
## Baseline - X3.months -0.400 0.166 29 -2.413 0.0224
## Baseline - X6.months -0.850 0.150 29 -5.676 <0.0001
## Baseline - X12.months -1.360 0.146 29 -9.298 <0.0001
## X3.months - X6.months -0.451 0.111 29 -4.054 0.0010
## X3.months - X12.months -0.960 0.138 29 -6.947 <0.0001
## X6.months - X12.months -0.509 0.130 29 -3.924 0.0010
##
## P value adjustment: holm method for 6 tests
friedman_test(d1,healthy~time|id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 healthy 30 44.6 3 0.00000000113 Friedman test
# 4. MIXED DESIGN ANOVA: treatment x time
# asumsi
dat_long |> group_by(treatment,time) |> shapiro_test(healthy)
## # A tibble: 12 × 5
## treatment time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 EM only Baseline healthy 0.970 0.550
## 2 EM only 3 months healthy 0.979 0.794
## 3 EM only 6 months healthy 0.957 0.252
## 4 EM only 12 months healthy 0.984 0.917
## 5 HC + EM Baseline healthy 0.964 0.400
## 6 HC + EM 3 months healthy 0.982 0.868
## 7 HC + EM 6 months healthy 0.988 0.974
## 8 HC + EM 12 months healthy 0.977 0.743
## 9 TM + EM Baseline healthy 0.975 0.678
## 10 TM + EM 3 months healthy 0.944 0.117
## 11 TM + EM 6 months healthy 0.978 0.773
## 12 TM + EM 12 months healthy 0.979 0.799
dat_long |> group_by(time) |> levene_test(healthy~treatment)
## # A tibble: 4 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 87 0.388 0.680
## 2 3 months 2 87 0.927 0.400
## 3 6 months 2 87 0.468 0.628
## 4 12 months 2 87 0.890 0.414
box_m(dat_wide[,c("healthy_baseline","healthy_3m","healthy_6m","healthy_12m")],
dat_wide$treatment)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 14.2 0.818 20 Box's M-test for Homogeneity of Covariance Matric…
aov2 <- aov_ez(id="id",dv="healthy",data=dat_long,
between="treatment",within="time",
anova_table=list(es=c("ges","pes"),correction="GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: healthy
## Effect df MSE F ges pes p.value
## 1 treatment 2, 87 1.58 2.13 .030 .047 .125
## 2 time 2.85, 248.36 0.33 54.13 *** .188 .384 <.001
## 3 treatment:time 5.71, 248.36 0.33 2.17 * .018 .048 .049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2) # termasuk Mauchly + GG/HF
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 4889.4 1 137.543 87 3092.6909 < 2e-16 ***
## treatment 6.7 2 137.543 87 2.1328 0.12467
## time 50.7 3 81.524 261 54.1294 < 2e-16 ***
## treatment:time 4.1 6 81.524 261 2.1712 0.04616 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.92444 0.24116
## treatment:time 0.92444 0.24116
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.95159 < 2e-16 ***
## treatment:time 0.95159 0.04941 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.987258 6.346626e-27
## treatment:time 0.987258 4.699463e-02
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.38 | [0.31, 1.00]
## treatment:time | 0.05 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# simple effects dan post hoc
joint_tests(aov2,by="treatment")
## treatment = EM only:
## model term df1 df2 F.ratio p.value
## time 3 87 7.039 0.0003
##
## treatment = HC + EM:
## model term df1 df2 F.ratio p.value
## time 3 87 15.672 <0.0001
##
## treatment = TM + EM:
## model term df1 df2 F.ratio p.value
## time 3 87 29.389 <0.0001
joint_tests(aov2,by="time")
## time = Baseline:
## model term df1 df2 F.ratio p.value
## treatment 2 87 0.000 1.0000
##
## time = X3.months:
## model term df1 df2 F.ratio p.value
## treatment 2 87 0.477 0.6223
##
## time = X6.months:
## model term df1 df2 F.ratio p.value
## treatment 2 87 2.253 0.1111
##
## time = X12.months:
## model term df1 df2 F.ratio p.value
## treatment 2 87 4.785 0.0107
em2 <- emmeans(aov2,~time|treatment)
contrast(em2,"trt.vs.ctrl",ref=1,adjust="holm")
## treatment = EM only:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline 0.201 0.150 87 1.342 0.1831
## X6.months - Baseline 0.400 0.157 87 2.543 0.0255
## X12.months - Baseline 0.670 0.153 87 4.373 0.0001
##
## treatment = HC + EM:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline 0.300 0.150 87 2.005 0.0481
## X6.months - Baseline 0.650 0.157 87 4.130 0.0002
## X12.months - Baseline 0.990 0.153 87 6.460 <0.0001
##
## treatment = TM + EM:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline 0.400 0.150 87 2.668 0.0091
## X6.months - Baseline 0.850 0.157 87 5.406 <0.0001
## X12.months - Baseline 1.360 0.153 87 8.869 <0.0001
##
## P value adjustment: holm method for 3 tests
em2b <- emmeans(aov2,~treatment|time)
pairs(em2b,adjust="tukey")
## time = Baseline:
## contrast estimate SE df t.ratio p.value
## EM only - (HC + EM) 0.000000 0.176 87 0.000 1.0000
## EM only - (TM + EM) -0.000667 0.176 87 -0.004 1.0000
## (HC + EM) - (TM + EM) -0.000667 0.176 87 -0.004 1.0000
##
## time = X3.months:
## contrast estimate SE df t.ratio p.value
## EM only - (HC + EM) -0.099333 0.204 87 -0.487 0.8778
## EM only - (TM + EM) -0.199333 0.204 87 -0.977 0.5935
## (HC + EM) - (TM + EM) -0.100000 0.204 87 -0.490 0.8763
##
## time = X6.months:
## contrast estimate SE df t.ratio p.value
## EM only - (HC + EM) -0.249667 0.213 87 -1.173 0.4724
## EM only - (TM + EM) -0.451000 0.213 87 -2.119 0.0919
## (HC + EM) - (TM + EM) -0.201333 0.213 87 -0.946 0.6129
##
## time = X12.months:
## contrast estimate SE df t.ratio p.value
## EM only - (HC + EM) -0.320000 0.223 87 -1.433 0.3284
## EM only - (TM + EM) -0.690000 0.223 87 -3.091 0.0075
## (HC + EM) - (TM + EM) -0.370000 0.223 87 -1.657 0.2275
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# perubahan baseline -> 12 bulan antar-treatment
change_dat <- dat_wide |>
mutate(change12=healthy_12m-healthy_baseline)
anova_test(data=change_dat,dv=change12,between=treatment,effect.size="pes")
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 treatment 2 87 5.063 0.008 * 0.104
change_dat |> tukey_hsd(change12~treatment)
## # A tibble: 3 × 9
## term group1 group2 null.value estimate conf.low conf.high p.adj
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 treatment EM only HC + EM 0 0.32 -0.197 0.837 0.307
## 2 treatment EM only TM + EM 0 0.689 0.172 1.21 0.00574
## 3 treatment HC + EM TM + EM 0 0.369 -0.148 0.886 0.210
## # ℹ 1 more variable: p.adj.signif <chr>
# 5. LINEAR MIXED MODEL
lmm <- lmer(healthy ~ 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 1.332 0.6662 2 87 2.1328 0.12467
## time 50.722 16.9074 3 261 54.1294 < 2e-16 ***
## treatment:time 4.069 0.6782 6 261 2.1712 0.04616 *
## ---
## 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: healthy ~ treatment * time + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 781.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.60778 -0.59199 0.01982 0.58991 2.34618
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.3172 0.5632
## Residual 0.3124 0.5589
## Number of obs: 360, groups: id, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.685e+00 6.627e-02 8.700e+01 55.612 < 2e-16 ***
## treatment1 -1.675e-01 9.372e-02 8.700e+01 -1.787 0.077376 .
## treatment2 -2.500e-04 9.372e-02 8.700e+01 -0.003 0.997878
## time1 -4.851e-01 5.102e-02 2.610e+02 -9.508 < 2e-16 ***
## time2 -1.848e-01 5.102e-02 2.610e+02 -3.622 0.000352 ***
## time3 1.482e-01 5.102e-02 2.610e+02 2.905 0.003984 **
## treatment1:time1 1.673e-01 7.215e-02 2.610e+02 2.318 0.021200 *
## treatment2:time1 2.778e-05 7.215e-02 2.610e+02 0.000 0.999693
## treatment1:time2 6.794e-02 7.215e-02 2.610e+02 0.942 0.347223
## treatment2:time2 2.778e-05 7.215e-02 2.610e+02 0.000 0.999693
## treatment1:time3 -6.606e-02 7.215e-02 2.610e+02 -0.916 0.360770
## treatment2:time3 1.636e-02 7.215e-02 2.610e+02 0.227 0.820788
## ---
## 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.504
## Unadjusted ICC: 0.396