# ================================================================
# BIOSTATISTIKA INTERMEDIATE
# REPEATED MEASURE ANALYSIS DENGAN R
# NAMA : IMAM FATHONI
# NIM : 261018020
# Data simulasi berdasarkan Sun et al. (2026), Journal of Medical Internet Research
# 3 treatment x 3 repeated measurements
# Outcome: MVPA (minutes/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_BLEND_MVPA_Data.csv")
dat_wide$treatment <- factor(dat_wide$treatment,
levels=c("Control","Web-based","Blended"))
dat_wide$id <- factor(dat_wide$id)
dat_long <- dat_wide |>
pivot_longer(c(mvpa_baseline,mvpa_12w,mvpa_24w),
names_to="time",values_to="mvpa") |>
mutate(time=factor(time,
levels=c("mvpa_baseline","mvpa_12w","mvpa_24w"),
labels=c("Baseline","12 weeks","24 weeks")))
# 2. EKSPLORASI
desk <- dat_long |> group_by(treatment,time) |>
get_summary_stats(mvpa,type="mean_sd")
print(desk)
## # A tibble: 9 × 6
## treatment time variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Control Baseline mvpa 30 159. 35.5
## 2 Control 12 weeks mvpa 30 184. 39.1
## 3 Control 24 weeks mvpa 30 195. 41.6
## 4 Web-based Baseline mvpa 30 163. 45.5
## 5 Web-based 12 weeks mvpa 30 186. 42.3
## 6 Web-based 24 weeks mvpa 30 198. 39.1
## 7 Blended Baseline mvpa 30 163. 32.0
## 8 Blended 12 weeks mvpa 30 245. 33.8
## 9 Blended 24 weeks mvpa 30 268 33.5
ggplot(dat_long,aes(time,mvpa,colour=treatment,group=treatment))+
stat_summary(fun=mean,geom="line",linewidth=1)+
stat_summary(fun=mean,geom="point",size=2.5)+
labs(y="MVPA (min/week)",x="Time")
ggplot(dat_long,aes(time,mvpa,group=id,colour=treatment))+
geom_line(alpha=.18)+facet_wrap(~treatment)+
stat_summary(aes(group=treatment),fun=mean,geom="line",linewidth=1.2)
# 3. REPEATED MEASURE ANOVA SATU ARAH: Blended
d1 <- droplevels(filter(dat_long,treatment=="Blended"))
d1 |> group_by(time) |> shapiro_test(mvpa)
## # A tibble: 3 × 4
## time variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline mvpa 0.974 0.648
## 2 12 weeks mvpa 0.950 0.170
## 3 24 weeks mvpa 0.956 0.240
aov1_rs <- anova_test(data=d1,dv=mvpa,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 2 58 334.176 1.47e-32 * 0.92
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 time 0.955 0.524
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 time 0.957 1.91, 55.5 2.94e-31 * 1.023 2.05, 59.32 1.47e-32
## 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 2 58 334.176 1.47e-32 * 0.92
aov1 <- aov_ez(id="id",dv="mvpa",data=d1,within="time",
anova_table=list(es=c("ges","pes"),correction="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) 4569354 1 79662 29 1663.43 < 2.2e-16 ***
## time 181912 2 15786 58 334.18 < 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.95486 0.52381
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.95681 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 1.022815 1.466062e-32
em1 <- emmeans(aov1,~time)
pairs(em1,adjust="holm")
## contrast estimate SE df t.ratio p.value
## Baseline - X12.weeks -81.6 4.09 29 -19.945 <0.0001
## Baseline - X24.weeks -104.8 3.97 29 -26.431 <0.0001
## X12.weeks - X24.weeks -23.2 4.69 29 -4.951 <0.0001
##
## P value adjustment: holm method for 3 tests
friedman_test(d1,mvpa~time|id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 mvpa 30 49.3 2 2.00e-11 Friedman test
# 4. MIXED DESIGN ANOVA: treatment x time
dat_long |> group_by(treatment,time) |> shapiro_test(mvpa)
## # A tibble: 9 × 5
## treatment time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Control Baseline mvpa 0.972 0.583
## 2 Control 12 weeks mvpa 0.975 0.694
## 3 Control 24 weeks mvpa 0.967 0.459
## 4 Web-based Baseline mvpa 0.966 0.425
## 5 Web-based 12 weeks mvpa 0.959 0.298
## 6 Web-based 24 weeks mvpa 0.977 0.734
## 7 Blended Baseline mvpa 0.974 0.648
## 8 Blended 12 weeks mvpa 0.950 0.170
## 9 Blended 24 weeks mvpa 0.956 0.240
dat_long |> group_by(time) |> levene_test(mvpa~treatment)
## # A tibble: 3 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 87 2.29 0.107
## 2 12 weeks 2 87 0.993 0.375
## 3 24 weeks 2 87 0.795 0.455
box_m(dat_wide[,c("mvpa_baseline","mvpa_12w","mvpa_24w")],dat_wide$treatment)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 8.87 0.714 12 Box's M-test for Homogeneity of Covariance Matric…
aov2 <- aov_ez(id="id",dv="mvpa",data=dat_long,
between="treatment",within="time",
anova_table=list(es=c("ges","pes"),correction="GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: mvpa
## Effect df MSE F ges pes p.value
## 1 treatment 2, 87 3733.77 16.03 *** .238 .269 <.001
## 2 time 1.97, 171.72 336.52 246.11 *** .299 .739 <.001
## 3 treatment:time 3.95, 171.72 336.52 42.35 *** .128 .493 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2) # Mauchly + GG/HF
## 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) 10329789 1 324838 87 2766.586 < 2.2e-16 ***
## treatment 119738 2 324838 87 16.035 1.18e-06 ***
## time 163470 2 57788 174 246.106 < 2.2e-16 ***
## treatment:time 56257 4 57788 174 42.347 < 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.98672 0.56267
## treatment:time 0.98672 0.56267
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.98689 < 2.2e-16 ***
## treatment:time 0.98689 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 1.009645 1.877247e-51
## treatment:time 1.009645 9.054597e-25
eta_squared(aov2,partial=TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## treatment | 0.27 | [0.14, 1.00]
## time | 0.74 | [0.69, 1.00]
## treatment:time | 0.49 | [0.40, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# simple effects dan post hoc
joint_tests(aov2,by="treatment")
## treatment = Control:
## model term df1 df2 F.ratio p.value
## time 2 87 32.011 <0.0001
##
## treatment = Web-based:
## model term df1 df2 F.ratio p.value
## time 2 87 29.185 <0.0001
##
## treatment = Blended:
## model term df1 df2 F.ratio p.value
## time 2 87 300.114 <0.0001
joint_tests(aov2,by="time")
## time = Baseline:
## model term df1 df2 F.ratio p.value
## treatment 2 87 0.118 0.8888
##
## time = X12.weeks:
## model term df1 df2 F.ratio p.value
## treatment 2 87 24.304 <0.0001
##
## time = X24.weeks:
## model term df1 df2 F.ratio p.value
## treatment 2 87 35.392 <0.0001
em2 <- emmeans(aov2,~time|treatment)
contrast(em2,"trt.vs.ctrl",ref=1,adjust="holm")
## treatment = Control:
## contrast estimate SE df t.ratio p.value
## X12.weeks - Baseline 24.5 4.49 87 5.461 <0.0001
## X24.weeks - Baseline 35.4 4.66 87 7.587 <0.0001
##
## treatment = Web-based:
## contrast estimate SE df t.ratio p.value
## X12.weeks - Baseline 22.3 4.49 87 4.968 <0.0001
## X24.weeks - Baseline 34.2 4.66 87 7.342 <0.0001
##
## treatment = Blended:
## contrast estimate SE df t.ratio p.value
## X12.weeks - Baseline 81.6 4.49 87 18.180 <0.0001
## X24.weeks - Baseline 104.8 4.66 87 22.486 <0.0001
##
## P value adjustment: holm method for 2 tests
em2b <- emmeans(aov2,~treatment|time)
pairs(em2b,adjust="tukey")
## time = Baseline:
## contrast estimate SE df t.ratio p.value
## Control - (Web-based) -4.281 9.84 87 -0.435 0.9010
## Control - Blended -3.980 9.84 87 -0.405 0.9138
## (Web-based) - Blended 0.301 9.84 87 0.031 0.9995
##
## time = X12.weeks:
## contrast estimate SE df t.ratio p.value
## Control - (Web-based) -2.070 9.95 87 -0.208 0.9764
## Control - Blended -61.090 9.95 87 -6.139 <0.0001
## (Web-based) - Blended -59.021 9.95 87 -5.931 <0.0001
##
## time = X24.weeks:
## contrast estimate SE df t.ratio p.value
## Control - (Web-based) -3.141 9.87 87 -0.318 0.9458
## Control - Blended -73.441 9.87 87 -7.440 <0.0001
## (Web-based) - Blended -70.300 9.87 87 -7.122 <0.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# perubahan baseline -> 24 minggu
change_dat <- dat_wide |>
mutate(change24=mvpa_24w-mvpa_baseline)
anova_test(data=change_dat,dv=change24,between=treatment,effect.size="pes")
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 treatment 2 87 75.23 1.07e-19 * 0.634
change_dat |> tukey_hsd(change24~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 Control Web-based 0 -1.14 -16.9 14.6 9.84e- 1
## 2 treatment Control Blended 0 69.5 53.7 85.2 3.09e-10
## 3 treatment Web-based Blended 0 70.6 54.9 86.3 3.09e-10
## # ℹ 1 more variable: p.adj.signif <chr>
# 5. LINEAR MIXED MODEL
lmm <- lmer(mvpa ~ 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 10650 5325 2 87 16.035 1.18e-06 ***
## time 163470 81735 2 174 246.106 < 2.2e-16 ***
## treatment:time 56257 14064 4 174 42.347 < 2.2e-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: mvpa ~ treatment * time + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 2510.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07924 -0.58062 0.02088 0.58084 1.99674
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1133.9 33.67
## Residual 332.1 18.22
## Number of obs: 270, groups: id, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 195.598 3.719 87.000 52.598 < 2e-16 ***
## treatment1 -16.445 5.259 87.000 -3.127 0.0024 **
## treatment2 -13.281 5.259 87.000 -2.525 0.0134 *
## time1 -33.654 1.568 174.000 -21.457 < 2e-16 ***
## time2 9.166 1.568 174.000 5.844 2.46e-08 ***
## treatment1:time1 13.691 2.218 174.000 6.172 4.60e-09 ***
## treatment2:time1 14.808 2.218 174.000 6.676 3.17e-10 ***
## treatment1:time2 -4.609 2.218 174.000 -2.078 0.0392 *
## treatment2:time2 -5.703 2.218 174.000 -2.571 0.0110 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn1 trtmn2 time1 time2 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.500
## trtmnt1:tm1 0.000 0.000 0.000 0.000 0.000
## trtmnt2:tm1 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.500 0.250
## trtmnt2:tm2 0.000 0.000 0.000 0.000 0.000 0.250 -0.500 -0.500
performance::icc(lmm)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.773
## Unadjusted ICC: 0.416