1. Pendahuluan

1.1 Latar belakang dan tujuan

Pemberian Makanan Tambahan (PMT) merupakan salah satu intervensi untuk memperbaiki status gizi balita. Laporan ini membandingkan perubahan berat badan balita selama tiga bulan pada tiga kelompok:

Kontrol (tanpa PMT),
PMT Biskuit,
PMT Lokal.

Berat badan ditimbang empat kali (bulan ke-0, 1, 2, dan 3) pada balita yang sama, sehingga datanya berupa pengukuran berulang (repeated measures).

Tujuan analisis:

menilai apakah berat badan berubah selama tiga bulan (efek waktu);
menilai apakah pola perubahan berat badan berbeda antarkelompok (efek interaksi kelompok × waktu);
mengetahui kelompok mana yang berbeda dan pada waktu mana (uji lanjut);
memastikan hasil tetap konsisten setelah usia dan jenis kelamin dikontrol (linear mixed model, ANCOVA).

1.2 Desain dan variabel

Komponen	Keterangan
Subjek	90 balita (30 per kelompok)
Faktor antar-subjek (between)	Kelompok: Kontrol, PMT Biskuit, PMT Lokal
Faktor dalam-subjek (within)	Waktu: bulan 0, 1, 2, 3
Variabel terikat	Berat badan (`BB_kg`, kg)
Kovariat	Usia (`Usia_bulan`, 12–59 bulan) dan jenis kelamin
Uji utama	Mixed (two-way) ANOVA 3 × 4, dengan koreksi Greenhouse-Geisser
Uji pendukung	RM ANOVA satu arah, Friedman, ANCOVA, linear mixed model (LMM)

1.3 Hipotesis uji utama

H0 (interaksi): perubahan berat badan dari bulan 0 sampai bulan 3 sama pada ketiga kelompok.
H1 (interaksi): perubahan berat badan berbeda pada sedikitnya satu kelompok.

2. Ringkasan Temuan

Pertanyaan	Hasil utama
Apakah berat badan berubah seiring waktu?	Ya. Efek waktu signifikan, F(1,98; 172,38) = 803,10; p < 0,001; η² parsial = 0,902.
Apakah pola perubahan berbeda antarkelompok?	Ya. Interaksi kelompok × waktu signifikan, F(3,96; 172,38) = 70,39; p < 0,001; η² parsial = 0,618.
Berapa kenaikan berat badan bulan 3 − bulan 0?	Kontrol 0,34 kg; PMT Biskuit 0,95 kg; PMT Lokal 1,13 kg.
Apakah kenaikan antarkelompok berbeda?	Ya. Biskuit − Kontrol 0,61 kg; Lokal − Kontrol 0,79 kg; Lokal − Biskuit 0,18 kg (semua p ≤ 0,0014, Holm).
Apakah efek utama kelompok signifikan?	Tidak (p = 0,687). Berat badan awal ketiga kelompok mirip, sehingga perbedaan baru tampak pada laju kenaikan.
Apakah hasil tahan terhadap koreksi usia dan jenis kelamin?	Ya. Pada ANCOVA dan LMM, perbedaan antarkelompok tetap signifikan (p < 0,001).

3. Persiapan Paket dan Data

3.1 Paket

paket <- c("readxl", "dplyr", "tidyr", "ggplot2", "afex", "emmeans",
           "rstatix", "car", "effectsize", "lme4", "lmerTest", "pbkrtest",
           "performance", "ggpubr", "Hmisc", "knitr")
baru <- paket[!paket %in% rownames(installed.packages())]
if (length(baru) > 0) install.packages(baru)

suppressPackageStartupMessages({
  library(readxl)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(afex)
  library(emmeans)
  library(rstatix)
  library(car)
  library(effectsize)
  library(lme4)
  library(lmerTest)
  library(performance)
  library(ggpubr)
})

options(contrasts = c("contr.sum", "contr.poly"))
afex_options(emmeans_model = "multivariate")
theme_set(theme_bw(base_size = 12))

Pengaturan contr.sum dipakai agar jumlah kuadrat tipe III benar. emmeans_model = "multivariate" membuat uji lanjut tetap sahih walaupun asumsi sferisitas dilanggar.

3.2 Membaca data

File Excel dicari otomatis di folder yang sama dengan file .Rmd ini (nama boleh berawalan angka atau berakhiran (2)). Jika ada beberapa file, dipilih yang memiliki sheet Data_Long dan Data_Wide.

kandidat <- list.files(pattern = "Data_PMT_Balita.*\\.xlsx$", ignore.case = TRUE)
kandidat <- kandidat[!startsWith(kandidat, "~$")]   # abaikan file sementara Excel

if (length(kandidat) == 0) {
  stop("File Excel 'Data_PMT_Balita...xlsx' tidak ditemukan di folder: ", getwd(),
       "\nLetakkan file Excel satu folder dengan file .Rmd ini.")
}

skor <- vapply(
  kandidat,
  function(f) {
    s <- excel_sheets(f)
    2 * ("Data_Long" %in% s) + 1 * ("Data_Wide" %in% s)
  },
  numeric(1)
)
if (max(skor) < 2) stop("Tidak ada file Excel dengan sheet 'Data_Long'.")

file_data <- kandidat[which.max(skor)]
message("Membaca file: ", file_data)

# LONG: satu baris = satu penimbangan
data_long <- as.data.frame(read_excel(file_data, sheet = "Data_Long"))

# WIDE: satu baris = satu balita (dari sheet Data_Wide bila ada, jika tidak dibentuk dari LONG)
if ("Data_Wide" %in% excel_sheets(file_data)) {
  data_wide <- as.data.frame(read_excel(file_data, sheet = "Data_Wide"))
  names(data_wide) <- sub("^BB_Bulan", "BB_B", names(data_wide))
} else {
  data_wide <- data_long %>%
    mutate(Waktu_lab = paste0("B", Waktu)) %>%
    dplyr::select(ID, Kelompok, Usia_bulan, Jenis_Kelamin, Waktu_lab, BB_kg) %>%
    pivot_wider(names_from = Waktu_lab, values_from = BB_kg, names_prefix = "BB_") %>%
    as.data.frame()
}

head(data_long)

head(data_wide)

str(data_long)

## 'data.frame':    360 obs. of  6 variables:
##  $ ID           : chr  "BSK-01" "BSK-01" "BSK-01" "BSK-01" ...
##  $ Kelompok     : chr  "PMT Biskuit" "PMT Biskuit" "PMT Biskuit" "PMT Biskuit" ...
##  $ Usia_bulan   : num  24 24 24 24 19 19 19 19 39 39 ...
##  $ Jenis_Kelamin: chr  "Laki-laki" "Laki-laki" "Laki-laki" "Laki-laki" ...
##  $ Waktu        : num  0 1 2 3 0 1 2 3 0 1 ...
##  $ BB_kg        : num  10.2 10.5 10.7 11.1 9.1 9.1 9.5 9.9 15.6 16 ...

str(data_wide)

## 'data.frame':    90 obs. of  8 variables:
##  $ ID           : chr  "KON-01" "KON-02" "KON-03" "KON-04" ...
##  $ Kelompok     : chr  "Kontrol" "Kontrol" "Kontrol" "Kontrol" ...
##  $ Usia_bulan   : num  52 51 19 42 52 59 58 40 15 42 ...
##  $ Jenis_Kelamin: chr  "Laki-laki" "Perempuan" "Perempuan" "Perempuan" ...
##  $ BB_B0        : num  15.4 14.3 9.6 12.7 15.3 15.7 17.4 12.6 8.3 15.4 ...
##  $ BB_B1        : num  15.3 14.4 9.7 12.9 15.4 16 17.6 12.6 8.3 15.4 ...
##  $ BB_B2        : num  15.6 14.5 9.8 13 15.5 16 17.8 12.8 8.6 15.8 ...
##  $ BB_B3        : num  15.8 14.6 9.9 12.9 15.7 16.1 17.8 13 8.8 16.3 ...

3.3 Validasi struktur data

stopifnot(
  all(c("ID", "Kelompok", "Usia_bulan", "Jenis_Kelamin", "Waktu", "BB_kg") %in% names(data_long)),
  all(c("ID", "Kelompok", "BB_B0", "BB_B1", "BB_B2", "BB_B3") %in% names(data_wide))
)

# Jumlah balita (harus 90) dan per kelompok (harus 30)
n_distinct(data_long$ID)

## [1] 90

n_distinct(data_wide$ID)

## [1] 90

data_wide %>% count(Kelompok)

# Jumlah penimbangan tiap balita (harus 4)
data_long %>% count(ID) %>% count(n, name = "jumlah_balita")

# Duplikasi ID x Waktu (harus 0 baris)
data_long %>% count(ID, Waktu) %>% filter(n != 1)

# Data hilang
colSums(is.na(data_long))

##            ID      Kelompok    Usia_bulan Jenis_Kelamin         Waktu 
##             0             0             0             0             0 
##         BB_kg 
##             0

colSums(is.na(data_wide))

##            ID      Kelompok    Usia_bulan Jenis_Kelamin         BB_B0 
##             0             0             0             0             0 
##         BB_B1         BB_B2         BB_B3 
##             0             0             0

# Ringkasan outcome
summary(data_long$BB_kg)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.80   10.80   12.80   12.88   14.72   18.30

