# ANALISIS PENGUKURAN BERULANG KADAR BERAT BADAN: 2 KELOMPOK x 3 WAKTU
# Mahasiswa: Ika Wulan Sari
# NIM: 2611018005
# Jalankan seluruh script melalui Source di RStudio.
# Letakkan script ini dan journal.pone.0257326.s004 (1).xlsx dalam folder yang sama.
# Laporan HTML dibuat otomatis setelah seluruh analisis selesai.
#
# **Nama:** Ika Wulan Sari
# **NIM:** 2611018005
# **Mata Kuliah:** Biostatistika Intermediate
# **Program Studi:** Magister Kesehatan Masyarakat Universitas Mulawarman
#
# ## Metode
# Data mencakup dua kelompok dan tiga pengukuran berat badan per peserta.
# Analisis meliputi repeated-measures ANOVA satu arah pada kelompok intervention, mixed ANOVA kelompok x waktu, koreksi
# Greenhouse-Geisser, pendekatan multivariat, post hoc, kontras interaksi,
# alternatif nonparametrik, dan linear mixed model (LMM).
# Baseline, bulan ke-3, dan bulan ke-6 berjarak sama.
# Outcome utama: berat badan dari WeightKg, Weightkg.M3, dan Weightkg.M6.
# ANOVA/nonparametrik: peserta lengkap; LMM: seluruh pengukuran tersedia.
# Sumber: journal.pone.0257326.s004 (1).xlsx, sheet data_set.
# Satuan berat badan: kg berdasarkan nama kolom. Kode -999 dan 999 pada berat badan/usia menjadi NA.
# ID dibuat dari urutan baris karena tidak ada kolom ID pada workbook.
# ======================= PARAMETER ANALISIS =========================
FILE_DATA <- "journal.pone.0257326.s004 (1).xlsx"
SHEET_DATA <- "data_set"
FOLDER_HASIL <- "Hasil_Analisis_Berat_Badan_Ika_Wulan_Sari"
INSTALL_PAKET_OTOMATIS <- TRUE
BUKA_HTML_OTOMATIS <- TRUE
analisis_berat_badan <- function() {
paket <- c("readxl", "dplyr", "tidyr", "ggplot2", "afex", "emmeans", "rstatix",
"car", "effectsize", "lme4", "lmerTest", "pbkrtest",
"performance", "ggpubr", "WRS2", "htmltools", "base64enc")
belum <- paket[!vapply(paket, requireNamespace, logical(1), quietly = TRUE)]
if (length(belum)) {
if (!INSTALL_PAKET_OTOMATIS)
stop("Paket belum tersedia: ", paste(belum, collapse = ", "))
install.packages(belum, repos = "https://cloud.r-project.org")
}
gagal <- paket[!vapply(paket, requireNamespace, logical(1), quietly = TRUE)]
if (length(gagal)) stop("Paket belum tersedia: ", paste(gagal, collapse = ", "))
suppressPackageStartupMessages({
library(dplyr); library(tidyr); library(ggplot2)
})
opsi_lama <- options(contrasts = c("contr.sum", "contr.poly"), width = 120)
on.exit(options(opsi_lama), add = TRUE)
afex::afex_options(emmeans_model = "multivariate")
theme_set(theme_bw(base_size = 12))
set.seed(2026)
# Cari data di folder script dan direktori kerja.
sumber <- unlist(lapply(sys.frames(), function(fr) {
if (!is.null(fr$ofile)) as.character(fr$ofile) else character(0)
}))
if (requireNamespace("rstudioapi", quietly = TRUE) && rstudioapi::isAvailable()) {
aktif <- tryCatch(rstudioapi::getSourceEditorContext()$path,
error = function(e) "")
if (nzchar(aktif)) sumber <- c(sumber, aktif)
}
nama_data <- c(FILE_DATA, "01-journal.pone.0257326.s004-1-.xlsx")
kandidat <- unique(c(unlist(lapply(dirname(sumber), function(d) file.path(d, nama_data))), nama_data))
kandidat <- kandidat[file.exists(kandidat)]
if (!length(kandidat)) stop(
"Berkas Excel tidak ditemukan. Letakkan data dan script dalam folder yang sama, ",
"atau isi FILE_DATA dengan lokasi lengkap berkas. Nama yang dicari: ", FILE_DATA)
path_data <- normalizePath(kandidat[1], winslash = "/", mustWork = TRUE)
if (!SHEET_DATA %in% readxl::excel_sheets(path_data))
stop("Sheet tidak ditemukan: ", SHEET_DATA)
folder <- file.path(dirname(path_data), FOLDER_HASIL)
message("Data: ", path_data, "\nSheet: ", SHEET_DATA, "\nHasil: ", folder)
dir.create(folder, recursive = TRUE, showWarnings = FALSE)
hasil <- list(); isi <- list(); grafik <- list()
teks <- function(judul, x) {
cat("\n", judul, "\n", sep = "")
cat(x, "\n")
isi[[length(isi) + 1L]] <<- htmltools::tagList(
htmltools::tags$h2(judul), htmltools::tags$p(x))
}
tampil <- function(nama, x) {
hasil[[nama]] <<- x
cat("\n", gsub("_", " ", nama), "\n", sep = "")
print(x)
cetak <- capture.output(print(x))
writeLines(cetak, file.path(folder, paste0(nama, "_Ika_Wulan_Sari.txt")))
if (is.data.frame(x)) write.csv(x, file.path(folder, paste0(nama, "_Ika_Wulan_Sari.csv")), row.names = FALSE)
isi[[length(isi) + 1L]] <<- htmltools::tagList(
htmltools::tags$h2(gsub("_", " ", nama)),
htmltools::tags$pre(paste(cetak, collapse = "\n")))
invisible(x)
}
gambar <- function(nama, p, lebar = 11, tinggi = 6) {
print(p)
berkas <- file.path(folder, paste0(nama, "_Ika_Wulan_Sari.png"))
ggsave(berkas, plot = p, width = lebar, height = tinggi, dpi = 300, bg = "white")
grafik[[nama]] <<- p
isi[[length(isi) + 1L]] <<- htmltools::tagList(
htmltools::tags$h2(gsub("_", " ", nama)),
htmltools::tags$img(src = base64enc::dataURI(file = berkas, mime = "image/png"),
style = "max-width:100%;height:auto;"))
invisible(p)
}
opsional <- function(nama, ekspresi) {
tryCatch(tampil(nama, ekspresi), error = function(e) {
teks(nama, paste("Analisis tidak dapat diestimasi:", conditionMessage(e)))
invisible(NULL)
})
}
fp <- function(p) ifelse(is.na(p), "NA", ifelse(p < .001, "< 0,001",
paste0("= ", formatC(p, digits = 3, format = "f", decimal.mark = ","))))
asli <- as.data.frame(readxl::read_excel(path_data, sheet = SHEET_DATA))
names(asli) <- trimws(names(asli))
wajib <- c("Group", "Ageyears", "Gender", "WeightKg", "Weightkg.M3", "Weightkg.M6")
if (!all(wajib %in% names(asli))) stop("Kolom tidak ditemukan: ", paste(setdiff(wajib, names(asli)), collapse = ", "))
angka <- function(x) {
x <- trimws(as.character(x))
x[x %in% c("", "-999", "999", "#NULL!", "NA")] <- NA_character_
y <- suppressWarnings(as.numeric(x))
if (any(!is.na(x) & is.na(y))) stop("Nilai tidak numerik pada kolom analisis.")
if (any(!is.na(y) & !is.finite(y))) stop("Nilai tak hingga pada kolom analisis.")
y
}
if (!nrow(asli)) stop("Tidak ada peserta.")
dat_wide <- data.frame(id = sprintf("ROW%03d", seq_len(nrow(asli))),
baris_excel = seq_len(nrow(asli)) + 1L,
kelompok = trimws(as.character(asli$Group)), usia = angka(asli$Ageyears),
jk = trimws(as.character(asli$Gender)), BB_M0 = angka(asli$WeightKg),
BB_M3 = angka(asli$Weightkg.M3), BB_M6 = angka(asli$Weightkg.M6))
dat_wide$jk[dat_wide$jk %in% c("9", "999", "-999", "")] <- NA_character_
kelompok <- c("control", "intervention")
if (anyNA(dat_wide$kelompok) || !setequal(unique(dat_wide$kelompok), kelompok))
stop("Group harus berisi control dan intervention.")
kolom_bb <- c("BB_M0", "BB_M3", "BB_M6")
if (any(as.matrix(dat_wide[kolom_bb]) <= 0, na.rm = TRUE)) stop("Berat badan nonpositif ditemukan.")
bulan <- c(0, 3, 6)
dat_wide <- dat_wide |> mutate(id = factor(id), kelompok = factor(kelompok, levels = c("control", "intervention")))
long_semua <- dat_wide |>
pivot_longer(all_of(kolom_bb), names_to = "waktu", values_to = "berat_badan") |>
mutate(waktu = factor(waktu, levels = kolom_bb, labels = c("M0", "M3", "M6")),
bulan = bulan[as.integer(waktu)]) |> arrange(id, waktu)
lengkap <- complete.cases(dat_wide[kolom_bb])
wide_lengkap <- droplevels(dat_wide[lengkap, , drop = FALSE])
dat_long <- droplevels(filter(long_semua, id %in% wide_lengkap$id))
dat_lmm <- droplevels(filter(long_semua, !is.na(berat_badan)))
if (any(table(wide_lengkap$kelompok) < 3)) stop("Peserta lengkap per kelompok kurang dari tiga.")
teks("Identitas", paste("Nama: Ika Wulan Sari | NIM: 2611018005 |",
"Mata Kuliah: Biostatistika Intermediate | Magister Kesehatan Masyarakat Universitas Mulawarman"))
teks("Desain", paste(nrow(dat_wide), "peserta; dua kelompok; tiga waktu (M0, M3, M6).",
sum(lengkap), "peserta lengkap untuk ANOVA;", sum(is.na(long_semua$berat_badan)),
"pengukuran hilang; LMM menggunakan", nrow(dat_lmm), "pengukuran tersedia. Outcome: berat badan (kg)."))
tampil("00_Kelengkapan", long_semua |> group_by(kelompok, waktu) |>
summarise(peserta = n(), tersedia = sum(!is.na(berat_badan)), hilang = sum(is.na(berat_badan)), .groups = "drop"))
teks("Penanganan data", "ID ROW dibuat dari urutan baris, dengan baris_excel untuk penelusuran; setiap baris diasumsikan satu peserta sesuai struktur wide. Kode -999 dan 999 pada berat badan/usia dibaca sebagai NA. Nilai berat badan 999 pada bulan 6 diperlakukan sebagai kode hilang berdasarkan pola kode pada baris yang sama (usia 999, gender 9, aktivitas -999). Tidak dilakukan imputasi. Kolom turunan overall_PA dan delta_M6_BL tidak digunakan. ANOVA, Friedman, Wilcoxon, dan robust RM memakai peserta lengkap; LMM memakai seluruh pengukuran berat badan tersedia. Hasil kedua pendekatan dapat berbeda karena sampel berbeda.")
tampil("01_Jumlah_Peserta", dat_wide |> count(kelompok, name = "n"))
tampil("01a_Usia", dat_wide |> group_by(kelompok) |> summarise(n = n(), mean = mean(usia, na.rm = TRUE), sd = sd(usia, na.rm = TRUE), min = min(usia, na.rm = TRUE), max = max(usia, na.rm = TRUE), .groups = "drop"))
tampil("01b_Jenis_Kelamin", dat_wide |> count(kelompok, jk))
teks("Berkas data yang digunakan", paste(path_data, "| Sheet:", SHEET_DATA))
teks("Sumber data", "Workbook lampiran journal.pone.0257326.s004 (1).xlsx, sheet data_set. Berat badan menggunakan WeightKg, Weightkg.M3, dan Weightkg.M6; kelompok memakai Group. Analisis ini adalah adaptasi metode latihan, bukan reproduksi seluruh analisis artikel.")
tampil("02_Data_Wide", head(dat_wide))
tampil("03_Data_Long", head(dat_long, 12))
desk <- dat_lmm |> group_by(kelompok, waktu, bulan) |>
summarise(n = n(), mean = mean(berat_badan), sd = sd(berat_badan), median = median(berat_badan),
min = min(berat_badan), max = max(berat_badan), se = sd / sqrt(n),
lower = mean - qt(.975, n - 1) * se,
upper = mean + qt(.975, n - 1) * se, .groups = "drop")
tampil("04_Deskriptif", desk)
teks("Sampel deskriptif", "Tabel deskriptif dan grafik profil rerata menggunakan seluruh berat badan yang tersedia pada tiap waktu. Tabel ANOVA serta grafik profil mixed ANOVA menggunakan peserta dengan tiga pengukuran lengkap.")
tampil("05_Kovarians", cov(wide_lengkap[kolom_bb]))
tampil("06_Korelasi", cor(wide_lengkap[kolom_bb]))
pasangan <- combn(kolom_bb, 2)
tampil("07_Varians_Selisih", data.frame(
pasangan = apply(pasangan, 2, paste, collapse = " - "),
varians = apply(pasangan, 2, function(p) var(wide_lengkap[[p[1]]] - wide_lengkap[[p[2]]]))))
gambar("G01_Profil_Rerata_CI95", ggplot(desk, aes(bulan, mean, color = kelompok, group = kelompok)) +
geom_line(linewidth = 1) + geom_point(size = 3) +
geom_errorbar(aes(ymin = lower, ymax = upper), width = .5) +
scale_x_continuous(breaks = bulan) +
labs(title = "Profil rerata berat badan dan interval kepercayaan 95%", x = "Bulan", y = "Berat badan (kg)", color = "Kelompok") +
theme(legend.position = "bottom"))
gambar("G02_Lintasan_Individu", ggplot(long_semua, aes(bulan, berat_badan, group = id)) +
geom_line(alpha = .25) + stat_summary(aes(group = kelompok), fun = mean,
geom = "line", color = "firebrick", linewidth = 1.2) +
facet_wrap(~ kelompok) + scale_x_continuous(breaks = bulan) +
labs(title = "Lintasan individu dan rerata kelompok", x = "Bulan", y = "Berat badan (kg)"))
gambar("G03_Boxplot", ggplot(dat_long, aes(waktu, berat_badan, fill = kelompok)) +
geom_boxplot() + labs(title = "Distribusi berat badan berdasarkan kelompok dan waktu", x = "Waktu", y = "Berat badan (kg)", fill = "Kelompok") +
theme(legend.position = "bottom"))
d1 <- droplevels(filter(dat_long, kelompok == "intervention"))
tampil("08_Outlier_Satu_Kelompok", d1 |> group_by(waktu) |> rstatix::identify_outliers(berat_badan))
tampil("09_Shapiro_Satu_Kelompok", d1 |> group_by(waktu) |> rstatix::shapiro_test(berat_badan))
gambar("G04_QQ_Satu_Kelompok", ggpubr::ggqqplot(d1, "berat_badan", facet.by = "waktu"))
aov1_rs <- rstatix::anova_test(data = d1, dv = berat_badan, wid = id, within = waktu, effect.size = "pes")
tampil("10_ANOVA_Rstatix_Satu_Kelompok", aov1_rs)
aov1 <- afex::aov_ez(id = "id", dv = "berat_badan", data = d1, within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
tampil("11_RM_ANOVA", aov1)
tampil("12_Sferisitas_dan_Koreksi_RM", summary(aov1))
tampil("13_Multivariat_RM", aov1$Anova)
tampil("14_Eta_Kuadrat_RM", effectsize::eta_squared(aov1, partial = TRUE))
em1 <- emmeans::emmeans(aov1, ~ waktu)
tampil("15_EMM_RM", as.data.frame(em1))
tampil("16_Posthoc_RM_Bonferroni", as.data.frame(pairs(em1, adjust = "bonferroni")))
tampil("17_RM_vs_Baseline_Holm", as.data.frame(emmeans::contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm")))
tampil("18_Tren_Polinomial_RM", as.data.frame(emmeans::contrast(em1, "poly")))
tampil("19_Friedman", rstatix::friedman_test(d1, berat_badan ~ waktu | id))
tampil("20_Kendall_W", rstatix::friedman_effsize(d1, berat_badan ~ waktu | id))
tampil("21_Wilcoxon_Berpasangan", d1 |> arrange(waktu, id) |>
rstatix::wilcox_test(berat_badan ~ waktu, paired = TRUE, p.adjust.method = "bonferroni"))
opsional("22_Robust_RM_Trim20", WRS2::rmanova(d1$berat_badan, d1$waktu, d1$id, tr = .2))
tampil("23_Outlier_Mixed", dat_long |> group_by(kelompok, waktu) |> rstatix::identify_outliers(berat_badan))
tampil("24_Shapiro_Mixed", dat_long |> group_by(kelompok, waktu) |> rstatix::shapiro_test(berat_badan))
gambar("G05_QQ_Mixed", ggpubr::ggqqplot(dat_long, "berat_badan", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok), tinggi = 9)
tampil("25_Levene", dat_long |> group_by(waktu) |> rstatix::levene_test(berat_badan ~ kelompok))
tampil("26_Box_M", rstatix::box_m(wide_lengkap[kolom_bb], wide_lengkap$kelompok))
aov2 <- afex::aov_ez(id = "id", dv = "berat_badan", data = dat_long, between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
tampil("27_Mixed_ANOVA", aov2)
tampil("28_Sferisitas_dan_Koreksi_Mixed", summary(aov2))
tampil("29_Multivariat_Mixed", aov2$Anova)
tampil("30_Eta_Kuadrat_Mixed", effectsize::eta_squared(aov2, partial = TRUE))
gambar("G06_Profil_Mixed_ANOVA", afex::afex_plot(aov2, x = "waktu", trace = "kelompok",
error = "within", mapping = c("colour", "shape", "linetype")) +
labs(title = "Profil mixed ANOVA", x = "Waktu", y = "Berat badan (kg)"))
teks("Interval grafik mixed ANOVA", "Interval within-subject pada grafik afex berbeda dari interval kepercayaan rerata biasa pada grafik profil deskriptif.")
tampil("31_Efek_Waktu_per_Kelompok", as.data.frame(emmeans::joint_tests(aov2, by = "kelompok")))
tampil("32_Efek_Kelompok_per_Waktu", as.data.frame(emmeans::joint_tests(aov2, by = "waktu")))
em2 <- emmeans::emmeans(aov2, ~ waktu | kelompok)
em2b <- emmeans::emmeans(aov2, ~ kelompok | waktu)
tampil("33_EMM_Waktu_per_Kelompok", as.data.frame(em2))
tampil("34_Waktu_vs_Baseline_Holm", as.data.frame(emmeans::contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm")))
tampil("35_Antarkelompok_Tukey", as.data.frame(pairs(em2b, adjust = "tukey")))
em_full <- emmeans::emmeans(aov2, ~ waktu * kelompok)
ki <- emmeans::contrast(em_full, interaction = list(
waktu = list("M6-M0" = c(-1, 0, 1)), kelompok = "pairwise"), adjust = "holm")
tampil("36_Kontras_Perbedaan_Perubahan", as.data.frame(ki))
tren <- as.data.frame(emmeans::contrast(em_full,
interaction = c(waktu = "poly", kelompok = "pairwise"), adjust = "none"))
tren_lin <- subset(tren, waktu_poly == "linear")
tren_lin$p.holm <- p.adjust(tren_lin$p.value, "holm")
tampil("37_Kontras_Tren_Linear", tren_lin)
teks("Koreksi perbandingan", "Holm untuk waktu versus baseline berlaku dalam tiap kelompok; Tukey berlaku dalam tiap waktu. Kontras perbedaan perubahan dan tren linear memakai Holm pada masing-masing keluarga kontras.")
tabel_gg <- as.data.frame(afex::nice(aov2, correction = "GG", es = "pes"))
tampil("38_Tabel_ANOVA_GG", tabel_gg)
tab_num <- as.data.frame(aov2$anova_table)
if ("kelompok:waktu" %in% rownames(tab_num)) {
r <- tab_num["kelompok:waktu", , drop = FALSE]
teks("Interpretasi interaksi", paste0("Mixed ANOVA dengan koreksi Greenhouse-Geisser menghasilkan interaksi kelompok x waktu: F(",
round(r[["num Df"]], 2), ", ", round(r[["den Df"]], 2), ") = ", round(r[["F"]], 3),
"; p ", fp(r[["Pr(>F)"]]), ". ",
if (r[["Pr(>F)"]] < .05) "Pola perubahan berat badan berbeda secara statistik antar kelompok."
else "Belum ditemukan bukti statistik perbedaan pola perubahan berat badan antar kelompok.",
" Arah dan besarnya perbedaan dinilai dari rerata serta kontras perubahan M6-M0."))
}
kontrol <- lme4::lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 200000))
lmm1 <- lmerTest::lmer(berat_badan ~ kelompok * waktu + (1 | id), data = dat_lmm, REML = TRUE, control = kontrol)
lmm2 <- lmerTest::lmer(berat_badan ~ kelompok * waktu + (1 + bulan | id), data = dat_lmm, REML = TRUE, control = kontrol)
tampil("39_LMM_Intersep", summary(lmm1))
tampil("40_LMM_Intersep_Slope", summary(lmm2))
tampil("41_Perbandingan_LMM_REML", anova(lmm1, lmm2, refit = FALSE))
teks("Perbandingan model", "Kedua model mempunyai fixed effects yang sama. Perbandingan REML mengikuti materi; nilai p likelihood-ratio bersifat pendekatan karena varians acak diuji pada batas ruang parameter.")
tampil("42_Singularitas_LMM", data.frame(model = c("Intersep", "Intersep dan slope"),
singular = c(lme4::isSingular(lmm1), lme4::isSingular(lmm2))))
opsional("43_LMM_Kenward_Roger", anova(lmm2, ddf = "Kenward-Roger"))
opsional("44_ICC", performance::icc(lmm1))
rd <- data.frame(prediksi = fitted(lmm2), residual = resid(lmm2))
gambar("G07_QQ_Residual_LMM", ggplot(rd, aes(sample = residual)) + stat_qq() + stat_qq_line() +
labs(title = "Q-Q residual LMM", x = "Kuantil teoretis", y = "Kuantil residual"))
ra <- data.frame(intersep = lme4::ranef(lmm2)$id[, 1])
gambar("G08_QQ_Intersep_Acak", ggplot(ra, aes(sample = intersep)) + stat_qq() + stat_qq_line() +
labs(title = "Q-Q intersep acak LMM", x = "Kuantil teoretis", y = "Kuantil intersep acak"))
gambar("G09_Residual_vs_Prediksi", ggplot(rd, aes(prediksi, residual)) + geom_point(alpha = .5) +
geom_hline(yintercept = 0, linetype = 2) + labs(title = "Residual versus nilai prediksi LMM", x = "Prediksi berat badan", y = "Residual"))
em_lmm <- emmeans::emmeans(lmm2, ~ waktu | kelompok, lmer.df = "kenward-roger")
tampil("45_EMM_LMM", as.data.frame(em_lmm))
# Data hilang yang nyata menggantikan demonstrasi penghapusan buatan.
dat_miss <- long_semua
lmm_miss <- lmm2
teks("LMM dan data hilang", "Model LMM utama sudah memakai data hilang yang nyata: baris dengan berat badan NA dikeluarkan, pengukuran tersedia dari peserta tidak lengkap tetap digunakan. Tidak dibuat penghapusan pengukuran secara buatan. Interpretasi LMM mengandalkan asumsi missing at random bersyarat pada model; asumsi ini tidak terbukti hanya dari data ini.")
tampil("47_Ringkasan_Data_Hilang", data.frame(pengukuran_hilang = sum(is.na(long_semua$berat_badan)),
peserta_tidak_lengkap = sum(!lengkap), peserta_lengkap = sum(lengkap),
pengukuran_dianalisis_LMM = nobs(lmm2)))
tampil("48_Informasi_Sesi", sessionInfo())
write.csv(dat_wide, file.path(folder, "Data_Berat_Badan_Ika_Wulan_Sari_Wide.csv"), row.names = FALSE)
write.csv(long_semua, file.path(folder, "Data_Berat_Badan_Ika_Wulan_Sari_Long.csv"), row.names = FALSE)
write.csv(dat_miss, file.path(folder, "Data_Berat_Badan_Ika_Wulan_Sari_Missing_Asli.csv"), row.names = FALSE)
saveRDS(list(hasil = hasil, data_wide = dat_wide, data_long = long_semua, data_anova = dat_long, data_lmm = dat_lmm,
aov1 = aov1, aov2 = aov2, lmm1 = lmm1, lmm2 = lmm2, lmm_miss = lmm_miss),
file.path(folder, "Objek_Analisis_Berat_Badan_Ika_Wulan_Sari.rds"))
laporan <- htmltools::tags$html(lang = "id",
htmltools::tags$head(htmltools::tags$meta(charset = "UTF-8"),
htmltools::tags$title("Analisis Berat Badan Ika Wulan Sari: 2 x 3"),
htmltools::tags$style("body{font-family:Arial,sans-serif;max-width:1100px;margin:40px auto;padding:20px;line-height:1.6;color:#17342e}h1,h2{color:#0b6b4f}h2{border-bottom:1px solid #d9e4e0;padding-top:24px}pre{background:#f3f7f5;padding:16px;overflow:auto;font-size:13px}img{display:block;margin:15px auto}")),
htmltools::tags$body(htmltools::tags$h1("Analisis Pengukuran Berulang Berat Badan Ika Wulan Sari: 2 x 3"), isi))
path_html <- file.path(folder, "Laporan_Analisis_Berat_Badan_Ika_Wulan_Sari.html")
htmltools::save_html(laporan, path_html)
if (BUKA_HTML_OTOMATIS && interactive())
try(utils::browseURL(path_html), silent = TRUE)
cat("\nHasil analisis: ", normalizePath(folder), "\n", sep = "")
invisible(list(hasil = hasil, grafik = grafik, data_wide = dat_wide,
data_long = long_semua, data_anova = dat_long, data_lmm = dat_lmm, aov1 = aov1, aov2 = aov2,
lmm1 = lmm1, lmm2 = lmm2, lmm_miss = lmm_miss))
}
hasil_berat_badan <- analisis_berat_badan()
## Registered S3 method overwritten by 'lme4':
## method from
## na.action.merMod car
## Data: /Users/nazmy/Downloads/journal.pone.0257326.s004 (1).xlsx
## Sheet: data_set
## Hasil: /Users/nazmy/Downloads/Hasil_Analisis_Berat_Badan_Ika_Wulan_Sari
##
## Identitas
## Nama: Ika Wulan Sari | NIM: 2611018005 | Mata Kuliah: Biostatistika Intermediate | Magister Kesehatan Masyarakat Universitas Mulawarman
##
## Desain
## 166 peserta; dua kelompok; tiga waktu (M0, M3, M6). 162 peserta lengkap untuk ANOVA; 5 pengukuran hilang; LMM menggunakan 493 pengukuran tersedia. Outcome: berat badan (kg).
##
## 00 Kelengkapan
## # A tibble: 6 × 5
## kelompok waktu peserta tersedia hilang
## <fct> <fct> <int> <int> <int>
## 1 control M0 81 81 0
## 2 control M3 81 81 0
## 3 control M6 81 79 2
## 4 intervention M0 85 85 0
## 5 intervention M3 85 84 1
## 6 intervention M6 85 83 2
##
## Penanganan data
## ID ROW dibuat dari urutan baris, dengan baris_excel untuk penelusuran; setiap baris diasumsikan satu peserta sesuai struktur wide. Kode -999 dan 999 pada berat badan/usia dibaca sebagai NA. Nilai berat badan 999 pada bulan 6 diperlakukan sebagai kode hilang berdasarkan pola kode pada baris yang sama (usia 999, gender 9, aktivitas -999). Tidak dilakukan imputasi. Kolom turunan overall_PA dan delta_M6_BL tidak digunakan. ANOVA, Friedman, Wilcoxon, dan robust RM memakai peserta lengkap; LMM memakai seluruh pengukuran berat badan tersedia. Hasil kedua pendekatan dapat berbeda karena sampel berbeda.
##
## 01 Jumlah Peserta
## kelompok n
## 1 control 81
## 2 intervention 85
##
## 01a Usia
## # A tibble: 2 × 6
## kelompok n mean sd min max
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 control 81 71.2 5.04 65 81
## 2 intervention 85 70.4 4.61 65 85
##
## 01b Jenis Kelamin
## kelompok jk n
## 1 control F 47
## 2 control M 34
## 3 intervention F 50
## 4 intervention M 35
##
## Berkas data yang digunakan
## /Users/nazmy/Downloads/journal.pone.0257326.s004 (1).xlsx | Sheet: data_set
##
## Sumber data
## Workbook lampiran journal.pone.0257326.s004 (1).xlsx, sheet data_set. Berat badan menggunakan WeightKg, Weightkg.M3, dan Weightkg.M6; kelompok memakai Group. Analisis ini adalah adaptasi metode latihan, bukan reproduksi seluruh analisis artikel.
##
## 02 Data Wide
## id baris_excel kelompok usia jk BB_M0 BB_M3 BB_M6
## 1 ROW001 2 intervention 78 M 84 84 84
## 2 ROW002 3 intervention 65 M 97 NA NA
## 3 ROW003 4 intervention 74 M 92 92 90
## 4 ROW004 5 intervention 69 F 62 62 62
## 5 ROW005 6 intervention 68 F 60 60 NA
## 6 ROW006 7 intervention 75 M 83 83 81
##
## 03 Data Long
## # A tibble: 12 × 8
## id baris_excel kelompok usia jk waktu berat_badan bulan
## <fct> <int> <fct> <dbl> <chr> <fct> <dbl> <dbl>
## 1 ROW001 2 intervention 78 M M0 84 0
## 2 ROW001 2 intervention 78 M M3 84 3
## 3 ROW001 2 intervention 78 M M6 84 6
## 4 ROW003 4 intervention 74 M M0 92 0
## 5 ROW003 4 intervention 74 M M3 92 3
## 6 ROW003 4 intervention 74 M M6 90 6
## 7 ROW004 5 intervention 69 F M0 62 0
## 8 ROW004 5 intervention 69 F M3 62 3
## 9 ROW004 5 intervention 69 F M6 62 6
## 10 ROW006 7 intervention 75 M M0 83 0
## 11 ROW006 7 intervention 75 M M3 83 3
## 12 ROW006 7 intervention 75 M M6 81 6
##
## 04 Deskriptif
## # A tibble: 6 × 12
## kelompok waktu bulan n mean sd median min max se lower upper
## <fct> <fct> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 control M0 0 81 78.0 15.4 75 52 122 1.71 74.6 81.4
## 2 control M3 3 81 78.0 15.5 75 49 122 1.72 74.6 81.4
## 3 control M6 6 79 77.2 15.4 74 52 125 1.73 73.8 80.7
## 4 intervention M0 0 85 77.9 15.1 77 50 115 1.64 74.7 81.2
## 5 intervention M3 3 84 77.2 15.0 76.5 49 115 1.64 74.0 80.5
## 6 intervention M6 6 83 76.7 14.8 76 50 117 1.63 73.5 80.0
##
## Sampel deskriptif
## Tabel deskriptif dan grafik profil rerata menggunakan seluruh berat badan yang tersedia pada tiap waktu. Tabel ANOVA serta grafik profil mixed ANOVA menggunakan peserta dengan tiga pengukuran lengkap.
##
## 05 Kovarians
## BB_M0 BB_M3 BB_M6
## BB_M0 228.2643 226.4029 220.7430
## BB_M3 226.4029 230.2750 223.9765
## BB_M6 220.7430 223.9765 226.6943
##
## 06 Korelasi
## BB_M0 BB_M3 BB_M6
## BB_M0 1.0000000 0.9875058 0.9703933
## BB_M3 0.9875058 1.0000000 0.9802996
## BB_M6 0.9703933 0.9802996 1.0000000
##
## 07 Varians Selisih
## pasangan varians
## 1 BB_M0 - BB_M3 5.733456
## 2 BB_M0 - BB_M6 13.472471
## 3 BB_M3 - BB_M6 9.016218
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 5 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 5 rows containing missing values or values outside the scale
## range (`geom_line()`).
##
## 08 Outlier Satu Kelompok
## # A tibble: 1 × 10
## waktu id baris_excel kelompok usia jk berat_badan bulan is.outlier is.extreme
## <fct> <fct> <int> <fct> <dbl> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 M6 ROW049 50 intervention 67 M 117 6 TRUE FALSE
##
## 09 Shapiro Satu Kelompok
## # A tibble: 3 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 berat_badan 0.980 0.222
## 2 M3 berat_badan 0.980 0.222
## 3 M6 berat_badan 0.980 0.229
##
## 10 ANOVA Rstatix Satu Kelompok
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2 164 4.3 0.015 * 0.05
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.735 3.76e-06 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05
## 1 waktu 0.79 1.58, 129.6 0.023 * 0.803 1.61, 131.72 0.022 *
##
##
## 11 RM ANOVA
## Anova Table (Type 3 tests)
##
## Response: berat_badan
## Effect df MSE F ges pes p.value
## 1 waktu 1.58, 129.60 8.52 4.30 * .001 .050 .023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
##
## 12 Sferisitas dan Koreksi RM
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1490059 1 53785 82 2271.7469 < 2e-16 ***
## waktu 58 2 1105 164 4.3004 0.01513 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.7346 3.7587e-06
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.79026 0.02307 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8031868 0.02247642
##
## 13 Multivariat RM
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.96516 2271.75 1 82 < 2e-16 ***
## waktu 1 0.06571 2.85 2 81 0.06376 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 14 Eta Kuadrat RM
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.05 | [0.01, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
## 15 EMM RM
## waktu emmean SE df lower.CL upper.CL
## M0 77.89157 1.644411 82 74.62031 81.16282
## M3 77.45783 1.645062 82 74.18528 80.73038
## M6 76.72289 1.629268 82 73.48176 79.96402
##
## Confidence level used: 0.95
##
## 16 Posthoc RM Bonferroni
## contrast estimate SE df t.ratio p.value
## M0 - M3 0.4337349 0.3346889 82 1.296 0.5959
## M0 - M6 1.1686747 0.4952965 82 2.360 0.0620
## M3 - M6 0.7349398 0.3600227 82 2.041 0.1333
##
## P value adjustment: bonferroni method for 3 tests
##
## 17 RM vs Baseline Holm
## contrast estimate SE df t.ratio p.value
## M3 - M0 -0.4337349 0.3346889 82 -1.296 0.1986
## M6 - M0 -1.1686747 0.4952965 82 -2.360 0.0413
##
## P value adjustment: holm method for 2 tests
##
## 18 Tren Polinomial RM
## contrast estimate SE df t.ratio p.value
## linear -1.1686747 0.4952965 82 -2.360 0.0207
## quadratic -0.3012048 0.4877985 82 -0.617 0.5386
##
##
## 19 Friedman
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 berat_badan 83 5.17 2 0.0753 Friedman test
##
## 20 Kendall W
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 berat_badan 83 0.0312 Kendall W small
##
## 21 Wilcoxon Berpasangan
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 berat_badan M0 M3 83 83 58 0.128 0.383 ns
## 2 berat_badan M0 M6 83 83 1524 0.0304 0.0913 ns
## 3 berat_badan M3 M6 83 83 1410. 0.0504 0.151 ns
##
## 22 Robust RM Trim20
## Call:
## WRS2::rmanova(y = d1$berat_badan, groups = d1$waktu, blocks = d1$id,
## tr = 0.2)
##
## Test statistic: F = 3.2291
## Degrees of freedom 1: 1.68
## Degrees of freedom 2: 83.89
## p-value: 0.05298
##
##
## 23 Outlier Mixed
## # A tibble: 8 × 10
## kelompok waktu id baris_excel usia jk berat_badan bulan is.outlier is.extreme
## <fct> <fct> <fct> <int> <dbl> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 control M0 ROW140 141 68 M 122 0 TRUE FALSE
## 2 control M3 ROW140 141 68 M 122 3 TRUE FALSE
## 3 control M6 ROW092 93 65 F 112 6 TRUE FALSE
## 4 control M6 ROW127 128 72 M 111 6 TRUE FALSE
## 5 control M6 ROW140 141 68 M 125 6 TRUE FALSE
## 6 control M6 ROW144 145 69 M 107 6 TRUE FALSE
## 7 control M6 ROW163 164 78 F 110 6 TRUE FALSE
## 8 intervention M6 ROW049 50 67 M 117 6 TRUE FALSE
##
## 24 Shapiro Mixed
## # A tibble: 6 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 control M0 berat_badan 0.956 0.00786
## 2 control M3 berat_badan 0.965 0.0294
## 3 control M6 berat_badan 0.954 0.00662
## 4 intervention M0 berat_badan 0.980 0.222
## 5 intervention M3 berat_badan 0.980 0.222
## 6 intervention M6 berat_badan 0.980 0.229
##
## 25 Levene
## # A tibble: 3 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 1 160 0.0311 0.860
## 2 M3 1 160 0.0169 0.897
## 3 M6 1 160 0.0813 0.776
##
## 26 Box M
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 58.5 9.10e-11 6 Box's M-test for Homogeneity of Covariance Matrices
##
## 27 Mixed ANOVA
## Anova Table (Type 3 tests)
##
## Response: berat_badan
## Effect df MSE F ges pes p.value
## 1 kelompok 1, 160 680.04 0.00 <.001 <.001 .954
## 2 waktu 1.63, 261.28 5.75 5.84 ** <.001 .035 .006
## 3 kelompok:waktu 1.63, 261.28 5.75 1.18 <.001 .007 .302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
##
## 28 Sferisitas dan Koreksi Mixed
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 2911656 1 108806 160 4281.6198 < 2.2e-16 ***
## kelompok 2 1 108806 160 0.0033 0.954163
## waktu 55 2 1504 320 5.8397 0.003229 **
## kelompok:waktu 11 2 1504 320 1.1783 0.309124
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.77527 1.6271e-09
## kelompok:waktu 0.77527 1.6271e-09
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.8165 0.005903 **
## kelompok:waktu 0.8165 0.302380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8237657 0.005763773
## kelompok:waktu 0.8237657 0.302710119
##
## 29 Multivariat Mixed
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.96398 4281.6 1 160 < 2e-16 ***
## kelompok 1 0.00002 0.0 1 160 0.95416
## waktu 1 0.04721 3.9 2 159 0.02139 *
## kelompok:waktu 1 0.01127 0.9 2 159 0.40601
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 30 Eta Kuadrat Mixed
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 2.07e-05 | [0.00, 1.00]
## waktu | 0.04 | [0.01, 1.00]
## kelompok:waktu | 7.31e-03 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
## Warning: Panel(s) show a mixed within-between-design.
## Error bars do not allow comparisons across all means.
## Suppress error bars with: error = "none"
##
## Interval grafik mixed ANOVA
## Interval within-subject pada grafik afex berbeda dari interval kepercayaan rerata biasa pada grafik profil deskriptif.
##
## 31 Efek Waktu per Kelompok
## kelompok = control:
## model term df1 df2 F.ratio p.value
## waktu 2 160 0.738 0.4797
##
## kelompok = intervention:
## model term df1 df2 F.ratio p.value
## waktu 2 160 4.224 0.0163
##
##
## 32 Efek Kelompok per Waktu
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 1 160 0.011 0.9179
##
## waktu = M3:
## model term df1 df2 F.ratio p.value
## kelompok 1 160 0.005 0.9460
##
## waktu = M6:
## model term df1 df2 F.ratio p.value
## kelompok 1 160 0.043 0.8360
##
##
## 33 EMM Waktu per Kelompok
## kelompok = control:
## waktu emmean SE df lower.CL upper.CL
## M0 77.64557 1.705077 160 74.27821 81.01293
## M3 77.62025 1.712602 160 74.23803 81.00247
## M6 77.21519 1.699031 160 73.85977 80.57061
##
## kelompok = intervention:
## waktu emmean SE df lower.CL upper.CL
## M0 77.89157 1.663483 160 74.60635 81.17678
## M3 77.45783 1.670825 160 74.15812 80.75755
## M6 76.72289 1.657585 160 73.44932 79.99646
##
## Confidence level used: 0.95
##
## 34 Waktu vs Baseline Holm
## kelompok = control:
## contrast estimate SE df t.ratio p.value
## M3 - M0 -0.0253165 0.2692487 160 -0.094 0.9252
## M6 - M0 -0.4303797 0.4121384 160 -1.044 0.5959
##
## kelompok = intervention:
## contrast estimate SE df t.ratio p.value
## M3 - M0 -0.4337349 0.2626806 160 -1.651 0.1007
## M6 - M0 -1.1686747 0.4020847 160 -2.907 0.0083
##
## P value adjustment: holm method for 2 tests
##
## 35 Antarkelompok Tukey
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## control - intervention -0.2459966 2.382113 160 -0.103 0.9179
##
## waktu = M3:
## contrast estimate SE df t.ratio p.value
## control - intervention 0.1624218 2.392627 160 0.068 0.9460
##
## waktu = M6:
## contrast estimate SE df t.ratio p.value
## control - intervention 0.4922983 2.373667 160 0.207 0.8360
##
##
## 36 Kontras Perbedaan Perubahan
## waktu_custom kelompok_pairwise estimate SE df t.ratio p.value
## M6-M0 control - intervention 0.738295 0.5757866 160 1.282 0.2016
##
##
## 37 Kontras Tren Linear
## waktu_poly kelompok_pairwise estimate SE df t.ratio p.value p.holm
## 1 linear control - intervention 0.738295 0.5757866 160 1.282237 0.2016141 0.2016141
##
## Koreksi perbandingan
## Holm untuk waktu versus baseline berlaku dalam tiap kelompok; Tukey berlaku dalam tiap waktu. Kontras perbedaan perubahan dan tren linear memakai Holm pada masing-masing keluarga kontras.
##
## 38 Tabel ANOVA GG
## Effect df MSE F pes p.value
## 1 kelompok 1, 160 680.04 0.00 <.001 .954
## 2 waktu 1.63, 261.28 5.75 5.84 ** .035 .006
## 3 kelompok:waktu 1.63, 261.28 5.75 1.18 .007 .302
##
## Interpretasi interaksi
## Mixed ANOVA dengan koreksi Greenhouse-Geisser menghasilkan interaksi kelompok x waktu: F(1.63, 261.28) = 1.178; p = 0,302. Belum ditemukan bukti statistik perbedaan pola perubahan berat badan antar kelompok. Arah dan besarnya perbedaan dinilai dari rerata serta kontras perubahan M6-M0.
##
## 39 LMM Intersep
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: berat_badan ~ kelompok * waktu + (1 | id)
## Data: dat_lmm
## Control: kontrol
##
## REML criterion at convergence: 2982.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.9171 -0.2753 -0.0353 0.3056 8.2191
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 226.983 15.066
## Residual 4.655 2.158
## Number of obs: 493, groups: id, 166
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 77.61576 1.17376 163.97873 66.126 <2e-16 ***
## kelompok1 0.24376 1.17376 163.97873 0.208 0.836
## waktu1 0.34336 0.13749 323.03899 2.497 0.013 *
## waktu2 0.11444 0.13749 323.00216 0.832 0.406
## kelompok1:waktu1 -0.19053 0.13749 323.03899 -1.386 0.167
## kelompok1:waktu2 0.01369 0.13749 323.00216 0.100 0.921
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) klmpk1 waktu1 waktu2 klm1:1
## kelompok1 0.024
## waktu1 -0.001 0.000
## waktu2 0.000 -0.001 -0.493
## klmpk1:wkt1 0.000 -0.001 0.020 -0.007
## klmpk1:wkt2 -0.001 0.000 -0.007 0.020 -0.493
##
## 40 LMM Intersep Slope
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: berat_badan ~ kelompok * waktu + (1 + bulan | id)
## Data: dat_lmm
## Control: kontrol
##
## REML criterion at convergence: 2949.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4050 -0.1686 -0.0384 0.1796 5.9556
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 231.0333 15.200
## bulan 0.2219 0.471 -0.13
## Residual 2.6604 1.631
## Number of obs: 493, groups: id, 166
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 77.61539 1.17365 163.97447 66.132 <2e-16 ***
## kelompok1 0.24375 1.17365 163.97447 0.208 0.8357
## waktu1 0.34373 0.15149 204.84543 2.269 0.0243 *
## waktu2 0.11448 0.10413 163.87485 1.099 0.2732
## kelompok1:waktu1 -0.19052 0.15149 204.84543 -1.258 0.2099
## kelompok1:waktu2 0.01404 0.10413 163.87485 0.135 0.8929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) klmpk1 waktu1 waktu2 klm1:1
## kelompok1 0.024
## waktu1 0.023 0.001
## waktu2 0.000 -0.001 -0.335
## klmpk1:wkt1 0.001 0.023 0.020 -0.004
## klmpk1:wkt2 -0.001 0.000 -0.004 0.021 -0.335
##
## 41 Perbandingan LMM REML
## Data: dat_lmm
## Models:
## lmm1: berat_badan ~ kelompok * waktu + (1 | id)
## lmm2: berat_badan ~ kelompok * waktu + (1 + bulan | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 8 2998.5 3032.1 -1491.2 2982.5
## lmm2 10 2969.3 3011.3 -1474.7 2949.3 33.209 2 6.15e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Perbandingan model
## Kedua model mempunyai fixed effects yang sama. Perbandingan REML mengikuti materi; nilai p likelihood-ratio bersifat pendekatan karena varians acak diuji pada batas ruang parameter.
##
## 42 Singularitas LMM
## model singular
## 1 Intersep FALSE
## 2 Intersep dan slope FALSE
##
## 43 LMM 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.1148 0.1148 1 164.00 0.0431 0.8357
## waktu 24.0500 12.0250 2 215.34 4.5061 0.0121 *
## kelompok:waktu 4.4532 2.2266 2 215.34 0.8344 0.4356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 44 ICC
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.980
## Unadjusted ICC: 0.979
##
## 45 EMM LMM
## kelompok = control:
## waktu emmean SE df lower.CL upper.CL
## M0 78.01235 1.698561 164.61 74.65857 81.36612
## M3 77.98765 1.686043 166.49 74.65887 81.31643
## M6 77.57742 1.688780 164.64 74.24296 80.91188
##
## kelompok = intervention:
## waktu emmean SE df lower.CL upper.CL
## M0 77.90588 1.658113 164.61 74.63197 81.17979
## M3 77.47208 1.646201 166.61 74.22198 80.72219
## M6 76.73695 1.648781 164.72 73.48148 79.99242
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## LMM dan data hilang
## Model LMM utama sudah memakai data hilang yang nyata: baris dengan berat badan NA dikeluarkan, pengukuran tersedia dari peserta tidak lengkap tetap digunakan. Tidak dibuat penghapusan pengukuran secara buatan. Interpretasi LMM mengandalkan asumsi missing at random bersyarat pada model; asumsi ini tidak terbukti hanya dari data ini.
##
## 47 Ringkasan Data Hilang
## pengukuran_hilang peserta_tidak_lengkap peserta_lengkap pengukuran_dianalisis_LMM
## 1 5 4 162 493
##
## 48 Informasi Sesi
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Asia/Singapore
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_4.0.3 tidyr_1.3.2 dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 farver_2.1.2 S7_0.2.2 fastmap_1.2.0 reshape_0.8.10
## [6] bayestestR_0.19.0 digest_0.6.39 estimability_2.0.0 lifecycle_1.0.5 magrittr_2.0.5
## [11] compiler_4.6.1 rlang_1.3.0 sass_0.4.10 tools_4.6.1 utf8_1.2.6
## [16] yaml_2.3.12 knitr_1.52 ggsignif_0.6.4 labeling_0.4.3 plyr_1.8.9
## [21] RColorBrewer_1.1-3 abind_1.4-8 withr_3.0.3 purrr_1.2.2 numDeriv_2016.8-1.1
## [26] grid_4.6.1 afex_1.5-1 datawizard_1.4.0 ggpubr_1.0.0 emmeans_2.0.4
## [31] scales_1.4.0 MASS_7.3-66 insight_1.5.4 cli_3.6.6 mvtnorm_1.4-2
## [36] rmarkdown_2.32 ragg_1.5.2 reformulas_0.4.4 generics_0.1.4 otel_0.2.0
## [41] rstudioapi_0.19.0 performance_0.18.2 reshape2_1.4.5 parameters_0.29.3 readxl_1.5.0.1
## [46] minqa_1.2.8 cachem_1.1.0 stringr_1.6.0 splines_4.6.1 parallel_4.6.1
## [51] effectsize_1.0.3 cellranger_1.1.0 WRS2_1.1-7 base64enc_0.1-6 vctrs_0.7.3
## [56] boot_1.3-32 Matrix_1.7-6 jsonlite_2.0.0 carData_3.0-6 car_3.1-5
## [61] pbkrtest_0.5.5 rstatix_1.1.0 Formula_1.2-6 systemfonts_1.3.2 jquerylib_0.1.4
## [66] glue_1.8.1 nloptr_2.2.1 stringi_1.8.9 gtable_0.3.6 lme4_2.0-6
## [71] lmerTest_3.2-1 tibble_3.3.1 pillar_1.11.1 htmltools_0.5.9 R6_2.6.1
## [76] textshaping_1.0.5 Rdpack_2.6.6 evaluate_1.0.5 lattice_0.23-1 rbibutils_2.4.1
## [81] backports_1.5.1 broom_1.0.13 bslib_0.12.0 Rcpp_1.1.2 nlme_3.1-171
## [86] xfun_0.61 pkgconfig_2.0.3
##
## Hasil analisis: /Users/nazmy/Downloads/Hasil_Analisis_Berat_Badan_Ika_Wulan_Sari