# =============================================================================
# Nama : Syalmitha Auralia Nur
# NIM : 2611018004
# Sumber Rujukan : Rossen et al. (2021)
# Link DOI : https://doi.org/10.1186/s12966-021-01193-w
# REPEATED MEASURE ANALYSIS DENGAN R
# Intervensi Pemantauan Langkah terhadap Kadar HbA1c
#
# DATA RANCANGAN
# Kelompok : Kontrol, Langkah, Langkah+Konseling
# Responden : 90 orang (30 per kelompok)
# Pengukuran : M0, M6, M12, M18, M24
# Satuan : HbA1c (mmol/mol)
# =============================================================================
# 0. PAKET & PENGATURAN
# Instal sekali apabila ada paket yang belum tersedia:
# install.packages(c("dplyr", "tidyr", "ggplot2", "afex",
# "emmeans", "rstatix", "car", "effectsize",
# "lme4", "lmerTest", "performance", "ggpubr"))
suppressPackageStartupMessages({
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 grafik R Markdown untuk output Word
if (requireNamespace("knitr", quietly = TRUE)) {
knitr::opts_chunk$set(
dev = "jpeg",
fig.width = 8,
fig.height = 6,
dpi = 160,
fig.path = "figure-hba1c/",
echo = TRUE,
message = FALSE
)
}
# Lokasi data
setwd("C:/Users/HP/Downloads")
getwd()
## [1] "C:/Users/HP/Downloads"
list.files(pattern = "\\.csv$")
## [1] "data_balita 2.csv"
## [2] "data_balita.csv"
## [3] "data_hba1c_diabetes_long.csv"
## [4] "data_hba1c_diabetes_wide.csv"
## [5] "data_tds_hipertensi_long.csv"
## [6] "data_tds_hipertensi_wide.csv"
## [7] "mahasiswa-S2-Kesehatan-Masyarakat-angkatan-2026.csv"
# Folder hasil analisis
folder_hasil <- "C:/Users/HP/Downloads/hasil_analisis_HbA1c"
if (!dir.exists(folder_hasil)) {
dir.create(
folder_hasil,
recursive = TRUE,
showWarnings = FALSE
)
}
stopifnot(dir.exists(folder_hasil))
# 1. MEMANGGIL DATA (FORMAT WIDE DAN LONG)
n_per <- 30
kel_lab <- c(
"Kontrol",
"Langkah",
"Langkah+Konseling"
)
bulan <- c(0, 6, 12, 18, 24)
# Memanggil data wide
file_wide <- "C:/Users/HP/Downloads/data_hba1c_diabetes_wide.csv"
if (!file.exists(file_wide)) {
stop("File CSV tidak ditemukan. Periksa nama file di Downloads.")
}
# Menentukan pemisah CSV
baris_awal <- readLines(
file_wide,
n = 1,
warn = FALSE
)
pemisah <- if (
grepl(";", baris_awal, fixed = TRUE)
) ";" else ","
dat_wide <- read.table(
file_wide,
header = TRUE,
sep = pemisah,
dec = ".",
fileEncoding = "UTF-8-BOM",
stringsAsFactors = FALSE,
check.names = FALSE
)
# Pengaturan tipe data
dat_wide$id <- factor(dat_wide$id)
dat_wide$kelompok <- factor(
dat_wide$kelompok,
levels = kel_lab
)
dat_wide$jk <- factor(
dat_wide$jk,
levels = c("L", "P")
)
dat_wide$diagnosis <- factor(
dat_wide$diagnosis,
levels = c("Prediabetes", "DM tipe 2")
)
kolom_ukur <- paste0("HbA1c_M", bulan)
dat_wide[kolom_ukur] <- lapply(
dat_wide[kolom_ukur],
as.numeric
)
# Mengubah data wide menjadi long
dat_long <- dat_wide |>
pivot_longer(
all_of(kolom_ukur),
names_to = "waktu",
values_to = "hba1c"
) |>
mutate(
waktu = factor(
waktu,
levels = kolom_ukur,
labels = paste0("M", bulan)
),
bulan = as.numeric(
sub("M", "", as.character(waktu))
),
bulan6 = bulan / 6
)
# Pemeriksaan struktur data
stopifnot(
nrow(dat_wide) == 90,
nrow(dat_long) == 450,
n_distinct(dat_long$id) == 90,
all(table(dat_long$id) == 5),
!anyNA(dat_long$hba1c),
all(table(dat_wide$kelompok) == 30)
)
print(head(dat_wide))
## id kelompok usia jk diagnosis HbA1c_M0 HbA1c_M6 HbA1c_M12 HbA1c_M18
## 1 P001 Kontrol 60 P DM tipe 2 66.9 62.5 61.0 60.9
## 2 P002 Kontrol 76 L DM tipe 2 55.0 53.3 59.7 55.7
## 3 P003 Kontrol 68 P DM tipe 2 63.6 51.9 55.4 52.7
## 4 P004 Kontrol 64 L DM tipe 2 43.0 51.5 57.3 46.8
## 5 P005 Kontrol 72 P DM tipe 2 56.4 47.0 52.6 52.9
## 6 P006 Kontrol 54 P DM tipe 2 55.3 56.7 58.8 57.6
## HbA1c_M24
## 1 70.4
## 2 58.7
## 3 55.4
## 4 50.2
## 5 54.9
## 6 63.0
print(head(dat_long))
## # A tibble: 6 × 9
## id kelompok usia jk diagnosis waktu hba1c bulan bulan6
## <fct> <fct> <int> <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 P001 Kontrol 60 P DM tipe 2 M0 66.9 0 0
## 2 P001 Kontrol 60 P DM tipe 2 M6 62.5 6 1
## 3 P001 Kontrol 60 P DM tipe 2 M12 61 12 2
## 4 P001 Kontrol 60 P DM tipe 2 M18 60.9 18 3
## 5 P001 Kontrol 60 P DM tipe 2 M24 70.4 24 4
## 6 P002 Kontrol 76 L DM tipe 2 M0 55 0 0
str(dat_long)
## tibble [450 × 9] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 90 levels "P001","P002",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ kelompok : Factor w/ 3 levels "Kontrol","Langkah",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : int [1:450] 60 60 60 60 60 76 76 76 76 76 ...
## $ jk : Factor w/ 2 levels "L","P": 2 2 2 2 2 1 1 1 1 1 ...
## $ diagnosis: Factor w/ 2 levels "Prediabetes",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ waktu : Factor w/ 5 levels "M0","M6","M12",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ hba1c : num [1:450] 66.9 62.5 61 60.9 70.4 55 53.3 59.7 55.7 58.7 ...
## $ bulan : num [1:450] 0 6 12 18 24 0 6 12 18 24 ...
## $ bulan6 : num [1:450] 0 1 2 3 4 0 1 2 3 4 ...
print(table(dat_long$kelompok, dat_long$waktu))
##
## M0 M6 M12 M18 M24
## Kontrol 30 30 30 30 30
## Langkah 30 30 30 30 30
## Langkah+Konseling 30 30 30 30 30
# 2. EKSPLORASI DATA
# Statistik deskriptif
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(
hba1c,
type = "mean_sd"
)
print(desk)
## # A tibble: 15 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 hba1c 30 50.5 10.8
## 2 Kontrol M6 hba1c 30 50.4 9.41
## 3 Kontrol M12 hba1c 30 51.9 10.6
## 4 Kontrol M18 hba1c 30 51.6 10.0
## 5 Kontrol M24 hba1c 30 52.6 11.2
## 6 Langkah M0 hba1c 30 50.3 10.5
## 7 Langkah M6 hba1c 30 51.4 9.11
## 8 Langkah M12 hba1c 30 52.9 9.79
## 9 Langkah M18 hba1c 30 53.0 10.1
## 10 Langkah M24 hba1c 30 57.1 9.89
## 11 Langkah+Konseling M0 hba1c 30 50.0 9.27
## 12 Langkah+Konseling M6 hba1c 30 51.2 10.3
## 13 Langkah+Konseling M12 hba1c 30 52.6 9.74
## 14 Langkah+Konseling M18 hba1c 30 53.6 11.2
## 15 Langkah+Konseling M24 hba1c 30 53.1 10.8
# Matriks kovarians dan korelasi antarwaktu
S <- cov(dat_wide[, kolom_ukur])
R <- cor(dat_wide[, kolom_ukur])
print(round(S, 1))
## HbA1c_M0 HbA1c_M6 HbA1c_M12 HbA1c_M18 HbA1c_M24
## HbA1c_M0 102.1 85.3 87.8 90.6 86.8
## HbA1c_M6 85.3 90.7 85.9 89.6 88.7
## HbA1c_M12 87.8 85.9 98.7 95.2 95.8
## HbA1c_M18 90.6 89.6 95.2 107.8 98.7
## HbA1c_M24 86.8 88.7 95.8 98.7 114.7
print(round(R, 2))
## HbA1c_M0 HbA1c_M6 HbA1c_M12 HbA1c_M18 HbA1c_M24
## HbA1c_M0 1.00 0.89 0.87 0.86 0.80
## HbA1c_M6 0.89 1.00 0.91 0.91 0.87
## HbA1c_M12 0.87 0.91 1.00 0.92 0.90
## HbA1c_M18 0.86 0.91 0.92 1.00 0.89
## HbA1c_M24 0.80 0.87 0.90 0.89 1.00
# Varians selisih antarpasangan waktu
pasangan <- combn(kolom_ukur, 2)
var_selisih <- apply(
pasangan,
2,
function(p) {
var(dat_wide[[p[1]]] - dat_wide[[p[2]]])
}
)
names(var_selisih) <- apply(
pasangan,
2,
paste,
collapse = " - "
)
print(round(var_selisih, 1))
## HbA1c_M0 - HbA1c_M6 HbA1c_M0 - HbA1c_M12 HbA1c_M0 - HbA1c_M18
## 22.2 25.2 28.6
## HbA1c_M0 - HbA1c_M24 HbA1c_M6 - HbA1c_M12 HbA1c_M6 - HbA1c_M18
## 43.3 17.7 19.3
## HbA1c_M6 - HbA1c_M24 HbA1c_M12 - HbA1c_M18 HbA1c_M12 - HbA1c_M24
## 28.0 16.0 21.8
## HbA1c_M18 - HbA1c_M24
## 25.0
# Profile plot: rerata +/- 95% CI
rata_ci95 <- function(x) {
m <- mean(x, na.rm = TRUE)
n <- sum(!is.na(x))
se <- sd(x, na.rm = TRUE) / sqrt(n)
margin <- qt(0.975, df = n - 1) * se
data.frame(
y = m,
ymin = m - margin,
ymax = m + margin
)
}
p_profil <- ggplot(
dat_long,
aes(
bulan,
hba1c,
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 = rata_ci95,
geom = "errorbar",
width = 0.6
) +
scale_x_continuous(breaks = bulan) +
labs(
x = "Bulan ke-",
y = "HbA1c (mmol/mol)",
colour = "Kelompok",
title = "Profil Rerata HbA1c (95% CI)"
) +
theme(legend.position = "bottom")
print(p_profil)
# Spaghetti plot: lintasan tiap responden
p_spag <- ggplot(
dat_long,
aes(
bulan,
hba1c,
group = id
)
) +
geom_line(alpha = 0.3) +
stat_summary(
aes(group = kelompok),
fun = mean,
geom = "line",
colour = "firebrick",
linewidth = 1.2
) +
facet_wrap(~ kelompok) +
scale_x_continuous(breaks = bulan) +
labs(
x = "Bulan ke-",
y = "HbA1c (mmol/mol)",
title = "Lintasan Individu dan Rerata Kelompok"
)
print(p_spag)
# 3. REPEATED MEASURE ANOVA SATU ARAH
# WITHIN-SUBJECT: WAKTU
# Perubahan HbA1c selama 24 bulan pada kelompok Langkah+Konseling
d1 <- droplevels(
filter(
dat_long,
kelompok == "Langkah+Konseling"
)
)
d1w <- filter(
dat_wide,
kelompok == "Langkah+Konseling"
)
# 3a. UJI ASUMSI
# (i) Outlier per waktu
print(
d1 |>
group_by(waktu) |>
identify_outliers(hba1c)
)
## [1] waktu id kelompok usia jk diagnosis
## [7] hba1c bulan bulan6 is.outlier is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas Shapiro-Wilk
print(
d1 |>
group_by(waktu) |>
shapiro_test(hba1c)
)
## # A tibble: 5 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 hba1c 0.849 0.000577
## 2 M6 hba1c 0.922 0.0309
## 3 M12 hba1c 0.930 0.0485
## 4 M18 hba1c 0.923 0.0317
## 5 M24 hba1c 0.943 0.111
# Q-Q plot
p_qq1 <- ggpubr::ggqqplot(
d1,
"hba1c",
facet.by = "waktu"
)
print(p_qq1)
# (iii) Sferisitas Mauchly
aov1_rs <- anova_test(
data = d1,
dv = hba1c,
wid = id,
within = waktu,
effect.size = "pes"
)
print(aov1_rs)
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 4 116 5.293 0.000596 * 0.154
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.6 0.123
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05
## 1 waktu 0.775 3.1, 89.88 0.002 * 0.878 3.51, 101.86 0.001 *
print(
get_anova_table(
aov1_rs,
correction = "auto"
)
)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 4 116 5.293 0.000596 * 0.154
# 3b. REPEATED MEASURE ANOVA DENGAN AFEX
aov1 <- aov_ez(
id = "id",
dv = "hba1c",
data = d1,
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
print(aov1)
## Anova Table (Type 3 tests)
##
## Response: hba1c
## Effect df MSE F ges pes p.value
## 1 waktu 3.10, 89.88 15.71 5.29 ** .017 .154 .002
## ---
## 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) 407255 1 13931.0 29 847.7795 < 2.2e-16 ***
## waktu 258 4 1411.7 116 5.2929 0.0005958 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.60038 0.1234
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.77481 0.001871 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8781361 0.001104726
# Ukuran efek
eta_squared(
aov1,
partial = TRUE
)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.15 | [0.05, 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.01 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 3c. PENDEKATAN MULTIVARIAT (MANOVA)
print(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.96692 847.78 1 29 < 2e-16 ***
## waktu 1 0.33654 3.30 4 26 0.02594 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3d. POST HOC DAN KONTRAS TREN
em1 <- emmeans(aov1, ~ waktu)
print(em1)
## waktu emmean SE df lower.CL upper.CL
## M0 50.0 1.69 29 46.6 53.5
## M6 51.2 1.88 29 47.3 55.0
## M12 52.6 1.78 29 49.0 56.3
## M18 53.6 2.04 29 49.4 57.8
## M24 53.1 1.97 29 49.1 57.1
##
## Confidence level used: 0.95
# Semua pasangan waktu
pairs(
em1,
adjust = "bonferroni"
)
## contrast estimate SE df t.ratio p.value
## M0 - M6 -1.167 0.919 29 -1.269 1.0000
## M0 - M12 -2.600 0.992 29 -2.622 0.1379
## M0 - M18 -3.547 1.060 29 -3.352 0.0224
## M0 - M24 -3.083 1.200 29 -2.564 0.1579
## M6 - M12 -1.433 0.752 29 -1.905 0.6676
## M6 - M18 -2.380 0.767 29 -3.105 0.0423
## M6 - M24 -1.917 0.861 29 -2.225 0.3398
## M12 - M18 -0.947 0.704 29 -1.345 1.0000
## M12 - M24 -0.483 0.746 29 -0.648 1.0000
## M18 - M24 0.463 0.878 29 0.528 1.0000
##
## P value adjustment: bonferroni method for 10 tests
# Setiap waktu vs baseline
contrast(
em1,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## contrast estimate SE df t.ratio p.value
## M6 - M0 1.17 0.919 29 1.269 0.2145
## M12 - M0 2.60 0.992 29 2.622 0.0414
## M18 - M0 3.55 1.060 29 3.352 0.0090
## M24 - M0 3.08 1.200 29 2.564 0.0414
##
## P value adjustment: holm method for 4 tests
# Tren polinomial
contrast(
em1,
"poly"
)
## contrast estimate SE df t.ratio p.value
## linear 8.55 2.62 29 3.263 0.0028
## quadratic -3.75 2.15 29 -1.743 0.0919
## cubic -1.68 1.82 29 -0.922 0.3642
## quartic -0.17 4.39 29 -0.039 0.9694
# 3e. ALTERNATIF NONPARAMETRIK: FRIEDMAN
friedman_test(
d1,
hba1c ~ waktu | id
)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 hba1c 30 15.4 4 0.00398 Friedman test
friedman_effsize(
d1,
hba1c ~ waktu | id
)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 hba1c 30 0.128 Kendall W small
# Post hoc Wilcoxon berpasangan
d1 |>
wilcox_test(
hba1c ~ waktu,
paired = TRUE,
p.adjust.method = "bonferroni"
)
## # A tibble: 10 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 hba1c M0 M6 30 30 162 0.150 1 ns
## 2 hba1c M0 M12 30 30 118 0.0172 0.172 ns
## 3 hba1c M0 M18 30 30 82.5 0.00137 0.0137 *
## 4 hba1c M0 M24 30 30 120 0.0197 0.197 ns
## 5 hba1c M6 M12 30 30 136. 0.0478 0.478 ns
## 6 hba1c M6 M18 30 30 102. 0.00588 0.0588 ns
## 7 hba1c M6 M24 30 30 100 0.00538 0.0538 ns
## 8 hba1c M12 M18 30 30 171 0.213 1 ns
## 9 hba1c M12 M24 30 30 190. 0.396 1 ns
## 10 hba1c M18 M24 30 30 251 0.704 1 ns
# 4. MIXED DESIGN ANOVA
# BETWEEN: KELOMPOK x WITHIN: WAKTU
# 4a. UJI ASUMSI
# (i) Outlier per kelompok dan waktu
dat_long |>
group_by(kelompok, waktu) |>
identify_outliers(hba1c)
## # A tibble: 1 × 11
## kelompok waktu id usia jk diagnosis hba1c bulan bulan6 is.outlier
## <fct> <fct> <fct> <int> <fct> <fct> <dbl> <dbl> <dbl> <lgl>
## 1 Langkah M0 P045 77 L DM tipe 2 75.5 0 0 TRUE
## # ℹ 1 more variable: is.extreme <lgl>
# (ii) Normalitas per kelompok dan waktu
dat_long |>
group_by(kelompok, waktu) |>
shapiro_test(hba1c)
## # A tibble: 15 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 hba1c 0.918 0.0245
## 2 Kontrol M6 hba1c 0.946 0.136
## 3 Kontrol M12 hba1c 0.933 0.0605
## 4 Kontrol M18 hba1c 0.946 0.128
## 5 Kontrol M24 hba1c 0.943 0.108
## 6 Langkah M0 hba1c 0.909 0.0143
## 7 Langkah M6 hba1c 0.925 0.0363
## 8 Langkah M12 hba1c 0.950 0.174
## 9 Langkah M18 hba1c 0.954 0.212
## 10 Langkah M24 hba1c 0.953 0.203
## 11 Langkah+Konseling M0 hba1c 0.849 0.000577
## 12 Langkah+Konseling M6 hba1c 0.922 0.0309
## 13 Langkah+Konseling M12 hba1c 0.930 0.0485
## 14 Langkah+Konseling M18 hba1c 0.923 0.0317
## 15 Langkah+Konseling M24 hba1c 0.943 0.111
# Q-Q Plot
p_qq2 <- ggpubr::ggqqplot(
dat_long,
"hba1c",
ggtheme = theme_bw()
) +
facet_grid(waktu ~ kelompok)
print(p_qq2)
# (iii) Homogenitas varians Levene
dat_long |>
group_by(waktu) |>
levene_test(hba1c ~ kelompok)
## # A tibble: 5 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 0.416 0.661
## 2 M6 2 87 0.137 0.872
## 3 M12 2 87 0.401 0.671
## 4 M18 2 87 0.363 0.697
## 5 M24 2 87 0.761 0.470
# (iv) Homogenitas matriks kovarians Box's M
box_m(
dat_wide[, kolom_ukur],
dat_wide$kelompok
)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 38.8 0.129 30 Box's M-test for Homogeneity of Covariance Matric…
# 4b. MIXED DESIGN ANOVA
aov2 <- aov_ez(
id = "id",
dv = "hba1c",
data = dat_long,
between = "kelompok",
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
print(aov2)
## Anova Table (Type 3 tests)
##
## Response: hba1c
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 473.29 0.19 .004 .004 .830
## 2 waktu 3.33, 289.68 14.23 18.39 *** .019 .175 <.001
## 3 kelompok:waktu 6.66, 289.68 14.23 2.90 ** .006 .063 .007
## ---
## 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) 1223611 1 41176 87 2585.3210 < 2.2e-16 ***
## kelompok 177 2 41176 87 0.1865 0.83020
## waktu 872 4 4122 348 18.3941 1.006e-13 ***
## kelompok:waktu 275 8 4122 348 2.9002 0.00384 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.71585 0.00077217
## kelompok:waktu 0.71585 0.00077217
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.83242 8.378e-12 ***
## kelompok:waktu 0.83242 0.006973 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8695202 3.144467e-12
## kelompok:waktu 0.8695202 6.105600e-03
print(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.96744 2585.32 1 87 < 2.2e-16 ***
## kelompok 2 0.00427 0.19 2 87 0.83020
## waktu 1 0.34140 10.89 4 84 3.698e-07 ***
## kelompok:waktu 2 0.20159 2.38 8 170 0.01858 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA versi rstatix (Tipe III)
aov2_rs <- anova_test(
data = dat_long,
dv = hba1c,
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 0.186 8.30e-01 0.004
## 2 waktu 3.33 289.68 18.394 8.38e-12 * 0.175
## 3 kelompok:waktu 6.66 289.68 2.900 7.00e-03 * 0.063
# Ukuran efek
eta_squared(
aov2,
partial = TRUE
)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 4.27e-03 | [0.00, 1.00]
## waktu | 0.17 | [0.11, 1.00]
## kelompok:waktu | 0.06 | [0.01, 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.00 | [0.00, 1.00]
## waktu | 0.02 | [0.00, 1.00]
## kelompok:waktu | 3.91e-03 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Plot interaksi model
p_afex <- afex_plot(
aov2,
x = "waktu",
trace = "kelompok",
error = "within",
mapping = c("colour", "shape", "linetype")
) +
labs(
y = "HbA1c (mmol/mol)",
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"
print(p_afex)
# 4b.1 LINE PLOT DENGAN ERROR BARS (MEAN +/- SE)
# Menghitung summary statistics untuk plot
df_summary <- dat_long %>%
group_by(waktu, kelompok) %>%
summarise(
mean_hba1c = mean(hba1c, na.rm = TRUE),
sd_hba1c = sd(hba1c, na.rm = TRUE),
n = sum(!is.na(hba1c)),
se_hba1c = sd_hba1c / sqrt(n),
.groups = "drop"
)
print(df_summary)
## # A tibble: 15 × 6
## waktu kelompok mean_hba1c sd_hba1c n se_hba1c
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 M0 Kontrol 50.5 10.8 30 1.97
## 2 M0 Langkah 50.3 10.5 30 1.92
## 3 M0 Langkah+Konseling 50.0 9.27 30 1.69
## 4 M6 Kontrol 50.4 9.41 30 1.72
## 5 M6 Langkah 51.4 9.11 30 1.66
## 6 M6 Langkah+Konseling 51.2 10.3 30 1.88
## 7 M12 Kontrol 51.9 10.6 30 1.93
## 8 M12 Langkah 52.9 9.79 30 1.79
## 9 M12 Langkah+Konseling 52.6 9.74 30 1.78
## 10 M18 Kontrol 51.6 10.0 30 1.83
## 11 M18 Langkah 53.0 10.1 30 1.85
## 12 M18 Langkah+Konseling 53.6 11.2 30 2.04
## 13 M24 Kontrol 52.6 11.2 30 2.04
## 14 M24 Langkah 57.1 9.89 30 1.81
## 15 M24 Langkah+Konseling 53.1 10.8 30 1.97
# Membuat Line Plot dengan Error Bars
p_interaksi <- ggplot(
df_summary,
aes(
x = waktu,
y = mean_hba1c,
group = kelompok,
color = kelompok
)
) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
geom_errorbar(
aes(
ymin = mean_hba1c - se_hba1c,
ymax = mean_hba1c + se_hba1c
),
width = 0.1
) +
scale_color_manual(
values = c(
"Kontrol" = "#E46726",
"Langkah" = "#2C82C9",
"Langkah+Konseling" = "#379467"
)
) +
labs(
title = "Perubahan Kadar HbA1c Seiring Waktu",
subtitle = "Interaksi antara Kelompok dan Waktu (Mixed ANOVA)",
x = "Waktu Pengukuran",
y = "Rata-rata HbA1c (mmol/mol)",
color = "Kelompok"
) +
theme_minimal() +
theme(
plot.title = element_text(
face = "bold",
size = 14
),
legend.position = "bottom"
)
print(p_interaksi)
# 4c. ANALISIS EFEK SEDERHANA DAN POST HOC
em2 <- emmeans(
aov2,
~ waktu | kelompok
)
# Efek waktu di dalam tiap kelompok
joint_tests(
aov2,
by = "kelompok"
)
## kelompok = Kontrol:
## model term df1 df2 F.ratio p.value
## waktu 4 87 1.592 0.1837
##
## kelompok = Langkah:
## model term df1 df2 F.ratio p.value
## waktu 4 87 10.961 <0.0001
##
## kelompok = Langkah+Konseling:
## model term df1 df2 F.ratio p.value
## waktu 4 87 3.860 0.0062
# Efek kelompok pada tiap waktu
joint_tests(
aov2,
by = "waktu"
)
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.015 0.9852
##
## waktu = M6:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.089 0.9147
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.088 0.9157
##
## waktu = M18:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.267 0.7663
##
## waktu = M24:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 1.569 0.2141
# Post hoc tiap waktu vs baseline
contrast(
em2,
"trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M6 - M0 -0.0567 0.863 87 -0.066 0.9478
## M12 - M0 1.3800 0.921 87 1.498 0.4130
## M18 - M0 1.1667 0.971 87 1.202 0.4653
## M24 - M0 2.1367 1.160 87 1.845 0.2738
##
## kelompok = Langkah:
## contrast estimate SE df t.ratio p.value
## M6 - M0 1.0900 0.863 87 1.263 0.2100
## M12 - M0 2.5767 0.921 87 2.798 0.0190
## M18 - M0 2.6300 0.971 87 2.709 0.0190
## M24 - M0 6.7267 1.160 87 5.808 <0.0001
##
## kelompok = Langkah+Konseling:
## contrast estimate SE df t.ratio p.value
## M6 - M0 1.1667 0.863 87 1.352 0.1800
## M12 - M0 2.6000 0.921 87 2.823 0.0177
## M18 - M0 3.5467 0.971 87 3.654 0.0018
## M24 - M0 3.0833 1.160 87 2.662 0.0185
##
## P value adjustment: holm method for 4 tests
# Perbandingan antarkelompok pada setiap waktu
em2b <- emmeans(
aov2,
~ kelompok | waktu
)
pairs(
em2b,
adjust = "tukey"
)
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Kontrol - Langkah 0.147 2.64 87 0.056 0.9983
## Kontrol - (Langkah+Konseling) 0.447 2.64 87 0.169 0.9843
## Langkah - (Langkah+Konseling) 0.300 2.64 87 0.114 0.9929
##
## waktu = M6:
## contrast estimate SE df t.ratio p.value
## Kontrol - Langkah -1.000 2.49 87 -0.402 0.9147
## Kontrol - (Langkah+Konseling) -0.777 2.49 87 -0.313 0.9476
## Langkah - (Langkah+Konseling) 0.223 2.49 87 0.090 0.9956
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - Langkah -1.050 2.59 87 -0.405 0.9136
## Kontrol - (Langkah+Konseling) -0.773 2.59 87 -0.298 0.9521
## Langkah - (Langkah+Konseling) 0.277 2.59 87 0.107 0.9937
##
## waktu = M18:
## contrast estimate SE df t.ratio p.value
## Kontrol - Langkah -1.317 2.70 87 -0.487 0.8776
## Kontrol - (Langkah+Konseling) -1.933 2.70 87 -0.715 0.7551
## Langkah - (Langkah+Konseling) -0.617 2.70 87 -0.228 0.9717
##
## waktu = M24:
## contrast estimate SE df t.ratio p.value
## Kontrol - Langkah -4.443 2.75 87 -1.617 0.2440
## Kontrol - (Langkah+Konseling) -0.500 2.75 87 -0.182 0.9819
## Langkah - (Langkah+Konseling) 3.943 2.75 87 1.435 0.3276
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4d. KONTRAS INTERAKSI
# Perbedaan perubahan M24 - M0 antarkelompok
em_full <- emmeans(
aov2,
~ waktu * kelompok
)
kontras_akhir <- contrast(
em_full,
interaction = list(
waktu = list(
"M24-M0" = c(-1, 0, 0, 0, 1)
),
kelompok = "pairwise"
),
adjust = "holm"
)
print(kontras_akhir)
## waktu_custom kelompok_pairwise estimate SE df t.ratio p.value
## M24-M0 Kontrol - Langkah -4.590 1.64 87 -2.802 0.0188
## M24-M0 Kontrol - (Langkah+Konseling) -0.947 1.64 87 -0.578 0.5648
## M24-M0 Langkah - (Langkah+Konseling) 3.643 1.64 87 2.224 0.0574
##
## P value adjustment: holm method for 3 tests
# Tren polinomial per kelompok
tren_kelompok <- as.data.frame(
contrast(em2, "poly")
)
subset(
tren_kelompok,
contrast == "linear"
)
## contrast kelompok estimate SE df t.ratio p.value
## 1 linear Kontrol 5.496667 2.582513 87 2.128418 3.612998e-02
## 5 linear Langkah 14.993333 2.582513 87 5.805715 1.023210e-07
## 9 linear Langkah+Konseling 8.546667 2.582513 87 3.309438 1.361411e-03
# Perbandingan tren linear antarkelompok
tren_int <- as.data.frame(
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,
method = "holm"
)
print(tren_lin)
## waktu_poly kelompok_pairwise estimate SE df t.ratio
## 1 linear Kontrol - Langkah -9.496667 3.652225 87 -2.6002416
## 5 linear Kontrol - (Langkah+Konseling) -3.050000 3.652225 87 -0.8351074
## 9 linear Langkah - (Langkah+Konseling) 6.446667 3.652225 87 1.7651342
## p.value p.holm
## 1 0.01094446 0.03283337
## 5 0.40594470 0.40594470
## 9 0.08105026 0.16210053
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Pengaturan optimizer
kontrol_lmm <- lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 2e5)
)
# Model 1: Random intercept
lmm1 <- lmerTest::lmer(
hba1c ~ kelompok * waktu + (1 | id),
data = dat_long,
REML = TRUE,
control = kontrol_lmm
)
# Model 2: Random intercept dan random slope
lmm2 <- lmerTest::lmer(
hba1c ~ kelompok * waktu + (1 + bulan6 | id),
data = dat_long,
REML = TRUE,
control = kontrol_lmm
)
# Perbandingan struktur efek acak
anova(
lmm1,
lmm2,
refit = FALSE
)
## Data: dat_long
## Models:
## lmm1: hba1c ~ kelompok * waktu + (1 | id)
## lmm2: hba1c ~ kelompok * waktu + (1 + bulan6 | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 17 2736.3 2806.1 -1351.1 2702.3
## lmm2 19 2716.1 2794.1 -1339.0 2678.1 24.207 2 5.54e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji F Tipe III efek tetap
anova(
lmm2,
type = 3,
ddf = "Satterthwaite"
)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 3.40 1.702 2 87.00 0.1865 0.830198
## waktu 412.91 103.228 4 229.98 11.3133 2.131e-08 ***
## kelompok:waktu 198.11 24.763 8 229.98 2.7139 0.007109 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Korelasi intrakelas
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.886
## Unadjusted ICC: 0.862
# Pemeriksaan model singular
isSingular(lmm2)
## [1] FALSE
# 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 30 nilai hilang pada pengukuran pasca-baseline
set.seed(1)
dat_miss <- dat_long
idx_miss <- sample(
which(dat_miss$waktu != "M0"),
30
)
dat_miss$hba1c[idx_miss] <- NA
# LMM dengan data hilang
lmm_miss <- lmerTest::lmer(
hba1c ~ kelompok * waktu + (1 | id),
data = dat_miss,
REML = TRUE,
na.action = na.omit,
control = kontrol_lmm
)
# Uji F Tipe III pada data hilang
anova(
lmm_miss,
type = 3,
ddf = "Satterthwaite"
)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 5.27 2.635 2 87.06 0.2192 0.803617
## waktu 866.54 216.634 4 318.45 18.0169 2.374e-13 ***
## kelompok:waktu 256.94 32.117 8 318.45 2.6711 0.007485 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Menghitung jumlah subjek dengan data hilang
n_distinct(
dat_miss$id[is.na(dat_miss$hba1c)]
)
## [1] 25
# 6. MENYIMPAN DATA, GRAFIK, DAN RINGKASAN HASIL
folder_hasil <- "C:/Users/HP/Downloads/hasil_analisis_HbA1c"
if (!dir.exists(folder_hasil)) {
dir.create(
folder_hasil,
recursive = TRUE,
showWarnings = FALSE
)
}
stopifnot(dir.exists(folder_hasil))
# Menyimpan data wide
write.table(
dat_wide,
file.path(folder_hasil, "hasil_data_wide.csv"),
sep = ";",
dec = ".",
row.names = FALSE,
quote = FALSE
)
# Menyimpan data long
write.table(
dat_long |> select(-bulan6),
file.path(folder_hasil, "hasil_data_long.csv"),
sep = ";",
dec = ".",
row.names = FALSE,
quote = FALSE
)
# Menyimpan statistik deskriptif
write.table(
desk,
file.path(folder_hasil, "statistik_deskriptif.csv"),
sep = ";",
dec = ".",
row.names = FALSE,
quote = FALSE
)
# Menyimpan summary statistics Line Plot
write.table(
df_summary,
file.path(folder_hasil, "summary_line_plot_HbA1c.csv"),
sep = ";",
dec = ".",
row.names = FALSE,
quote = FALSE
)
# Menyimpan grafik dalam format PDF
simpan_pdf <- function(plot, nama_file, lebar, tinggi) {
lokasi <- file.path(folder_hasil, nama_file)
grDevices::pdf(
file = lokasi,
width = lebar,
height = tinggi,
onefile = TRUE,
useDingbats = FALSE
)
tryCatch(
print(plot),
finally = grDevices::dev.off()
)
message("Grafik tersimpan: ", lokasi)
invisible(lokasi)
}
# Menyimpan profile plot
simpan_pdf(
p_profil,
"profile_plot_HbA1c.pdf",
lebar = 9,
tinggi = 5.6
)
# Menyimpan spaghetti plot
simpan_pdf(
p_spag,
"spaghetti_plot_HbA1c.pdf",
lebar = 10,
tinggi = 5.3
)
# Menyimpan Line Plot dengan Error Bars
simpan_pdf(
p_interaksi,
"line_plot_interaksi_HbA1c.pdf",
lebar = 9,
tinggi = 5.6
)
# Menyimpan plot interaksi Mixed ANOVA
simpan_pdf(
p_afex,
"plot_mixed_anova_HbA1c.pdf",
lebar = 9,
tinggi = 5.6
)
# Menyimpan Q-Q plot
simpan_pdf(
p_qq2,
"qqplot_normalitas_HbA1c.pdf",
lebar = 11,
tinggi = 8
)
# Menyimpan ringkasan hasil analisis
capture.output(
{
cat("HASIL ANALISIS DATA SIMULASI HbA1c\n\n")
cat("Repeated Measure ANOVA - Langkah+Konseling\n")
print(aov1)
cat("\nMixed Design ANOVA - Semua Kelompok\n")
print(aov2)
cat("\nRingkasan Rerata dan SD\n")
print(desk)
cat("\nSummary Statistics Line Plot (Mean +/- SE)\n")
print(df_summary)
cat("\nKontras Perubahan Baseline-Akhir\n")
print(kontras_akhir)
cat("\nLinear Mixed Model\n")
print(
anova(
lmm2,
type = 3,
ddf = "Satterthwaite"
)
)
},
file = file.path(
folder_hasil,
"ringkasan_hasil_analisis.txt"
)
)
# Memeriksa seluruh file hasil
print(list.files(folder_hasil))
## [1] "fig_knit_unnamed-chunk-10-1.pdf" "fig_knit_unnamed-chunk-17-1.pdf"
## [3] "fig_knit_unnamed-chunk-18-1.pdf" "fig_knit_unnamed-chunk-20-1.pdf"
## [5] "fig_knit_unnamed-chunk-24-1.pdf" "fig_knit_unnamed-chunk-7-1.pdf"
## [7] "fig_knit_unnamed-chunk-7-2.pdf" "hasil_data_long.csv"
## [9] "hasil_data_wide.csv" "line_plot_interaksi_HbA1c.pdf"
## [11] "plot_mixed_anova_HbA1c.pdf" "profile_plot_HbA1c.pdf"
## [13] "profile_plot_HbA1c.png" "qqplot_normalitas_HbA1c.pdf"
## [15] "ringkasan_hasil_analisis.txt" "spaghetti_plot_HbA1c.pdf"
## [17] "spaghetti_plot_HbA1c.png" "statistik_deskriptif.csv"
## [19] "summary_line_plot_HbA1c.csv"
message(
"Selesai. Hasil tersimpan di: ",
normalizePath(folder_hasil)
)
# =============================================================================