# NAMA : SILVIA ANDRYANI
# NIM : 261108045
# TUGAS : BIOSTATIK INTERMEDIATE 5
# DOSEN : Dr. M. FATHURAHMAN, S.Si., M.Si.
# =============================================================================
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: eksperimen diabetes, gula darah puasa (GDP)
# Data menggunakan format long dari file CSV
#
# Isi:
# 0. Paket & pengaturan
# 1. Membaca data CSV 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, sfierisitas (Mauchly)
# 3b. ANOVA + koreksi Greenhouse-Geisser / Huynh-Feldt
# 3c. Pendekatan multivariat (MANOVA) sebagai pembanding
# 3d. Post hoc berpasangan & kontras polinomial (tren)
# 3e. Alternatif nonparametrik: uji Friedman
# 4. Mixed Design ANOVA (between: Kelompok x within: Waktu)
# 4a. Uji asumsi: outlier, normalitas, Levene, Box's M, Mauchly
# 4b. ANOVA campuran + ukuran efek
# 4c. Analisis efek sederhana (simple effects) & post hoc
# 4d. Kontras interaksi (perubahan 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("tidyverse", "afex", "emmeans", "rstatix", "car",
# "effectsize", "lme4", "lmerTest", "performance",
# "ggpubr", "Hmisc"))
suppressPackageStartupMessages({
library(dplyr) # manipulasi data
library(tidyr) # format panjang <-> lebar
library(ggplot2) # grafik
library(afex) # ANOVA within/mixed, koreksi GG/HF, MANOVA
library(emmeans) # rerata marginal, simple effects, post hoc, kontras
library(rstatix) # uji asumsi: outlier, Shapiro, Levene, Box's M
library(car) # leveneTest, Anova
library(effectsize) # eta kuadrat parsial, omega kuadrat
library(lme4) # linear mixed model
library(lmerTest) # uji F/t dengan Satterthwaite / Kenward-Roger
library(performance) # diagnostik model
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 CSV DAN MENYIAPKAN DATA
# Data sudah tersedia, sehingga bagian simulasi pada contoh diganti menjadi
# membaca file CSV. Struktur data:
# - id : identitas subjek
# - kelompok : kelompok perlakuan
# - usia : usia subjek
# - jk : jenis kelamin
# - waktu : M0, M4, M8, M12
# - gdp : gula darah puasa
# - minggu : 0, 4, 8, 12
file_csv <- "D:/Documents/S2 KESMAS UNMUL/MATERI S2 KESMAS/BIOSTATISTIK/BIOSTATIK 5/MATERI DAN TUGAS BIOSTATISTIK 5/DATA BUATAN/data_sintetis_diabetes_gdp_long_60_format_contoh.csv"
dat_long <- read.csv(file_csv, stringsAsFactors = FALSE)
dat_long <- dat_long |>
mutate(
id = factor(id),
kelompok = factor(kelompok,
levels = c("Standar",
"Standar_Diet",
"Standar_Diet_Aktivitas")),
usia = as.numeric(usia),
jk = factor(jk, levels = c("L", "P")),
waktu = factor(waktu, levels = c("M0", "M4", "M8", "M12")),
gdp = as.numeric(gdp),
minggu = as.numeric(minggu)
) |>
arrange(kelompok, id, minggu)
head(dat_long)
## id kelompok usia jk waktu gdp minggu
## 1 D001 Standar 48 L M0 168.9 0
## 2 D001 Standar 48 L M4 163.6 4
## 3 D001 Standar 48 L M8 159.8 8
## 4 D001 Standar 48 L M12 147.3 12
## 5 D002 Standar 43 P M0 154.2 0
## 6 D002 Standar 43 P M4 143.3 4
str(dat_long)
## 'data.frame': 240 obs. of 7 variables:
## $ id : Factor w/ 60 levels "D001","D002",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok: Factor w/ 3 levels "Standar","Standar_Diet",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num 48 48 48 48 43 43 43 43 55 55 ...
## $ jk : Factor w/ 2 levels "L","P": 1 1 1 1 2 2 2 2 1 1 ...
## $ waktu : Factor w/ 4 levels "M0","M4","M8",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ gdp : num 169 164 160 147 154 ...
## $ minggu : num 0 4 8 12 0 4 8 12 0 4 ...
# Cek kelengkapan data
nrow(dat_long)
## [1] 240
n_distinct(dat_long$id)
## [1] 60
table(dat_long$kelompok)
##
## Standar Standar_Diet Standar_Diet_Aktivitas
## 80 80 80
table(dat_long$waktu)
##
## M0 M4 M8 M12
## 60 60 60 60
table(dat_long$kelompok, dat_long$waktu)
##
## M0 M4 M8 M12
## Standar 20 20 20 20
## Standar_Diet 20 20 20 20
## Standar_Diet_Aktivitas 20 20 20 20
colSums(is.na(dat_long))
## id kelompok usia jk waktu gdp minggu
## 0 0 0 0 0 0 0
# Cek apakah setiap subjek punya 4 kali pengukuran
dat_long |>
count(id, name = "jumlah_pengukuran") |>
count(jumlah_pengukuran)
## jumlah_pengukuran n
## 1 4 60
# Format lebar (satu baris = satu subjek)
dat_wide <- dat_long |>
pivot_wider(
id_cols = c(id, kelompok, usia, jk),
names_from = waktu,
values_from = gdp,
names_prefix = "GDP_"
) |>
arrange(kelompok, id)
head(dat_wide)
## # A tibble: 6 × 8
## id kelompok usia jk GDP_M0 GDP_M4 GDP_M8 GDP_M12
## <fct> <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 D001 Standar 48 L 169. 164. 160. 147.
## 2 D002 Standar 43 P 154. 143. 128. 118.
## 3 D003 Standar 55 L 135. 136. 129. 133.
## 4 D004 Standar 47 L 203. 190. 197. 190.
## 5 D005 Standar 44 L 169. 165. 169. 164
## 6 D006 Standar 56 L 150. 138. 134. 122.
str(dat_wide)
## tibble [60 × 8] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 60 levels "D001","D002",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ kelompok: Factor w/ 3 levels "Standar","Standar_Diet",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : num [1:60] 48 43 55 47 44 56 62 38 41 70 ...
## $ jk : Factor w/ 2 levels "L","P": 1 2 1 1 1 1 1 2 2 2 ...
## $ GDP_M0 : num [1:60] 169 154 135 203 169 ...
## $ GDP_M4 : num [1:60] 164 143 136 190 165 ...
## $ GDP_M8 : num [1:60] 160 128 129 197 169 ...
## $ GDP_M12 : num [1:60] 147 118 133 190 164 ...
# 2. EKSPLORASI DATA
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(gdp, type = "mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Standar M0 gdp 20 177. 20.2
## 2 Standar M4 gdp 20 173. 21.1
## 3 Standar M8 gdp 20 168. 22.3
## 4 Standar M12 gdp 20 166. 25.6
## 5 Standar_Diet M0 gdp 20 177. 15.4
## 6 Standar_Diet M4 gdp 20 164. 15.9
## 7 Standar_Diet M8 gdp 20 156. 17.5
## 8 Standar_Diet M12 gdp 20 142. 18.8
## 9 Standar_Diet_Aktivitas M0 gdp 20 171. 23.1
## 10 Standar_Diet_Aktivitas M4 gdp 20 156. 22.9
## 11 Standar_Diet_Aktivitas M8 gdp 20 142. 25.3
## 12 Standar_Diet_Aktivitas M12 gdp 20 128. 25.7
# Matriks kovarians & korelasi antarwaktu
minggu <- c(0, 4, 8, 12)
kolom_gdp <- paste0("GDP_M", minggu)
S <- cov(dat_wide[, kolom_gdp])
R <- cor(dat_wide[, kolom_gdp])
round(S, 1)
## GDP_M0 GDP_M4 GDP_M8 GDP_M12
## GDP_M0 387.4 368.6 386.8 398.7
## GDP_M4 368.6 440.4 473.7 521.7
## GDP_M8 386.8 473.7 578.4 638.1
## GDP_M12 398.7 521.7 638.1 784.4
round(R, 2)
## GDP_M0 GDP_M4 GDP_M8 GDP_M12
## GDP_M0 1.00 0.89 0.82 0.72
## GDP_M4 0.89 1.00 0.94 0.89
## GDP_M8 0.82 0.94 1.00 0.95
## GDP_M12 0.72 0.89 0.95 1.00
# Varians selisih antarpasangan waktu (inti asumsi sfierisitas)
pasangan <- combn(kolom_gdp, 2)
var_selisih <- apply(pasangan, 2, function(p) {
var(dat_wide[[p[1]]] - dat_wide[[p[2]]])
})
names(var_selisih) <- apply(pasangan, 2, paste, collapse = " - ")
round(var_selisih, 1)
## GDP_M0 - GDP_M4 GDP_M0 - GDP_M8 GDP_M0 - GDP_M12 GDP_M4 - GDP_M8
## 90.7 192.2 374.4 71.4
## GDP_M4 - GDP_M12 GDP_M8 - GDP_M12
## 181.5 86.6
# Profile plot: rerata +/- 95% CI per kelompok
p_profil <- ggplot(dat_long, aes(minggu, gdp,
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 = "Gula darah puasa (GDP)",
colour = "Kelompok",
title = "Profil rerata GDP (+/- 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.
# Spaghetti plot: lintasan tiap subjek
p_spag <- ggplot(dat_long, aes(minggu, gdp, 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 = "GDP",
title = "Lintasan individu dan rerata kelompok")
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan: apakah GDP berubah selama 12 minggu pada kelompok
# Standar_Diet_Aktivitas?
d1 <- droplevels(filter(dat_long, kelompok == "Standar_Diet_Aktivitas"))
d1w <- filter(dat_wide, kelompok == "Standar_Diet_Aktivitas")
## 3a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per waktu
d1 |>
group_by(waktu) |>
identify_outliers(gdp)
## # A tibble: 3 × 9
## waktu id kelompok usia jk gdp minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 M0 D048 Standar_Diet_Aktiv… 64 P 103. 0 TRUE FALSE
## 2 M4 D048 Standar_Diet_Aktiv… 64 P 87.7 4 TRUE FALSE
## 3 M8 D048 Standar_Diet_Aktiv… 64 P 76.3 8 TRUE FALSE
# (ii) Normalitas per waktu (Shapiro-Wilk) dan Q-Q plot
d1 |>
group_by(waktu) |>
shapiro_test(gdp)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 gdp 0.895 0.0331
## 2 M4 gdp 0.893 0.0312
## 3 M8 gdp 0.954 0.433
## 4 M12 gdp 0.966 0.659
ggpubr::ggqqplot(d1, "gdp", facet.by = "waktu")
# (iii) Sfierisitas: Mauchly
aov1_rs <- anova_test(data = d1, dv = gdp, 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 57 98.59 1.54e-22 * 0.838
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.229 8.67e-05 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.522 1.57, 29.78 7.94e-13 * 0.561 1.68, 31.98 1.29e-13
## p[HF]<.05
## 1 *
get_anova_table(aov1_rs, correction = "auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 1.57 29.78 98.59 7.94e-13 * 0.838
## 3b. ANOVA dengan afex ------------------------------------------------------
aov1 <- aov_ez(id = "id", dv = "gdp", data = d1, within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov1
## Anova Table (Type 3 tests)
##
## Response: gdp
## Effect df MSE F ges pes p.value
## 1 waktu 1.57, 29.78 136.02 98.59 *** .319 .838 <.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) 1780792 1 40803 19 829.24 < 2.2e-16 ***
## waktu 21017 3 4050 57 98.59 < 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.22882 8.6723e-05
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.52242 7.942e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.5611078 1.289203e-13
# Ukuran efek tambahan
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.84 | [0.77, 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.31 | [0.13, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
## 3c. 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.97760 829.24 1 19 < 2.2e-16 ***
## waktu 1 0.87464 39.54 3 17 6.985e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 3d. Post hoc & kontras tren ------------------------------------------------
em1 <- emmeans(aov1, ~ waktu)
em1
## waktu emmean SE df lower.CL upper.CL
## M0 171 5.16 19 160 182
## M4 156 5.13 19 145 167
## M8 142 5.66 19 130 154
## M12 128 5.75 19 116 140
##
## Confidence level used: 0.95
pairs(em1, adjust = "bonferroni")
## contrast estimate SE df t.ratio p.value
## M0 - M4 15.0 2.11 19 7.104 <0.0001
## M0 - M8 29.4 3.00 19 9.780 <0.0001
## M0 - M12 43.5 3.83 19 11.370 <0.0001
## M4 - M8 14.4 2.03 19 7.082 <0.0001
## M4 - M12 28.6 2.81 19 10.174 <0.0001
## M8 - M12 14.2 1.59 19 8.913 <0.0001
##
## P value adjustment: bonferroni method for 6 tests
contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm")
## contrast estimate SE df t.ratio p.value
## M4 - M0 -15.0 2.11 19 -7.104 <0.0001
## M8 - M0 -29.4 3.00 19 -9.780 <0.0001
## M12 - M0 -43.5 3.83 19 -11.370 <0.0001
##
## P value adjustment: holm method for 3 tests
contrast(em1, "poly")
## contrast estimate SE df t.ratio p.value
## linear -144.96 12.90 19 -11.258 <0.0001
## quadratic 0.79 2.25 19 0.351 0.7294
## cubic -0.32 4.70 19 -0.068 0.9464
## 3e. Alternatif nonparametrik -----------------------------------------------
friedman_test(d1, gdp ~ waktu | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 gdp 20 55.3 3 5.87e-12 Friedman test
friedman_effsize(d1, gdp ~ waktu | id)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 gdp 20 0.922 Kendall W large
d1 |>
wilcox_test(gdp ~ 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 gdp M0 M4 20 20 209 0.00000381 0.0000229 ****
## 2 gdp M0 M8 20 20 210 0.00000191 0.0000114 ****
## 3 gdp M0 M12 20 20 210 0.00000191 0.0000114 ****
## 4 gdp M4 M8 20 20 206 0.0000134 0.0000801 ****
## 5 gdp M4 M12 20 20 210 0.00000191 0.0000114 ****
## 6 gdp M8 M12 20 20 209 0.00000381 0.0000229 ****
# (Opsional) ANOVA robust berbasis trimmed mean -- paket WRS2
if (requireNamespace("WRS2", quietly = TRUE)) {
print(WRS2::rmanova(d1$gdp, d1$waktu, d1$id, tr = 0.2))
}
## Call:
## WRS2::rmanova(y = d1$gdp, groups = d1$waktu, blocks = d1$id,
## tr = 0.2)
##
## Test statistic: F = 65.3718
## Degrees of freedom 1: 2.05
## Degrees of freedom 2: 22.56
## p-value: 0
# 4. MIXED DESIGN ANOVA (Kelompok [between] x Waktu [within])
# Pertanyaan: apakah pola perubahan GDP berbeda antarkelompok perlakuan?
## 4a. Uji asumsi -------------------------------------------------------------
# (i) Outlier per sel
dat_long |>
group_by(kelompok, waktu) |>
identify_outliers(gdp)
## # A tibble: 5 × 9
## kelompok waktu id usia jk gdp minggu is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 Standar_Diet M8 D030 51 P 120. 8 TRUE FALSE
## 2 Standar_Diet M8 D036 62 P 194. 8 TRUE FALSE
## 3 Standar_Diet_Aktiv… M0 D048 64 P 103. 0 TRUE FALSE
## 4 Standar_Diet_Aktiv… M4 D048 64 P 87.7 4 TRUE FALSE
## 5 Standar_Diet_Aktiv… M8 D048 64 P 76.3 8 TRUE FALSE
# (ii) Normalitas per sel (3 x 4 = 12 sel) dan residual model
dat_long |>
group_by(kelompok, waktu) |>
shapiro_test(gdp)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Standar M0 gdp 0.973 0.810
## 2 Standar M4 gdp 0.966 0.667
## 3 Standar M8 gdp 0.940 0.241
## 4 Standar M12 gdp 0.958 0.504
## 5 Standar_Diet M0 gdp 0.976 0.872
## 6 Standar_Diet M4 gdp 0.969 0.743
## 7 Standar_Diet M8 gdp 0.963 0.602
## 8 Standar_Diet M12 gdp 0.940 0.237
## 9 Standar_Diet_Aktivitas M0 gdp 0.895 0.0331
## 10 Standar_Diet_Aktivitas M4 gdp 0.893 0.0312
## 11 Standar_Diet_Aktivitas M8 gdp 0.954 0.433
## 12 Standar_Diet_Aktivitas M12 gdp 0.966 0.659
ggpubr::ggqqplot(dat_long, "gdp", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas varians antarkelompok pada tiap waktu
dat_long |>
group_by(waktu) |>
levene_test(gdp ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 57 0.635 0.533
## 2 M4 2 57 0.797 0.456
## 3 M8 2 57 0.877 0.422
## 4 M12 2 57 0.607 0.549
# (iv) Homogenitas matriks kovarians antarkelompok
box_m(dat_wide[, kolom_gdp], dat_wide$kelompok)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 19.8 0.471 20 Box's M-test for Homogeneity of Covariance Matric…
# (v) Sfierisitas: Mauchly (dari summary model di bawah)
## 4b. ANOVA campuran ---------------------------------------------------------
aov2 <- aov_ez(id = "id", dv = "gdp", data = dat_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: gdp
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 57 1685.46 5.48 ** .150 .161 .007
## 2 waktu 1.89, 107.71 80.55 194.90 *** .221 .774 <.001
## 3 kelompok:waktu 3.78, 107.71 80.55 19.80 *** .054 .410 <.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) 6135619 1 96071 57 3640.3291 < 2e-16 ***
## kelompok 18476 2 96071 57 5.4811 0.00665 **
## waktu 29667 3 8676 171 194.9002 < 2e-16 ***
## kelompok:waktu 6029 6 8676 171 19.8033 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.43072 5.9188e-09
## kelompok:waktu 0.43072 5.9188e-09
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.62989 < 2.2e-16 ***
## kelompok:waktu 0.62989 7.778e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.6508229 1.106445e-36
## kelompok:waktu 0.6508229 3.710862e-12
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.98458 3640.3 1 57 < 2.2e-16 ***
## kelompok 2 0.16130 5.5 2 57 0.00665 **
## waktu 1 0.83354 91.8 3 55 < 2.2e-16 ***
## kelompok:waktu 2 0.58289 7.7 6 112 6.388e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix (format ringkas untuk laporan)
aov2_rs <- anova_test(data = dat_long, dv = gdp, 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 57.00 5.481 7.00e-03 * 0.161
## 2 waktu 1.89 107.71 194.900 1.38e-35 * 0.774
## 3 kelompok:waktu 3.78 107.71 19.803 7.78e-12 * 0.410
# Ukuran efek
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.16 | [0.03, 1.00]
## waktu | 0.77 | [0.73, 1.00]
## kelompok:waktu | 0.41 | [0.31, 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.13 | [0.01, 1.00]
## waktu | 0.22 | [0.12, 1.00]
## kelompok:waktu | 0.05 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Plot interaksi dari model
afex_plot(aov2, x = "waktu", trace = "kelompok", error = "within",
mapping = c("colour", "shape", "linetype")) +
labs(y = "GDP", 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 ---------------------------------------------------------
em2 <- emmeans(aov2, ~ waktu | kelompok)
# Efek WAKTU di dalam tiap kelompok
joint_tests(aov2, by = "kelompok")
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## kelompok = Standar:
## model term df1 df2 F.ratio p.value
## waktu 3 57 5.202 0.0030
##
## kelompok = Standar_Diet:
## model term df1 df2 F.ratio p.value
## waktu 3 57 44.071 <0.0001
##
## kelompok = Standar_Diet_Aktivitas:
## model term df1 df2 F.ratio p.value
## waktu 3 57 67.410 <0.0001
# Efek KELOMPOK pada tiap waktu
joint_tests(aov2, by = "waktu")
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 57 0.619 0.5418
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 57 3.351 0.0421
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 57 6.899 0.0021
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 57 13.178 <0.0001
# Post hoc: tiap waktu vs baseline di dalam tiap kelompok
contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm")
## kelompok = Standar:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -4.20 1.87 57 -2.247 0.0285
## M8 - M0 -9.36 2.52 57 -3.706 0.0014
## M12 - M0 -11.37 3.10 57 -3.669 0.0014
##
## kelompok = Standar_Diet:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -13.57 1.87 57 -7.261 <0.0001
## M8 - M0 -21.78 2.52 57 -8.628 <0.0001
## M12 - M0 -35.24 3.10 57 -11.375 <0.0001
##
## kelompok = Standar_Diet_Aktivitas:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -14.96 1.87 57 -8.002 <0.0001
## M8 - M0 -29.36 2.52 57 -11.629 <0.0001
## M12 - M0 -43.52 3.10 57 -14.045 <0.0001
##
## P value adjustment: holm method for 3 tests
# Post hoc: perbandingan antarkelompok pada tiap waktu
em2b <- emmeans(aov2, ~ kelompok | waktu)
pairs(em2b, adjust = "tukey")
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Standar - Standar_Diet -0.515 6.26 57 -0.082 0.9963
## Standar - Standar_Diet_Aktivitas 5.765 6.26 57 0.920 0.6299
## Standar_Diet - Standar_Diet_Aktivitas 6.280 6.26 57 1.002 0.5784
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Standar - Standar_Diet 8.855 6.39 57 1.386 0.3547
## Standar - Standar_Diet_Aktivitas 16.520 6.39 57 2.587 0.0324
## Standar_Diet - Standar_Diet_Aktivitas 7.665 6.39 57 1.200 0.4580
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Standar - Standar_Diet 11.910 6.94 57 1.715 0.2083
## Standar - Standar_Diet_Aktivitas 25.765 6.94 57 3.711 0.0013
## Standar_Diet - Standar_Diet_Aktivitas 13.855 6.94 57 1.996 0.1225
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Standar - Standar_Diet 23.360 7.45 57 3.135 0.0075
## Standar - Standar_Diet_Aktivitas 37.915 7.45 57 5.088 <0.0001
## Standar_Diet - Standar_Diet_Aktivitas 14.555 7.45 57 1.953 0.1333
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## 4d. Kontras interaksi ------------------------------------------------------
# Apakah penurunan (M12 - M0) berbeda antarkelompok?
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
## M12-M0 Standar - Standar_Diet 23.88 4.38 57 5.448
## M12-M0 Standar - Standar_Diet_Aktivitas 32.15 4.38 57 7.337
## M12-M0 Standar_Diet - Standar_Diet_Aktivitas 8.28 4.38 57 1.888
## p.value
## <0.0001
## <0.0001
## 0.0641
##
## P value adjustment: holm method for 3 tests
# Tren linear per kelompok dan perbandingannya
contrast(em2, "poly")[c(1, 4, 7)]
## contrast kelompok estimate SE df t.ratio p.value
## linear Standar -39.3 10.2 57 -3.831 0.0003
## linear Standar_Diet -113.9 10.2 57 -11.117 <0.0001
## linear Standar_Diet_Aktivitas -145.0 10.2 57 -14.143 <0.0001
tren_int <- summary(
contrast(em_full,
interaction = c(waktu = "poly", kelompok = "pairwise"),
adjust = "none")
)
tren_lin <- subset(tren_int, waktu_poly == "linear")
tren_lin$p.holm <- p.adjust(tren_lin$p.value, "holm")
tren_lin
## waktu_poly kelompok_pairwise estimate SE df
## 1 linear Standar - Standar_Diet 74.680 14.49466 57
## 4 linear Standar - Standar_Diet_Aktivitas 105.695 14.49466 57
## 7 linear Standar_Diet - Standar_Diet_Aktivitas 31.015 14.49466 57
## t.ratio p.value p.holm
## 1 5.152241 3.342437e-06 6.684875e-06
## 4 7.291995 1.037498e-09 3.112495e-09
## 7 2.139753 3.666865e-02 3.666865e-02
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Tidak mensyaratkan sfierisitas, menampung data hilang (MAR) dan waktu
# pengukuran yang tidak seragam.
lmm1 <- lmer(gdp ~ kelompok * waktu + (1 | id),
data = dat_long, REML = TRUE)
lmm2 <- lmer(gdp ~ kelompok * waktu + (1 + minggu | id),
data = dat_long, REML = TRUE,
control = lmerControl(optimizer = "bobyqa"))
anova(lmm1, lmm2, refit = FALSE)
## Data: dat_long
## Models:
## lmm1: gdp ~ kelompok * waktu + (1 | id)
## lmm2: gdp ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 1823.0 1871.8 -897.53 1795.0
## lmm2 16 1775.4 1831.0 -871.68 1743.4 51.69 2 5.966e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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 258.5 129.26 2 57.00 5.4811 0.00665 **
## waktu 6684.5 2228.16 3 121.08 93.8206 < 2.2e-16 ***
## kelompok:waktu 1450.5 241.75 6 134.40 10.1625 2.961e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.890
## Unadjusted ICC: 0.596
# 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))
# Simulasi nilai hilang untuk menunjukkan keunggulan LMM
set.seed(1)
dat_miss <- dat_long
dat_miss$gdp[sample(which(dat_miss$waktu != "M0"), 20)] <- NA
lmm_miss <- lmer(gdp ~ kelompok * waktu + (1 + minggu | id),
data = dat_miss,
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 277.1 138.55 2 57.01 5.5473 0.006291 **
## waktu 7232.2 2410.74 3 109.16 95.8187 < 2.2e-16 ***
## kelompok:waktu 1546.3 257.72 6 120.17 10.2259 3.98e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# RM ANOVA akan membuang subjek yang punya minimal 1 nilai hilang
n_distinct(dat_miss$id[is.na(dat_miss$gdp)])
## [1] 19
# 6. MENYIMPAN DATA & RINGKASAN HASIL
dir.create("hasil_repeated_measures_diabetes", showWarnings = FALSE)
write.csv(dat_wide,
"hasil_repeated_measures_diabetes/data_diabetes_gdp_wide.csv",
row.names = FALSE)
write.csv(dat_long,
"hasil_repeated_measures_diabetes/data_diabetes_gdp_long.csv",
row.names = FALSE)
write.csv(desk,
"hasil_repeated_measures_diabetes/deskriptif_gdp.csv",
row.names = FALSE)
write.csv(as.data.frame(get_anova_table(aov1_rs, correction = "auto")),
"hasil_repeated_measures_diabetes/one_way_rm_anova.csv",
row.names = FALSE)
write.csv(as.data.frame(get_anova_table(aov2_rs, correction = "GG")),
"hasil_repeated_measures_diabetes/mixed_design_anova.csv",
row.names = FALSE)
ggsave("hasil_repeated_measures_diabetes/profile_plot_gdp.png",
p_profil, width = 8, height = 5, dpi = 300)
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
ggsave("hasil_repeated_measures_diabetes/spaghetti_plot_gdp.png",
p_spag, width = 8, height = 5, dpi = 300)
sink("hasil_repeated_measures_diabetes/session_info.txt")
sessionInfo()
sink()
list.files("hasil_repeated_measures_diabetes")
## [1] "data_diabetes_gdp_long.csv" "data_diabetes_gdp_wide.csv"
## [3] "deskriptif_gdp.csv" "mixed_design_anova.csv"
## [5] "one_way_rm_anova.csv" "profile_plot_gdp.png"
## [7] "session_info.txt" "spaghetti_plot_gdp.png"