# =============================================================================
# Nama : Dian Wahyu Andini Putri
# NIM : 2611018016
# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan:
# Analisis perubahan kadar hemoglobin (Hb) pada remaja dengan anemia
# berdasarkan kelompok intervensi selama 12 minggu
#
# DATA SUMBER:
# data_hb_anemia_remaja.xlsx
#
# Struktur analisis mengikuti script contoh:
# 0. Paket & pengaturan
# 1. Membaca data Excel dan menyiapkan format panjang & lebar
# 2. Eksplorasi data: statistik deskriptif, profile plot, spaghetti plot
# 3. Repeated Measure ANOVA satu arah (within-subject: Waktu)
# 3a. Uji asumsi: outlier, normalitas, sferisitas (Mauchly)
# 3b. ANOVA + koreksi Greenhouse-Geisser
# 3c. Pendekatan multivariat (MANOVA)
# 3d. Post hoc berpasangan & kontras polinomial
# 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. Efek sederhana & post hoc
# 4d. Kontras interaksi (perubahan Hb dari baseline antarkelompok)
# 5. Pembanding: Linear Mixed Model (LMM)
# 6. Menyimpan data & ringkasan hasil
# =============================================================================
# 0. PAKET & PENGATURAN
# Jalankan sekali bila paket belum terpasang:
# install.packages(c("readxl", "writexl", "tidyverse", "afex", "emmeans",
# "rstatix", "car", "effectsize", "lme4", "lmerTest",
# "performance", "ggpubr"))
suppressPackageStartupMessages({
library(readxl) # membaca Excel
library(writexl) # menyimpan Excel
library(dplyr) # manipulasi data
library(tidyr) # format panjang <-> lebar
library(ggplot2) # grafik
library(afex) # ANOVA within/mixed
library(emmeans) # rerata marginal, post hoc, kontras
library(rstatix) # uji asumsi dan uji nonparametrik
library(car) # leveneTest / Anova
library(effectsize) # ukuran efek
library(lme4) # linear mixed model
library(lmerTest) # uji F/t pada LMM
library(performance) # ICC dan diagnostik
library(ggpubr) # Q-Q plot
})
options(contrasts = c("contr.sum", "contr.poly"))
afex_options(emmeans_model = "multivariate")
theme_set(theme_bw(base_size = 12))
# 1. MEMBACA DATA EXCEL & MENYIAPKAN DATA
# Jika file berada di working directory R:
file_excel <- "data_hb_anemia_remaja.xlsx"
# Jika R tidak menemukan file, gunakan path lengkap, contoh:
# file_excel <- "C:/Users/Nama/Documents/data_hb_anemia_remaja.xlsx"
# Lihat nama sheet
excel_sheets(file_excel)
## [1] "Keterangan" "Data_wide" "Data_long" "Ringkasan"
# Pilih sheet data utama
# Script menggunakan sheet "Data_wide" sebagai sumber utama.
# Data_wide berisi satu baris untuk satu responden dengan pengukuran Hb
# pada minggu 0, 4, 8, dan 12.
dat_wide <- read_excel(file_excel, sheet = "Data_wide")
# Lihat struktur awal data
head(dat_wide)
## # A tibble: 6 × 7
## id kelompok usia Hb_W0 Hb_W4 Hb_W8 Hb_W12
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 R001 Kontrol 13 10.3 9.8 10 10.4
## 2 R002 Kontrol 16 10.7 10.8 10.4 11.4
## 3 R003 Kontrol 16 9.9 9.9 9.3 9.2
## 4 R004 Kontrol 18 10.2 9.3 9.2 9
## 5 R005 Kontrol 17 10.5 10.6 10.1 10.3
## 6 R006 Kontrol 13 9.2 10.3 10.5 9.9
str(dat_wide)
## tibble [90 × 7] (S3: tbl_df/tbl/data.frame)
## $ id : chr [1:90] "R001" "R002" "R003" "R004" ...
## $ kelompok: chr [1:90] "Kontrol" "Kontrol" "Kontrol" "Kontrol" ...
## $ usia : num [1:90] 13 16 16 18 17 13 12 13 12 14 ...
## $ Hb_W0 : num [1:90] 10.3 10.7 9.9 10.2 10.5 9.2 10 9.6 10.5 11.1 ...
## $ Hb_W4 : num [1:90] 9.8 10.8 9.9 9.3 10.6 10.3 9.5 10 10.3 10.8 ...
## $ Hb_W8 : num [1:90] 10 10.4 9.3 9.2 10.1 10.5 9.5 10.1 10.4 11.1 ...
## $ Hb_W12 : num [1:90] 10.4 11.4 9.2 9 10.3 9.9 9.7 10.1 10.6 10.5 ...
names(dat_wide)
## [1] "id" "kelompok" "usia" "Hb_W0" "Hb_W4" "Hb_W8" "Hb_W12"
# Jika kolom kelompok memiliki nama berbeda, sesuaikan bagian ini.
# Berdasarkan struktur data, variabel yang digunakan:
# id, kelompok, Hb_W0, Hb_W4, Hb_W8, Hb_W12
# Pastikan ID menjadi faktor
dat_wide$id <- factor(dat_wide$id)
# Pastikan kelompok menjadi faktor.
# Urutan level mengikuti urutan yang terdapat pada data.
dat_wide$kelompok <- factor(dat_wide$kelompok)
# Pastikan variabel Hb bersifat numerik
hb_cols <- c("Hb_W0", "Hb_W4", "Hb_W8", "Hb_W12")
dat_wide[hb_cols] <- lapply(dat_wide[hb_cols], as.numeric)
# Waktu pengukuran
minggu <- c(0, 4, 8, 12)
# Format panjang
# Satu baris = satu pengukuran Hb pada satu responden.
dat_long <- dat_wide |>
pivot_longer(
cols = all_of(hb_cols),
names_to = "waktu",
values_to = "Hb"
) |>
mutate(
waktu = factor(
waktu,
levels = hb_cols,
labels = paste0("M", minggu)
),
minggu = as.numeric(sub("M", "", as.character(waktu)))
)
# Periksa data
head(dat_wide)
## # A tibble: 6 × 7
## id kelompok usia Hb_W0 Hb_W4 Hb_W8 Hb_W12
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 R001 Kontrol 13 10.3 9.8 10 10.4
## 2 R002 Kontrol 16 10.7 10.8 10.4 11.4
## 3 R003 Kontrol 16 9.9 9.9 9.3 9.2
## 4 R004 Kontrol 18 10.2 9.3 9.2 9
## 5 R005 Kontrol 17 10.5 10.6 10.1 10.3
## 6 R006 Kontrol 13 9.2 10.3 10.5 9.9
head(dat_long)
## # A tibble: 6 × 6
## id kelompok usia waktu Hb minggu
## <fct> <fct> <dbl> <fct> <dbl> <dbl>
## 1 R001 Kontrol 13 M0 10.3 0
## 2 R001 Kontrol 13 M4 9.8 4
## 3 R001 Kontrol 13 M8 10 8
## 4 R001 Kontrol 13 M12 10.4 12
## 5 R002 Kontrol 16 M0 10.7 0
## 6 R002 Kontrol 16 M4 10.8 4
str(dat_long)
## tibble [360 × 6] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 90 levels "R001","R002",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok: Factor w/ 3 levels "Kontrol","TTD",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num [1:360] 13 13 13 13 16 16 16 16 16 16 ...
## $ waktu : Factor w/ 4 levels "M0","M4","M8",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Hb : num [1:360] 10.3 9.8 10 10.4 10.7 10.8 10.4 11.4 9.9 9.9 ...
## $ minggu : num [1:360] 0 4 8 12 0 4 8 12 0 4 ...
# Jumlah responden
n_distinct(dat_wide$id)
## [1] 90
# Jumlah responden menurut kelompok
table(dat_wide$kelompok)
##
## Kontrol TTD TTD+VitC
## 30 30 30
# Jumlah observasi menurut kelompok dan waktu
table(dat_long$kelompok, dat_long$waktu)
##
## M0 M4 M8 M12
## Kontrol 30 30 30 30
## TTD 30 30 30 30
## TTD+VitC 30 30 30 30
# 2. EKSPLORASI DATA
# 2.1 Statistik deskriptif
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(Hb, type = "mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 Hb 30 10.6 0.625
## 2 Kontrol M4 Hb 30 10.8 0.772
## 3 Kontrol M8 Hb 30 10.8 0.879
## 4 Kontrol M12 Hb 30 10.8 0.916
## 5 TTD M0 Hb 30 10.6 0.625
## 6 TTD M4 Hb 30 11.1 0.637
## 7 TTD M8 Hb 30 11.5 0.826
## 8 TTD M12 Hb 30 11.8 0.83
## 9 TTD+VitC M0 Hb 30 10.6 0.636
## 10 TTD+VitC M4 Hb 30 11.3 0.646
## 11 TTD+VitC M8 Hb 30 11.8 0.694
## 12 TTD+VitC M12 Hb 30 12.1 0.746
# Statistik yang lebih lengkap
desk_lengkap <- dat_long |>
group_by(kelompok, waktu) |>
summarise(
n = sum(!is.na(Hb)),
mean = mean(Hb, na.rm = TRUE),
sd = sd(Hb, na.rm = TRUE),
median = median(Hb, na.rm = TRUE),
min = min(Hb, na.rm = TRUE),
max = max(Hb, na.rm = TRUE),
.groups = "drop"
)
desk_lengkap
## # A tibble: 12 × 8
## kelompok waktu n mean sd median min max
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kontrol M0 30 10.6 0.625 10.7 9.2 12.3
## 2 Kontrol M4 30 10.8 0.772 10.8 9.3 12.8
## 3 Kontrol M8 30 10.8 0.879 10.5 9.2 13.4
## 4 Kontrol M12 30 10.8 0.916 10.9 9 13.5
## 5 TTD M0 30 10.6 0.625 10.8 9.5 12.1
## 6 TTD M4 30 11.1 0.637 11.3 9.8 12.5
## 7 TTD M8 30 11.5 0.826 11.6 9.9 13
## 8 TTD M12 30 11.8 0.830 11.9 10.2 12.9
## 9 TTD+VitC M0 30 10.6 0.636 10.8 9.4 11.6
## 10 TTD+VitC M4 30 11.3 0.646 11.4 9.5 12.4
## 11 TTD+VitC M8 30 11.8 0.694 11.9 10.1 12.8
## 12 TTD+VitC M12 30 12.1 0.746 12.1 10.4 13.4
# 2.2 Matriks kovarians dan korelasi antarwaktu
S <- cov(
dat_wide[, hb_cols],
use = "pairwise.complete.obs"
)
R <- cor(
dat_wide[, hb_cols],
use = "pairwise.complete.obs"
)
round(S, 3)
## Hb_W0 Hb_W4 Hb_W8 Hb_W12
## Hb_W0 0.387 0.328 0.387 0.400
## Hb_W4 0.328 0.510 0.549 0.591
## Hb_W8 0.387 0.549 0.817 0.792
## Hb_W12 0.400 0.591 0.792 0.970
round(R, 3)
## Hb_W0 Hb_W4 Hb_W8 Hb_W12
## Hb_W0 1.000 0.739 0.689 0.653
## Hb_W4 0.739 1.000 0.850 0.840
## Hb_W8 0.689 0.850 1.000 0.890
## Hb_W12 0.653 0.840 0.890 1.000
# 2.3 Varians selisih antarpasangan waktu
# Berkaitan dengan asumsi sferisitas
pasangan <- combn(hb_cols, 2)
var_selisih <- apply(
pasangan,
2,
function(p) {
x <- dat_wide[[p[1]]]
y <- dat_wide[[p[2]]]
var(x - y, na.rm = TRUE)
}
)
names(var_selisih) <- apply(
pasangan,
2,
paste,
collapse = " - "
)
round(var_selisih, 3)
## Hb_W0 - Hb_W4 Hb_W0 - Hb_W8 Hb_W0 - Hb_W12 Hb_W4 - Hb_W8 Hb_W4 - Hb_W12
## 0.240 0.429 0.557 0.229 0.299
## Hb_W8 - Hb_W12
## 0.203
# 2.4 Profile plot
# Rerata Hb +/- 95% CI
p_profil <- ggplot(
dat_long,
aes(
x = minggu,
y = Hb,
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 = .6
) +
scale_x_continuous(breaks = minggu) +
labs(
x = "Minggu ke-",
y = "Kadar hemoglobin (g/dL)",
colour = "Kelompok",
title = "Profil rerata kadar Hb (± 95% CI)"
) +
theme(
legend.position = "bottom"
)
p_profil
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
# 2.5 Spaghetti plot
# Lintasan kadar Hb setiap responden
p_spag <- ggplot(
dat_long,
aes(
x = minggu,
y = Hb,
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 = minggu) +
labs(
x = "Minggu ke-",
y = "Kadar Hb (g/dL)",
title = "Lintasan individu dan rerata kadar Hb menurut kelompok"
)
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan:
# Apakah kadar Hb berubah selama 12 minggu pada kelompok intervensi
# yang dipilih?
#
# Pada script contoh, analisis satu arah dilakukan pada satu kelompok.
# Di sini kita menggunakan kelompok intervensi "TTD+VitC".
#
# Jika nama kelompok pada Excel berbeda, cek:
# levels(dat_wide$kelompok)
# lalu ubah nilai berikut.
kelompok_rm <- "TTD+VitC"
# Periksa apakah kelompok tersedia
if (!kelompok_rm %in% levels(dat_wide$kelompok)) {
stop(
paste0(
"Kelompok '", kelompok_rm,
"' tidak ditemukan. Jalankan levels(dat_wide$kelompok) ",
"untuk melihat nama kelompok yang tersedia."
)
)
}
d1 <- droplevels(
filter(dat_long, kelompok == kelompok_rm)
)
d1w <- filter(
dat_wide,
kelompok == kelompok_rm
)
# 3a. UJI ASUMSI
# (i) Outlier per waktu
d1 |>
group_by(waktu) |>
identify_outliers(Hb)
## [1] waktu id kelompok usia Hb minggu is.outlier
## [8] is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas per waktu: Shapiro-Wilk
d1 |>
group_by(waktu) |>
shapiro_test(Hb)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 Hb 0.941 0.0976
## 2 M4 Hb 0.966 0.439
## 3 M8 Hb 0.958 0.278
## 4 M12 Hb 0.977 0.756
# Q-Q plot
ggqqplot(
d1,
"Hb",
facet.by = "waktu"
)
# (iii) Sferisitas: Mauchly
# Dilaporkan otomatis oleh anova_test dan afex
aov1_rs <- anova_test(
data = d1,
dv = Hb,
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 143.204 1.52e-33 * 0.832
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.948 0.917
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.964 2.89, 83.86 2.03e-32 * 1.082 3.25, 94.14 1.52e-33
## p[HF]<.05
## 1 *
# Tabel ANOVA dengan koreksi otomatis
# Jika Mauchly p < 0.05, Greenhouse-Geisser digunakan.
get_anova_table(
aov1_rs,
correction = "auto"
)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 87 143.204 1.52e-33 * 0.832
# 3b. REPEATED MEASURE ANOVA DENGAN AFEX
aov1 <- aov_ez(
id = "id",
dv = "Hb",
data = d1,
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov1
## Anova Table (Type 3 tests)
##
## Response: Hb
## Effect df MSE F ges pes p.value
## 1 waktu 2.89, 83.86 0.09 143.20 *** .416 .832 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Ringkasan:
# Mauchly, epsilon GG/HF, dan hasil ANOVA
summary(aov1)
## Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): HF eps > 1
## treated as 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 15720.9 1 46.140 29 9880.8 < 2.2e-16 ***
## waktu 38.4 3 7.783 87 143.2 < 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.94836 0.91659
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.96388 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 1.082018 1.523274e-33
# Ukuran efek
eta_squared(
aov1,
partial = TRUE
)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.83 | [0.78, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
omega_squared(
aov1,
partial = TRUE
)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## -------------------------------------------
## waktu | 0.41 | [0.27, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 3c. PENDEKATAN MULTIVARIAT
# Pendekatan ini tidak memerlukan asumsi sferisitas.
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.99707 9880.8 1 29 < 2.2e-16 ***
## waktu 1 0.93980 140.5 3 27 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3d. POST HOC & KONTRAS TREN
# Estimated marginal means
em1 <- emmeans(
aov1,
~ waktu
)
em1
## waktu emmean SE df lower.CL upper.CL
## M0 10.6 0.116 29 10.4 10.8
## M4 11.3 0.118 29 11.0 11.5
## M8 11.8 0.127 29 11.5 12.0
## M12 12.1 0.136 29 11.8 12.4
##
## Confidence level used: 0.95
# Semua pasangan waktu
pairs(
em1,
adjust = "bonferroni"
)
## contrast estimate SE df t.ratio p.value
## M0 - M4 -0.670 0.0778 29 -8.614 <0.0001
## M0 - M8 -1.173 0.0732 29 -16.034 <0.0001
## M0 - M12 -1.500 0.0744 29 -20.152 <0.0001
## M4 - M8 -0.503 0.0771 29 -6.530 <0.0001
## M4 - M12 -0.830 0.0867 29 -9.571 <0.0001
## M8 - M12 -0.327 0.0733 29 -4.455 0.0007
##
## P value adjustment: bonferroni method for 6 tests
# Setiap waktu dibandingkan dengan baseline M0
contrast(
em1,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## contrast estimate SE df t.ratio p.value
## M4 - M0 0.67 0.0778 29 8.614 <0.0001
## M8 - M0 1.17 0.0732 29 16.034 <0.0001
## M12 - M0 1.50 0.0744 29 20.152 <0.0001
##
## P value adjustment: holm method for 3 tests
# Tren linear, kuadratik, dan kubik
contrast(
em1,
"poly"
)
## contrast estimate SE df t.ratio p.value
## linear 5.003 0.245 29 20.401 <0.0001
## quadratic -0.343 0.113 29 -3.032 0.0051
## cubic -0.010 0.234 29 -0.043 0.9662
# 3e. ALTERNATIF NONPARAMETRIK: FRIEDMAN
friedman_test(
d1,
Hb ~ waktu | id
)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 Hb 30 75.2 3 3.24e-16 Friedman test
friedman_effsize(
d1,
Hb ~ waktu | id
)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 Hb 30 0.836 Kendall W large
# Perbandingan berpasangan
d1 |>
wilcox_test(
Hb ~ 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 Hb M0 M4 30 30 5 0.0000000205 1.23e-7 ****
## 2 Hb M0 M8 30 30 0 0.00000000186 1.12e-8 ****
## 3 Hb M0 M12 30 30 0 0.00000000186 1.12e-8 ****
## 4 Hb M4 M8 30 30 24 0.00000180 1.08e-5 ****
## 5 Hb M4 M12 30 30 2 0.00000000559 3.35e-8 ****
## 6 Hb M8 M12 30 30 48 0.0000729 4.37e-4 ***
# 4. MIXED DESIGN ANOVA
# Kelompok [between] x Waktu [within]
#
# Pertanyaan utama:
# Apakah pola perubahan kadar Hb dari waktu ke waktu berbeda
# antar kelompok intervensi?
#
# Fokus utama biasanya adalah INTERAKSI:
# kelompok x waktu
#
# Jika interaksi signifikan, berarti perubahan Hb dari waktu ke waktu
# tidak sama antar kelompok.
# 4a. UJI ASUMSI
# (i) Outlier pada setiap kombinasi kelompok x waktu
dat_long |>
group_by(kelompok, waktu) |>
identify_outliers(Hb)
## # A tibble: 3 × 8
## kelompok waktu id usia Hb minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol M0 R022 16 12.3 0 TRUE FALSE
## 2 Kontrol M8 R022 16 13.4 8 TRUE FALSE
## 3 Kontrol M12 R022 16 13.5 12 TRUE FALSE
# (ii) Normalitas pada setiap sel
dat_long |>
group_by(kelompok, waktu) |>
shapiro_test(Hb)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 Hb 0.976 0.712
## 2 Kontrol M4 Hb 0.980 0.833
## 3 Kontrol M8 Hb 0.949 0.156
## 4 Kontrol M12 Hb 0.960 0.319
## 5 TTD M0 Hb 0.954 0.220
## 6 TTD M4 Hb 0.972 0.609
## 7 TTD M8 Hb 0.969 0.516
## 8 TTD M12 Hb 0.929 0.0461
## 9 TTD+VitC M0 Hb 0.941 0.0976
## 10 TTD+VitC M4 Hb 0.966 0.439
## 11 TTD+VitC M8 Hb 0.958 0.278
## 12 TTD+VitC M12 Hb 0.977 0.756
# Q-Q plot per kelompok dan waktu
ggqqplot(
dat_long,
"Hb",
ggtheme = theme_bw()
) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antar kelompok pada setiap waktu
dat_long |>
group_by(waktu) |>
levene_test(Hb ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 0.147 0.863
## 2 M4 2 87 0.413 0.663
## 3 M8 2 87 0.562 0.572
## 4 M12 2 87 0.402 0.670
# (iv) Homogenitas matriks kovarians antar kelompok
# Box's M
box_m(
dat_wide[, hb_cols],
dat_wide$kelompok
)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 15.9 0.723 20 Box's M-test for Homogeneity of Covariance Matric…
# (v) Sferisitas:
# Mauchly diperiksa melalui summary(aov2) di bawah.
# 4b. MIXED ANOVA
aov2 <- aov_ez(
id = "id",
dv = "Hb",
data = dat_long,
between = "kelompok",
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
aov2
## Anova Table (Type 3 tests)
##
## Response: Hb
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 1.88 8.51 *** .143 .164 <.001
## 2 waktu 2.82, 245.28 0.12 139.53 *** .194 .616 <.001
## 3 kelompok:waktu 5.64, 245.28 0.12 22.17 *** .071 .338 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Ringkasan model:
# termasuk Mauchly, epsilon GG/HF, dan p-value terkoreksi
summary(aov2)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 44778 1 163.339 87 23850.5989 < 2.2e-16 ***
## kelompok 32 2 163.339 87 8.5072 0.0004222 ***
## waktu 46 3 28.834 261 139.5322 < 2.2e-16 ***
## kelompok:waktu 15 6 28.834 261 22.1709 < 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.91392 0.17264
## kelompok:waktu 0.91392 0.17264
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.93976 < 2.2e-16 ***
## kelompok:waktu 0.93976 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.9744767 1.192469e-52
## kelompok:waktu 0.9744767 1.391737e-20
# Pendekatan multivariat
aov2$Anova
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99637 23850.6 1 87 < 2.2e-16 ***
## kelompok 2 0.16358 8.5 2 87 0.0004222 ***
## waktu 1 0.78126 101.2 3 85 < 2.2e-16 ***
## kelompok:waktu 2 0.53466 10.5 6 172 7.023e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix
aov2_rs <- anova_test(
data = dat_long,
dv = Hb,
wid = id,
between = kelompok,
within = waktu,
effect.size = "pes",
type = 3
)
get_anova_table(
aov2_rs,
correction = "GG"
)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 kelompok 2.00 87.00 8.507 4.22e-04 * 0.164
## 2 waktu 2.82 245.28 139.532 7.14e-51 * 0.616
## 3 kelompok:waktu 5.64 245.28 22.171 6.22e-20 * 0.338
# Ukuran efek
eta_squared(
aov2,
partial = TRUE
)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.16 | [0.05, 1.00]
## waktu | 0.62 | [0.56, 1.00]
## kelompok:waktu | 0.34 | [0.25, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
omega_squared(
aov2,
partial = TRUE
)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## ------------------------------------------------
## kelompok | 0.14 | [0.04, 1.00]
## waktu | 0.19 | [0.12, 1.00]
## kelompok:waktu | 0.07 | [0.01, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 4b.1 PLOT INTERAKSI
afex_plot(
aov2,
x = "waktu",
trace = "kelompok",
error = "within",
mapping = c("colour", "shape", "linetype")
) +
labs(
y = "Kadar Hb (g/dL)",
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 & POST HOC
# Estimated marginal means:
# waktu di dalam setiap kelompok
em2 <- emmeans(
aov2,
~ waktu | kelompok
)
em2
## kelompok = Kontrol:
## waktu emmean SE df lower.CL upper.CL
## M0 10.6 0.115 87 10.4 10.8
## M4 10.8 0.126 87 10.5 11.0
## M8 10.8 0.147 87 10.5 11.1
## M12 10.8 0.152 87 10.5 11.1
##
## kelompok = TTD:
## waktu emmean SE df lower.CL upper.CL
## M0 10.6 0.115 87 10.4 10.9
## M4 11.1 0.126 87 10.9 11.4
## M8 11.5 0.147 87 11.2 11.8
## M12 11.8 0.152 87 11.5 12.1
##
## kelompok = TTD+VitC:
## waktu emmean SE df lower.CL upper.CL
## M0 10.6 0.115 87 10.4 10.8
## M4 11.3 0.126 87 11.0 11.5
## M8 11.8 0.147 87 11.5 12.1
## M12 12.1 0.152 87 11.8 12.4
##
## Confidence level used: 0.95
# Efek WAKTU di dalam masing-masing kelompok
joint_tests(
aov2,
by = "kelompok"
)
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## kelompok = Kontrol:
## model term df1 df2 F.ratio p.value
## waktu 3 87 2.280 0.0850
##
## kelompok = TTD:
## model term df1 df2 F.ratio p.value
## waktu 3 87 48.324 <0.0001
##
## kelompok = TTD+VitC:
## model term df1 df2 F.ratio p.value
## waktu 3 87 86.213 <0.0001
# Efek KELOMPOK pada masing-masing waktu
joint_tests(
aov2,
by = "waktu"
)
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.044 0.9571
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 4.485 0.0140
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 12.746 <0.0001
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 18.635 <0.0001
# Setiap waktu dibandingkan dengan baseline M0
# di dalam masing-masing kelompok
contrast(
em2,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M4 - M0 0.170 0.0818 87 2.078 0.0813
## M8 - M0 0.177 0.0928 87 1.903 0.0813
## M12 - M0 0.243 0.0971 87 2.507 0.0421
##
## kelompok = TTD:
## contrast estimate SE df t.ratio p.value
## M4 - M0 0.500 0.0818 87 6.112 <0.0001
## M8 - M0 0.877 0.0928 87 9.445 <0.0001
## M12 - M0 1.123 0.0971 87 11.574 <0.0001
##
## kelompok = TTD+VitC:
## contrast estimate SE df t.ratio p.value
## M4 - M0 0.670 0.0818 87 8.190 <0.0001
## M8 - M0 1.173 0.0928 87 12.641 <0.0001
## M12 - M0 1.500 0.0971 87 15.455 <0.0001
##
## P value adjustment: holm method for 3 tests
# Perbandingan antarkelompok pada masing-masing waktu
em2b <- emmeans(
aov2,
~ kelompok | waktu
)
pairs(
em2b,
adjust = "tukey"
)
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Kontrol - TTD -0.0467 0.162 87 -0.287 0.9555
## Kontrol - (TTD+VitC) -0.0133 0.162 87 -0.082 0.9963
## TTD - (TTD+VitC) 0.0333 0.162 87 0.205 0.9770
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - TTD -0.3767 0.178 87 -2.122 0.0913
## Kontrol - (TTD+VitC) -0.5133 0.178 87 -2.892 0.0133
## TTD - (TTD+VitC) -0.1367 0.178 87 -0.770 0.7224
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - TTD -0.7467 0.208 87 -3.598 0.0015
## Kontrol - (TTD+VitC) -1.0100 0.208 87 -4.867 <0.0001
## TTD - (TTD+VitC) -0.2633 0.208 87 -1.269 0.4165
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - TTD -0.9267 0.215 87 -4.306 0.0001
## Kontrol - (TTD+VitC) -1.2700 0.215 87 -5.901 <0.0001
## TTD - (TTD+VitC) -0.3433 0.215 87 -1.595 0.2531
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4d. KONTRAS INTERAKSI
# Pertanyaan:
# Apakah perubahan Hb dari baseline (M0) ke minggu ke-12 (M12)
# berbeda antar kelompok?
#
# Kontras M12 - M0:
# (-1, 0, 0, 1)
em_full <- emmeans(
aov2,
~ waktu * kelompok
)
contrast(
em_full,
interaction = list(
waktu = list(
"M12-M0" = c(-1, 0, 0, 1)
),
kelompok = "pairwise"
),
adjust = "holm"
)
## waktu_custom kelompok_pairwise estimate SE df t.ratio p.value
## M12-M0 Kontrol - TTD -0.880 0.137 87 -6.411 <0.0001
## M12-M0 Kontrol - (TTD+VitC) -1.257 0.137 87 -9.155 <0.0001
## M12-M0 TTD - (TTD+VitC) -0.377 0.137 87 -2.744 0.0074
##
## P value adjustment: holm method for 3 tests
# 4d.1 TREN LINEAR PER KELOMPOK
contrast(
em2,
"poly"
)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## linear 0.73667 0.312 87 2.361 0.0205
## quadratic -0.10333 0.113 87 -0.914 0.3632
## cubic 0.22333 0.244 87 0.914 0.3633
##
## kelompok = TTD:
## contrast estimate SE df t.ratio p.value
## linear 3.74667 0.312 87 12.009 <0.0001
## quadratic -0.25333 0.113 87 -2.241 0.0276
## cubic -0.00667 0.244 87 -0.027 0.9783
##
## kelompok = TTD+VitC:
## contrast estimate SE df t.ratio p.value
## linear 5.00333 0.312 87 16.037 <0.0001
## quadratic -0.34333 0.113 87 -3.037 0.0032
## cubic -0.01000 0.244 87 -0.041 0.9674
tren_int <- summary(
contrast(
em_full,
interaction = c(
waktu = "poly",
kelompok = "pairwise"
),
adjust = "none"
)
)
# Ambil tren linear
tren_lin <- subset(
tren_int,
waktu_poly == "linear"
)
# Koreksi Holm untuk tiga perbandingan antarkelompok
tren_lin$p.holm <- p.adjust(
tren_lin$p.value,
"holm"
)
tren_lin
## waktu_poly kelompok_pairwise estimate SE df t.ratio p.value
## 1 linear Kontrol - TTD -3.010000 0.4412226 87 -6.821953 1.138307e-09
## 4 linear Kontrol - (TTD+VitC) -4.266667 0.4412226 87 -9.670100 1.907740e-15
## 7 linear TTD - (TTD+VitC) -1.256667 0.4412226 87 -2.848147 5.486643e-03
## p.holm
## 1 2.276615e-09
## 4 5.723219e-15
## 7 5.486643e-03
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
#
# LMM berguna apabila:
# - terdapat data Hb yang hilang;
# - jumlah pengukuran tidak lengkap;
# - ingin memodelkan korelasi intra-subjek;
# - ingin memasukkan variasi individual antar responden.
#
# Model 1: random intercept
# Model 2: random intercept + random slope waktu
lmm1 <- lmer(
Hb ~ kelompok * waktu + (1 | id),
data = dat_long,
REML = TRUE
)
lmm2 <- lmer(
Hb ~ kelompok * waktu + (1 + minggu | id),
data = dat_long,
REML = TRUE
)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0156912 (tol = 0.002, component 1)
## See ?lme4::convergence and ?lme4::troubleshooting.
# Perbandingan struktur random effect
anova(
lmm1,
lmm2,
refit = FALSE
)
## Data: dat_long
## Models:
## lmm1: Hb ~ kelompok * waktu + (1 | id)
## lmm2: Hb ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 553.33 607.74 -262.67 525.33
## lmm2 16 529.10 591.28 -248.55 497.10 28.233 2 7.401e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji efek tetap
anova(
lmm2,
ddf = "Kenward-Roger"
)
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 1.5730 0.7865 2 87.00 8.4915 0.0004278 ***
## waktu 29.7959 9.9320 3 185.37 106.7398 < 2.2e-16 ***
## kelompok:waktu 9.4242 1.5707 6 206.40 16.8618 8.17e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Intraclass correlation
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.800
## Unadjusted ICC: 0.545
# Ringkasan model
summary(lmm2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Hb ~ kelompok * waktu + (1 + minggu | id)
## Data: dat_long
##
## REML criterion at convergence: 497.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.54392 -0.57566 -0.04789 0.63841 2.14608
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.3050293 0.55229
## minggu 0.0006671 0.02583 0.69
## Residual 0.0926217 0.30434
## Number of obs: 360, groups: id, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.15278 0.07228 86.76358 154.294 < 2e-16 ***
## kelompok1 -0.40861 0.10222 86.76358 -3.997 0.000134 ***
## kelompok2 0.11556 0.10222 86.76358 1.130 0.261411
## waktu1 -0.53611 0.03223 161.70347 -16.635 < 2e-16 ***
## waktu2 -0.08944 0.02831 210.13257 -3.159 0.001814 **
## waktu3 0.20611 0.02831 210.13256 7.280 6.56e-12 ***
## kelompok1:waktu1 0.38861 0.04558 161.70347 8.526 1.02e-14 ***
## kelompok2:waktu1 -0.08889 0.04558 161.70347 -1.950 0.052873 .
## kelompok1:waktu2 0.11194 0.04004 210.13256 2.796 0.005654 **
## kelompok2:waktu2 -0.03556 0.04004 210.13256 -0.888 0.375524
## kelompok1:waktu3 -0.17694 0.04004 210.13256 -4.419 1.58e-05 ***
## kelompok2:waktu3 0.04556 0.04004 210.13256 1.138 0.256489
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) klmpk1 klmpk2 waktu1 waktu2 waktu3 klm1:1 klm2:1 klm1:2
## kelompok1 0.000
## kelompok2 0.000 -0.500
## waktu1 -0.396 0.000 0.000
## waktu2 -0.150 0.000 0.000 -0.185
## waktu3 0.150 0.000 0.000 -0.379 -0.358
## klmpk1:wkt1 0.000 -0.396 0.198 0.000 0.000 0.000
## klmpk2:wkt1 0.000 0.198 -0.396 0.000 0.000 0.000 -0.500
## klmpk1:wkt2 0.000 -0.150 0.075 0.000 0.000 0.000 -0.185 0.092
## klmpk2:wkt2 0.000 0.075 -0.150 0.000 0.000 0.000 0.092 -0.185 -0.500
## klmpk1:wkt3 0.000 0.150 -0.075 0.000 0.000 0.000 -0.379 0.190 -0.358
## klmpk2:wkt3 0.000 -0.075 0.150 0.000 0.000 0.000 0.190 -0.379 0.179
## klm2:2 klm1:3
## kelompok1
## kelompok2
## waktu1
## waktu2
## waktu3
## klmpk1:wkt1
## klmpk2:wkt1
## klmpk1:wkt2
## klmpk2:wkt2
## klmpk1:wkt3 0.179
## klmpk2:wkt3 -0.358 -0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0156912 (tol = 0.002, component 1)
## See ?lme4::convergence and ?lme4::troubleshooting.
# 5.1 DIAGNOSTIK RESIDUAL LMM
par(mfrow = c(1, 3))
qqnorm(
resid(lmm2),
main = "Q-Q Residual"
)
qqline(resid(lmm2))
qqnorm(
ranef(lmm2)$id[, 1],
main = "Q-Q 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))
# 5.2 CONTOH ANALISIS LMM DENGAN DATA HILANG
# Bagian ini opsional.
# Digunakan untuk menunjukkan bahwa LMM dapat mempertahankan responden
# yang memiliki sebagian pengukuran Hb yang hilang.
set.seed(1)
dat_miss <- dat_long
# Hanya membuat beberapa pengukuran setelah baseline menjadi NA
idx_miss <- sample(
which(
dat_miss$waktu != "M0" &
!is.na(dat_miss$Hb)
),
size = min(
30,
sum(
dat_miss$waktu != "M0" &
!is.na(dat_miss$Hb)
)
)
)
dat_miss$Hb[idx_miss] <- NA
# LMM dengan data hilang
lmm_miss <- lmer(
Hb ~ kelompok * waktu + (1 + minggu | id),
data = dat_miss,
REML = TRUE,
na.action = na.exclude,
control = lmerControl(
optimizer = "bobyqa"
)
)
anova(
lmm_miss,
ddf = "Kenward-Roger"
)
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 1.4926 0.7463 2 86.989 8.8102 0.000328 ***
## waktu 26.2863 8.7621 3 168.596 102.9432 < 2.2e-16 ***
## kelompok:waktu 9.5428 1.5905 6 186.327 18.6644 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Jumlah responden yang memiliki minimal satu Hb hilang
n_distinct(
dat_miss$id[is.na(dat_miss$Hb)]
)
## [1] 25
# 6. RINGKASAN & PENYIMPANAN HASIL
# 6.1 Simpan data panjang
write.csv(
dat_long,
"data_Hb_anemia_format_long.csv",
row.names = FALSE
)
# 6.2 Simpan data lebar
write.csv(
dat_wide,
"data_Hb_anemia_format_wide.csv",
row.names = FALSE
)
# 6.3 Simpan statistik deskriptif
write.csv(
desk_lengkap,
"deskriptif_Hb_menurut_kelompok_dan_waktu.csv",
row.names = FALSE
)
# 6.4 Simpan varians selisih antarwaktu
write.csv(
data.frame(
pasangan_waktu = names(var_selisih),
varians_selisih = as.numeric(var_selisih)
),
"varians_selisih_Hb.csv",
row.names = FALSE
)
# 6.5 Simpan grafik
ggsave(
"profile_plot_Hb.png",
p_profil,
width = 9,
height = 6,
dpi = 300
)
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
ggsave(
"spaghetti_plot_Hb.png",
p_spag,
width = 10,
height = 7,
dpi = 300
)
# 6.6 Tampilkan kembali hasil utama
cat("\n=============================================\n")
##
## =============================================
cat("HASIL UTAMA REPEATED MEASURE ANALYSIS Hb\n")
## HASIL UTAMA REPEATED MEASURE ANALYSIS Hb
cat("=============================================\n\n")
## =============================================
cat("Statistik deskriptif:\n")
## Statistik deskriptif:
print(desk_lengkap)
## # A tibble: 12 × 8
## kelompok waktu n mean sd median min max
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kontrol M0 30 10.6 0.625 10.7 9.2 12.3
## 2 Kontrol M4 30 10.8 0.772 10.8 9.3 12.8
## 3 Kontrol M8 30 10.8 0.879 10.5 9.2 13.4
## 4 Kontrol M12 30 10.8 0.916 10.9 9 13.5
## 5 TTD M0 30 10.6 0.625 10.8 9.5 12.1
## 6 TTD M4 30 11.1 0.637 11.3 9.8 12.5
## 7 TTD M8 30 11.5 0.826 11.6 9.9 13
## 8 TTD M12 30 11.8 0.830 11.9 10.2 12.9
## 9 TTD+VitC M0 30 10.6 0.636 10.8 9.4 11.6
## 10 TTD+VitC M4 30 11.3 0.646 11.4 9.5 12.4
## 11 TTD+VitC M8 30 11.8 0.694 11.9 10.1 12.8
## 12 TTD+VitC M12 30 12.1 0.746 12.1 10.4 13.4
cat("\n\nRepeated Measure ANOVA satu arah:\n")
##
##
## Repeated Measure ANOVA satu arah:
print(get_anova_table(aov1_rs, correction = "auto"))
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 87 143.204 1.52e-33 * 0.832
cat("\n\nMixed Design ANOVA:\n")
##
##
## Mixed Design ANOVA:
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 87.00 8.507 4.22e-04 * 0.164
## 2 waktu 2.82 245.28 139.532 7.14e-51 * 0.616
## 3 kelompok:waktu 5.64 245.28 22.171 6.22e-20 * 0.338
cat("\n\nLMM:\n")
##
##
## LMM:
print(anova(lmm2, ddf = "Kenward-Roger"))
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 1.5730 0.7865 2 87.00 8.4915 0.0004278 ***
## waktu 29.7959 9.9320 3 185.37 106.7398 < 2.2e-16 ***
## kelompok:waktu 9.4242 1.5707 6 206.40 16.8618 8.17e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n\nAnalisis selesai.\n")
##
##
## Analisis selesai.