# =============================================================================
# Nama : [Sundari]
# NIM : [2611018024]
# =============================================================================
# REPEATED MEASURES ANALYSIS DENGAN R
# Topik : Perubahan berat badan balita pada 3 kelompok selama 3 bulan
# (Kontrol, PMT Biskuit, PMT Lokal)
# DATA : Data_PMT_Balita_-long.xlsx (sheet Data_Long)
#
# Struktur data:
# - 90 balita (30 per kelompok)
# - 3 kelompok: Kontrol, PMT Biskuit, PMT Lokal
# - 4 kali penimbangan: bulan 0, 1, 2, 3
# - Outcome : BB_kg (berat badan, kg)
# - Kovariat : Usia_bulan (12-59 bulan), Jenis_Kelamin
#
# Alur analisis mengikuti format syntax terlampir:
# 0. Paket & pengaturan
# 1. Membaca data Excel (LONG), membentuk WIDE + validasi
# 2. Eksplorasi data + grafik
# 3. Repeated Measures ANOVA satu arah (contoh PMT Lokal)
# 3a. Outlier, normalitas, sphericity/Mauchly
# 3b. ANOVA + Greenhouse-Geisser/Huynh-Feldt
# 3c. Pendekatan multivariat
# 3d. Post-hoc dan tren
# 3e. Friedman
# 4. Mixed Design ANOVA (Kelompok x Waktu)
# 4a. Asumsi: outlier, normalitas, Levene, Box's M, Mauchly
# 4b. ANOVA + effect size
# 4c-4g. Ukuran efek, plot interaksi, simple effects, post-hoc
# 4h. Kontras perubahan bulan 0 -> bulan 3
# 4i. Tren linear antarkelompok
# 4j. Pembanding: selisih BB dan ANCOVA
# 5. Linear Mixed Model (LMM)
# 6. Opsional: simulasi missing value
# 7. Export hasil
# 8. Ringkasan akhir
# =============================================================================
# 0. PAKET & PENGATURAN
# Paket yang belum terpasang akan dipasang otomatis (perlu internet).
paket <- c("readxl", "dplyr", "tidyr", "ggplot2", "afex", "emmeans",
"rstatix", "car", "effectsize", "lme4", "lmerTest", "pbkrtest",
"performance", "ggpubr", "Hmisc")
baru <- paket[!paket %in% rownames(installed.packages())]
if (length(baru) > 0) install.packages(baru)
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 MEMBENTUK WIDE
# Pastikan file Excel berada di Working Directory.
# Cek dengan:
getwd()
## [1] "C:/Users/ASUS/Downloads/ndari R studio"
list.files()
## [1] "Data_PMT_Balita(2).xlsx"
## [2] "dataset_PMT_balita_LONG.csv"
## [3] "dataset_PMT_balita_WIDE.csv"
## [4] "grafik_01_profile_BB.png"
## [5] "grafik_02_spaghetti_BB.png"
## [6] "grafik_03_interaksi_Kelompok_x_Waktu.png"
## [7] "hasil_01_statistik_deskriptif_BB.csv"
## [8] "hasil_02_shapiro_Kelompok_x_Waktu.csv"
## [9] "hasil_03_levene_per_Waktu.csv"
## [10] "hasil_04_mixed_ANOVA_GG.csv"
## [11] "hasil_05_posthoc_antar_kelompok_Holm.csv"
## [12] "hasil_06_posthoc_dalam_kelompok_Holm.csv"
## [13] "hasil_07_estimated_marginal_means.csv"
## [14] "RStudio_Repeated_Measures_BB_PMT_Balita.docx"
## [15] "RStudio_Repeated_Measures_BB_PMT_Balita.html"
## [16] "RStudio_Repeated_Measures_BB_PMT_Balita.R"
## [17] "RStudio_Repeated_Measures_BB_PMT_Balita.spin.R"
## [18] "RStudio_Repeated_Measures_BB_PMT_Balita.spin.Rmd"
# Jika file belum ditemukan, atur Working Directory melalui:
# Session -> Set Working Directory -> Choose Directory...
# (atau: Session -> Set Working Directory -> To Source File Location)
# 1A. Baca data LONG dari Excel
# File dicari otomatis (nama boleh berawalan angka, mis.
# "1791007307248_Data_PMT_Balita_-long.xlsx"). Jika tidak ada, muncul
# jendela untuk memilih file secara manual.
kandidat <- list.files(pattern = "Data_PMT_Balita.*\\.xlsx$", ignore.case = TRUE)
kandidat <- kandidat[!startsWith(kandidat, "~$")] # abaikan file sementara Excel
file_data <- if (length(kandidat) > 0) kandidat[1] else file.choose()
message("Membaca file: ", file_data)
## Membaca file: Data_PMT_Balita(2).xlsx
data_long <- as.data.frame(read_excel(file_data, sheet = "Data_Long"))
# 1B. Bentuk WIDE dari LONG
# File Excel hanya berisi format LONG, sehingga WIDE dibentuk dari LONG.
# (Jika file Excel juga memiliki sheet "Data_Wide", sheet itu dipakai
# pada langkah 1E untuk pengecekan konsistensi.)
data_wide <- data_long %>%
mutate(Waktu_lab = paste0("B", Waktu)) %>%
dplyr::select(ID, Kelompok, Usia_bulan, Jenis_Kelamin, Waktu_lab, BB_kg) %>%
pivot_wider(
names_from = Waktu_lab,
values_from = BB_kg,
names_prefix = "BB_"
) %>%
as.data.frame()
# Tampilkan data
View(data_long)
View(data_wide)
# Struktur dan dimensi
head(data_long)
## ID Kelompok Usia_bulan Jenis_Kelamin Waktu BB_kg
## 1 BSK-01 PMT Biskuit 24 Laki-laki 0 10.2
## 2 BSK-01 PMT Biskuit 24 Laki-laki 1 10.5
## 3 BSK-01 PMT Biskuit 24 Laki-laki 2 10.7
## 4 BSK-01 PMT Biskuit 24 Laki-laki 3 11.1
## 5 BSK-02 PMT Biskuit 19 Laki-laki 0 9.1
## 6 BSK-02 PMT Biskuit 19 Laki-laki 1 9.1
head(data_wide)
## ID Kelompok Usia_bulan Jenis_Kelamin BB_B0 BB_B1 BB_B2 BB_B3
## 1 BSK-01 PMT Biskuit 24 Laki-laki 10.2 10.5 10.7 11.1
## 2 BSK-02 PMT Biskuit 19 Laki-laki 9.1 9.1 9.5 9.9
## 3 BSK-03 PMT Biskuit 39 Laki-laki 15.6 16.0 16.2 16.8
## 4 BSK-04 PMT Biskuit 33 Perempuan 12.2 12.6 12.8 13.3
## 5 BSK-05 PMT Biskuit 25 Laki-laki 11.1 11.3 11.8 12.0
## 6 BSK-06 PMT Biskuit 22 Perempuan 10.9 11.1 11.4 11.9
str(data_long)
## 'data.frame': 360 obs. of 6 variables:
## $ ID : chr "BSK-01" "BSK-01" "BSK-01" "BSK-01" ...
## $ Kelompok : chr "PMT Biskuit" "PMT Biskuit" "PMT Biskuit" "PMT Biskuit" ...
## $ Usia_bulan : num 24 24 24 24 19 19 19 19 39 39 ...
## $ Jenis_Kelamin: chr "Laki-laki" "Laki-laki" "Laki-laki" "Laki-laki" ...
## $ Waktu : num 0 1 2 3 0 1 2 3 0 1 ...
## $ BB_kg : num 10.2 10.5 10.7 11.1 9.1 9.1 9.5 9.9 15.6 16 ...
str(data_wide)
## 'data.frame': 90 obs. of 8 variables:
## $ ID : chr "BSK-01" "BSK-02" "BSK-03" "BSK-04" ...
## $ Kelompok : chr "PMT Biskuit" "PMT Biskuit" "PMT Biskuit" "PMT Biskuit" ...
## $ Usia_bulan : num 24 19 39 33 25 22 45 12 46 37 ...
## $ Jenis_Kelamin: chr "Laki-laki" "Laki-laki" "Laki-laki" "Perempuan" ...
## $ BB_B0 : num 10.2 9.1 15.6 12.2 11.1 10.9 14.5 9 13.4 11.3 ...
## $ BB_B1 : num 10.5 9.1 16 12.6 11.3 11.1 14.9 9.3 13.8 11.6 ...
## $ BB_B2 : num 10.7 9.5 16.2 12.8 11.8 11.4 15.1 9.2 14.1 11.9 ...
## $ BB_B3 : num 11.1 9.9 16.8 13.3 12 11.9 15.5 9.6 14.4 12.2 ...
dim(data_long)
## [1] 360 6
dim(data_wide)
## [1] 90 8
names(data_long)
## [1] "ID" "Kelompok" "Usia_bulan" "Jenis_Kelamin"
## [5] "Waktu" "BB_kg"
names(data_wide)
## [1] "ID" "Kelompok" "Usia_bulan" "Jenis_Kelamin"
## [5] "BB_B0" "BB_B1" "BB_B2" "BB_B3"
# 1C. VALIDASI STRUKTUR DATA
# Kolom yang diharapkan pada LONG
stopifnot(
all(c("ID", "Kelompok", "Usia_bulan", "Jenis_Kelamin", "Waktu", "BB_kg")
%in% names(data_long))
)
# Kolom yang diharapkan pada WIDE
stopifnot(
all(c("ID", "Kelompok", "BB_B0", "BB_B1", "BB_B2", "BB_B3")
%in% names(data_wide))
)
# Pastikan jumlah balita 90
n_distinct(data_long$ID)
## [1] 90
n_distinct(data_wide$ID)
## [1] 90
# Jumlah balita tiap kelompok
data_wide %>% count(Kelompok)
## Kelompok n
## 1 Kontrol 30
## 2 PMT Biskuit 30
## 3 PMT Lokal 30
data_long %>%
distinct(ID, Kelompok) %>%
count(Kelompok)
## Kelompok n
## 1 Kontrol 30
## 2 PMT Biskuit 30
## 3 PMT Lokal 30
# Jumlah penimbangan tiap balita (harus 4)
data_long %>%
count(ID) %>%
count(n, name = "jumlah_balita")
## n jumlah_balita
## 1 4 90
# Cek duplikasi ID x Waktu
cek_duplikasi <- data_long %>%
count(ID, Waktu) %>%
filter(n != 1)
cek_duplikasi
## [1] ID Waktu n
## <0 rows> (or 0-length row.names)
# Cek missing value
colSums(is.na(data_long))
## ID Kelompok Usia_bulan Jenis_Kelamin Waktu
## 0 0 0 0 0
## BB_kg
## 0
colSums(is.na(data_wide))
## ID Kelompok Usia_bulan Jenis_Kelamin BB_B0
## 0 0 0 0 0
## BB_B1 BB_B2 BB_B3
## 0 0 0
# Ringkasan outcome
summary(data_long$BB_kg)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.80 10.80 12.80 12.88 14.72 18.30
# 1D. RECODING VARIABEL
# LONG digunakan sebagai basis analisis.
# Waktu_num tetap numerik untuk grafik/LMM,
# sedangkan Waktu menjadi faktor untuk repeated-measures ANOVA.
kel_lab <- c("Kontrol", "PMT Biskuit", "PMT Lokal")
data_long <- data_long %>%
mutate(
ID = factor(ID),
Kelompok = factor(Kelompok, levels = kel_lab),
Jenis_Kelamin = factor(Jenis_Kelamin),
Waktu_num = as.numeric(Waktu),
Waktu = factor(
Waktu,
levels = c(0, 1, 2, 3),
labels = c("B0", "B1", "B2", "B3")
)
)
# WIDE juga diberi factor untuk keperluan Box's M / pengecekan.
data_wide <- data_wide %>%
mutate(
ID = factor(ID),
Kelompok = factor(Kelompok, levels = kel_lab),
Jenis_Kelamin = factor(Jenis_Kelamin)
)
# Cek kembali level
levels(data_long$Kelompok)
## [1] "Kontrol" "PMT Biskuit" "PMT Lokal"
levels(data_long$Waktu)
## [1] "B0" "B1" "B2" "B3"
# 1E. CEK KONSISTENSI LONG VS WIDE
# Buat WIDE dari LONG untuk membandingkan dengan WIDE.
wide_from_long <- data_long %>%
dplyr::select(ID, Kelompok, Waktu, BB_kg) %>%
pivot_wider(
names_from = Waktu,
values_from = BB_kg,
names_prefix = "BB_"
)
# Bila ingin melihat hasil transformasi:
View(wide_from_long)
# Pengecekan sederhana jumlah baris
nrow(wide_from_long)
## [1] 90
nrow(data_wide)
## [1] 90
# Bila file Excel memiliki sheet "Data_Wide", cocokkan nilainya
if ("Data_Wide" %in% excel_sheets(file_data)) {
wide_file <- read_excel(file_data, sheet = "Data_Wide")
kol_file <- c("BB_Bulan0", "BB_Bulan1", "BB_Bulan2", "BB_Bulan3")
urut <- match(as.character(data_wide$ID), wide_file$ID)
print(all.equal(
as.matrix(wide_file[urut, kol_file]),
as.matrix(data_wide[, c("BB_B0", "BB_B1", "BB_B2", "BB_B3")]),
check.attributes = FALSE
))
}
## [1] TRUE
# 1F. KARAKTERISTIK AWAL (kesetaraan kelompok)
karakteristik <- data_wide %>%
group_by(Kelompok) %>%
summarise(
n = n(),
usia_mean = mean(Usia_bulan),
usia_sd = sd(Usia_bulan),
laki = sum(Jenis_Kelamin == "Laki-laki"),
perempuan = sum(Jenis_Kelamin == "Perempuan"),
bb0_mean = mean(BB_B0),
bb0_sd = sd(BB_B0),
.groups = "drop"
)
karakteristik
## # A tibble: 3 × 8
## Kelompok n usia_mean usia_sd laki perempuan bb0_mean bb0_sd
## <fct> <int> <dbl> <dbl> <int> <int> <dbl> <dbl>
## 1 Kontrol 30 41.9 13.9 12 18 13.0 2.36
## 2 PMT Biskuit 30 34.5 13.6 19 11 12.1 2.39
## 3 PMT Lokal 30 35.9 14.3 19 11 12.3 2.47
chisq.test(table(data_wide$Kelompok, data_wide$Jenis_Kelamin)) # jenis kelamin
##
## Pearson's Chi-squared test
##
## data: table(data_wide$Kelompok, data_wide$Jenis_Kelamin)
## X-squared = 4.41, df = 2, p-value = 0.1103
summary(aov(Usia_bulan ~ Kelompok, data = data_wide)) # usia
## Df Sum Sq Mean Sq F value Pr(>F)
## Kelompok 2 927 463.6 2.39 0.0976 .
## Residuals 87 16875 194.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(BB_B0 ~ Kelompok, data = data_wide)) # BB awal
## Df Sum Sq Mean Sq F value Pr(>F)
## Kelompok 2 11.6 5.802 1 0.372
## Residuals 87 504.6 5.800
# 2. EKSPLORASI DATA
# 2A. Statistik deskriptif per kelompok dan waktu
deskriptif <- data_long %>%
group_by(Kelompok, Waktu) %>%
summarise(
n = n(),
mean = mean(BB_kg, na.rm = TRUE),
sd = sd(BB_kg, na.rm = TRUE),
median = median(BB_kg, na.rm = TRUE),
min = min(BB_kg, na.rm = TRUE),
max = max(BB_kg, na.rm = TRUE),
.groups = "drop"
)
deskriptif
## # A tibble: 12 × 8
## Kelompok Waktu n mean sd median min max
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kontrol B0 30 13.0 2.36 13 8.3 17.4
## 2 Kontrol B1 30 13.1 2.37 13.2 8.3 17.6
## 3 Kontrol B2 30 13.2 2.40 13.2 8.6 17.8
## 4 Kontrol B3 30 13.3 2.42 13.4 8.8 17.8
## 5 PMT Biskuit B0 30 12.1 2.39 12.0 7.8 16
## 6 PMT Biskuit B1 30 12.5 2.46 12.4 7.9 16.7
## 7 PMT Biskuit B2 30 12.7 2.47 12.6 8.3 17.1
## 8 PMT Biskuit B3 30 13.1 2.47 13 8.6 17.4
## 9 PMT Lokal B0 30 12.3 2.47 12.2 8.3 17.3
## 10 PMT Lokal B1 30 12.7 2.43 12.6 8.7 17.5
## 11 PMT Lokal B2 30 13.1 2.44 12.7 9.2 17.9
## 12 PMT Lokal B3 30 13.4 2.44 13.0 9.6 18.3
# Versi rstatix
data_long %>%
group_by(Kelompok, Waktu) %>%
get_summary_stats(BB_kg, type = "mean_sd")
## # A tibble: 12 × 6
## Kelompok Waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol B0 BB_kg 30 13.0 2.36
## 2 Kontrol B1 BB_kg 30 13.1 2.37
## 3 Kontrol B2 BB_kg 30 13.2 2.40
## 4 Kontrol B3 BB_kg 30 13.3 2.42
## 5 PMT Biskuit B0 BB_kg 30 12.1 2.39
## 6 PMT Biskuit B1 BB_kg 30 12.5 2.46
## 7 PMT Biskuit B2 BB_kg 30 12.7 2.47
## 8 PMT Biskuit B3 BB_kg 30 13.1 2.47
## 9 PMT Lokal B0 BB_kg 30 12.3 2.47
## 10 PMT Lokal B1 BB_kg 30 12.7 2.43
## 11 PMT Lokal B2 BB_kg 30 13.1 2.44
## 12 PMT Lokal B3 BB_kg 30 13.4 2.44
# 2B. Matriks kovarians dan korelasi antar waktu
vars_waktu <- c("BB_B0", "BB_B1", "BB_B2", "BB_B3")
S <- cov(data_wide[, vars_waktu], use = "complete.obs")
R <- cor(data_wide[, vars_waktu], use = "complete.obs")
round(S, 2)
## BB_B0 BB_B1 BB_B2 BB_B3
## BB_B0 5.80 5.79 5.79 5.75
## BB_B1 5.79 5.80 5.81 5.78
## BB_B2 5.79 5.81 5.86 5.84
## BB_B3 5.75 5.78 5.84 5.85
round(R, 3)
## BB_B0 BB_B1 BB_B2 BB_B3
## BB_B0 1.000 0.998 0.993 0.986
## BB_B1 0.998 1.000 0.998 0.993
## BB_B2 0.993 0.998 1.000 0.997
## BB_B3 0.986 0.993 0.997 1.000
# 2C. Varians selisih antar waktu
pasangan <- combn(vars_waktu, 2)
var_selisih <- apply(
pasangan,
2,
function(p) var(data_wide[[p[1]]] - data_wide[[p[2]]], na.rm = TRUE)
)
names(var_selisih) <- apply(
pasangan,
2,
paste,
collapse = " - "
)
round(var_selisih, 3)
## BB_B0 - BB_B1 BB_B0 - BB_B2 BB_B0 - BB_B3 BB_B1 - BB_B2 BB_B1 - BB_B3
## 0.027 0.078 0.158 0.028 0.087
## BB_B2 - BB_B3
## 0.031
# 2D. Profile plot: rerata +/- 95% CI
p_profil <- ggplot(
data_long,
aes(
x = Waktu_num,
y = BB_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 = .12
) +
scale_x_continuous(breaks = c(0, 1, 2, 3)) +
labs(
x = "Bulan penimbangan",
y = "Berat badan (kg)",
colour = "Kelompok",
title = "Profil rerata berat badan balita selama 3 bulan"
) +
theme(legend.position = "bottom")
p_profil
# 2E. Spaghetti plot
p_spag <- ggplot(
data_long,
aes(
x = Waktu_num,
y = BB_kg,
group = ID
)
) +
geom_line(alpha = .25) +
stat_summary(
aes(group = Kelompok),
fun = mean,
geom = "line",
linewidth = 1.3,
colour = "firebrick"
) +
facet_wrap(~ Kelompok) +
scale_x_continuous(breaks = c(0, 1, 2, 3)) +
labs(
x = "Bulan",
y = "Berat badan (kg)",
title = "Lintasan individu dan rerata kelompok"
)
p_spag
# 3. REPEATED MEASURES ANOVA SATU ARAH
# Contoh: kelompok PMT Lokal
# Pertanyaan:
# Apakah berat badan berubah selama 3 bulan pada kelompok PMT Lokal?
# (ganti kelompok_fokus dengan "Kontrol" atau "PMT Biskuit" bila diperlukan)
kelompok_fokus <- "PMT Lokal"
d1 <- data_long %>%
filter(Kelompok == kelompok_fokus) %>%
droplevels()
d1w <- data_wide %>%
filter(Kelompok == kelompok_fokus) %>%
droplevels()
# 3A. UJI ASUMSI
# (i) Outlier per waktu
outlier_1way <- d1 %>%
group_by(Waktu) %>%
identify_outliers(BB_kg)
outlier_1way
## [1] Waktu ID Kelompok Usia_bulan Jenis_Kelamin
## [6] BB_kg Waktu_num is.outlier is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas Shapiro-Wilk per waktu
shapiro_1way <- d1 %>%
group_by(Waktu) %>%
shapiro_test(BB_kg)
shapiro_1way
## # A tibble: 4 × 4
## Waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 B0 BB_kg 0.961 0.328
## 2 B1 BB_kg 0.962 0.353
## 3 B2 BB_kg 0.958 0.275
## 4 B3 BB_kg 0.955 0.232
# Q-Q plot
qqplot_1way <- ggpubr::ggqqplot(
d1,
"BB_kg",
facet.by = "Waktu"
)
qqplot_1way
# (iii) Mauchly's Test
# anova_test melaporkan Mauchly + GG/HF
aov1_rs <- anova_test(
data = d1,
dv = BB_kg,
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 87 602.823 4.51e-58 * 0.954
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 Waktu 0.323 8.07e-06 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 Waktu 0.637 1.91, 55.4 7.14e-38 * 0.681 2.04, 59.21 2.6e-40
## 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.91 55.4 602.823 7.14e-38 * 0.954
# 3B. REPEATED MEASURES ANOVA DENGAN AFEX
aov1 <- aov_ez(
id = "ID",
dv = "BB_kg",
data = d1,
within = "Waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov1
## Anova Table (Type 3 tests)
##
## Response: BB_kg
## Effect df MSE F ges pes p.value
## 1 Waktu 1.91, 55.40 0.02 602.82 *** .030 .954 <.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) 19902.2 1 693.01 29 832.83 < 2.2e-16 ***
## Waktu 21.2 3 1.02 87 602.82 < 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.32253 8.0703e-06
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Waktu 0.63682 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Waktu 0.6806248 2.602525e-40
nice(aov1, correction = "GG", es = c("ges", "pes"))
## Anova Table (Type 3 tests)
##
## Response: BB_kg
## Effect df MSE F ges pes p.value
## 1 Waktu 1.91, 55.40 0.02 602.82 *** .030 .954 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## Waktu | 0.95 | [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.03 | [0.00, 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.96635 832.83 1 29 < 2.2e-16 ***
## Waktu 1 0.97729 387.22 3 27 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3D. POST-HOC DAN TREN
em1 <- emmeans(aov1, ~ Waktu)
em1
## Waktu emmean SE df lower.CL upper.CL
## B0 12.3 0.452 29 11.4 13.2
## B1 12.7 0.443 29 11.8 13.6
## B2 13.1 0.446 29 12.2 14.0
## B3 13.4 0.445 29 12.5 14.4
##
## Confidence level used: 0.95
# Semua pasangan waktu
pairs(em1, adjust = "holm")
## contrast estimate SE df t.ratio p.value
## B0 - B1 -0.363 0.0200 29 -18.123 <0.0001
## B0 - B2 -0.743 0.0341 29 -21.777 <0.0001
## B0 - B3 -1.127 0.0346 29 -32.607 <0.0001
## B1 - B2 -0.380 0.0232 29 -16.384 <0.0001
## B1 - B3 -0.763 0.0305 29 -25.022 <0.0001
## B2 - B3 -0.383 0.0215 29 -17.840 <0.0001
##
## P value adjustment: holm method for 6 tests
# Masing-masing waktu dibandingkan dengan baseline B0
contrast(
em1,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## contrast estimate SE df t.ratio p.value
## B1 - B0 0.363 0.0200 29 18.123 <0.0001
## B2 - B0 0.743 0.0341 29 21.777 <0.0001
## B3 - B0 1.127 0.0346 29 32.607 <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 3.7600 0.1220 29 30.721 <0.0001
## quadratic 0.0200 0.0350 29 0.571 0.5725
## cubic -0.0133 0.0484 29 -0.276 0.7847
# 3E. ALTERNATIF NONPARAMETRIK: FRIEDMAN
friedman_test(d1, BB_kg ~ Waktu | ID)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 BB_kg 30 90 3 2.19e-19 Friedman test
friedman_effsize(d1, BB_kg ~ Waktu | ID)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 BB_kg 30 1 Kendall W large
# Wilcoxon berpasangan sebagai post-hoc alternatif
# bila diperlukan
d1 %>%
wilcox_test(
BB_kg ~ Waktu,
paired = TRUE,
p.adjust.method = "holm"
)
## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 BB_kg B0 B1 30 30 0 0.00000000186 1.12e-8 ****
## 2 BB_kg B0 B2 30 30 0 0.00000000186 1.12e-8 ****
## 3 BB_kg B0 B3 30 30 0 0.00000000186 1.12e-8 ****
## 4 BB_kg B1 B2 30 30 0 0.00000000186 1.12e-8 ****
## 5 BB_kg B1 B3 30 30 0 0.00000000186 1.12e-8 ****
## 6 BB_kg B2 B3 30 30 0 0.00000000186 1.12e-8 ****
# 4. MIXED DESIGN ANOVA
# Kelompok (between) x Waktu (within)
# Pertanyaan utama:
# Apakah perubahan berat badan dari bulan 0 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, Waktu) %>%
identify_outliers(BB_kg)
outlier_mixed
## [1] Kelompok Waktu ID Usia_bulan Jenis_Kelamin
## [6] BB_kg Waktu_num is.outlier is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas Shapiro-Wilk per Kelompok x Waktu
shapiro_mixed <- data_long %>%
group_by(Kelompok, Waktu) %>%
shapiro_test(BB_kg)
shapiro_mixed
## # A tibble: 12 × 5
## Kelompok Waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol B0 BB_kg 0.958 0.270
## 2 Kontrol B1 BB_kg 0.960 0.315
## 3 Kontrol B2 BB_kg 0.959 0.291
## 4 Kontrol B3 BB_kg 0.956 0.243
## 5 PMT Biskuit B0 BB_kg 0.957 0.260
## 6 PMT Biskuit B1 BB_kg 0.965 0.409
## 7 PMT Biskuit B2 BB_kg 0.964 0.401
## 8 PMT Biskuit B3 BB_kg 0.964 0.386
## 9 PMT Lokal B0 BB_kg 0.961 0.328
## 10 PMT Lokal B1 BB_kg 0.962 0.353
## 11 PMT Lokal B2 BB_kg 0.958 0.275
## 12 PMT Lokal B3 BB_kg 0.955 0.232
# Q-Q plot per sel
qqplot_mixed <- ggpubr::ggqqplot(
data_long,
"BB_kg"
) +
facet_grid(Waktu ~ Kelompok)
qqplot_mixed
# (iii) Homogenitas varians pada setiap waktu: Levene
levene_mixed <- data_long %>%
group_by(Waktu) %>%
levene_test(BB_kg ~ Kelompok)
levene_mixed
## # A tibble: 4 × 5
## Waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 B0 2 87 0.169 0.845
## 2 B1 2 87 0.0557 0.946
## 3 B2 2 87 0.0514 0.950
## 4 B3 2 87 0.0383 0.962
# (iv) Homogenitas matriks kovarians: Box's M
box_m_result <- box_m(
data_wide[, vars_waktu],
data_wide$Kelompok
)
box_m_result
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 34.3 0.0245 20 Box's M-test for Homogeneity of Covariance Matric…
# (v) Sphericity: Mauchly dilihat pada summary(aov2)
# 4B. MIXED DESIGN REPEATED MEASURES ANOVA
aov2 <- aov_ez(
id = "ID",
dv = "BB_kg",
data = data_long,
between = "Kelompok",
within = "Waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov2
## Anova Table (Type 3 tests)
##
## Response: BB_kg
## Effect df MSE F ges pes p.value
## 1 Kelompok 2, 87 23.53 0.38 .009 .009 .687
## 2 Waktu 1.98, 172.38 0.02 803.10 *** .015 .902 <.001
## 3 Kelompok:Waktu 3.96, 172.38 0.02 70.39 *** .003 .618 <.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) 59678 1 2046.98 87 2536.4165 <2e-16 ***
## Kelompok 18 2 2046.98 87 0.3776 0.6866
## Waktu 32 3 3.48 261 803.0987 <2e-16 ***
## Kelompok:Waktu 6 6 3.48 261 70.3855 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## Waktu 0.42937 3.1623e-14
## Kelompok:Waktu 0.42937 3.1623e-14
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Waktu 0.66046 < 2.2e-16 ***
## Kelompok:Waktu 0.66046 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Waktu 0.6757837 9.223254e-90
## Kelompok:Waktu 0.6757837 8.423267e-36
nice(aov2, correction = "GG", es = c("ges", "pes"))
## Anova Table (Type 3 tests)
##
## Response: BB_kg
## Effect df MSE F ges pes p.value
## 1 Kelompok 2, 87 23.53 0.38 .009 .009 .687
## 2 Waktu 1.98, 172.38 0.02 803.10 *** .015 .902 <.001
## 3 Kelompok:Waktu 3.96, 172.38 0.02 70.39 *** .003 .618 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Alternatif rstatix
aov2_rs <- anova_test(
data = data_long,
dv = BB_kg,
wid = ID,
between = Kelompok,
within = Waktu,
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.378 6.87e-01 0.009
## 2 Waktu 3 261 803.099 1.97e-131 * 0.902
## 3 Kelompok:Waktu 6 261 70.386 9.55e-52 * 0.618
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 Waktu 0.429 3.16e-14 *
## 2 Kelompok:Waktu 0.429 3.16e-14 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF]
## 1 Waktu 0.66 1.98, 172.38 8.60e-88 * 0.676 2.03, 176.38
## 2 Kelompok:Waktu 0.66 3.96, 172.38 4.78e-35 * 0.676 4.05, 176.38
## p[HF] p[HF]<.05
## 1 9.22e-90 *
## 2 8.42e-36 *
get_anova_table(aov2_rs, correction = "GG")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 Kelompok 2.00 87.00 0.378 6.87e-01 0.009
## 2 Waktu 1.98 172.38 803.099 8.60e-88 * 0.902
## 3 Kelompok:Waktu 3.96 172.38 70.386 4.78e-35 * 0.618
# 4C. UKURAN EFEK
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## Kelompok | 8.61e-03 | [0.00, 1.00]
## Waktu | 0.90 | [0.89, 1.00]
## Kelompok:Waktu | 0.62 | [0.56, 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.00 | [0.00, 1.00]
## Waktu | 0.02 | [0.00, 1.00]
## Kelompok:Waktu | 2.67e-03 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 4D. PLOT INTERAKSI
p_interaksi <- afex_plot(
aov2,
x = "Waktu",
trace = "Kelompok",
error = "within",
mapping = c("colour", "shape", "linetype")
) +
labs(
y = "Berat badan (kg)",
x = "Waktu",
title = "Interaksi Kelompok x Waktu pada berat badan balita"
) +
theme(legend.position = "bottom")
## Warning: Panel(s) show a mixed within-between-design.
## Error bars do not allow comparisons across all means.
## Suppress error bars with: error = "none"
p_interaksi
# 4E. SIMPLE EFFECTS
# Estimated marginal means waktu dalam setiap kelompok
em2 <- emmeans(aov2, ~ Waktu | Kelompok)
em2
## Kelompok = Kontrol:
## Waktu emmean SE df lower.CL upper.CL
## B0 13.0 0.440 87 12.1 13.8
## B1 13.1 0.442 87 12.2 14.0
## B2 13.2 0.445 87 12.3 14.1
## B3 13.3 0.446 87 12.4 14.2
##
## Kelompok = PMT Biskuit:
## Waktu emmean SE df lower.CL upper.CL
## B0 12.1 0.440 87 11.3 13.0
## B1 12.5 0.442 87 11.6 13.3
## B2 12.7 0.445 87 11.8 13.6
## B3 13.1 0.446 87 12.2 14.0
##
## Kelompok = PMT Lokal:
## Waktu emmean SE df lower.CL upper.CL
## B0 12.3 0.440 87 11.4 13.2
## B1 12.7 0.442 87 11.8 13.6
## B2 13.1 0.445 87 12.2 13.9
## B3 13.4 0.446 87 12.6 14.3
##
## 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 87 26.397 <0.0001
##
## Kelompok = PMT Biskuit:
## model term df1 df2 F.ratio p.value
## Waktu 3 87 211.382 <0.0001
##
## Kelompok = PMT Lokal:
## model term df1 df2 F.ratio p.value
## Waktu 3 87 290.387 <0.0001
# Efek kelompok pada masing-masing waktu
joint_tests(aov2, by = "Waktu")
## Waktu = B0:
## model term df1 df2 F.ratio p.value
## Kelompok 2 87 1.000 0.3719
##
## Waktu = B1:
## model term df1 df2 F.ratio p.value
## Kelompok 2 87 0.532 0.5892
##
## Waktu = B2:
## model term df1 df2 F.ratio p.value
## Kelompok 2 87 0.297 0.7437
##
## Waktu = B3:
## model term df1 df2 F.ratio p.value
## Kelompok 2 87 0.170 0.8439
# 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
## B0 - B1 -0.117 0.0233 87 -5.014 <0.0001
## B0 - B2 -0.233 0.0330 87 -7.061 <0.0001
## B0 - B3 -0.340 0.0385 87 -8.835 <0.0001
## B1 - B2 -0.117 0.0236 87 -4.937 <0.0001
## B1 - B3 -0.223 0.0338 87 -6.602 <0.0001
## B2 - B3 -0.107 0.0228 87 -4.681 <0.0001
##
## Kelompok = PMT Biskuit:
## contrast estimate SE df t.ratio p.value
## B0 - B1 -0.317 0.0233 87 -13.610 <0.0001
## B0 - B2 -0.597 0.0330 87 -18.056 <0.0001
## B0 - B3 -0.947 0.0385 87 -24.599 <0.0001
## B1 - B2 -0.280 0.0236 87 -11.848 <0.0001
## B1 - B3 -0.630 0.0338 87 -18.625 <0.0001
## B2 - B3 -0.350 0.0228 87 -15.359 <0.0001
##
## Kelompok = PMT Lokal:
## contrast estimate SE df t.ratio p.value
## B0 - B1 -0.363 0.0233 87 -15.615 <0.0001
## B0 - B2 -0.743 0.0330 87 -22.495 <0.0001
## B0 - B3 -1.127 0.0385 87 -29.277 <0.0001
## B1 - B2 -0.380 0.0236 87 -16.080 <0.0001
## B1 - B3 -0.763 0.0338 87 -22.567 <0.0001
## B2 - B3 -0.383 0.0228 87 -16.822 <0.0001
##
## P value adjustment: holm method for 6 tests
# Tiap waktu vs baseline B0
contrast(
em2,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## Kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## B1 - B0 0.117 0.0233 87 5.014 <0.0001
## B2 - B0 0.233 0.0330 87 7.061 <0.0001
## B3 - B0 0.340 0.0385 87 8.835 <0.0001
##
## Kelompok = PMT Biskuit:
## contrast estimate SE df t.ratio p.value
## B1 - B0 0.317 0.0233 87 13.610 <0.0001
## B2 - B0 0.597 0.0330 87 18.056 <0.0001
## B3 - B0 0.947 0.0385 87 24.599 <0.0001
##
## Kelompok = PMT Lokal:
## contrast estimate SE df t.ratio p.value
## B1 - B0 0.363 0.0233 87 15.615 <0.0001
## B2 - B0 0.743 0.0330 87 22.495 <0.0001
## B3 - B0 1.127 0.0385 87 29.277 <0.0001
##
## P value adjustment: holm method for 3 tests
# 4G. POST-HOC ANTARKELOMPOK PADA SETIAP WAKTU
em2b <- emmeans(aov2, ~ Kelompok | Waktu)
em2b
## Waktu = B0:
## Kelompok emmean SE df lower.CL upper.CL
## Kontrol 13.0 0.440 87 12.1 13.8
## PMT Biskuit 12.1 0.440 87 11.3 13.0
## PMT Lokal 12.3 0.440 87 11.4 13.2
##
## Waktu = B1:
## Kelompok emmean SE df lower.CL upper.CL
## Kontrol 13.1 0.442 87 12.2 14.0
## PMT Biskuit 12.5 0.442 87 11.6 13.3
## PMT Lokal 12.7 0.442 87 11.8 13.6
##
## Waktu = B2:
## Kelompok emmean SE df lower.CL upper.CL
## Kontrol 13.2 0.445 87 12.3 14.1
## PMT Biskuit 12.7 0.445 87 11.8 13.6
## PMT Lokal 13.1 0.445 87 12.2 13.9
##
## Waktu = B3:
## Kelompok emmean SE df lower.CL upper.CL
## Kontrol 13.3 0.446 87 12.4 14.2
## PMT Biskuit 13.1 0.446 87 12.2 14.0
## PMT Lokal 13.4 0.446 87 12.6 14.3
##
## Confidence level used: 0.95
# Tukey
pairs(em2b, adjust = "tukey")
## Waktu = B0:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit 0.837 0.622 87 1.345 0.3741
## Kontrol - PMT Lokal 0.653 0.622 87 1.051 0.5472
## PMT Biskuit - PMT Lokal -0.183 0.622 87 -0.295 0.9532
##
## Waktu = B1:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit 0.637 0.625 87 1.019 0.5672
## Kontrol - PMT Lokal 0.407 0.625 87 0.651 0.7925
## PMT Biskuit - PMT Lokal -0.230 0.625 87 -0.368 0.9281
##
## Waktu = B2:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit 0.473 0.630 87 0.752 0.7335
## Kontrol - PMT Lokal 0.143 0.630 87 0.228 0.9719
## PMT Biskuit - PMT Lokal -0.330 0.630 87 -0.524 0.8598
##
## Waktu = B3:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit 0.230 0.630 87 0.365 0.9293
## Kontrol - PMT Lokal -0.133 0.630 87 -0.212 0.9756
## PMT Biskuit - PMT Lokal -0.363 0.630 87 -0.576 0.8330
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# Holm sebagai alternatif koreksi
pairs(em2b, adjust = "holm")
## Waktu = B0:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit 0.837 0.622 87 1.345 0.5459
## Kontrol - PMT Lokal 0.653 0.622 87 1.051 0.5927
## PMT Biskuit - PMT Lokal -0.183 0.622 87 -0.295 0.7688
##
## Waktu = B1:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit 0.637 0.625 87 1.019 0.9336
## Kontrol - PMT Lokal 0.407 0.625 87 0.651 1.0000
## PMT Biskuit - PMT Lokal -0.230 0.625 87 -0.368 1.0000
##
## Waktu = B2:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit 0.473 0.630 87 0.752 1.0000
## Kontrol - PMT Lokal 0.143 0.630 87 0.228 1.0000
## PMT Biskuit - PMT Lokal -0.330 0.630 87 -0.524 1.0000
##
## Waktu = B3:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit 0.230 0.630 87 0.365 1.0000
## Kontrol - PMT Lokal -0.133 0.630 87 -0.212 1.0000
## PMT Biskuit - PMT Lokal -0.363 0.630 87 -0.576 1.0000
##
## P value adjustment: holm method for 3 tests
# 4H. KONTRAS PERUBAHAN BULAN 0 -> BULAN 3
em_full <- emmeans(aov2, ~ Waktu * Kelompok)
em_full
## Waktu Kelompok emmean SE df lower.CL upper.CL
## B0 Kontrol 13.0 0.440 87 12.1 13.8
## B1 Kontrol 13.1 0.442 87 12.2 14.0
## B2 Kontrol 13.2 0.445 87 12.3 14.1
## B3 Kontrol 13.3 0.446 87 12.4 14.2
## B0 PMT Biskuit 12.1 0.440 87 11.3 13.0
## B1 PMT Biskuit 12.5 0.442 87 11.6 13.3
## B2 PMT Biskuit 12.7 0.445 87 11.8 13.6
## B3 PMT Biskuit 13.1 0.446 87 12.2 14.0
## B0 PMT Lokal 12.3 0.440 87 11.4 13.2
## B1 PMT Lokal 12.7 0.442 87 11.8 13.6
## B2 PMT Lokal 13.1 0.445 87 12.2 13.9
## B3 PMT Lokal 13.4 0.446 87 12.6 14.3
##
## Confidence level used: 0.95
# Pastikan urutan grid: B0-B3 untuk Kontrol, lalu PMT Biskuit, lalu PMT Lokal
stopifnot(
identical(as.character(em_full@grid$Waktu), rep(c("B0", "B1", "B2", "B3"), 3)),
identical(as.character(em_full@grid$Kelompok), rep(kel_lab, each = 4))
)
kc <- c(-1, 0, 0, 1) # kontras B3 - B0
z <- rep(0, 4)
# (i) Kenaikan BB (B3 - B0) di dalam tiap kelompok
contrast(
em_full,
list(
"Kontrol: B3-B0" = c(kc, z, z),
"PMT Biskuit: B3-B0" = c(z, kc, z),
"PMT Lokal: B3-B0" = c(z, z, kc)
),
adjust = "none"
)
## contrast estimate SE df t.ratio p.value
## Kontrol: B3-B0 0.340 0.0385 87 8.835 <0.0001
## PMT Biskuit: B3-B0 0.947 0.0385 87 24.599 <0.0001
## PMT Lokal: B3-B0 1.127 0.0385 87 29.277 <0.0001
# (ii) Perbedaan perubahan B3 - B0 antar kelompok
contrast(
em_full,
list(
"Biskuit - Kontrol" = c(-kc, kc, z),
"Lokal - Kontrol" = c(-kc, z, kc),
"Lokal - Biskuit" = c(z, -kc, kc)
),
adjust = "holm"
)
## contrast estimate SE df t.ratio p.value
## Biskuit - Kontrol 0.607 0.0544 87 11.147 <0.0001
## Lokal - Kontrol 0.787 0.0544 87 14.454 <0.0001
## Lokal - Biskuit 0.180 0.0544 87 3.307 0.0014
##
## P value adjustment: holm method for 3 tests
# 4I. TREN LINEAR ANTARKELOMPOK
tren_int <- summary(
contrast(
em_full,
interaction = c("poly", "pairwise"),
adjust = "none"
)
)
tren_int
## Waktu_poly Kelompok_pairwise estimate SE df t.ratio p.value
## linear Kontrol - PMT Biskuit -1.98333 0.1870 87 -10.628 <0.0001
## quadratic Kontrol - PMT Biskuit -0.04333 0.0501 87 -0.864 0.3899
## cubic Kontrol - PMT Biskuit -0.11667 0.0772 87 -1.511 0.1344
## linear Kontrol - PMT Lokal -2.62333 0.1870 87 -14.057 <0.0001
## quadratic Kontrol - PMT Lokal -0.03000 0.0501 87 -0.598 0.5512
## cubic Kontrol - PMT Lokal 0.00333 0.0772 87 0.043 0.9657
## linear PMT Biskuit - PMT Lokal -0.64000 0.1870 87 -3.429 0.0009
## quadratic PMT Biskuit - PMT Lokal 0.01333 0.0501 87 0.266 0.7910
## cubic PMT Biskuit - PMT Lokal 0.12000 0.0772 87 1.554 0.1238
# Nama kolom dapat berbeda menurut versi emmeans.
# Cek names(tren_int) terlebih dahulu.
names(tren_int)
## [1] "Waktu_poly" "Kelompok_pairwise" "estimate"
## [4] "SE" "df" "t.ratio"
## [7] "p.value"
# Biasanya komponen tren Waktu dapat dipisahkan seperti berikut:
if ("Waktu_poly" %in% names(tren_int)) {
tren_lin <- subset(tren_int, Waktu_poly == "linear")
} else if ("Waktu.poly" %in% names(tren_int)) {
tren_lin <- subset(tren_int, Waktu.poly == "linear")
} else {
tren_lin <- tren_int
}
if ("p.value" %in% names(tren_lin)) {
tren_lin$p_holm <- p.adjust(tren_lin$p.value, method = "holm")
}
tren_lin
## Waktu_poly Kelompok_pairwise estimate SE df t.ratio
## 1 linear Kontrol - PMT Biskuit -1.983333 0.1866208 87 -10.627610
## 4 linear Kontrol - PMT Lokal -2.623333 0.1866208 87 -14.057024
## 7 linear PMT Biskuit - PMT Lokal -0.640000 0.1866208 87 -3.429414
## p.value p_holm
## 1 2.141143e-17 4.282286e-17
## 4 4.148663e-24 1.244599e-23
## 7 9.267448e-04 9.267448e-04
# 4J. PEMBANDING: SELISIH BB (B3 - B0) DAN ANCOVA
data_wide <- data_wide %>%
mutate(Selisih = BB_B3 - BB_B0)
data_wide %>%
group_by(Kelompok) %>%
get_summary_stats(Selisih, type = "mean_sd")
## # A tibble: 3 × 5
## Kelompok variable n mean sd
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol Selisih 30 0.34 0.213
## 2 PMT Biskuit Selisih 30 0.947 0.229
## 3 PMT Lokal Selisih 30 1.13 0.189
# Asumsi
data_wide %>% group_by(Kelompok) %>% shapiro_test(Selisih)
## # A tibble: 3 × 4
## Kelompok variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Kontrol Selisih 0.942 0.104
## 2 PMT Biskuit Selisih 0.947 0.137
## 3 PMT Lokal Selisih 0.951 0.179
levene_test(data_wide, Selisih ~ Kelompok)
## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 87 0.0852 0.918
# ANOVA satu arah + Tukey
aov_sel <- aov(Selisih ~ Kelompok, data = data_wide)
summary(aov_sel)
## Df Sum Sq Mean Sq F value Pr(>F)
## Kelompok 2 10.193 5.096 114.7 <2e-16 ***
## Residuals 87 3.865 0.044
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov_sel)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Selisih ~ Kelompok, data = data_wide)
##
## $Kelompok
## diff lwr upr p adj
## PMT Biskuit-Kontrol 0.6066667 0.47689442 0.7364389 0.0000000
## PMT Lokal-Kontrol 0.7866667 0.65689442 0.9164389 0.0000000
## PMT Lokal-PMT Biskuit 0.1800000 0.05022775 0.3097722 0.0038808
# Jika asumsi tidak terpenuhi:
# kruskal.test(Selisih ~ Kelompok, data = data_wide)
# Paired t-test B0 vs B3 pada tiap kelompok
data_wide %>%
pivot_longer(c(BB_B0, BB_B3), names_to = "waktu_uji", values_to = "bb_uji") %>%
group_by(Kelompok) %>%
t_test(bb_uji ~ waktu_uji, paired = TRUE)
## # A tibble: 3 × 9
## Kelompok .y. group1 group2 n1 n2 statistic df p
## * <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 Kontrol bb_uji BB_B0 BB_B3 30 30 -8.76 29 1.23e- 9
## 2 PMT Biskuit bb_uji BB_B0 BB_B3 30 30 -22.7 29 5.24e-20
## 3 PMT Lokal bb_uji BB_B0 BB_B3 30 30 -32.6 29 2.10e-24
# ANCOVA: BB bulan 3 antarkelompok, dikontrol BB awal, usia, dan jenis kelamin
ancova <- lm(BB_B3 ~ BB_B0 + Usia_bulan + Jenis_Kelamin + Kelompok,
data = data_wide)
car::Anova(ancova, type = 3)
## Anova Table (Type III tests)
##
## Response: BB_B3
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.289 1 6.6414 0.01171 *
## BB_B0 78.306 1 1796.5835 < 2e-16 ***
## Usia_bulan 0.152 1 3.4879 0.06531 .
## Jenis_Kelamin 0.000 1 0.0032 0.95512
## Kelompok 8.871 2 101.7674 < 2e-16 ***
## Residuals 3.661 84
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(ancova, pairwise ~ Kelompok, adjust = "tukey")
## $emmeans
## Kelompok emmean SE df lower.CL upper.CL
## Kontrol 12.8 0.0393 84 12.7 12.9
## PMT Biskuit 13.4 0.0389 84 13.3 13.5
## PMT Lokal 13.6 0.0387 84 13.5 13.7
##
## Results are averaged over the levels of: Jenis_Kelamin
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit -0.592 0.0565 84 -10.471 <0.0001
## Kontrol - PMT Lokal -0.774 0.0560 84 -13.811 <0.0001
## PMT Biskuit - PMT Lokal -0.182 0.0540 84 -3.364 0.0033
##
## Results are averaged over the levels of: Jenis_Kelamin
## P value adjustment: tukey method for comparing a family of 3 estimates
# 5. LINEAR MIXED MODEL (LMM)
# LMM random intercept
lmm1 <- lmer(
BB_kg ~ Kelompok * Waktu + (1 | ID),
data = data_long,
REML = TRUE
)
summary(lmm1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BB_kg ~ Kelompok * Waktu + (1 | ID)
## Data: data_long
##
## REML criterion at convergence: 193.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8798 -0.5037 0.0005 0.5458 3.5407
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 5.87880 2.4246
## Residual 0.01334 0.1155
## Number of obs: 360, groups: ID, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.875278 0.255650 86.999999 50.363 < 2e-16 ***
## Kelompok1 0.270556 0.361544 87.000014 0.748 0.456
## Kelompok2 -0.273611 0.361544 87.000014 -0.757 0.451
## Waktu1 -0.398611 0.010544 261.000000 -37.805 < 2e-16 ***
## Waktu2 -0.133056 0.010544 261.000000 -12.619 < 2e-16 ***
## Waktu3 0.125833 0.010544 261.000000 11.934 < 2e-16 ***
## Kelompok1:Waktu1 0.226111 0.014911 261.000000 15.164 < 2e-16 ***
## Kelompok2:Waktu1 -0.066389 0.014911 261.000000 -4.452 1.26e-05 ***
## Kelompok1:Waktu2 0.077222 0.014911 261.000000 5.179 4.47e-07 ***
## Kelompok2:Waktu2 -0.015278 0.014911 261.000000 -1.025 0.307
## Kelompok1:Waktu3 -0.065000 0.014911 261.000000 -4.359 1.88e-05 ***
## Kelompok2:Waktu3 0.005833 0.014911 261.000000 0.391 0.696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Klmpk1 Klmpk2 Waktu1 Waktu2 Waktu3 Kl1:W1 Kl2:W1 Kl1:W2
## Kelompok1 0.000
## Kelompok2 0.000 -0.500
## Waktu1 0.000 0.000 0.000
## Waktu2 0.000 0.000 0.000 -0.333
## Waktu3 0.000 0.000 0.000 -0.333 -0.333
## Klmpk1:Wkt1 0.000 0.000 0.000 0.000 0.000 0.000
## Klmpk2:Wkt1 0.000 0.000 0.000 0.000 0.000 0.000 -0.500
## Klmpk1:Wkt2 0.000 0.000 0.000 0.000 0.000 0.000 -0.333 0.167
## Klmpk2:Wkt2 0.000 0.000 0.000 0.000 0.000 0.000 0.167 -0.333 -0.500
## Klmpk1:Wkt3 0.000 0.000 0.000 0.000 0.000 0.000 -0.333 0.167 -0.333
## Klmpk2:Wkt3 0.000 0.000 0.000 0.000 0.000 0.000 0.167 -0.333 0.167
## Kl2:W2 Kl1:W3
## Kelompok1
## Kelompok2
## Waktu1
## Waktu2
## Waktu3
## Klmpk1:Wkt1
## Klmpk2:Wkt1
## Klmpk1:Wkt2
## Klmpk2:Wkt2
## Klmpk1:Wkt3 0.167
## Klmpk2:Wkt3 -0.333 -0.500
# LMM random intercept + random slope
lmm2 <- lmer(
BB_kg ~ Kelompok * Waktu + (1 + Waktu_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: BB_kg ~ Kelompok * Waktu + (1 + Waktu_num | ID)
## Data: data_long
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 136.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.34694 -0.45022 -0.01438 0.48019 2.26498
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 5.803933 2.40914
## Waktu_num 0.003834 0.06192 0.15
## Residual 0.006951 0.08337
## Number of obs: 360, groups: ID, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.875278 0.255652 86.998052 50.362 < 2e-16 ***
## Kelompok1 0.270556 0.361547 86.998151 0.748 0.456281
## Kelompok2 -0.273611 0.361547 86.998151 -0.757 0.451227
## Waktu1 -0.398611 0.012400 118.737351 -32.145 < 2e-16 ***
## Waktu2 -0.133056 0.008281 244.689332 -16.068 < 2e-16 ***
## Waktu3 0.125833 0.008281 244.689332 15.196 < 2e-16 ***
## Kelompok1:Waktu1 0.226111 0.017537 118.737353 12.893 < 2e-16 ***
## Kelompok2:Waktu1 -0.066389 0.017537 118.737353 -3.786 0.000242 ***
## Kelompok1:Waktu2 0.077222 0.011711 244.689332 6.594 2.62e-10 ***
## Kelompok2:Waktu2 -0.015278 0.011711 244.689332 -1.305 0.193262
## Kelompok1:Waktu3 -0.065000 0.011711 244.689332 -5.550 7.40e-08 ***
## Kelompok2:Waktu3 0.005833 0.011711 244.689332 0.498 0.618853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Klmpk1 Klmpk2 Waktu1 Waktu2 Waktu3 Kl1:W1 Kl2:W1 Kl1:W2
## Kelompok1 0.000
## Kelompok2 0.000 -0.500
## Waktu1 -0.149 0.000 0.000
## Waktu2 -0.075 0.000 0.000 0.123
## Waktu3 0.075 0.000 0.000 -0.499 -0.437
## Klmpk1:Wkt1 0.000 -0.149 0.075 0.000 0.000 0.000
## Klmpk2:Wkt1 0.000 0.075 -0.149 0.000 0.000 0.000 -0.500
## Klmpk1:Wkt2 0.000 -0.075 0.037 0.000 0.000 0.000 0.123 -0.062
## Klmpk2:Wkt2 0.000 0.037 -0.075 0.000 0.000 0.000 -0.062 0.123 -0.500
## Klmpk1:Wkt3 0.000 0.075 -0.037 0.000 0.000 0.000 -0.499 0.250 -0.437
## Klmpk2:Wkt3 0.000 -0.037 0.075 0.000 0.000 0.000 0.250 -0.499 0.218
## Kl2:W2 Kl1:W3
## Kelompok1
## Kelompok2
## Waktu1
## Waktu2
## Waktu3
## Klmpk1:Wkt1
## Klmpk2:Wkt1
## Klmpk1:Wkt2
## Klmpk2:Wkt2
## Klmpk1:Wkt3 0.218
## Klmpk2:Wkt3 -0.437 -0.500
# Bandingkan struktur random effect
anova(lmm1, lmm2, refit = FALSE)
## Data: data_long
## Models:
## lmm1: BB_kg ~ Kelompok * Waktu + (1 | ID)
## lmm2: BB_kg ~ Kelompok * Waktu + (1 + Waktu_num | ID)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 221.55 275.95 -96.773 193.55
## lmm2 16 168.25 230.42 -68.123 136.25 57.3 2 3.61e-13 ***
## ---
## 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.0052 0.00262 2 87.00 0.3776 0.6866
## Waktu 8.5589 2.85298 3 185.37 408.5763 <2e-16 ***
## Kelompok:Waktu 1.5149 0.25248 6 206.40 36.1183 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.998
## Unadjusted ICC: 0.972
# LMM dengan kovariat usia dan jenis kelamin (Waktu_num sebagai numerik)
lmm3 <- lmer(
BB_kg ~ Kelompok * Waktu_num + Usia_bulan + Jenis_Kelamin + (1 | ID),
data = data_long,
REML = TRUE
)
anova(lmm3, 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.019 0.009 2 86.084 0.7222 0.4886
## Waktu_num 32.133 32.133 1 267.000 2443.3082 <2e-16 ***
## Usia_bulan 6.462 6.462 1 85.000 491.3483 <2e-16 ***
## Jenis_Kelamin 0.024 0.024 1 85.000 1.8062 0.1825
## Kelompok:Waktu_num 5.613 2.806 2 267.000 213.3784 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmm3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BB_kg ~ Kelompok * Waktu_num + Usia_bulan + Jenis_Kelamin + (1 |
## ID)
## Data: data_long
##
## REML criterion at convergence: 1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8760 -0.5403 -0.0130 0.5250 3.6046
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.88462 0.9405
## Residual 0.01315 0.1147
## Number of obs: 360, groups: ID, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.429612 0.290550 85.132562 22.129 < 2e-16 ***
## Kelompok1 -0.175643 0.147694 86.033282 -1.189 0.238
## Kelompok2 0.110239 0.143289 86.098344 0.769 0.444
## Waktu_num 0.267222 0.005406 267.000000 49.430 < 2e-16 ***
## Usia_bulan 0.161073 0.007267 85.000002 22.166 < 2e-16 ***
## Jenis_Kelamin1 0.137981 0.102668 85.000000 1.344 0.183
## Kelompok1:Waktu_num -0.153556 0.007645 267.000000 -20.085 < 2e-16 ***
## Kelompok2:Waktu_num 0.044778 0.007645 267.000000 5.857 1.38e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Klmpk1 Klmpk2 Wkt_nm Us_bln Jns_K1 Kl1:W_
## Kelompok1 0.186
## Kelompok2 -0.129 -0.523
## Waktu_num -0.028 0.000 0.000
## Usia_bulan -0.939 -0.207 0.142 0.000
## Jenis_Klmn1 -0.095 0.203 -0.103 0.000 0.059
## Klmpk1:Wkt_ 0.000 -0.078 0.040 0.000 0.000 0.000
## Klmpk2:Wkt_ 0.000 0.039 -0.080 0.000 0.000 0.000 -0.500
# Laju kenaikan BB per bulan (kg/bulan) tiap kelompok dan perbandingannya
emtrends(lmm3, pairwise ~ Kelompok, var = "Waktu_num", adjust = "tukey")
## $emtrends
## Kelompok Waktu_num.trend SE df lower.CL upper.CL
## Kontrol 0.114 0.00936 267 0.0952 0.132
## PMT Biskuit 0.312 0.00936 267 0.2936 0.330
## PMT Lokal 0.376 0.00936 267 0.3576 0.394
##
## Results are averaged over the levels of: Jenis_Kelamin
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT Biskuit -0.198 0.0132 267 -14.977 <0.0001
## Kontrol - PMT Lokal -0.262 0.0132 267 -19.810 <0.0001
## PMT Biskuit - PMT Lokal -0.064 0.0132 267 -4.833 <0.0001
##
## Results are averaged over the levels of: Jenis_Kelamin
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
# 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.9941, p-value = 0.177
# 6. OPSIONAL: SIMULASI MISSING VALUE UNTUK DEMONSTRASI LMM
# Bagian ini TIDAK mengubah data utama.
# Hanya digunakan bila dosen meminta demonstrasi keunggulan LMM.
set.seed(2026)
data_miss <- data_long
idx_miss <- sample(
which(data_miss$Waktu != "B0"),
30,
replace = FALSE
)
data_miss$BB_kg[idx_miss] <- NA
lmm_miss <- lmer(
BB_kg ~ Kelompok * Waktu + (1 + Waktu_num | ID),
data = data_miss,
REML = TRUE,
na.action = na.exclude,
control = lmerControl(optimizer = "bobyqa")
)
anova(lmm_miss, ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Kelompok 0.0048 0.00241 2 87.00 0.3678 0.6933
## Waktu 7.2185 2.40618 3 167.59 365.1244 <2e-16 ***
## Kelompok:Waktu 1.2838 0.21397 6 185.16 32.4319 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Jumlah balita yang kehilangan >= 1 nilai (akan dibuang oleh RM ANOVA)
n_distinct(
data_miss$ID[is.na(data_miss$BB_kg)]
)
## [1] 28
# 7. EXPORT HASIL
# Data format LONG dan WIDE (CSV)
write.csv(
data_long,
"dataset_PMT_balita_LONG.csv",
row.names = FALSE
)
write.csv(
data_wide,
"dataset_PMT_balita_WIDE.csv",
row.names = FALSE
)
# Statistik deskriptif
write.csv(
deskriptif,
"hasil_01_statistik_deskriptif_BB.csv",
row.names = FALSE
)
# Normalitas
write.csv(
shapiro_mixed,
"hasil_02_shapiro_Kelompok_x_Waktu.csv",
row.names = FALSE
)
# Levene
write.csv(
levene_mixed,
"hasil_03_levene_per_Waktu.csv",
row.names = FALSE
)
# ANOVA utama
hasil_anova <- as.data.frame(
nice(
aov2,
correction = "GG",
es = c("ges", "pes")
)
)
write.csv(
hasil_anova,
"hasil_04_mixed_ANOVA_GG.csv",
row.names = FALSE
)
# Post-hoc antarkelompok
posthoc_group <- as.data.frame(
pairs(em2b, adjust = "holm")
)
write.csv(
posthoc_group,
"hasil_05_posthoc_antar_kelompok_Holm.csv",
row.names = FALSE
)
# Post-hoc waktu dalam kelompok
posthoc_time <- as.data.frame(
pairs(em2, adjust = "holm")
)
write.csv(
posthoc_time,
"hasil_06_posthoc_dalam_kelompok_Holm.csv",
row.names = FALSE
)
# EMM
emm_group_time <- as.data.frame(em2b)
write.csv(
emm_group_time,
"hasil_07_estimated_marginal_means.csv",
row.names = FALSE
)
# Gambar
ggsave(
"grafik_01_profile_BB.png",
p_profil,
width = 8,
height = 5.5,
dpi = 300
)
ggsave(
"grafik_02_spaghetti_BB.png",
p_spag,
width = 9,
height = 6,
dpi = 300
)
ggsave(
"grafik_03_interaksi_Kelompok_x_Waktu.png",
p_interaksi,
width = 8,
height = 5.5,
dpi = 300
)
# 8. RINGKASAN AKHIR DI CONSOLE
cat("\n============================================================\n")
##
## ============================================================
cat("ANALISIS REPEATED MEASURES BERAT BADAN BALITA SELESAI\n")
## ANALISIS REPEATED MEASURES BERAT BADAN BALITA SELESAI
cat("============================================================\n")
## ============================================================
cat("Jumlah balita :", n_distinct(data_long$ID), "\n")
## Jumlah balita : 90
cat("Kelompok :", n_distinct(data_long$Kelompok), "\n")
## Kelompok : 3
cat("Waktu penimbangan:", n_distinct(data_long$Waktu), "\n")
## Waktu penimbangan: 4
cat("Outcome : BB_kg\n")
## Outcome : BB_kg
cat("\nSemua file hasil tersimpan di:\n")
##
## Semua file hasil tersimpan di:
cat(getwd(), "\n")
## C:/Users/ASUS/Downloads/ndari R studio
cat("============================================================\n")
## ============================================================
# =============================================================================
# SELESAI
# =============================================================================