# ============================================================================
# Nama : Andi Tiawarman
# NIM : 2611018049
# REPEATED MEASURE ANALYSIS DENGAN R
# Topik: Insomnia Severity Index (ISI) pada 3 treatment x 5 pengukuran
# Data : SIMULASI berdasarkan protokol RCT Furihata et al. (2026)
# Sumber: Healthcare. 2026;14(10):1386. DOI:10.3390/healthcare14101386
# ============================================================================
# Jalankan sekali bila paket belum tersedia:
# 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. MEMBACA DATA -------------------------------------------------------------
dat_wide <- read.csv("FINAL_DBBTI_ISI_Data.csv", check.names = FALSE)
dat_wide$id <- factor(dat_wide$id)
dat_wide$treatment <- factor(dat_wide$treatment,
levels = c("Waitlist Control", "Digital BBT-I", "Digital BBT-I + LT"))
# Format panjang: baseline + minggu 1-4 = lima pengukuran utama latihan
# Protokol asli juga memiliki follow-up 3 bulan, tetapi latihan 3x5 ini mengikuti
# model primer artikel: baseline dan minggu 1-4.
dat_long <- dat_wide |>
pivot_longer(cols = starts_with("ISI_"), names_to = "time", values_to = "isi") |>
mutate(time = factor(time,
levels = c("ISI_Baseline","ISI_Week1","ISI_Week2","ISI_Week3","ISI_Week4"),
labels = c("Baseline","Week1","Week2","Week3","Week4")),
week = c(0,1,2,3,4)[match(time, levels(time))])
head(dat_wide)
## id treatment ISI_Baseline ISI_Week1 ISI_Week2 ISI_Week3 ISI_Week4
## 1 C01 Waitlist Control 17 16 19 16 15
## 2 C02 Waitlist Control 16 12 11 13 13
## 3 C03 Waitlist Control 14 16 18 15 12
## 4 C04 Waitlist Control 18 20 18 21 21
## 5 C05 Waitlist Control 16 13 14 17 18
## 6 C06 Waitlist Control 13 14 15 10 14
head(dat_long)
## # A tibble: 6 × 5
## id treatment time isi week
## <fct> <fct> <fct> <int> <dbl>
## 1 C01 Waitlist Control Baseline 17 0
## 2 C01 Waitlist Control Week1 16 1
## 3 C01 Waitlist Control Week2 19 2
## 4 C01 Waitlist Control Week3 16 3
## 5 C01 Waitlist Control Week4 15 4
## 6 C02 Waitlist Control Baseline 16 0
str(dat_long)
## tibble [450 × 5] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 90 levels "B01","B02","B03",..: 31 31 31 31 31 32 32 32 32 32 ...
## $ treatment: Factor w/ 3 levels "Waitlist Control",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ time : Factor w/ 5 levels "Baseline","Week1",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ isi : int [1:450] 17 16 19 16 15 16 12 11 13 13 ...
## $ week : num [1:450] 0 1 2 3 4 0 1 2 3 4 ...
# 2. EKSPLORASI DATA ----------------------------------------------------------
desk <- dat_long |>
group_by(treatment, time) |>
get_summary_stats(isi, type = "mean_sd")
desk
## # A tibble: 15 × 6
## treatment time variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Waitlist Control Baseline isi 30 15.6 3.29
## 2 Waitlist Control Week1 isi 30 15.9 3.47
## 3 Waitlist Control Week2 isi 30 15.3 3.80
## 4 Waitlist Control Week3 isi 30 15.4 3.56
## 5 Waitlist Control Week4 isi 30 15.3 3.45
## 6 Digital BBT-I Baseline isi 30 16.2 3.15
## 7 Digital BBT-I Week1 isi 30 16.2 2.98
## 8 Digital BBT-I Week2 isi 30 15.4 3.39
## 9 Digital BBT-I Week3 isi 30 14.6 2.67
## 10 Digital BBT-I Week4 isi 30 13.7 3.38
## 11 Digital BBT-I + LT Baseline isi 30 15.7 3.52
## 12 Digital BBT-I + LT Week1 isi 30 15.1 3.44
## 13 Digital BBT-I + LT Week2 isi 30 13.8 3.72
## 14 Digital BBT-I + LT Week3 isi 30 12.7 3.44
## 15 Digital BBT-I + LT Week4 isi 30 12.1 3.43
# Profile plot berwarna
cols <- c("Waitlist Control"="#4C78A8", "Digital BBT-I"="#F58518",
"Digital BBT-I + LT"="#54A24B")
ggplot(dat_long, aes(x = time, y = isi, 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 = .15) +
scale_colour_manual(values = cols) +
labs(x = "Waktu", y = "Insomnia Severity Index (ISI)", colour = "Treatment",
title = "Profil rerata ISI (95% CI)") +
theme(legend.position = "bottom")
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
# Spaghetti plot
ggplot(dat_long, aes(x = time, y = isi, group = id)) +
geom_line(alpha = .20) +
stat_summary(aes(group = treatment, colour = treatment), fun = mean,
geom = "line", linewidth = 1.2) +
facet_wrap(~ treatment) +
scale_colour_manual(values = cols) +
labs(x = "Waktu", y = "ISI", title = "Lintasan individu dan rerata treatment") +
theme(legend.position = "none")
# 3. REPEATED MEASURE ANOVA SATU ARAH ---------------------------------------
# Contoh fokus: apakah ISI berubah pada kelompok Digital BBT-I + LT?
d1 <- droplevels(filter(dat_long, treatment == "Digital BBT-I + LT"))
# 3a. Uji asumsi
# Outlier
d1 |> group_by(time) |> identify_outliers(isi)
## # A tibble: 3 × 7
## time id treatment isi week is.outlier is.extreme
## <fct> <fct> <fct> <int> <dbl> <lgl> <lgl>
## 1 Baseline L28 Digital BBT-I + LT 26 0 TRUE FALSE
## 2 Week2 L12 Digital BBT-I + LT 23 2 TRUE FALSE
## 3 Week3 L28 Digital BBT-I + LT 21 3 TRUE FALSE
# Normalitas per waktu
d1 |> group_by(time) |> shapiro_test(isi)
## # A tibble: 5 × 4
## time variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline isi 0.957 0.263
## 2 Week1 isi 0.959 0.292
## 3 Week2 isi 0.971 0.558
## 4 Week3 isi 0.964 0.389
## 5 Week4 isi 0.958 0.274
ggpubr::ggqqplot(d1, "isi", facet.by = "time")
# Sphericity + RM ANOVA
rm1 <- anova_test(data = d1, dv = isi, wid = id, within = time,
effect.size = "pes")
rm1
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 time 4 116 22.753 6.82e-14 * 0.44
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 time 0.641 0.203
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 time 0.826 3.3, 95.79 7.45e-12 * 0.945 3.78, 109.58 3.03e-13
## p[HF]<.05
## 1 *
get_anova_table(rm1, correction = "auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 time 4 116 22.753 6.82e-14 * 0.44
# 3b. Model afex + effect size
aov1 <- aov_ez(id = "id", dv = "isi", data = d1, within = "time")
nice(aov1, correction = "none", es = c("ges","pes"))
## Anova Table (Type 3 tests)
##
## Response: isi
## Effect df MSE F ges pes p.value
## 1 time 4, 116 3.14 22.75 *** .138 .440 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
summary(aov1)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 28925.9 1 1423.07 29 589.465 < 2.2e-16 ***
## time 285.8 4 364.23 116 22.753 6.821e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.64097 0.20335
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.82579 7.451e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.9446145 3.028674e-13
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## time | 0.44 | [0.32, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Post hoc: setiap minggu vs baseline
em1 <- emmeans(aov1, ~ time)
contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm")
## contrast estimate SE df t.ratio p.value
## Week1 - Baseline -0.567 0.380 29 -1.493 0.1463
## Week2 - Baseline -1.867 0.557 29 -3.354 0.0045
## Week3 - Baseline -3.033 0.473 29 -6.408 <0.0001
## Week4 - Baseline -3.600 0.409 29 -8.812 <0.0001
##
## P value adjustment: holm method for 4 tests
# 3c. Alternatif nonparametrik
friedman_test(d1, isi ~ time | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 isi 30 58.4 4 6.40e-12 Friedman test
friedman_effsize(d1, isi ~ time | id)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 isi 30 0.486 Kendall W moderate
d1 |> wilcox_test(isi ~ time, paired = TRUE, p.adjust.method = "holm")
## # A tibble: 10 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 isi Baseline Week1 30 30 291 0.0924 1.85e-1 ns
## 2 isi Baseline Week2 30 30 380 0.00102 6.15e-3 **
## 3 isi Baseline Week3 30 30 446 0.000000387 2.71e-6 ****
## 4 isi Baseline Week4 30 30 459 0.0000000149 1.49e-7 ****
## 5 isi Week1 Week2 30 30 363 0.00346 1.73e-2 *
## 6 isi Week1 Week3 30 30 437 0.000000238 1.91e-6 ****
## 7 isi Week1 Week4 30 30 453 0.0000000969 8.72e-7 ****
## 8 isi Week2 Week3 30 30 336 0.0193 5.80e-2 ns
## 9 isi Week2 Week4 30 30 359 0.00349 1.73e-2 *
## 10 isi Week3 Week4 30 30 287 0.253 2.53e-1 ns
# 4. MIXED DESIGN ANOVA: TREATMENT x TIME -----------------------------------
# 4a. Uji asumsi
# Normalitas per sel
dat_long |> group_by(treatment, time) |> shapiro_test(isi)
## # A tibble: 15 × 5
## treatment time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Waitlist Control Baseline isi 0.959 0.291
## 2 Waitlist Control Week1 isi 0.966 0.427
## 3 Waitlist Control Week2 isi 0.934 0.0623
## 4 Waitlist Control Week3 isi 0.977 0.752
## 5 Waitlist Control Week4 isi 0.957 0.258
## 6 Digital BBT-I Baseline isi 0.965 0.422
## 7 Digital BBT-I Week1 isi 0.963 0.361
## 8 Digital BBT-I Week2 isi 0.964 0.401
## 9 Digital BBT-I Week3 isi 0.966 0.437
## 10 Digital BBT-I Week4 isi 0.968 0.490
## 11 Digital BBT-I + LT Baseline isi 0.957 0.263
## 12 Digital BBT-I + LT Week1 isi 0.959 0.292
## 13 Digital BBT-I + LT Week2 isi 0.971 0.558
## 14 Digital BBT-I + LT Week3 isi 0.964 0.389
## 15 Digital BBT-I + LT Week4 isi 0.958 0.274
# Homogenitas varians per waktu
dat_long |> group_by(time) |> levene_test(isi ~ treatment)
## # A tibble: 5 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 87 0.113 0.893
## 2 Week1 2 87 0.436 0.648
## 3 Week2 2 87 0.314 0.731
## 4 Week3 2 87 1.06 0.349
## 5 Week4 2 87 0.00268 0.997
# Homogenitas matriks kovarians (Box's M)
box_m(dat_wide[, c("ISI_Baseline","ISI_Week1","ISI_Week2","ISI_Week3","ISI_Week4")],
dat_wide$treatment)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 31.7 0.383 30 Box's M-test for Homogeneity of Covariance Matric…
# Mauchly + ANOVA campuran melalui rstatix
aov2_rs <- anova_test(data = dat_long, dv = isi, wid = id,
between = treatment, within = time,
effect.size = "pes", type = 3)
aov2_rs
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 treatment 2 87 2.485 8.9e-02 0.054
## 2 time 4 348 22.847 8.9e-17 * 0.208
## 3 treatment:time 8 348 4.407 4.3e-05 * 0.092
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 time 0.875 0.25
## 2 treatment:time 0.875 0.25
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF]
## 1 time 0.935 3.74, 325.53 7.60e-16 * 0.983 3.93, 341.98
## 2 treatment:time 0.935 7.48, 325.53 7.02e-05 * 0.983 7.86, 341.98
## p[HF] p[HF]<.05
## 1 1.58e-16 *
## 2 4.90e-05 *
get_anova_table(aov2_rs, correction = "auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 treatment 2 87 2.485 8.9e-02 0.054
## 2 time 4 348 22.847 8.9e-17 * 0.208
## 3 treatment:time 8 348 4.407 4.3e-05 * 0.092
# afex sebagai model utama pembanding
aov2 <- aov_ez(id = "id", dv = "isi", data = dat_long,
between = "treatment", within = "time")
nice(aov2, correction = "none", es = c("ges","pes"))
## Anova Table (Type 3 tests)
##
## Response: isi
## Effect df MSE F ges pes p.value
## 1 treatment 2, 87 43.79 2.49 + .042 .054 .089
## 2 time 4, 348 3.42 22.85 *** .059 .208 <.001
## 3 treatment:time 8, 348 3.42 4.41 *** .024 .092 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
summary(aov2)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 99309 1 3809.9 87 2267.7304 < 2.2e-16 ***
## treatment 218 2 3809.9 87 2.4853 0.0892 .
## time 313 4 1190.7 348 22.8471 < 2.2e-16 ***
## treatment:time 121 8 1190.7 348 4.4066 4.297e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.87519 0.25017
## treatment:time 0.87519 0.25017
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.93544 7.602e-16 ***
## treatment:time 0.93544 7.025e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.9827148 1.580701e-16
## treatment:time 0.9827148 4.900641e-05
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.21 | [0.14, 1.00]
## treatment:time | 0.09 | [0.03, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Plot interaksi
afex_plot(aov2, x = "time", trace = "treatment", error = "within",
mapping = c("colour","shape","linetype")) +
scale_colour_manual(values = cols) +
labs(x = "Waktu", y = "ISI")
## Warning: Panel(s) show a mixed within-between-design.
## Error bars do not allow comparisons across all means.
## Suppress error bars with: error = "none"
# 4c. Simple effects + post hoc
joint_tests(aov2, by = "treatment")
## treatment = Waitlist Control:
## model term df1 df2 F.ratio p.value
## time 4 87 0.693 0.5990
##
## treatment = Digital BBT-I:
## model term df1 df2 F.ratio p.value
## time 4 87 11.085 <0.0001
##
## treatment = Digital BBT-I + LT:
## model term df1 df2 F.ratio p.value
## time 4 87 22.445 <0.0001
joint_tests(aov2, by = "time")
## time = Baseline:
## model term df1 df2 F.ratio p.value
## treatment 2 87 0.248 0.7807
##
## time = Week1:
## model term df1 df2 F.ratio p.value
## treatment 2 87 0.775 0.4636
##
## time = Week2:
## model term df1 df2 F.ratio p.value
## treatment 2 87 1.707 0.1873
##
## time = Week3:
## model term df1 df2 F.ratio p.value
## treatment 2 87 5.697 0.0047
##
## time = Week4:
## model term df1 df2 F.ratio p.value
## treatment 2 87 6.433 0.0025
em2 <- emmeans(aov2, ~ time | treatment)
contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm")
## treatment = Waitlist Control:
## contrast estimate SE df t.ratio p.value
## Week1 - Baseline 0.267 0.477 87 0.559 1.0000
## Week2 - Baseline -0.333 0.539 87 -0.619 1.0000
## Week3 - Baseline -0.167 0.501 87 -0.333 1.0000
## Week4 - Baseline -0.333 0.445 87 -0.749 1.0000
##
## treatment = Digital BBT-I:
## contrast estimate SE df t.ratio p.value
## Week1 - Baseline 0.000 0.477 87 0.000 1.0000
## Week2 - Baseline -0.767 0.539 87 -1.423 0.3166
## Week3 - Baseline -1.600 0.501 87 -3.196 0.0058
## Week4 - Baseline -2.500 0.445 87 -5.615 <0.0001
##
## treatment = Digital BBT-I + LT:
## contrast estimate SE df t.ratio p.value
## Week1 - Baseline -0.567 0.477 87 -1.189 0.2378
## Week2 - Baseline -1.867 0.539 87 -3.465 0.0017
## Week3 - Baseline -3.033 0.501 87 -6.058 <0.0001
## Week4 - Baseline -3.600 0.445 87 -8.086 <0.0001
##
## P value adjustment: holm method for 4 tests
em2b <- emmeans(aov2, ~ treatment | time)
pairs(em2b, adjust = "tukey")
## time = Baseline:
## contrast estimate SE df t.ratio p.value
## Waitlist Control - (Digital BBT-I) -0.567 0.858 87 -0.660 0.7871
## Waitlist Control - (Digital BBT-I + LT) -0.100 0.858 87 -0.116 0.9925
## (Digital BBT-I) - (Digital BBT-I + LT) 0.467 0.858 87 0.544 0.8500
##
## time = Week1:
## contrast estimate SE df t.ratio p.value
## Waitlist Control - (Digital BBT-I) -0.300 0.854 87 -0.351 0.9342
## Waitlist Control - (Digital BBT-I + LT) 0.733 0.854 87 0.859 0.6674
## (Digital BBT-I) - (Digital BBT-I + LT) 1.033 0.854 87 1.210 0.4503
##
## time = Week2:
## contrast estimate SE df t.ratio p.value
## Waitlist Control - (Digital BBT-I) -0.133 0.940 87 -0.142 0.9890
## Waitlist Control - (Digital BBT-I + LT) 1.433 0.940 87 1.525 0.2844
## (Digital BBT-I) - (Digital BBT-I + LT) 1.567 0.940 87 1.667 0.2239
##
## time = Week3:
## contrast estimate SE df t.ratio p.value
## Waitlist Control - (Digital BBT-I) 0.867 0.838 87 1.034 0.5578
## Waitlist Control - (Digital BBT-I + LT) 2.767 0.838 87 3.300 0.0040
## (Digital BBT-I) - (Digital BBT-I + LT) 1.900 0.838 87 2.266 0.0661
##
## time = Week4:
## contrast estimate SE df t.ratio p.value
## Waitlist Control - (Digital BBT-I) 1.600 0.883 87 1.812 0.1716
## Waitlist Control - (Digital BBT-I + LT) 3.167 0.883 87 3.587 0.0016
## (Digital BBT-I) - (Digital BBT-I + LT) 1.567 0.883 87 1.774 0.1842
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4d. Kontras perubahan Week4 - Baseline antar-treatment
em_full <- emmeans(aov2, ~ time * treatment)
contrast(em_full,
interaction = list(time = list("Week4-Baseline" = c(-1,0,0,0,1)),
treatment = "pairwise"),
adjust = "holm")
## time_custom treatment_pairwise estimate SE df
## Week4-Baseline Waitlist Control - (Digital BBT-I) 2.17 0.63 87
## Week4-Baseline Waitlist Control - (Digital BBT-I + LT) 3.27 0.63 87
## Week4-Baseline (Digital BBT-I) - (Digital BBT-I + LT) 1.10 0.63 87
## t.ratio p.value
## 3.441 0.0018
## 5.188 <0.0001
## 1.747 0.0842
##
## P value adjustment: holm method for 3 tests
# 5. LINEAR MIXED MODEL -------------------------------------------------------
# Pembanding yang lebih fleksibel terhadap data hilang.
lmm <- lmer(isi ~ 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 17.007 8.503 2 87 2.4853 0.0892 .
## time 312.689 78.172 4 348 22.8471 < 2.2e-16 ***
## treatment:time 120.618 15.077 8 348 4.4066 4.297e-05 ***
## ---
## 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: isi ~ treatment * time + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 2063
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1836 -0.5842 -0.0309 0.6476 3.3044
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 8.074 2.842
## Residual 3.422 1.850
## Number of obs: 450, groups: id, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.855556 0.311956 86.999995 47.621 < 2e-16 ***
## treatment1 0.631111 0.441172 86.999995 1.431 0.156146
## treatment2 0.337778 0.441172 86.999995 0.766 0.445965
## time1 0.966667 0.174395 348.000002 5.543 5.88e-08 ***
## time2 0.866667 0.174395 348.000002 4.970 1.05e-06 ***
## time3 -0.022222 0.174395 348.000002 -0.127 0.898678
## time4 -0.633333 0.174395 348.000002 -3.632 0.000324 ***
## treatment1:time1 -0.853333 0.246632 348.000002 -3.460 0.000607 ***
## treatment2:time1 0.006667 0.246632 348.000002 0.027 0.978451
## treatment1:time2 -0.486667 0.246632 348.000002 -1.973 0.049258 *
## treatment2:time2 0.106667 0.246632 348.000002 0.432 0.665651
## treatment1:time3 -0.197778 0.246632 348.000002 -0.802 0.423149
## treatment2:time3 0.228889 0.246632 348.000002 0.928 0.354020
## treatment1:time4 0.580000 0.246632 348.000002 2.352 0.019246 *
## treatment2:time4 0.006667 0.246632 348.000002 0.027 0.978451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
performance::icc(lmm)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.702
## Unadjusted ICC: 0.624
# Diagnostik residual LMM
par(mfrow = c(1,2))
qqnorm(resid(lmm), main = "Q-Q residual LMM"); qqline(resid(lmm))
plot(fitted(lmm), resid(lmm), xlab = "Nilai prediksi", ylab = "Residual",
main = "Residual vs prediksi"); abline(h = 0, lty = 2)
par(mfrow = c(1,1))
# 6. CATATAN AKADEMIK ---------------------------------------------------------
# Artikel sumber adalah PROTOKOL RCT. Karena hasil uji klinis belum dilaporkan,
# angka individual, mean tiap minggu, F, p, effect size, dan post hoc pada latihan
# ini SELURUHNYA berasal dari data simulasi. Yang diadopsi dari artikel ialah:
# desain 3-arm RCT, outcome utama ISI, rencana n=90 (30 per group), dan struktur
# pengukuran primer baseline + minggu 1-4. Artikel juga merencanakan follow-up 3 bulan.