Hasil. Data lengkap: 90 balita (30 per kelompok), masing-masing 4 penimbangan (360 baris pada format LONG), tanpa duplikasi dan tanpa data hilang. Berat badan berkisar 7,8–18,3 kg dengan rerata 12,88 kg. Struktur ini seimbang (balanced), sehingga memenuhi syarat untuk RM ANOVA.

3.4 Pengkodean ulang variabel

Waktu_num tetap numerik (untuk grafik dan LMM), sedangkan Waktu dijadikan faktor (untuk RM ANOVA).

kel_lab <- c("Kontrol", "PMT Biskuit", "PMT Lokal")

data_long <- data_long %>%
  mutate(
    ID = factor(ID),
    Kelompok = factor(Kelompok, levels = kel_lab),
    Jenis_Kelamin = factor(Jenis_Kelamin),
    Waktu_num = as.numeric(Waktu),
    Waktu = factor(Waktu, levels = c(0, 1, 2, 3), labels = c("B0", "B1", "B2", "B3"))
  )

data_wide <- data_wide %>%
  mutate(
    ID = factor(ID),
    Kelompok = factor(Kelompok, levels = kel_lab),
    Jenis_Kelamin = factor(Jenis_Kelamin)
  )

levels(data_long$Kelompok)

## [1] "Kontrol"     "PMT Biskuit" "PMT Lokal"

levels(data_long$Waktu)

## [1] "B0" "B1" "B2" "B3"

3.5 Konsistensi format LONG dan WIDE

wide_from_long <- data_long %>%
  dplyr::select(ID, Kelompok, Waktu, BB_kg) %>%
  pivot_wider(names_from = Waktu, values_from = BB_kg, names_prefix = "BB_")

nrow(wide_from_long)

## [1] 90

nrow(data_wide)

## [1] 90

kol_bb <- c("BB_B0", "BB_B1", "BB_B2", "BB_B3")
urut   <- match(as.character(data_wide$ID), as.character(wide_from_long$ID))
stopifnot(!anyNA(urut))

# Hasil yang diharapkan: TRUE
all.equal(
  as.matrix(wide_from_long[urut, kol_bb]),
  as.matrix(data_wide[, kol_bb]),
  check.attributes = FALSE
)

## [1] TRUE

Nilai pada format WIDE identik dengan format LONG (TRUE), sehingga kedua format dapat dipakai bergantian.

3.6 Karakteristik awal: apakah kelompok setara?

karakteristik <- data_wide %>%
  group_by(Kelompok) %>%
  summarise(
    n = n(),
    usia_mean = mean(Usia_bulan), usia_sd = sd(Usia_bulan),
    laki = sum(Jenis_Kelamin == "Laki-laki"),
    perempuan = sum(Jenis_Kelamin == "Perempuan"),
    bb0_mean = mean(BB_B0), bb0_sd = sd(BB_B0),
    .groups = "drop"
  )
knitr::kable(karakteristik, digits = 2,
             caption = "Karakteristik awal balita menurut kelompok")

Karakteristik awal balita menurut kelompok
Kelompok	n	usia_mean	usia_sd	laki	perempuan	bb0_mean	bb0_sd
Kontrol	30	41.9	13.90	12	18	12.97	2.36
PMT Biskuit	30	34.5	13.59	19	11	12.14	2.39
PMT Lokal	30	35.9	14.28	19	11	12.32	2.47

chisq.test(table(data_wide$Kelompok, data_wide$Jenis_Kelamin))   # jenis kelamin

## 
##  Pearson's Chi-squared test
## 
## data:  table(data_wide$Kelompok, data_wide$Jenis_Kelamin)
## X-squared = 4.41, df = 2, p-value = 0.1103

summary(aov(Usia_bulan ~ Kelompok, data = data_wide))             # usia

##             Df Sum Sq Mean Sq F value Pr(>F)  
## Kelompok     2    927   463.6    2.39 0.0976 .
## Residuals   87  16875   194.0                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(BB_B0 ~ Kelompok, data = data_wide))                  # BB awal

##             Df Sum Sq Mean Sq F value Pr(>F)
## Kelompok     2   11.6   5.802       1  0.372
## Residuals   87  504.6   5.800

Interpretasi. Ketiga kelompok relatif setara pada awal pengamatan: jenis kelamin (χ² = 4,41; p = 0,110), usia (F = 2,39; p = 0,098), dan berat badan awal (F = 1,00; p = 0,372), semuanya p > 0,05. Usia rata-rata kelompok Kontrol (41,9 bulan) agak lebih tinggi dibandingkan PMT Biskuit (34,5) dan PMT Lokal (35,9). Karena p usia mendekati batas 0,05, usia dan jenis kelamin tetap dikontrol pada ANCOVA dan LMM (bagian 7 dan 8).

4. Eksplorasi Data

4.1 Statistik deskriptif

deskriptif <- data_long %>%
  group_by(Kelompok, Waktu) %>%
  summarise(
    n = n(),
    mean = mean(BB_kg, na.rm = TRUE),
    sd = sd(BB_kg, na.rm = TRUE),
    median = median(BB_kg, na.rm = TRUE),
    min = min(BB_kg, na.rm = TRUE),
    max = max(BB_kg, na.rm = TRUE),
    .groups = "drop"
  )
knitr::kable(deskriptif, digits = 2,
             caption = "Berat badan (kg) per kelompok dan waktu")

Berat badan (kg) per kelompok dan waktu
Kelompok	Waktu	n	mean	sd	median	min	max
Kontrol	B0	30	12.97	2.36	13.00	8.3	17.4
Kontrol	B1	30	13.09	2.37	13.20	8.3	17.6
Kontrol	B2	30	13.21	2.40	13.25	8.6	17.8
Kontrol	B3	30	13.31	2.42	13.45	8.8	17.8
PMT Biskuit	B0	30	12.14	2.39	12.05	7.8	16.0
PMT Biskuit	B1	30	12.45	2.46	12.40	7.9	16.7
PMT Biskuit	B2	30	12.73	2.47	12.60	8.3	17.1
PMT Biskuit	B3	30	13.08	2.47	13.00	8.6	17.4
PMT Lokal	B0	30	12.32	2.47	12.20	8.3	17.3
PMT Lokal	B1	30	12.68	2.43	12.55	8.7	17.5
PMT Lokal	B2	30	13.06	2.44	12.70	9.2	17.9
PMT Lokal	B3	30	13.45	2.44	13.05	9.6	18.3

Rerata berat badan pada bulan 0 sampai bulan 3: Kontrol 12,97 → 13,31 kg; PMT Biskuit 12,14 → 13,08 kg; PMT Lokal 12,32 → 13,45 kg. Simpangan baku sekitar 2,4 kg pada semua sel, yang mencerminkan perbedaan usia dan ukuran tubuh antarbalita.

4.2 Matriks kovarians, korelasi, dan varians selisih antarwaktu

vars_waktu <- c("BB_B0", "BB_B1", "BB_B2", "BB_B3")

S <- cov(data_wide[, vars_waktu], use = "complete.obs")
R <- cor(data_wide[, vars_waktu], use = "complete.obs")
round(S, 2)

##       BB_B0 BB_B1 BB_B2 BB_B3
## BB_B0  5.80  5.79  5.79  5.75
## BB_B1  5.79  5.80  5.81  5.78
## BB_B2  5.79  5.81  5.86  5.84
## BB_B3  5.75  5.78  5.84  5.85

round(R, 3)

##       BB_B0 BB_B1 BB_B2 BB_B3
## BB_B0 1.000 0.998 0.993 0.986
## BB_B1 0.998 1.000 0.998 0.993
## BB_B2 0.993 0.998 1.000 0.997
## BB_B3 0.986 0.993 0.997 1.000

pasangan <- combn(vars_waktu, 2)
var_selisih <- apply(
  pasangan, 2,
  function(p) var(data_wide[[p[1]]] - data_wide[[p[2]]], na.rm = TRUE)
)
names(var_selisih) <- apply(pasangan, 2, paste, collapse = " - ")
round(var_selisih, 3)

## BB_B0 - BB_B1 BB_B0 - BB_B2 BB_B0 - BB_B3 BB_B1 - BB_B2 BB_B1 - BB_B3 
##         0.027         0.078         0.158         0.028         0.087 
## BB_B2 - BB_B3 
##         0.031

Interpretasi. Berat badan antarwaktu berkorelasi sangat tinggi (r = 0,986–0,998), artinya balita yang berat sejak awal tetap berat pada bulan berikutnya. Varians selisih antarpasangan waktu tidak seragam: 0,027 (B0−B1) sampai 0,158 (B0−B3), dan membesar seiring jarak waktu. Ini petunjuk awal bahwa asumsi sferisitas dilanggar, dan akan dikonfirmasi oleh uji Mauchly pada bagian 5 dan 6.

4.3 Profile plot

p_profil <- ggplot(
  data_long,
  aes(x = Waktu_num, y = BB_kg, colour = Kelompok, group = Kelompok)
) +
  stat_summary(fun = mean, geom = "line", linewidth = 1) +
  stat_summary(fun = mean, geom = "point", size = 2.8) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = .12) +
  scale_x_continuous(breaks = c(0, 1, 2, 3)) +
  labs(
    x = "Bulan penimbangan", y = "Berat badan (kg)", colour = "Kelompok",
    title = "Profil rerata berat badan balita selama 3 bulan"
  ) +
  theme(legend.position = "bottom")

p_profil

Ketiga garis naik, tetapi kemiringannya berbeda: PMT Lokal paling curam, PMT Biskuit di tengah, dan Kontrol paling landai. Garis yang tidak sejajar ini menunjukkan adanya interaksi kelompok × waktu.

