# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R - DATA EXCEL 3 KELOMPOK x 4 PENGUKURAN
#
# NAMA : VINA RAHMADANI
# NIM : 2611018027
#
# Data : data_stunting_tb_3kelompok.xlsx, sheet "Data"
# (DATA SIMULASI untuk latihan, bukan data penelitian nyata)
# Desain : 75 balita stunting dalam 3 kelompok (Kontrol, PMT, PMT+Edukasi;
# n = 25 per kelompok), tinggi badan (TB, cm) diukur 4 kali (B0-B3)
# Kolom : id, kelompok, usia_bln, jk, TB_B0, TB_B1, TB_B2, TB_B3
#
# Isi:
# 0. Paket & pengaturan
# 1. Impor data Excel (format lebar -> panjang)
# 2. Eksplorasi data: statistik deskriptif, profile plot, spaghetti plot
# 3. Repeated Measure ANOVA satu arah (within-subject: Waktu)
# 3a. Uji asumsi: outlier, normalitas, sfierisitas (Mauchly)
# 3b. ANOVA + koreksi Greenhouse-Geisser / Huynh-Feldt
# 3c. Pendekatan multivariat (MANOVA) sebagai pembanding
# 3d. Post hoc berpasangan & kontras polinomial (tren)
# 3e. Alternatif nonparametrik: uji Friedman
# 4. Mixed Design ANOVA (between: Kelompok x within: Waktu)
# 4a. Uji asumsi: outlier, normalitas, Levene, Box's M, Mauchly
# 4b. ANOVA campuran + ukuran efek
# 4c. Analisis efek sederhana (simple effects) & post hoc
# 4d. Kontras interaksi (pertambahan TB antarkelompok)
# 5. Pembanding: Linear Mixed Model (LMM) + kovariat usia & jenis kelamin
# 6. Analisis pertambahan TB (B3 - B0) dengan ANCOVA
# 7. Menyimpan output
# =============================================================================
# 0. PAKET & PENGATURAN
paket <- c("dplyr", "tidyr", "ggplot2", "afex", "emmeans", "rstatix", "car",
"effectsize", "lme4", "lmerTest", "pbkrtest", "performance", "ggpubr",
"Hmisc", "readxl")
baru <- paket[!paket %in% rownames(installed.packages())]
if (length(baru) > 0) install.packages(baru, repos = "https://cloud.r-project.org")
suppressPackageStartupMessages({
library(dplyr) # manipulasi data (menyediakan pipe %>%)
library(tidyr) # format panjang <-> lebar
library(ggplot2) # grafik
library(afex) # ANOVA within/mixed (tipe III, koreksi GG/HF, MANOVA)
library(emmeans) # rerata marginal, efek sederhana, post hoc, kontras
library(rstatix) # uji asumsi yang ramah pipe
library(car) # leveneTest, Anova
library(effectsize) # eta kuadrat parsial, omega kuadrat
library(lme4) # linear mixed model
library(lmerTest) # uji F/t dengan derajat bebas Kenward-Roger
})
options(contrasts = c("contr.sum", "contr.poly")) # kontras jumlah-nol untuk SS tipe III
afex_options(emmeans_model = "multivariate") # post hoc tahan terhadap non-sfierisitas
theme_set(theme_bw(base_size = 12))
# Pembungkus: jalankan langkah tambahan; jika galat, tampilkan pesan lalu lanjut.
coba <- function(expr) {
hasil <- try(expr, silent = TRUE)
if (inherits(hasil, "try-error")) {
message(" [dilewati] ", conditionMessage(attr(hasil, "condition")))
} else {
print(hasil)
}
invisible(hasil)
}
# 1. IMPOR DATA EXCEL
# Letakkan file Excel di working directory (cek dengan getwd(); ubah dengan setwd("folder_anda"))
# atau isi file_data dengan path lengkap, mis. "C:/Users/nama/Documents/data_stunting_tb_3kelompok.xlsx"
calon_file <- c("data_stunting_tb_3kelompok.xlsx", "data_stunting_tb_3kelompok (1).xlsx")
file_data <- calon_file[file.exists(calon_file)][1]
if (is.na(file_data)) {
if (interactive()) {
message("File Excel tidak ditemukan di: ", getwd(), "\nSilakan pilih file .xlsx secara manual.")
file_data <- file.choose()
} else {
stop("File Excel tidak ditemukan di: ", getwd(),
"\nPindahkan file ke folder tersebut atau gunakan setwd().")
}
}
message("Membaca data dari: ", file_data)
## Membaca data dari: data_stunting_tb_3kelompok.xlsx
kel_lab <- c("Kontrol", "PMT", "PMT+Edukasi")
bln_vec <- c(0, 1, 2, 3) # ASUMSI: B0-B3 = bulan ke-0 s.d. 3; ubah sesuai desain
kol_tb <- paste0("TB_B", 0:3) # nama kolom TB di file
dat_wide <- as.data.frame(readxl::read_excel(file_data, sheet = "Data"))
names(dat_wide) <- trimws(names(dat_wide))
# Pemeriksaan struktur & data hilang
stopifnot(all(c("id", "kelompok", "usia_bln", "jk", kol_tb) %in% names(dat_wide)))
# Pastikan kolom numerik terbaca sebagai angka (jaga-jaga bila desimal memakai koma)
for (k in c("usia_bln", kol_tb)) {
if (!is.numeric(dat_wide[[k]])) dat_wide[[k]] <- as.numeric(gsub(",", ".", dat_wide[[k]]))
}
dat_wide$kelompok <- trimws(dat_wide$kelompok)
dat_wide$jk <- trimws(dat_wide$jk)
print(colSums(is.na(dat_wide)))
## id kelompok usia_bln jk TB_B0 TB_B1 TB_B2 TB_B3
## 0 0 0 0 0 0 0 0
dat_wide$kelompok <- factor(dat_wide$kelompok, levels = kel_lab)
stopifnot(!any(is.na(dat_wide$kelompok))) # berhenti jika ada label kelompok tak cocok
dat_wide$jk <- factor(dat_wide$jk, levels = c("L", "P"))
dat_wide$id <- factor(dat_wide$id)
print(table(dat_wide$kelompok))
##
## Kontrol PMT PMT+Edukasi
## 25 25 25
# Format panjang (satu baris = satu pengukuran), dibutuhkan afex/rstatix/lme4
dat_long <- dat_wide %>%
pivot_longer(all_of(kol_tb), names_to = "waktu", values_to = "tb") %>%
mutate(waktu = factor(waktu, levels = kol_tb, labels = paste0("B", 0:3)),
bulan = bln_vec[as.integer(waktu)])
print(head(dat_wide))
## id kelompok usia_bln jk TB_B0 TB_B1 TB_B2 TB_B3
## 1 S001 Kontrol 12 L 68.4 69.0 70.5 70.5
## 2 S002 Kontrol 13 P 66.9 69.5 69.6 70.3
## 3 S003 Kontrol 13 L 71.0 71.0 72.2 72.6
## 4 S004 Kontrol 14 P 70.0 70.7 71.3 71.5
## 5 S005 Kontrol 15 P 70.4 71.2 72.5 73.0
## 6 S006 Kontrol 16 L 73.4 75.3 75.6 76.1
print(head(dat_long))
## # A tibble: 6 × 7
## id kelompok usia_bln jk waktu tb bulan
## <fct> <fct> <dbl> <fct> <fct> <dbl> <dbl>
## 1 S001 Kontrol 12 L B0 68.4 0
## 2 S001 Kontrol 12 L B1 69 1
## 3 S001 Kontrol 12 L B2 70.5 2
## 4 S001 Kontrol 12 L B3 70.5 3
## 5 S002 Kontrol 13 P B0 66.9 0
## 6 S002 Kontrol 13 P B1 69.5 1
str(dat_long)
## tibble [300 × 7] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 75 levels "S001","S002",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok: Factor w/ 3 levels "Kontrol","PMT",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia_bln: num [1:300] 12 12 12 12 13 13 13 13 13 13 ...
## $ jk : Factor w/ 2 levels "L","P": 1 1 1 1 2 2 2 2 1 1 ...
## $ waktu : Factor w/ 4 levels "B0","B1","B2",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ tb : num [1:300] 68.4 69 70.5 70.5 66.9 69.5 69.6 70.3 71 71 ...
## $ bulan : num [1:300] 0 1 2 3 0 1 2 3 0 1 ...
# 2. EKSPLORASI DATA
desk <- dat_long %>%
group_by(kelompok, waktu) %>%
get_summary_stats(tb, type = "mean_sd")
print(desk)
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol B0 tb 25 76.9 5.48
## 2 Kontrol B1 tb 25 77.7 5.18
## 3 Kontrol B2 tb 25 78.5 5.30
## 4 Kontrol B3 tb 25 79.2 5.32
## 5 PMT B0 tb 25 77.2 5.56
## 6 PMT B1 tb 25 78.2 5.38
## 7 PMT B2 tb 25 79.1 5.48
## 8 PMT B3 tb 25 80.1 5.33
## 9 PMT+Edukasi B0 tb 25 77.8 5.50
## 10 PMT+Edukasi B1 tb 25 78.8 5.52
## 11 PMT+Edukasi B2 tb 25 79.9 5.40
## 12 PMT+Edukasi B3 tb 25 81.0 5.46
# Kesetaraan kelompok pada baseline (usia, jenis kelamin, TB awal)
print(dat_wide %>% group_by(kelompok) %>%
get_summary_stats(usia_bln, TB_B0, type = "mean_sd"))
## # A tibble: 6 × 5
## kelompok variable n mean sd
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol usia_bln 25 23.2 7.63
## 2 Kontrol TB_B0 25 76.9 5.48
## 3 PMT usia_bln 25 23.1 7.56
## 4 PMT TB_B0 25 77.2 5.56
## 5 PMT+Edukasi usia_bln 25 23.3 7.62
## 6 PMT+Edukasi TB_B0 25 77.8 5.50
print(table(dat_wide$kelompok, dat_wide$jk))
##
## L P
## Kontrol 13 12
## PMT 8 17
## PMT+Edukasi 14 11
print(anova_test(data = dat_wide, dv = usia_bln, between = kelompok))
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 kelompok 2 72 0.003 0.997 8.33e-05
print(anova_test(data = dat_wide, dv = TB_B0, between = kelompok))
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 kelompok 2 72 0.181 0.835 0.005
# Matriks kovarians & korelasi antarwaktu (seluruh subjek)
S <- cov(dat_wide[, kol_tb])
R <- cor(dat_wide[, kol_tb])
print(round(S, 1)); print(round(R, 3))
## TB_B0 TB_B1 TB_B2 TB_B3
## TB_B0 29.7 28.8 29.0 28.9
## TB_B1 28.8 28.2 28.3 28.2
## TB_B2 29.0 28.3 28.7 28.5
## TB_B3 28.9 28.2 28.5 28.6
## TB_B0 TB_B1 TB_B2 TB_B3
## TB_B0 1.000 0.996 0.995 0.990
## TB_B1 0.996 1.000 0.997 0.994
## TB_B2 0.995 0.997 1.000 0.996
## TB_B3 0.990 0.994 0.996 1.000
# Varians selisih antarpasangan waktu (inti asumsi sfierisitas)
pasangan <- combn(kol_tb, 2)
var_selisih <- apply(pasangan, 2, function(p) var(dat_wide[[p[1]]] - dat_wide[[p[2]]]))
names(var_selisih) <- apply(pasangan, 2, paste, collapse = " - ")
print(round(var_selisih, 3))
## TB_B0 - TB_B1 TB_B0 - TB_B2 TB_B0 - TB_B3 TB_B1 - TB_B2 TB_B1 - TB_B3
## 0.271 0.319 0.578 0.192 0.361
## TB_B2 - TB_B3
## 0.241
# Profile plot: rerata dan 95% CI per kelompok (mean_cl_normal memerlukan paket Hmisc)
p_profil <- ggplot(dat_long, aes(bulan, tb, colour = kelompok, group = kelompok)) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", size = 2.5) +
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = .1) +
scale_x_continuous(breaks = bln_vec, labels = paste0("B", 0:3)) +
labs(x = "Waktu pengukuran", y = "Tinggi badan (cm)", colour = "Kelompok",
title = "Profil rerata TB (dengan 95% CI)") +
theme(legend.position = "bottom")
print(p_profil)
# Spaghetti plot: lintasan tiap anak
p_spag <- ggplot(dat_long, aes(bulan, tb, group = id)) +
geom_line(alpha = .3) +
stat_summary(aes(group = kelompok), fun = mean, geom = "line",
colour = "firebrick", linewidth = 1.2) +
facet_wrap(~ kelompok) +
scale_x_continuous(breaks = bln_vec, labels = paste0("B", 0:3)) +
labs(x = "Waktu pengukuran", y = "TB (cm)", title = "Lintasan individu dan rerata kelompok")
print(p_spag)
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan: apakah TB berubah selama pengukuran pada kelompok PMT+Edukasi?
# (ganti dengan "PMT" atau "Kontrol" untuk kelompok lain)
d1 <- droplevels(dplyr::filter(dat_long, kelompok == "PMT+Edukasi"))
d1w <- dplyr::filter(dat_wide, kelompok == "PMT+Edukasi")
## 3a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per waktu (ekstrem = di luar Q1-3IQR / Q3+3IQR)
print(d1 %>% group_by(waktu) %>% identify_outliers(tb))
## [1] waktu id kelompok usia_bln jk tb bulan
## [8] is.outlier is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas per waktu (Shapiro-Wilk) dan Q-Q plot
print(d1 %>% group_by(waktu) %>% shapiro_test(tb))
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 B0 tb 0.951 0.262
## 2 B1 tb 0.951 0.259
## 3 B2 tb 0.958 0.370
## 4 B3 tb 0.952 0.285
print(ggpubr::ggqqplot(d1, "tb", facet.by = "waktu"))
# (iii) Sfierisitas: Mauchly (dilaporkan otomatis oleh anova_test & afex)
aov1_rs <- anova_test(data = d1, dv = tb, wid = id, within = waktu,
effect.size = "pes")
print(aov1_rs) # ANOVA, Mauchly, koreksi GG & HF
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 72 353.278 5.61e-43 * 0.936
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.481 0.005 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.656 1.97, 47.26 5.05e-29 * 0.715 2.14, 51.45 2.18e-31
## p[HF]<.05
## 1 *
print(get_anova_table(aov1_rs, correction = "auto")) # GG dipakai jika Mauchly p < .05
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 1.97 47.26 353.278 5.05e-29 * 0.936
## 3b. ANOVA dengan afex (sumber utama laporan) -------------------------------
aov1 <- aov_ez(id = "id", dv = "tb", data = d1, within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
print(aov1)
## Anova Table (Type 3 tests)
##
## Response: tb
## Effect df MSE F ges pes p.value
## 1 waktu 1.97, 47.26 0.20 353.28 *** .047 .936 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
print(summary(aov1)) # univariat tanpa koreksi + Mauchly + epsilon GG & HF
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 629976 1 2865.36 24 5276.62 < 2.2e-16 ***
## waktu 140 3 9.53 72 353.28 < 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.48114 0.0053176
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.65638 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.7145845 2.175399e-31
# Ukuran efek tambahan
coba(eta_squared(aov1, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.94 | [0.91, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
coba(omega_squared(aov1, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## -------------------------------------------
## waktu | 0.04 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
## 3c. Pendekatan multivariat (tidak memerlukan sfierisitas) ------------------
coba(aov1$Anova) # Pillai, Wilks, Hotelling-Lawley, Roy
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99547 5276.6 1 24 < 2.2e-16 ***
## waktu 1 0.95850 169.4 3 22 2.383e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 3d. Post hoc & kontras tren ------------------------------------------------
em1 <- emmeans(aov1, ~ waktu)
print(em1)
## waktu emmean SE df lower.CL upper.CL
## B0 77.8 1.10 24 75.6 80.1
## B1 78.8 1.10 24 76.5 81.1
## B2 79.9 1.08 24 77.7 82.1
## B3 81.0 1.09 24 78.7 83.2
##
## Confidence level used: 0.95
print(pairs(em1, adjust = "bonferroni")) # semua pasangan waktu (6)
## contrast estimate SE df t.ratio p.value
## B0 - B1 -0.944 0.0887 24 -10.642 <0.0001
## B0 - B2 -2.052 0.1120 24 -18.400 <0.0001
## B0 - B3 -3.160 0.1420 24 -22.270 <0.0001
## B1 - B2 -1.108 0.0816 24 -13.573 <0.0001
## B1 - B3 -2.216 0.1000 24 -22.062 <0.0001
## B2 - B3 -1.108 0.0798 24 -13.889 <0.0001
##
## P value adjustment: bonferroni method for 6 tests
print(contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm")) # tiap waktu vs B0
## contrast estimate SE df t.ratio p.value
## B1 - B0 0.944 0.0887 24 10.642 <0.0001
## B2 - B0 2.052 0.1120 24 18.400 <0.0001
## B3 - B0 3.160 0.1420 24 22.270 <0.0001
##
## P value adjustment: holm method for 3 tests
print(contrast(em1, "poly")) # tren linear, kuadratik, kubik
## contrast estimate SE df t.ratio p.value
## linear 10.588 0.4610 24 22.955 <0.0001
## quadratic 0.164 0.0998 24 1.643 0.1134
## cubic -0.164 0.2350 24 -0.698 0.4920
## 3e. Alternatif nonparametrik ------------------------------------------------
print(friedman_test(d1, tb ~ waktu | id))
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 tb 25 73.8 3 6.40e-16 Friedman test
print(friedman_effsize(d1, tb ~ waktu | id)) # Kendall's W
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 tb 25 0.985 Kendall W large
print(d1 %>% wilcox_test(tb ~ waktu, paired = TRUE, p.adjust.method = "bonferroni"))
## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 tb B0 B1 25 25 1 0.000000119 7.15e-7 ****
## 2 tb B0 B2 25 25 0 0.0000000596 3.58e-7 ****
## 3 tb B0 B3 25 25 0 0.0000000596 3.58e-7 ****
## 4 tb B1 B2 25 25 0 0.0000000596 3.58e-7 ****
## 5 tb B1 B3 25 25 0 0.0000000596 3.58e-7 ****
## 6 tb B2 B3 25 25 0 0.0000000596 3.58e-7 ****
# (Opsional) ANOVA robust berbasis trimmed mean -- paket WRS2
if (requireNamespace("WRS2", quietly = TRUE)) {
coba(WRS2::rmanova(d1$tb, d1$waktu, d1$id, tr = 0.2))
}
## Call:
## WRS2::rmanova(y = d1$tb, groups = d1$waktu, blocks = d1$id, tr = 0.2)
##
## Test statistic: F = 121.4877
## Degrees of freedom 1: 1.65
## Degrees of freedom 2: 23.06
## p-value: 0
# 4. MIXED DESIGN ANOVA (Kelompok [between] x Waktu [within])
# Pertanyaan: apakah pola pertambahan TB berbeda antarkelompok intervensi?
## 4a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per sel
print(dat_long %>% group_by(kelompok, waktu) %>% identify_outliers(tb))
## [1] kelompok waktu id usia_bln jk tb bulan
## [8] is.outlier is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas per sel (3 x 4 = 12 sel) dan Q-Q plot
print(dat_long %>% group_by(kelompok, waktu) %>% shapiro_test(tb))
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol B0 tb 0.950 0.252
## 2 Kontrol B1 tb 0.939 0.144
## 3 Kontrol B2 tb 0.937 0.126
## 4 Kontrol B3 tb 0.946 0.205
## 5 PMT B0 tb 0.967 0.559
## 6 PMT B1 tb 0.965 0.525
## 7 PMT B2 tb 0.964 0.506
## 8 PMT B3 tb 0.969 0.619
## 9 PMT+Edukasi B0 tb 0.951 0.262
## 10 PMT+Edukasi B1 tb 0.951 0.259
## 11 PMT+Edukasi B2 tb 0.958 0.370
## 12 PMT+Edukasi B3 tb 0.952 0.285
print(ggpubr::ggqqplot(dat_long, "tb", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok))
# (iii) Homogenitas varians antarkelompok pada TIAP waktu (Levene)
print(dat_long %>% group_by(waktu) %>% levene_test(tb ~ kelompok))
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 B0 2 72 0.00155 0.998
## 2 B1 2 72 0.0395 0.961
## 3 B2 2 72 0.00856 0.991
## 4 B3 2 72 0.00783 0.992
# (iv) Homogenitas matriks kovarians antarkelompok (Box's M; uji pada alpha = .001)
coba(box_m(dat_wide[, kol_tb], dat_wide$kelompok))
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 16.6 0.678 20 Box's M-test for Homogeneity of Covariance Matric…
# (v) Sfierisitas: Mauchly (dari summary model di bawah)
## 4b. ANOVA campuran ---------------------------------------------------------
aov2 <- aov_ez(id = "id", dv = "tb", data = dat_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
print(aov2)
## Anova Table (Type 3 tests)
##
## Response: tb
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 72 116.72 0.37 .010 .010 .689
## 2 waktu 2.49, 179.01 0.17 752.58 *** .036 .913 <.001
## 3 kelompok:waktu 4.97, 179.01 0.17 7.21 *** <.001 .167 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
print(summary(aov2)) # Mauchly, epsilon GG/HF, p terkoreksi
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1857934 1 8403.8 72 15917.8949 < 2.2e-16 ***
## kelompok 87 2 8403.8 72 0.3748 0.6888
## waktu 316 3 30.2 216 752.5845 < 2.2e-16 ***
## kelompok:waktu 6 6 30.2 216 7.2059 5.056e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.74759 0.0009762
## kelompok:waktu 0.74759 0.0009762
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.82876 < 2.2e-16 ***
## kelompok:waktu 0.82876 3.775e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8607309 1.493665e-98
## kelompok:waktu 0.8607309 2.591159e-06
coba(aov2$Anova) # uji multivariat untuk efek within & interaksi
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99550 15917.9 1 72 < 2.2e-16 ***
## kelompok 2 0.01030 0.4 2 72 0.6887510
## waktu 1 0.95270 470.0 3 70 < 2.2e-16 ***
## kelompok:waktu 2 0.30488 4.3 6 142 0.0005658 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix (hasil identik; format ringkas untuk laporan)
aov2_rs <- anova_test(data = dat_long, dv = tb, wid = id,
between = kelompok, within = waktu, effect.size = "pes",
type = 3)
print(get_anova_table(aov2_rs, correction = "GG"))
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 kelompok 2.00 72.00 0.375 6.89e-01 0.010
## 2 waktu 2.49 179.01 752.584 5.37e-95 * 0.913
## 3 kelompok:waktu 4.97 179.01 7.206 3.77e-06 * 0.167
# Ukuran efek
coba(eta_squared(aov2, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.01 | [0.00, 1.00]
## waktu | 0.91 | [0.90, 1.00]
## kelompok:waktu | 0.17 | [0.08, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
coba(omega_squared(aov2, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## ------------------------------------------------
## kelompok | 0.00 | [0.00, 1.00]
## waktu | 0.04 | [0.00, 1.00]
## kelompok:waktu | 6.09e-04 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Plot interaksi dari model (data_geom = geom_point agar tidak butuh ggbeeswarm)
coba(afex_plot(aov2, x = "waktu", trace = "kelompok", error = "within",
mapping = c("colour", "shape", "linetype"),
data_geom = ggplot2::geom_point, data_alpha = .25) +
labs(y = "TB (cm)", x = "Waktu"))
## Warning: Panel(s) show a mixed within-between-design.
## Error bars do not allow comparisons across all means.
## Suppress error bars with: error = "none"
## 4c. Efek sederhana (bila interaksi signifikan) -----------------------------
em2 <- emmeans(aov2, ~ waktu | kelompok)
# Efek WAKTU di dalam tiap kelompok (uji F gabungan per kelompok)
coba(joint_tests(aov2, by = "kelompok"))
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## kelompok = Kontrol:
## model term df1 df2 F.ratio p.value
## waktu 3 72 111.449 <0.0001
##
## kelompok = PMT:
## model term df1 df2 F.ratio p.value
## waktu 3 72 165.173 <0.0001
##
## kelompok = PMT+Edukasi:
## model term df1 df2 F.ratio p.value
## waktu 3 72 216.700 <0.0001
# Efek KELOMPOK pada tiap waktu
coba(joint_tests(aov2, by = "waktu"))
## waktu = B0:
## model term df1 df2 F.ratio p.value
## kelompok 2 72 0.181 0.8346
##
## waktu = B1:
## model term df1 df2 F.ratio p.value
## kelompok 2 72 0.273 0.7619
##
## waktu = B2:
## model term df1 df2 F.ratio p.value
## kelompok 2 72 0.427 0.6540
##
## waktu = B3:
## model term df1 df2 F.ratio p.value
## kelompok 2 72 0.727 0.4869
# Post hoc: tiap waktu vs B0 di dalam tiap kelompok
print(contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm"))
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## B1 - B0 0.752 0.103 72 7.282 <0.0001
## B2 - B0 1.572 0.107 72 14.662 <0.0001
## B3 - B0 2.256 0.134 72 16.836 <0.0001
##
## kelompok = PMT:
## contrast estimate SE df t.ratio p.value
## B1 - B0 1.004 0.103 72 9.723 <0.0001
## B2 - B0 1.860 0.107 72 17.348 <0.0001
## B3 - B0 2.836 0.134 72 21.165 <0.0001
##
## kelompok = PMT+Edukasi:
## contrast estimate SE df t.ratio p.value
## B1 - B0 0.944 0.103 72 9.142 <0.0001
## B2 - B0 2.052 0.107 72 19.138 <0.0001
## B3 - B0 3.160 0.134 72 23.583 <0.0001
##
## P value adjustment: holm method for 3 tests
# Post hoc: perbandingan antarkelompok pada tiap waktu
em2b <- emmeans(aov2, ~ kelompok | waktu)
print(pairs(em2b, adjust = "tukey"))
## waktu = B0:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT -0.340 1.56 72 -0.218 0.9742
## Kontrol - (PMT+Edukasi) -0.928 1.56 72 -0.595 0.8233
## PMT - (PMT+Edukasi) -0.588 1.56 72 -0.377 0.9247
##
## waktu = B1:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT -0.592 1.52 72 -0.390 0.9196
## Kontrol - (PMT+Edukasi) -1.120 1.52 72 -0.738 0.7415
## PMT - (PMT+Edukasi) -0.528 1.52 72 -0.348 0.9354
##
## waktu = B2:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT -0.628 1.53 72 -0.411 0.9110
## Kontrol - (PMT+Edukasi) -1.408 1.53 72 -0.922 0.6279
## PMT - (PMT+Edukasi) -0.780 1.53 72 -0.511 0.8662
##
## waktu = B3:
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT -0.920 1.52 72 -0.606 0.8176
## Kontrol - (PMT+Edukasi) -1.832 1.52 72 -1.206 0.4537
## PMT - (PMT+Edukasi) -0.912 1.52 72 -0.600 0.8204
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## 4d. Kontras interaksi -------------------------------------------------------
# Apakah pertambahan TB (B3 - B0) berbeda antarkelompok?
# Langkah 1: pertambahan TB (B3 - B0) di tiap kelompok.
gain_k <- contrast(em2, list("B3-B0" = c(-1, 0, 0, 1)))
print(gain_k)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## B3-B0 2.26 0.134 72 16.836 <0.0001
##
## kelompok = PMT:
## contrast estimate SE df t.ratio p.value
## B3-B0 2.84 0.134 72 21.165 <0.0001
##
## kelompok = PMT+Edukasi:
## contrast estimate SE df t.ratio p.value
## B3-B0 3.16 0.134 72 23.583 <0.0001
print(confint(gain_k))
## kelompok = Kontrol:
## contrast estimate SE df lower.CL upper.CL
## B3-B0 2.26 0.134 72 1.99 2.52
##
## kelompok = PMT:
## contrast estimate SE df lower.CL upper.CL
## B3-B0 2.84 0.134 72 2.57 3.10
##
## kelompok = PMT+Edukasi:
## contrast estimate SE df lower.CL upper.CL
## B3-B0 3.16 0.134 72 2.89 3.43
##
## Confidence level used: 0.95
# Langkah 2: bandingkan pertambahan tersebut antarkelompok (koreksi Holm).
print(pairs(gain_k, by = NULL, adjust = "holm"))
## contrast estimate SE df t.ratio p.value
## (B3-B0 Kontrol) - (B3-B0 PMT) -0.580 0.189 72 -3.061 0.0062
## (B3-B0 Kontrol) - (B3-B0 PMT+Edukasi) -0.904 0.189 72 -4.770 <0.0001
## (B3-B0 PMT) - (B3-B0 PMT+Edukasi) -0.324 0.189 72 -1.710 0.0916
##
## P value adjustment: holm method for 3 tests
# Laju pertambahan TB (kemiringan linear, cm/bulan; ASUMSI jarak waktu sama = 1 bulan)
lin_k <- contrast(em2, list("slope_cm_per_bulan" = c(-.3, -.1, .1, .3)))
print(lin_k)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## slope_cm_per_bulan 0.759 0.0424 72 17.914 <0.0001
##
## kelompok = PMT:
## contrast estimate SE df t.ratio p.value
## slope_cm_per_bulan 0.936 0.0424 72 22.107 <0.0001
##
## kelompok = PMT+Edukasi:
## contrast estimate SE df t.ratio p.value
## slope_cm_per_bulan 1.059 0.0424 72 24.997 <0.0001
print(pairs(lin_k, by = NULL, adjust = "holm")) # apakah laju berbeda antarkelompok?
## contrast estimate SE
## slope_cm_per_bulan Kontrol - slope_cm_per_bulan PMT -0.178 0.0599
## slope_cm_per_bulan Kontrol - (slope_cm_per_bulan PMT+Edukasi) -0.300 0.0599
## slope_cm_per_bulan PMT - (slope_cm_per_bulan PMT+Edukasi) -0.122 0.0599
## df t.ratio p.value
## 72 -2.965 0.0082
## 72 -5.008 <0.0001
## 72 -2.043 0.0447
##
## P value adjustment: holm method for 3 tests
# Tren polinomial per kelompok (linear, kuadratik, kubik)
print(contrast(em2, "poly"))
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## linear 7.588 0.424 72 17.914 <0.0001
## quadratic -0.068 0.130 72 -0.523 0.6025
## cubic -0.204 0.268 72 -0.760 0.4497
##
## kelompok = PMT:
## contrast estimate SE df t.ratio p.value
## linear 9.364 0.424 72 22.107 <0.0001
## quadratic -0.028 0.130 72 -0.215 0.8301
## cubic 0.268 0.268 72 0.999 0.3213
##
## kelompok = PMT+Edukasi:
## contrast estimate SE df t.ratio p.value
## linear 10.588 0.424 72 24.997 <0.0001
## quadratic 0.164 0.130 72 1.261 0.2112
## cubic -0.164 0.268 72 -0.611 0.5431
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Tidak mensyaratkan sfierisitas dan menampung data hilang (MAR).
# lmm3 menambah usia (dipusatkan) dan jenis kelamin sebagai kovariat, karena
# TB sangat dipengaruhi usia dan kelompok bisa berbeda usia pada baseline.
dat_long$usia_c <- dat_long$usia_bln - mean(dat_wide$usia_bln) # usia dipusatkan
lmm1 <- lmer(tb ~ kelompok * waktu + (1 | id), data = dat_long, REML = TRUE)
lmm2 <- lmer(tb ~ kelompok * waktu + (1 + bulan | id), data = dat_long, REML = TRUE,
control = lmerControl(optimizer = "bobyqa"))
coba(anova(lmm1, lmm2, refit = FALSE)) # uji rasio kemungkinan struktur acak
## Data: dat_long
## Models:
## lmm1: tb ~ kelompok * waktu + (1 | id)
## lmm2: tb ~ kelompok * waktu + (1 + bulan | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 819.03 870.88 -395.52 791.03
## lmm2 16 802.86 862.12 -385.43 770.86 20.167 2 4.176e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coba(anova(lmm2, ddf = "Kenward-Roger")) # uji F tipe III efek tetap
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 0.073 0.037 2 72.00 0.3748 0.6887508
## waktu 137.868 45.956 3 153.23 467.1471 < 2.2e-16 ***
## kelompok:waktu 2.845 0.474 6 170.40 4.8134 0.0001457 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coba(performance::icc(lmm1)) # korelasi intrakelas
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.995
## Unadjusted ICC: 0.951
lmm3 <- lmer(tb ~ kelompok * waktu + usia_c + jk + (1 + bulan | id),
data = dat_long, REML = TRUE,
control = lmerControl(optimizer = "bobyqa"))
coba(anova(lmm3, ddf = "Kenward-Roger"))
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 1.315 0.658 2 70.258 6.7227 0.0021308 **
## waktu 137.868 45.956 3 153.227 467.1471 < 2.2e-16 ***
## usia_c 137.801 137.801 1 70.000 1408.5559 < 2.2e-16 ***
## jk 3.773 3.773 1 70.000 38.5687 3.319e-08 ***
## kelompok:waktu 2.845 0.474 6 170.399 4.8134 0.0001457 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(lmm3))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: tb ~ kelompok * waktu + usia_c + jk + (1 + bulan | id)
## Data: dat_long
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 558.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.83439 -0.46382 0.00244 0.50915 2.70801
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 1.28095 1.1318
## bulan 0.02529 0.1590 0.21
## Residual 0.09783 0.3128
## Number of obs: 300, groups: id, 75
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 78.75421 0.14067 69.56597 559.860 < 2e-16 ***
## kelompok1 -0.72162 0.19906 69.58189 -3.625 0.000546 ***
## kelompok2 0.27311 0.20257 70.01500 1.348 0.181935
## waktu1 -1.36967 0.04168 115.09906 -32.864 < 2e-16 ***
## waktu2 -0.46967 0.03260 185.09393 -14.408 < 2e-16 ***
## waktu3 0.45833 0.03260 185.09393 14.060 < 2e-16 ***
## usia_c 0.68687 0.01805 69.99998 38.063 < 2e-16 ***
## jk1 0.86812 0.13783 69.99998 6.298 2.31e-08 ***
## kelompok1:waktu1 0.22467 0.05894 115.09906 3.812 0.000223 ***
## kelompok2:waktu1 -0.05533 0.05894 115.09906 -0.939 0.349794
## kelompok1:waktu2 0.07667 0.04610 185.09393 1.663 0.097994 .
## kelompok2:waktu2 0.04867 0.04610 185.09393 1.056 0.292492
## kelompok1:waktu3 -0.03133 0.04610 185.09393 -0.680 0.497555
## kelompok2:waktu3 -0.02333 0.04610 185.09393 -0.506 0.613356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(summary(lmm3), correlation=TRUE) or
## vcov(summary(lmm3)) if you need it
# Diagnostik residual LMM (normalitas & homogenitas)
par(mfrow = c(1, 3))
qqnorm(resid(lmm2), main = "Q-Q residual"); qqline(resid(lmm2))
qqnorm(ranef(lmm2)$id[, 1], main = "Q-Q intersep acak"); qqline(ranef(lmm2)$id[, 1])
plot(fitted(lmm2), resid(lmm2), xlab = "Nilai prediksi", ylab = "Residual",
main = "Residual vs prediksi"); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# 6. PERTAMBAHAN TB (B3 - B0) DENGAN ANCOVA
# CATATAN: TB_B0 dan usia_bln berkorelasi sangat tinggi (r ~ 0,96 pada data ini), jadi jangan
# dimasukkan bersamaan dalam satu model (multikolinearitas). Dua model terpisah:
dat_wide$gain <- dat_wide$TB_B3 - dat_wide$TB_B0
print(dat_wide %>% group_by(kelompok) %>% get_summary_stats(gain, type = "mean_sd"))
## # A tibble: 3 × 5
## kelompok variable n mean sd
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol gain 25 2.26 0.638
## 2 PMT gain 25 2.84 0.66
## 3 PMT+Edukasi gain 25 3.16 0.709
print(cor(dat_wide$usia_bln, dat_wide$TB_B0))
## [1] 0.9631404
# Model A: koreksi TB baseline
m_a <- lm(gain ~ TB_B0 + jk + kelompok, data = dat_wide)
print(car::Anova(m_a, type = 3))
## Anova Table (Type III tests)
##
## Response: gain
## Sum Sq Df F value Pr(>F)
## (Intercept) 9.0833 1 21.3492 1.696e-05 ***
## TB_B0 1.9785 1 4.6503 0.03448 *
## jk 0.1617 1 0.3801 0.53957
## kelompok 11.0879 2 13.0305 1.546e-05 ***
## Residuals 29.7823 70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(pairs(emmeans(m_a, ~ kelompok), adjust = "tukey"))
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT -0.571 0.187 70 -3.046 0.0091
## Kontrol - (PMT+Edukasi) -0.937 0.185 70 -5.063 <0.0001
## PMT - (PMT+Edukasi) -0.366 0.188 70 -1.943 0.1343
##
## Results are averaged over the levels of: jk
## P value adjustment: tukey method for comparing a family of 3 estimates
# Model B: koreksi usia
m_b <- lm(gain ~ usia_bln + jk + kelompok, data = dat_wide)
print(car::Anova(m_b, type = 3))
## Anova Table (Type III tests)
##
## Response: gain
## Sum Sq Df F value Pr(>F)
## (Intercept) 74.396 1 175.6857 < 2.2e-16 ***
## usia_bln 2.118 1 5.0026 0.0285 *
## jk 0.400 1 0.9450 0.3343
## kelompok 10.555 2 12.4633 2.343e-05 ***
## Residuals 29.642 70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(pairs(emmeans(m_b, ~ kelompok), adjust = "tukey"))
## contrast estimate SE df t.ratio p.value
## Kontrol - PMT -0.549 0.187 70 -2.942 0.0121
## Kontrol - (PMT+Edukasi) -0.913 0.184 70 -4.956 <0.0001
## PMT - (PMT+Edukasi) -0.364 0.188 70 -1.937 0.1359
##
## Results are averaged over the levels of: jk
## P value adjustment: tukey method for comparing a family of 3 estimates
# Diagnostik residual model A
par(mfrow = c(1, 2))
qqnorm(resid(m_a), main = "Q-Q residual"); qqline(resid(m_a))
plot(fitted(m_a), resid(m_a), xlab = "Nilai prediksi", ylab = "Residual",
main = "Residual vs prediksi"); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# 7. SIMPAN OUTPUT & SESSION INFO
write.csv(desk, "Ringkasan_Deskriptif_TB.csv", row.names = FALSE)
write.csv(dat_long, "data_tb_long.csv", row.names = FALSE)
ggsave("Profil_TB.png", p_profil, width = 7, height = 4.5, dpi = 300)
ggsave("Spaghetti_TB.png", p_spag, width = 9, height = 4.5, dpi = 300)
message("Selesai. Output disimpan di: ", getwd())
## Selesai. Output disimpan di: /Users/cang/Downloads/untitled folder
print(sessionInfo())
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Ventura 13.3.1
##
## 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/Makassar
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lmerTest_3.2-1 effectsize_1.0.3 car_3.1-5 carData_3.0-6
## [5] rstatix_1.1.0 emmeans_2.0.4 afex_1.5-1 lme4_2.0-6
## [9] Matrix_1.7-5 ggplot2_4.0.3 tidyr_1.3.2 dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 farver_2.1.2 S7_0.2.2
## [4] fastmap_1.2.0 reshape_0.8.10 bayestestR_0.19.0
## [7] digest_0.6.39 rpart_4.1.27 estimability_2.0.0
## [10] lifecycle_1.0.5 cluster_2.1.8.2 magrittr_2.0.5
## [13] compiler_4.6.1 rlang_1.3.0 Hmisc_5.3-0
## [16] sass_0.4.10 tools_4.6.1 utf8_1.2.6
## [19] yaml_2.3.12 data.table_1.18.6.1 ggsignif_0.6.4
## [22] knitr_1.51 labeling_0.4.3 htmlwidgets_1.6.4
## [25] plyr_1.8.9 RColorBrewer_1.1-3 abind_1.4-8
## [28] withr_3.0.3 foreign_0.8-91 purrr_1.2.2
## [31] numDeriv_2016.8-1.1 nnet_7.3-20 grid_4.6.1
## [34] datawizard_1.4.0 ggpubr_1.0.0 colorspace_2.1-3
## [37] scales_1.4.0 MASS_7.3-65 insight_1.5.4
## [40] cli_3.6.6 mvtnorm_1.4-2 rmarkdown_2.31
## [43] reformulas_0.4.4 generics_0.1.4 performance_0.18.2
## [46] rstudioapi_0.19.0 reshape2_1.4.5 parameters_0.29.3
## [49] readxl_1.5.0 minqa_1.2.8 cachem_1.1.0
## [52] stringr_1.6.0 splines_4.6.1 parallel_4.6.1
## [55] WRS2_1.1-7 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 Formula_1.2-6 htmlTable_2.5.0
## [64] jquerylib_0.1.4 glue_1.8.1 nloptr_2.2.1
## [67] stringi_1.8.9 gtable_0.3.6 tibble_3.3.1
## [70] pillar_1.11.1 htmltools_0.5.9 R6_2.6.1
## [73] Rdpack_2.6.6 evaluate_1.0.5 lattice_0.22-9
## [76] rbibutils_2.4.1 backports_1.5.1 broom_1.0.13
## [79] bslib_0.12.0 Rcpp_1.1.2 gridExtra_2.3.1
## [82] nlme_3.1-169 checkmate_2.3.4 xfun_0.60
## [85] pkgconfig_2.0.3