Analisis Repeated Measures Data Balita 3 Kelompok Selama 3 Bulan
Nurul Adha
2026-10-06
Deskripsi Analisis
Analisis ini mengevaluasi perubahan berat badan balita pada tiga kelompok (Kontrol, PMT Biskuit, dan PMT Lokal) selama tiga bulan dengan pendekatan repeated-measures ANOVA dan linear mixed model. Data terdiri dari 90 balita, dengan tiga pengukuran berat badan per balita.
Catatan: File Excel
Data_Balita_3Kelompok_3Penimbangan_RStudio.xlsxharus berada di working directory yang sama dengan file R Markdown saat dilakukan knitting.
0. PAKET & PENGATURAN
# Jalankan sekali bila paket belum tersedia:
# install.packages(c("readxl", "dplyr", "tidyr", "ggplot2", "afex",
# "emmeans", "rstatix", "car", "effectsize", "lme4",
# "lmerTest", "performance", "ggpubr"))
suppressPackageStartupMessages({
library(readxl)
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 EXCEL LONG DAN WIDE
# Pastikan file Excel berada di Working Directory.
# Cek dengan:
getwd()## [1] "C:/Users/Hype AMD/Downloads/New folder"list.files()## [1] "Analisis_Balita_3Kelompok_3Bulan_RStudio--2-.Rmd"
## [2] "Analisis_Balita_3Kelompok_3Bulan_RStudio (2).Rmd"
## [3] "Analisis_Balita_3Kelompok_3Bulan_RStudio.R"
## [4] "Data_Balita_3Kelompok_3Penimbangan_RStudio.xlsx"# Jika file belum ditemukan, atur Working Directory melalui:
# Session -> Set Working Directory -> Choose Directory...
file_excel <- "Data_Balita_3Kelompok_3Penimbangan_RStudio.xlsx"
# Jika nama/path file berbeda:
# file_excel <- file.choose()
# 1A. Baca DATA WIDE dari sheet Data_Wide
data_wide <- read_excel(
file_excel,
sheet = "Data_Wide"
)
# 1B. Baca DATA LONG dari sheet Data_Long_R
data_long <- read_excel(
file_excel,
sheet = "Data_Long_R"
)
# Tampilkan data
# View(data_long)
# View(data_wide)
# Struktur dan dimensi
head(data_long)head(data_wide)str(data_long)## tibble [270 × 7] (S3: tbl_df/tbl/data.frame)
## $ ID : chr [1:270] "B001" "B001" "B001" "B002" ...
## $ Kelompok : chr [1:270] "Kontrol" "Kontrol" "Kontrol" "Kontrol" ...
## $ Usia_Bulan : num [1:270] 20 20 20 11 11 11 14 14 14 8 ...
## $ Jenis_Kelamin : chr [1:270] "Perempuan" "Perempuan" "Perempuan" "Perempuan" ...
## $ Waktu : chr [1:270] "BB_Bulan_1_kg" "BB_Bulan_2_kg" "BB_Bulan_3_kg" "BB_Bulan_1_kg" ...
## $ Berat_Badan_kg: num [1:270] 10.9 11 11.2 8.9 9 9.1 9.2 9.3 9.6 9.8 ...
## $ Bulan : num [1:270] 1 2 3 1 2 3 1 2 3 1 ...str(data_wide)## tibble [90 × 8] (S3: tbl_df/tbl/data.frame)
## $ ID : chr [1:90] "B001" "B002" "B003" "B004" ...
## $ Kelompok : chr [1:90] "Kontrol" "Kontrol" "Kontrol" "Kontrol" ...
## $ Usia_Bulan : num [1:90] 20 11 14 8 14 7 12 19 20 6 ...
## $ Jenis_Kelamin : chr [1:90] "Perempuan" "Perempuan" "Perempuan" "Laki-laki" ...
## $ BB_Bulan_1_kg : num [1:90] 10.9 8.9 9.2 9.8 10.5 8.1 9.3 9.5 10.7 8.7 ...
## $ BB_Bulan_2_kg : num [1:90] 11 9 9.3 10 10.6 8.5 9.4 9.8 11 9 ...
## $ BB_Bulan_3_kg : num [1:90] 11.2 9.1 9.6 10.2 10.9 8.7 9.5 10.1 11.2 9.1 ...
## $ Kenaikan_BB_1_3_Bulan_kg: num [1:90] 0.3 0.2 0.4 0.4 0.4 0.6 0.2 0.6 0.5 0.4 ...dim(data_long)## [1] 270 7dim(data_wide)## [1] 90 8names(data_long)## [1] "ID" "Kelompok" "Usia_Bulan" "Jenis_Kelamin"
## [5] "Waktu" "Berat_Badan_kg" "Bulan"names(data_wide)## [1] "ID" "Kelompok"
## [3] "Usia_Bulan" "Jenis_Kelamin"
## [5] "BB_Bulan_1_kg" "BB_Bulan_2_kg"
## [7] "BB_Bulan_3_kg" "Kenaikan_BB_1_3_Bulan_kg"# 1C. VALIDASI STRUKTUR DATA
# Kolom yang diharapkan pada LONG
stopifnot(
all(c("ID", "Kelompok", "Usia_Bulan", "Jenis_Kelamin",
"Waktu", "Berat_Badan_kg", "Bulan") %in% names(data_long))
)
# Kolom yang diharapkan pada WIDE
stopifnot(
all(c("ID", "Kelompok", "Usia_Bulan", "Jenis_Kelamin",
"BB_Bulan_1_kg", "BB_Bulan_2_kg", "BB_Bulan_3_kg") %in% names(data_wide))
)
# Pastikan jumlah responden 90
n_distinct(data_long$ID)## [1] 90n_distinct(data_wide$ID)## [1] 90# Jumlah responden tiap kelompok
data_wide %>% count(Kelompok)data_long %>%
distinct(ID, Kelompok) %>%
count(Kelompok)# Jumlah pengukuran tiap responden
data_long %>%
count(ID) %>%
count(n, name = "jumlah_responden")# Cek duplikasi ID x Bulan
cek_duplikasi <- data_long %>%
count(ID, Bulan) %>%
filter(n != 1)
cek_duplikasi# Cek missing value
colSums(is.na(data_long))## ID Kelompok Usia_Bulan Jenis_Kelamin Waktu
## 0 0 0 0 0
## Berat_Badan_kg Bulan
## 0 0colSums(is.na(data_wide))## ID Kelompok Usia_Bulan
## 0 0 0
## Jenis_Kelamin BB_Bulan_1_kg BB_Bulan_2_kg
## 0 0 0
## BB_Bulan_3_kg Kenaikan_BB_1_3_Bulan_kg
## 0 0# Ringkasan outcome
summary(data_long$Berat_Badan_kg)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.50 9.20 10.10 10.15 11.00 14.40# 1D. RECODING VARIABEL
# LONG digunakan sebagai basis analisis.
# Bulan_num tetap numerik untuk grafik/LMM,
# sedangkan Bulan menjadi faktor untuk repeated-measures ANOVA.
data_long <- data_long %>%
mutate(
ID = factor(ID),
Kelompok = factor(
Kelompok,
levels = c("Kontrol", "PMT Biskuit", "PMT Lokal")
),
Jenis_Kelamin = factor(
Jenis_Kelamin,
levels = c("Laki-laki", "Perempuan")
),
Bulan_num = as.numeric(Bulan),
Bulan = factor(
Bulan,
levels = c(1, 2, 3),
labels = c("Bulan_1", "Bulan_2", "Bulan_3")
)
)
# WIDE juga diberi factor
data_wide <- data_wide %>%
mutate(
ID = factor(ID),
Kelompok = factor(
Kelompok,
levels = c("Kontrol", "PMT Biskuit", "PMT Lokal")
),
Jenis_Kelamin = factor(
Jenis_Kelamin,
levels = c("Laki-laki", "Perempuan")
)
)
# Cek kembali level
levels(data_long$Kelompok)## [1] "Kontrol" "PMT Biskuit" "PMT Lokal"levels(data_long$Bulan)## [1] "Bulan_1" "Bulan_2" "Bulan_3"levels(data_long$Jenis_Kelamin)## [1] "Laki-laki" "Perempuan"# 1E. CEK KONSISTENSI LONG VS WIDE
wide_from_long <- data_long %>%
select(ID, Kelompok, Bulan, Berat_Badan_kg) %>%
pivot_wider(
names_from = Bulan,
values_from = Berat_Badan_kg,
names_prefix = "BB_"
)
# View(wide_from_long)
nrow(wide_from_long)## [1] 90nrow(data_wide)## [1] 902. EKSPLORASI DATA
# 2A. Statistik deskriptif per kelompok dan waktu
deskriptif <- data_long %>%
group_by(Kelompok, Bulan) %>%
summarise(
n = n(),
mean = mean(Berat_Badan_kg, na.rm = TRUE),
sd = sd(Berat_Badan_kg, na.rm = TRUE),
median = median(Berat_Badan_kg, na.rm = TRUE),
min = min(Berat_Badan_kg, na.rm = TRUE),
max = max(Berat_Badan_kg, na.rm = TRUE),
.groups = "drop"
)
deskriptif# Versi rstatix
data_long %>%
group_by(Kelompok, Bulan) %>%
get_summary_stats(Berat_Badan_kg, type = "mean_sd")# 2B. Deskriptif usia dan jenis kelamin
deskriptif_usia <- data_wide %>%
group_by(Kelompok) %>%
summarise(
n = n(),
mean_usia = mean(Usia_Bulan, na.rm = TRUE),
sd_usia = sd(Usia_Bulan, na.rm = TRUE),
median_usia = median(Usia_Bulan, na.rm = TRUE),
min_usia = min(Usia_Bulan, na.rm = TRUE),
max_usia = max(Usia_Bulan, na.rm = TRUE),
.groups = "drop"
)
deskriptif_usiatabel_jenis_kelamin <- data_wide %>%
count(Kelompok, Jenis_Kelamin) %>%
group_by(Kelompok) %>%
mutate(persen = round(100 * n / sum(n), 1))
tabel_jenis_kelamin# 2C. Matriks kovarians dan korelasi antar waktu
vars_waktu <- c(
"BB_Bulan_1_kg", "BB_Bulan_2_kg", "BB_Bulan_3_kg"
)
S <- cov(data_wide[, vars_waktu], use = "complete.obs")
R <- cor(data_wide[, vars_waktu], use = "complete.obs")
round(S, 2)## BB_Bulan_1_kg BB_Bulan_2_kg BB_Bulan_3_kg
## BB_Bulan_1_kg 1.31 1.31 1.32
## BB_Bulan_2_kg 1.31 1.33 1.35
## BB_Bulan_3_kg 1.32 1.35 1.38round(R, 3)## BB_Bulan_1_kg BB_Bulan_2_kg BB_Bulan_3_kg
## BB_Bulan_1_kg 1.000 0.993 0.982
## BB_Bulan_2_kg 0.993 1.000 0.994
## BB_Bulan_3_kg 0.982 0.994 1.000# 2D. 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, 3)## BB_Bulan_1_kg - BB_Bulan_2_kg BB_Bulan_1_kg - BB_Bulan_3_kg
## 0.019 0.050
## BB_Bulan_2_kg - BB_Bulan_3_kg
## 0.018# 2E. Profile plot: rerata +/- 95% CI
p_profil <- ggplot(
data_long,
aes(
x = Bulan_num,
y = Berat_Badan_kg,
colour = Kelompok,
group = Kelompok
)
) +
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 = .08
) +
scale_x_continuous(
breaks = c(1, 2, 3),
labels = c("Bulan 1", "Bulan 2", "Bulan 3")
) +
labs(
x = "Waktu pengukuran",
y = "Berat badan (kg)",
colour = "Kelompok",
title = "Profil rerata berat badan selama 3 bulan"
) +
theme(legend.position = "bottom")
p_profil
# 2F. Spaghetti plot
p_spag <- ggplot(
data_long,
aes(
x = Bulan_num,
y = Berat_Badan_kg,
group = ID
)
) +
geom_line(alpha = .25) +
stat_summary(
aes(group = Kelompok),
fun = mean,
geom = "line",
linewidth = 1.3
) +
facet_wrap(~ Kelompok) +
scale_x_continuous(
breaks = c(1, 2, 3),
labels = c("Bulan 1", "Bulan 2", "Bulan 3")
) +
labs(
x = "Bulan",
y = "Berat badan (kg)",
title = "Lintasan individu dan rerata kelompok"
)
p_spag
3. REPEATED MEASURES ANOVA SATU ARAH
# Pertanyaan:
# Apakah berat badan berubah selama 3 bulan pada kelompok PMT Biskuit?
d1 <- data_long %>%
filter(Kelompok == "PMT Biskuit") %>%
droplevels()
d1w <- data_wide %>%
filter(Kelompok == "PMT Biskuit") %>%
droplevels()
# 3A. UJI ASUMSI
# (i) Outlier per waktu
outlier_1way <- d1 %>%
group_by(Bulan) %>%
identify_outliers(Berat_Badan_kg)
outlier_1way# (ii) Normalitas Shapiro-Wilk per waktu
shapiro_1way <- d1 %>%
group_by(Bulan) %>%
shapiro_test(Berat_Badan_kg)
shapiro_1way# Q-Q plot
qqplot_1way <- ggpubr::ggqqplot(
d1,
"Berat_Badan_kg",
facet.by = "Bulan"
)
qqplot_1way
# (iii) Mauchly's Test
# anova_test melaporkan Mauchly + GG/HF
aov1_rs <- anova_test(
data = d1,
dv = Berat_Badan_kg,
wid = ID,
within = Bulan,
effect.size = "pes"
)
aov1_rs## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 Bulan 2 58 430.106 1.64e-35 * 0.937
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 Bulan 0.883 0.176
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 Bulan 0.896 1.79, 51.94 4.71e-32 * 0.951 1.9, 55.14 7.06e-34
## p[HF]<.05
## 1 *get_anova_table(aov1_rs, correction = "auto")# 3B. REPEATED MEASURES ANOVA DENGAN AFEX
aov1 <- aov_ez(
id = "ID",
dv = "Berat_Badan_kg",
data = d1,
within = "Bulan",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov1## Anova Table (Type 3 tests)
##
## Response: Berat_Badan_kg
## Effect df MSE F ges pes p.value
## 1 Bulan 1.79, 51.94 0.01 430.11 *** .051 .937 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GGsummary(aov1)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 9379.9 1 122.186 29 2226.26 < 2.2e-16 ***
## Bulan 6.6 2 0.445 58 430.11 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## Bulan 0.88331 0.17602
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Bulan 0.8955 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Bulan 0.9506071 7.059862e-34nice(aov1, correction = "GG", es = c("ges", "pes"))# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)omega_squared(aov1, partial = TRUE)# 3C. POST-HOC DAN TREN
em1 <- emmeans(aov1, ~ Bulan)
em1## Bulan emmean SE df lower.CL upper.CL
## Bulan_1 9.88 0.217 29 9.44 10.3
## Bulan_2 10.20 0.216 29 9.76 10.6
## Bulan_3 10.54 0.217 29 10.10 11.0
##
## Confidence level used: 0.95# Semua pasangan waktu
pairs(em1, adjust = "holm")## contrast estimate SE df t.ratio p.value
## Bulan_1 - Bulan_2 -0.323 0.0218 29 -14.828 <0.0001
## Bulan_1 - Bulan_3 -0.663 0.0260 29 -25.479 <0.0001
## Bulan_2 - Bulan_3 -0.340 0.0195 29 -17.405 <0.0001
##
## P value adjustment: holm method for 3 tests# Masing-masing waktu dibandingkan dengan Bulan 1
contrast(
em1,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)## contrast estimate SE df t.ratio p.value
## Bulan_2 - Bulan_1 0.323 0.0218 29 14.828 <0.0001
## Bulan_3 - Bulan_1 0.663 0.0260 29 25.479 <0.0001
##
## P value adjustment: holm method for 2 tests# Tren linear dan kuadratik
contrast(em1, "poly")## contrast estimate SE df t.ratio p.value
## linear 0.6633 0.0260 29 25.479 <0.0001
## quadratic 0.0167 0.0322 29 0.518 0.6086# 3D. ALTERNATIF NONPARAMETRIK: FRIEDMAN
friedman_test(d1, Berat_Badan_kg ~ Bulan | ID)friedman_effsize(d1, Berat_Badan_kg ~ Bulan | ID)# Wilcoxon berpasangan sebagai post-hoc alternatif
d1 %>%
wilcox_test(
Berat_Badan_kg ~ Bulan,
paired = TRUE,
p.adjust.method = "holm"
)4. MIXED DESIGN ANOVA
# Pertanyaan utama:
# Apakah perubahan berat badan dari Bulan 1 sampai Bulan 3 berbeda
# antar kelompok Kontrol, PMT Biskuit, dan PMT Lokal?
# 4A. UJI ASUMSI
# (i) Outlier per sel
outlier_mixed <- data_long %>%
group_by(Kelompok, Bulan) %>%
identify_outliers(Berat_Badan_kg)
outlier_mixed# (ii) Normalitas Shapiro-Wilk per Kelompok x Bulan
shapiro_mixed <- data_long %>%
group_by(Kelompok, Bulan) %>%
shapiro_test(Berat_Badan_kg)
shapiro_mixed# Q-Q plot per sel
qqplot_mixed <- ggpubr::ggqqplot(
data_long,
"Berat_Badan_kg"
) +
facet_grid(Bulan ~ Kelompok)
qqplot_mixed
# (iii) Homogenitas varians pada setiap waktu: Levene
levene_mixed <- data_long %>%
group_by(Bulan) %>%
levene_test(Berat_Badan_kg ~ Kelompok)
levene_mixed# (iv) Homogenitas matriks kovarians: Box's M
box_m_result <- box_m(
data_wide[, vars_waktu],
data_wide$Kelompok
)
box_m_result# (v) Sphericity: Mauchly dilihat pada summary(aov2)
# 4B. MIXED DESIGN REPEATED MEASURES ANOVA
aov2 <- aov_ez(
id = "ID",
dv = "Berat_Badan_kg",
data = data_long,
between = "Kelompok",
within = "Bulan",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov2## Anova Table (Type 3 tests)
##
## Response: Berat_Badan_kg
## Effect df MSE F ges pes p.value
## 1 Kelompok 2, 87 4.02 0.65 .015 .015 .526
## 2 Bulan 1.39, 120.61 0.02 807.69 *** .046 .903 <.001
## 3 Kelompok:Bulan 2.77, 120.61 0.02 17.31 *** .002 .285 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GGsummary(aov2)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 27814.0 1 350.06 87 6912.5129 < 2.2e-16 ***
## Kelompok 5.2 2 350.06 87 0.6475 0.5258
## Bulan 17.1 2 1.84 174 807.6924 < 2.2e-16 ***
## Kelompok:Bulan 0.7 4 1.84 174 17.3057 5.706e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## Bulan 0.55733 1.2098e-11
## Kelompok:Bulan 0.55733 1.2098e-11
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Bulan 0.69316 < 2.2e-16 ***
## Kelompok:Bulan 0.69316 5.827e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Bulan 0.7007975 4.517553e-63
## Kelompok:Bulan 0.7007975 4.900643e-09nice(aov2, correction = "GG", es = c("ges", "pes"))# Alternatif rstatix
aov2_rs <- anova_test(
data = data_long,
dv = Berat_Badan_kg,
wid = ID,
between = Kelompok,
within = Bulan,
effect.size = "pes",
type = 3
)
aov2_rs## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 Kelompok 2 87 0.648 5.26e-01 0.015
## 2 Bulan 2 174 807.692 8.76e-89 * 0.903
## 3 Kelompok:Bulan 4 174 17.306 5.71e-12 * 0.285
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 Bulan 0.557 1.21e-11 *
## 2 Kelompok:Bulan 0.557 1.21e-11 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF]
## 1 Bulan 0.693 1.39, 120.61 2.05e-62 * 0.701 1.4, 121.94
## 2 Kelompok:Bulan 0.693 2.77, 120.61 5.83e-09 * 0.701 2.8, 121.94
## p[HF] p[HF]<.05
## 1 4.52e-63 *
## 2 4.90e-09 *get_anova_table(aov2_rs, correction = "GG")# 4C. UKURAN EFEK
eta_squared(aov2, partial = TRUE)omega_squared(aov2, partial = TRUE)# 4D. PLOT INTERAKSI
p_interaksi <- afex_plot(
aov2,
x = "Bulan",
trace = "Kelompok",
error = "within",
mapping = c("colour", "shape", "linetype")
) +
labs(
y = "Berat badan (kg)",
x = "Waktu",
title = "Interaksi Kelompok x Waktu pada berat badan"
) +
theme(legend.position = "bottom")
p_interaksi
# 4E. SIMPLE EFFECTS
# Estimated marginal means waktu dalam setiap kelompok
em2 <- emmeans(aov2, ~ Bulan | Kelompok)
em2## Kelompok = Kontrol:
## Bulan emmean SE df lower.CL upper.CL
## Bulan_1 9.74 0.211 87 9.32 10.2
## Bulan_2 9.96 0.212 87 9.54 10.4
## Bulan_3 10.18 0.214 87 9.76 10.6
##
## Kelompok = PMT Biskuit:
## Bulan emmean SE df lower.CL upper.CL
## Bulan_1 9.88 0.211 87 9.46 10.3
## Bulan_2 10.20 0.212 87 9.78 10.6
## Bulan_3 10.54 0.214 87 10.12 11.0
##
## Kelompok = PMT Lokal:
## Bulan emmean SE df lower.CL upper.CL
## Bulan_1 9.93 0.211 87 9.51 10.3
## Bulan_2 10.25 0.212 87 9.83 10.7
## Bulan_3 10.67 0.214 87 10.24 11.1
##
## Confidence level used: 0.95# Efek waktu di dalam masing-masing kelompok
joint_tests(aov2, by = "Kelompok")# Efek kelompok pada masing-masing waktu
joint_tests(aov2, by = "Bulan")# 4F. POST-HOC DALAM KELOMPOK
# Semua pasangan waktu dalam masing-masing kelompok
pairs(
em2,
adjust = "holm"
)## Kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## Bulan_1 - Bulan_2 -0.220 0.0237 87 -9.301 <0.0001
## Bulan_1 - Bulan_3 -0.443 0.0341 87 -13.000 <0.0001
## Bulan_2 - Bulan_3 -0.223 0.0198 87 -11.295 <0.0001
##
## Kelompok = PMT Biskuit:
## contrast estimate SE df t.ratio p.value
## Bulan_1 - Bulan_2 -0.323 0.0237 87 -13.669 <0.0001
## Bulan_1 - Bulan_3 -0.663 0.0341 87 -19.451 <0.0001
## Bulan_2 - Bulan_3 -0.340 0.0198 87 -17.196 <0.0001
##
## Kelompok = PMT Lokal:
## contrast estimate SE df t.ratio p.value
## Bulan_1 - Bulan_2 -0.327 0.0237 87 -13.810 <0.0001
## Bulan_1 - Bulan_3 -0.740 0.0341 87 -21.699 <0.0001
## Bulan_2 - Bulan_3 -0.413 0.0198 87 -20.905 <0.0001
##
## P value adjustment: holm method for 3 tests# Tiap waktu vs Bulan 1
contrast(
em2,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)## Kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## Bulan_2 - Bulan_1 0.220 0.0237 87 9.301 <0.0001
## Bulan_3 - Bulan_1 0.443 0.0341 87 13.000 <0.0001
##
## Kelompok = PMT Biskuit:
## contrast estimate SE df t.ratio p.value
## Bulan_2 - Bulan_1 0.323 0.0237 87 13.669 <0.0001
## Bulan_3 - Bulan_1 0.663 0.0341 87 19.451 <0.0001
##
## Kelompok = PMT Lokal:
## contrast estimate SE df t.ratio p.value
## Bulan_2 - Bulan_1 0.327 0.0237 87 13.810 <0.0001
## Bulan_3 - Bulan_1 0.740 0.0341 87 21.699 <0.0001
##
## P value adjustment: holm method for 2 tests# 4G. POST-HOC ANTARKELOMPOK PADA SETIAP WAKTU
em2b <- emmeans(aov2, ~ Kelompok | Bulan)
em2b## Bulan = Bulan_1:
## Kelompok emmean SE df lower.CL upper.CL
## Kontrol 9.74 0.211 87 9.32 10.2
## PMT Biskuit 9.88 0.211 87 9.46 10.3
## PMT Lokal 9.93 0.211 87 9.51 10.3
##
## Bulan = Bulan_2:
## Kelompok emmean SE df lower.CL upper.CL
## Kontrol 9.96 0.212 87 9.54 10.4
## PMT Biskuit 10.20 0.212 87 9.78 10.6
## PMT Lokal 10.25 0.212 87 9.83 10.7
##
## Bulan = Bulan_3:
## Kelompok emmean SE df lower.CL upper.CL
## Kontrol 10.18 0.214 87 9.76 10.6
## PMT Biskuit 10.54 0.214 87 10.12 11.0
## PMT Lokal 10.67 0.214 87 10.24 11.1
##
## Confidence level used: 0.95# Tukey
pairs(em2b, adjust = "tukey")## Bulan = Bulan_1:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit -0.1433 0.298 87 -0.481 0.8804
## Kontrol - PMT Lokal -0.1900 0.298 87 -0.638 0.7997
## PMT Biskuit - PMT Lokal -0.0467 0.298 87 -0.157 0.9866
##
## Bulan = Bulan_2:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit -0.2467 0.299 87 -0.824 0.6893
## Kontrol - PMT Lokal -0.2967 0.299 87 -0.991 0.5846
## PMT Biskuit - PMT Lokal -0.0500 0.299 87 -0.167 0.9847
##
## Bulan = Bulan_3:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit -0.3633 0.302 87 -1.203 0.4549
## Kontrol - PMT Lokal -0.4867 0.302 87 -1.611 0.2466
## PMT Biskuit - PMT Lokal -0.1233 0.302 87 -0.408 0.9123
##
## P value adjustment: tukey method for comparing a family of 3 estimates# Holm sebagai alternatif koreksi
pairs(em2b, adjust = "holm")## Bulan = Bulan_1:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit -0.1433 0.298 87 -0.481 1.0000
## Kontrol - PMT Lokal -0.1900 0.298 87 -0.638 1.0000
## PMT Biskuit - PMT Lokal -0.0467 0.298 87 -0.157 1.0000
##
## Bulan = Bulan_2:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit -0.2467 0.299 87 -0.824 0.9735
## Kontrol - PMT Lokal -0.2967 0.299 87 -0.991 0.9735
## PMT Biskuit - PMT Lokal -0.0500 0.299 87 -0.167 0.9735
##
## Bulan = Bulan_3:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit -0.3633 0.302 87 -1.203 0.4648
## Kontrol - PMT Lokal -0.4867 0.302 87 -1.611 0.3326
## PMT Biskuit - PMT Lokal -0.1233 0.302 87 -0.408 0.6841
##
## P value adjustment: holm method for 3 tests# 4H. KONTRAS PERUBAHAN BULAN 1 -> BULAN 3
em_full <- emmeans(aov2, ~ Bulan * Kelompok)
em_full## Bulan Kelompok emmean SE df lower.CL upper.CL
## Bulan_1 Kontrol 9.74 0.211 87 9.32 10.2
## Bulan_2 Kontrol 9.96 0.212 87 9.54 10.4
## Bulan_3 Kontrol 10.18 0.214 87 9.76 10.6
## Bulan_1 PMT Biskuit 9.88 0.211 87 9.46 10.3
## Bulan_2 PMT Biskuit 10.20 0.212 87 9.78 10.6
## Bulan_3 PMT Biskuit 10.54 0.214 87 10.12 11.0
## Bulan_1 PMT Lokal 9.93 0.211 87 9.51 10.3
## Bulan_2 PMT Lokal 10.25 0.212 87 9.83 10.7
## Bulan_3 PMT Lokal 10.67 0.214 87 10.24 11.1
##
## Confidence level used: 0.95# Perbedaan perubahan Bulan 3 - Bulan 1 antar kelompok
contrast(
em_full,
interaction = list(
Bulan = list("Bulan3-Bulan1" = c(-1, 0, 1)),
Kelompok = "pairwise"
),
adjust = "holm"
)## Bulan_custom Kelompok_pairwise estimate SE df t.ratio p.value
## Bulan3-Bulan1 Kontrol - PMT Biskuit -0.2200 0.0482 87 -4.562 <0.0001
## Bulan3-Bulan1 Kontrol - PMT Lokal -0.2967 0.0482 87 -6.151 <0.0001
## Bulan3-Bulan1 PMT Biskuit - PMT Lokal -0.0767 0.0482 87 -1.590 0.1155
##
## P value adjustment: holm method for 3 tests# 4I. TREN LINEAR ANTARKELOMPOK
tren_int <- summary(
contrast(
em_full,
interaction = c(
Bulan = "poly",
Kelompok = "pairwise"
),
adjust = "none"
)
)
tren_int# Nama kolom dapat berbeda menurut versi emmeans.
# Cek names(tren_int) terlebih dahulu.
names(tren_int)## [1] "Bulan_poly" "Kelompok_pairwise" "estimate"
## [4] "SE" "df" "t.ratio"
## [7] "p.value"if ("Bulan_poly" %in% names(tren_int)) {
tren_lin <- subset(tren_int, Bulan_poly == "linear")
} else if ("Bulan.poly" %in% names(tren_int)) {
tren_lin <- subset(tren_int, Bulan.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_lin5. LINEAR MIXED MODEL (LMM)
# Model random intercept
lmm1 <- lmer(
Berat_Badan_kg ~ Kelompok * Bulan + (1 | ID),
data = data_long,
REML = TRUE
)
summary(lmm1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Berat_Badan_kg ~ Kelompok * Bulan + (1 | ID)
## Data: data_long
##
## REML criterion at convergence: 113.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.49430 -0.54808 0.01916 0.58303 2.09479
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.33772 1.1566
## Residual 0.01057 0.1028
## Number of obs: 270, groups: ID, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.149630 0.122077 86.999999 83.142 < 2e-16 ***
## Kelompok1 -0.191852 0.172642 86.999999 -1.111 0.2695
## Kelompok2 0.059259 0.172642 86.999999 0.343 0.7322
## Bulan1 -0.301852 0.008847 174.000000 -34.118 < 2e-16 ***
## Bulan2 -0.011852 0.008847 174.000000 -1.340 0.1821
## Kelompok1:Bulan1 0.080741 0.012512 174.000000 6.453 1.05e-09 ***
## Kelompok2:Bulan1 -0.027037 0.012512 174.000000 -2.161 0.0321 *
## Kelompok1:Bulan2 0.010741 0.012512 174.000000 0.858 0.3918
## Kelompok2:Bulan2 0.006296 0.012512 174.000000 0.503 0.6154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Klmpk1 Klmpk2 Bulan1 Bulan2 Kl1:B1 Kl2:B1 Kl1:B2
## Kelompok1 0.000
## Kelompok2 0.000 -0.500
## Bulan1 0.000 0.000 0.000
## Bulan2 0.000 0.000 0.000 -0.500
## Klmpk1:Bln1 0.000 0.000 0.000 0.000 0.000
## Klmpk2:Bln1 0.000 0.000 0.000 0.000 0.000 -0.500
## Klmpk1:Bln2 0.000 0.000 0.000 0.000 0.000 -0.500 0.250
## Klmpk2:Bln2 0.000 0.000 0.000 0.000 0.000 0.250 -0.500 -0.500# Model random intercept + random slope waktu numerik
lmm2 <- lmer(
Berat_Badan_kg ~ Kelompok * Bulan_num + (1 + Bulan_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: Berat_Badan_kg ~ Kelompok * Bulan_num + (1 + Bulan_num | ID)
## Data: data_long
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 50.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.13804 -0.43265 0.00516 0.48912 1.56640
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.326812 1.15187
## Bulan_num 0.006722 0.08199 -0.04
## Residual 0.004000 0.06325
## Number of obs: 270, groups: ID, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.534074 0.121844 87.000058 78.248 < 2e-16 ***
## Kelompok1 -0.019630 0.172314 87.000070 -0.114 0.9096
## Kelompok2 0.011481 0.172314 87.000070 0.067 0.9470
## Bulan_num 0.307778 0.009844 87.000016 31.264 < 2e-16 ***
## Kelompok1:Bulan_num -0.086111 0.013922 87.000012 -6.185 1.96e-08 ***
## Kelompok2:Bulan_num 0.023889 0.013922 87.000012 1.716 0.0897 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Klmpk1 Klmpk2 Bln_nm Kl1:B_
## Kelompok1 0.000
## Kelompok2 0.000 -0.500
## Bulan_num -0.069 0.000 0.000
## Klmpk1:Bln_ 0.000 -0.069 0.034 0.000
## Klmpk2:Bln_ 0.000 0.034 -0.069 0.000 -0.500# Bandingkan struktur random effect
anova(lmm1, lmm2, refit = FALSE)# Uji efek tetap
anova(lmm2, ddf = "Kenward-Roger")# ICC
performance::icc(lmm1)# 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.97828, p-value = 0.00038916. EXPORT HASIL
# Statistik deskriptif
write.csv(
deskriptif,
"hasil_01_statistik_deskriptif_BB.csv",
row.names = FALSE
)
# Deskriptif usia
write.csv(
deskriptif_usia,
"hasil_02_statistik_usia.csv",
row.names = FALSE
)
# Distribusi jenis kelamin
write.csv(
tabel_jenis_kelamin,
"hasil_03_jenis_kelamin.csv",
row.names = FALSE
)
# Normalitas
write.csv(
shapiro_mixed,
"hasil_04_shapiro_Group_x_Bulan.csv",
row.names = FALSE
)
# Levene
write.csv(
levene_mixed,
"hasil_05_levene_per_Bulan.csv",
row.names = FALSE
)
# ANOVA utama
hasil_anova <- as.data.frame(
nice(
aov2,
correction = "GG",
es = c("ges", "pes")
)
)
write.csv(
hasil_anova,
"hasil_06_mixed_ANOVA_GG.csv",
row.names = FALSE
)
# Post-hoc antarkelompok
posthoc_group <- as.data.frame(
pairs(em2b, adjust = "holm")
)
write.csv(
posthoc_group,
"hasil_07_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_08_posthoc_dalam_kelompok_Holm.csv",
row.names = FALSE
)
# Estimated marginal means
emm_group_time <- as.data.frame(em2b)
write.csv(
emm_group_time,
"hasil_09_estimated_marginal_means.csv",
row.names = FALSE
)
# Kontras perubahan Bulan 1 -> Bulan 3
kontras_bulan1_bulan3 <- as.data.frame(
contrast(
em_full,
interaction = list(
Bulan = list("Bulan3-Bulan1" = c(-1, 0, 1)),
Kelompok = "pairwise"
),
adjust = "holm"
)
)
write.csv(
kontras_bulan1_bulan3,
"hasil_10_kontras_perubahan_Bulan1_Bulan3.csv",
row.names = FALSE
)
# Gambar
ggsave(
"grafik_01_profile_Berat_Badan.png",
p_profil,
width = 8,
height = 5.5,
dpi = 300
)
ggsave(
"grafik_02_spaghetti_Berat_Badan.png",
p_spag,
width = 9,
height = 6,
dpi = 300
)
ggsave(
"grafik_03_interaksi_Kelompok_x_Bulan.png",
p_interaksi,
width = 8,
height = 5.5,
dpi = 300
)7. RINGKASAN AKHIR DI CONSOLE
cat("\n============================================================\n")##
## ============================================================cat("ANALISIS REPEATED MEASURES DATA BALITA SELESAI\n")## ANALISIS REPEATED MEASURES DATA BALITA SELESAIcat("============================================================\n")## ============================================================cat("Jumlah responden :", n_distinct(data_long$ID), "\n")## Jumlah responden : 90cat("Kelompok :", n_distinct(data_long$Kelompok), "\n")## Kelompok : 3cat("Waktu pengukuran :", n_distinct(data_long$Bulan), "\n")## Waktu pengukuran : 3cat("Outcome : Berat badan (kg)\n")## Outcome : Berat badan (kg)cat("Kelompok : Kontrol, PMT Biskuit, PMT Lokal\n")## Kelompok : Kontrol, PMT Biskuit, PMT Lokalcat("\nSemua file hasil tersimpan di:\n")##
## Semua file hasil tersimpan di:cat(getwd(), "\n")## C:/Users/Hype AMD/Downloads/New foldercat("============================================================\n")## ============================================================# =============================================================================
# SELESAI
# =============================================================================Catatan Interpretasi
Hasil statistik dan grafik akan dihitung ulang saat dokumen di-knit. Nilai pada laporan HTML/Word sumber tidak ditanam sebagai output statis, sehingga laporan R Markdown tetap dapat direproduksi dari data Excel yang sama.