4.4 Spaghetti plot

p_spag <- ggplot(data_long, aes(x = Waktu_num, y = BB_kg, group = ID)) +
  geom_line(alpha = .25) +
  stat_summary(
    aes(group = Kelompok), fun = mean, geom = "line",
    linewidth = 1.3, colour = "firebrick"
  ) +
  facet_wrap(~ Kelompok) +
  scale_x_continuous(breaks = c(0, 1, 2, 3)) +
  labs(x = "Bulan", y = "Berat badan (kg)",
       title = "Lintasan individu dan rerata kelompok")

p_spag

Lintasan individu hampir sejajar satu sama lain (perbedaan terutama pada level awal), dan hampir semua balita naik. Hal ini sejalan dengan ICC yang sangat tinggi pada bagian 8.

5. RM ANOVA Satu Arah: Contoh Kelompok PMT Lokal

Pertanyaan: apakah berat badan balita kelompok PMT Lokal berubah selama tiga bulan? (Ganti kelompok_fokus dengan “Kontrol” atau “PMT Biskuit” untuk kelompok lain.)

kelompok_fokus <- "PMT Lokal"

d1  <- data_long %>% filter(Kelompok == kelompok_fokus) %>% droplevels()
d1w <- data_wide %>% filter(Kelompok == kelompok_fokus) %>% droplevels()

5.1 Uji asumsi

# (i) Outlier per waktu
outlier_1way <- d1 %>% group_by(Waktu) %>% identify_outliers(BB_kg)
outlier_1way

# (ii) Normalitas Shapiro-Wilk per waktu
shapiro_1way <- d1 %>% group_by(Waktu) %>% shapiro_test(BB_kg)
shapiro_1way

# Q-Q plot
ggpubr::ggqqplot(d1, "BB_kg", facet.by = "Waktu")

# (iii) Sferisitas (Mauchly) + koreksi GG/HF
aov1_rs <- anova_test(
  data = d1, dv = BB_kg, wid = ID, within = Waktu, effect.size = "pes"
)
aov1_rs

## ANOVA Table (type III tests)
## 
## $ANOVA
##   Effect DFn DFd       F        p p<.05   pes
## 1  Waktu   3  87 602.823 4.51e-58     * 0.954
## 
## $`Mauchly's Test for Sphericity`
##   Effect     W        p p<.05
## 1  Waktu 0.323 8.07e-06     *
## 
## $`Sphericity Corrections`
##   Effect   GGe     DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]   p[HF]
## 1  Waktu 0.637 1.91, 55.4 7.14e-38         * 0.681 2.04, 59.21 2.6e-40
##   p[HF]<.05
## 1         *

get_anova_table(aov1_rs, correction = "auto")

Interpretasi. Tidak ada outlier, dan data normal pada keempat waktu (Shapiro-Wilk p = 0,232–0,353). Uji Mauchly signifikan (W = 0,323; p < 0,001), sehingga sferisitas tidak terpenuhi dan dipakai koreksi Greenhouse-Geisser (ε = 0,637; Huynh-Feldt ε = 0,681).

5.2 RM ANOVA dengan afex

aov1 <- aov_ez(
  id = "ID", dv = "BB_kg", data = d1, within = "Waktu",
  anova_table = list(es = c("ges", "pes"), correction = "GG")
)
aov1

## Anova Table (Type 3 tests)
## 
## Response: BB_kg
##   Effect          df  MSE          F  ges  pes p.value
## 1  Waktu 1.91, 55.40 0.02 602.82 *** .030 .954   <.001
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
## 
## Sphericity correction method: GG

summary(aov1)

## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##              Sum Sq num Df Error SS den Df F value    Pr(>F)    
## (Intercept) 19902.2      1   693.01     29  832.83 < 2.2e-16 ***
## Waktu          21.2      3     1.02     87  602.82 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##       Test statistic    p-value
## Waktu        0.32253 8.0703e-06
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##        GG eps Pr(>F[GG])    
## Waktu 0.63682  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          HF eps   Pr(>F[HF])
## Waktu 0.6806248 2.602525e-40

nice(aov1, correction = "GG", es = c("ges", "pes"))

# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)

omega_squared(aov1, partial = TRUE)

Interpretasi. Berat badan berubah bermakna selama tiga bulan pada kelompok PMT Lokal: F(1,91; 55,40) = 602,82; p < 0,001; η² parsial = 0,954 (efek sangat besar).

5.3 Pendekatan multivariat

aov1$Anova

## 
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
##             Df test stat approx F num Df den Df    Pr(>F)    
## (Intercept)  1   0.96635   832.83      1     29 < 2.2e-16 ***
## Waktu        1   0.97729   387.22      3     27 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Uji multivariat (Pillai) tidak membutuhkan sferisitas dan memberi kesimpulan yang sama: F(3, 27) = 387,22; p < 0,001.

5.4 Uji lanjut dan tren

em1 <- emmeans(aov1, ~ Waktu)
em1

##  Waktu emmean    SE df lower.CL upper.CL
##  B0      12.3 0.452 29     11.4     13.2
##  B1      12.7 0.443 29     11.8     13.6
##  B2      13.1 0.446 29     12.2     14.0
##  B3      13.4 0.445 29     12.5     14.4
## 
## Confidence level used: 0.95

# Semua pasangan waktu
pairs(em1, adjust = "holm")

##  contrast estimate     SE df t.ratio p.value
##  B0 - B1    -0.363 0.0200 29 -18.123 <0.0001
##  B0 - B2    -0.743 0.0341 29 -21.777 <0.0001
##  B0 - B3    -1.127 0.0346 29 -32.607 <0.0001
##  B1 - B2    -0.380 0.0232 29 -16.384 <0.0001
##  B1 - B3    -0.763 0.0305 29 -25.022 <0.0001
##  B2 - B3    -0.383 0.0215 29 -17.840 <0.0001
## 
## P value adjustment: holm method for 6 tests

# Tiap waktu vs baseline B0
contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm")

##  contrast estimate     SE df t.ratio p.value
##  B1 - B0     0.363 0.0200 29  18.123 <0.0001
##  B2 - B0     0.743 0.0341 29  21.777 <0.0001
##  B3 - B0     1.127 0.0346 29  32.607 <0.0001
## 
## P value adjustment: holm method for 3 tests

# Tren linear, kuadratik, kubik
contrast(em1, "poly")

##  contrast  estimate     SE df t.ratio p.value
##  linear      3.7600 0.1220 29  30.721 <0.0001
##  quadratic   0.0200 0.0350 29   0.571  0.5725
##  cubic      -0.0133 0.0484 29  -0.276  0.7847

Interpretasi. Semua pasangan waktu berbeda bermakna (Holm, p < 0,0001). Dibandingkan bulan 0, berat badan naik 0,36 kg (bulan 1), 0,74 kg (bulan 2), dan 1,13 kg (bulan 3). Polanya linear (kontras linear = 3,76; p < 0,0001), sedangkan komponen kuadratik (p = 0,573) dan kubik (p = 0,785) tidak signifikan, artinya kenaikan berat badan stabil sekitar 0,38 kg per bulan.

5.5 Alternatif nonparametrik: Friedman

friedman_test(d1, BB_kg ~ Waktu | ID)

friedman_effsize(d1, BB_kg ~ Waktu | ID)

d1 %>% wilcox_test(BB_kg ~ Waktu, paired = TRUE, p.adjust.method = "holm")

Uji Friedman memberi kesimpulan yang sama: χ²(3) = 90,0; p < 0,001; Kendall’s W = 1,00 (kesepakatan sempurna: seluruh balita naik berurutan dari bulan ke bulan). Seluruh pasangan Wilcoxon juga signifikan setelah koreksi Holm.

6. Mixed Design ANOVA: Kelompok × Waktu

Ini adalah uji utama untuk menjawab: apakah perubahan berat badan dari bulan 0 sampai bulan 3 berbeda antara Kontrol, PMT Biskuit, dan PMT Lokal?

6.1 Uji asumsi

# (i) Outlier per sel
outlier_mixed <- data_long %>% group_by(Kelompok, Waktu) %>% identify_outliers(BB_kg)
outlier_mixed

# (ii) Normalitas per Kelompok x Waktu
shapiro_mixed <- data_long %>% group_by(Kelompok, Waktu) %>% shapiro_test(BB_kg)
shapiro_mixed

ggpubr::ggqqplot(data_long, "BB_kg") + facet_grid(Waktu ~ Kelompok)

# (iii) Homogenitas varians pada tiap waktu: Levene
levene_mixed <- data_long %>% group_by(Waktu) %>% levene_test(BB_kg ~ Kelompok)
levene_mixed

# (iv) Homogenitas matriks kovarians: Box's M
box_m_result <- box_m(data_wide[, vars_waktu], data_wide$Kelompok)
box_m_result

Asumsi	Hasil	Kesimpulan
Outlier	Tidak ada	Terpenuhi
Normalitas (12 sel)	Shapiro-Wilk p = 0,232–0,409	Terpenuhi
Homogenitas varians (Levene)	p = 0,845; 0,946; 0,950; 0,962	Terpenuhi
Homogenitas kovarians (Box’s M)	M = 34,3; p = 0,0245	Terpenuhi (ambang uji Box’s M umumnya p < 0,001)
Sferisitas (Mauchly)	Lihat 6.2	Dilanggar, pakai koreksi GG

