# =============================================================================
# NAMA : ANDINA OKTAVIANTY
# NIM : 2611018012
# DOSEN PENGAMPU : DR.M.FATHURAHMAN,S.SI.,M.SI
# MATA KULIAH : BIOSTATISTIKA INTERMEDIATE
# =============================================================================
# =============================================================================
# ISI:
# 0. Paket & pengaturan
# 1. Simulasi data (format panjang & lebar)
# 2. Eksplorasi data: statistik deskriptif, profile plot, spaghetti plot
# 3. Repeated Measure ANOVA satu arah (within-subject: Waktu)
# 3a. Uji asumsi: outlier, normalitas, sfierisitas (Mauchly)
# 3b. ANOVA + koreksi Greenhouse-Geisser / Huynh-Feldt
# 3c. Pendekatan multivariat (MANOVA) sebagai pembanding
# 3d. Post hoc berpasangan & kontras polinomial (tren)
# 3e. Alternatif nonparametrik: uji Friedman
# 4. Mixed Design ANOVA (between: Kelompok x within: Waktu)
# 4a. Uji asumsi: outlier, normalitas, Levene, Box's M, Mauchly
# 4b. ANOVA campuran + ukuran efek
# 4c. Analisis efek sederhana (simple effects) & post hoc
# 4d. Kontras interaksi (perubahan dari baseline antarkelompok)
# 5. Pembanding: Linear Mixed Model (LMM)
# 6. Menyimpan data & ringkasan hasil
# =============================================================================
# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# KASUS: PROGRAM PENATALAKSANAAN OBESITAS DAN INDEKS MASSA TUBUH (IMT/BMI)
# DATA SIMULASI
#
# Skenario:
# 120 orang dewasa dengan obesitas diacak menjadi 3 kelompok (n = 40/kelompok):
# - Kontrol : edukasi standar
# - Diet : intervensi diet rendah kalori
# - Diet+Aktif : intervensi diet rendah kalori + aktivitas fisik
#
# IMT (kg/m^2) diukur pada minggu ke-0, 4, 8, dan 12.
#
# Pertanyaan analisis:
# 1. Apakah IMT berubah dari waktu ke waktu?
# 2. Apakah pola perubahan IMT berbeda antar kelompok?
#
# Struktur analisis sengaja dibuat sama seperti latihan sebelumnya:
# - Simulasi data wide & long
# - Statistik deskriptif
# - Profile plot & spaghetti plot
# - Repeated Measure ANOVA satu arah
# - Outlier, Shapiro-Wilk, Mauchly, Greenhouse-Geisser/Huynh-Feldt
# - Pendekatan multivariat (MANOVA)
# - Post hoc & kontras tren
# - Alternatif nonparametrik Friedman
# - Mixed Design ANOVA
# - Simple effects & post hoc
# - Linear Mixed Model (LMM)
# - Simulasi data hilang
# - Menyimpan data
# =============================================================================
# 0. PAKET & PENGATURAN
# Jalankan sekali bila paket belum terpasang:
# 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)
})
options(contrasts = c("contr.sum", "contr.poly"))
afex_options(emmeans_model = "multivariate")
theme_set(theme_bw(base_size = 12))
# 1. SIMULASI DATA
set.seed(2026)
n_per <- 40
kel_lab <- c("Kontrol", "Diet", "Diet+Aktif")
minggu <- c(0, 4, 8, 12)
# Rerata populasi IMT (kg/m^2) per kelompok x waktu
# Kontrol : perubahan kecil
# Diet : penurunan sedang
# Diet+Aktif : penurunan lebih besar
mu <- rbind(
"Kontrol" = c(31.5, 31.3, 31.1, 30.9),
"Diet" = c(31.5, 30.6, 29.8, 29.1),
"Diet+Aktif" = c(31.5, 30.3, 29.2, 28.3)
)
# Komponen variasi antar individu dan galat pengukuran
sd_int <- 1.7
sd_slope <- 0.07
sd_eps <- 0.45
dat_wide <- lapply(seq_along(kel_lab), function(g) {
id <- (g - 1) * n_per + seq_len(n_per)
# Perbedaan IMT awal antar peserta
b0 <- rnorm(n_per, 0, sd_int)
# Perbedaan laju perubahan antar peserta
b1 <- rnorm(n_per, 0, sd_slope)
y <- sapply(seq_along(minggu), function(t)
mu[g, t] +
b0 +
b1 * minggu[t] +
rnorm(n_per, 0, sd_eps)
)
colnames(y) <- paste0("IMT_M", minggu)
data.frame(
id = sprintf("O%03d", id),
kelompok = kel_lab[g],
usia = round(runif(n_per, 25, 60)),
jk = sample(
c("L", "P"),
n_per,
replace = TRUE,
prob = c(.40, .60)
),
round(y, 2)
)
}) |> bind_rows()
dat_wide$kelompok <- factor(dat_wide$kelompok, levels = kel_lab)
dat_wide$id <- factor(dat_wide$id)
# Format panjang
dat_long <- dat_wide |>
pivot_longer(
starts_with("IMT_M"),
names_to = "waktu",
values_to = "imt"
) |>
mutate(
waktu = factor(
waktu,
levels = paste0("IMT_M", minggu),
labels = paste0("M", minggu)
),
minggu = as.numeric(sub("M", "", waktu))
)
head(dat_wide)
## id kelompok usia jk IMT_M0 IMT_M4 IMT_M8 IMT_M12
## 1 O001 Kontrol 46 L 31.86 32.72 31.97 30.62
## 2 O002 Kontrol 43 L 30.34 30.64 29.31 29.92
## 3 O003 Kontrol 44 L 31.22 31.30 30.27 30.02
## 4 O004 Kontrol 40 L 30.66 31.65 31.48 30.18
## 5 O005 Kontrol 56 P 30.34 29.15 29.32 28.94
## 6 O006 Kontrol 55 L 27.60 27.63 27.52 28.61
head(dat_long)
## # A tibble: 6 × 7
## id kelompok usia jk waktu imt minggu
## <fct> <fct> <dbl> <chr> <fct> <dbl> <dbl>
## 1 O001 Kontrol 46 L M0 31.9 0
## 2 O001 Kontrol 46 L M4 32.7 4
## 3 O001 Kontrol 46 L M8 32.0 8
## 4 O001 Kontrol 46 L M12 30.6 12
## 5 O002 Kontrol 43 L M0 30.3 0
## 6 O002 Kontrol 43 L M4 30.6 4
str(dat_long)
## tibble [480 × 7] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 120 levels "O001","O002",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok: Factor w/ 3 levels "Kontrol","Diet",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num [1:480] 46 46 46 46 43 43 43 43 44 44 ...
## $ jk : chr [1:480] "L" "L" "L" "L" ...
## $ waktu : Factor w/ 4 levels "M0","M4","M8",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ imt : num [1:480] 31.9 32.7 32 30.6 30.3 ...
## $ minggu : num [1:480] 0 4 8 12 0 4 8 12 0 4 ...
# 2. EKSPLORASI DATA
# Statistik deskriptif: mean dan SD
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(imt, type = "mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 imt 40 31.4 1.72
## 2 Kontrol M4 imt 40 31.3 1.56
## 3 Kontrol M8 imt 40 31.0 1.73
## 4 Kontrol M12 imt 40 30.8 1.67
## 5 Diet M0 imt 40 31.3 1.94
## 6 Diet M4 imt 40 30.5 1.91
## 7 Diet M8 imt 40 29.8 1.94
## 8 Diet M12 imt 40 29.1 1.98
## 9 Diet+Aktif M0 imt 40 31.7 1.61
## 10 Diet+Aktif M4 imt 40 30.6 1.70
## 11 Diet+Aktif M8 imt 40 29.5 1.70
## 12 Diet+Aktif M12 imt 40 28.5 1.93
# Matriks kovarians
S <- cov(dat_wide[, paste0("IMT_M", minggu)])
round(S, 2)
## IMT_M0 IMT_M4 IMT_M8 IMT_M12
## IMT_M0 3.09 2.76 2.70 2.76
## IMT_M4 2.76 3.07 3.07 3.23
## IMT_M8 2.70 3.07 3.57 3.68
## IMT_M12 2.76 3.23 3.68 4.34
# Matriks korelasi
R <- cor(dat_wide[, paste0("IMT_M", minggu)])
round(R, 2)
## IMT_M0 IMT_M4 IMT_M8 IMT_M12
## IMT_M0 1.00 0.90 0.81 0.75
## IMT_M4 0.90 1.00 0.93 0.88
## IMT_M8 0.81 0.93 1.00 0.93
## IMT_M12 0.75 0.88 0.93 1.00
# Varians selisih antar pasangan waktu
pasangan <- combn(paste0("IMT_M", minggu), 2)
var_selisih <- apply(
pasangan,
2,
function(p)
var(dat_wide[[p[1]]] - dat_wide[[p[2]]])
)
names(var_selisih) <- apply(
pasangan,
2,
paste,
collapse = " - "
)
round(var_selisih, 2)
## IMT_M0 - IMT_M4 IMT_M0 - IMT_M8 IMT_M0 - IMT_M12 IMT_M4 - IMT_M8
## 0.64 1.25 1.90 0.50
## IMT_M4 - IMT_M12 IMT_M8 - IMT_M12
## 0.96 0.56
# 2a. PROFILE PLOT
p_profil <- ggplot(
dat_long,
aes(
minggu,
imt,
colour = kelompok,
group = kelompok
)
) +
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 = .6
) +
scale_x_continuous(
breaks = minggu
) +
labs(
x = "Minggu ke-",
y = "Indeks Massa Tubuh (kg/m²)",
colour = "Kelompok",
title = "Profil rerata IMT (± 95% CI)"
) +
theme(
legend.position = "bottom"
)
p_profil
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
# 2b. SPAGHETTI PLOT
p_spag <- ggplot(
dat_long,
aes(
minggu,
imt,
group = id
)
) +
geom_line(alpha = .3) +
stat_summary(
aes(group = kelompok),
fun = mean,
geom = "line",
colour = "firebrick",
linewidth = 1.2
) +
facet_wrap(~ kelompok) +
scale_x_continuous(
breaks = minggu
) +
labs(
x = "Minggu ke-",
y = "IMT (kg/m²)",
title = "Lintasan individu dan rerata kelompok"
)
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan:
# Apakah IMT berubah selama 12 minggu pada kelompok Diet+Aktif?
d1 <- droplevels(
filter(dat_long, kelompok == "Diet+Aktif")
)
d1w <- filter(
dat_wide,
kelompok == "Diet+Aktif"
)
# 3a. UJI ASUMSI
# (i) Outlier per waktu
d1 |>
group_by(waktu) |>
identify_outliers(imt)
## # A tibble: 5 × 9
## waktu id kelompok usia jk imt minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 M0 O105 Diet+Aktif 60 P 36.0 0 TRUE FALSE
## 2 M4 O105 Diet+Aktif 60 P 35.4 4 TRUE FALSE
## 3 M8 O105 Diet+Aktif 60 P 34.4 8 TRUE FALSE
## 4 M12 O091 Diet+Aktif 37 P 33.1 12 TRUE FALSE
## 5 M12 O105 Diet+Aktif 60 P 33.6 12 TRUE FALSE
# (ii) Normalitas per waktu
d1 |>
group_by(waktu) |>
shapiro_test(imt)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 imt 0.978 0.599
## 2 M4 imt 0.962 0.196
## 3 M8 imt 0.966 0.264
## 4 M12 imt 0.970 0.371
# Q-Q plot
ggpubr::ggqqplot(
d1,
"imt",
facet.by = "waktu"
)
# (iii) Sphericity: Mauchly
aov1_rs <- anova_test(
data = d1,
dv = imt,
wid = id,
within = waktu,
effect.size = "pes"
)
aov1_rs
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 117 349.407 3.31e-58 * 0.9
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.71 0.024 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.8 2.4, 93.59 4.48e-47 * 0.856 2.57, 100.15 3.39e-50
## 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 waktu 2.4 93.59 349.407 4.48e-47 * 0.9
# 3b. REPEATED MEASURE ANOVA DENGAN AFEX
aov1 <- aov_ez(
id = "id",
dv = "imt",
data = d1,
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov1
## Anova Table (Type 3 tests)
##
## Response: imt
## Effect df MSE F ges pes p.value
## 1 waktu 2.40, 93.59 0.27 349.41 *** .327 .900 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov1)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 144717 1 447.04 39 12625.30 < 2.2e-16 ***
## waktu 229 3 25.58 117 349.41 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.71039 0.024391
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.79991 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8560008 3.389893e-50
# Ukuran efek
eta_squared(
aov1,
partial = TRUE
)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.90 | [0.87, 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
## -------------------------------------------
## waktu | 0.32 | [0.20, 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.99692 12625.3 1 39 < 2.2e-16 ***
## waktu 1 0.94116 197.3 3 37 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3d. POST HOC & KONTRAS TREN
em1 <- emmeans(
aov1,
~ waktu
)
em1
## waktu emmean SE df lower.CL upper.CL
## M0 31.7 0.255 39 31.2 32.2
## M4 30.6 0.268 39 30.1 31.2
## M8 29.5 0.269 39 29.0 30.0
## M12 28.5 0.306 39 27.9 29.1
##
## Confidence level used: 0.95
# Semua pasangan waktu
pairs(
em1,
adjust = "bonferroni"
)
## contrast estimate SE df t.ratio p.value
## M0 - M4 1.07 0.0868 39 12.378 <0.0001
## M0 - M8 2.19 0.1070 39 20.429 <0.0001
## M0 - M12 3.20 0.1300 39 24.514 <0.0001
## M4 - M8 1.11 0.0901 39 12.336 <0.0001
## M4 - M12 2.12 0.1130 39 18.750 <0.0001
## M8 - M12 1.01 0.0931 39 10.864 <0.0001
##
## P value adjustment: bonferroni method for 6 tests
# Tiap waktu dibandingkan baseline M0
contrast(
em1,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## contrast estimate SE df t.ratio p.value
## M4 - M0 -1.07 0.0868 39 -12.378 <0.0001
## M8 - M0 -2.19 0.1070 39 -20.429 <0.0001
## M12 - M0 -3.20 0.1300 39 -24.514 <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 -10.7037 0.431 39 -24.855 <0.0001
## quadratic 0.0633 0.124 39 0.511 0.6123
## cubic 0.1388 0.257 39 0.540 0.5919
# 3e. ALTERNATIF NONPARAMETRIK
friedman_test(
d1,
imt ~ waktu | id
)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 imt 40 115. 3 7.86e-25 Friedman test
friedman_effsize(
d1,
imt ~ waktu | id
)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 imt 40 0.961 Kendall W large
d1 |>
wilcox_test(
imt ~ waktu,
paired = TRUE,
p.adjust.method = "bonferroni"
)
## # 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 imt M0 M4 40 40 820 1.82e-12 1.09e-11 ****
## 2 imt M0 M8 40 40 820 1.82e-12 1.09e-11 ****
## 3 imt M0 M12 40 40 820 1.82e-12 1.09e-11 ****
## 4 imt M4 M8 40 40 817 9.09e-12 5.46e-11 ****
## 5 imt M4 M12 40 40 820 1.82e-12 1.09e-11 ****
## 6 imt M8 M12 40 40 816 1.27e-11 7.64e-11 ****
# 4. MIXED DESIGN ANOVA
# BETWEEN = KELOMPOK
# WITHIN = WAKTU
# Pertanyaan:
# Apakah pola perubahan IMT berbeda antar kelompok?
# 4a. UJI ASUMSI
# (i) Outlier per sel
dat_long |>
group_by(kelompok, waktu) |>
identify_outliers(imt)
## # A tibble: 15 × 9
## kelompok waktu id usia jk imt minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol M0 O015 30 L 26.4 0 TRUE FALSE
## 2 Kontrol M4 O006 55 L 27.6 4 TRUE FALSE
## 3 Kontrol M4 O015 30 L 27.0 4 TRUE FALSE
## 4 Kontrol M8 O006 55 L 27.5 8 TRUE FALSE
## 5 Kontrol M8 O015 30 L 26.4 8 TRUE FALSE
## 6 Kontrol M8 O026 27 P 34.5 8 TRUE FALSE
## 7 Kontrol M12 O015 30 L 25.9 12 TRUE FALSE
## 8 Diet M4 O076 41 P 25.6 4 TRUE FALSE
## 9 Diet M8 O075 60 P 24.7 8 TRUE FALSE
## 10 Diet M12 O076 41 P 24.5 12 TRUE FALSE
## 11 Diet+Aktif M0 O105 60 P 36.0 0 TRUE FALSE
## 12 Diet+Aktif M4 O105 60 P 35.4 4 TRUE FALSE
## 13 Diet+Aktif M8 O105 60 P 34.4 8 TRUE FALSE
## 14 Diet+Aktif M12 O091 37 P 33.1 12 TRUE FALSE
## 15 Diet+Aktif M12 O105 60 P 33.6 12 TRUE FALSE
# (ii) Normalitas per sel
dat_long |>
group_by(kelompok, waktu) |>
shapiro_test(imt)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 imt 0.975 0.505
## 2 Kontrol M4 imt 0.967 0.285
## 3 Kontrol M8 imt 0.979 0.640
## 4 Kontrol M12 imt 0.981 0.714
## 5 Diet M0 imt 0.961 0.179
## 6 Diet M4 imt 0.967 0.278
## 7 Diet M8 imt 0.963 0.218
## 8 Diet M12 imt 0.979 0.645
## 9 Diet+Aktif M0 imt 0.978 0.599
## 10 Diet+Aktif M4 imt 0.962 0.196
## 11 Diet+Aktif M8 imt 0.966 0.264
## 12 Diet+Aktif M12 imt 0.970 0.371
# Q-Q plot per kelompok dan waktu
ggpubr::ggqqplot(
dat_long,
"imt",
ggtheme = theme_bw()
) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antarkelompok pada tiap waktu
dat_long |>
group_by(waktu) |>
levene_test(imt ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 117 0.898 0.410
## 2 M4 2 117 0.753 0.473
## 3 M8 2 117 0.466 0.629
## 4 M12 2 117 0.687 0.505
# (iv) Homogenitas matriks kovarians
box_m(
dat_wide[, paste0("IMT_M", minggu)],
dat_wide$kelompok
)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 33.1 0.0329 20 Box's M-test for Homogeneity of Covariance Matric…
# (v) Sphericity diperiksa dari summary model
# 4b. MIXED DESIGN ANOVA
aov2 <- aov_ez(
id = "id",
dv = "imt",
data = dat_long,
between = "kelompok",
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov2
## Anova Table (Type 3 tests)
##
## Response: imt
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 117 11.96 4.28 * .064 .068 .016
## 2 waktu 2.61, 305.91 0.32 322.21 *** .153 .734 <.001
## 3 kelompok:waktu 5.23, 305.91 0.32 44.35 *** .047 .431 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 445364 1 1398.75 117 37252.9424 < 2e-16 ***
## kelompok 102 2 1398.75 117 4.2837 0.01602 *
## waktu 271 3 98.36 351 322.2149 < 2e-16 ***
## kelompok:waktu 75 6 98.36 351 44.3468 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.80121 0.00010448
## kelompok:waktu 0.80121 0.00010448
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.87153 < 2.2e-16 ***
## kelompok:waktu 0.87153 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8932393 4.915785e-90
## kelompok:waktu 0.8932393 3.101891e-36
aov2$Anova
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99687 37253 1 117 < 2.2e-16 ***
## kelompok 2 0.06823 4 2 117 0.01602 *
## waktu 1 0.84888 215 3 115 < 2.2e-16 ***
## kelompok:waktu 2 0.61410 17 6 232 2.247e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix
aov2_rs <- anova_test(
data = dat_long,
dv = imt,
wid = id,
between = kelompok,
within = waktu,
effect.size = "pes",
type = 3
)
get_anova_table(
aov2_rs,
correction = "GG"
)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 kelompok 2.00 117.00 4.284 1.60e-02 * 0.068
## 2 waktu 2.61 305.91 322.215 6.43e-88 * 0.734
## 3 kelompok:waktu 5.23 305.91 44.347 2.04e-35 * 0.431
# Ukuran efek
eta_squared(
aov2,
partial = TRUE
)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.07 | [0.01, 1.00]
## waktu | 0.73 | [0.70, 1.00]
## kelompok:waktu | 0.43 | [0.36, 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
## ------------------------------------------------
## kelompok | 0.05 | [0.00, 1.00]
## waktu | 0.15 | [0.09, 1.00]
## kelompok:waktu | 0.05 | [0.01, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 4c. PLOT INTERAKSI
afex_plot(
aov2,
x = "waktu",
trace = "kelompok",
error = "within",
mapping = c("colour", "shape", "linetype")
) +
labs(
y = "IMT (kg/m²)",
x = "Waktu"
)
## Warning: Panel(s) show a mixed within-between-design.
## Error bars do not allow comparisons across all means.
## Suppress error bars with: error = "none"
# 4d. EFEK SEDERHANA
em2 <- emmeans(
aov2,
~ waktu | kelompok
)
# Efek WAKTU dalam tiap kelompok
joint_tests(
aov2,
by = "kelompok"
)
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## kelompok = Kontrol:
## model term df1 df2 F.ratio p.value
## waktu 3 117 8.957 <0.0001
##
## kelompok = Diet:
## model term df1 df2 F.ratio p.value
## waktu 3 117 85.698 <0.0001
##
## kelompok = Diet+Aktif:
## model term df1 df2 F.ratio p.value
## waktu 3 117 184.390 <0.0001
# Efek KELOMPOK pada tiap waktu
joint_tests(
aov2,
by = "waktu"
)
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 117 0.544 0.5819
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 117 2.523 0.0846
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 117 7.617 0.0008
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 117 15.746 <0.0001
# Post hoc: setiap waktu vs baseline
contrast(
em2,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -0.0432 0.107 117 -0.402 0.6881
## M8 - M0 -0.3730 0.133 117 -2.810 0.0116
## M12 - M0 -0.6110 0.139 117 -4.400 <0.0001
##
## kelompok = Diet:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -0.7718 0.107 117 -7.180 <0.0001
## M8 - M0 -1.4690 0.133 117 -11.067 <0.0001
## M12 - M0 -2.1963 0.139 117 -15.815 <0.0001
##
## kelompok = Diet+Aktif:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -1.0742 0.107 117 -9.994 <0.0001
## M8 - M0 -2.1862 0.133 117 -16.471 <0.0001
## M12 - M0 -3.1972 0.139 117 -23.023 <0.0001
##
## P value adjustment: holm method for 3 tests
# Perbandingan antarkelompok pada tiap waktu
em2b <- emmeans(
aov2,
~ kelompok | waktu
)
pairs(
em2b,
adjust = "tukey"
)
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 0.0605 0.395 117 0.153 0.9871
## Kontrol - (Diet+Aktif) -0.3222 0.395 117 -0.817 0.6934
## Diet - (Diet+Aktif) -0.3827 0.395 117 -0.970 0.5972
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 0.7890 0.387 117 2.041 0.1071
## Kontrol - (Diet+Aktif) 0.7087 0.387 117 1.833 0.1633
## Diet - (Diet+Aktif) -0.0803 0.387 117 -0.208 0.9765
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 1.1565 0.401 117 2.885 0.0129
## Kontrol - (Diet+Aktif) 1.4910 0.401 117 3.719 0.0009
## Diet - (Diet+Aktif) 0.3345 0.401 117 0.834 0.6825
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - Diet 1.6458 0.417 117 3.946 0.0004
## Kontrol - (Diet+Aktif) 2.2640 0.417 117 5.429 <0.0001
## Diet - (Diet+Aktif) 0.6182 0.417 117 1.482 0.3031
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4e. KONTRAS INTERAKSI
# Apakah perubahan IMT M12 - M0 berbeda antar kelompok?
em_full <- emmeans(
aov2,
~ waktu * kelompok
)
contrast(
em_full,
interaction = list(
waktu = list(
"M12-M0" = c(-1, 0, 0, 1)
),
kelompok = "pairwise"
),
adjust = "holm"
)
## waktu_custom kelompok_pairwise estimate SE df t.ratio p.value
## M12-M0 Kontrol - Diet 1.59 0.196 117 8.072 <0.0001
## M12-M0 Kontrol - (Diet+Aktif) 2.59 0.196 117 13.168 <0.0001
## M12-M0 Diet - (Diet+Aktif) 1.00 0.196 117 5.097 <0.0001
##
## P value adjustment: holm method for 3 tests
# 4f. TREN LINEAR ANTARKELOMPOK
contrast(
em2,
"poly"
)[c(1, 4, 7)]
## contrast kelompok estimate SE df t.ratio p.value
## linear Kontrol -2.16 0.457 117 -4.735 <0.0001
## linear Diet -7.29 0.457 117 -15.952 <0.0001
## linear Diet+Aktif -10.70 0.457 117 -23.434 <0.0001
tren_int <- summary(
contrast(
em_full,
interaction = c(
waktu = "poly",
kelompok = "pairwise"
),
adjust = "none"
)
)
tren_lin <- subset(
tren_int,
waktu_poly == "linear"
)
tren_lin$p.holm <- p.adjust(
tren_lin$p.value,
"holm"
)
tren_lin
## waktu_poly kelompok_pairwise estimate SE df t.ratio
## 1 linear Kontrol - Diet 5.12325 0.6459466 117 7.931383
## 4 linear Kontrol - (Diet+Aktif) 8.54100 0.6459466 117 13.222455
## 7 linear Diet - (Diet+Aktif) 3.41775 0.6459466 117 5.291072
## p.value p.holm
## 1 1.435931e-12 2.871862e-12
## 4 5.678700e-25 1.703610e-24
## 7 5.748809e-07 5.748809e-07
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
lmm1 <- lmer(
imt ~ kelompok * waktu + (1 | id),
data = dat_long,
REML = TRUE
)
lmm2 <- lmer(
imt ~ kelompok * waktu + (1 + minggu | id),
data = dat_long,
REML = TRUE
)
# Perbandingan struktur efek acak
anova(
lmm1,
lmm2,
refit = FALSE
)
## Data: dat_long
## Models:
## lmm1: imt ~ kelompok * waktu + (1 | id)
## lmm2: imt ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 1261.3 1319.7 -616.64 1233.3
## lmm2 16 1244.0 1310.7 -605.98 1212.0 21.316 2 2.351e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 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)
## kelompok 1.814 0.907 2 117.00 4.2836 0.01602 *
## waktu 137.527 45.842 3 249.65 215.7770 < 2e-16 ***
## kelompok:waktu 38.156 6.359 6 278.39 29.9079 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Intraclass correlation
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.912
## Unadjusted ICC: 0.706
# 5a. 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 intersep acak"
)
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))
# 5b. SIMULASI DATA HILANG
set.seed(1)
dat_miss <- dat_long
dat_miss$imt[
sample(
which(dat_miss$waktu != "M0"),
30
)
] <- NA
lmm_miss <- lmer(
imt ~ kelompok * waktu + (1 + minggu | id),
data = dat_miss,
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)
## kelompok 1.842 0.921 2 116.99 4.2967 0.01582 *
## waktu 136.381 45.460 3 232.06 211.3945 < 2e-16 ***
## kelompok:waktu 37.861 6.310 6 257.42 29.3180 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Jumlah ID yang memiliki minimal satu data hilang
n_distinct(
dat_miss$id[
is.na(dat_miss$imt)
]
)
## [1] 27
# 6. SIMPAN DATA
write.csv(
dat_wide,
"data_imt_obesitas_wide.csv",
row.names = FALSE
)
write.csv(
dat_long,
"data_imt_obesitas_long.csv",
row.names = FALSE
)
# 7. SESSION INFO
sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_Indonesia.utf8 LC_CTYPE=English_Indonesia.utf8
## [3] LC_MONETARY=English_Indonesia.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_Indonesia.utf8
##
## time zone: Asia/Makassar
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lmerTest_3.2-1 effectsize_1.0.3 car_3.1-5 carData_3.0-6
## [5] rstatix_1.1.0 emmeans_2.0.4 afex_1.5-1 lme4_2.0-6
## [9] Matrix_1.7-5 ggplot2_4.0.3 tidyr_1.3.2 dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.61 bslib_0.12.0
## [4] bayestestR_0.19.0 insight_1.5.4 lattice_0.22-9
## [7] numDeriv_2016.8-1.1 vctrs_0.7.3 tools_4.6.1
## [10] Rdpack_2.6.6 generics_0.1.4 pbkrtest_0.5.5
## [13] parallel_4.6.1 datawizard_1.4.0 tibble_3.3.1
## [16] pkgconfig_2.0.3 RColorBrewer_1.1-3 S7_0.2.2
## [19] lifecycle_1.0.5 compiler_4.6.1 farver_2.1.2
## [22] stringr_1.6.0 htmltools_0.5.9 sass_0.4.10
## [25] yaml_2.3.12 Formula_1.2-6 ggpubr_1.0.0
## [28] pillar_1.11.1 nloptr_2.2.1 jquerylib_0.1.4
## [31] MASS_7.3-65 cachem_1.1.0 reformulas_0.4.4
## [34] boot_1.3-32 abind_1.4-8 nlme_3.1-169
## [37] tidyselect_1.2.1 digest_0.6.39 performance_0.18.2
## [40] mvtnorm_1.4-2 stringi_1.8.9 reshape2_1.4.5
## [43] purrr_1.2.2 labeling_0.4.3 splines_4.6.1
## [46] fastmap_1.2.0 grid_4.6.1 cli_3.6.6
## [49] magrittr_2.0.5 utf8_1.2.6 broom_1.0.13
## [52] withr_3.0.3 scales_1.4.0 backports_1.5.1
## [55] estimability_2.0.0 rmarkdown_2.32 ggsignif_0.6.4
## [58] evaluate_1.0.5 knitr_1.52 parameters_0.29.3
## [61] rbibutils_2.4.1 rlang_1.3.0 Rcpp_1.1.2
## [64] glue_1.8.1 minqa_1.2.8 jsonlite_2.0.0
## [67] R6_2.6.1 plyr_1.8.9