# =============================================================================
# Nama : Rizky Katherine
# NIM : 2611018011
# =============================================================================
# REPEATED MEASURES ANALYSIS DENGAN R
# Topik : Perubahan skor FMA-UE pada 3 kelompok selama 12 minggu
# DATA : CSV LONG dan CSV WIDE
#
# Struktur data:
# - 90 responden (30 per kelompok)
# - 3 kelompok: CRT, BWS_TCY, RAT
# - 4 kali pengukuran: minggu 0, 4, 8, 12
# - Outcome: FMA_UE (0-66)
#
# Catatan:
# Data merupakan DATA SIMULASI untuk latihan Biostatistik dan bukan
# data individu asli dari artikel Zhang et al. (2025).
#
# Alur analisis mengikuti format syntax terlampir:
# 0. Paket & pengaturan
# 1. Membaca data CSV LONG & WIDE + validasi
# 2. Eksplorasi data + grafik
# 3. Repeated Measures ANOVA satu arah (contoh BWS_TCY)
# 3a. Outlier, normalitas, sphericity/Mauchly
# 3b. ANOVA + Greenhouse-Geisser/Huynh-Feldt
# 3c. Pendekatan multivariat
# 3d. Post-hoc dan tren
# 3e. Friedman
# 4. Mixed Design ANOVA (Group x Week)
# 4a. Asumsi: outlier, normalitas, Levene, Box's M, Mauchly
# 4b. ANOVA + effect size
# 4c. Simple effects + post-hoc
# 4d. Kontras perubahan baseline -> minggu 12
# 5. Linear Mixed Model (LMM)
# 6. Export hasil
# =============================================================================
# 0. PAKET & PENGATURAN
# Jalankan sekali bila paket belum tersedia:
# install.packages(c("dplyr", "tidyr", "ggplot2", "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 CSV LONG DAN WIDE
# Pastikan kedua file CSV berada di Working Directory.
# Cek dengan:
getwd()
## [1] "C:/Users/kathe_5vwnplr/Downloads/DAFTAR S2/SEMESTER 1/BIOSTATISTIK/R STUDIO-REPEATED MEASUREMENT ANOVA"
list.files()
## [1] "dataset_repeated_measures_3group_4time_LONG.csv"
## [2] "dataset_repeated_measures_3group_4time_WIDE.csv"
## [3] "grafik_01_profile_FMA_UE.png"
## [4] "grafik_02_spaghetti_FMA_UE.png"
## [5] "grafik_03_interaksi_Group_x_Week.png"
## [6] "hasil_01_statistik_deskriptif_FMA_UE.csv"
## [7] "hasil_02_shapiro_Group_x_Week.csv"
## [8] "hasil_03_levene_per_Week.csv"
## [9] "hasil_04_mixed_ANOVA_GG.csv"
## [10] "hasil_05_posthoc_antar_kelompok_Holm.csv"
## [11] "hasil_06_posthoc_dalam_kelompok_Holm.csv"
## [12] "hasil_07_estimated_marginal_means.csv"
## [13] "Rizky Katherine_2611018011_Repeated Measurement ANOVA.docx"
## [14] "RStudio_Repeated_Measures_FMA_UE_FINAL_CSV_LONG_WIDE.docx"
## [15] "RStudio_Repeated_Measures_FMA_UE_FINAL_CSV_LONG_WIDE.R"
## [16] "RStudio_Repeated_Measures_FMA_UE_FINAL_CSV_LONG_WIDE.spin.R"
## [17] "RStudio_Repeated_Measures_FMA_UE_FINAL_CSV_LONG_WIDE.spin.Rmd"
# Jika CSV belum ditemukan, atur Working Directory melalui:
# Session -> Set Working Directory -> Choose Directory...
# 1A. Baca CSV LONG
data_long <- read.csv(
"dataset_repeated_measures_3group_4time_LONG.csv",
header = TRUE,
stringsAsFactors = FALSE
)
# 1B. Baca CSV WIDE
data_wide <- read.csv(
"dataset_repeated_measures_3group_4time_WIDE.csv",
header = TRUE,
stringsAsFactors = FALSE
)
# Jika nama/path file berbeda, gunakan:
# data_long <- read.csv(file.choose(), header = TRUE)
# data_wide <- read.csv(file.choose(), header = TRUE)
# Tampilkan data
View(data_long)
View(data_wide)
# Struktur dan dimensi
head(data_long)
## ID Group Week FMA_UE
## 1 P001 CRT 0 35.5
## 2 P001 CRT 4 37.7
## 3 P001 CRT 8 38.7
## 4 P001 CRT 12 39.5
## 5 P002 CRT 0 34.8
## 6 P002 CRT 4 37.6
head(data_wide)
## ID Group FMA_UE_W0 FMA_UE_W4 FMA_UE_W8 FMA_UE_W12
## 1 P001 CRT 35.5 37.7 38.7 39.5
## 2 P002 CRT 34.8 37.6 36.5 46.3
## 3 P003 CRT 28.8 33.5 39.4 38.2
## 4 P004 CRT 36.1 46.3 42.6 48.3
## 5 P005 CRT 24.2 32.5 37.7 34.0
## 6 P006 CRT 37.1 39.1 37.1 46.2
str(data_long)
## 'data.frame': 360 obs. of 4 variables:
## $ ID : chr "P001" "P001" "P001" "P001" ...
## $ Group : chr "CRT" "CRT" "CRT" "CRT" ...
## $ Week : int 0 4 8 12 0 4 8 12 0 4 ...
## $ FMA_UE: num 35.5 37.7 38.7 39.5 34.8 37.6 36.5 46.3 28.8 33.5 ...
str(data_wide)
## 'data.frame': 90 obs. of 6 variables:
## $ ID : chr "P001" "P002" "P003" "P004" ...
## $ Group : chr "CRT" "CRT" "CRT" "CRT" ...
## $ FMA_UE_W0 : num 35.5 34.8 28.8 36.1 24.2 37.1 29.6 26 36.4 27 ...
## $ FMA_UE_W4 : num 37.7 37.6 33.5 46.3 32.5 39.1 35.9 36.3 36.5 36 ...
## $ FMA_UE_W8 : num 38.7 36.5 39.4 42.6 37.7 37.1 32.1 36.3 44.6 43.6 ...
## $ FMA_UE_W12: num 39.5 46.3 38.2 48.3 34 46.2 43.7 36.7 38.4 39.3 ...
dim(data_long)
## [1] 360 4
dim(data_wide)
## [1] 90 6
names(data_long)
## [1] "ID" "Group" "Week" "FMA_UE"
names(data_wide)
## [1] "ID" "Group" "FMA_UE_W0" "FMA_UE_W4" "FMA_UE_W8"
## [6] "FMA_UE_W12"
# 1C. VALIDASI STRUKTUR DATA
# Kolom yang diharapkan pada LONG
stopifnot(
all(c("ID", "Group", "Week", "FMA_UE") %in% names(data_long))
)
# Kolom yang diharapkan pada WIDE
stopifnot(
all(c("ID", "Group", "FMA_UE_W0", "FMA_UE_W4",
"FMA_UE_W8", "FMA_UE_W12") %in% names(data_wide))
)
# Pastikan jumlah responden 90
n_distinct(data_long$ID)
## [1] 90
n_distinct(data_wide$ID)
## [1] 90
# Jumlah responden tiap kelompok
data_wide %>% count(Group)
## Group n
## 1 BWS_TCY 30
## 2 CRT 30
## 3 RAT 30
data_long %>%
distinct(ID, Group) %>%
count(Group)
## Group n
## 1 BWS_TCY 30
## 2 CRT 30
## 3 RAT 30
# Jumlah pengukuran tiap responden
data_long %>%
count(ID) %>%
count(n, name = "jumlah_responden")
## n jumlah_responden
## 1 4 90
# Cek duplikasi ID x Week
cek_duplikasi <- data_long %>%
count(ID, Week) %>%
filter(n != 1)
cek_duplikasi
## [1] ID Week n
## <0 rows> (or 0-length row.names)
# Cek missing value
colSums(is.na(data_long))
## ID Group Week FMA_UE
## 0 0 0 0
colSums(is.na(data_wide))
## ID Group FMA_UE_W0 FMA_UE_W4 FMA_UE_W8 FMA_UE_W12
## 0 0 0 0 0 0
# Ringkasan outcome
summary(data_long$FMA_UE)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 22.80 36.17 41.15 41.65 46.90 63.40
# 1D. RECODING VARIABEL
# LONG digunakan sebagai basis analisis.
# Week_num tetap numerik untuk grafik/LMM,
# sedangkan Week menjadi faktor untuk repeated-measures ANOVA.
data_long <- data_long %>%
mutate(
ID = factor(ID),
Group = factor(
Group,
levels = c("CRT", "BWS_TCY", "RAT")
),
Week_num = as.numeric(Week),
Week = factor(
Week,
levels = c(0, 4, 8, 12),
labels = c("W0", "W4", "W8", "W12")
)
)
# WIDE juga diberi factor untuk keperluan Box's M / pengecekan.
data_wide <- data_wide %>%
mutate(
ID = factor(ID),
Group = factor(
Group,
levels = c("CRT", "BWS_TCY", "RAT")
)
)
# Cek kembali level
levels(data_long$Group)
## [1] "CRT" "BWS_TCY" "RAT"
levels(data_long$Week)
## [1] "W0" "W4" "W8" "W12"
# 1E. CEK KONSISTENSI LONG VS WIDE
# Buat WIDE dari LONG untuk membandingkan dengan file WIDE.
wide_from_long <- data_long %>%
select(ID, Group, Week, FMA_UE) %>%
pivot_wider(
names_from = Week,
values_from = FMA_UE,
names_prefix = "FMA_UE_"
)
# Bila ingin melihat hasil transformasi:
View(wide_from_long)
# Pengecekan sederhana jumlah baris
nrow(wide_from_long)
## [1] 90
nrow(data_wide)
## [1] 90
# 2. EKSPLORASI DATA
# 2A. Statistik deskriptif per kelompok dan waktu
deskriptif <- data_long %>%
group_by(Group, Week) %>%
summarise(
n = n(),
mean = mean(FMA_UE, na.rm = TRUE),
sd = sd(FMA_UE, na.rm = TRUE),
median = median(FMA_UE, na.rm = TRUE),
min = min(FMA_UE, na.rm = TRUE),
max = max(FMA_UE, na.rm = TRUE),
.groups = "drop"
)
deskriptif
## # A tibble: 12 × 8
## Group Week n mean sd median min max
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 CRT W0 30 32.8 5.35 33.7 22.8 41.2
## 2 CRT W4 30 36.9 4.85 37.0 27.6 47.6
## 3 CRT W8 30 40.1 4.41 41.1 32.1 49.8
## 4 CRT W12 30 42.2 4.88 41.8 33.2 53.7
## 5 BWS_TCY W0 30 35.6 4.83 35.5 24.6 46.3
## 6 BWS_TCY W4 30 41.1 4.99 40.9 30 51.4
## 7 BWS_TCY W8 30 46.6 4.91 46.0 35.5 54.6
## 8 BWS_TCY W12 30 53.1 4.69 52.4 43.5 59.9
## 9 RAT W0 30 36.0 5.82 34.6 26.6 49.9
## 10 RAT W4 30 39.4 5.28 39.4 28.7 50.7
## 11 RAT W8 30 45.9 5.34 46.8 35.8 57.5
## 12 RAT W12 30 50.0 5.72 49.8 38.9 63.4
# Versi rstatix
data_long %>%
group_by(Group, Week) %>%
get_summary_stats(FMA_UE, type = "mean_sd")
## # A tibble: 12 × 6
## Group Week variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 CRT W0 FMA_UE 30 32.8 5.34
## 2 CRT W4 FMA_UE 30 36.9 4.85
## 3 CRT W8 FMA_UE 30 40.1 4.41
## 4 CRT W12 FMA_UE 30 42.2 4.88
## 5 BWS_TCY W0 FMA_UE 30 35.6 4.83
## 6 BWS_TCY W4 FMA_UE 30 41.1 4.99
## 7 BWS_TCY W8 FMA_UE 30 46.6 4.91
## 8 BWS_TCY W12 FMA_UE 30 53.1 4.69
## 9 RAT W0 FMA_UE 30 36.0 5.82
## 10 RAT W4 FMA_UE 30 39.4 5.28
## 11 RAT W8 FMA_UE 30 45.9 5.34
## 12 RAT W12 FMA_UE 30 50.0 5.72
# 2B. Matriks kovarians dan korelasi antar waktu
vars_waktu <- c(
"FMA_UE_W0", "FMA_UE_W4", "FMA_UE_W8", "FMA_UE_W12"
)
S <- cov(data_wide[, vars_waktu], use = "complete.obs")
R <- cor(data_wide[, vars_waktu], use = "complete.obs")
round(S, 2)
## FMA_UE_W0 FMA_UE_W4 FMA_UE_W8 FMA_UE_W12
## FMA_UE_W0 30.08 22.52 23.73 25.80
## FMA_UE_W4 22.52 27.79 22.37 27.28
## FMA_UE_W8 23.73 22.37 32.03 30.60
## FMA_UE_W12 25.80 27.28 30.60 46.71
round(R, 3)
## FMA_UE_W0 FMA_UE_W4 FMA_UE_W8 FMA_UE_W12
## FMA_UE_W0 1.000 0.779 0.765 0.688
## FMA_UE_W4 0.779 1.000 0.750 0.757
## FMA_UE_W8 0.765 0.750 1.000 0.791
## FMA_UE_W12 0.688 0.757 0.791 1.000
# 2C. Varians selisih antar waktu
pasangan <- combn(vars_waktu, 2)
var_selisih <- apply(
pasangan,
2,
function(p) var(data_wide[[p[1]]] - data_wide[[p[2]]], na.rm = TRUE)
)
names(var_selisih) <- apply(
pasangan,
2,
paste,
collapse = " - "
)
round(var_selisih, 2)
## FMA_UE_W0 - FMA_UE_W4 FMA_UE_W0 - FMA_UE_W8 FMA_UE_W0 - FMA_UE_W12
## 12.82 14.65 25.19
## FMA_UE_W4 - FMA_UE_W8 FMA_UE_W4 - FMA_UE_W12 FMA_UE_W8 - FMA_UE_W12
## 15.09 19.94 17.54
# 2D. Profile plot: rerata +/- 95% CI
p_profil <- ggplot(
data_long,
aes(
x = Week_num,
y = FMA_UE,
colour = Group,
group = Group
)
) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", size = 2.8) +
stat_summary(
fun.data = mean_cl_normal,
geom = "errorbar",
width = .5
) +
scale_x_continuous(breaks = c(0, 4, 8, 12)) +
labs(
x = "Minggu pengukuran",
y = "Skor FMA-UE",
colour = "Kelompok",
title = "Profil rerata FMA-UE selama 12 minggu"
) +
theme(legend.position = "bottom")
p_profil
# 2E. Spaghetti plot
p_spag <- ggplot(
data_long,
aes(
x = Week_num,
y = FMA_UE,
group = ID
)
) +
geom_line(alpha = .25) +
stat_summary(
aes(group = Group),
fun = mean,
geom = "line",
linewidth = 1.3
) +
facet_wrap(~ Group) +
scale_x_continuous(breaks = c(0, 4, 8, 12)) +
labs(
x = "Minggu",
y = "FMA-UE",
title = "Lintasan individu dan rerata kelompok"
)
p_spag
# 3. REPEATED MEASURES ANOVA SATU ARAH
# Contoh: kelompok BWS_TCY
# Pertanyaan:
# Apakah skor FMA-UE berubah selama 12 minggu pada kelompok BWS_TCY?
d1 <- data_long %>%
filter(Group == "BWS_TCY") %>%
droplevels()
d1w <- data_wide %>%
filter(Group == "BWS_TCY") %>%
droplevels()
# 3A. UJI ASUMSI
# (i) Outlier per waktu
outlier_1way <- d1 %>%
group_by(Week) %>%
identify_outliers(FMA_UE)
outlier_1way
## [1] Week ID Group FMA_UE Week_num is.outlier is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas Shapiro-Wilk per waktu
shapiro_1way <- d1 %>%
group_by(Week) %>%
shapiro_test(FMA_UE)
shapiro_1way
## # A tibble: 4 × 4
## Week variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 W0 FMA_UE 0.990 0.990
## 2 W4 FMA_UE 0.977 0.748
## 3 W8 FMA_UE 0.963 0.373
## 4 W12 FMA_UE 0.948 0.147
# Q-Q plot
qqplot_1way <- ggpubr::ggqqplot(
d1,
"FMA_UE",
facet.by = "Week"
)
qqplot_1way
# (iii) Mauchly's Test
# anova_test melaporkan Mauchly + GG/HF
aov1_rs <- anova_test(
data = d1,
dv = FMA_UE,
wid = ID,
within = Week,
effect.size = "pes"
)
aov1_rs
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 Week 3 87 293.332 2.28e-45 * 0.91
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 Week 0.962 0.956
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 Week 0.975 2.93, 84.86 2.6e-44 * 1.097 3.29, 95.41 2.28e-45
## 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 Week 3 87 293.332 2.28e-45 * 0.91
# 3B. REPEATED MEASURES ANOVA DENGAN AFEX
aov1 <- aov_ez(
id = "ID",
dv = "FMA_UE",
data = d1,
within = "Week",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov1
## Anova Table (Type 3 tests)
##
## Response: FMA_UE
## Effect df MSE F ges pes p.value
## 1 Week 2.93, 84.86 5.88 293.33 *** .648 .910 <.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) 233483 1 2236.74 29 3027.17 < 2.2e-16 ***
## Week 5047 3 498.93 87 293.33 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## Week 0.96171 0.95568
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Week 0.97538 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Week 1.096687 2.283027e-45
nice(aov1, correction = "GG", es = c("ges", "pes"))
## Anova Table (Type 3 tests)
##
## Response: FMA_UE
## Effect df MSE F ges pes p.value
## 1 Week 2.93, 84.86 5.88 293.33 *** .648 .910 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## Week | 0.91 | [0.88, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
omega_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## -------------------------------------------
## Week | 0.64 | [0.54, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 3C. PENDEKATAN MULTIVARIAT
aov1$Anova
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99051 3027.17 1 29 < 2.2e-16 ***
## Week 1 0.96743 267.29 3 27 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3D. POST-HOC DAN TREN
em1 <- emmeans(aov1, ~ Week)
em1
## Week emmean SE df lower.CL upper.CL
## W0 35.6 0.881 29 33.8 37.4
## W4 41.1 0.912 29 39.2 43.0
## W8 46.6 0.897 29 44.7 48.4
## W12 53.1 0.856 29 51.4 54.9
##
## Confidence level used: 0.95
# Semua pasangan waktu
pairs(em1, adjust = "holm")
## contrast estimate SE df t.ratio p.value
## W0 - W4 -5.44 0.617 29 -8.820 <0.0001
## W0 - W8 -10.92 0.591 29 -18.495 <0.0001
## W0 - W12 -17.49 0.649 29 -26.939 <0.0001
## W4 - W8 -5.48 0.559 29 -9.804 <0.0001
## W4 - W12 -12.04 0.639 29 -18.858 <0.0001
## W8 - W12 -6.56 0.650 29 -10.097 <0.0001
##
## P value adjustment: holm method for 6 tests
# Masing-masing waktu dibandingkan dengan baseline W0
contrast(
em1,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## contrast estimate SE df t.ratio p.value
## W4 - W0 5.44 0.617 29 8.820 <0.0001
## W8 - W0 10.92 0.591 29 18.495 <0.0001
## W12 - W0 17.49 0.649 29 26.939 <0.0001
##
## P value adjustment: holm method for 3 tests
# Tren linear, kuadratik, kubik
contrast(em1, "poly")
## contrast estimate SE df t.ratio p.value
## linear 57.94 1.990 29 29.101 <0.0001
## quadratic 1.12 0.909 29 1.232 0.2278
## cubic 1.05 1.840 29 0.570 0.5732
# 3E. ALTERNATIF NONPARAMETRIK: FRIEDMAN
friedman_test(d1, FMA_UE ~ Week | ID)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 FMA_UE 30 86.5 3 1.22e-18 Friedman test
friedman_effsize(d1, FMA_UE ~ Week | ID)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 FMA_UE 30 0.961 Kendall W large
# Wilcoxon berpasangan sebagai post-hoc alternatif
# bila diperlukan
d1 %>%
wilcox_test(
FMA_UE ~ Week,
paired = TRUE,
p.adjust.method = "holm"
)
## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 FMA_UE W0 W4 30 30 5 0.0000000168 1.68e-8 ****
## 2 FMA_UE W0 W8 30 30 0 0.00000000186 1.12e-8 ****
## 3 FMA_UE W0 W12 30 30 0 0.00000000186 1.12e-8 ****
## 4 FMA_UE W4 W8 30 30 2 0.00000000559 1.12e-8 ****
## 5 FMA_UE W4 W12 30 30 0 0.00000000186 1.12e-8 ****
## 6 FMA_UE W8 W12 30 30 0 0.00000000186 1.12e-8 ****
# 4. MIXED DESIGN ANOVA
# Group (between) x Week (within)
# Pertanyaan utama:
# Apakah perubahan FMA-UE dari minggu 0 sampai minggu 12 berbeda
# antar kelompok CRT, BWS_TCY, dan RAT?
# 4A. UJI ASUMSI
# (i) Outlier per sel
outlier_mixed <- data_long %>%
group_by(Group, Week) %>%
identify_outliers(FMA_UE)
outlier_mixed
## # A tibble: 1 × 7
## Group Week ID FMA_UE Week_num is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 RAT W12 P088 63.4 12 TRUE FALSE
# (ii) Normalitas Shapiro-Wilk per Group x Week
shapiro_mixed <- data_long %>%
group_by(Group, Week) %>%
shapiro_test(FMA_UE)
shapiro_mixed
## # A tibble: 12 × 5
## Group Week variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 CRT W0 FMA_UE 0.951 0.182
## 2 CRT W4 FMA_UE 0.972 0.599
## 3 CRT W8 FMA_UE 0.974 0.668
## 4 CRT W12 FMA_UE 0.970 0.534
## 5 BWS_TCY W0 FMA_UE 0.990 0.990
## 6 BWS_TCY W4 FMA_UE 0.977 0.748
## 7 BWS_TCY W8 FMA_UE 0.963 0.373
## 8 BWS_TCY W12 FMA_UE 0.948 0.147
## 9 RAT W0 FMA_UE 0.958 0.268
## 10 RAT W4 FMA_UE 0.991 0.995
## 11 RAT W8 FMA_UE 0.984 0.920
## 12 RAT W12 FMA_UE 0.979 0.794
# Q-Q plot per sel
qqplot_mixed <- ggpubr::ggqqplot(
data_long,
"FMA_UE"
) +
facet_grid(Week ~ Group)
qqplot_mixed
# (iii) Homogenitas varians pada setiap waktu: Levene
levene_mixed <- data_long %>%
group_by(Week) %>%
levene_test(FMA_UE ~ Group)
levene_mixed
## # A tibble: 4 × 5
## Week df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 W0 2 87 0.448 0.640
## 2 W4 2 87 0.361 0.698
## 3 W8 2 87 0.302 0.740
## 4 W12 2 87 0.359 0.699
# (iv) Homogenitas matriks kovarians: Box's M
box_m_result <- box_m(
data_wide[, vars_waktu],
data_wide$Group
)
box_m_result
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 12.0 0.915 20 Box's M-test for Homogeneity of Covariance Matric…
# (v) Sphericity: Mauchly dilihat pada summary(aov2)
# 4B. MIXED DESIGN REPEATED MEASURES ANOVA
aov2 <- aov_ez(
id = "ID",
dv = "FMA_UE",
data = data_long,
between = "Group",
within = "Week",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov2
## Anova Table (Type 3 tests)
##
## Response: FMA_UE
## Effect df MSE F ges pes p.value
## 1 Group 2, 87 84.39 14.67 *** .214 .252 <.001
## 2 Week 2.94, 255.76 6.74 480.23 *** .512 .847 <.001
## 3 Group:Week 5.88, 255.76 6.74 15.54 *** .064 .263 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2)
## 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) 624458 1 7341.7 87 7399.913 < 2.2e-16 ***
## Group 2475 2 7341.7 87 14.666 3.244e-06 ***
## Week 9522 3 1725.1 261 480.225 < 2.2e-16 ***
## Group:Week 616 6 1725.1 261 15.538 3.086e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## Week 0.96997 0.75933
## Group:Week 0.96997 0.75933
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Week 0.97994 < 2.2e-16 ***
## Group:Week 0.97994 5.63e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Week 1.017933 6.573898e-106
## Group:Week 1.017933 3.086179e-15
nice(aov2, correction = "GG", es = c("ges", "pes"))
## Anova Table (Type 3 tests)
##
## Response: FMA_UE
## Effect df MSE F ges pes p.value
## 1 Group 2, 87 84.39 14.67 *** .214 .252 <.001
## 2 Week 2.94, 255.76 6.74 480.23 *** .512 .847 <.001
## 3 Group:Week 5.88, 255.76 6.74 15.54 *** .064 .263 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Alternatif rstatix
aov2_rs <- anova_test(
data = data_long,
dv = FMA_UE,
wid = ID,
between = Group,
within = Week,
effect.size = "pes",
type = 3
)
aov2_rs
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 Group 2 87 14.666 3.24e-06 * 0.252
## 2 Week 3 261 480.225 6.57e-106 * 0.847
## 3 Group:Week 6 261 15.538 3.09e-15 * 0.263
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 Week 0.97 0.759
## 2 Group:Week 0.97 0.759
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 Week 0.98 2.94, 255.76 7.66e-104 * 1.018 3.05, 265.68 6.57e-106
## 2 Group:Week 0.98 5.88, 255.76 5.63e-15 * 1.018 6.11, 265.68 3.09e-15
## p[HF]<.05
## 1 *
## 2 *
get_anova_table(aov2_rs, correction = "GG")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 Group 2.00 87.00 14.666 3.24e-06 * 0.252
## 2 Week 2.94 255.76 480.225 7.66e-104 * 0.847
## 3 Group:Week 5.88 255.76 15.538 5.63e-15 * 0.263
# 4C. UKURAN EFEK
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ------------------------------------------
## Group | 0.25 | [0.12, 1.00]
## Week | 0.85 | [0.82, 1.00]
## Group:Week | 0.26 | [0.18, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
omega_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## --------------------------------------------
## Group | 0.23 | [0.11, 1.00]
## Week | 0.51 | [0.44, 1.00]
## Group:Week | 0.06 | [0.01, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 4D. PLOT INTERAKSI
p_interaksi <- afex_plot(
aov2,
x = "Week",
trace = "Group",
error = "within",
mapping = c("colour", "shape", "linetype")
) +
labs(
y = "Skor FMA-UE",
x = "Waktu",
title = "Interaksi Kelompok x Waktu pada skor FMA-UE"
) +
theme(legend.position = "bottom")
## Warning: Panel(s) show a mixed within-between-design.
## Error bars do not allow comparisons across all means.
## Suppress error bars with: error = "none"
p_interaksi
# 4E. SIMPLE EFFECTS
# Estimated marginal means waktu dalam setiap kelompok
em2 <- emmeans(aov2, ~ Week | Group)
em2
## Group = CRT:
## Week emmean SE df lower.CL upper.CL
## W0 32.8 0.976 87 30.8 34.7
## W4 36.9 0.921 87 35.1 38.8
## W8 40.1 0.895 87 38.3 41.9
## W12 42.2 0.934 87 40.4 44.1
##
## Group = BWS_TCY:
## Week emmean SE df lower.CL upper.CL
## W0 35.6 0.976 87 33.7 37.6
## W4 41.1 0.921 87 39.3 42.9
## W8 46.6 0.895 87 44.8 48.3
## W12 53.1 0.934 87 51.3 55.0
##
## Group = RAT:
## Week emmean SE df lower.CL upper.CL
## W0 36.0 0.976 87 34.1 38.0
## W4 39.4 0.921 87 37.5 41.2
## W8 45.9 0.895 87 44.1 47.7
## W12 50.0 0.934 87 48.1 51.8
##
## Confidence level used: 0.95
# Efek waktu di dalam masing-masing kelompok
joint_tests(aov2, by = "Group")
## Group = CRT:
## model term df1 df2 F.ratio p.value
## Week 3 87 74.947 <0.0001
##
## Group = BWS_TCY:
## model term df1 df2 F.ratio p.value
## Week 3 87 245.246 <0.0001
##
## Group = RAT:
## model term df1 df2 F.ratio p.value
## Week 3 87 177.657 <0.0001
# Efek kelompok pada masing-masing waktu
joint_tests(aov2, by = "Week")
## Week = W0:
## model term df1 df2 F.ratio p.value
## Group 2 87 3.357 0.0394
##
## Week = W4:
## model term df1 df2 F.ratio p.value
## Group 2 87 5.134 0.0078
##
## Week = W8:
## model term df1 df2 F.ratio p.value
## Group 2 87 15.785 <0.0001
##
## Week = W12:
## model term df1 df2 F.ratio p.value
## Group 2 87 35.898 <0.0001
# 4F. POST-HOC DALAM KELOMPOK
# Semua pasangan waktu dalam masing-masing kelompok
pairs(
em2,
adjust = "holm"
)
## Group = CRT:
## contrast estimate SE df t.ratio p.value
## W0 - W4 -4.17 0.641 87 -6.509 <0.0001
## W0 - W8 -7.35 0.649 87 -11.318 <0.0001
## W0 - W12 -9.49 0.700 87 -13.557 <0.0001
## W4 - W8 -3.18 0.667 87 -4.760 <0.0001
## W4 - W12 -5.31 0.625 87 -8.495 <0.0001
## W8 - W12 -2.14 0.696 87 -3.069 0.0029
##
## Group = BWS_TCY:
## contrast estimate SE df t.ratio p.value
## W0 - W4 -5.44 0.641 87 -8.490 <0.0001
## W0 - W8 -10.92 0.649 87 -16.821 <0.0001
## W0 - W12 -17.49 0.700 87 -24.989 <0.0001
## W4 - W8 -5.48 0.667 87 -8.211 <0.0001
## W4 - W12 -12.04 0.625 87 -19.254 <0.0001
## W8 - W12 -6.56 0.696 87 -9.428 <0.0001
##
## Group = RAT:
## contrast estimate SE df t.ratio p.value
## W0 - W4 -3.32 0.641 87 -5.178 <0.0001
## W0 - W8 -9.89 0.649 87 -15.230 <0.0001
## W0 - W12 -13.92 0.700 87 -19.897 <0.0001
## W4 - W8 -6.57 0.667 87 -9.844 <0.0001
## W4 - W12 -10.60 0.625 87 -16.952 <0.0001
## W8 - W12 -4.03 0.696 87 -5.794 <0.0001
##
## P value adjustment: holm method for 6 tests
# Tiap waktu vs baseline W0
contrast(
em2,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## Group = CRT:
## contrast estimate SE df t.ratio p.value
## W4 - W0 4.17 0.641 87 6.509 <0.0001
## W8 - W0 7.35 0.649 87 11.318 <0.0001
## W12 - W0 9.49 0.700 87 13.557 <0.0001
##
## Group = BWS_TCY:
## contrast estimate SE df t.ratio p.value
## W4 - W0 5.44 0.641 87 8.490 <0.0001
## W8 - W0 10.92 0.649 87 16.821 <0.0001
## W12 - W0 17.49 0.700 87 24.989 <0.0001
##
## Group = RAT:
## contrast estimate SE df t.ratio p.value
## W4 - W0 3.32 0.641 87 5.178 <0.0001
## W8 - W0 9.89 0.649 87 15.230 <0.0001
## W12 - W0 13.92 0.700 87 19.897 <0.0001
##
## P value adjustment: holm method for 3 tests
# 4G. POST-HOC ANTARKELOMPOK PADA SETIAP WAKTU
em2b <- emmeans(aov2, ~ Group | Week)
em2b
## Week = W0:
## Group emmean SE df lower.CL upper.CL
## CRT 32.8 0.976 87 30.8 34.7
## BWS_TCY 35.6 0.976 87 33.7 37.6
## RAT 36.0 0.976 87 34.1 38.0
##
## Week = W4:
## Group emmean SE df lower.CL upper.CL
## CRT 36.9 0.921 87 35.1 38.8
## BWS_TCY 41.1 0.921 87 39.3 42.9
## RAT 39.4 0.921 87 37.5 41.2
##
## Week = W8:
## Group emmean SE df lower.CL upper.CL
## CRT 40.1 0.895 87 38.3 41.9
## BWS_TCY 46.6 0.895 87 44.8 48.3
## RAT 45.9 0.895 87 44.1 47.7
##
## Week = W12:
## Group emmean SE df lower.CL upper.CL
## CRT 42.2 0.934 87 40.4 44.1
## BWS_TCY 53.1 0.934 87 51.3 55.0
## RAT 50.0 0.934 87 48.1 51.8
##
## Confidence level used: 0.95
# Tukey
pairs(em2b, adjust = "tukey")
## Week = W0:
## contrast estimate SE df t.ratio p.value
## CRT - BWS_TCY -2.883 1.38 87 -2.089 0.0979
## CRT - RAT -3.273 1.38 87 -2.372 0.0515
## BWS_TCY - RAT -0.390 1.38 87 -0.283 0.9569
##
## Week = W4:
## contrast estimate SE df t.ratio p.value
## CRT - BWS_TCY -4.153 1.30 87 -3.190 0.0056
## CRT - RAT -2.420 1.30 87 -1.859 0.1570
## BWS_TCY - RAT 1.733 1.30 87 1.331 0.3818
##
## Week = W8:
## contrast estimate SE df t.ratio p.value
## CRT - BWS_TCY -6.457 1.27 87 -5.100 <0.0001
## CRT - RAT -5.813 1.27 87 -4.592 <0.0001
## BWS_TCY - RAT 0.643 1.27 87 0.508 0.8676
##
## Week = W12:
## contrast estimate SE df t.ratio p.value
## CRT - BWS_TCY -10.883 1.32 87 -8.238 <0.0001
## CRT - RAT -7.710 1.32 87 -5.836 <0.0001
## BWS_TCY - RAT 3.173 1.32 87 2.402 0.0479
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# Holm sebagai alternatif koreksi
pairs(em2b, adjust = "holm")
## Week = W0:
## contrast estimate SE df t.ratio p.value
## CRT - BWS_TCY -2.883 1.38 87 -2.089 0.0792
## CRT - RAT -3.273 1.38 87 -2.372 0.0597
## BWS_TCY - RAT -0.390 1.38 87 -0.283 0.7781
##
## Week = W4:
## contrast estimate SE df t.ratio p.value
## CRT - BWS_TCY -4.153 1.30 87 -3.190 0.0059
## CRT - RAT -2.420 1.30 87 -1.859 0.1329
## BWS_TCY - RAT 1.733 1.30 87 1.331 0.1866
##
## Week = W8:
## contrast estimate SE df t.ratio p.value
## CRT - BWS_TCY -6.457 1.27 87 -5.100 <0.0001
## CRT - RAT -5.813 1.27 87 -4.592 <0.0001
## BWS_TCY - RAT 0.643 1.27 87 0.508 0.6126
##
## Week = W12:
## contrast estimate SE df t.ratio p.value
## CRT - BWS_TCY -10.883 1.32 87 -8.238 <0.0001
## CRT - RAT -7.710 1.32 87 -5.836 <0.0001
## BWS_TCY - RAT 3.173 1.32 87 2.402 0.0184
##
## P value adjustment: holm method for 3 tests
# 4H. KONTRAS PERUBAHAN BASELINE -> MINGGU 12
em_full <- emmeans(aov2, ~ Week * Group)
em_full
## Week Group emmean SE df lower.CL upper.CL
## W0 CRT 32.8 0.976 87 30.8 34.7
## W4 CRT 36.9 0.921 87 35.1 38.8
## W8 CRT 40.1 0.895 87 38.3 41.9
## W12 CRT 42.2 0.934 87 40.4 44.1
## W0 BWS_TCY 35.6 0.976 87 33.7 37.6
## W4 BWS_TCY 41.1 0.921 87 39.3 42.9
## W8 BWS_TCY 46.6 0.895 87 44.8 48.3
## W12 BWS_TCY 53.1 0.934 87 51.3 55.0
## W0 RAT 36.0 0.976 87 34.1 38.0
## W4 RAT 39.4 0.921 87 37.5 41.2
## W8 RAT 45.9 0.895 87 44.1 47.7
## W12 RAT 50.0 0.934 87 48.1 51.8
##
## Confidence level used: 0.95
# Perbedaan perubahan W12 - W0 antar kelompok
contrast(
em_full,
interaction = list(
Week = list("W12-W0" = c(-1, 0, 0, 1)),
Group = "pairwise"
),
adjust = "holm"
)
## Week_custom Group_pairwise estimate SE df t.ratio p.value
## W12-W0 CRT - BWS_TCY -8.00 0.99 87 -8.084 <0.0001
## W12-W0 CRT - RAT -4.44 0.99 87 -4.483 <0.0001
## W12-W0 BWS_TCY - RAT 3.56 0.99 87 3.601 0.0005
##
## P value adjustment: holm method for 3 tests
# 4I. TREN LINEAR ANTARKELOMPOK
tren_int <- summary(
contrast(
em_full,
interaction = c(
Week = "poly",
Group = "pairwise"
),
adjust = "none"
)
)
tren_int
## Week_poly Group_pairwise estimate SE df t.ratio p.value
## linear CRT - BWS_TCY -26.303 3.03 87 -8.668 <0.0001
## quadratic CRT - BWS_TCY -3.157 1.24 87 -2.538 0.0129
## cubic CRT - BWS_TCY -1.090 3.08 87 -0.354 0.7244
## linear CRT - RAT -16.703 3.03 87 -5.505 <0.0001
## quadratic CRT - RAT -2.750 1.24 87 -2.211 0.0297
## cubic CRT - RAT 5.743 3.08 87 1.864 0.0657
## linear BWS_TCY - RAT 9.600 3.03 87 3.164 0.0021
## quadratic BWS_TCY - RAT 0.407 1.24 87 0.327 0.7445
## cubic BWS_TCY - RAT 6.833 3.08 87 2.218 0.0292
# Nama kolom dapat berbeda menurut versi emmeans.
# Cek names(tren_int) terlebih dahulu.
names(tren_int)
## [1] "Week_poly" "Group_pairwise" "estimate" "SE"
## [5] "df" "t.ratio" "p.value"
# Biasanya komponen tren Week dapat dipisahkan seperti berikut:
if ("Week_poly" %in% names(tren_int)) {
tren_lin <- subset(tren_int, Week_poly == "linear")
} else if ("Week.poly" %in% names(tren_int)) {
tren_lin <- subset(tren_int, Week.poly == "linear")
} else {
tren_lin <- tren_int
}
if ("p.value" %in% names(tren_lin)) {
tren_lin$p_holm <- p.adjust(tren_lin$p.value, method = "holm")
}
tren_lin
## Week_poly Group_pairwise estimate SE df t.ratio p.value
## 1 linear CRT - BWS_TCY -26.30333 3.034478 87 -8.668158 2.146274e-13
## 4 linear CRT - RAT -16.70333 3.034478 87 -5.504516 3.685402e-07
## 7 linear BWS_TCY - RAT 9.60000 3.034478 87 3.163641 2.146560e-03
## p_holm
## 1 6.438821e-13
## 4 7.370804e-07
## 7 2.146560e-03
# 5. LINEAR MIXED MODEL (LMM)
# LMM random intercept
lmm1 <- lmer(
FMA_UE ~ Group * Week + (1 | ID),
data = data_long,
REML = TRUE
)
summary(lmm1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: FMA_UE ~ Group * Week + (1 | ID)
## Data: data_long
##
## REML criterion at convergence: 1924.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.17262 -0.56319 0.00375 0.55887 2.31646
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 19.444 4.410
## Residual 6.609 2.571
## Number of obs: 360, groups: ID, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 41.64861 0.48416 87.00000 86.023 < 2e-16 ***
## Group1 -3.63278 0.68470 87.00000 -5.306 8.45e-07 ***
## Group2 2.46139 0.68470 87.00000 3.595 0.000538 ***
## Week1 -6.83306 0.23469 261.00000 -29.115 < 2e-16 ***
## Week2 -2.52083 0.23469 261.00000 -10.741 < 2e-16 ***
## Week3 2.55472 0.23469 261.00000 10.886 < 2e-16 ***
## Group1:Week1 1.58056 0.33190 261.00000 4.762 3.18e-06 ***
## Group2:Week1 -1.63028 0.33190 261.00000 -4.912 1.59e-06 ***
## Group1:Week2 1.44167 0.33190 261.00000 4.344 2.01e-05 ***
## Group2:Week2 -0.49917 0.33190 261.00000 -1.504 0.133798
## Group1:Week3 -0.45722 0.33190 261.00000 -1.378 0.169509
## Group2:Week3 -0.09472 0.33190 261.00000 -0.285 0.775568
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Group1 Group2 Week1 Week2 Week3 Gr1:W1 Gr2:W1 Gr1:W2
## Group1 0.000
## Group2 0.000 -0.500
## Week1 0.000 0.000 0.000
## Week2 0.000 0.000 0.000 -0.333
## Week3 0.000 0.000 0.000 -0.333 -0.333
## Group1:Wek1 0.000 0.000 0.000 0.000 0.000 0.000
## Group2:Wek1 0.000 0.000 0.000 0.000 0.000 0.000 -0.500
## Group1:Wek2 0.000 0.000 0.000 0.000 0.000 0.000 -0.333 0.167
## Group2:Wek2 0.000 0.000 0.000 0.000 0.000 0.000 0.167 -0.333 -0.500
## Group1:Wek3 0.000 0.000 0.000 0.000 0.000 0.000 -0.333 0.167 -0.333
## Group2:Wek3 0.000 0.000 0.000 0.000 0.000 0.000 0.167 -0.333 0.167
## Gr2:W2 Gr1:W3
## Group1
## Group2
## Week1
## Week2
## Week3
## Group1:Wek1
## Group2:Wek1
## Group1:Wek2
## Group2:Wek2
## Group1:Wek3 0.167
## Group2:Wek3 -0.333 -0.500
# LMM random intercept + random slope
lmm2 <- lmer(
FMA_UE ~ Group * Week + (1 + Week_num | ID),
data = data_long,
REML = TRUE,
control = lmerControl(optimizer = "bobyqa")
)
summary(lmm2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: FMA_UE ~ Group * Week + (1 + Week_num | ID)
## Data: data_long
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 1923.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.18768 -0.58588 0.00148 0.55450 2.17603
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 21.225918 4.60716
## Week_num 0.005561 0.07457 -0.47
## Residual 6.461131 2.54188
## Number of obs: 360, groups: ID, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 41.64861 0.48416 87.00009 86.023 < 2e-16 ***
## Group1 -3.63278 0.68470 87.00009 -5.306 8.45e-07 ***
## Group2 2.46139 0.68470 87.00009 3.595 0.000538 ***
## Week1 -6.83306 0.23679 192.02147 -28.858 < 2e-16 ***
## Week2 -2.52083 0.23257 199.25996 -10.839 < 2e-16 ***
## Week3 2.55472 0.23257 199.25996 10.985 < 2e-16 ***
## Group1:Week1 1.58056 0.33487 192.02147 4.720 4.54e-06 ***
## Group2:Week1 -1.63028 0.33487 192.02147 -4.868 2.34e-06 ***
## Group1:Week2 1.44167 0.32891 199.25996 4.383 1.89e-05 ***
## Group2:Week2 -0.49917 0.32891 199.25996 -1.518 0.130687
## Group1:Week3 -0.45722 0.32891 199.25996 -1.390 0.166042
## Group2:Week3 -0.09472 0.32891 199.25996 -0.288 0.773653
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Group1 Group2 Week1 Week2 Week3 Gr1:W1 Gr2:W1 Gr1:W2
## Group1 0.000
## Group2 0.000 -0.500
## Week1 0.075 0.000 0.000
## Week2 0.025 0.000 0.000 -0.312
## Week3 -0.025 0.000 0.000 -0.339 -0.336
## Group1:Wek1 0.000 0.075 -0.037 0.000 0.000 0.000
## Group2:Wek1 0.000 -0.037 0.075 0.000 0.000 0.000 -0.500
## Group1:Wek2 0.000 0.025 -0.013 0.000 0.000 0.000 -0.312 0.156
## Group2:Wek2 0.000 -0.013 0.025 0.000 0.000 0.000 0.156 -0.312 -0.500
## Group1:Wek3 0.000 -0.025 0.013 0.000 0.000 0.000 -0.339 0.170 -0.336
## Group2:Wek3 0.000 0.013 -0.025 0.000 0.000 0.000 0.170 -0.339 0.168
## Gr2:W2 Gr1:W3
## Group1
## Group2
## Week1
## Week2
## Week3
## Group1:Wek1
## Group2:Wek1
## Group1:Wek2
## Group2:Wek2
## Group1:Wek3 0.168
## Group2:Wek3 -0.336 -0.500
# Bandingkan struktur random effect
anova(lmm1, lmm2, refit = FALSE)
## Data: data_long
## Models:
## lmm1: FMA_UE ~ Group * Week + (1 | ID)
## lmm2: FMA_UE ~ Group * Week + (1 + Week_num | ID)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 1952.3 2006.7 -962.14 1924.3
## lmm2 16 1955.4 2017.5 -961.68 1923.4 0.9242 2 0.63
# Uji efek tetap
anova(lmm2, ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Group 189.5 94.76 2 87.00 14.666 3.244e-06 ***
## Week 8909.3 2969.78 3 185.37 457.530 < 2.2e-16 ***
## Group:Week 581.9 96.99 6 206.40 14.926 3.693e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.746
## Unadjusted ICC: 0.318
# Diagnostik residual LMM
par(mfrow = c(1, 3))
qqnorm(
resid(lmm2),
main = "Q-Q residual"
)
qqline(resid(lmm2))
qqnorm(
ranef(lmm2)$ID[, 1],
main = "Q-Q random intercept"
)
qqline(ranef(lmm2)$ID[, 1])
plot(
fitted(lmm2),
resid(lmm2),
xlab = "Nilai prediksi",
ylab = "Residual",
main = "Residual vs prediksi"
)
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# Shapiro residual
shapiro.test(resid(lmm2))
##
## Shapiro-Wilk normality test
##
## data: resid(lmm2)
## W = 0.99571, p-value = 0.4315
# 6. OPSIONAL: SIMULASI MISSING VALUE UNTUK DEMONSTRASI LMM
# Bagian ini TIDAK mengubah data utama.
# Hanya digunakan bila dosen meminta demonstrasi keunggulan LMM.
set.seed(2026)
data_miss <- data_long
idx_miss <- sample(
which(data_miss$Week != "W0"),
30,
replace = FALSE
)
data_miss$FMA_UE[idx_miss] <- NA
lmm_miss <- lmer(
FMA_UE ~ Group * Week + (1 + Week_num | ID),
data = data_miss,
REML = TRUE,
na.action = na.exclude,
control = lmerControl(optimizer = "bobyqa")
)
anova(lmm_miss, ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Group 193.3 96.67 2 86.94 14.383 4.016e-06 ***
## Week 8194.1 2731.38 3 168.36 404.434 < 2.2e-16 ***
## Group:Week 587.7 97.96 6 186.28 14.487 1.531e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
n_distinct(
data_miss$ID[is.na(data_miss$FMA_UE)]
)
## [1] 28
# 7. EXPORT HASIL
#Statistik deskriptif
write.csv(
deskriptif,
"hasil_01_statistik_deskriptif_FMA_UE.csv",
row.names = FALSE
)
#Normalitas
write.csv(
shapiro_mixed,
"hasil_02_shapiro_Group_x_Week.csv",
row.names = FALSE
)
#Levene
write.csv(
levene_mixed,
"hasil_03_levene_per_Week.csv",
row.names = FALSE
)
#ANOVA utama
hasil_anova <- as.data.frame(
nice(
aov2,
correction = "GG",
es = c("ges", "pes")
)
)
write.csv(
hasil_anova,
"hasil_04_mixed_ANOVA_GG.csv",
row.names = FALSE
)
#Post-hoc antarkelompok
posthoc_group <- as.data.frame(
pairs(em2b, adjust = "holm")
)
write.csv(
posthoc_group,
"hasil_05_posthoc_antar_kelompok_Holm.csv",
row.names = FALSE
)
# Post-hoc waktu dalam kelompok
posthoc_time <- as.data.frame(
pairs(em2, adjust = "holm")
)
write.csv(
posthoc_time,
"hasil_06_posthoc_dalam_kelompok_Holm.csv",
row.names = FALSE
)
#EMM
emm_group_time <- as.data.frame(em2b)
write.csv(
emm_group_time,
"hasil_07_estimated_marginal_means.csv",
row.names = FALSE
)
#Gambar
ggsave(
"grafik_01_profile_FMA_UE.png",
p_profil,
width = 8,
height = 5.5,
dpi = 300
)
ggsave(
"grafik_02_spaghetti_FMA_UE.png",
p_spag,
width = 9,
height = 6,
dpi = 300
)
ggsave(
"grafik_03_interaksi_Group_x_Week.png",
p_interaksi,
width = 8,
height = 5.5,
dpi = 300
)
# 8. RINGKASAN AKHIR DI CONSOLE
cat("\n============================================================\n")
##
## ============================================================
cat("ANALISIS REPEATED MEASURES FMA-UE SELESAI\n")
## ANALISIS REPEATED MEASURES FMA-UE SELESAI
cat("============================================================\n")
## ============================================================
cat("Jumlah responden :", n_distinct(data_long$ID), "\n")
## Jumlah responden : 90
cat("Kelompok :", n_distinct(data_long$Group), "\n")
## Kelompok : 3
cat("Waktu pengukuran :", n_distinct(data_long$Week), "\n")
## Waktu pengukuran : 4
cat("Outcome : FMA_UE\n")
## Outcome : FMA_UE
cat("\nSemua file hasil tersimpan di:\n")
##
## Semua file hasil tersimpan di:
cat(getwd(), "\n")
## C:/Users/kathe_5vwnplr/Downloads/DAFTAR S2/SEMESTER 1/BIOSTATISTIK/R STUDIO-REPEATED MEASUREMENT ANOVA
cat("============================================================\n")
## ============================================================
# =============================================================================
# SELESAI
# =============================================================================