6.2 Mixed ANOVA

aov2 <- aov_ez(
  id = "ID", dv = "BB_kg", data = data_long,
  between = "Kelompok", within = "Waktu",
  anova_table = list(es = c("ges", "pes"), correction = "GG")
)
aov2

## Anova Table (Type 3 tests)
## 
## Response: BB_kg
##           Effect           df   MSE          F  ges  pes p.value
## 1       Kelompok        2, 87 23.53       0.38 .009 .009    .687
## 2          Waktu 1.98, 172.38  0.02 803.10 *** .015 .902   <.001
## 3 Kelompok:Waktu 3.96, 172.38  0.02  70.39 *** .003 .618   <.001
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
## 
## Sphericity correction method: GG

summary(aov2)

## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                Sum Sq num Df Error SS den Df   F value Pr(>F)    
## (Intercept)     59678      1  2046.98     87 2536.4165 <2e-16 ***
## Kelompok           18      2  2046.98     87    0.3776 0.6866    
## Waktu              32      3     3.48    261  803.0987 <2e-16 ***
## Kelompok:Waktu      6      6     3.48    261   70.3855 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##                Test statistic    p-value
## Waktu                 0.42937 3.1623e-14
## Kelompok:Waktu        0.42937 3.1623e-14
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                 GG eps Pr(>F[GG])    
## Waktu          0.66046  < 2.2e-16 ***
## Kelompok:Waktu 0.66046  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   HF eps   Pr(>F[HF])
## Waktu          0.6757837 9.223254e-90
## Kelompok:Waktu 0.6757837 8.423267e-36

nice(aov2, correction = "GG", es = c("ges", "pes"))

# Versi rstatix (hasil setara)
aov2_rs <- anova_test(
  data = data_long, dv = BB_kg, wid = ID,
  between = Kelompok, within = Waktu, effect.size = "pes", type = 3
)
aov2_rs

## ANOVA Table (type III tests)
## 
## $ANOVA
##           Effect DFn DFd       F         p p<.05   pes
## 1       Kelompok   2  87   0.378  6.87e-01       0.009
## 2          Waktu   3 261 803.099 1.97e-131     * 0.902
## 3 Kelompok:Waktu   6 261  70.386  9.55e-52     * 0.618
## 
## $`Mauchly's Test for Sphericity`
##           Effect     W        p p<.05
## 1          Waktu 0.429 3.16e-14     *
## 2 Kelompok:Waktu 0.429 3.16e-14     *
## 
## $`Sphericity Corrections`
##           Effect  GGe       DF[GG]    p[GG] p[GG]<.05   HFe       DF[HF]
## 1          Waktu 0.66 1.98, 172.38 8.60e-88         * 0.676 2.03, 176.38
## 2 Kelompok:Waktu 0.66 3.96, 172.38 4.78e-35         * 0.676 4.05, 176.38
##      p[HF] p[HF]<.05
## 1 9.22e-90         *
## 2 8.42e-36         *

get_anova_table(aov2_rs, correction = "GG")

Efek	F (df, GG)	p	η² parsial
Kelompok	0,38 (2; 87)	0,687	0,009
Waktu	803,10 (1,98; 172,38)	< 0,001	0,902
Kelompok × Waktu	70,39 (3,96; 172,38)	< 0,001	0,618

Interpretasi.

Sferisitas dilanggar (Mauchly W = 0,429; p < 0,001; ε GG = 0,660), jadi tabel yang dilaporkan adalah yang terkoreksi Greenhouse-Geisser.
Efek waktu signifikan: berat badan berubah selama tiga bulan.
Interaksi kelompok × waktu signifikan, sehingga H0 ditolak: laju kenaikan berat badan berbeda antarkelompok.
Efek utama kelompok tidak signifikan (p = 0,687). Ini bukan berarti PMT tidak berpengaruh: efek utama merata-ratakan seluruh waktu, sedangkan ketiga kelompok mulai dari berat badan yang mirip dan baru terpisah seiring waktu. Karena interaksi signifikan, interpretasi dilanjutkan dengan simple effects.

6.3 Ukuran efek

eta_squared(aov2, partial = TRUE)

omega_squared(aov2, partial = TRUE)

η² parsial untuk waktu (0,90) dan interaksi (0,62) termasuk besar. Perhatikan bahwa η² generalized (ges) jauh lebih kecil (waktu 0,015; interaksi 0,003) karena variasi berat badan antarbalita sangat besar (ICC 0,97). Jangan membandingkan angka partial dan generalized secara langsung, dan sebutkan jenis ukuran efek yang dilaporkan.

6.4 Plot interaksi

p_interaksi <- afex_plot(
  aov2, x = "Waktu", trace = "Kelompok", error = "within",
  mapping = c("colour", "shape", "linetype")
) +
  labs(
    y = "Berat badan (kg)", x = "Waktu",
    title = "Interaksi Kelompok x Waktu pada berat badan balita"
  ) +
  theme(legend.position = "bottom")

p_interaksi

Catatan: afex_plot memberi peringatan bahwa error bar pada desain campuran tidak dapat dipakai untuk membandingkan seluruh rerata. Gunakan error bar hanya sebagai gambaran dan andalkan uji lanjut untuk inferensi.

6.5 Simple effects

# Rerata marginal waktu dalam tiap kelompok
em2 <- emmeans(aov2, ~ Waktu | Kelompok)
em2

## Kelompok = Kontrol:
##  Waktu emmean    SE df lower.CL upper.CL
##  B0      13.0 0.440 87     12.1     13.8
##  B1      13.1 0.442 87     12.2     14.0
##  B2      13.2 0.445 87     12.3     14.1
##  B3      13.3 0.446 87     12.4     14.2
## 
## Kelompok = PMT Biskuit:
##  Waktu emmean    SE df lower.CL upper.CL
##  B0      12.1 0.440 87     11.3     13.0
##  B1      12.5 0.442 87     11.6     13.3
##  B2      12.7 0.445 87     11.8     13.6
##  B3      13.1 0.446 87     12.2     14.0
## 
## Kelompok = PMT Lokal:
##  Waktu emmean    SE df lower.CL upper.CL
##  B0      12.3 0.440 87     11.4     13.2
##  B1      12.7 0.442 87     11.8     13.6
##  B2      13.1 0.445 87     12.2     13.9
##  B3      13.4 0.446 87     12.6     14.3
## 
## Confidence level used: 0.95

# Efek waktu di dalam masing-masing kelompok
joint_tests(aov2, by = "Kelompok")

# Efek kelompok pada masing-masing waktu
joint_tests(aov2, by = "Waktu")

Efek waktu dalam kelompok: signifikan pada ketiga kelompok (p < 0,0001), dengan kekuatan yang berbeda: F = 26,4 (Kontrol), 211,4 (PMT Biskuit), dan 290,4 (PMT Lokal).
Efek kelompok pada tiap waktu: tidak signifikan pada semua waktu (p = 0,372; 0,589; 0,744; 0,844). Artinya, berat badan absolut antarkelompok tidak berbeda. Perbedaan terletak pada besarnya kenaikan, bukan pada level berat badan.

6.6 Uji lanjut dalam kelompok

# Semua pasangan waktu dalam tiap kelompok
pairs(em2, adjust = "holm")

## Kelompok = Kontrol:
##  contrast estimate     SE df t.ratio p.value
##  B0 - B1    -0.117 0.0233 87  -5.014 <0.0001
##  B0 - B2    -0.233 0.0330 87  -7.061 <0.0001
##  B0 - B3    -0.340 0.0385 87  -8.835 <0.0001
##  B1 - B2    -0.117 0.0236 87  -4.937 <0.0001
##  B1 - B3    -0.223 0.0338 87  -6.602 <0.0001
##  B2 - B3    -0.107 0.0228 87  -4.681 <0.0001
## 
## Kelompok = PMT Biskuit:
##  contrast estimate     SE df t.ratio p.value
##  B0 - B1    -0.317 0.0233 87 -13.610 <0.0001
##  B0 - B2    -0.597 0.0330 87 -18.056 <0.0001
##  B0 - B3    -0.947 0.0385 87 -24.599 <0.0001
##  B1 - B2    -0.280 0.0236 87 -11.848 <0.0001
##  B1 - B3    -0.630 0.0338 87 -18.625 <0.0001
##  B2 - B3    -0.350 0.0228 87 -15.359 <0.0001
## 
## Kelompok = PMT Lokal:
##  contrast estimate     SE df t.ratio p.value
##  B0 - B1    -0.363 0.0233 87 -15.615 <0.0001
##  B0 - B2    -0.743 0.0330 87 -22.495 <0.0001
##  B0 - B3    -1.127 0.0385 87 -29.277 <0.0001
##  B1 - B2    -0.380 0.0236 87 -16.080 <0.0001
##  B1 - B3    -0.763 0.0338 87 -22.567 <0.0001
##  B2 - B3    -0.383 0.0228 87 -16.822 <0.0001
## 
## P value adjustment: holm method for 6 tests

# Tiap waktu vs baseline B0
contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm")

