# TUGAS BIOSTATISTIKA INTERMEDIATE
# Nama : Nisa Nur Fitriana
# Nim : 2611018026
# REPEATED MEASURE ANALYSIS - PERCEIVED SUSCEPTIBILITY, 3 TREATMENT x 3 WAKTU
# Sumber: Bahrami A, Rahimzadeh M, Safari-Moradabadi A. (2025)
# Preventive Medicine Reports, 57, 103198. DOI: 10.1016/j.pmedr.2025.103198
# PMID: 40838176; PMCID: PMC12362373.
#
# CATATAN: Data individual pada CSV adalah DATA SIMULASI untuk latihan statistik.
# Mean dan SD setiap sel dikalibrasi pada Tabel 3 artikel asli untuk konstruk
# perceived susceptibility. Statistik inferensial di bawah berasal dari data simulasi,
# bukan reproduksi output statistik artikel asli.
# 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_PAP_HBM_Data.csv")
dat_wide$id <- factor(dat_wide$id)
dat_wide$treatment <- factor(dat_wide$treatment,
levels=c("Control","In-Person","Mobile-Based"))
dat_long <- dat_wide |>
pivot_longer(c(susceptibility_baseline,
susceptibility_immediate,
susceptibility_3months),
names_to="time",values_to="susceptibility") |>
mutate(time=factor(time,
levels=c("susceptibility_baseline","susceptibility_immediate","susceptibility_3months"),
labels=c("Baseline","Immediate post","3 months")))
head(dat_wide)
## id treatment susceptibility_baseline susceptibility_immediate
## 1 C01 Control 25.870 28.213
## 2 C02 Control 22.517 28.209
## 3 C03 Control 18.664 17.195
## 4 C04 Control 17.766 19.127
## 5 C05 Control 26.151 23.855
## 6 C06 Control 29.986 25.205
## susceptibility_3months
## 1 23.668
## 2 25.822
## 3 20.268
## 4 20.381
## 5 30.006
## 6 21.566
str(dat_long)
## tibble [405 × 4] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 135 levels "C01","C02","C03",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ treatment : Factor w/ 3 levels "Control","In-Person",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ time : Factor w/ 3 levels "Baseline","Immediate post",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ susceptibility: num [1:405] 25.9 28.2 23.7 22.5 28.2 ...
# 2. EKSPLORASI DATA
dat_long |> group_by(treatment,time) |>
get_summary_stats(susceptibility,type="mean_sd")
## # A tibble: 9 × 6
## treatment time variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Control Baseline susceptibility 45 22.5 5.19
## 2 Control Immediate post susceptibility 45 22.8 5.15
## 3 Control 3 months susceptibility 45 22.8 5.28
## 4 In-Person Baseline susceptibility 45 22.5 5.25
## 5 In-Person Immediate post susceptibility 45 26.4 4.19
## 6 In-Person 3 months susceptibility 45 25.2 4.6
## 7 Mobile-Based Baseline susceptibility 45 24.3 3.96
## 8 Mobile-Based Immediate post susceptibility 45 27.5 3.62
## 9 Mobile-Based 3 months susceptibility 45 25.6 3.76
ggplot(dat_long,aes(time,susceptibility,colour=treatment,group=treatment)) +
stat_summary(fun=mean,geom="line",linewidth=1) +
stat_summary(fun=mean,geom="point",size=2.5) +
labs(y="Perceived susceptibility score",x="Time",colour="Treatment")
ggplot(dat_long,aes(time,susceptibility,group=id,colour=treatment)) +
geom_line(alpha=.20) + facet_wrap(~treatment) +
stat_summary(aes(group=1),fun=mean,geom="line",linewidth=1.2) +
labs(y="Perceived susceptibility score",x="Time")
# 3. REPEATED MEASURE ANOVA SATU ARAH - Mobile-Based
d1 <- droplevels(filter(dat_long,treatment=="Mobile-Based"))
d1 |> group_by(time) |> identify_outliers(susceptibility)
## # A tibble: 1 × 6
## time id treatment susceptibility is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <lgl> <lgl>
## 1 3 months M28 Mobile-Based 15.2 TRUE FALSE
d1 |> group_by(time) |> shapiro_test(susceptibility)
## # A tibble: 3 × 4
## time variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline susceptibility 0.962 0.146
## 2 Immediate post susceptibility 0.975 0.435
## 3 3 months susceptibility 0.967 0.234
ggpubr::ggqqplot(d1,"susceptibility",facet.by="time")
aov1_rs <- anova_test(data=d1,dv=susceptibility,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 88 23.687 5.88e-09 * 0.35
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 time 0.994 0.872
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 time 0.994 1.99, 87.44 6.49e-09 * 1.041 2.08, 91.57 5.88e-09
## 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 88 23.687 5.88e-09 * 0.35
aov1 <- aov_ez(id="id",dv="susceptibility",data=d1,within="time",
anova_table=list(es=c("ges","pes"),correction="GG"))
aov1
## Anova Table (Type 3 tests)
##
## Response: susceptibility
## Effect df MSE F ges pes p.value
## 1 time 1.99, 87.44 4.93 23.69 *** .109 .350 <.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) 89815 1 1457.71 44 2710.994 < 2.2e-16 ***
## time 232 2 430.93 88 23.687 5.883e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.99363 0.87166
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.99367 6.493e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 1.040525 5.882585e-09
em1 <- emmeans(aov1,~time)
pairs(em1,adjust="holm")
## contrast estimate SE df t.ratio p.value
## Baseline - Immediate.post -3.20 0.478 44 -6.697 <0.0001
## Baseline - X3.months -1.37 0.473 44 -2.894 0.0059
## Immediate.post - X3.months 1.83 0.448 44 4.088 0.0004
##
## P value adjustment: holm method for 3 tests
friedman_test(d1,susceptibility ~ time | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 susceptibility 45 28.3 2 0.000000712 Friedman test
d1 |> wilcox_test(susceptibility ~ time,paired=TRUE,p.adjust.method="holm")
## # 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 susceptibili… Basel… Immed… 45 45 88 9.62e-8 2.89e-7 ****
## 2 susceptibili… Basel… 3 mon… 45 45 286 8.21e-3 8.21e-3 **
## 3 susceptibili… Immed… 3 mon… 45 45 839 1.63e-4 3.27e-4 ***
# 4. MIXED DESIGN ANOVA: treatment (between) x time (within)
dat_long |> group_by(treatment,time) |> shapiro_test(susceptibility)
## # A tibble: 9 × 5
## treatment time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Control Baseline susceptibility 0.962 0.152
## 2 Control Immediate post susceptibility 0.985 0.813
## 3 Control 3 months susceptibility 0.978 0.531
## 4 In-Person Baseline susceptibility 0.974 0.407
## 5 In-Person Immediate post susceptibility 0.968 0.253
## 6 In-Person 3 months susceptibility 0.971 0.326
## 7 Mobile-Based Baseline susceptibility 0.962 0.146
## 8 Mobile-Based Immediate post susceptibility 0.975 0.435
## 9 Mobile-Based 3 months susceptibility 0.967 0.234
dat_long |> group_by(time) |> levene_test(susceptibility ~ treatment)
## # A tibble: 3 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 132 0.665 0.516
## 2 Immediate post 2 132 1.61 0.203
## 3 3 months 2 132 2.02 0.136
box_m(dat_wide[,c("susceptibility_baseline","susceptibility_immediate","susceptibility_3months")],
dat_wide$treatment)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 12.2 0.432 12 Box's M-test for Homogeneity of Covariance Matric…
aov2 <- aov_ez(id="id",dv="susceptibility",data=dat_long,
between="treatment",within="time",
anova_table=list(es=c("ges","pes"),correction="GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: susceptibility
## Effect df MSE F ges pes p.value
## 1 treatment 2, 132 53.71 6.23 ** .074 .086 .003
## 2 time 1.95, 257.30 5.03 42.59 *** .047 .244 <.001
## 3 treatment:time 3.90, 257.30 5.03 9.38 *** .021 .124 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2) # termasuk Mauchly dan koreksi GG/HF
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 240989 1 7089.3 132 4487.1076 < 2.2e-16 ***
## treatment 669 2 7089.3 132 6.2255 0.002608 **
## time 418 2 1294.5 264 42.5906 < 2.2e-16 ***
## treatment:time 184 4 1294.5 264 9.3762 4.226e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.97398 0.17782
## treatment:time 0.97398 0.17782
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.97464 < 2.2e-16 ***
## treatment:time 0.97464 5.688e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.9890518 1.333939e-16
## treatment:time 0.9890518 4.804119e-07
eta_squared(aov2,partial=TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## treatment | 0.09 | [0.02, 1.00]
## time | 0.24 | [0.17, 1.00]
## treatment:time | 0.12 | [0.06, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Simple effects
joint_tests(aov2,by="treatment")
## treatment = Control:
## model term df1 df2 F.ratio p.value
## time 2 132 0.245 0.7834
##
## treatment = In-Person:
## model term df1 df2 F.ratio p.value
## time 2 132 34.317 <0.0001
##
## treatment = Mobile-Based:
## model term df1 df2 F.ratio p.value
## time 2 132 23.753 <0.0001
joint_tests(aov2,by="time")
## time = Baseline:
## model term df1 df2 F.ratio p.value
## treatment 2 132 2.043 0.1337
##
## time = Immediate.post:
## model term df1 df2 F.ratio p.value
## treatment 2 132 14.454 <0.0001
##
## time = X3.months:
## model term df1 df2 F.ratio p.value
## treatment 2 132 4.891 0.0089
# Post hoc waktu pada tiap treatment
em2 <- emmeans(aov2,~time|treatment)
contrast(em2,"trt.vs.ctrl",ref=1,adjust="holm")
## treatment = Control:
## contrast estimate SE df t.ratio p.value
## Immediate.post - Baseline 0.27 0.478 132 0.565 1.0000
## X3.months - Baseline 0.33 0.492 132 0.671 1.0000
##
## treatment = In-Person:
## contrast estimate SE df t.ratio p.value
## Immediate.post - Baseline 3.95 0.478 132 8.265 <0.0001
## X3.months - Baseline 2.70 0.492 132 5.491 <0.0001
##
## treatment = Mobile-Based:
## contrast estimate SE df t.ratio p.value
## Immediate.post - Baseline 3.20 0.478 132 6.696 <0.0001
## X3.months - Baseline 1.37 0.492 132 2.786 0.0061
##
## P value adjustment: holm method for 2 tests
# Perbandingan treatment pada tiap waktu
em2b <- emmeans(aov2,~treatment|time)
pairs(em2b,adjust="tukey")
## time = Baseline:
## contrast estimate SE df t.ratio p.value
## Control - (In-Person) 0.00996 1.020 132 0.010 0.9999
## Control - (Mobile-Based) -1.78000 1.020 132 -1.746 0.1922
## (In-Person) - (Mobile-Based) -1.78996 1.020 132 -1.755 0.1888
##
## time = Immediate.post:
## contrast estimate SE df t.ratio p.value
## Control - (In-Person) -3.67000 0.920 132 -3.987 0.0003
## Control - (Mobile-Based) -4.71000 0.920 132 -5.117 <0.0001
## (In-Person) - (Mobile-Based) -1.04000 0.920 132 -1.130 0.4974
##
## time = X3.months:
## contrast estimate SE df t.ratio p.value
## Control - (In-Person) -2.35996 0.967 132 -2.439 0.0421
## Control - (Mobile-Based) -2.81998 0.967 132 -2.915 0.0116
## (In-Person) - (Mobile-Based) -0.46002 0.967 132 -0.476 0.8830
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# Kontras perubahan: immediate-baseline dan 3months-baseline antarkelompok
em_full <- emmeans(aov2,~time*treatment)
contrast(em_full,
interaction=list(time=list("Immediate-Baseline"=c(-1,1,0)),treatment="pairwise"),
adjust="holm")
## time_custom treatment_pairwise estimate SE df t.ratio
## Immediate-Baseline Control - (In-Person) -3.68 0.676 132 -5.445
## Immediate-Baseline Control - (Mobile-Based) -2.93 0.676 132 -4.335
## Immediate-Baseline (In-Person) - (Mobile-Based) 0.75 0.676 132 1.110
## p.value
## <0.0001
## <0.0001
## 0.2692
##
## P value adjustment: holm method for 3 tests
contrast(em_full,
interaction=list(time=list("3months-Baseline"=c(-1,0,1)),treatment="pairwise"),
adjust="holm")
## time_custom treatment_pairwise estimate SE df t.ratio
## 3months-Baseline Control - (In-Person) -2.37 0.695 132 -3.408
## 3months-Baseline Control - (Mobile-Based) -1.04 0.695 132 -1.496
## 3months-Baseline (In-Person) - (Mobile-Based) 1.33 0.695 132 1.913
## p.value
## 0.0026
## 0.1371
## 0.1159
##
## P value adjustment: holm method for 3 tests
# 5. PEMBANDING: LINEAR MIXED MODEL
lmm <- lmer(susceptibility ~ 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 61.05 30.526 2 132 6.2255 0.002608 **
## time 417.69 208.843 2 264 42.5906 < 2.2e-16 ***
## treatment:time 183.90 45.976 4 264 9.3762 4.226e-07 ***
## ---
## 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: susceptibility ~ treatment * time + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 2116.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02243 -0.54999 -0.03114 0.55121 2.40500
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 16.268 4.033
## Residual 4.903 2.214
## Number of obs: 405, groups: id, 135
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 24.3933 0.3642 132.0000 66.986 < 2e-16 ***
## treatment1 -1.7033 0.5150 132.0000 -3.307 0.00121 **
## treatment2 0.3033 0.5150 132.0000 0.589 0.55686
## time1 -1.3133 0.1556 264.0000 -8.440 2.13e-15 ***
## time2 1.1600 0.1556 264.0000 7.455 1.29e-12 ***
## treatment1:time1 1.1133 0.2201 264.0000 5.059 7.90e-07 ***
## treatment2:time1 -0.9033 0.2201 264.0000 -4.105 5.40e-05 ***
## treatment1:time2 -1.0900 0.2201 264.0000 -4.953 1.31e-06 ***
## treatment2:time2 0.5733 0.2201 264.0000 2.605 0.00970 **
## ---
## 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.768
## Unadjusted ICC: 0.669
# 6. SIMPAN RINGKASAN
write.csv(dat_long,"FINAL_PAP_HBM_Data_LONG.csv",row.names=FALSE)