# Nama : Dhea Fara Dhina
# NIM : 2611018015
#=============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: intervensi gaya hidup dan berat badan
# pada subjek obesitas (DATA SIMULASI/MODIFIKASI)
#
# Dasar modifikasi:
# Bhutani S, Klempel MC, Kroeger CM, Trepanowski JF, Varady KA. (2013).
# Alternate day fasting and endurance exercise combine to reduce body weight
# and favorably alter plasma lipids in obese humans.
# Obesity (Silver Spring), 21(7), 1370-1379.
# DOI: 10.1002/oby.20353
#
# Catatan:
# Data yang digunakan adalah DATA SIMULASI yang dimodifikasi berdasarkan
# desain penelitian dan pola perubahan berat badan dari artikel tersebut.
# Nilai Bulan 1 dan Bulan 2 dibuat untuk kebutuhan latihan repeated measure,
# bukan merupakan data bulanan asli dari artikel.
#
# Isi:
# 0. Paket & pengaturan
# 1. Import data utama CSV + pembentukan data long
# 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 ke bulan 3)
# 5. Pembanding: Linear Mixed Model (LMM)
# 6. Menyimpan data & ringkasan hasil
# =============================================================================
# 0. PAKET & PENGATURAN
# Jalankan sekali bila paket belum terpasang:
# install.packages(c("dplyr", "tidyr", "ggplot2", "afex", "emmeans",
# "rstatix", "car", "effectsize", "lme4", "lmerTest",
# "performance", "ggpubr", "pbkrtest"))
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. IMPORT DATA
# Skenario penelitian. 64 subjek dengan obesitas dibagi ke empat kelompok:
# - Kontrol : kelompok kontrol
# - ADF : alternate day fasting
# - Olahraga : endurance exercise
# - ADF + Olahraga : alternate day fasting + endurance exercise
# Berat badan diukur pada Bulan ke-0 (baseline), 1, 2, dan 3.
file_data <- "DATA_UTAMA_Obesitas_64.csv"
# Jika file tidak berada di working directory, pilih file secara manual:
if (!file.exists(file_data)) {
file_data <- file.choose()
}
dat_wide_raw <- read.csv(
file_data,
stringsAsFactors = FALSE,
check.names = FALSE
)
# Cek data awal
head(dat_wide_raw)
## ID Kelompok Usia JK BB_Bulan_0_kg BB_Bulan_1_kg BB_Bulan_2_kg
## 1 P001 ADF + Olahraga 56 P 83.3 80.6 77.4
## 2 P002 ADF + Olahraga 55 P 87.5 85.4 83.1
## 3 P003 Kontrol 46 P 88.7 88.8 88.5
## 4 P004 Kontrol 57 P 82.1 82.1 82.1
## 5 P005 Kontrol 35 P 95.1 95.0 95.0
## 6 P006 Olahraga 53 P 95.0 94.7 94.2
## BB_Bulan_3_kg Perubahan_BB_3_Bulan
## 1 74.5 -8.8
## 2 81.3 -6.2
## 3 88.5 -0.2
## 4 82.1 0.0
## 5 95.2 0.1
## 6 94.1 -0.9
str(dat_wide_raw)
## 'data.frame': 64 obs. of 9 variables:
## $ ID : chr "P001" "P002" "P003" "P004" ...
## $ Kelompok : chr "ADF + Olahraga" "ADF + Olahraga" "Kontrol" "Kontrol" ...
## $ Usia : int 56 55 46 57 35 53 59 49 53 54 ...
## $ JK : chr "P" "P" "P" "P" ...
## $ BB_Bulan_0_kg : num 83.3 87.5 88.7 82.1 95.1 ...
## $ BB_Bulan_1_kg : num 80.6 85.4 88.8 82.1 95 ...
## $ BB_Bulan_2_kg : num 77.4 83.1 88.5 82.1 95 ...
## $ BB_Bulan_3_kg : num 74.5 81.3 88.5 82.1 95.2 ...
## $ Perubahan_BB_3_Bulan: num -8.8 -6.2 -0.2 0 0.1 -0.9 -3.8 0.3 -2.9 -6.2 ...
dim(dat_wide_raw)
## [1] 64 9
# Menghapus kolom index jika ada
if ("index" %in% names(dat_wide_raw)) {
dat_wide_raw <- dat_wide_raw |> select(-index)
}
names(dat_wide_raw)
## [1] "ID" "Kelompok" "Usia"
## [4] "JK" "BB_Bulan_0_kg" "BB_Bulan_1_kg"
## [7] "BB_Bulan_2_kg" "BB_Bulan_3_kg" "Perubahan_BB_3_Bulan"
# Faktor
kel_lab <- c("Kontrol", "ADF", "Olahraga", "ADF + Olahraga")
bulan <- c(0, 1, 2, 3)
dat_wide <- dat_wide_raw |>
mutate(
id = factor(ID),
kelompok = factor(Kelompok, levels = kel_lab),
usia = Usia,
jk = factor(JK)
) |>
select(
id, kelompok, usia, jk,
BB_Bulan_0_kg, BB_Bulan_1_kg, BB_Bulan_2_kg, BB_Bulan_3_kg,
Perubahan_BB_3_Bulan
)
# Format panjang (satu baris = satu pengukuran)
dat_long <- dat_wide |>
pivot_longer(
cols = starts_with("BB_Bulan_"),
names_to = "waktu",
values_to = "bb"
) |>
mutate(
bulan = as.numeric(
sub("BB_Bulan_", "", sub("_kg", "", waktu))
),
waktu = factor(
waktu,
levels = paste0("BB_Bulan_", bulan, "_kg"),
labels = paste0("M", bulan)
)
) |>
arrange(id, bulan)
head(dat_wide)
## id kelompok usia jk BB_Bulan_0_kg BB_Bulan_1_kg BB_Bulan_2_kg
## 1 P001 ADF + Olahraga 56 P 83.3 80.6 77.4
## 2 P002 ADF + Olahraga 55 P 87.5 85.4 83.1
## 3 P003 Kontrol 46 P 88.7 88.8 88.5
## 4 P004 Kontrol 57 P 82.1 82.1 82.1
## 5 P005 Kontrol 35 P 95.1 95.0 95.0
## 6 P006 Olahraga 53 P 95.0 94.7 94.2
## BB_Bulan_3_kg Perubahan_BB_3_Bulan
## 1 74.5 -8.8
## 2 81.3 -6.2
## 3 88.5 -0.2
## 4 82.1 0.0
## 5 95.2 0.1
## 6 94.1 -0.9
head(dat_long)
## # A tibble: 6 × 8
## id kelompok usia jk Perubahan_BB_3_Bulan waktu bb bulan
## <fct> <fct> <int> <fct> <dbl> <fct> <dbl> <dbl>
## 1 P001 ADF + Olahraga 56 P -8.8 M0 83.3 0
## 2 P001 ADF + Olahraga 56 P -8.8 M1 80.6 1
## 3 P001 ADF + Olahraga 56 P -8.8 M2 77.4 2
## 4 P001 ADF + Olahraga 56 P -8.8 M3 74.5 3
## 5 P002 ADF + Olahraga 55 P -6.2 M0 87.5 0
## 6 P002 ADF + Olahraga 55 P -6.2 M1 85.4 1
str(dat_long)
## tibble [256 × 8] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 64 levels "P001","P002",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok : Factor w/ 4 levels "Kontrol","ADF",..: 4 4 4 4 4 4 4 4 1 1 ...
## $ usia : int [1:256] 56 56 56 56 55 55 55 55 46 46 ...
## $ jk : Factor w/ 2 levels "L","P": 2 2 2 2 2 2 2 2 2 2 ...
## $ Perubahan_BB_3_Bulan: num [1:256] -8.8 -8.8 -8.8 -8.8 -6.2 -6.2 -6.2 -6.2 -0.2 -0.2 ...
## $ waktu : Factor w/ 4 levels "M0","M1","M2",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ bb : num [1:256] 83.3 80.6 77.4 74.5 87.5 85.4 83.1 81.3 88.7 88.8 ...
## $ bulan : num [1:256] 0 1 2 3 0 1 2 3 0 1 ...
# Opsional: simpan hasil format long ke CSV
write.csv(
dat_long,
"DATA_LONG_Obesitas_64.csv",
row.names = FALSE
)
# 2. EKSPLORASI DATA
# Statistik deskriptif per kelompok dan waktu
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(bb, type = "mean_sd")
desk
## # A tibble: 16 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 bb 16 94.2 8.35
## 2 Kontrol M1 bb 16 94.1 8.36
## 3 Kontrol M2 bb 16 94.1 8.43
## 4 Kontrol M3 bb 16 94.1 8.42
## 5 ADF M0 bb 16 94.8 7.40
## 6 ADF M1 bb 16 93.7 7.44
## 7 ADF M2 bb 16 92.8 7.46
## 8 ADF M3 bb 16 91.8 7.50
## 9 Olahraga M0 bb 16 94.1 6.39
## 10 Olahraga M1 bb 16 93.8 6.39
## 11 Olahraga M2 bb 16 93.4 6.38
## 12 Olahraga M3 bb 16 93.2 6.33
## 13 ADF + Olahraga M0 bb 16 91.7 5.80
## 14 ADF + Olahraga M1 bb 16 89.8 5.86
## 15 ADF + Olahraga M2 bb 16 87.4 5.92
## 16 ADF + Olahraga M3 bb 16 85.4 6.05
# Matriks kovarians & korelasi antarwaktu
S <- cov(dat_wide[, paste0("BB_Bulan_", bulan, "_kg")])
R <- cor(dat_wide[, paste0("BB_Bulan_", bulan, "_kg")])
round(S, 2)
## BB_Bulan_0_kg BB_Bulan_1_kg BB_Bulan_2_kg BB_Bulan_3_kg
## BB_Bulan_0_kg 48.85 49.60 50.66 51.39
## BB_Bulan_1_kg 49.60 50.96 52.71 54.06
## BB_Bulan_2_kg 50.66 52.71 55.35 57.44
## BB_Bulan_3_kg 51.39 54.06 57.44 60.20
round(R, 2)
## BB_Bulan_0_kg BB_Bulan_1_kg BB_Bulan_2_kg BB_Bulan_3_kg
## BB_Bulan_0_kg 1.00 0.99 0.97 0.95
## BB_Bulan_1_kg 0.99 1.00 0.99 0.98
## BB_Bulan_2_kg 0.97 0.99 1.00 1.00
## BB_Bulan_3_kg 0.95 0.98 1.00 1.00
# Varians selisih antarpasangan waktu
pasangan <- combn(paste0("BB_Bulan_", bulan, "_kg"), 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)
## BB_Bulan_0_kg - BB_Bulan_1_kg BB_Bulan_0_kg - BB_Bulan_2_kg
## 0.61 2.89
## BB_Bulan_0_kg - BB_Bulan_3_kg BB_Bulan_1_kg - BB_Bulan_2_kg
## 6.27 0.89
## BB_Bulan_1_kg - BB_Bulan_3_kg BB_Bulan_2_kg - BB_Bulan_3_kg
## 3.04 0.67
# Profile plot: rerata +/- 95% CI per kelompok
p_profil <- ggplot(
dat_long,
aes(bulan, bb, 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 = .08
) +
scale_x_continuous(breaks = bulan) +
labs(
x = "Bulan ke-",
y = "Berat badan (kg)",
colour = "Kelompok",
title = "Profil rerata berat badan (± 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.
# Spaghetti plot: lintasan tiap subjek
p_spag <- ggplot(
dat_long,
aes(bulan, bb, 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 = bulan) +
labs(
x = "Bulan ke-",
y = "Berat badan (kg)",
title = "Lintasan individu dan rerata kelompok"
)
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan:
# apakah berat badan berubah selama 3 bulan pada kelompok ADF + Olahraga?
d1 <- droplevels(
filter(dat_long, kelompok == "ADF + Olahraga")
)
d1w <- filter(
dat_wide,
kelompok == "ADF + Olahraga"
)
## 3a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per waktu
d1 |>
group_by(waktu) |>
identify_outliers(bb)
## [1] waktu id kelompok
## [4] usia jk Perubahan_BB_3_Bulan
## [7] bb bulan is.outlier
## [10] is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas per waktu (Shapiro-Wilk)
sw1 <- d1 |>
group_by(waktu) |>
shapiro_test(bb)
sw1
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 bb 0.936 0.298
## 2 M1 bb 0.935 0.296
## 3 M2 bb 0.953 0.544
## 4 M3 bb 0.958 0.630
# Q-Q plot
ggpubr::ggqqplot(
d1,
"bb",
facet.by = "waktu"
)
# (iii) Sfierisitas: Mauchly
aov1_rs <- anova_test(
data = d1,
dv = bb,
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 45 337.53 7.53e-31 * 0.957
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.002 8.48e-17 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.345 1.04, 15.53 4.9e-12 * 0.348 1.04, 15.65 4.13e-12
## 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 1.04 15.53 337.53 4.9e-12 * 0.957
## 3b. ANOVA dengan afex -------------------------------------------------------
aov1 <- aov_ez(
id = "id",
dv = "bb",
data = d1,
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov1
## Anova Table (Type 3 tests)
##
## Response: bb
## Effect df MSE F ges pes p.value
## 1 waktu 1.04, 15.53 1.02 337.53 *** .145 .957 <.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) 502132 1 2078.48 15 3623.79 < 2.2e-16 ***
## waktu 356 3 15.83 45 337.53 < 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.0019453 8.4755e-17
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.34514 4.902e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.3477097 4.132171e-12
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.96 | [0.94, 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.14 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
## 3c. Pendekatan multivariat (MANOVA) ----------------------------------------
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.99588 3623.8 1 15 < 2.2e-16 ***
## waktu 1 0.95986 103.6 3 13 2.502e-09 ***
## ---
## 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 91.7 1.45 15 88.6 94.8
## M1 89.8 1.47 15 86.7 92.9
## M2 87.4 1.48 15 84.3 90.6
## M3 85.4 1.51 15 82.2 88.6
##
## Confidence level used: 0.95
# Semua pasangan waktu
pairs(
em1,
adjust = "bonferroni"
)
## contrast estimate SE df t.ratio p.value
## M0 - M1 1.89 0.111 15 17.061 <0.0001
## M0 - M2 4.24 0.235 15 18.030 <0.0001
## M0 - M3 6.24 0.336 15 18.561 <0.0001
## M1 - M2 2.35 0.133 15 17.722 <0.0001
## M1 - M3 4.36 0.233 15 18.718 <0.0001
## M2 - M3 2.01 0.107 15 18.813 <0.0001
##
## P value adjustment: bonferroni method for 6 tests
# Setiap waktu dibandingkan dengan baseline
contrast(
em1,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## contrast estimate SE df t.ratio p.value
## M1 - M0 -1.89 0.111 15 -17.061 <0.0001
## M2 - M0 -4.24 0.235 15 -18.030 <0.0001
## M3 - M0 -6.24 0.336 15 -18.561 <0.0001
##
## P value adjustment: holm method for 3 tests
# Tren linear, kuadratik, dan kubik
contrast(
em1,
"poly"
)
## contrast estimate SE df t.ratio p.value
## linear -21.081 1.1400 15 -18.536 <0.0001
## quadratic -0.119 0.0476 15 -2.493 0.0248
## cubic 0.806 0.1180 15 6.805 <0.0001
## 3e. Alternatif nonparametrik ----------------------------------------------
friedman_test(
d1,
bb ~ waktu | id
)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 bb 16 48 3 2.13e-10 Friedman test
friedman_effsize(
d1,
bb ~ waktu | id
)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 bb 16 1 Kendall W large
d1 |>
wilcox_test(
bb ~ 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 bb M0 M1 16 16 136 0.0000305 0.000183 ***
## 2 bb M0 M2 16 16 136 0.0000305 0.000183 ***
## 3 bb M0 M3 16 16 136 0.0000305 0.000183 ***
## 4 bb M1 M2 16 16 136 0.0000305 0.000183 ***
## 5 bb M1 M3 16 16 136 0.0000305 0.000183 ***
## 6 bb M2 M3 16 16 136 0.0000305 0.000183 ***
# (Opsional) ANOVA robust berbasis trimmed mean -- paket WRS2
if (requireNamespace("WRS2", quietly = TRUE)) {
print(
WRS2::rmanova(
d1$bb,
d1$waktu,
d1$id,
tr = 0.2
)
)
}
## Call:
## WRS2::rmanova(y = d1$bb, groups = d1$waktu, blocks = d1$id, tr = 0.2)
##
## Test statistic: F = 268.6431
## Degrees of freedom 1: 1.1
## Degrees of freedom 2: 9.87
## p-value: 0
# 4. MIXED DESIGN ANOVA (Kelompok [between] x Waktu [within])
## 4a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per sel
dat_long |>
group_by(kelompok, waktu) |>
identify_outliers(bb)
## # A tibble: 8 × 10
## kelompok waktu id usia jk Perubahan_BB_3_Bulan bb bulan is.outlier
## <fct> <fct> <fct> <int> <fct> <dbl> <dbl> <dbl> <lgl>
## 1 ADF M0 P041 25 P -2.9 111. 0 TRUE
## 2 ADF M0 P060 57 P -3.1 80.4 0 TRUE
## 3 ADF M1 P041 25 P -2.9 110. 1 TRUE
## 4 ADF M1 P060 57 P -3.1 79.2 1 TRUE
## 5 ADF M2 P041 25 P -2.9 109. 2 TRUE
## 6 ADF M2 P060 57 P -3.1 78.3 2 TRUE
## 7 ADF M3 P041 25 P -2.9 108. 3 TRUE
## 8 ADF M3 P060 57 P -3.1 77.3 3 TRUE
## # ℹ 1 more variable: is.extreme <lgl>
# (ii) Normalitas per sel
dat_long |>
group_by(kelompok, waktu) |>
shapiro_test(bb)
## # A tibble: 16 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 bb 0.959 0.636
## 2 Kontrol M1 bb 0.959 0.645
## 3 Kontrol M2 bb 0.957 0.616
## 4 Kontrol M3 bb 0.958 0.617
## 5 ADF M0 bb 0.963 0.717
## 6 ADF M1 bb 0.964 0.728
## 7 ADF M2 bb 0.968 0.811
## 8 ADF M3 bb 0.970 0.833
## 9 Olahraga M0 bb 0.899 0.0788
## 10 Olahraga M1 bb 0.902 0.0862
## 11 Olahraga M2 bb 0.896 0.0699
## 12 Olahraga M3 bb 0.899 0.0778
## 13 ADF + Olahraga M0 bb 0.936 0.298
## 14 ADF + Olahraga M1 bb 0.935 0.296
## 15 ADF + Olahraga M2 bb 0.953 0.544
## 16 ADF + Olahraga M3 bb 0.958 0.630
# Q-Q plot per kelompok dan waktu
ggpubr::ggqqplot(
dat_long,
"bb",
ggtheme = theme_bw()
) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antarkelompok pada tiap waktu
dat_long |>
group_by(waktu) |>
levene_test(bb ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 3 60 0.652 0.585
## 2 M1 3 60 0.607 0.613
## 3 M2 3 60 0.580 0.630
## 4 M3 3 60 0.528 0.665
# (iv) Homogenitas matriks kovarians antarkelompok (Box's M)
box_m(
dat_wide[, paste0("BB_Bulan_", bulan, "_kg")],
dat_wide$kelompok
)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 84.7 0.000000406 30 Box's M-test for Homogeneity of Covariance Ma…
## 4b. ANOVA campuran ---------------------------------------------------------
aov2 <- aov_ez(
id = "id",
dv = "bb",
data = dat_long,
between = "kelompok",
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov2
## Anova Table (Type 3 tests)
##
## Response: bb
## Effect df MSE F ges pes p.value
## 1 kelompok 3, 60 201.17 2.11 .095 .095 .109
## 2 waktu 1.11, 66.39 0.28 756.83 *** .019 .927 <.001
## 3 kelompok:waktu 3.32, 66.39 0.28 219.61 *** .017 .917 <.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) 2185630 1 12070.2 60 10864.5716 <2e-16 ***
## kelompok 1271 3 12070.2 60 2.1053 0.109
## waktu 238 3 18.9 180 756.8253 <2e-16 ***
## kelompok:waktu 207 9 18.9 180 219.6091 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.015709 1.3786e-50
## kelompok:waktu 0.015709 1.3786e-50
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.36883 < 2.2e-16 ***
## kelompok:waktu 0.36883 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.370703 1.991850e-39
## kelompok:waktu 0.370703 1.157069e-35
# Pendekatan multivariat
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.99451 10864.6 1 60 <2e-16 ***
## kelompok 3 0.09524 2.1 3 60 0.109
## waktu 1 0.93162 263.4 3 58 <2e-16 ***
## kelompok:waktu 3 1.43814 18.4 9 180 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix
aov2_rs <- anova_test(
data = dat_long,
dv = bb,
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 3.00 60.00 2.105 1.09e-01 0.095
## 2 waktu 1.11 66.39 756.825 3.06e-39 * 0.927
## 3 kelompok:waktu 3.32 66.39 219.609 1.70e-35 * 0.917
# Ukuran efek
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.10 | [0.00, 1.00]
## waktu | 0.93 | [0.91, 1.00]
## kelompok:waktu | 0.92 | [0.90, 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.02 | [0.00, 1.00]
## kelompok:waktu | 0.02 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Plot interaksi
afex_plot(
aov2,
x = "waktu",
trace = "kelompok",
error = "within",
mapping = c("colour", "shape", "linetype")
) +
labs(
y = "Berat badan (kg)",
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"
## 4c. Efek sederhana & post hoc ---------------------------------------------
em2 <- emmeans(
aov2,
~ waktu | kelompok
)
em2
## kelompok = Kontrol:
## waktu emmean SE df lower.CL upper.CL
## M0 94.2 1.76 60 90.6 97.7
## M1 94.1 1.77 60 90.6 97.7
## M2 94.1 1.78 60 90.6 97.7
## M3 94.1 1.78 60 90.5 97.7
##
## kelompok = ADF:
## waktu emmean SE df lower.CL upper.CL
## M0 94.8 1.76 60 91.3 98.4
## M1 93.7 1.77 60 90.2 97.2
## M2 92.8 1.78 60 89.2 96.4
## M3 91.8 1.78 60 88.2 95.3
##
## kelompok = Olahraga:
## waktu emmean SE df lower.CL upper.CL
## M0 94.1 1.76 60 90.6 97.6
## M1 93.8 1.77 60 90.2 97.3
## M2 93.4 1.78 60 89.9 97.0
## M3 93.2 1.78 60 89.6 96.7
##
## kelompok = ADF + Olahraga:
## waktu emmean SE df lower.CL upper.CL
## M0 91.7 1.76 60 88.1 95.2
## M1 89.8 1.77 60 86.2 93.3
## M2 87.4 1.78 60 83.9 91.0
## M3 85.4 1.78 60 81.9 89.0
##
## Confidence level used: 0.95
# Efek WAKTU di dalam masing-masing 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 60 0.053 0.9836
##
## kelompok = ADF:
## model term df1 df2 F.ratio p.value
## waktu 3 60 118.045 <0.0001
##
## kelompok = Olahraga:
## model term df1 df2 F.ratio p.value
## waktu 3 60 8.813 <0.0001
##
## kelompok = ADF + Olahraga:
## model term df1 df2 F.ratio p.value
## waktu 3 60 398.003 <0.0001
# Efek KELOMPOK pada tiap waktu
joint_tests(
aov2,
by = "waktu"
)
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 3 60 0.621 0.6041
##
## waktu = M1:
## model term df1 df2 F.ratio p.value
## kelompok 3 60 1.347 0.2678
##
## waktu = M2:
## model term df1 df2 F.ratio p.value
## kelompok 3 60 2.954 0.0396
##
## waktu = M3:
## model term df1 df2 F.ratio p.value
## kelompok 3 60 4.800 0.0046
# Setiap waktu dibandingkan dengan baseline
contrast(
em2,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M1 - M0 -0.0187 0.0623 60 -0.301 1.0000
## M2 - M0 -0.0437 0.1270 60 -0.345 1.0000
## M3 - M0 -0.0688 0.1810 60 -0.379 1.0000
##
## kelompok = ADF:
## contrast estimate SE df t.ratio p.value
## M1 - M0 -1.1500 0.0623 60 -18.446 <0.0001
## M2 - M0 -2.0500 0.1270 60 -16.177 <0.0001
## M3 - M0 -3.0688 0.1810 60 -16.938 <0.0001
##
## kelompok = Olahraga:
## contrast estimate SE df t.ratio p.value
## M1 - M0 -0.3063 0.0623 60 -4.912 <0.0001
## M2 - M0 -0.6438 0.1270 60 -5.080 <0.0001
## M3 - M0 -0.9187 0.1810 60 -5.071 <0.0001
##
## kelompok = ADF + Olahraga:
## contrast estimate SE df t.ratio p.value
## M1 - M0 -1.8875 0.0623 60 -30.276 <0.0001
## M2 - M0 -4.2375 0.1270 60 -33.439 <0.0001
## M3 - M0 -6.2438 0.1810 60 -34.462 <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 - ADF -0.6813 2.49 60 -0.273 0.9928
## Kontrol - Olahraga 0.0813 2.49 60 0.033 1.0000
## Kontrol - (ADF + Olahraga) 2.4937 2.49 60 1.000 0.7499
## ADF - Olahraga 0.7625 2.49 60 0.306 0.9900
## ADF - (ADF + Olahraga) 3.1750 2.49 60 1.273 0.5833
## Olahraga - (ADF + Olahraga) 2.4125 2.49 60 0.967 0.7683
##
## waktu = M1:
## contrast estimate SE df t.ratio p.value
## Kontrol - ADF 0.4500 2.50 60 0.180 0.9979
## Kontrol - Olahraga 0.3688 2.50 60 0.147 0.9988
## Kontrol - (ADF + Olahraga) 4.3625 2.50 60 1.743 0.3111
## ADF - Olahraga -0.0813 2.50 60 -0.032 1.0000
## ADF - (ADF + Olahraga) 3.9125 2.50 60 1.563 0.4073
## Olahraga - (ADF + Olahraga) 3.9937 2.50 60 1.595 0.3889
##
## waktu = M2:
## contrast estimate SE df t.ratio p.value
## Kontrol - ADF 1.3250 2.52 60 0.527 0.9523
## Kontrol - Olahraga 0.6813 2.52 60 0.271 0.9930
## Kontrol - (ADF + Olahraga) 6.6875 2.52 60 2.658 0.0481
## ADF - Olahraga -0.6438 2.52 60 -0.256 0.9941
## ADF - (ADF + Olahraga) 5.3625 2.52 60 2.131 0.1550
## Olahraga - (ADF + Olahraga) 6.0062 2.52 60 2.387 0.0905
##
## waktu = M3:
## contrast estimate SE df t.ratio p.value
## Kontrol - ADF 2.3188 2.52 60 0.919 0.7950
## Kontrol - Olahraga 0.9313 2.52 60 0.369 0.9827
## Kontrol - (ADF + Olahraga) 8.6687 2.52 60 3.434 0.0058
## ADF - Olahraga -1.3875 2.52 60 -0.550 0.9463
## ADF - (ADF + Olahraga) 6.3500 2.52 60 2.516 0.0676
## Olahraga - (ADF + Olahraga) 7.7375 2.52 60 3.065 0.0167
##
## P value adjustment: tukey method for comparing a family of 4 estimates
## 4d. Kontras interaksi -------------------------------------------------------
# Fokus: apakah perubahan dari Bulan 0 ke Bulan 3 berbeda antarkelompok?
em_full <- emmeans(
aov2,
~ waktu * kelompok
)
change_m0_m3 <- contrast(
em_full,
interaction = c("revpairwise", "revpairwise"),
adjust = "holm"
)
change_m0_m3
## waktu_revpairwise kelompok_revpairwise estimate SE df t.ratio
## M1 - M0 ADF - Kontrol -1.131 0.0882 60 -12.831
## M2 - M0 ADF - Kontrol -2.006 0.1790 60 -11.195
## M2 - M1 ADF - Kontrol -0.875 0.1050 60 -8.337
## M3 - M0 ADF - Kontrol -3.000 0.2560 60 -11.708
## M3 - M1 ADF - Kontrol -1.869 0.1800 60 -10.386
## M3 - M2 ADF - Kontrol -0.994 0.0921 60 -10.787
## M1 - M0 Olahraga - Kontrol -0.287 0.0882 60 -3.261
## M2 - M0 Olahraga - Kontrol -0.600 0.1790 60 -3.348
## M2 - M1 Olahraga - Kontrol -0.312 0.1050 60 -2.977
## M3 - M0 Olahraga - Kontrol -0.850 0.2560 60 -3.317
## M3 - M1 Olahraga - Kontrol -0.562 0.1800 60 -3.126
## M3 - M2 Olahraga - Kontrol -0.250 0.0921 60 -2.714
## M1 - M0 Olahraga - ADF 0.844 0.0882 60 9.570
## M2 - M0 Olahraga - ADF 1.406 0.1790 60 7.847
## M2 - M1 Olahraga - ADF 0.562 0.1050 60 5.359
## M3 - M0 Olahraga - ADF 2.150 0.2560 60 8.391
## M3 - M1 Olahraga - ADF 1.306 0.1800 60 7.259
## M3 - M2 Olahraga - ADF 0.744 0.0921 60 8.073
## M1 - M0 (ADF + Olahraga) - Kontrol -1.869 0.0882 60 -21.196
## M2 - M0 (ADF + Olahraga) - Kontrol -4.194 0.1790 60 -23.401
## M2 - M1 (ADF + Olahraga) - Kontrol -2.325 0.1050 60 -22.152
## M3 - M0 (ADF + Olahraga) - Kontrol -6.175 0.2560 60 -24.100
## M3 - M1 (ADF + Olahraga) - Kontrol -4.306 0.1800 60 -23.932
## M3 - M2 (ADF + Olahraga) - Kontrol -1.981 0.0921 60 -21.506
## M1 - M0 (ADF + Olahraga) - ADF -0.738 0.0882 60 -8.365
## M2 - M0 (ADF + Olahraga) - ADF -2.188 0.1790 60 -12.206
## M2 - M1 (ADF + Olahraga) - ADF -1.450 0.1050 60 -13.815
## M3 - M0 (ADF + Olahraga) - ADF -3.175 0.2560 60 -12.391
## M3 - M1 (ADF + Olahraga) - ADF -2.438 0.1800 60 -13.546
## M3 - M2 (ADF + Olahraga) - ADF -0.988 0.0921 60 -10.719
## M1 - M0 (ADF + Olahraga) - Olahraga -1.581 0.0882 60 -17.935
## M2 - M0 (ADF + Olahraga) - Olahraga -3.594 0.1790 60 -20.053
## M2 - M1 (ADF + Olahraga) - Olahraga -2.013 0.1050 60 -19.175
## M3 - M0 (ADF + Olahraga) - Olahraga -5.325 0.2560 60 -20.783
## M3 - M1 (ADF + Olahraga) - Olahraga -3.744 0.1800 60 -20.806
## M3 - M2 (ADF + Olahraga) - Olahraga -1.731 0.0921 60 -18.792
## p.value
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## 0.0085
## 0.0085
## 0.0085
## 0.0085
## 0.0085
## 0.0087
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
## <0.0001
##
## P value adjustment: holm method for 36 tests
# Alternatif: perubahan aktual tiap subjek dari baseline ke Bulan 3
change_data <- dat_wide |>
mutate(
perubahan_M3_M0 = BB_Bulan_3_kg - BB_Bulan_0_kg
)
change_data |>
group_by(kelompok) |>
get_summary_stats(
perubahan_M3_M0,
type = "mean_sd"
)
## # A tibble: 4 × 5
## kelompok variable n mean sd
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol perubahan_M3_M0 16 -0.069 0.166
## 2 ADF perubahan_M3_M0 16 -3.07 0.436
## 3 Olahraga perubahan_M3_M0 16 -0.919 0.269
## 4 ADF + Olahraga perubahan_M3_M0 16 -6.24 1.35
# Uji beda perubahan M3-M0 antar kelompok
anova_change <- aov(
perubahan_M3_M0 ~ kelompok,
data = change_data
)
summary(anova_change)
## Df Sum Sq Mean Sq F value Pr(>F)
## kelompok 3 363.6 121.22 230.8 <2e-16 ***
## Residuals 60 31.5 0.53
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Post hoc perubahan antar kelompok
TukeyHSD(anova_change)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = perubahan_M3_M0 ~ kelompok, data = change_data)
##
## $kelompok
## diff lwr upr p adj
## ADF-Kontrol -3.000 -3.677079 -2.3229211 0.0000000
## Olahraga-Kontrol -0.850 -1.527079 -0.1729211 0.0082036
## ADF + Olahraga-Kontrol -6.175 -6.852079 -5.4979211 0.0000000
## Olahraga-ADF 2.150 1.472921 2.8270789 0.0000000
## ADF + Olahraga-ADF -3.175 -3.852079 -2.4979211 0.0000000
## ADF + Olahraga-Olahraga -5.325 -6.002079 -4.6479211 0.0000000
# Tren linear per kelompok
tren_kelompok <- emmeans(
aov2,
~ waktu | kelompok
)
contrast(
tren_kelompok,
"poly"
)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## linear -0.23125 0.6110 60 -0.378 0.7064
## quadratic -0.00625 0.0452 60 -0.138 0.8905
## cubic 0.00625 0.1000 60 0.062 0.9505
##
## kelompok = ADF:
## contrast estimate SE df t.ratio p.value
## linear -10.10625 0.6110 60 -16.541 <0.0001
## quadratic 0.13125 0.0452 60 2.903 0.0052
## cubic -0.36875 0.1000 60 -3.679 0.0005
##
## kelompok = Olahraga:
## contrast estimate SE df t.ratio p.value
## linear -3.09375 0.6110 60 -5.064 <0.0001
## quadratic 0.03125 0.0452 60 0.691 0.4921
## cubic 0.09375 0.1000 60 0.935 0.3533
##
## kelompok = ADF + Olahraga:
## contrast estimate SE df t.ratio p.value
## linear -21.08125 0.6110 60 -34.504 <0.0001
## quadratic -0.11875 0.0452 60 -2.626 0.0109
## cubic 0.80625 0.1000 60 8.045 <0.0001
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Model dengan random intercept
lmm1 <- lmer(
bb ~ kelompok * waktu + (1 | id),
data = dat_long,
REML = TRUE
)
summary(lmm1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: bb ~ kelompok * waktu + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 660.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0530 -0.2838 0.0157 0.2841 3.8672
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 50.2663 7.090
## Residual 0.1049 0.324
## Number of obs: 256, groups: id, 64
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 92.39922 0.88647 60.00000 104.233 < 2e-16 ***
## kelompok1 1.73047 1.53540 60.00000 1.127 0.264
## kelompok2 0.87734 1.53540 60.00000 0.571 0.570
## kelompok3 1.21484 1.53540 60.00000 0.791 0.432
## waktu1 1.28984 0.03507 180.00000 36.780 < 2e-16 ***
## waktu2 0.44922 0.03507 180.00000 12.809 < 2e-16 ***
## waktu3 -0.45391 0.03507 180.00000 -12.943 < 2e-16 ***
## kelompok1:waktu1 -1.25703 0.06074 180.00000 -20.695 < 2e-16 ***
## kelompok2:waktu1 0.27734 0.06074 180.00000 4.566 9.19e-06 ***
## kelompok3:waktu1 -0.82266 0.06074 180.00000 -13.543 < 2e-16 ***
## kelompok1:waktu2 -0.43516 0.06074 180.00000 -7.164 1.94e-11 ***
## kelompok2:waktu2 -0.03203 0.06074 180.00000 -0.527 0.599
## kelompok3:waktu2 -0.28828 0.06074 180.00000 -4.746 4.22e-06 ***
## kelompok1:waktu3 0.44297 0.06074 180.00000 7.293 9.31e-12 ***
## kelompok2:waktu3 -0.02891 0.06074 180.00000 -0.476 0.635
## kelompok3:waktu3 0.27734 0.06074 180.00000 4.566 9.19e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
anova(lmm1, 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 0.663 0.221 3 60 2.1053 0.109
## waktu 238.282 79.427 3 180 756.8253 <2e-16 ***
## kelompok:waktu 207.428 23.048 9 180 219.6091 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Model dengan random intercept + random slope waktu
lmm2 <- lmer(
bb ~ kelompok * waktu + (1 + bulan | id),
data = dat_long,
REML = TRUE
)
summary(lmm2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: bb ~ kelompok * waktu + (1 + bulan | id)
## Data: dat_long
##
## REML criterion at convergence: 414.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15424 -0.45495 0.02092 0.42375 2.32736
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 49.791513 7.05631
## bulan 0.058104 0.24105 0.07
## Residual 0.008106 0.09004
## Number of obs: 256, groups: id, 64
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 92.39922 0.88645 60.00291 104.235 < 2e-16 ***
## kelompok1 1.73047 1.53538 60.00301 1.127 0.264202
## kelompok2 0.87734 1.53538 60.00301 0.571 0.569849
## kelompok3 1.21484 1.53538 60.00301 0.791 0.431923
## waktu1 1.28984 0.04624 62.18323 27.897 < 2e-16 ***
## waktu2 0.44922 0.01794 106.57779 25.035 < 2e-16 ***
## waktu3 -0.45391 0.01794 106.57779 -25.297 < 2e-16 ***
## kelompok1:waktu1 -1.25703 0.08008 62.18323 -15.697 < 2e-16 ***
## kelompok2:waktu1 0.27734 0.08008 62.18323 3.463 0.000972 ***
## kelompok3:waktu1 -0.82266 0.08008 62.18323 -10.273 5.05e-15 ***
## kelompok1:waktu2 -0.43516 0.03108 106.57779 -14.002 < 2e-16 ***
## kelompok2:waktu2 -0.03203 0.03108 106.57779 -1.031 0.305041
## kelompok3:waktu2 -0.28828 0.03108 106.57779 -9.276 2.29e-15 ***
## kelompok1:waktu3 0.44297 0.03108 106.57779 14.253 < 2e-16 ***
## kelompok2:waktu3 -0.02891 0.03108 106.57779 -0.930 0.354424
## kelompok3:waktu3 0.27734 0.03108 106.57779 8.924 1.42e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
# Perbandingan struktur random effect
anova(
lmm1,
lmm2,
refit = FALSE
)
## Data: dat_long
## Models:
## lmm1: bb ~ kelompok * waktu + (1 | id)
## lmm2: bb ~ kelompok * waktu + (1 + bulan | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 18 696.11 759.92 -330.05 660.11
## lmm2 20 454.67 525.58 -207.34 414.67 245.44 2 < 2.2e-16 ***
## ---
## 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 0.0512 0.01707 3 60.00 2.1054 0.109
## waktu 6.5260 2.17533 3 127.51 266.5624 <2e-16 ***
## kelompok:waktu 6.3177 0.70197 9 148.38 85.8264 <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.998
## Unadjusted ICC: 0.880
# 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))
# Simulasi 30 nilai hilang (MCAR) untuk latihan LMM
set.seed(1)
dat_miss <- dat_long
dat_miss$bb[
sample(which(dat_miss$waktu != "M0"), 30)
] <- NA
lmm_miss <- lmer(
bb ~ kelompok * waktu + (1 + bulan | 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 0.0511 0.01704 3 60.00 2.1112 0.1082
## waktu 6.1274 2.04246 3 107.46 251.1404 <2e-16 ***
## kelompok:waktu 5.9547 0.66163 9 123.20 81.1580 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Jumlah subjek yang memiliki setidaknya satu nilai hilang
n_distinct(
dat_miss$id[is.na(dat_miss$bb)]
)
## [1] 25
# 6. MENYIMPAN DATA & SESSION INFO
write.csv(
dat_wide,
"data_obesitas_wide_RStudio.csv",
row.names = FALSE
)
write.csv(
dat_long,
"data_obesitas_long_RStudio.csv",
row.names = FALSE
)
write.csv(
desk,
"deskriptif_berat_badan.csv",
row.names = FALSE
)
ringkasan_perubahan <- change_data |>
group_by(kelompok) |>
summarise(
n = n(),
mean_baseline = mean(BB_Bulan_0_kg, na.rm = TRUE),
mean_bulan3 = mean(BB_Bulan_3_kg, na.rm = TRUE),
mean_perubahan = mean(perubahan_M3_M0, na.rm = TRUE),
sd_perubahan = sd(perubahan_M3_M0, na.rm = TRUE),
.groups = "drop"
)
ringkasan_perubahan
## # A tibble: 4 × 6
## kelompok n mean_baseline mean_bulan3 mean_perubahan sd_perubahan
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Kontrol 16 94.2 94.1 -0.0688 0.166
## 2 ADF 16 94.8 91.8 -3.07 0.436
## 3 Olahraga 16 94.1 93.2 -0.919 0.269
## 4 ADF + Olahraga 16 91.7 85.4 -6.24 1.35
write.csv(
ringkasan_perubahan,
"ringkasan_perubahan_berat_badan.csv",
row.names = FALSE
)
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.60 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 WRS2_1.1-7 RColorBrewer_1.1-3
## [19] S7_0.2.2 lifecycle_1.0.5 compiler_4.6.1
## [22] farver_2.1.2 stringr_1.6.0 htmltools_0.5.9
## [25] sass_0.4.10 yaml_2.3.12 Formula_1.2-6
## [28] ggpubr_1.0.0 pillar_1.11.1 nloptr_2.2.1
## [31] jquerylib_0.1.4 MASS_7.3-65 cachem_1.1.0
## [34] reformulas_0.4.4 boot_1.3-32 abind_1.4-8
## [37] nlme_3.1-169 tidyselect_1.2.1 digest_0.6.39
## [40] performance_0.18.2 mvtnorm_1.4-2 stringi_1.8.9
## [43] reshape2_1.4.5 purrr_1.2.2 labeling_0.4.3
## [46] splines_4.6.1 fastmap_1.2.0 grid_4.6.1
## [49] cli_3.6.6 magrittr_2.0.5 utf8_1.2.6
## [52] broom_1.0.13 withr_3.0.3 scales_1.4.0
## [55] backports_1.5.1 estimability_2.0.0 rmarkdown_2.31
## [58] ggsignif_0.6.4 evaluate_1.0.5 knitr_1.51
## [61] parameters_0.29.3 rbibutils_2.4.1 rlang_1.3.0
## [64] Rcpp_1.1.2 glue_1.8.1 reshape_0.8.10
## [67] rstudioapi_0.19.0 minqa_1.2.8 jsonlite_2.0.0
## [70] R6_2.6.1 plyr_1.8.9