## Kelompok = Kontrol:
##  contrast estimate     SE df t.ratio p.value
##  B1 - B0     0.117 0.0233 87   5.014 <0.0001
##  B2 - B0     0.233 0.0330 87   7.061 <0.0001
##  B3 - B0     0.340 0.0385 87   8.835 <0.0001
## 
## Kelompok = PMT Biskuit:
##  contrast estimate     SE df t.ratio p.value
##  B1 - B0     0.317 0.0233 87  13.610 <0.0001
##  B2 - B0     0.597 0.0330 87  18.056 <0.0001
##  B3 - B0     0.947 0.0385 87  24.599 <0.0001
## 
## Kelompok = PMT Lokal:
##  contrast estimate     SE df t.ratio p.value
##  B1 - B0     0.363 0.0233 87  15.615 <0.0001
##  B2 - B0     0.743 0.0330 87  22.495 <0.0001
##  B3 - B0     1.127 0.0385 87  29.277 <0.0001
## 
## P value adjustment: holm method for 3 tests

Pada ketiga kelompok, setiap pasangan waktu berbeda bermakna (Holm, p < 0,0001). Kenaikan bulan 3 dibandingkan bulan 0: Kontrol 0,34 kg, PMT Biskuit 0,95 kg, PMT Lokal 1,13 kg.

6.7 Uji lanjut antarkelompok pada tiap waktu

em2b <- emmeans(aov2, ~ Kelompok | Waktu)
em2b

## Waktu = B0:
##  Kelompok    emmean    SE df lower.CL upper.CL
##  Kontrol       13.0 0.440 87     12.1     13.8
##  PMT Biskuit   12.1 0.440 87     11.3     13.0
##  PMT Lokal     12.3 0.440 87     11.4     13.2
## 
## Waktu = B1:
##  Kelompok    emmean    SE df lower.CL upper.CL
##  Kontrol       13.1 0.442 87     12.2     14.0
##  PMT Biskuit   12.5 0.442 87     11.6     13.3
##  PMT Lokal     12.7 0.442 87     11.8     13.6
## 
## Waktu = B2:
##  Kelompok    emmean    SE df lower.CL upper.CL
##  Kontrol       13.2 0.445 87     12.3     14.1
##  PMT Biskuit   12.7 0.445 87     11.8     13.6
##  PMT Lokal     13.1 0.445 87     12.2     13.9
## 
## Waktu = B3:
##  Kelompok    emmean    SE df lower.CL upper.CL
##  Kontrol       13.3 0.446 87     12.4     14.2
##  PMT Biskuit   13.1 0.446 87     12.2     14.0
##  PMT Lokal     13.4 0.446 87     12.6     14.3
## 
## Confidence level used: 0.95

pairs(em2b, adjust = "tukey")

## Waktu = B0:
##  contrast                estimate    SE df t.ratio p.value
##  Kontrol - PMT Biskuit      0.837 0.622 87   1.345  0.3741
##  Kontrol - PMT Lokal        0.653 0.622 87   1.051  0.5472
##  PMT Biskuit - PMT Lokal   -0.183 0.622 87  -0.295  0.9532
## 
## Waktu = B1:
##  contrast                estimate    SE df t.ratio p.value
##  Kontrol - PMT Biskuit      0.637 0.625 87   1.019  0.5672
##  Kontrol - PMT Lokal        0.407 0.625 87   0.651  0.7925
##  PMT Biskuit - PMT Lokal   -0.230 0.625 87  -0.368  0.9281
## 
## Waktu = B2:
##  contrast                estimate    SE df t.ratio p.value
##  Kontrol - PMT Biskuit      0.473 0.630 87   0.752  0.7335
##  Kontrol - PMT Lokal        0.143 0.630 87   0.228  0.9719
##  PMT Biskuit - PMT Lokal   -0.330 0.630 87  -0.524  0.8598
## 
## Waktu = B3:
##  contrast                estimate    SE df t.ratio p.value
##  Kontrol - PMT Biskuit      0.230 0.630 87   0.365  0.9293
##  Kontrol - PMT Lokal       -0.133 0.630 87  -0.212  0.9756
##  PMT Biskuit - PMT Lokal   -0.363 0.630 87  -0.576  0.8330
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

pairs(em2b, adjust = "holm")

## Waktu = B0:
##  contrast                estimate    SE df t.ratio p.value
##  Kontrol - PMT Biskuit      0.837 0.622 87   1.345  0.5459
##  Kontrol - PMT Lokal        0.653 0.622 87   1.051  0.5927
##  PMT Biskuit - PMT Lokal   -0.183 0.622 87  -0.295  0.7688
## 
## Waktu = B1:
##  contrast                estimate    SE df t.ratio p.value
##  Kontrol - PMT Biskuit      0.637 0.625 87   1.019  0.9336
##  Kontrol - PMT Lokal        0.407 0.625 87   0.651  1.0000
##  PMT Biskuit - PMT Lokal   -0.230 0.625 87  -0.368  1.0000
## 
## Waktu = B2:
##  contrast                estimate    SE df t.ratio p.value
##  Kontrol - PMT Biskuit      0.473 0.630 87   0.752  1.0000
##  Kontrol - PMT Lokal        0.143 0.630 87   0.228  1.0000
##  PMT Biskuit - PMT Lokal   -0.330 0.630 87  -0.524  1.0000
## 
## Waktu = B3:
##  contrast                estimate    SE df t.ratio p.value
##  Kontrol - PMT Biskuit      0.230 0.630 87   0.365  1.0000
##  Kontrol - PMT Lokal       -0.133 0.630 87  -0.212  1.0000
##  PMT Biskuit - PMT Lokal   -0.363 0.630 87  -0.576  1.0000
## 
## P value adjustment: holm method for 3 tests

Tidak ada perbedaan bermakna antarkelompok pada B0, B1, B2, maupun B3 (Tukey p = 0,374–0,976). Hasil ini konsisten dengan efek kelompok yang tidak signifikan pada 6.5.

6.8 Kontras perubahan bulan 0 → bulan 3

Inilah jawaban langsung atas pertanyaan penelitian: apakah besarnya kenaikan berat badan berbeda antarkelompok?

em_full <- emmeans(aov2, ~ Waktu * Kelompok)
em_full

##  Waktu Kelompok    emmean    SE df lower.CL upper.CL
##  B0    Kontrol       13.0 0.440 87     12.1     13.8
##  B1    Kontrol       13.1 0.442 87     12.2     14.0
##  B2    Kontrol       13.2 0.445 87     12.3     14.1
##  B3    Kontrol       13.3 0.446 87     12.4     14.2
##  B0    PMT Biskuit   12.1 0.440 87     11.3     13.0
##  B1    PMT Biskuit   12.5 0.442 87     11.6     13.3
##  B2    PMT Biskuit   12.7 0.445 87     11.8     13.6
##  B3    PMT Biskuit   13.1 0.446 87     12.2     14.0
##  B0    PMT Lokal     12.3 0.440 87     11.4     13.2
##  B1    PMT Lokal     12.7 0.442 87     11.8     13.6
##  B2    PMT Lokal     13.1 0.445 87     12.2     13.9
##  B3    PMT Lokal     13.4 0.446 87     12.6     14.3
## 
## Confidence level used: 0.95

# Pastikan urutan grid: B0-B3 untuk Kontrol, lalu PMT Biskuit, lalu PMT Lokal
stopifnot(
  identical(as.character(em_full@grid$Waktu), rep(c("B0", "B1", "B2", "B3"), 3)),
  identical(as.character(em_full@grid$Kelompok), rep(kel_lab, each = 4))
)

kc <- c(-1, 0, 0, 1)   # kontras B3 - B0
z  <- rep(0, 4)

# (i) Kenaikan BB (B3 - B0) di dalam tiap kelompok
contrast(
  em_full,
  list(
    "Kontrol: B3-B0"     = c(kc, z, z),
    "PMT Biskuit: B3-B0" = c(z, kc, z),
    "PMT Lokal: B3-B0"   = c(z, z, kc)
  ),
  adjust = "none"
)

##  contrast           estimate     SE df t.ratio p.value
##  Kontrol: B3-B0        0.340 0.0385 87   8.835 <0.0001
##  PMT Biskuit: B3-B0    0.947 0.0385 87  24.599 <0.0001
##  PMT Lokal: B3-B0      1.127 0.0385 87  29.277 <0.0001

# (ii) Perbedaan kenaikan antarkelompok
contrast(
  em_full,
  list(
    "Biskuit - Kontrol" = c(-kc, kc, z),
    "Lokal - Kontrol"   = c(-kc, z, kc),
    "Lokal - Biskuit"   = c(z, -kc, kc)
  ),
  adjust = "holm"
)

##  contrast          estimate     SE df t.ratio p.value
##  Biskuit - Kontrol    0.607 0.0544 87  11.147 <0.0001
##  Lokal - Kontrol      0.787 0.0544 87  14.454 <0.0001
##  Lokal - Biskuit      0.180 0.0544 87   3.307  0.0014
## 
## P value adjustment: holm method for 3 tests

Perbandingan kenaikan BB (B3 − B0)	Selisih (kg)	p (Holm)
PMT Biskuit − Kontrol	0,607	< 0,0001
PMT Lokal − Kontrol	0,787	< 0,0001
PMT Lokal − PMT Biskuit	0,180	0,0014

Interpretasi. Balita yang mendapat PMT Biskuit naik 0,61 kg lebih banyak daripada Kontrol, dan PMT Lokal naik 0,79 kg lebih banyak daripada Kontrol. PMT Lokal juga lebih unggul 0,18 kg dibandingkan PMT Biskuit. Ketiga perbedaan signifikan setelah koreksi Holm.

6.9 Tren linear antarkelompok

tren_int <- summary(
  contrast(em_full, interaction = c("poly", "pairwise"), adjust = "none")
)
tren_int

