# ============================================================
# TUGAS BIOSTATISTIKA INTERMEDIATE
# Nama : Zulkaeni Septia Rini
# Nim : 2611018051
# FINAL - REPEATED MEASURES ANALYSIS DENGAN R
# Data simulasi berdasarkan Chikama et al. (2026)
# Outcome: Athens Insomnia Scale (AIS)
# Desain utama: 3 treatment x 3 kali pengukuran
# ============================================================
# Jika paket belum tersedia, jalankan satu kali:
# install.packages(c("dplyr","tidyr","ggplot2","afex","emmeans",
# "rstatix","car","effectsize","lme4","lmerTest",
# "performance","ggpubr"))
required_packages <- c(
"dplyr","tidyr","ggplot2","afex","emmeans","rstatix","car",
"effectsize","lme4","lmerTest","performance","ggpubr"
)
missing_packages <- required_packages[!vapply(required_packages, requireNamespace,
quietly = TRUE, FUN.VALUE = logical(1))]
## Registered S3 method overwritten by 'lme4':
## method from
## na.action.merMod car
if (length(missing_packages) > 0) {
stop(
"Paket berikut belum terpasang: ", paste(missing_packages, collapse = ", "),
"\nJalankan install.packages() seperti pada bagian atas script, lalu ulangi analisis."
)
}
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 DAN VALIDASI STRUKTUR
data_file <- "FINAL_SLEEP_AIS_Data.csv"
if (!file.exists(data_file)) {
stop(
"File data tidak ditemukan: ", data_file,
"\nPastikan file .R dan FINAL_SLEEP_AIS_Data.csv berada pada folder kerja yang sama."
)
}
dat_wide <- read.csv(data_file, stringsAsFactors = FALSE)
required_columns <- c("id", "treatment", "AIS_baseline", "AIS_3m", "AIS_6m")
if (!all(required_columns %in% names(dat_wide))) {
stop("Kolom data tidak lengkap. Kolom wajib: ", paste(required_columns, collapse = ", "))
}
if (anyNA(dat_wide[, required_columns])) {
stop("Terdapat missing value pada variabel utama. Periksa dataset sebelum analisis.")
}
if (anyDuplicated(dat_wide$id) > 0) {
stop("ID peserta tidak unik pada data format wide.")
}
# Urutan kelompok dibuat konsisten dengan laporan: A, B, C.
dat_wide$treatment <- factor(
dat_wide$treatment,
levels = c(
"A - Guidance 6 months",
"B - Guidance 3 months",
"C - Report only"
)
)
dat_wide$id <- factor(dat_wide$id)
if (nlevels(dat_wide$treatment) < 3) {
stop("Desain belum memenuhi minimal 3 treatment.")
}
dat_long <- dat_wide |>
pivot_longer(
cols = c(AIS_baseline, AIS_3m, AIS_6m),
names_to = "time",
values_to = "AIS"
) |>
mutate(
time = factor(
time,
levels = c("AIS_baseline", "AIS_3m", "AIS_6m"),
labels = c("Baseline", "3 months", "6 months")
)
)
if (nlevels(dat_long$time) < 3) {
stop("Desain belum memenuhi minimal 3 kali pengukuran.")
}
cat("\n=== VALIDASI DESAIN ===\n")
##
## === VALIDASI DESAIN ===
cat("Jumlah peserta :", n_distinct(dat_wide$id), "\n")
## Jumlah peserta : 90
cat("Jumlah treatment:", nlevels(dat_wide$treatment), "\n")
## Jumlah treatment: 3
cat("Jumlah waktu :", nlevels(dat_long$time), "\n")
## Jumlah waktu : 3
print(table(dat_wide$treatment))
##
## A - Guidance 6 months B - Guidance 3 months C - Report only
## 30 30 30
# 2. EKSPLORASI DATA
cat("\n=== STATISTIK DESKRIPTIF ===\n")
##
## === STATISTIK DESKRIPTIF ===
descriptive <- dat_long |>
group_by(treatment, time) |>
get_summary_stats(AIS, type = "mean_sd")
print(descriptive)
## # A tibble: 9 × 6
## treatment time variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 A - Guidance 6 months Baseline AIS 30 6.67 2.25
## 2 A - Guidance 6 months 3 months AIS 30 4.93 2.21
## 3 A - Guidance 6 months 6 months AIS 30 5.13 2.28
## 4 B - Guidance 3 months Baseline AIS 30 6.27 2.02
## 5 B - Guidance 3 months 3 months AIS 30 4 2.32
## 6 B - Guidance 3 months 6 months AIS 30 4.87 2.14
## 7 C - Report only Baseline AIS 30 6.23 2.05
## 8 C - Report only 3 months AIS 30 5.93 2.16
## 9 C - Report only 6 months AIS 30 6.2 2.31
p_profile <- ggplot(dat_long, aes(x = time, y = AIS, 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 = 0.08) +
labs(
title = "Mean AIS profile (95% CI)",
x = "Time",
y = "Athens Insomnia Scale (AIS)",
colour = "Treatment"
)
print(p_profile)
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
# 3. ANALISIS TAMBAHAN: REPEATED-MEASURES ANOVA PADA GROUP A
# Bagian ini bukan analisis utama desain 3 x 3.
cat("\n=== ANALISIS TAMBAHAN: GROUP A ===\n")
##
## === ANALISIS TAMBAHAN: GROUP A ===
d1 <- droplevels(filter(dat_long, treatment == "A - Guidance 6 months"))
cat("\nNormalitas per waktu - Group A:\n")
##
## Normalitas per waktu - Group A:
print(d1 |> group_by(time) |> shapiro_test(AIS))
## # A tibble: 3 × 4
## time variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline AIS 0.939 0.0847
## 2 3 months AIS 0.960 0.303
## 3 6 months AIS 0.953 0.209
aov1 <- aov_ez(
id = "id",
dv = "AIS",
data = d1,
within = "time",
anova_table = list(es = c("ges", "pes"), correction = "GG")
)
cat("\nRepeated-measures ANOVA - Group A:\n")
##
## Repeated-measures ANOVA - Group A:
print(aov1)
## Anova Table (Type 3 tests)
##
## Response: AIS
## Effect df MSE F ges pes p.value
## 1 time 1.84, 53.47 1.38 21.13 *** .109 .422 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
print(summary(aov1)) # termasuk Mauchly/sphericity bila tersedia
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 2800.04 1 365.96 29 221.888 4.019e-15 ***
## time 53.96 2 74.04 58 21.132 1.277e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.91519 0.28917
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.92182 3.458e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.9815234 1.615764e-07
cat("\nPost hoc waktu - Group A (Holm):\n")
##
## Post hoc waktu - Group A (Holm):
print(emmeans(aov1, ~ time) |> pairs(adjust = "holm"))
## contrast estimate SE df t.ratio p.value
## Baseline - X3.months 1.73 0.318 29 5.454 <0.0001
## Baseline - X6.months 1.53 0.306 29 5.011 <0.0001
## X3.months - X6.months -0.20 0.246 29 -0.812 0.4235
##
## P value adjustment: holm method for 3 tests
cat("\nAlternatif nonparametrik - Friedman:\n")
##
## Alternatif nonparametrik - Friedman:
print(friedman_test(d1, AIS ~ time | id))
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 AIS 30 20.4 2 0.0000377 Friedman test
# 4. ANALISIS UTAMA: MIXED REPEATED-MEASURES ANOVA 3 x 3
# Between-subject factor : treatment (3 kelompok)
# Within-subject factor : time (3 pengukuran)
cat("\n=== ANALISIS UTAMA: MIXED ANOVA 3 x 3 ===\n")
##
## === ANALISIS UTAMA: MIXED ANOVA 3 x 3 ===
cat("\nNormalitas per treatment x waktu:\n")
##
## Normalitas per treatment x waktu:
print(dat_long |> group_by(treatment, time) |> shapiro_test(AIS))
## # A tibble: 9 × 5
## treatment time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 A - Guidance 6 months Baseline AIS 0.939 0.0847
## 2 A - Guidance 6 months 3 months AIS 0.960 0.303
## 3 A - Guidance 6 months 6 months AIS 0.953 0.209
## 4 B - Guidance 3 months Baseline AIS 0.957 0.254
## 5 B - Guidance 3 months 3 months AIS 0.941 0.0996
## 6 B - Guidance 3 months 6 months AIS 0.956 0.249
## 7 C - Report only Baseline AIS 0.935 0.0650
## 8 C - Report only 3 months AIS 0.951 0.175
## 9 C - Report only 6 months AIS 0.978 0.767
cat("\nHomogenitas varians (Levene) pada setiap waktu:\n")
##
## Homogenitas varians (Levene) pada setiap waktu:
print(dat_long |> group_by(time) |> levene_test(AIS ~ treatment))
## # A tibble: 3 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 87 0.0861 0.918
## 2 3 months 2 87 0.0786 0.925
## 3 6 months 2 87 0.0246 0.976
cat("\nHomogenitas matriks kovarians (Box's M):\n")
##
## Homogenitas matriks kovarians (Box's M):
print(box_m(
dat_wide[, c("AIS_baseline", "AIS_3m", "AIS_6m")],
dat_wide$treatment
))
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 7.12 0.850 12 Box's M-test for Homogeneity of Covariance Matric…
aov2 <- aov_ez(
id = "id",
dv = "AIS",
data = dat_long,
between = "treatment",
within = "time",
anova_table = list(es = c("ges", "pes"), correction = "GG")
)
cat("\nMixed repeated-measures ANOVA:\n")
##
## Mixed repeated-measures ANOVA:
print(aov2)
## Anova Table (Type 3 tests)
##
## Response: AIS
## Effect df MSE F ges pes p.value
## 1 treatment 2, 87 12.42 2.10 .040 .046 .128
## 2 time 1.94, 168.95 1.06 47.20 *** .071 .352 <.001
## 3 treatment:time 3.88, 168.95 1.06 9.05 *** .029 .172 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
print(summary(aov2))
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 8411.3 1 1080.8 87 677.0957 < 2.2e-16 ***
## treatment 52.3 2 1080.8 87 2.1040 0.1281
## time 96.9 2 178.6 174 47.2003 < 2.2e-16 ***
## treatment:time 37.2 4 178.6 174 9.0533 1.149e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.97008 0.27084
## treatment:time 0.97008 0.27084
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.97095 < 2.2e-16 ***
## treatment:time 0.97095 1.572e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.9927741 5.361951e-17
## treatment:time 0.9927741 1.241850e-06
cat("\nPartial eta-squared:\n")
##
## Partial eta-squared:
print(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.35 | [0.26, 1.00]
## treatment:time | 0.17 | [0.08, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Simple effects
cat("\nSimple effect waktu pada setiap treatment:\n")
##
## Simple effect waktu pada setiap treatment:
print(joint_tests(aov2, by = "treatment"))
## treatment = A - Guidance 6 months:
## model term df1 df2 F.ratio p.value
## time 2 87 23.675 <0.0001
##
## treatment = B - Guidance 3 months:
## model term df1 df2 F.ratio p.value
## time 2 87 37.715 <0.0001
##
## treatment = C - Report only:
## model term df1 df2 F.ratio p.value
## time 2 87 0.922 0.4017
cat("\nSimple effect treatment pada setiap waktu:\n")
##
## Simple effect treatment pada setiap waktu:
print(joint_tests(aov2, by = "time"))
## time = Baseline:
## model term df1 df2 F.ratio p.value
## treatment 2 87 0.393 0.6760
##
## time = X3.months:
## model term df1 df2 F.ratio p.value
## treatment 2 87 5.625 0.0050
##
## time = X6.months:
## model term df1 df2 F.ratio p.value
## treatment 2 87 2.955 0.0574
# Post hoc / contrasts
cat("\nPerbandingan waktu di setiap treatment (baseline sebagai referensi; Holm):\n")
##
## Perbandingan waktu di setiap treatment (baseline sebagai referensi; Holm):
em_time <- emmeans(aov2, ~ time | treatment)
print(contrast(em_time, "trt.vs.ctrl", ref = 1, adjust = "holm"))
## treatment = A - Guidance 6 months:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline -1.7333 0.261 87 -6.637 <0.0001
## X6.months - Baseline -1.5333 0.281 87 -5.463 <0.0001
##
## treatment = B - Guidance 3 months:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline -2.2667 0.261 87 -8.679 <0.0001
## X6.months - Baseline -1.4000 0.281 87 -4.988 <0.0001
##
## treatment = C - Report only:
## contrast estimate SE df t.ratio p.value
## X3.months - Baseline -0.3000 0.261 87 -1.149 0.5077
## X6.months - Baseline -0.0333 0.281 87 -0.119 0.9057
##
## P value adjustment: holm method for 2 tests
cat("\nPerbandingan treatment pada setiap waktu (Tukey):\n")
##
## Perbandingan treatment pada setiap waktu (Tukey):
em_group <- emmeans(aov2, ~ treatment | time)
print(pairs(em_group, adjust = "tukey"))
## time = Baseline:
## contrast estimate SE df t.ratio
## (A - Guidance 6 months) - (B - Guidance 3 months) 0.4000 0.544 87 0.736
## (A - Guidance 6 months) - (C - Report only) 0.4333 0.544 87 0.797
## (B - Guidance 3 months) - (C - Report only) 0.0333 0.544 87 0.061
## p.value
## 0.7431
## 0.7059
## 0.9979
##
## time = X3.months:
## contrast estimate SE df t.ratio
## (A - Guidance 6 months) - (B - Guidance 3 months) 0.9333 0.577 87 1.619
## (A - Guidance 6 months) - (C - Report only) -1.0000 0.577 87 -1.735
## (B - Guidance 3 months) - (C - Report only) -1.9333 0.577 87 -3.354
## p.value
## 0.2431
## 0.1982
## 0.0034
##
## time = X6.months:
## contrast estimate SE df t.ratio
## (A - Guidance 6 months) - (B - Guidance 3 months) 0.2667 0.580 87 0.459
## (A - Guidance 6 months) - (C - Report only) -1.0667 0.580 87 -1.838
## (B - Guidance 3 months) - (C - Report only) -1.3333 0.580 87 -2.297
## p.value
## 0.8903
## 0.1635
## 0.0615
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# Perubahan baseline -> 3 bulan dan baseline -> 6 bulan dibandingkan antarkelompok
cat("\nKontras difference-in-change antartreatment (Holm):\n")
##
## Kontras difference-in-change antartreatment (Holm):
em_full <- emmeans(aov2, ~ time * treatment)
change_contrasts <- contrast(
em_full,
interaction = list(
time = list(
"3m-base" = c(-1, 1, 0),
"6m-base" = c(-1, 0, 1)
),
treatment = "pairwise"
),
adjust = "holm"
)
print(change_contrasts)
## time_custom treatment_pairwise estimate SE
## 3m-base (A - Guidance 6 months) - (B - Guidance 3 months) 0.533 0.369
## 6m-base (A - Guidance 6 months) - (B - Guidance 3 months) -0.133 0.397
## 3m-base (A - Guidance 6 months) - (C - Report only) -1.433 0.369
## 6m-base (A - Guidance 6 months) - (C - Report only) -1.500 0.397
## 3m-base (B - Guidance 3 months) - (C - Report only) -1.967 0.369
## 6m-base (B - Guidance 3 months) - (C - Report only) -1.367 0.397
## df t.ratio p.value
## 87 1.444 0.3047
## 87 -0.336 0.7378
## 87 -3.881 0.0010
## 87 -3.779 0.0012
## 87 -5.325 <0.0001
## 87 -3.443 0.0027
##
## P value adjustment: holm method for 6 tests
# 5. LINEAR MIXED MODEL (PEMBANDING)
cat("\n=== PEMBANDING: LINEAR MIXED MODEL ===\n")
##
## === PEMBANDING: LINEAR MIXED MODEL ===
lmm <- lmer(AIS ~ treatment * time + (1 | id), data = dat_long, REML = TRUE)
print(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 4.319 2.160 2 87 2.1040 0.1281
## time 96.896 48.448 2 174 47.2003 < 2.2e-16 ***
## treatment:time 37.170 9.293 4 174 9.0533 1.149e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(lmm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: AIS ~ treatment * time + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 1008.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.13333 -0.51716 0.01363 0.51692 2.22260
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 3.799 1.949
## Residual 1.026 1.013
## Number of obs: 270, groups: id, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.581481 0.214499 86.999996 26.021 < 2e-16 ***
## treatment1 -0.003704 0.303347 86.999996 -0.012 0.990286
## treatment2 -0.537037 0.303347 86.999996 -1.770 0.080168 .
## time1 0.807407 0.087197 174.000001 9.260 < 2e-16 ***
## time2 -0.625926 0.087197 174.000001 -7.178 1.96e-11 ***
## treatment1:time1 0.281481 0.123315 174.000001 2.283 0.023662 *
## treatment2:time1 0.414815 0.123315 174.000001 3.364 0.000945 ***
## treatment1:time2 -0.018519 0.123315 174.000001 -0.150 0.880802
## treatment2:time2 -0.418519 0.123315 174.000001 -3.394 0.000853 ***
## ---
## 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
print(performance::icc(lmm))
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.787
## Unadjusted ICC: 0.688
cat("\nAnalisis selesai. Analisis utama tugas adalah Mixed Repeated-Measures ANOVA 3 treatment x 3 pengukuran.\n")
##
## Analisis selesai. Analisis utama tugas adalah Mixed Repeated-Measures ANOVA 3 treatment x 3 pengukuran.