# TUGAS BIOSTATISTIKA INTERMEDIETE
# ANALISIS PENGUKURAN BERULANG TEKANAN DARAH SISTOLIK
# Nama : Muhammad Ridhan Nazmy
# NIM : 2611018008
# Rpubs : r_nazmy
# Mata Kuliah : Biostatistika Intermediete
# Magister Kesehatan Masyarakat Universitas Mulawarman
#
# Sumber data: James dkk. (2019), PLOS Medicine, S1 Data.
# Artikel: https://doi.org/10.1371/journal.pmed.1002870
# Data: https://doi.org/10.1371/journal.pmed.1002870.s007
# Analisis sekunder menggunakan RM ANOVA, Mixed Design ANOVA, dan LMM.
# ======================= PARAMETER ANALISIS =========================
FILE_DATA <- "pmed.1002870.s007.xls"
KELOMPOK_FOKUS <- "Drink powder"
FOLDER_HASIL <- "hasil_analisis_tds"
INSTALL_PAKET_OTOMATIS <- TRUE
SIMULASI_MISSING_TAMBAHAN <- FALSE # Data asli sudah memiliki nilai hilang.
BUKA_HTML_OTOMATIS <- TRUE
TAMPILKAN_GRAFIK_RSTUDIO <- TRUE # Tampilan grafik pada perangkat grafis aktif.
# ======================== PROGRAM UTAMA ==============================
jalankan_analisis <- function() {
waktu_mulai <- Sys.time()
paket <- c("readxl", "dplyr", "tidyr", "ggplot2", "afex", "emmeans",
"rstatix", "car", "effectsize", "lme4", "lmerTest", "pbkrtest",
"performance", "ggpubr", "knitr", "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 = ", "))
message("Menginstal paket: ", paste(belum, collapse = ", "))
install.packages(belum, repos = "https://cloud.r-project.org")
}
gagal_paket <- paket[!vapply(paket, requireNamespace, logical(1), quietly = TRUE)]
if (length(gagal_paket)) stop("Instalasi gagal: ", paste(gagal_paket, collapse = ", "))
suppressPackageStartupMessages({
library(dplyr); library(tidyr); library(ggplot2)
library(afex); library(emmeans); library(rstatix)
library(car); library(effectsize); library(lme4); library(lmerTest)
})
opsi_lama <- options(contrasts = c("contr.sum", "contr.poly"))
on.exit(options(opsi_lama), add = TRUE)
# Resolusi lokasi berkas data berdasarkan direktori 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)
}
kandidat <- unique(c(FILE_DATA, file.path(dirname(sumber), FILE_DATA)))
kandidat <- kandidat[file.exists(kandidat)]
if (!length(kandidat)) stop(
"Berkas data tidak ditemukan pada direktori script atau direktori kerja.\n",
"Nama file yang dicari: ", FILE_DATA)
path_data <- normalizePath(kandidat[1], winslash = "/", mustWork = TRUE)
folder <- file.path(dirname(path_data), FOLDER_HASIL)
dir.create(folder, recursive = TRUE, showWarnings = FALSE)
folder <- normalizePath(folder, winslash = "/", mustWork = TRUE)
message("Data: ", path_data, "\nHasil: ", folder)
# Penyusunan laporan HTML, tabel, dan visualisasi.
isi <- character(0); status <- data.frame()
tambah <- function(html) isi <<- c(isi, html)
escape <- function(x) as.character(htmltools::htmlEscape(paste(x, collapse = "\n")))
narasi <- function(x) tambah(paste0("<p>", escape(x), "</p>"))
output <- function(x) {
teks <- capture.output(print(x))
cat(paste(teks, collapse = "\n"), "\n")
tambah(paste0("<pre>", escape(teks), "</pre>"))
invisible(x)
}
tabel <- function(x, nama) {
x <- as.data.frame(x)
write.csv(x, file.path(folder, paste0(nama, ".csv")), row.names = FALSE,
na = "", fileEncoding = "UTF-8")
if (nrow(x)) {
tambah(paste0("<div class='table'>",
as.character(knitr::kable(x, format = "html", digits = 4)), "</div>"))
} else narasi("Tidak ada baris untuk tabel ini.")
print(x)
invisible(x)
}
gambar <- function(p, nama, lebar = 10, tinggi = 6) {
path <- file.path(folder, paste0(nama, ".png"))
ggplot2::ggsave(path, plot = p, width = lebar, height = tinggi, dpi = 180)
# Pencetakan eksplisit objek ggplot pada perangkat grafis aktif.
# Grafik diekspor ke PNG dan disematkan pada laporan HTML.
if (TAMPILKAN_GRAFIK_RSTUDIO) {
tryCatch(print(p), error = function(e) {
message("Tampilan Plots gagal, PNG tetap tersimpan: ", conditionMessage(e))
})
}
tambah(paste0("<figure><img src='", base64enc::dataURI(file = path, mime = "image/png"),
"' alt='", escape(nama), "'><figcaption>", escape(nama), "</figcaption></figure>"))
message("Grafik tersimpan: ", path)
}
bagian <- function(judul, f) {
message("\n>>> ", judul)
tambah(paste0("<h2>", escape(judul), "</h2>"))
warnings <- character(0)
ok <- tryCatch(withCallingHandlers({f(); TRUE}, warning = function(w) {
warnings <<- c(warnings, conditionMessage(w)); invokeRestart("muffleWarning")
}), error = function(e) {
narasi(paste("BAGIAN GAGAL:", conditionMessage(e)))
message("Bagian gagal: ", conditionMessage(e)); FALSE
})
if (length(warnings)) {
narasi(paste("Peringatan estimasi:", paste(unique(warnings), collapse = "; ")))
}
status <<- rbind(status, data.frame(bagian = judul,
status = if (ok) "Selesai" else "Gagal",
peringatan = paste(unique(warnings), collapse = "; ")))
invisible(ok)
}
fp <- function(p) {
if (is.na(p)) return("tidak tersedia")
if (p < .001) "< 0,001" else paste0("= ", formatC(p, digits = 3, format = "f", decimal.mark = ","))
}
ringkas_interaksi <- function(tab, pcol, model) {
idx <- which(grepl("kelompok:waktu|waktu:kelompok", rownames(tab)))
if (!length(idx)) {narasi("Efek interaksi tidak teridentifikasi secara otomatis pada tabel model."); return(invisible(NULL))}
p <- as.numeric(tab[idx[1], pcol])
narasi(paste0(model, ": interaksi kelompok × waktu memiliki p ", fp(p), ". ",
if (is.na(p)) "Pengujian tidak dapat disimpulkan." else if (p < .05)
"Terdapat bukti bahwa pola perubahan TDS berbeda antarkelompok. Perbedaan perubahan antarkelompok dievaluasi melalui kontras interaksi."
else "Belum ditemukan bukti bahwa pola perubahan TDS berbeda antarkelompok. Ini tidak membuktikan kelompok setara."))
}
# ========================= 1. MEMBACA DATA ==========================
raw <- as.data.frame(readxl::read_excel(path_data, sheet = 1))
names(raw) <- trimws(names(raw))
kolom_tds <- c("Systolic BP baseline", "Systolic BP midline", "Systolic BP endline")
wajib <- c("Study ID", "Intervention", "Age at baseline (years)", "BMI Baseline", kolom_tds)
if (!all(wajib %in% names(raw))) stop("Kolom tidak ditemukan: ",
paste(setdiff(wajib, names(raw)), collapse = ", "))
level_kelompok <- c("Control", "UNIMMAP", "Drink powder")
label_waktu <- c("M0", "M5", "M12")
minggu_nilai <- c(0, 5, 12)
kel <- trimws(as.character(raw[["Intervention"]]))
if (anyNA(kel) || any(!kel %in% level_kelompok)) stop("Label kelompok tidak dikenali.")
if (anyNA(raw[["Study ID"]]) || anyDuplicated(raw[["Study ID"]])) stop("Study ID hilang atau duplikat.")
if (!KELOMPOK_FOKUS %in% level_kelompok) stop("KELOMPOK_FOKUS tidak dikenali.")
if (!all(vapply(raw[kolom_tds], is.numeric, logical(1)))) stop("TDS harus numerik.")
if (any(!is.finite(as.matrix(raw[kolom_tds])) & !is.na(as.matrix(raw[kolom_tds]))))
stop("TDS memiliki nilai tak hingga.")
dat_wide <- data.frame(id = factor(raw[["Study ID"]]),
kelompok = factor(kel, levels = level_kelompok),
usia = raw[["Age at baseline (years)"]], bmi_baseline = raw[["BMI Baseline"]],
TDS_M0 = raw[[kolom_tds[1]]], TDS_M5 = raw[[kolom_tds[2]]], TDS_M12 = raw[[kolom_tds[3]]])
tds_names <- c("TDS_M0", "TDS_M5", "TDS_M12")
dat_wide$lengkap <- complete.cases(dat_wide[tds_names])
dat_long <- dat_wide |>
pivot_longer(all_of(tds_names), names_to = "waktu", values_to = "tds") |>
mutate(waktu = factor(waktu, levels = tds_names, labels = label_waktu),
minggu = minggu_nilai[as.integer(waktu)]) |>
arrange(id, waktu)
wide_cc <- droplevels(dat_wide[dat_wide$lengkap, ])
long_cc <- droplevels(filter(dat_long, lengkap))
long_obs <- droplevels(filter(dat_long, !is.na(tds)))
d1 <- droplevels(filter(long_cc, kelompok == KELOMPOK_FOKUS))
if (nrow(wide_cc) < 10 || n_distinct(d1$id) < 4) stop("Peserta lengkap terlalu sedikit.")
write.csv(dat_wide, file.path(folder, "data_wide.csv"), row.names = FALSE, na = "")
write.csv(dat_long, file.path(folder, "data_long.csv"), row.names = FALSE, na = "")
write.csv(wide_cc, file.path(folder, "data_complete_case.csv"), row.names = FALSE, na = "")
bagian("1. Sumber, desain, dan kelengkapan data", function() {
narasi("James dkk. (2019), PLOS Medicine, DOI 10.1371/journal.pmed.1002870. S1 Data: 10.1371/journal.pmed.1002870.s007. Analisis ini menggunakan luaran sekunder tekanan darah sistolik (mmHg). Populasi penelitian adalah perempuan usia reproduktif.")
narasi("Desain analisis terdiri atas 3 kelompok dan 3 waktu pengukuran (minggu 0, 5, dan 12). RM ANOVA, Mixed Design ANOVA, dan LMM diterapkan sebagai analisis sekunder. Artikel sumber menggunakan regresi dengan penyesuaian baseline; perbedaan spesifikasi model dapat menghasilkan estimasi dan nilai p yang berbeda.")
narasi(paste(nrow(dat_wide), "peserta;", nrow(wide_cc), "peserta lengkap untuk ANOVA;",
nrow(long_obs), "pengukuran tersedia untuk LMM. Nilai hilang asli tidak diimputasi dan tidak dihapus dari berkas sumber."))
tabel(dat_long |> group_by(kelompok, waktu) |>
summarise(n_total = n(), n_tersedia = sum(!is.na(tds)),
n_hilang = sum(is.na(tds)), .groups = "drop"), "01_kelengkapan")
tabel(dat_wide |> group_by(kelompok) |>
summarise(n_total = n(), n_lengkap = sum(lengkap),
n_tidak_lengkap = sum(!lengkap), .groups = "drop"), "01_peserta")
narasi("ANOVA memakai peserta dengan ketiga TDS tersedia. LMM memakai semua TDS yang teramati, dengan asumsi mekanisme missing yang dapat diabaikan (misalnya MAR bersyarat pada model); asumsi ini tidak otomatis terbukti dari dataset.")
})
# ======================= 2. EKSPLORASI ==============================
bagian("2. Statistik deskriptif dan visualisasi", function() {
desk <- function(d) d |> group_by(kelompok, waktu) |>
summarise(n = n(), mean = mean(tds), sd = sd(tds),
median = median(tds), min = min(tds), max = max(tds),
se = sd/sqrt(n), lower_ci = mean - qt(.975, n-1)*se,
upper_ci = mean + qt(.975, n-1)*se, .groups = "drop")
tabel(desk(long_obs), "02_deskriptif_semua_tersedia")
tabel(desk(long_cc), "02_deskriptif_complete_case")
gambar(ggplot(long_obs, aes(minggu, tds, colour = kelompok, group = kelompok)) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", size = 2) +
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = .3) +
scale_x_continuous(breaks = minggu_nilai) + theme_bw(base_size = 12) +
theme(legend.position = "bottom") +
labs(title = "TDS: semua pengukuran tersedia", x = "Minggu", y = "TDS (mmHg)", colour = "Kelompok"),
"02_profil_tds")
narasi("Interval kepercayaan pada grafik profil adalah CI rerata per sel, bukan CI selisih berpasangan. Peserta yang menyumbang rerata dapat berbeda antarwaktu.")
gambar(ggplot(long_cc, aes(minggu, tds, group = id)) + geom_line(alpha = .15) +
stat_summary(aes(group = kelompok), fun = mean, geom = "line", colour = "firebrick", linewidth = 1.1) +
facet_wrap(~kelompok) + scale_x_continuous(breaks = minggu_nilai) + theme_bw() +
labs(title = "Lintasan individu: peserta lengkap", x = "Minggu", y = "TDS (mmHg)"), "02_lintasan")
for (g in level_kelompok) {
z <- wide_cc[wide_cc$kelompok == g, tds_names]
narasi(paste("Kovarians dan korelasi dalam kelompok:", g))
output(round(cov(z), 3)); output(round(cor(z), 3))
pasangan <- combn(tds_names, 2)
tabel(data.frame(kelompok = g,
pasangan = apply(pasangan, 2, paste, collapse = " - "),
var_selisih = apply(pasangan, 2, function(p) var(z[[p[1]]] - z[[p[2]]]))),
paste0("02_varians_selisih_", make.names(g)))
}
narasi("Matriks dihitung per kelompok agar perbedaan rerata intervensi tidak tercampur dengan kovarians dalam kelompok. Ketidaksamaan varians selisih merupakan petunjuk; uji Mauchly ada pada bagian ANOVA.")
})
# ======================== 3. ASUMSI =================================
bagian("3. Asumsi ANOVA pada peserta lengkap", function() {
tabel(long_cc |> group_by(kelompok, waktu) |> identify_outliers(tds), "03_outlier")
tabel(long_cc |> group_by(kelompok, waktu) |> shapiro_test(tds), "03_shapiro")
gambar(ggpubr::ggqqplot(long_cc, "tds", facet.by = c("kelompok", "waktu")),
"03_qq_per_sel", 12, 9)
tabel(long_cc |> group_by(waktu) |> levene_test(tds ~ kelompok), "03_levene")
tabel(rstatix::box_m(wide_cc[tds_names], wide_cc$kelompok), "03_box_m")
narasi("p > 0,05 tidak membuktikan normalitas atau kesamaan varians. Observasi yang teridentifikasi sebagai outlier dipertahankan dalam analisis; identifikasi statistik tidak dengan sendirinya membuktikan kesalahan pengukuran. Kelompok lengkap tidak sama besar. Jika Levene/Box M mengindikasikan heterogenitas kuat, ANOVA dan model lmer dengan residual homogen perlu ditinjau; koreksi GG tidak menyelesaikan heterogenitas antarkelompok.")
})
# ====================== 4. RM ANOVA SATU ARAH ========================
aov1 <- NULL; aov2 <- NULL
bagian("4. RM ANOVA satu kelompok", function() {
narasi(paste("Kelompok fokus:", KELOMPOK_FOKUS, "—", n_distinct(d1$id), "peserta lengkap."))
gambar(ggpubr::ggqqplot(d1, "tds", facet.by = "waktu") +
labs(title = paste("Q–Q TDS per waktu:", KELOMPOK_FOKUS)), "04_qq_kelompok_fokus", 10, 5)
output(rstatix::anova_test(data = d1, dv = tds, wid = id, within = waktu, effect.size = "pes"))
aov1 <<- afex::aov_ez(id = "id", dv = "tds", data = d1, within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
output(aov1); output(summary(aov1)); output(aov1$Anova)
tabel(as.data.frame(aov1$anova_table) |> tibble::rownames_to_column("efek"), "04_anova_gg")
output(effectsize::eta_squared(aov1, partial = TRUE))
p <- as.numeric(aov1$anova_table[1, "Pr(>F)"])
narasi(paste0("Efek waktu dengan koreksi GG: p ", fp(p), ". ",
if (!is.na(p) && p < .05) "Ada bukti perbedaan rerata TDS antarwaktu dalam kelompok fokus."
else "Belum ada bukti perbedaan rerata TDS antarwaktu dalam kelompok fokus."))
narasi("Tabel utama memakai GG secara konsisten. Ringkasan memuat hasil tanpa koreksi, Mauchly, serta GG/HF. Perubahan dalam satu kelompok saja tidak membuktikan efek intervensi dibanding kontrol.")
})
bagian("4b. Grafik RM ANOVA satu kelompok", function() {
if (is.null(aov1)) stop("Model RM ANOVA belum tersedia.")
gambar(afex::afex_plot(aov1, x = "waktu", error = "within") +
labs(title = paste("RM ANOVA:", KELOMPOK_FOKUS), x = "Waktu", y = "TDS (mmHg)"),
"04_profil_rm_anova")
narasi("Error bar within-subject menggambarkan ketidakpastian untuk perbandingan dalam peserta. Interval tersebut bukan CI perbedaan antarkelompok.")
})
bagian("5. Post hoc satu kelompok dan tren waktu", function() {
if (is.null(aov1)) stop("Model RM ANOVA belum tersedia.")
em <- emmeans::emmeans(aov1, ~ waktu, model = "multivariate")
tabel(summary(em), "05_emmeans")
tabel(summary(pairs(em, adjust = "bonferroni"), infer = c(TRUE, TRUE)), "05_pasangan_bonferroni")
tabel(summary(contrast(em, "trt.vs.ctrl", ref = 1, adjust = "holm"),
infer = c(TRUE, TRUE)), "05_vs_baseline_holm")
# Waktu tidak berjarak sama: koefisien ortogonal berdasarkan minggu asli.
pol <- poly(minggu_nilai, degree = 2)
koef <- list(linear_minggu = pol[, 1], kuadratik_minggu = pol[, 2])
tabel(summary(contrast(em, method = koef, adjust = "holm"), infer = c(TRUE, TRUE)), "05_tren")
narasi("Tren memakai minggu 0, 5, 12. Koefisien ortogonal dinormalisasi, sehingga estimate linear bukan kemiringan mmHg per minggu. Hanya linear dan kuadratik dapat dihitung; tren kubik tidak tersedia. CI dan p disesuaikan sesuai metode yang tertera pada output.")
})
# ====================== 6. NONPARAMETRIK ============================
bagian("6. Pembanding nonparametrik satu kelompok", function() {
tabel(rstatix::friedman_test(d1, tds ~ waktu | id), "06_friedman")
tabel(rstatix::friedman_effsize(d1, tds ~ waktu | id), "06_kendall_w")
# Penyusunan pasangan pengukuran berdasarkan ID peserta.
w <- wide_cc[wide_cc$kelompok == KELOMPOK_FOKUS, ]
ij <- combn(seq_along(tds_names), 2)
hasil <- lapply(seq_len(ncol(ij)), function(j) {
a <- ij[1, j]; b <- ij[2, j]
wt <- wilcox.test(w[[tds_names[a]]], w[[tds_names[b]]], paired = TRUE, exact = FALSE)
data.frame(waktu_1 = label_waktu[a], waktu_2 = label_waktu[b],
n_pasangan = nrow(w), V = unname(wt$statistic), p = wt$p.value)
}) |> bind_rows() |> mutate(p_bonferroni = p.adjust(p, "bonferroni"))
tabel(hasil, "06_wilcoxon")
narasi("Friedman merupakan pembanding untuk efek waktu dalam satu kelompok; bukan pengganti pengujian interaksi kelompok × waktu. Wilcoxon signed-rank juga memiliki asumsi mengenai distribusi selisih. Pilihan metode perlu didasarkan pada desain dan diagnostik, bukan dipilih menurut nilai p yang lebih kecil.")
})
# ======================== 7. MIXED ANOVA ============================
bagian("7. Mixed Design ANOVA kelompok × waktu", function() {
aov2 <<- afex::aov_ez(id = "id", dv = "tds", data = long_cc,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
output(aov2); output(summary(aov2)); output(aov2$Anova)
tabel(as.data.frame(aov2$anova_table) |> tibble::rownames_to_column("efek"), "07_anova_gg")
output(effectsize::eta_squared(aov2, partial = TRUE))
ringkas_interaksi(aov2$anova_table, "Pr(>F)", "Mixed Design ANOVA")
narasi("Efek kelompok adalah perbedaan yang dirata-ratakan atas waktu; efek waktu dirata-ratakan atas kelompok. Interpretasi interaksi didasarkan pada profil kelompok dan kontras perubahan. Baseline p > 0,05 bukan bukti bahwa randomisasi berhasil atau kelompok benar-benar setara.")
})
bagian("7b. Visualisasi interaksi kelompok dan waktu", function() {
if (is.null(aov2)) stop("Model Mixed Design ANOVA belum tersedia.")
gambar(afex::afex_plot(aov2, x = "waktu", trace = "kelompok", error = "within",
mapping = c("colour", "shape", "linetype")) +
labs(title = "Mixed Design ANOVA: kelompok × waktu", x = "Waktu", y = "TDS (mmHg)"),
"07_interaksi_anova")
narasi("Visualisasi interaksi menggunakan afex_plot pada peserta dengan pengukuran lengkap. Grafik menyajikan rerata model, profil kelompok, dan error bar within-subject.")
})
# ======================= 8. EFEK SEDERHANA ==========================
bagian("8. Efek sederhana dan post hoc antarkelompok", function() {
if (is.null(aov2)) stop("Model Mixed Design ANOVA belum tersedia.")
jt1 <- as.data.frame(joint_tests(aov2, by = "kelompok", model = "multivariate"))
jt1$p_holm_antar_kelompok <- p.adjust(jt1$p.value, "holm")
tabel(jt1, "08_efek_waktu_per_kelompok")
jt2 <- as.data.frame(joint_tests(aov2, by = "waktu", model = "multivariate"))
jt2$p_holm_antar_waktu <- p.adjust(jt2$p.value, "holm")
tabel(jt2, "08_efek_kelompok_per_waktu")
ew <- emmeans(aov2, ~ waktu | kelompok, model = "multivariate")
tabel(summary(contrast(ew, "trt.vs.ctrl", ref = 1, adjust = "holm"),
infer = c(TRUE, TRUE)), "08_perubahan_vs_baseline")
ek <- emmeans(aov2, ~ kelompok | waktu, model = "multivariate")
tabel(summary(pairs(ek, adjust = "tukey"), infer = c(TRUE, TRUE)), "08_kelompok_per_waktu")
narasi("Holm untuk perubahan vs baseline berlaku pada dua perbandingan dalam masing-masing kelompok. Tukey berlaku pada tiga pasangan kelompok dalam masing-masing waktu. Keduanya bukan koreksi gabungan seluruh tabel. Efek sederhana di atas merupakan analisis lanjutan/eksploratif, terutama bila interaksi keseluruhan tidak signifikan.")
})
# ===================== 9. KONTRAS INTERAKSI =========================
bagian("9. Kontras perubahan baseline–akhir dan tren interaksi", function() {
if (is.null(aov2)) stop("Model Mixed Design ANOVA belum tersedia.")
e <- emmeans(aov2, ~ waktu * kelompok, model = "multivariate")
perubahan <- contrast(e, interaction = list(
waktu = list("M12-M0" = c(-1, 0, 1)), kelompok = "pairwise"), adjust = "holm")
tabel(summary(perubahan, infer = c(TRUE, TRUE)), "09_selisih_perubahan")
lin <- poly(minggu_nilai, degree = 2)[, 1]
tren <- contrast(e, interaction = list(
waktu = list(linear_minggu = lin), kelompok = "pairwise"), adjust = "holm")
tabel(summary(tren, infer = c(TRUE, TRUE)), "09_tren_linear_interaksi")
narasi("Kontras pertama membandingkan (TDS M12 − TDS M0) antarkelompok. Perubahan negatif berarti penurunan TDS. Estimate bertanda positif untuk Control − intervensi berarti perubahan kelompok intervensi lebih negatif (penurunannya lebih besar); interpretasi didasarkan pada CI dan nilai p Holm. Perbedaan signifikansi dalam masing-masing kelompok tidak dengan sendirinya membuktikan perbedaan efek antarkelompok.")
})
# ========================== 10. LMM ================================
lmm <- NULL; lmm_ri <- NULL; lmm_rs <- NULL
bagian("10. LMM dengan semua pengukuran tersedia", function() {
# REML cocok untuk estimasi dan KR. Perbandingan AIC di sini mempunyai fixed effects sama.
lmm_ri <<- lmerTest::lmer(tds ~ kelompok * waktu + (1 | id), data = long_obs,
REML = TRUE, control = lmerControl(optimizer = "bobyqa"))
lmm_rs <<- tryCatch(lmerTest::lmer(tds ~ kelompok * waktu + (1 + minggu | id), data = long_obs,
REML = TRUE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 200000))),
error = function(e) {
narasi(paste("Random slope gagal diestimasi; memakai random intercept:", conditionMessage(e)))
NULL
})
info_model <- function(m, nama) data.frame(model = nama, n_observasi = nobs(m),
AIC = AIC(m), BIC = BIC(m), singular = lme4::isSingular(m),
konvergensi = paste(m@optinfo$conv$lme4$messages, collapse = "; "))
info_rs <- if (!is.null(lmm_rs)) info_model(lmm_rs, "random_intercept_slope") else NULL
tabel(bind_rows(info_model(lmm_ri, "random_intercept"), info_rs), "10_perbandingan_model")
layak <- !is.null(lmm_rs) && !isSingular(lmm_rs) && is.null(lmm_rs@optinfo$conv$lme4$messages)
pakai_rs <- layak && AIC(lmm_rs) < AIC(lmm_ri)
lmm <<- if (pakai_rs) lmm_rs else lmm_ri
narasi(paste("Model yang dipakai:", if (pakai_rs) "random intercept + slope minggu." else "random intercept."))
narasi("Pemilihan ini memakai AIC dan diagnostik estimasi, bukan uji signifikansi slope acak. AIC merupakan pembanding relatif, bukan bukti bahwa asumsi model terpenuhi. Uji likelihood ratio standar untuk variance component memiliki persoalan batas parameter; tidak dipakai sebagai bukti formal di sini.")
output(summary(lmm)); output(VarCorr(lmm))
hasil <- anova(lmm, type = 3, ddf = "Kenward-Roger")
tabel(as.data.frame(hasil) |> tibble::rownames_to_column("efek"), "10_lmm_kr")
ringkas_interaksi(hasil, "Pr(>F)", "LMM seluruh pengukuran")
output(performance::icc(lmm_ri))
narasi("ICC di atas milik model random intercept. ICC pada model dengan slope acak bergantung pada waktu dan tidak boleh disamakan dengan satu persentase konstan. Faktor waktu pada fixed effects tetap kategorik; minggu numerik hanya digunakan pada random slope.")
e <- emmeans(lmm, ~ waktu * kelompok, lmer.df = "kenward-roger")
tabel(summary(contrast(e, interaction = list(waktu = list("M12-M0" = c(-1, 0, 1)),
kelompok = "pairwise"), adjust = "holm"), infer = c(TRUE, TRUE)), "10_lmm_selisih_perubahan")
tabel(summary(emmeans(lmm, ~ kelompok | waktu, lmer.df = "kenward-roger")), "10_lmm_rerata")
narasi("LMM mencakup kelompok, waktu, dan interaksi kelompok × waktu tanpa penyesuaian usia, BMI, status kehamilan, atau status puasa. Estimasi yang diperoleh merupakan hasil model tanpa pengendalian kovariat tersebut.")
})
bagian("11. Diagnostik LMM", function() {
if (is.null(lmm)) stop("LMM belum tersedia.")
diag <- data.frame(fitted = fitted(lmm), residual = resid(lmm), kelompok = long_obs$kelompok,
waktu = long_obs$waktu)
gambar(ggplot(diag, aes(sample = residual)) + stat_qq() + stat_qq_line() + theme_bw() +
labs(title = "Q–Q residual LMM"), "11_qq_residual")
gambar(ggplot(diag, aes(fitted, residual, colour = kelompok)) + geom_point(alpha = .35) +
geom_hline(yintercept = 0, linetype = 2) + theme_bw() +
labs(x = "Prediksi", y = "Residual", title = "Residual terhadap prediksi"), "11_residual_prediksi")
gambar(ggplot(diag, aes(waktu, residual, fill = kelompok)) + geom_boxplot() + theme_bw() +
labs(title = "Residual menurut kelompok dan waktu"), "11_residual_sel")
re <- ranef(lmm)$id
for (nm in names(re)) gambar(ggplot(data.frame(nilai = re[[nm]]), aes(sample = nilai)) +
stat_qq() + stat_qq_line() + theme_bw() + labs(title = paste("Q–Q efek acak", nm)),
paste0("11_qq_acak_", make.names(nm)))
output(performance::check_singularity(lmm))
output(performance::check_convergence(lmm))
narasi("Diagnostik mencakup bentuk distribusi residual, observasi ekstrem, dan variasi residual antarsel. Model lmer mengasumsikan residual homogen dan independen bersyarat pada efek acak. Model tidak mensyaratkan sferisitas, tetapi tetap memerlukan evaluasi asumsi residual dan efek acak.")
})
bagian("11b. Visualisasi diagnostik LMM tiga panel", function() {
if (is.null(lmm)) stop("LMM belum tersedia.")
# Gambar base-R tiga panel: Q-Q residual, Q-Q intercept acak, residual vs prediksi.
plot_tiga <- function() {
par_lama <- par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
on.exit(par(par_lama), add = TRUE)
qqnorm(resid(lmm), main = "Q-Q residual"); qqline(resid(lmm))
qqnorm(ranef(lmm)$id[, 1], main = "Q-Q intercept acak")
qqline(ranef(lmm)$id[, 1])
plot(fitted(lmm), resid(lmm), xlab = "Nilai prediksi", ylab = "Residual",
main = "Residual vs prediksi", pch = 16, col = grDevices::adjustcolor("black", .3))
abline(h = 0, lty = 2)
}
path <- file.path(folder, "11_diagnostik_lmm_tiga_panel.png")
grDevices::png(path, width = 2400, height = 800, res = 180)
tryCatch(plot_tiga(), finally = grDevices::dev.off())
if (TAMPILKAN_GRAFIK_RSTUDIO) plot_tiga()
tambah(paste0("<figure><img src='", base64enc::dataURI(file = path, mime = "image/png"),
"' alt='Diagnostik LMM tiga panel'><figcaption>Diagnostik LMM tiga panel</figcaption></figure>"))
message("Grafik tersimpan: ", path)
})
bagian("11c. Grafik rerata prediksi LMM", function() {
if (is.null(lmm)) stop("LMM belum tersedia.")
pred <- as.data.frame(summary(emmeans(lmm, ~ kelompok | waktu,
lmer.df = "kenward-roger")))
pred$minggu <- minggu_nilai[as.integer(pred$waktu)]
gambar(ggplot(pred, aes(minggu, emmean, colour = kelompok, group = kelompok)) +
geom_line(linewidth = 1) + geom_point(size = 2.5) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .3) +
scale_x_continuous(breaks = minggu_nilai) + theme_bw(base_size = 12) +
theme(legend.position = "bottom") +
labs(title = "Rerata prediksi LMM dan 95% CI", x = "Minggu", y = "TDS (mmHg)",
colour = "Kelompok"), "11_profil_prediksi_lmm")
narasi("Grafik LMM memakai seluruh pengukuran tersedia. Error bar merupakan CI rerata prediksi per sel, bukan CI selisih antar kelompok.")
})
bagian("12. LMM pada peserta lengkap: pembanding ANOVA", function() {
if (is.null(lmm)) stop("LMM belum tersedia.")
mcc <- update(lmm, data = long_cc)
output(summary(mcc))
tab <- anova(mcc, type = 3, ddf = "Kenward-Roger")
tabel(as.data.frame(tab) |> tibble::rownames_to_column("efek"), "12_lmm_complete_case")
ringkas_interaksi(tab, "Pr(>F)", "LMM peserta lengkap")
narasi("ANOVA dan LMM ini memakai peserta yang sama. Perbedaan dengan LMM semua pengukuran dapat berasal dari tambahan data dan asumsi kovarians. Tabel bukan pengujian formal atas selisih nilai p antar model.")
if (SIMULASI_MISSING_TAMBAHAN) {
set.seed(2026)
dm <- long_cc
pos <- which(dm$waktu != "M0")
dm$tds[sample(pos, min(30L, length(pos)))] <- NA_real_
ms <- update(lmm, data = dm, na.action = na.omit)
output(anova(ms, type = 3, ddf = "Kenward-Roger"))
narasi("Output tambahan ini SIMULASI kehilangan pengukuran, bukan hasil data asli. Seed = 2026; maksimal 30 pengukuran pasca-baseline dihilangkan dari peserta lengkap.")
} else narasi("Simulasi missing tambahan dimatikan karena data jurnal sudah memiliki missing asli.")
})
bagian("13. Informasi sesi dan batas interpretasi", function() {
narasi("Taraf signifikansi ditetapkan pada α = 0,05. Interpretasi hasil mempertimbangkan nilai p, ukuran efek, interval kepercayaan, dan arah perubahan. Nilai p yang ditampilkan sebagai 0 akibat pembulatan dilaporkan sebagai p < 0,001.")
narasi("Analisis menggunakan tiga waktu pengukuran asli tanpa penambahan pengukuran atau rekonstruksi data. Tren kubik tidak diestimasi karena hanya tersedia tiga waktu pengukuran. Kelayakan model ANOVA dan LMM dievaluasi berdasarkan diagnostik dan asumsi masing-masing model.")
output(sessionInfo())
})
# ======================= 14. MENYIMPAN LAPORAN ======================
tambah("<h2>Status setiap bagian</h2>")
tabel(status, "status_analisis")
tambah("<details><summary>Syntax yang digunakan</summary>")
file_script <- sumber[file.exists(sumber) & grepl("\\.R$", sumber, ignore.case = TRUE)]
if (length(file_script)) tambah(paste0("<pre>", escape(readLines(tail(file_script, 1), warn = FALSE)), "</pre>"))
tambah("</details>")
css <- "body{font:16px/1.6 Arial,sans-serif;max-width:1100px;margin:40px auto;padding:0 24px;color:#18332e}h1,h2{color:#0b6b4f}h2{border-top:1px solid #d9e4e0;padding-top:24px}pre{background:#f0f5f3;padding:16px;overflow:auto;font-size:13px;color:#182822}.table{overflow:auto}table{border-collapse:collapse;font-size:13px;width:100%}td,th{border:1px solid #d9e4e0;padding:8px;text-align:left}th{background:#e7f2ee}img{max-width:100%}summary{cursor:pointer;font-weight:bold}"
selesai <- sum(status$status == "Selesai")
header <- paste0("<!doctype html><html lang='id'><head><meta charset='UTF-8'>",
"<meta name='viewport' content='width=device-width,initial-scale=1'>",
"<title>Analisis TDS jurnal — pengukuran berulang</title><style>", css, "</style></head><body>",
"<h1>Analisis Pengukuran Berulang: Tekanan Darah Sistolik</h1>",
"<table class='identitas'><tbody>",
"<tr><th>Nama</th><td>Muhammad Ridhan Nazmy</td></tr>",
"<tr><th>NIM</th><td>2611018008</td></tr>",
"<tr><th>Rpubs</th><td>r_nazmy</td></tr>",
"<tr><th>Mata Kuliah</th><td>Biostatistika Intermediete</td></tr>",
"<tr><th>Program Studi / Institusi</th><td>Magister Kesehatan Masyarakat Universitas Mulawarman</td></tr>",
"</tbody></table>",
"<p>Data James dkk. (2019), PLOS Medicine · 3 kelompok × 3 waktu</p>",
"<p>Dibuat: ", escape(format(Sys.time())), " · Bagian selesai: ", selesai, "/", nrow(status), "</p>")
path_html <- file.path(folder, "laporan_analisis_tds.html")
writeLines(c(header, isi, "</body></html>"), path_html, useBytes = TRUE)
saveRDS(list(data_wide = dat_wide, data_long = dat_long, aov_satu = aov1,
aov_campuran = aov2, lmm = lmm, status = status), file.path(folder, "objek_analisis.rds"))
writeLines(capture.output(sessionInfo()), file.path(folder, "sessionInfo.txt"))
message("\nLaporan tersedia: ", path_html,
"\nDurasi: ", round(as.numeric(difftime(Sys.time(), waktu_mulai, units = "secs")), 1), " detik.")
if (any(status$status == "Gagal")) message("Terdapat bagian analisis yang gagal; rinciannya tercatat dalam status_analisis.csv dan laporan HTML.")
if (BUKA_HTML_OTOMATIS && interactive()) try(utils::browseURL(path_html), silent = TRUE)
invisible(list(folder = folder, status = status))
}
# Eksekusi analisis.
hasil_analisis <- jalankan_analisis()
## Registered S3 method overwritten by 'lme4':
## method from
## na.action.merMod car
## Data: /Users/nazmy/Downloads/paket_analisis_tds/pmed.1002870.s007.xls
## Hasil: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds
##
## >>> 1. Sumber, desain, dan kelengkapan data
## kelompok waktu n_total n_tersedia n_hilang
## 1 Control M0 100 100 0
## 2 Control M5 100 93 7
## 3 Control M12 100 82 18
## 4 UNIMMAP M0 105 104 1
## 5 UNIMMAP M5 105 92 13
## 6 UNIMMAP M12 105 84 21
## 7 Drink powder M0 93 93 0
## 8 Drink powder M5 93 88 5
## 9 Drink powder M12 93 77 16
## kelompok n_total n_lengkap n_tidak_lengkap
## 1 Control 100 80 20
## 2 UNIMMAP 105 81 24
## 3 Drink powder 93 77 16
##
## >>> 2. Statistik deskriptif dan visualisasi
## kelompok waktu n mean sd median min max se
## 1 Control M0 100 115.2733 14.49512 114.0000 87.33334 155.6667 1.449512
## 2 Control M5 93 107.7061 12.28819 106.6667 81.33334 157.0000 1.274226
## 3 Control M12 82 107.1870 11.99189 105.6667 84.33334 144.3333 1.324283
## 4 UNIMMAP M0 104 114.3365 13.43051 113.8333 88.33334 173.3333 1.316970
## 5 UNIMMAP M5 92 109.2464 11.27824 108.5000 86.66666 160.6667 1.175838
## 6 UNIMMAP M12 84 107.4524 10.73418 106.6667 89.33334 139.0000 1.171195
## 7 Drink powder M0 93 115.7670 13.71666 113.6667 85.33334 163.6667 1.422352
## 8 Drink powder M5 88 108.4091 11.07558 107.3333 88.00000 142.6667 1.180661
## 9 Drink powder M12 77 105.6234 10.11590 104.3333 89.00000 138.0000 1.152814
## lower_ci upper_ci
## 1 112.3972 118.1495
## 2 105.1754 110.2368
## 3 104.5521 109.8219
## 4 111.7246 116.9484
## 5 106.9107 111.5820
## 6 105.1229 109.7818
## 7 112.9421 118.5919
## 8 106.0624 110.7558
## 9 103.3273 107.9194
## kelompok waktu n mean sd median min max se
## 1 Control M0 80 115.2458 13.12953 114.5000 88.33334 155.6667 1.467927
## 2 Control M5 80 107.3375 11.65834 106.6667 81.33334 157.0000 1.303443
## 3 Control M12 80 107.0708 12.07078 105.6667 84.33334 144.3333 1.349554
## 4 UNIMMAP M0 81 113.7901 13.86347 112.0000 88.33334 173.3333 1.540385
## 5 UNIMMAP M5 81 109.1646 11.31544 108.0000 86.66666 160.6667 1.257271
## 6 UNIMMAP M12 81 107.1811 10.74374 106.0000 89.33334 139.0000 1.193749
## 7 Drink powder M0 77 116.2078 14.16053 113.6667 85.33334 163.6667 1.613742
## 8 Drink powder M5 77 108.4026 11.41000 107.3333 88.00000 142.6667 1.300290
## 9 Drink powder M12 77 105.6234 10.11590 104.3333 89.00000 138.0000 1.152814
## lower_ci upper_ci
## 1 112.3240 118.1677
## 2 104.7431 109.9319
## 3 104.3846 109.7571
## 4 110.7247 116.8556
## 5 106.6626 111.6667
## 6 104.8054 109.5567
## 7 112.9937 119.4218
## 8 105.8128 110.9924
## 9 103.3273 107.9194
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/02_profil_tds.png
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/02_lintasan.png
## TDS_M0 TDS_M5 TDS_M12
## TDS_M0 172.385 106.420 106.566
## TDS_M5 106.420 135.917 103.717
## TDS_M12 106.566 103.717 145.704
## TDS_M0 TDS_M5 TDS_M12
## TDS_M0 1.000 0.695 0.672
## TDS_M5 0.695 1.000 0.737
## TDS_M12 0.672 0.737 1.000
## kelompok pasangan var_selisih
## 1 Control TDS_M0 - TDS_M5 95.46266
## 2 Control TDS_M0 - TDS_M12 104.95634
## 3 Control TDS_M5 - TDS_M12 74.18677
## TDS_M0 TDS_M5 TDS_M12
## TDS_M0 192.196 120.700 109.179
## TDS_M5 120.700 128.039 95.649
## TDS_M12 109.179 95.649 115.428
## TDS_M0 TDS_M5 TDS_M12
## TDS_M0 1.000 0.769 0.733
## TDS_M5 0.769 1.000 0.787
## TDS_M12 0.733 0.787 1.000
## kelompok pasangan var_selisih
## 1 UNIMMAP TDS_M0 - TDS_M5 78.83439
## 2 UNIMMAP TDS_M0 - TDS_M12 89.26609
## 3 UNIMMAP TDS_M5 - TDS_M12 52.16917
## TDS_M0 TDS_M5 TDS_M12
## TDS_M0 200.521 105.048 92.531
## TDS_M5 105.048 130.188 88.488
## TDS_M12 92.531 88.488 102.331
## TDS_M0 TDS_M5 TDS_M12
## TDS_M0 1.000 0.650 0.646
## TDS_M5 0.650 1.000 0.767
## TDS_M12 0.646 0.767 1.000
## kelompok pasangan var_selisih
## 1 Drink powder TDS_M0 - TDS_M5 120.61215
## 2 Drink powder TDS_M0 - TDS_M12 117.78994
## 3 Drink powder TDS_M5 - TDS_M12 55.54272
##
## >>> 3. Asumsi ANOVA pada peserta lengkap
## kelompok waktu id usia bmi_baseline lengkap tds minggu
## 1 Control M0 219 27.80561 25.05134 TRUE 155.6667 0
## 2 Control M0 251 44.75565 16.35249 TRUE 147.6667 0
## 3 Control M5 219 27.80561 25.05134 TRUE 157.0000 5
## 4 Control M12 219 27.80561 25.05134 TRUE 144.3333 12
## 5 UNIMMAP M0 17 41.75496 20.26891 TRUE 173.3333 0
## 6 UNIMMAP M5 17 41.75496 20.26891 TRUE 160.6667 5
## 7 UNIMMAP M5 78 26.74333 29.67046 TRUE 135.0000 5
## 8 UNIMMAP M5 100 20.23819 27.16084 TRUE 135.6667 5
## 9 UNIMMAP M12 17 41.75496 20.26891 TRUE 139.0000 12
## 10 Drink powder M0 139 39.36208 27.02350 TRUE 159.6667 0
## 11 Drink powder M0 178 43.75359 22.89886 TRUE 149.3333 0
## 12 Drink powder M0 193 42.73785 26.27128 TRUE 163.6667 0
## 13 Drink powder M5 193 42.73785 26.27128 TRUE 142.6667 5
## 14 Drink powder M12 193 42.73785 26.27128 TRUE 138.0000 12
## is.outlier is.extreme
## 1 TRUE FALSE
## 2 TRUE FALSE
## 3 TRUE FALSE
## 4 TRUE FALSE
## 5 TRUE FALSE
## 6 TRUE TRUE
## 7 TRUE FALSE
## 8 TRUE FALSE
## 9 TRUE FALSE
## 10 TRUE FALSE
## 11 TRUE FALSE
## 12 TRUE FALSE
## 13 TRUE FALSE
## 14 TRUE FALSE
## kelompok waktu variable statistic p
## 1 Control M0 tds 0.9812348 0.2895775831
## 2 Control M5 tds 0.9527453 0.0049658648
## 3 Control M12 tds 0.9775393 0.1714387091
## 4 UNIMMAP M0 tds 0.9381109 0.0007017987
## 5 UNIMMAP M5 tds 0.9281315 0.0002143989
## 6 UNIMMAP M12 tds 0.9735596 0.0918928831
## 7 Drink powder M0 tds 0.9467618 0.0028473693
## 8 Drink powder M5 tds 0.9749280 0.1316866543
## 9 Drink powder M12 tds 0.9625858 0.0230943713
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/03_qq_per_sel.png
## waktu df1 df2 statistic p
## 1 M0 2 235 0.02160737 0.9786263
## 2 M5 2 235 0.27905977 0.7567450
## 3 M12 2 235 1.25261352 0.2876578
## statistic p.value parameter
## 1 11.32126 0.5016055 12
## method
## 1 Box's M-test for Homogeneity of Covariance Matrices
##
## >>> 4. RM ANOVA satu kelompok
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/04_qq_kelompok_fokus.png
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2 152 47.329 1.05e-16 * 0.384
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.812 0.000408 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.842 1.68, 127.96 1.89e-14 * 0.859 1.72, 130.53 1.08e-14
## p[HF]<.05
## 1 *
##
## Anova Table (Type 3 tests)
##
## Response: tds
## Effect df MSE F ges pes p.value
## 1 waktu 1.68, 127.96 58.20 47.33 *** .124 .384 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 2799061 1 25464.5 76 8353.947 < 2.2e-16 ***
## waktu 4637 2 7446.6 152 47.328 < 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.81212 0.00040811
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.84184 1.892e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8587805 1.084174e-14
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99098 8353.9 1 76 < 2.2e-16 ***
## waktu 1 0.49163 36.3 2 75 9.592e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## efek num Df den Df MSE F pes ges Pr(>F)
## 1 waktu 1.683673 127.9591 58.19515 47.32852 0.3837598 0.1235027 1.892223e-14
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.38 | [0.28, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
##
## >>> 4b. Grafik RM ANOVA satu kelompok
## Bagian gagal: package ggbeeswarm is required.
##
## >>> 5. Post hoc satu kelompok dan tren waktu
## waktu emmean SE df lower.CL upper.CL
## M0 116.2078 1.613742 76 112.9938 119.4218
## M5 108.4026 1.300290 76 105.8128 110.9924
## M12 105.6234 1.152814 76 103.3273 107.9194
##
## Confidence level used: 0.95
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## M0 - M5 7.805195 1.2515556 76 4.741232 10.869157 6.236 <0.0001
## M0 - M12 10.584416 1.2368264 76 7.556512 13.612319 8.558 <0.0001
## M5 - M12 2.779221 0.8493139 76 0.699996 4.858446 3.272 0.0048
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## M5 - M0 -7.805195 1.251556 76 -10.66710 -4.943291 -6.236 <0.0001
## M12 - M0 -10.584416 1.236826 76 -13.41264 -7.756194 -8.558 <0.0001
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: holm method for 2 tests
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## linear_minggu -7.253370 0.8461635 76 -9.188273 -5.318468 -8.572 <0.0001
## kuadratik_minggu 2.759278 0.7459864 76 1.053449 4.465108 3.699 0.0004
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: holm method for 2 tests
##
## >>> 6. Pembanding nonparametrik satu kelompok
## .y. n statistic df p method
## 1 tds 77 42.18241 2 6.921591e-10 Friedman test
## .y. n effsize method magnitude
## 1 tds 77 0.2739118 Kendall W small
## waktu_1 waktu_2 n_pasangan V p p_bonferroni
## 1 M0 M5 77 2535.0 1.560334e-07 4.681003e-07
## 2 M0 M12 77 2808.5 3.267120e-11 9.801360e-11
## 3 M5 M12 77 2015.0 4.296053e-03 1.288816e-02
##
## >>> 7. Mixed Design ANOVA kelompok × waktu
## Anova Table (Type 3 tests)
##
## Response: tds
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 235 353.53 0.01 <.001 <.001 .993
## 2 waktu 1.83, 429.09 47.90 109.06 *** .084 .317 <.001
## 3 kelompok:waktu 3.65, 429.09 47.90 2.76 * .005 .023 .032
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 8635803 1 83079 235 24427.5864 < 2e-16 ***
## kelompok 5 2 83079 235 0.0072 0.99283
## waktu 9538 2 20552 470 109.0614 < 2e-16 ***
## kelompok:waktu 483 4 20552 470 2.7635 0.02714 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.90466 8.1088e-06
## kelompok:waktu 0.90466 8.1088e-06
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.91296 < 2e-16 ***
## kelompok:waktu 0.91296 0.03156 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.9197368 1.106992e-36
## kelompok:waktu 0.9197368 3.118931e-02
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99047 24427.6 1 235 < 2e-16 ***
## kelompok 2 0.00006 0.0 2 235 0.99283
## waktu 1 0.41833 84.1 2 234 < 2e-16 ***
## kelompok:waktu 2 0.04347 2.6 4 470 0.03493 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## efek num Df den Df MSE F pes
## 1 kelompok 2.000000 235.0000 353.52665 7.195235e-03 6.123229e-05
## 2 waktu 1.825922 429.0917 47.89593 1.090614e+02 3.169823e-01
## 3 kelompok:waktu 3.651844 429.0917 47.89593 2.763450e+00 2.297831e-02
## ges Pr(>F)
## 1 4.908945e-05 9.928308e-01
## 2 8.428043e-02 1.965421e-36
## 3 4.642522e-03 3.155838e-02
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 6.12e-05 | [0.00, 1.00]
## waktu | 0.32 | [0.26, 1.00]
## kelompok:waktu | 0.02 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
##
## >>> 7b. Visualisasi interaksi kelompok dan waktu
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/07_interaksi_anova.png
##
## >>> 8. Efek sederhana dan post hoc antarkelompok
## model term kelompok df1 df2 F.ratio p.value p_holm_antar_kelompok
## 1 waktu Control 2 235 30.189 2.143183e-12 4.286366e-12
## 3 waktu UNIMMAP 2 235 17.063 1.204340e-07 1.204340e-07
## 5 waktu Drink powder 2 235 41.855 2.828708e-16 8.486125e-16
## model term waktu df1 df2 F.ratio p.value p_holm_antar_waktu
## 1 kelompok M0 2 235 0.624 0.5368033 1
## 3 kelompok M5 2 235 0.516 0.5978404 1
## 5 kelompok M12 2 235 0.487 0.6152642 1
## kelompok = Control:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## M5 - M0 -7.908333 1.106432 235 -10.404285 -5.412382 -7.148 <0.0001
## M12 - M0 -8.175000 1.138889 235 -10.744169 -5.605831 -7.178 <0.0001
##
## kelompok = UNIMMAP:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## M5 - M0 -4.625514 1.099581 235 -7.106011 -2.145018 -4.207 <0.0001
## M12 - M0 -6.609054 1.131837 235 -9.162314 -4.055793 -5.839 <0.0001
##
## kelompok = Drink powder:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## M5 - M0 -7.805195 1.127780 235 -10.349304 -5.261085 -6.921 <0.0001
## M12 - M0 -10.584416 1.160863 235 -13.203155 -7.965676 -9.118 <0.0001
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: holm method for 2 tests
## waktu = M0:
## contrast estimate SE df lower.CL upper.CL t.ratio
## Control - UNIMMAP 1.4557095 2.162558 235 -3.645031 6.556450 0.673
## Control - Drink powder -0.9619589 2.190291 235 -6.128112 4.204194 -0.439
## UNIMMAP - Drink powder -2.4176684 2.183650 235 -7.568157 2.732820 -1.107
## p.value
## 0.7793
## 0.8992
## 0.5106
##
## waktu = M5:
## contrast estimate SE df lower.CL upper.CL t.ratio
## Control - UNIMMAP -1.8271095 1.806734 235 -6.088582 2.434363 -1.011
## Control - Drink powder -1.0650977 1.829903 235 -5.381220 3.251024 -0.582
## UNIMMAP - Drink powder 0.7620118 1.824355 235 -3.541024 5.065047 0.418
## p.value
## 0.5705
## 0.8299
## 0.9084
##
## waktu = M12:
## contrast estimate SE df lower.CL upper.CL t.ratio
## Control - UNIMMAP -0.1102370 1.736527 235 -4.206116 3.985642 -0.063
## Control - Drink powder 1.4474567 1.758797 235 -2.700949 5.595862 0.823
## UNIMMAP - Drink powder 1.5576937 1.753464 235 -2.578134 5.693521 0.888
## p.value
## 0.9978
## 0.6892
## 0.6482
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimates
## P value adjustment: tukey method for comparing a family of 3 estimates
##
## >>> 9. Kontras perubahan baseline–akhir dan tren interaksi
## waktu_custom kelompok_pairwise estimate SE df lower.CL upper.CL
## M12-M0 Control - UNIMMAP -1.565946 1.605653 235 -5.437562 2.305669
## M12-M0 Control - Drink powder 2.409416 1.626244 235 -1.511850 6.330681
## M12-M0 UNIMMAP - Drink powder 3.975362 1.621314 235 0.065986 7.884738
## t.ratio p.value
## -0.975 0.3304
## 1.482 0.2796
## 2.452 0.0448
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: holm method for 3 tests
## waktu_custom kelompok_pairwise estimate SE df lower.CL
## linear_minggu Control - UNIMMAP -0.9066969 1.111020 235 -3.585633
## linear_minggu Control - Drink powder 1.7981626 1.125268 235 -0.915129
## linear_minggu UNIMMAP - Drink powder 2.7048595 1.121856 235 -0.000205
## upper.CL t.ratio p.value
## 1.772239 -0.816 0.4153
## 4.511454 1.598 0.2228
## 5.409924 2.411 0.0500
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: holm method for 3 tests
##
## >>> 10. LMM dengan semua pengukuran tersedia
## model n_observasi AIC BIC singular konvergensi
## 1 random_intercept 813 5998.540 6050.249 FALSE
## 2 random_intercept_slope 813 5978.954 6040.064 FALSE
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: tds ~ kelompok * waktu + (1 + minggu | id)
## Data: long_obs
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 5953
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4519 -0.5398 -0.0067 0.4443 3.6768
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 140.6573 11.8599
## minggu 0.1775 0.4213 -0.69
## Residual 38.1156 6.1738
## Number of obs: 813, groups: id, 298
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 110.1369 0.6413 294.1474 171.751 < 2e-16 ***
## kelompok1 -0.2348 0.9043 293.6072 -0.260 0.7953
## kelompok2 0.4430 0.8955 295.6780 0.495 0.6212
## waktu1 5.0354 0.3368 427.6135 14.949 < 2e-16 ***
## waktu2 -1.6282 0.3108 281.5937 -5.239 3.17e-07 ***
## kelompok1:waktu1 0.3358 0.4739 425.9012 0.709 0.4790
## kelompok2:waktu1 -1.1386 0.4731 431.9881 -2.407 0.0165 *
## kelompok1:waktu2 -0.8889 0.4375 282.2977 -2.032 0.0431 *
## kelompok2:waktu2 0.6981 0.4378 281.3613 1.595 0.1119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) klmpk1 klmpk2 waktu1 waktu2 klm1:1 klm2:1 klm1:2
## kelompok1 -0.008
## kelompok2 -0.036 -0.477
## waktu1 0.179 0.002 -0.015
## waktu2 0.026 -0.003 0.006 -0.386
## klmpk1:wkt1 0.002 0.182 -0.082 -0.015 0.009
## klmpk2:wkt1 -0.014 -0.081 0.172 -0.020 -0.001 -0.482
## klmpk1:wkt2 -0.003 0.025 -0.015 0.009 -0.013 -0.384 0.189
## klmpk2:wkt2 0.006 -0.015 0.031 -0.001 -0.011 0.189 -0.391 -0.487
## Groups Name Std.Dev. Corr
## id (Intercept) 11.85990
## minggu 0.42127 -0.695
## Residual 6.17378
## efek Sum Sq Mean Sq NumDF DenDF F value
## 1 kelompok 9.354443 4.677221 2 293.8515 0.1227115
## 2 waktu 8524.950471 4262.475236 2 356.4165 111.6258226
## 3 kelompok:waktu 369.277688 92.319422 4 399.2845 2.4162004
## Pr(>F)
## 1 8.845641e-01
## 2 2.283297e-38
## 3 4.825819e-02
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.710
## Unadjusted ICC: 0.650
## waktu_custom kelompok_pairwise estimate SE df lower.CL
## M12-M0 Control - UNIMMAP -1.361829 1.507993 270.73 -4.994505
## M12-M0 Control - Drink powder 2.013855 1.542200 268.30 -1.701435
## M12-M0 UNIMMAP - Drink powder 3.375684 1.534155 270.16 -0.320066
## upper.CL t.ratio p.value
## 2.270848 -0.903 0.3855
## 5.729145 1.306 0.3855
## 7.071433 2.200 0.0859
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: holm method for 3 tests
## waktu = M0:
## kelompok emmean SE df lower.CL upper.CL
## Control 115.2733 1.337060 320.91 112.6428 117.9038
## UNIMMAP 114.4767 1.307271 322.60 111.9048 117.0486
## Drink powder 115.7670 1.386467 320.91 113.0393 118.4947
##
## waktu = M5:
## kelompok emmean SE df lower.CL upper.CL
## Control 107.3850 1.237272 413.41 104.9529 109.8172
## UNIMMAP 109.6499 1.223874 427.58 107.2443 112.0554
## Drink powder 108.4914 1.278896 409.89 105.9774 111.0054
##
## waktu = M12:
## kelompok emmean SE df lower.CL upper.CL
## Control 107.0481 1.158294 302.20 104.7687 109.3274
## UNIMMAP 107.6133 1.140203 303.55 105.3696 109.8570
## Drink powder 105.5279 1.198209 300.59 103.1700 107.8859
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## >>> 11. Diagnostik LMM
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/11_qq_residual.png
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/11_residual_prediksi.png
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/11_residual_sel.png
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/11_qq_acak_X.Intercept..png
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/11_qq_acak_minggu.png
## [1] FALSE
## [1] TRUE
## attr(,"gradient")
## [1] 6.333805e-07
##
## >>> 11b. Visualisasi diagnostik LMM tiga panel
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/11_diagnostik_lmm_tiga_panel.png
##
## >>> 11c. Grafik rerata prediksi LMM
## Grafik tersimpan: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/11_profil_prediksi_lmm.png
##
## >>> 12. LMM pada peserta lengkap: pembanding ANOVA
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: tds ~ kelompok * waktu + (1 + minggu | id)
## Data: long_cc
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 5189.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4598 -0.5478 0.0054 0.4724 3.7376
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 135.1888 11.6271
## minggu 0.1639 0.4048 -0.66
## Residual 37.7729 6.1460
## Number of obs: 714, groups: id, 238
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 110.00264 0.70382 235.00008 156.293 < 2e-16 ***
## kelompok1 -0.11792 0.99316 235.00008 -0.119 0.90559
## kelompok2 0.04263 0.99011 235.00008 0.043 0.96569
## waktu1 5.07861 0.35774 393.11529 14.196 < 2e-16 ***
## waktu2 -1.70107 0.32582 240.72289 -5.221 3.85e-07 ***
## kelompok1:waktu1 0.28250 0.50480 393.11529 0.560 0.57606
## kelompok2:waktu1 -1.33376 0.50325 393.11529 -2.650 0.00837 **
## kelompok1:waktu2 -0.84615 0.45977 240.72289 -1.840 0.06694 .
## kelompok2:waktu2 0.82041 0.45836 240.72289 1.790 0.07473 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) klmpk1 klmpk2 waktu1 waktu2 klm1:1 klm2:1 klm1:2
## kelompok1 -0.006
## kelompok2 -0.015 -0.489
## waktu1 0.206 -0.001 -0.003
## waktu2 0.027 0.000 0.000 -0.432
## klmpk1:wkt1 -0.001 0.206 -0.101 -0.006 0.003
## klmpk2:wkt1 -0.003 -0.101 0.206 -0.015 0.006 -0.489
## klmpk1:wkt2 0.000 0.027 -0.013 0.003 -0.006 -0.432 0.211
## klmpk2:wkt2 0.000 -0.013 0.027 0.006 -0.015 0.211 -0.432 -0.489
## efek Sum Sq Mean Sq NumDF DenDF F value
## 1 kelompok 0.5435703 0.2717852 2 235.0000 7.195235e-03
## 2 waktu 7651.0757931 3825.5378966 2 312.6321 1.010618e+02
## 3 kelompok:waktu 411.9696384 102.9924096 4 351.0272 2.718905e+00
## Pr(>F)
## 1 9.928308e-01
## 2 1.402839e-34
## 3 2.963512e-02
##
## >>> 13. Informasi sesi dan batas interpretasi
## 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] lmerTest_3.2-1 effectsize_1.0.3 car_3.1-5 carData_3.0-6
## [5] rstatix_1.1.0 emmeans_2.0.4 afex_1.5-1 lme4_2.0-6
## [9] Matrix_1.7-6 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
## [4] fastmap_1.2.0 bayestestR_0.19.0 rpart_4.1.27
## [7] digest_0.6.39 estimability_2.0.0 lifecycle_1.0.5
## [10] cluster_2.1.8.3 magrittr_2.0.5 compiler_4.6.1
## [13] rlang_1.3.0 Hmisc_5.3-0 sass_0.4.10
## [16] tools_4.6.1 yaml_2.3.12 data.table_1.18.6.1
## [19] knitr_1.52 ggsignif_0.6.4 labeling_0.4.3
## [22] htmlwidgets_1.6.4 plyr_1.8.9 RColorBrewer_1.1-3
## [25] abind_1.4-8 foreign_0.8-91 withr_3.0.3
## [28] purrr_1.2.2 numDeriv_2016.8-1.1 stats4_4.6.1
## [31] nnet_7.3-21 grid_4.6.1 datawizard_1.4.0
## [34] ggpubr_1.0.0 colorspace_2.1-3 scales_1.4.0
## [37] MASS_7.3-66 insight_1.5.4 cli_3.6.6
## [40] mvtnorm_1.4-2 rmarkdown_2.32 ragg_1.5.2
## [43] reformulas_0.4.4 generics_0.1.4 otel_0.2.0
## [46] rstudioapi_0.19.0 performance_0.18.2 reshape2_1.4.5
## [49] parameters_0.29.3 readxl_1.5.0.1 minqa_1.2.8
## [52] cachem_1.1.0 stringr_1.6.0 splines_4.6.1
## [55] parallel_4.6.1 cellranger_1.1.0 base64enc_0.1-6
## [58] vctrs_0.7.3 boot_1.3-32 jsonlite_2.0.0
## [61] pbkrtest_0.5.5 htmlTable_2.5.0 Formula_1.2-6
## [64] systemfonts_1.3.2 jquerylib_0.1.4 glue_1.8.1
## [67] nloptr_2.2.1 stringi_1.8.9 gtable_0.3.6
## [70] tibble_3.3.1 pillar_1.11.1 htmltools_0.5.9
## [73] R6_2.6.1 textshaping_1.0.5 Rdpack_2.6.6
## [76] evaluate_1.0.5 lattice_0.23-1 rbibutils_2.4.1
## [79] backports_1.5.1 broom_1.0.13 bslib_0.12.0
## [82] Rcpp_1.1.2 checkmate_2.3.4 gridExtra_2.3.1
## [85] nlme_3.1-171 xfun_0.61 pkgconfig_2.0.3
## bagian status
## 1 1. Sumber, desain, dan kelengkapan data Selesai
## 2 2. Statistik deskriptif dan visualisasi Selesai
## 3 3. Asumsi ANOVA pada peserta lengkap Selesai
## 4 4. RM ANOVA satu kelompok Selesai
## 5 4b. Grafik RM ANOVA satu kelompok Gagal
## 6 5. Post hoc satu kelompok dan tren waktu Selesai
## 7 6. Pembanding nonparametrik satu kelompok Selesai
## 8 7. Mixed Design ANOVA kelompok × waktu Selesai
## 9 7b. Visualisasi interaksi kelompok dan waktu Selesai
## 10 8. Efek sederhana dan post hoc antarkelompok Selesai
## 11 9. Kontras perubahan baseline–akhir dan tren interaksi Selesai
## 12 10. LMM dengan semua pengukuran tersedia Selesai
## 13 11. Diagnostik LMM Selesai
## 14 11b. Visualisasi diagnostik LMM tiga panel Selesai
## 15 11c. Grafik rerata prediksi LMM Selesai
## 16 12. LMM pada peserta lengkap: pembanding ANOVA Selesai
## 17 13. Informasi sesi dan batas interpretasi Selesai
## peringatan
## 1
## 2
## 3
## 4
## 5
## 6
## 7
## 8
## 9 Panel(s) show a mixed within-between-design.\nError bars do not allow comparisons across all means.\nSuppress error bars with: error = "none"
## 10
## 11
## 12
## 13
## 14
## 15
## 16
## 17
##
## Laporan tersedia: /Users/nazmy/Downloads/paket_analisis_tds/hasil_analisis_tds/laporan_analisis_tds.html
## Durasi: 6.1 detik.
## Terdapat bagian analisis yang gagal; rinciannya tercatat dalam status_analisis.csv dan laporan HTML.