names(tren_int)

## [1] "Waktu_poly"        "Kelompok_pairwise" "estimate"         
## [4] "SE"                "df"                "t.ratio"          
## [7] "p.value"

# Nama kolom dapat berbeda menurut versi emmeans
if ("Waktu_poly" %in% names(tren_int)) {
  tren_lin <- subset(tren_int, Waktu_poly == "linear")
} else if ("Waktu.poly" %in% names(tren_int)) {
  tren_lin <- subset(tren_int, Waktu.poly == "linear")
} else {
  tren_lin <- tren_int
}

if ("p.value" %in% names(tren_lin)) {
  tren_lin$p_holm <- p.adjust(tren_lin$p.value, method = "holm")
}
tren_lin

Laju kenaikan linear berbeda bermakna pada ketiga pasangan kelompok (p Holm: Kontrol vs Biskuit < 0,0001; Kontrol vs Lokal < 0,0001; Biskuit vs Lokal = 0,0009). Komponen kuadratik dan kubik tidak signifikan, jadi perbedaan antarkelompok terutama pada kemiringan (kg per bulan), bukan pada bentuk kurva.

7. Pembanding: Selisih Berat Badan dan ANCOVA

Analisis ini menyederhanakan pertanyaan menjadi satu angka per balita (selisih B3 − B0), sehingga hasilnya mudah dibaca dan dilaporkan.

data_wide <- data_wide %>% mutate(Selisih = BB_B3 - BB_B0)

data_wide %>% group_by(Kelompok) %>% get_summary_stats(Selisih, type = "mean_sd")

# Asumsi
data_wide %>% group_by(Kelompok) %>% shapiro_test(Selisih)

levene_test(data_wide, Selisih ~ Kelompok)

# ANOVA satu arah + Tukey
aov_sel <- aov(Selisih ~ Kelompok, data = data_wide)
summary(aov_sel)

##             Df Sum Sq Mean Sq F value Pr(>F)    
## Kelompok     2 10.193   5.096   114.7 <2e-16 ***
## Residuals   87  3.865   0.044                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov_sel)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Selisih ~ Kelompok, data = data_wide)
## 
## $Kelompok
##                            diff        lwr       upr     p adj
## PMT Biskuit-Kontrol   0.6066667 0.47689442 0.7364389 0.0000000
## PMT Lokal-Kontrol     0.7866667 0.65689442 0.9164389 0.0000000
## PMT Lokal-PMT Biskuit 0.1800000 0.05022775 0.3097722 0.0038808

# Jika asumsi tidak terpenuhi:
# kruskal.test(Selisih ~ Kelompok, data = data_wide)

# Paired t-test B0 vs B3 pada tiap kelompok
data_wide %>%
  pivot_longer(c(BB_B0, BB_B3), names_to = "waktu_uji", values_to = "bb_uji") %>%
  group_by(Kelompok) %>%
  t_test(bb_uji ~ waktu_uji, paired = TRUE)

Interpretasi.

Selisih berat badan (rerata ± SD): Kontrol 0,34 ± 0,21 kg; PMT Biskuit 0,95 ± 0,23 kg; PMT Lokal 1,13 ± 0,19 kg.
Asumsi terpenuhi (Shapiro-Wilk per kelompok p = 0,104–0,179; Levene p = 0,918).
ANOVA: F(2, 87) = 114,7; p < 0,001. Tukey: semua pasangan berbeda (Biskuit−Kontrol 0,61 kg; Lokal−Kontrol 0,79 kg; Lokal−Biskuit 0,18 kg, p = 0,004).
Paired t-test menunjukkan kenaikan bermakna pada ketiga kelompok, termasuk Kontrol (t(29) = −8,76; p < 0,001), tetapi dengan besaran yang jauh lebih kecil.

7.1 ANCOVA: dikontrol berat badan awal, usia, dan jenis kelamin

ancova <- lm(BB_B3 ~ BB_B0 + Usia_bulan + Jenis_Kelamin + Kelompok, data = data_wide)
car::Anova(ancova, type = 3)

emmeans(ancova, pairwise ~ Kelompok, adjust = "tukey")

## $emmeans
##  Kelompok    emmean     SE df lower.CL upper.CL
##  Kontrol       12.8 0.0393 84     12.7     12.9
##  PMT Biskuit   13.4 0.0389 84     13.3     13.5
##  PMT Lokal     13.6 0.0387 84     13.5     13.7
## 
## Results are averaged over the levels of: Jenis_Kelamin 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate     SE df t.ratio p.value
##  Kontrol - PMT Biskuit     -0.592 0.0565 84 -10.471 <0.0001
##  Kontrol - PMT Lokal       -0.774 0.0560 84 -13.811 <0.0001
##  PMT Biskuit - PMT Lokal   -0.182 0.0540 84  -3.364  0.0033
## 
## Results are averaged over the levels of: Jenis_Kelamin 
## P value adjustment: tukey method for comparing a family of 3 estimates

Setelah berat badan awal, usia, dan jenis kelamin dikontrol, kelompok tetap berpengaruh terhadap berat badan bulan 3 (F(2, 84) = 101,77; p < 0,001). Rerata terkoreksi: Kontrol 12,8 kg; PMT Biskuit 13,4 kg; PMT Lokal 13,6 kg, dan semua pasangan berbeda (Tukey p ≤ 0,0033). Berat badan awal sangat berpengaruh (p < 0,001), usia mendekati signifikan (p = 0,065), dan jenis kelamin tidak berpengaruh (p = 0,955).

8. Pembanding: Linear Mixed Model (LMM)

LMM tidak mensyaratkan sferisitas, dapat menampung data hilang, dan memungkinkan usia serta jenis kelamin dimasukkan sebagai kovariat.

