##..............................................................................
#NAMA : WULANDARI
#NIM : 2611018022
##..............................................................................
#PERUBAHAN HbA1c (%) PADA PASIEN DIABETES TIPE 2
#DENGAN TIGA KELOMPOK INTERVENSI SELAMA 3 BULAN
##..............................................................................
suppressPackageStartupMessages({
library(readxl)
library(dplyr) # manipulasi data
library(tidyr) # format panjang <-> lebar
library(ggplot2) # grafik
library(afex) # ANOVA within/mixed (tipe III, koreksi GG/HF, MANOVA)
library(emmeans) # rerata marginal, efek sederhana, post hoc, kontras
library(rstatix) # uji asumsi yang ramah pipe (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 derajat bebas Satterthwaite / Kenward-Roger
library(performance)
library(ggpubr)
})
options(contrasts = c("contr.sum", "contr.poly")) # kontras jumlah-nol untuk SS tipe III
afex_options(emmeans_model = "multivariate") # post hoc memakai model multivariat (tahan terhadap non-sfierisitas)
theme_set(theme_bw(base_size = 12))
data_wide <- read_excel("D:/BIOSTATISTIK/Data_Eksperimen_DM_Wide.xlsx", sheet = "Dta DM_Wide")
View(data_wide)
str(data_wide)
## tibble [90 × 8] (S3: tbl_df/tbl/data.frame)
## $ ID_Responden : chr [1:90] "P001" "P002" "P003" "P004" ...
## $ Umur_Tahun : num [1:90] 46 59 54 50 47 60 46 65 58 62 ...
## $ JK : chr [1:90] "P" "P" "P" "L" ...
## $ Kelompok_Perlakuan: chr [1:90] "Kelompok A (Diet & Olahraga)" "Kelompok A (Diet & Olahraga)" "Kelompok A (Diet & Olahraga)" "Kelompok A (Diet & Olahraga)" ...
## $ Pemeriksaan_bln0 : num [1:90] 8.3 9.4 9 8.2 8.6 8.5 9.4 9.2 9 8.6 ...
## $ Pemeriksaan_bln1 : num [1:90] 8 9.2 8.8 7.8 8.4 8.3 9.2 8.8 8.7 8.4 ...
## $ Pemeriksaan_bln2 : num [1:90] 7.8 8.9 8.7 7.6 8.2 8.2 8.9 8.6 8.5 8.2 ...
## $ Pemeriksaan_bln3 : num [1:90] 7.5 8.7 8.4 7.5 8 7.9 8.7 8.5 8.2 8 ...
data_wide <- data_wide %>%
mutate(
id = factor(ID_Responden),
kelompok = factor(Kelompok_Perlakuan,
levels = c("Kelompok A (Diet & Olahraga)",
"Kelompok B (Metformin Standar)",
"Kelompok C (Metformin + Herbal)")),
umur = as.numeric(Umur_Tahun),
jk = factor(JK)
)
data_wide %>% count(kelompok)
## # A tibble: 3 × 2
## kelompok n
## <fct> <int>
## 1 Kelompok A (Diet & Olahraga) 30
## 2 Kelompok B (Metformin Standar) 30
## 3 Kelompok C (Metformin + Herbal) 30
data_long <- data_wide %>%
pivot_longer(
cols = starts_with("Pemeriksaan_bln"),
names_to = "waktu",
values_to = "pemeriksaan"
) %>%
mutate(
waktu = factor(waktu,
levels = c("Pemeriksaan_bln0", "Pemeriksaan_bln1",
"Pemeriksaan_bln2", "Pemeriksaan_bln3"),
labels = c("Bln-0", "Bln-1", "Bln-2", "Bln-3")),
bulan = as.numeric(sub("Bln", "", waktu))
)
head(data_long)
## # A tibble: 6 × 11
## ID_Responden Umur_Tahun JK Kelompok_Perlakuan id kelompok umur jk
## <chr> <dbl> <chr> <chr> <fct> <fct> <dbl> <fct>
## 1 P001 46 P Kelompok A (Diet & O… P001 Kelompo… 46 P
## 2 P001 46 P Kelompok A (Diet & O… P001 Kelompo… 46 P
## 3 P001 46 P Kelompok A (Diet & O… P001 Kelompo… 46 P
## 4 P001 46 P Kelompok A (Diet & O… P001 Kelompo… 46 P
## 5 P002 59 P Kelompok A (Diet & O… P002 Kelompo… 59 P
## 6 P002 59 P Kelompok A (Diet & O… P002 Kelompo… 59 P
## # ℹ 3 more variables: waktu <fct>, pemeriksaan <dbl>, bulan <dbl>
str(data_long)
## tibble [360 × 11] (S3: tbl_df/tbl/data.frame)
## $ ID_Responden : chr [1:360] "P001" "P001" "P001" "P001" ...
## $ Umur_Tahun : num [1:360] 46 46 46 46 59 59 59 59 54 54 ...
## $ JK : chr [1:360] "P" "P" "P" "P" ...
## $ Kelompok_Perlakuan: chr [1:360] "Kelompok A (Diet & Olahraga)" "Kelompok A (Diet & Olahraga)" "Kelompok A (Diet & Olahraga)" "Kelompok A (Diet & Olahraga)" ...
## $ id : Factor w/ 90 levels "P001","P002",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok : Factor w/ 3 levels "Kelompok A (Diet & Olahraga)",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ umur : num [1:360] 46 46 46 46 59 59 59 59 54 54 ...
## $ jk : Factor w/ 2 levels "L","P": 2 2 2 2 2 2 2 2 2 2 ...
## $ waktu : Factor w/ 4 levels "Bln-0","Bln-1",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ pemeriksaan : num [1:360] 8.3 8 7.8 7.5 9.4 9.2 8.9 8.7 9 8.8 ...
## $ bulan : num [1:360] 0 -1 -2 -3 0 -1 -2 -3 0 -1 ...
colSums(is.na(data_wide))
## ID_Responden Umur_Tahun JK Kelompok_Perlakuan
## 0 0 0 0
## Pemeriksaan_bln0 Pemeriksaan_bln1 Pemeriksaan_bln2 Pemeriksaan_bln3
## 0 0 0 0
## id kelompok umur jk
## 0 0 0 0
colSums(is.na(data_long))
## ID_Responden Umur_Tahun JK Kelompok_Perlakuan
## 0 0 0 0
## id kelompok umur jk
## 0 0 0 0
## waktu pemeriksaan bulan
## 0 0 0
data_long %>% count(id, waktu) %>% filter(n != 1)
## # A tibble: 0 × 3
## # ℹ 3 variables: id <fct>, waktu <fct>, n <int>
dim(data_wide)
## [1] 90 12
dim(data_long)
## [1] 360 11
table(data_long$kelompok, data_long$waktu)
##
## Bln-0 Bln-1 Bln-2 Bln-3
## Kelompok A (Diet & Olahraga) 30 30 30 30
## Kelompok B (Metformin Standar) 30 30 30 30
## Kelompok C (Metformin + Herbal) 30 30 30 30
# 1. EKSPLORASI DATA
# 1a. Statistik deskriptif per kelompok dan waktu
deskriptif <- data_long %>%
group_by(kelompok, waktu) %>%
get_summary_stats(pemeriksaan, type = "mean_sd")
deskriptif
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kelompok A (Diet & Olahraga) Bln-0 pemeriksaan 30 8.61 0.566
## 2 Kelompok A (Diet & Olahraga) Bln-1 pemeriksaan 30 8.34 0.57
## 3 Kelompok A (Diet & Olahraga) Bln-2 pemeriksaan 30 8.15 0.573
## 4 Kelompok A (Diet & Olahraga) Bln-3 pemeriksaan 30 7.92 0.594
## 5 Kelompok B (Metformin Standar) Bln-0 pemeriksaan 30 8.52 0.542
## 6 Kelompok B (Metformin Standar) Bln-1 pemeriksaan 30 7.95 0.567
## 7 Kelompok B (Metformin Standar) Bln-2 pemeriksaan 30 7.49 0.554
## 8 Kelompok B (Metformin Standar) Bln-3 pemeriksaan 30 7.18 0.561
## 9 Kelompok C (Metformin + Herbal) Bln-0 pemeriksaan 30 8.24 0.474
## 10 Kelompok C (Metformin + Herbal) Bln-1 pemeriksaan 30 7.45 0.439
## 11 Kelompok C (Metformin + Herbal) Bln-2 pemeriksaan 30 6.85 0.434
## 12 Kelompok C (Metformin + Herbal) Bln-3 pemeriksaan 30 6.40 0.443
# 1b.Rerata perubahan dari baseline
perubahan <- data_wide %>%
mutate(
perubahan_Bln3_Bln0 = Pemeriksaan_bln3 - Pemeriksaan_bln0,
penurunan_Bln0_Bln3 = Pemeriksaan_bln0 - Pemeriksaan_bln3
) %>%
group_by(kelompok) %>%
summarise(
n = n(),
mean_perubahan = mean(perubahan_Bln3_Bln0, na.rm = TRUE),
mean_penurunan = mean(penurunan_Bln0_Bln3, na.rm = TRUE),
sd_penurunan = sd(penurunan_Bln0_Bln3, na.rm = TRUE)
)
perubahan
## # A tibble: 3 × 5
## kelompok n mean_perubahan mean_penurunan sd_penurunan
## <fct> <int> <dbl> <dbl> <dbl>
## 1 Kelompok A (Diet & Olahraga) 30 -0.69 0.69 0.118
## 2 Kelompok B (Metformin Standa… 30 -1.35 1.35 0.161
## 3 Kelompok C (Metformin + Herb… 30 -1.84 1.84 0.192
# 1c. Matriks Kovarians dan korelasi
S <- cov(data_wide[, c("Pemeriksaan_bln0","Pemeriksaan_bln1",
"Pemeriksaan_bln2","Pemeriksaan_bln3")],
use = "complete.obs")
R <- cor(data_wide[, c("Pemeriksaan_bln0","Pemeriksaan_bln1",
"Pemeriksaan_bln2","Pemeriksaan_bln3")],
use = "complete.obs")
round(S, 3)
## Pemeriksaan_bln0 Pemeriksaan_bln1 Pemeriksaan_bln2
## Pemeriksaan_bln0 0.299 0.324 0.342
## Pemeriksaan_bln1 0.324 0.410 0.463
## Pemeriksaan_bln2 0.342 0.463 0.555
## Pemeriksaan_bln3 0.361 0.501 0.606
## Pemeriksaan_bln3
## Pemeriksaan_bln0 0.361
## Pemeriksaan_bln1 0.501
## Pemeriksaan_bln2 0.606
## Pemeriksaan_bln3 0.674
round(R, 3)
## Pemeriksaan_bln0 Pemeriksaan_bln1 Pemeriksaan_bln2
## Pemeriksaan_bln0 1.000 0.925 0.841
## Pemeriksaan_bln1 0.925 1.000 0.972
## Pemeriksaan_bln2 0.841 0.972 1.000
## Pemeriksaan_bln3 0.805 0.953 0.991
## Pemeriksaan_bln3
## Pemeriksaan_bln0 0.805
## Pemeriksaan_bln1 0.953
## Pemeriksaan_bln2 0.991
## Pemeriksaan_bln3 1.000
# 1d. Profile plot rerata ± 95% CI
p_profil <- ggplot(data_long, aes(bulan, pemeriksaan,
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 = .15) +
scale_x_continuous(breaks = 0:3) +
labs(x = "Bulan", y = "Nilai pemeriksaan HbA1c",
colour = "Kelompok", title = "Profil rerata nilai pemeriksaan") +
theme_bw()
p_profil
# 1e. Spaghetti plot
data_long$bulan <- as.numeric(gsub("[^0-9]", "", data_long$bulan))
p_spag <- ggplot(data_long, aes(bulan, pemeriksaan, group = id)) +
geom_line(alpha = .25) +
stat_summary(aes(group = kelompok), fun = mean, geom = "line", linewidth = 1.2) +
facet_wrap(~ kelompok) +
scale_x_continuous(breaks = 0:3, labels = c("0", "1", "2", "3")) +
labs(x = "Bulan ke-", y = "Nilai pemeriksaan HbA1c", title = "Lintasan individu dan rerata kelompok") +
theme_bw()
p_spag
# 2. REPEATED MEASURES ANOVA SATU ARAH
# 2a. Satu kelompok
d1 <- data_long %>%
filter(kelompok == "Kelompok A (Diet & Olahraga)") %>%
droplevels()
d1w <- data_wide %>%
filter(kelompok == "Kelompok A (Diet & Olahraga)") %>%
droplevels()
# 2b. Uji Asumsi ...............................................................
# (i) Uji Outlier Per waktu
d1 %>% group_by(waktu) %>% identify_outliers(pemeriksaan)
## [1] waktu ID_Responden Umur_Tahun JK
## [5] Kelompok_Perlakuan id kelompok umur
## [9] jk pemeriksaan bulan is.outlier
## [13] is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Uji Normalitas per waktu
d1 %>% group_by(waktu) %>% shapiro_test(pemeriksaan)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Bln-0 pemeriksaan 0.952 0.192
## 2 Bln-1 pemeriksaan 0.961 0.323
## 3 Bln-2 pemeriksaan 0.945 0.123
## 4 Bln-3 pemeriksaan 0.939 0.0847
ggpubr::ggqqplot(d1, "pemeriksaan", facet.by = "waktu")
# (iii) Uji sphericity + RM ANOVA
aov1_rs <- anova_test(
data = d1, dv = pemeriksaan, wid = id, within = waktu,
effect.size = "pes")
aov1_rs # berisi : ANOVA, Mauchly's Test, Koreksi GG & HF
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 87 602.751 4.53e-58 * 0.954
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.539 0.004 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.716 2.15, 62.29 2.83e-42 * 0.775 2.33, 67.43 1.46e-45
## p[HF]<.05
## 1 *
get_anova_table(aov1_rs, correction = "auto") # auto : GG dipakai jika Mauchly p<0.5
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2.15 62.29 602.751 2.83e-42 * 0.954
# 2b. Uji ANOVA dengan Afex
aov1 <- aov_ez(
id = "id", dv = "pemeriksaan", data = d1, within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG")
)
aov1
## Anova Table (Type 3 tests)
##
## Response: pemeriksaan
## Effect df MSE F ges pes p.value
## 1 waktu 2.15, 62.29 0.01 602.75 *** .167 .954 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov1)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 8182.4 1 38.086 29 6230.37 < 2.2e-16 ***
## waktu 7.7 3 0.371 87 602.75 < 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.53944 0.0043201
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.71595 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.7750481 1.456296e-45
# (i) Ukuran Efek
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.95 | [0.94, 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.16 | [0.04, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 2c. Pendekatan Multivariat
aov1$Anova # Pillai, Wilks, Hotelling-Lawley,Roy
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99537 6230.4 1 29 < 2.2e-16 ***
## waktu 1 0.97466 346.2 3 27 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 2d. Post Hoct
em1 <- emmeans(aov1, ~ waktu)
em1
## waktu emmean SE df lower.CL upper.CL
## Bln.0 8.61 0.103 29 8.40 8.82
## Bln.1 8.34 0.104 29 8.13 8.56
## Bln.2 8.15 0.105 29 7.94 8.36
## Bln.3 7.92 0.108 29 7.70 8.15
##
## Confidence level used: 0.95
pairs(em1, adjust = "bonferroni") # semua pasangan waktu (6 perbandingan)
## contrast estimate SE df t.ratio p.value
## Bln.0 - Bln.1 0.270 0.0174 29 15.529 <0.0001
## Bln.0 - Bln.2 0.463 0.0200 29 23.111 <0.0001
## Bln.0 - Bln.3 0.690 0.0216 29 31.902 <0.0001
## Bln.1 - Bln.2 0.193 0.0117 29 16.554 <0.0001
## Bln.1 - Bln.3 0.420 0.0155 29 27.163 <0.0001
## Bln.2 - Bln.3 0.227 0.0126 29 17.954 <0.0001
##
## P value adjustment: bonferroni method for 6 tests
contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm") # tiap waktu vs baseline
## contrast estimate SE df t.ratio p.value
## Bln.1 - Bln.0 -0.270 0.0174 29 -15.529 <0.0001
## Bln.2 - Bln.0 -0.463 0.0200 29 -23.111 <0.0001
## Bln.3 - Bln.0 -0.690 0.0216 29 -31.902 <0.0001
##
## P value adjustment: holm method for 3 tests
contrast(em1, "poly") # Tren linear, kuadratik, kubik
## contrast estimate SE df t.ratio p.value
## linear -2.2633 0.0699 29 -32.384 <0.0001
## quadratic 0.0433 0.0223 29 1.941 0.0620
## cubic -0.1100 0.0340 29 -3.233 0.0030
# 2e. Alternatif nonparametrik..................................................
friedman_test(d1, pemeriksaan ~ waktu | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 pemeriksaan 30 90 3 2.19e-19 Friedman test
friedman_effsize(d1, pemeriksaan ~ waktu | id) # Kendalll's W
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 pemeriksaan 30 1 Kendall W large
d1 %>% wilcox_test(pemeriksaan ~ 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 pemeriksaan Bln-0 Bln-1 30 30 465 1.86e-9 1.12e-8 ****
## 2 pemeriksaan Bln-0 Bln-2 30 30 465 1.86e-9 1.12e-8 ****
## 3 pemeriksaan Bln-0 Bln-3 30 30 465 1.86e-9 1.12e-8 ****
## 4 pemeriksaan Bln-1 Bln-2 30 30 465 1.86e-9 1.12e-8 ****
## 5 pemeriksaan Bln-1 Bln-3 30 30 465 1.86e-9 1.12e-8 ****
## 6 pemeriksaan Bln-2 Bln-3 30 30 465 1.86e-9 1.12e-8 ****
# 3. MIXED DESIGN ANOVA
# 3a. uji Asumsi................................................................
# (i) Uji Outlier per sel
data_long |> group_by(kelompok, waktu) |> identify_outliers(pemeriksaan)
## [1] kelompok waktu ID_Responden Umur_Tahun
## [5] JK Kelompok_Perlakuan id umur
## [9] jk pemeriksaan bulan is.outlier
## [13] is.extreme
## <0 rows> (or 0-length row.names)
# (ii) Normalitas per sel
data_long %>%
group_by(kelompok, waktu) %>%
shapiro_test(pemeriksaan)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kelompok A (Diet & Olahraga) Bln-0 pemeriksaan 0.952 0.192
## 2 Kelompok A (Diet & Olahraga) Bln-1 pemeriksaan 0.961 0.323
## 3 Kelompok A (Diet & Olahraga) Bln-2 pemeriksaan 0.945 0.123
## 4 Kelompok A (Diet & Olahraga) Bln-3 pemeriksaan 0.939 0.0847
## 5 Kelompok B (Metformin Standar) Bln-0 pemeriksaan 0.947 0.139
## 6 Kelompok B (Metformin Standar) Bln-1 pemeriksaan 0.937 0.0739
## 7 Kelompok B (Metformin Standar) Bln-2 pemeriksaan 0.958 0.276
## 8 Kelompok B (Metformin Standar) Bln-3 pemeriksaan 0.961 0.324
## 9 Kelompok C (Metformin + Herbal) Bln-0 pemeriksaan 0.943 0.109
## 10 Kelompok C (Metformin + Herbal) Bln-1 pemeriksaan 0.944 0.119
## 11 Kelompok C (Metformin + Herbal) Bln-2 pemeriksaan 0.960 0.318
## 12 Kelompok C (Metformin + Herbal) Bln-3 pemeriksaan 0.968 0.482
ggpubr::ggqqplot(data_long, "pemeriksaan", ggtheme = theme_bw()) +
facet_grid(waktu ~ kelompok)
# (iii) Homogenitas Varians
data_long %>% group_by(waktu) %>% levene_test(pemeriksaan ~ kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Bln-0 2 87 0.825 0.442
## 2 Bln-1 2 87 1.79 0.173
## 3 Bln-2 2 87 1.56 0.216
## 4 Bln-3 2 87 1.53 0.222
# (iv) Homogrnitas matriks kovarians antarkelompok
box_m(
data_wide[, c("Pemeriksaan_bln0","Pemeriksaan_bln1",
"Pemeriksaan_bln2","Pemeriksaan_bln3")],
data_wide$kelompok
)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 31.1 0.0537 20 Box's M-test for Homogeneity of Covariance Matric…
# 3b. Mixed Anova...............................................................
aov2 <- aov_ez(
id = "id", dv = "pemeriksaan", data = data_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges", "pes"), correction = "GG")
)
aov2
## Anova Table (Type 3 tests)
##
## Response: pemeriksaan
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 1.10 28.76 *** .393 .398 <.001
## 2 waktu 2.20, 191.49 0.01 3540.89 *** .463 .976 <.001
## 3 kelompok:waktu 4.40, 191.49 0.01 245.63 *** .107 .850 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2) # Mouchly, epsilon GG/HF, p terkoreksi
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 21669.0 1 95.554 87 19729.109 < 2.2e-16 ***
## kelompok 63.2 2 95.554 87 28.755 2.591e-10 ***
## waktu 84.3 3 2.070 261 3540.886 < 2.2e-16 ***
## kelompok:waktu 11.7 6 2.070 261 245.627 < 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.57292 4.0174e-09
## kelompok:waktu 0.57292 4.0174e-09
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.73369 < 2.2e-16 ***
## kelompok:waktu 0.73369 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.7535227 9.334283e-160
## kelompok:waktu 0.7535227 3.040800e-79
aov2$Anova # uji multivariat untuk efek within & interaksi
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99561 19729.1 1 87 < 2.2e-16 ***
## kelompok 2 0.39797 28.8 2 87 2.591e-10 ***
## waktu 1 0.98563 1943.8 3 85 < 2.2e-16 ***
## kelompok:waktu 2 0.97626 27.3 6 172 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Versi rstatix (hasil identik; format ringkas untuk laporan)
aov2_rs <- anova_test(data = data_long, dv = pemeriksaan, 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.0 87.00 28.755 2.59e-10 * 0.398
## 2 waktu 2.2 191.49 3540.886 1.25e-155 * 0.976
## 3 kelompok:waktu 4.4 191.49 245.627 3.15e-77 * 0.850
# Ukuran efek
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.40 | [0.26, 1.00]
## waktu | 0.98 | [0.97, 1.00]
## kelompok:waktu | 0.85 | [0.83, 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.38 | [0.25, 1.00]
## waktu | 0.46 | [0.39, 1.00]
## kelompok:waktu | 0.11 | [0.04, 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 = "pemeriksaan HbA1c (%)", 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"
# 3c. Efek sederhana (karena interaksi signifikan)..............................
em2 <- emmeans(aov2, ~ waktu | kelompok)
# Efek Waktu di dalam tiap kelompok (uji F gabungan per kelompok)
joint_tests(aov2, by = "kelompok")
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## kelompok = Kelompok A (Diet & Olahraga):
## model term df1 df2 F.ratio p.value
## waktu 3 87 191.223 <0.0001
##
## kelompok = Kelompok B (Metformin Standar):
## model term df1 df2 F.ratio p.value
## waktu 3 87 725.406 <0.0001
##
## kelompok = Kelompok C (Metformin + Herbal):
## model term df1 df2 F.ratio p.value
## waktu 3 87 1345.121 <0.0001
# Efek KELOMPOK pada tiap waktu
joint_tests(aov2, by = "waktu")
## waktu = Bln.0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 4.074 0.0204
##
## waktu = Bln.1:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 21.690 <0.0001
##
## waktu = Bln.2:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 46.334 <0.0001
##
## waktu = Bln.3:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 60.703 <0.0001
# Post hoc: tiap waktu vs baseline di dalam tiap kelompok
contrast(em2, "trt.vs.ctrl", ref = 1, adjust = "holm")
## kelompok = Kelompok A (Diet & Olahraga):
## contrast estimate SE df t.ratio p.value
## Bln.1 - Bln.0 -0.270 0.0222 87 -12.170 <0.0001
## Bln.2 - Bln.0 -0.463 0.0267 87 -17.366 <0.0001
## Bln.3 - Bln.0 -0.690 0.0293 87 -23.578 <0.0001
##
## kelompok = Kelompok B (Metformin Standar):
## contrast estimate SE df t.ratio p.value
## Bln.1 - Bln.0 -0.570 0.0222 87 -25.692 <0.0001
## Bln.2 - Bln.0 -1.037 0.0267 87 -38.854 <0.0001
## Bln.3 - Bln.0 -1.347 0.0293 87 -46.016 <0.0001
##
## kelompok = Kelompok C (Metformin + Herbal):
## contrast estimate SE df t.ratio p.value
## Bln.1 - Bln.0 -0.793 0.0222 87 -35.759 <0.0001
## Bln.2 - Bln.0 -1.393 0.0267 87 -52.222 <0.0001
## Bln.3 - Bln.0 -1.843 0.0293 87 -62.987 <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") # Tukey per waktu
## waktu = Bln.0:
## contrast estimate
## Kelompok A (Diet & Olahraga) - Kelompok B (Metformin Standar) 0.090
## Kelompok A (Diet & Olahraga) - (Kelompok C (Metformin + Herbal)) 0.373
## Kelompok B (Metformin Standar) - (Kelompok C (Metformin + Herbal)) 0.283
## SE df t.ratio p.value
## 0.137 87 0.659 0.7876
## 0.137 87 2.735 0.0205
## 0.137 87 2.075 0.1009
##
## waktu = Bln.1:
## contrast estimate
## Kelompok A (Diet & Olahraga) - Kelompok B (Metformin Standar) 0.390
## Kelompok A (Diet & Olahraga) - (Kelompok C (Metformin + Herbal)) 0.897
## Kelompok B (Metformin Standar) - (Kelompok C (Metformin + Herbal)) 0.507
## SE df t.ratio p.value
## 0.137 87 2.857 0.0147
## 0.137 87 6.568 <0.0001
## 0.137 87 3.711 0.0010
##
## waktu = Bln.2:
## contrast estimate
## Kelompok A (Diet & Olahraga) - Kelompok B (Metformin Standar) 0.663
## Kelompok A (Diet & Olahraga) - (Kelompok C (Metformin + Herbal)) 1.303
## Kelompok B (Metformin Standar) - (Kelompok C (Metformin + Herbal)) 0.640
## SE df t.ratio p.value
## 0.135 87 4.899 <0.0001
## 0.135 87 9.626 <0.0001
## 0.135 87 4.727 <0.0001
##
## waktu = Bln.3:
## contrast estimate
## Kelompok A (Diet & Olahraga) - Kelompok B (Metformin Standar) 0.747
## Kelompok A (Diet & Olahraga) - (Kelompok C (Metformin + Herbal)) 1.527
## Kelompok B (Metformin Standar) - (Kelompok C (Metformin + Herbal)) 0.780
## SE df t.ratio p.value
## 0.139 87 5.388 <0.0001
## 0.139 87 11.018 <0.0001
## 0.139 87 5.629 <0.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4d. Kontras interaksi -------------------------------------------------------
# Apakah penurunan (Bln3 - Bln0) berbeda antarkelompok? -- inti pertanyaan uji klinis
em_full <- emmeans(aov2, ~ waktu * kelompok)
contrast(em_full, interaction = list(waktu = list("Bln3-Bln0" = c(-1, 0, 0, 1)),
kelompok = "pairwise"),adjust = "holm")
## waktu_custom
## Bln3-Bln0
## Bln3-Bln0
## Bln3-Bln0
## kelompok_pairwise estimate
## Kelompok A (Diet & Olahraga) - Kelompok B (Metformin Standar) 0.657
## Kelompok A (Diet & Olahraga) - (Kelompok C (Metformin + Herbal)) 1.153
## Kelompok B (Metformin Standar) - (Kelompok C (Metformin + Herbal)) 0.497
## SE df t.ratio p.value
## 0.0414 87 15.866 <0.0001
## 0.0414 87 27.867 <0.0001
## 0.0414 87 12.000 <0.0001
##
## 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 Kelompok A (Diet & Olahraga) -2.26 0.0966 87 -23.422 <0.0001
## linear Kelompok B (Metformin Standar) -4.51 0.0966 87 -46.638 <0.0001
## linear Kelompok C (Metformin + Herbal) -6.13 0.0966 87 -63.437 <0.0001
tren_int <- summary(contrast(em_full, interaction = c(waktu = "poly", kelompok = "pairwise"),
adjust = "none"))
tren_lin <- subset(tren_int, waktu_poly == "linear") # apakah laju penurunan linear berbeda?
tren_lin$p.holm <- p.adjust(tren_lin$p.value, "holm") # koreksi Holm untuk 3 perbandingan
tren_lin
## waktu_poly kelompok_pairwise
## 1 linear Kelompok A (Diet & Olahraga) - Kelompok B (Metformin Standar)
## 4 linear Kelompok A (Diet & Olahraga) - (Kelompok C (Metformin + Herbal))
## 7 linear Kelompok B (Metformin Standar) - (Kelompok C (Metformin + Herbal))
## estimate SE df t.ratio p.value p.holm
## 1 2.243333 0.1366578 87 16.41570 2.217873e-28 4.435747e-28
## 4 3.866667 0.1366578 87 28.29452 1.188141e-45 3.564422e-45
## 7 1.623333 0.1366578 87 11.87882 6.643211e-20 6.643211e-20
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
# Tidak mensyaratkan sfierisitas, menampung data hilang (MAR) dan waktu
# pengukuran yang tidak seragam.
lmm1 <- lmer(
pemeriksaan ~ kelompok * waktu + (1 | id),
data = data_long, REML = TRUE
)
summary(lmm1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pemeriksaan ~ kelompok * waktu + (1 | id)
## Data: data_long
##
## REML criterion at convergence: -208.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.30820 -0.57911 0.05473 0.49325 2.12640
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.272599 0.52211
## Residual 0.007932 0.08906
## Number of obs: 360, groups: id, 90
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.75833 0.05523 87.00000 140.460 < 2e-16 ***
## kelompok1 0.49917 0.07811 87.00000 6.390 7.92e-09 ***
## kelompok2 0.02667 0.07811 87.00000 0.341 0.73364
## waktu1 0.70056 0.00813 261.00000 86.169 < 2e-16 ***
## waktu2 0.15611 0.00813 261.00000 19.202 < 2e-16 ***
## waktu3 -0.26389 0.00813 261.00000 -32.459 < 2e-16 ***
## kelompok1:waktu1 -0.34472 0.01150 261.00000 -29.982 < 2e-16 ***
## kelompok2:waktu1 0.03778 0.01150 261.00000 3.286 0.00116 **
## kelompok1:waktu2 -0.07028 0.01150 261.00000 -6.112 3.56e-09 ***
## kelompok2:waktu2 0.01222 0.01150 261.00000 1.063 0.28875
## kelompok1:waktu3 0.15639 0.01150 261.00000 13.602 < 2e-16 ***
## kelompok2:waktu3 -0.03444 0.01150 261.00000 -2.996 0.00300 **
## ---
## 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.000 0.000 0.000
## waktu2 0.000 0.000 0.000 -0.333
## waktu3 0.000 0.000 0.000 -0.333 -0.333
## klmpk1:wkt1 0.000 0.000 0.000 0.000 0.000 0.000
## klmpk2:wkt1 0.000 0.000 0.000 0.000 0.000 0.000 -0.500
## klmpk1:wkt2 0.000 0.000 0.000 0.000 0.000 0.000 -0.333 0.167
## klmpk2:wkt2 0.000 0.000 0.000 0.000 0.000 0.000 0.167 -0.333 -0.500
## klmpk1:wkt3 0.000 0.000 0.000 0.000 0.000 0.000 -0.333 0.167 -0.333
## klmpk2:wkt3 0.000 0.000 0.000 0.000 0.000 0.000 0.167 -0.333 0.167
## klm2:2 klm1:3
## kelompok1
## kelompok2
## waktu1
## waktu2
## waktu3
## klmpk1:wkt1
## klmpk2:wkt1
## klmpk1:wkt2
## klmpk2:wkt2
## klmpk1:wkt3 0.167
## klmpk2:wkt3 -0.333 -0.500
lmm2 <- lmer(pemeriksaan ~ kelompok * waktu + (1 + bulan | id), data = data_long, REML = TRUE)
anova(lmm1, lmm2, refit = FALSE) # uji rasio kemungkinan struktur acak
## Data: data_long
## Models:
## lmm1: pemeriksaan ~ kelompok * waktu + (1 | id)
## lmm2: pemeriksaan ~ kelompok * waktu + (1 + bulan | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 -180.77 -126.36 104.38 -208.77
## lmm2 16 -211.55 -149.38 121.78 -243.55 34.786 2 2.794e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm2, ddf = "Kenward-Roger") # uji F tipe III efek tetap
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 0.2815 0.1407 2 87.00 28.755 2.591e-10 ***
## waktu 30.1249 10.0416 3 185.37 2042.303 < 2.2e-16 ***
## kelompok:waktu 4.3318 0.7220 6 206.40 146.674 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1) # korelasi intrakelas
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.972
## Unadjusted ICC: 0.377
# Diagnostik residual LMM (normalitas & homogenitas)
par(mfrow = c(1, 3))
qqnorm(resid(lmm2), main = "Q-Q residual"); qqline(resid(lmm2))
qqnorm(ranef(lmm2)$id[, 1], main = "Q-Q intersep acak"); qqline(ranef(lmm2)$id[, 1])
plot(fitted(lmm2), resid(lmm2), xlab = "Nilai prediksi", ylab = "Residual",
main = "Residual vs prediksi"); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# (paket 'see' + performance::check_model(lmm2) memberi panel diagnostik lengkap)
# JIka ada missing value
# untuk menunjukkan keunggulan LMM
# contoh jika ada missing value
data_long_miss <- data_long
lmm_miss <- lmer(
pemeriksaan ~ kelompok * waktu + (1 + bulan | id),
data = data_long_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 0.2815 0.1407 2 87.00 28.755 2.591e-10 ***
## waktu 30.1245 10.0415 3 185.37 2042.284 < 2.2e-16 ***
## kelompok:waktu 4.3318 0.7220 6 206.40 146.673 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 4. SIMPAN DATA & SESSION INFO
write.csv(data_wide, "Data_Eksperimen_DM_Wide.csv", row.names = FALSE)
write.csv(data_long, "data_eksperimen_DM_long.csv", row.names = FALSE)