# Random intercept
lmm1 <- lmer(BB_kg ~ Kelompok * Waktu + (1 | ID), data = data_long, REML = TRUE)
summary(lmm1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BB_kg ~ Kelompok * Waktu + (1 | ID)
##    Data: data_long
## 
## REML criterion at convergence: 193.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8798 -0.5037  0.0005  0.5458  3.5407 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 5.87880  2.4246  
##  Residual             0.01334  0.1155  
## Number of obs: 360, groups:  ID, 90
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       12.875278   0.255650  86.999999  50.363  < 2e-16 ***
## Kelompok1          0.270556   0.361544  87.000014   0.748    0.456    
## Kelompok2         -0.273611   0.361544  87.000014  -0.757    0.451    
## Waktu1            -0.398611   0.010544 261.000000 -37.805  < 2e-16 ***
## Waktu2            -0.133056   0.010544 261.000000 -12.619  < 2e-16 ***
## Waktu3             0.125833   0.010544 261.000000  11.934  < 2e-16 ***
## Kelompok1:Waktu1   0.226111   0.014911 261.000000  15.164  < 2e-16 ***
## Kelompok2:Waktu1  -0.066389   0.014911 261.000000  -4.452 1.26e-05 ***
## Kelompok1:Waktu2   0.077222   0.014911 261.000000   5.179 4.47e-07 ***
## Kelompok2:Waktu2  -0.015278   0.014911 261.000000  -1.025    0.307    
## Kelompok1:Waktu3  -0.065000   0.014911 261.000000  -4.359 1.88e-05 ***
## Kelompok2:Waktu3   0.005833   0.014911 261.000000   0.391    0.696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Klmpk1 Klmpk2 Waktu1 Waktu2 Waktu3 Kl1:W1 Kl2:W1 Kl1:W2
## Kelompok1    0.000                                                        
## Kelompok2    0.000 -0.500                                                 
## Waktu1       0.000  0.000  0.000                                          
## Waktu2       0.000  0.000  0.000 -0.333                                   
## Waktu3       0.000  0.000  0.000 -0.333 -0.333                            
## Klmpk1:Wkt1  0.000  0.000  0.000  0.000  0.000  0.000                     
## Klmpk2:Wkt1  0.000  0.000  0.000  0.000  0.000  0.000 -0.500              
## Klmpk1:Wkt2  0.000  0.000  0.000  0.000  0.000  0.000 -0.333  0.167       
## Klmpk2:Wkt2  0.000  0.000  0.000  0.000  0.000  0.000  0.167 -0.333 -0.500
## Klmpk1:Wkt3  0.000  0.000  0.000  0.000  0.000  0.000 -0.333  0.167 -0.333
## Klmpk2:Wkt3  0.000  0.000  0.000  0.000  0.000  0.000  0.167 -0.333  0.167
##             Kl2:W2 Kl1:W3
## Kelompok1                
## Kelompok2                
## Waktu1                   
## Waktu2                   
## Waktu3                   
## Klmpk1:Wkt1              
## Klmpk2:Wkt1              
## Klmpk1:Wkt2              
## Klmpk2:Wkt2              
## Klmpk1:Wkt3  0.167       
## Klmpk2:Wkt3 -0.333 -0.500

# Random intercept + random slope
lmm2 <- lmer(
  BB_kg ~ Kelompok * Waktu + (1 + Waktu_num | ID),
  data = data_long, REML = TRUE,
  control = lmerControl(optimizer = "bobyqa")
)
summary(lmm2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BB_kg ~ Kelompok * Waktu + (1 + Waktu_num | ID)
##    Data: data_long
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 136.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.34694 -0.45022 -0.01438  0.48019  2.26498 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID       (Intercept) 5.803933 2.40914       
##           Waktu_num   0.003834 0.06192  0.15 
##  Residual             0.006951 0.08337       
## Number of obs: 360, groups:  ID, 90
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       12.875278   0.255652  86.998052  50.362  < 2e-16 ***
## Kelompok1          0.270556   0.361547  86.998151   0.748 0.456281    
## Kelompok2         -0.273611   0.361547  86.998151  -0.757 0.451227    
## Waktu1            -0.398611   0.012400 118.737351 -32.145  < 2e-16 ***
## Waktu2            -0.133056   0.008281 244.689332 -16.068  < 2e-16 ***
## Waktu3             0.125833   0.008281 244.689332  15.196  < 2e-16 ***
## Kelompok1:Waktu1   0.226111   0.017537 118.737353  12.893  < 2e-16 ***
## Kelompok2:Waktu1  -0.066389   0.017537 118.737353  -3.786 0.000242 ***
## Kelompok1:Waktu2   0.077222   0.011711 244.689332   6.594 2.62e-10 ***
## Kelompok2:Waktu2  -0.015278   0.011711 244.689332  -1.305 0.193262    
## Kelompok1:Waktu3  -0.065000   0.011711 244.689332  -5.550 7.40e-08 ***
## Kelompok2:Waktu3   0.005833   0.011711 244.689332   0.498 0.618853    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Klmpk1 Klmpk2 Waktu1 Waktu2 Waktu3 Kl1:W1 Kl2:W1 Kl1:W2
## Kelompok1    0.000                                                        
## Kelompok2    0.000 -0.500                                                 
## Waktu1      -0.149  0.000  0.000                                          
## Waktu2      -0.075  0.000  0.000  0.123                                   
## Waktu3       0.075  0.000  0.000 -0.499 -0.437                            
## Klmpk1:Wkt1  0.000 -0.149  0.075  0.000  0.000  0.000                     
## Klmpk2:Wkt1  0.000  0.075 -0.149  0.000  0.000  0.000 -0.500              
## Klmpk1:Wkt2  0.000 -0.075  0.037  0.000  0.000  0.000  0.123 -0.062       
## Klmpk2:Wkt2  0.000  0.037 -0.075  0.000  0.000  0.000 -0.062  0.123 -0.500
## Klmpk1:Wkt3  0.000  0.075 -0.037  0.000  0.000  0.000 -0.499  0.250 -0.437
## Klmpk2:Wkt3  0.000 -0.037  0.075  0.000  0.000  0.000  0.250 -0.499  0.218
##             Kl2:W2 Kl1:W3
## Kelompok1                
## Kelompok2                
## Waktu1                   
## Waktu2                   
## Waktu3                   
## Klmpk1:Wkt1              
## Klmpk2:Wkt1              
## Klmpk1:Wkt2              
## Klmpk2:Wkt2              
## Klmpk1:Wkt3  0.218       
## Klmpk2:Wkt3 -0.437 -0.500

# Bandingkan struktur random effect
anova(lmm1, lmm2, refit = FALSE)

# Uji efek tetap
anova(lmm2, ddf = "Kenward-Roger")

# ICC
performance::icc(lmm1)

Interpretasi.

Model dengan random slope lebih baik daripada random intercept saja (χ²(2) = 57,3; p < 0,001; AIC 168,3 vs 221,6), artinya laju kenaikan berat badan memang bervariasi antarbalita.
Uji efek tetap (Kenward-Roger) konsisten dengan mixed ANOVA: kelompok p = 0,687; waktu F(3; 185,4) = 408,58; interaksi F(6; 206,4) = 36,12 (keduanya p < 0,001).
ICC = 0,972 (tanpa penyesuaian) dan 0,998 (disesuaikan): hampir seluruh variasi berat badan berasal dari perbedaan antarbalita, bukan dari pengukuran ulang.

8.1 LMM dengan usia dan jenis kelamin sebagai kovariat

lmm3 <- lmer(
  BB_kg ~ Kelompok * Waktu_num + Usia_bulan + Jenis_Kelamin + (1 | ID),
  data = data_long, REML = TRUE
)
anova(lmm3, ddf = "Kenward-Roger")

summary(lmm3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BB_kg ~ Kelompok * Waktu_num + Usia_bulan + Jenis_Kelamin + (1 |  
##     ID)
##    Data: data_long
## 
## REML criterion at convergence: 1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8760 -0.5403 -0.0130  0.5250  3.6046 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.88462  0.9405  
##  Residual             0.01315  0.1147  
## Number of obs: 360, groups:  ID, 90
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           6.429612   0.290550  85.132562  22.129  < 2e-16 ***
## Kelompok1            -0.175643   0.147694  86.033282  -1.189    0.238    
## Kelompok2             0.110239   0.143289  86.098344   0.769    0.444    
## Waktu_num             0.267222   0.005406 267.000000  49.430  < 2e-16 ***
## Usia_bulan            0.161073   0.007267  85.000002  22.166  < 2e-16 ***
## Jenis_Kelamin1        0.137981   0.102668  85.000000   1.344    0.183    
## Kelompok1:Waktu_num  -0.153556   0.007645 267.000000 -20.085  < 2e-16 ***
## Kelompok2:Waktu_num   0.044778   0.007645 267.000000   5.857 1.38e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Klmpk1 Klmpk2 Wkt_nm Us_bln Jns_K1 Kl1:W_
## Kelompok1    0.186                                          
## Kelompok2   -0.129 -0.523                                   
## Waktu_num   -0.028  0.000  0.000                            
## Usia_bulan  -0.939 -0.207  0.142  0.000                     
## Jenis_Klmn1 -0.095  0.203 -0.103  0.000  0.059              
## Klmpk1:Wkt_  0.000 -0.078  0.040  0.000  0.000  0.000       
## Klmpk2:Wkt_  0.000  0.039 -0.080  0.000  0.000  0.000 -0.500

# Laju kenaikan BB per bulan (kg/bulan) tiap kelompok dan perbandingannya
emtrends(lmm3, pairwise ~ Kelompok, var = "Waktu_num", adjust = "tukey")

## $emtrends
##  Kelompok    Waktu_num.trend      SE  df lower.CL upper.CL
##  Kontrol               0.114 0.00936 267   0.0952    0.132
##  PMT Biskuit           0.312 0.00936 267   0.2936    0.330
##  PMT Lokal             0.376 0.00936 267   0.3576    0.394
## 
## Results are averaged over the levels of: Jenis_Kelamin 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate     SE  df t.ratio p.value
##  Kontrol - PMT Biskuit     -0.198 0.0132 267 -14.977 <0.0001
##  Kontrol - PMT Lokal       -0.262 0.0132 267 -19.810 <0.0001
##  PMT Biskuit - PMT Lokal   -0.064 0.0132 267  -4.833 <0.0001
## 
## Results are averaged over the levels of: Jenis_Kelamin 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

Interpretasi.

Usia berpengaruh sangat kuat terhadap berat badan (F(1; 85) = 491,35; p < 0,001), sedangkan jenis kelamin tidak (p = 0,183).
Interaksi kelompok × waktu tetap signifikan setelah usia dan jenis kelamin dikontrol (F(2; 267) = 213,38; p < 0,001).
Laju kenaikan berat badan: Kontrol 0,114 kg/bulan; PMT Biskuit 0,312 kg/bulan; PMT Lokal 0,376 kg/bulan. Semua pasangan berbeda bermakna (Tukey p < 0,0001).

8.2 Diagnostik residual

par(mfrow = c(1, 3))

qqnorm(resid(lmm2), main = "Q-Q residual")
qqline(resid(lmm2))

qqnorm(ranef(lmm2)$ID[, 1], main = "Q-Q random intercept")
qqline(ranef(lmm2)$ID[, 1])

plot(fitted(lmm2), resid(lmm2),
     xlab = "Nilai prediksi", ylab = "Residual", main = "Residual vs prediksi")
abline(h = 0, lty = 2)

par(mfrow = c(1, 1))

shapiro.test(resid(lmm2))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(lmm2)
## W = 0.9941, p-value = 0.177

Residual berdistribusi normal (Shapiro-Wilk W = 0,994; p = 0,177) dan titik-titik pada Q-Q plot mengikuti garis lurus, sehingga asumsi LMM terpenuhi.

9. Opsional: Simulasi Data Hilang untuk Demonstrasi Keunggulan LMM

Bagian ini tidak mengubah data utama. Sebanyak 30 nilai (di luar bulan 0) dihapus secara acak untuk menunjukkan bahwa LMM tetap dapat dipakai.

set.seed(2026)

data_miss <- data_long
idx_miss <- sample(which(data_miss$Waktu != "B0"), 30, replace = FALSE)
data_miss$BB_kg[idx_miss] <- NA

lmm_miss <- lmer(
  BB_kg ~ Kelompok * Waktu + (1 + Waktu_num | ID),
  data = data_miss, REML = TRUE, na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa")
)
anova(lmm_miss, ddf = "Kenward-Roger")

# Jumlah balita yang kehilangan >= 1 nilai (akan dibuang seluruhnya oleh RM ANOVA)
n_distinct(data_miss$ID[is.na(data_miss$BB_kg)])

## [1] 28

Dengan 30 nilai hilang, 28 dari 90 balita (31%) kehilangan sedikitnya satu penimbangan. RM ANOVA akan membuang seluruh 28 balita tersebut, sedangkan LMM tetap memakai semua data yang ada. Hasil LMM tetap konsisten dengan analisis lengkap (interaksi F(6; 185,2) = 32,43; p < 0,001).

10. Kesimpulan

Berat badan balita naik bermakna selama tiga bulan pada ketiga kelompok (efek waktu, p < 0,001).
Pola kenaikan berbeda antarkelompok (interaksi kelompok × waktu, F = 70,39; p < 0,001; η² parsial = 0,618), sehingga H0 ditolak.
Urutan kenaikan berat badan: PMT Lokal (1,13 kg) > PMT Biskuit (0,95 kg) > Kontrol (0,34 kg), dengan semua perbedaan signifikan setelah koreksi Holm atau Tukey.
Berat badan absolut antarkelompok tidak berbeda pada tiap waktu; perbedaan terletak pada laju kenaikan (0,376; 0,312; dan 0,114 kg per bulan).
Hasil kokoh terhadap pilihan metode: mixed ANOVA dengan koreksi GG, kontras perubahan, ANOVA selisih, ANCOVA, dan LMM dengan kovariat usia serta jenis kelamin memberi kesimpulan yang sama.

Catatan dan keterbatasan

Sferisitas dilanggar (Mauchly p < 0,001), jadi semua uji within memakai koreksi Greenhouse-Geisser; uji lanjut memakai model multivariat.
Usia rata-rata kelompok Kontrol lebih tinggi; walaupun belum berbeda bermakna (p = 0,098), pengaruh usia terhadap berat badan sangat kuat, sehingga hasil terkoreksi (ANCOVA dan LMM) sebaiknya turut dilaporkan.
Korelasi antarwaktu yang sangat tinggi (ICC 0,97) membuat η² parsial tampak sangat besar; laporkan juga ukuran efek dalam satuan kg (selisih rerata) agar mudah diinterpretasi.

Kesimpulan Analisis

Analisis mixed ANOVA menunjukkan interaksi kelompok × waktu yang signifikan setelah koreksi Greenhouse-Geisser, F(3,96; 172,38) = 70,39; p < 0,001; η² parsial = 0,618, yang berarti pola perubahan berat badan selama tiga bulan berbeda antarkelompok. Kenaikan berat badan tertinggi terdapat pada kelompok PMT Lokal (1,13 kg), diikuti PMT Biskuit (0,95 kg) dan Kontrol (0,34 kg); ketiga kelompok berbeda bermakna (uji kontras, koreksi Holm, p ≤ 0,0014).

11. Menyimpan Hasil

Chunk ini menyimpan tabel dan grafik ke folder hasil_laporan (dibuat otomatis) di samping file .Rmd.

dir.create("hasil_laporan", showWarnings = FALSE)
f <- function(x) file.path("hasil_laporan", x)

# Data format LONG dan WIDE (CSV)
write.csv(data_long, f("dataset_PMT_balita_LONG.csv"), row.names = FALSE)
write.csv(data_wide, f("dataset_PMT_balita_WIDE.csv"), row.names = FALSE)

# Tabel hasil
write.csv(deskriptif, f("hasil_01_statistik_deskriptif_BB.csv"), row.names = FALSE)
write.csv(shapiro_mixed, f("hasil_02_shapiro_Kelompok_x_Waktu.csv"), row.names = FALSE)
write.csv(levene_mixed, f("hasil_03_levene_per_Waktu.csv"), row.names = FALSE)

hasil_anova <- as.data.frame(nice(aov2, correction = "GG", es = c("ges", "pes")))
write.csv(hasil_anova, f("hasil_04_mixed_ANOVA_GG.csv"), row.names = FALSE)

posthoc_group <- as.data.frame(pairs(em2b, adjust = "holm"))
write.csv(posthoc_group, f("hasil_05_posthoc_antar_kelompok_Holm.csv"), row.names = FALSE)

posthoc_time <- as.data.frame(pairs(em2, adjust = "holm"))
write.csv(posthoc_time, f("hasil_06_posthoc_dalam_kelompok_Holm.csv"), row.names = FALSE)

emm_group_time <- as.data.frame(em2b)
write.csv(emm_group_time, f("hasil_07_estimated_marginal_means.csv"), row.names = FALSE)

# Grafik
ggsave(f("grafik_01_profile_BB.png"), p_profil, width = 8, height = 5.5, dpi = 300)
ggsave(f("grafik_02_spaghetti_BB.png"), p_spag, width = 9, height = 6, dpi = 300)
ggsave(f("grafik_03_interaksi_Kelompok_x_Waktu.png"), p_interaksi, width = 8, height = 5.5, dpi = 300)

cat("Jumlah balita    :", n_distinct(data_long$ID), "\n")

## Jumlah balita    : 90

cat("Kelompok         :", n_distinct(data_long$Kelompok), "\n")

## Kelompok         : 3

cat("Waktu penimbangan:", n_distinct(data_long$Waktu), "\n")

## Waktu penimbangan: 4

cat("Outcome          : BB_kg\n")

## Outcome          : BB_kg

cat("File hasil tersimpan di:", file.path(getwd(), "hasil_laporan"), "\n")

## File hasil tersimpan di: C:/Users/ASUS/Downloads/hasil_laporan

Informasi Sesi

sessionInfo()

## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Singapore
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggpubr_1.0.0       performance_0.18.2 lmerTest_3.2-1     effectsize_1.0.3  
##  [5] car_3.1-5          carData_3.0-6      rstatix_1.1.0      emmeans_2.0.4     
##  [9] afex_1.5-1         lme4_2.0-6         Matrix_1.7-5       ggplot2_4.0.3     
## [13] tidyr_1.3.2        dplyr_1.2.1        readxl_1.5.0.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   digest_0.6.39      
##  [7] rpart_4.1.27        estimability_2.0.0  lifecycle_1.0.5    
## [10] cluster_2.1.8.2     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         withr_3.0.3         foreign_0.8-91     
## [28] purrr_1.2.2         numDeriv_2016.8-1.1 nnet_7.3-20        
## [31] grid_4.6.1          datawizard_1.4.0    colorspace_2.1-3   
## [34] scales_1.4.0        MASS_7.3-65         insight_1.5.4      
## [37] cli_3.6.6           mvtnorm_1.4-2       rmarkdown_2.32     
## [40] ragg_1.5.2          reformulas_0.4.4    generics_0.1.4     
## [43] otel_0.2.0          rstudioapi_0.19.0   reshape2_1.4.5     
## [46] parameters_0.29.3   minqa_1.2.8         cachem_1.1.0       
## [49] stringr_1.6.0       splines_4.6.1       parallel_4.6.1     
## [52] cellranger_1.1.0    base64enc_0.1-6     vctrs_0.7.3        
## [55] boot_1.3-32         jsonlite_2.0.0      pbkrtest_0.5.5     
## [58] htmlTable_2.5.0     Formula_1.2-6       systemfonts_1.3.2  
## [61] jquerylib_0.1.4     glue_1.8.1          nloptr_2.2.1       
## [64] stringi_1.8.9       gtable_0.3.6        tibble_3.3.1       
## [67] pillar_1.11.1       htmltools_0.5.9     R6_2.6.1           
## [70] textshaping_1.0.5   Rdpack_2.6.6        evaluate_1.0.5     
## [73] lattice_0.22-9      rbibutils_2.4.1     backports_1.5.1    
## [76] broom_1.0.13        bslib_0.12.0        Rcpp_1.1.2         
## [79] checkmate_2.3.4     gridExtra_2.3.1     nlme_3.1-169       
## [82] xfun_0.61           pkgconfig_2.0.3

Laporan Analisis Repeated Measures: Pengaruh PMT Biskuit dan PMT Lokal terhadap Berat Badan Balita

Sundari (NIM 2611018024)

03 October 2026

1. Pendahuluan

1.1 Latar belakang dan tujuan

1.2 Desain dan variabel

1.3 Hipotesis uji utama

2. Ringkasan Temuan

3. Persiapan Paket dan Data

3.1 Paket

3.2 Membaca data

3.3 Validasi struktur data

3.4 Pengkodean ulang variabel

3.5 Konsistensi format LONG dan WIDE

3.6 Karakteristik awal: apakah kelompok setara?

4. Eksplorasi Data

4.1 Statistik deskriptif

4.2 Matriks kovarians, korelasi, dan varians selisih antarwaktu

4.3 Profile plot

4.4 Spaghetti plot

5. RM ANOVA Satu Arah: Contoh Kelompok PMT Lokal

5.1 Uji asumsi

5.2 RM ANOVA dengan afex

5.3 Pendekatan multivariat

5.4 Uji lanjut dan tren

5.5 Alternatif nonparametrik: Friedman

6. Mixed Design ANOVA: Kelompok × Waktu

6.1 Uji asumsi

6.2 Mixed ANOVA

6.3 Ukuran efek

6.4 Plot interaksi

6.5 Simple effects

6.6 Uji lanjut dalam kelompok

6.7 Uji lanjut antarkelompok pada tiap waktu

6.8 Kontras perubahan bulan 0 → bulan 3

6.9 Tren linear antarkelompok

7. Pembanding: Selisih Berat Badan dan ANCOVA

7.1 ANCOVA: dikontrol berat badan awal, usia, dan jenis kelamin

8. Pembanding: Linear Mixed Model (LMM)

8.1 LMM dengan usia dan jenis kelamin sebagai kovariat

8.2 Diagnostik residual

9. Opsional: Simulasi Data Hilang untuk Demonstrasi Keunggulan LMM

10. Kesimpulan

Catatan dan keterbatasan

Kesimpulan Analisis

11. Menyimpan Hasil

Informasi Sesi