paket <- c(
"tidyverse", "afex", "emmeans", "rstatix",
"car", "effectsize", "lme4", "lmerTest",
"performance", "ggpubr", "pbkrtest"
)
belum_terpasang <- setdiff(
paket,
rownames(installed.packages())
)
if (length(belum_terpasang) > 0) {
install.packages(
belum_terpasang,
repos = "https://cloud.r-project.org"
)
}
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(afex)
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Registered S3 method overwritten by 'car':
## method from
## na.action.merMod lme4
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - Get and set global package options with: afex_options()
## - Set sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
##
## Attaching package: 'afex'
##
## The following object is masked from 'package:lme4':
##
## lmer
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(rstatix)
##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filter
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(effectsize)
##
## Attaching package: 'effectsize'
##
## The following objects are masked from 'package:rstatix':
##
## cohens_d, eta_squared, omega_squared
library(lme4)
library(lmerTest)
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
library(performance)
library(ggpubr)
options(
contrasts = c("contr.sum", "contr.poly")
)
afex_options(
emmeans_model = "multivariate"
)
theme_set(
theme_bw(base_size = 12)
)
dat_wide <- read.csv(
"data_tds_hipertensi_wide_dummy.csv",
stringsAsFactors = FALSE
)
head(dat_wide)
## id kelompok usia jk TDS_M0 TDS_M4 TDS_M8 TDS_M12
## 1 P001 Kontrol 47 L 148.9 147.7 146.2 138.9
## 2 P002 Kontrol 43 L 155.0 149.4 157.9 152.6
## 3 P003 Kontrol 64 P 153.2 152.1 151.4 145.9
## 4 P004 Kontrol 61 P 146.4 141.8 144.8 132.5
## 5 P005 Kontrol 64 P 133.4 146.3 143.6 141.7
## 6 P006 Kontrol 53 P 158.2 161.5 140.5 147.8
str(dat_wide)
## 'data.frame': 90 obs. of 8 variables:
## $ id : chr "P001" "P002" "P003" "P004" ...
## $ kelompok: chr "Kontrol" "Kontrol" "Kontrol" "Kontrol" ...
## $ usia : int 47 43 64 61 64 53 60 64 46 65 ...
## $ jk : chr "L" "L" "P" "P" ...
## $ TDS_M0 : num 149 155 153 146 133 ...
## $ TDS_M4 : num 148 149 152 142 146 ...
## $ TDS_M8 : num 146 158 151 145 144 ...
## $ TDS_M12 : num 139 153 146 132 142 ...
names(dat_wide)
## [1] "id" "kelompok" "usia" "jk" "TDS_M0" "TDS_M4" "TDS_M8"
## [8] "TDS_M12"
dim(dat_wide)
## [1] 90 8
table(dat_wide$kelompok)
##
## DASH DASH+AF Kontrol
## 30 30 30
colSums(is.na(dat_wide))
## id kelompok usia jk TDS_M0 TDS_M4 TDS_M8 TDS_M12
## 0 0 0 0 0 0 0 0
anyDuplicated(dat_wide$id)
## [1] 0
kel_lab <- c("Kontrol", "DASH", "DASH+AF")
minggu <- c(0, 4, 8, 12)
dat_wide <- dat_wide %>%
mutate(
id = factor(id),
kelompok = factor(
kelompok,
levels = kel_lab
),
jk = factor(
jk,
levels = c("L", "P")
)
)
dat_long <- dat_wide %>%
pivot_longer(
cols = c(
TDS_M0, TDS_M4,
TDS_M8, TDS_M12
),
names_to = "waktu",
values_to = "tds"
) %>%
mutate(
waktu = factor(
waktu,
levels = c(
"TDS_M0", "TDS_M4",
"TDS_M8", "TDS_M12"
),
labels = c("M0", "M4", "M8", "M12")
),
minggu = as.numeric(
sub("M", "", as.character(waktu))
)
)
head(dat_long)
## # A tibble: 6 × 7
## id kelompok usia jk waktu tds minggu
## <fct> <fct> <int> <fct> <fct> <dbl> <dbl>
## 1 P001 Kontrol 47 L M0 149. 0
## 2 P001 Kontrol 47 L M4 148. 4
## 3 P001 Kontrol 47 L M8 146. 8
## 4 P001 Kontrol 47 L M12 139. 12
## 5 P002 Kontrol 43 L M0 155 0
## 6 P002 Kontrol 43 L M4 149. 4
nrow(dat_long)
## [1] 360
table(
dat_long$kelompok,
dat_long$waktu
)
##
## M0 M4 M8 M12
## Kontrol 30 30 30 30
## DASH 30 30 30 30
## DASH+AF 30 30 30 30
sum(is.na(dat_long$tds))
## [1] 0
deskriptif <- dat_long %>%
group_by(kelompok, waktu) %>%
summarise(
n = n(),
rerata = mean(tds),
sd = sd(tds),
minimum = min(tds),
maksimum = max(tds),
.groups = "drop"
)
print(deskriptif, n = 12)
## # A tibble: 12 × 7
## kelompok waktu n rerata sd minimum maksimum
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Kontrol M0 30 151. 9.71 133. 173.
## 2 Kontrol M4 30 149. 11.7 128 169.
## 3 Kontrol M8 30 148. 9.79 133 176.
## 4 Kontrol M12 30 146. 10.7 126. 164.
## 5 DASH M0 30 151. 9.54 135. 173.
## 6 DASH M4 30 147. 10.7 129. 166.
## 7 DASH M8 30 144. 12.6 122. 170.
## 8 DASH M12 30 141. 12.1 124. 170.
## 9 DASH+AF M0 30 151. 10.4 129. 171.
## 10 DASH+AF M4 30 144. 13.2 116. 170.
## 11 DASH+AF M8 30 141. 10.5 118. 160.
## 12 DASH+AF M12 30 136. 11.7 105. 158.
print(deskriptif, n = 12)
## # A tibble: 12 × 7
## kelompok waktu n rerata sd minimum maksimum
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Kontrol M0 30 151. 9.71 133. 173.
## 2 Kontrol M4 30 149. 11.7 128 169.
## 3 Kontrol M8 30 148. 9.79 133 176.
## 4 Kontrol M12 30 146. 10.7 126. 164.
## 5 DASH M0 30 151. 9.54 135. 173.
## 6 DASH M4 30 147. 10.7 129. 166.
## 7 DASH M8 30 144. 12.6 122. 170.
## 8 DASH M12 30 141. 12.1 124. 170.
## 9 DASH+AF M0 30 151. 10.4 129. 171.
## 10 DASH+AF M4 30 144. 13.2 116. 170.
## 11 DASH+AF M8 30 141. 10.5 118. 160.
## 12 DASH+AF M12 30 136. 11.7 105. 158.
# Hitung interval kepercayaan 95%
deskriptif <- deskriptif %>%
mutate(
minggu = as.numeric(
sub("M", "", as.character(waktu))
),
se = sd / sqrt(n),
batas_bawah = rerata - qt(0.975, df = n - 1) * se,
batas_atas = rerata + qt(0.975, df = n - 1) * se
)
# Buat grafik
grafik_tds <- ggplot(
deskriptif,
aes(
x = minggu,
y = rerata,
colour = kelompok,
group = kelompok
)
) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
geom_errorbar(
aes(
ymin = batas_bawah,
ymax = batas_atas
),
width = 0.5
) +
scale_x_continuous(
breaks = c(0, 4, 8, 12)
) +
scale_colour_manual(
values = c(
"Kontrol" = "#64748B",
"DASH" = "#2563EB",
"DASH+AF" = "#0B6B4F"
)
) +
labs(
title = "Perubahan Rerata Tekanan Darah Sistolik",
subtitle = "Data dummy; garis vertikal = interval kepercayaan 95%",
x = "Minggu ke-",
y = "Rerata TDS (mmHg)",
colour = "Kelompok"
) +
theme_bw(base_size = 12) +
theme(
legend.position = "bottom"
)
print(grafik_tds)
hasil_outlier <- dat_long %>%
group_by(kelompok, waktu) %>%
rstatix::identify_outliers(tds) %>%
ungroup()
print(hasil_outlier, n = Inf)
## # A tibble: 1 × 9
## kelompok waktu id usia jk tds minggu is.outlier is.extreme
## <fct> <fct> <fct> <int> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 DASH+AF M12 P090 42 P 105. 12 TRUE FALSE
hasil_normalitas <- dat_long %>%
group_by(kelompok, waktu) %>%
rstatix::shapiro_test(tds) %>%
ungroup()
print(hasil_normalitas, n = 12)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 tds 0.980 0.822
## 2 Kontrol M4 tds 0.969 0.521
## 3 Kontrol M8 tds 0.964 0.385
## 4 Kontrol M12 tds 0.969 0.510
## 5 DASH M0 tds 0.976 0.713
## 6 DASH M4 tds 0.963 0.367
## 7 DASH M8 tds 0.966 0.429
## 8 DASH M12 tds 0.945 0.122
## 9 DASH+AF M0 tds 0.984 0.918
## 10 DASH+AF M4 tds 0.987 0.963
## 11 DASH+AF M8 tds 0.984 0.918
## 12 DASH+AF M12 tds 0.949 0.156
grafik_qq <- ggplot(
dat_long,
aes(sample = tds)
) +
stat_qq(size = 1.5, alpha = 0.7) +
stat_qq_line(colour = "#0B6B4F") +
facet_grid(waktu ~ kelompok) +
labs(
title = "Pemeriksaan Normalitas TDS",
x = "Kuantil teoretis",
y = "Kuantil data"
) +
theme_bw(base_size = 11)
print(grafik_qq)
d1 <- dat_long %>%
filter(kelompok == "DASH+AF") %>%
droplevels()
table(d1$waktu)
##
## M0 M4 M8 M12
## 30 30 30 30
anova_satu <- rstatix::anova_test(
data = d1,
dv = tds,
wid = id,
within = waktu,
effect.size = "pes"
)
print(anova_satu)
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 3 87 51.037 4e-19 * 0.638
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.812 0.33
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.868 2.6, 75.53 6.6e-17 * 0.961 2.88, 83.64 1.78e-18
## p[HF]<.05
## 1 *
aov1 <- afex::aov_ez(
id = "id",
dv = "tds",
data = d1,
within = "waktu",
anova_table = list(
es = c("ges", "pes"),
correction = "none"
)
)
# Rerata estimasi pada setiap waktu
em1 <- emmeans::emmeans(
aov1,
~ waktu
)
print(em1)
## waktu emmean SE df lower.CL upper.CL
## M0 151 1.89 29 147 155
## M4 144 2.41 29 139 149
## M8 141 1.92 29 137 144
## M12 136 2.13 29 132 141
##
## Confidence level used: 0.95
# Bandingkan seluruh pasangan waktu
posthoc_waktu <- pairs(
em1,
adjust = "bonferroni"
)
print(posthoc_waktu)
## contrast estimate SE df t.ratio p.value
## M0 - M4 6.57 1.08 29 6.080 <0.0001
## M0 - M8 10.36 1.13 29 9.136 <0.0001
## M0 - M12 14.61 1.36 29 10.709 <0.0001
## M4 - M8 3.79 1.21 29 3.128 0.0239
## M4 - M12 8.04 1.44 29 5.593 <0.0001
## M8 - M12 4.26 1.09 29 3.907 0.0031
##
## P value adjustment: bonferroni method for 6 tests
hasil_levene <- dat_long %>%
group_by(waktu) %>%
rstatix::levene_test(tds ~ kelompok) %>%
ungroup()
print(hasil_levene)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 0.175 0.840
## 2 M4 2 87 0.283 0.754
## 3 M8 2 87 1.71 0.186
## 4 M12 2 87 0.482 0.619
hasil_boxm <- rstatix::box_m(
dat_wide[, c(
"TDS_M0", "TDS_M4",
"TDS_M8", "TDS_M12"
)],
dat_wide$kelompok
)
print(hasil_boxm)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 20.5 0.430 20 Box's M-test for Homogeneity of Covariance Matric…
anova_campuran <- rstatix::anova_test(
data = dat_long,
dv = tds,
wid = id,
between = kelompok,
within = waktu,
effect.size = "pes"
)
print(anova_campuran)
## ANOVA Table (type II tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 kelompok 2 87 2.137 1.24e-01 0.047
## 2 waktu 3 261 69.270 5.54e-33 * 0.443
## 3 kelompok:waktu 6 261 5.570 1.88e-05 * 0.114
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.877 0.047 *
## 2 kelompok:waktu 0.877 0.047 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF]
## 1 waktu 0.922 2.77, 240.68 1.28e-30 * 0.955 2.87, 249.37
## 2 kelompok:waktu 0.922 5.53, 240.68 3.58e-05 * 0.955 5.73, 249.37
## p[HF] p[HF]<.05
## 1 1.25e-31 *
## 2 2.72e-05 *
tabel_anova_campuran <- rstatix::get_anova_table(
anova_campuran,
correction = "auto"
)
print(tabel_anova_campuran)
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 kelompok 2.00 87.00 2.137 1.24e-01 0.047
## 2 waktu 2.77 240.68 69.270 1.28e-30 * 0.443
## 3 kelompok:waktu 5.53 240.68 5.570 3.58e-05 * 0.114
aov2 <- afex::aov_ez(
id = "id",
dv = "tds",
data = dat_long,
between = "kelompok",
within = "waktu",
type = 3,
anova_table = list(
es = c("ges", "pes"),
correction = "GG"
)
)
print(aov2)
## Anova Table (Type 3 tests)
##
## Response: tds
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 423.12 2.14 .040 .047 .124
## 2 waktu 2.77, 240.68 25.64 69.27 *** .103 .443 <.001
## 3 kelompok:waktu 5.53, 240.68 25.64 5.57 *** .018 .114 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
# Rerata estimasi kelompok pada setiap waktu
em_kelompok <- emmeans::emmeans(
aov2,
~ kelompok | waktu
)
# Perbandingan pasangan kelompok
posthoc_kelompok <- pairs(
em_kelompok,
adjust = "tukey"
)
print(posthoc_kelompok)
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Kontrol - DASH -0.197 2.55 87 -0.077 0.9967
## Kontrol - (DASH+AF) 0.100 2.55 87 0.039 0.9992
## DASH - (DASH+AF) 0.297 2.55 87 0.116 0.9926
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - DASH 1.497 3.08 87 0.486 0.8781
## Kontrol - (DASH+AF) 4.580 3.08 87 1.487 0.3020
## DASH - (DASH+AF) 3.083 3.08 87 1.001 0.5781
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - DASH 4.370 2.85 87 1.534 0.2803
## Kontrol - (DASH+AF) 7.750 2.85 87 2.720 0.0213
## DASH - (DASH+AF) 3.380 2.85 87 1.186 0.4645
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - DASH 5.100 2.97 87 1.717 0.2047
## Kontrol - (DASH+AF) 9.527 2.97 87 3.207 0.0053
## DASH - (DASH+AF) 4.427 2.97 87 1.490 0.3006
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# Perubahan M12 - M0 dalam setiap kelompok
em_waktu_kelompok <- emmeans::emmeans(
aov2,
~ waktu | kelompok
)
perubahan <- emmeans::contrast(
em_waktu_kelompok,
method = list(
"M12 - M0" = c(-1, 0, 0, 1)
)
)
print(perubahan)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M12 - M0 -5.19 1.42 87 -3.641 0.0005
##
## kelompok = DASH:
## contrast estimate SE df t.ratio p.value
## M12 - M0 -10.48 1.42 87 -7.359 <0.0001
##
## kelompok = DASH+AF:
## contrast estimate SE df t.ratio p.value
## M12 - M0 -14.61 1.42 87 -10.259 <0.0001
# Bandingkan perubahan tersebut antarkelompok
beda_perubahan <- pairs(
perubahan,
by = NULL,
adjust = "holm"
)
print(beda_perubahan)
## contrast estimate SE df t.ratio p.value
## (M12 - M0 Kontrol) - (M12 - M0 DASH) 5.30 2.01 87 2.629 0.0202
## (M12 - M0 Kontrol) - (M12 - M0 DASH+AF) 9.43 2.01 87 4.679 <0.0001
## (M12 - M0 DASH) - (M12 - M0 DASH+AF) 4.13 2.01 87 2.050 0.0434
##
## P value adjustment: holm method for 3 tests
perubahan_baseline <- emmeans::contrast(
em_waktu_kelompok,
method = "trt.vs.ctrl",
ref = 1,
adjust = "holm"
)
print(perubahan_baseline)
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -2.09 1.08 87 -1.940 0.0732
## M8 - M0 -2.71 1.27 87 -2.123 0.0732
## M12 - M0 -5.19 1.42 87 -3.641 0.0014
##
## kelompok = DASH:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -3.78 1.08 87 -3.512 0.0007
## M8 - M0 -7.27 1.27 87 -5.705 <0.0001
## M12 - M0 -10.48 1.42 87 -7.359 <0.0001
##
## kelompok = DASH+AF:
## contrast estimate SE df t.ratio p.value
## M4 - M0 -6.57 1.08 87 -6.099 <0.0001
## M8 - M0 -10.36 1.27 87 -8.124 <0.0001
## M12 - M0 -14.61 1.42 87 -10.259 <0.0001
##
## P value adjustment: holm method for 3 tests
# Model 1: intersep acak
lmm1 <- lmerTest::lmer(
tds ~ kelompok * waktu + (1 | id),
data = dat_long,
REML = TRUE
)
# Model 2: intersep dan kemiringan acak
lmm2 <- lmerTest::lmer(
tds ~ 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.0165734 (tol = 0.002, component 1)
## See ?lme4::convergence and ?lme4::troubleshooting.
dat_long <- dat_long %>%
mutate(
minggu_skala = (minggu - 6) / 4
)
lmm2 <- lmerTest::lmer(
tds ~ kelompok * waktu +
(1 + minggu_skala | id),
data = dat_long,
REML = TRUE,
control = lme4::lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 200000)
)
)
# Pesan pemeriksaan konvergensi
lmm2@optinfo$conv$lme4$messages
## NULL
# Kode optimizer
lmm2@optinfo$conv$opt
## [1] 0
# Pemeriksaan singularitas
lme4::isSingular(lmm2)
## [1] FALSE
hasil_lmm <- anova(
lmm2,
type = 3,
ddf = "Kenward-Roger"
)
print(hasil_lmm)
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 82.66 41.33 2 87.00 2.1366 0.1242076
## waktu 2963.30 987.77 3 185.37 50.8295 < 2.2e-16 ***
## kelompok:waktu 489.70 81.62 6 206.40 4.1952 0.0005271 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simpan pengaturan tampilan
pengaturan_lama <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
# 1. Normalitas residual
qqnorm(
resid(lmm2),
main = "Q-Q residual"
)
qqline(resid(lmm2), col = "red")
# 2. Residual terhadap nilai prediksi
plot(
fitted(lmm2),
resid(lmm2),
xlab = "TDS prediksi",
ylab = "Residual",
main = "Residual vs prediksi"
)
abline(h = 0, col = "red", lty = 2)
# 3. Normalitas intersep acak
efek_acak <- lme4::ranef(lmm2)$id
qqnorm(
efek_acak[, "(Intercept)"],
main = "Q-Q intersep acak"
)
qqline(
efek_acak[, "(Intercept)"],
col = "red"
)
# 4. Normalitas kemiringan acak
qqnorm(
efek_acak[, "minggu_skala"],
main = "Q-Q kemiringan acak"
)
qqline(
efek_acak[, "minggu_skala"],
col = "red"
)
# Pulihkan pengaturan
par(pengaturan_lama)
perbandingan_model <- anova(
lmm1,
lmm2,
refit = FALSE
)
print(perbandingan_model)
## Data: dat_long
## Models:
## lmm1: tds ~ kelompok * waktu + (1 | id)
## lmm2: tds ~ kelompok * waktu + (1 + minggu_skala | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 2425.2 2479.6 -1198.6 2397.2
## lmm2 16 2419.4 2481.6 -1193.7 2387.4 9.7383 2 0.00768 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hasil_icc <- performance::icc(lmm1)
print(hasil_icc)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.809
## Unadjusted ICC: 0.691
efek_waktu_per_kelompok <- emmeans::joint_tests(
aov2,
by = "kelompok"
)
print(efek_waktu_per_kelompok)
## kelompok = Kontrol:
## model term df1 df2 F.ratio p.value
## waktu 3 87 4.421 0.0061
##
## kelompok = DASH:
## model term df1 df2 F.ratio p.value
## waktu 3 87 19.128 <0.0001
##
## kelompok = DASH+AF:
## model term df1 df2 F.ratio p.value
## waktu 3 87 37.185 <0.0001
efek_kelompok_per_waktu <- emmeans::joint_tests(
aov2,
by = "waktu"
)
print(efek_kelompok_per_waktu)
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.007 0.9930
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 1.150 0.3214
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 3.719 0.0282
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 5.151 0.0077
efek_waktu_holm <- as.data.frame(
efek_waktu_per_kelompok
) %>%
mutate(
p_holm = p.adjust(p.value, method = "holm")
)
efek_kelompok_holm <- as.data.frame(
efek_kelompok_per_waktu
) %>%
mutate(
p_holm = p.adjust(p.value, method = "holm")
)
print(efek_waktu_holm)
## model term kelompok df1 df2 F.ratio p.value p_holm
## 1 waktu Kontrol 3 87 4.421 6.092613e-03 6.092613e-03
## 4 waktu DASH 3 87 19.128 1.294689e-09 2.589378e-09
## 7 waktu DASH+AF 3 87 37.185 1.463446e-15 4.390337e-15
print(efek_kelompok_holm)
## model term waktu df1 df2 F.ratio p.value p_holm
## 1 kelompok M0 2 87 0.007 0.993023152 0.99302315
## 3 kelompok M4 2 87 1.150 0.321388935 0.64277787
## 5 kelompok M8 2 87 3.719 0.028190942 0.08457283
## 7 kelompok M12 2 87 5.151 0.007690895 0.03076358
tren_dash_af <- emmeans::contrast(
em1,
method = "poly"
)
print(tren_dash_af)
## contrast estimate SE df t.ratio p.value
## linear -47.63 4.61 29 -10.339 <0.0001
## quadratic 2.31 1.54 29 1.500 0.1445
## cubic -3.25 3.47 29 -0.937 0.3565
# Uji multivariat untuk perubahan waktu pada DASH+AF
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.99457 5311.3 1 29 < 2.2e-16 ***
## waktu 1 0.81515 39.7 3 27 4.903e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji multivariat kelompok, waktu, dan interaksi
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.99521 18080.1 1 87 < 2.2e-16 ***
## kelompok 2 0.04682 2.1 2 87 0.124208
## waktu 1 0.64223 50.9 3 85 < 2.2e-16 ***
## kelompok:waktu 2 0.23328 3.8 6 172 0.001445 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji perubahan antarwaktu pada DASH+AF
hasil_friedman <- rstatix::friedman_test(
d1,
tds ~ waktu | id
)
print(hasil_friedman)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 tds 30 59.4 3 8.03e-13 Friedman test
# Ukuran efek Kendall's W
efek_friedman <- rstatix::friedman_effsize(
d1,
tds ~ waktu | id
)
print(efek_friedman)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 tds 30 0.660 Kendall W large
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.99457 5311.3 1 29 < 2.2e-16 ***
## waktu 1 0.81515 39.7 3 27 4.903e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Pastikan urutan responden konsisten pada setiap waktu
d1_urut <- d1 %>%
arrange(waktu, id)
posthoc_wilcoxon <- d1_urut %>%
rstatix::wilcox_test(
tds ~ waktu,
paired = TRUE,
p.adjust.method = "bonferroni"
)
print(posthoc_wilcoxon, n = 6)
## # 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 tds M0 M4 30 30 440 0.00000168 1.01e-5 ****
## 2 tds M0 M8 30 30 461 0.0000000112 6.71e-8 ****
## 3 tds M0 M12 30 30 465 0.00000000186 1.12e-8 ****
## 4 tds M4 M8 30 30 365 0.00510 3.06e-2 *
## 5 tds M4 M12 30 30 436 0.00000324 1.94e-5 ****
## 6 tds M8 M12 30 30 399 0.000305 1.83e-3 **
print(hasil_friedman)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 tds 30 59.4 3 8.03e-13 Friedman test
kolom_tds <- c(
"TDS_M0", "TDS_M4",
"TDS_M8", "TDS_M12"
)
# Matriks kovarians
matriks_kovarians <- cov(
dat_wide[, kolom_tds]
)
print(round(matriks_kovarians, 2))
## TDS_M0 TDS_M4 TDS_M8 TDS_M12
## TDS_M0 95.43 100.35 83.50 82.80
## TDS_M4 100.35 142.75 111.59 118.61
## TDS_M8 83.50 111.59 129.22 114.94
## TDS_M12 82.80 118.61 114.94 144.72
# Matriks korelasi
matriks_korelasi <- cor(
dat_wide[, kolom_tds]
)
print(round(matriks_korelasi, 2))
## TDS_M0 TDS_M4 TDS_M8 TDS_M12
## TDS_M0 1.00 0.86 0.75 0.70
## TDS_M4 0.86 1.00 0.82 0.83
## TDS_M8 0.75 0.82 1.00 0.84
## TDS_M12 0.70 0.83 0.84 1.00
pasangan <- combn(kolom_tds, 2)
varians_selisih <- apply(
pasangan,
2,
function(x) {
var(dat_wide[[x[1]]] - dat_wide[[x[2]]])
}
)
names(varians_selisih) <- apply(
pasangan,
2,
paste,
collapse = " - "
)
print(round(varians_selisih, 2))
## TDS_M0 - TDS_M4 TDS_M0 - TDS_M8 TDS_M0 - TDS_M12 TDS_M4 - TDS_M8
## 37.48 57.65 74.56 48.79
## TDS_M4 - TDS_M12 TDS_M8 - TDS_M12
## 50.26 44.07
grafik_individu <- ggplot(
dat_long,
aes(x = minggu, y = tds, group = id)
) +
geom_line(alpha = 0.25, colour = "grey40") +
stat_summary(
aes(group = kelompok),
fun = mean,
geom = "line",
colour = "#0B6B4F",
linewidth = 1.2
) +
facet_wrap(~ kelompok) +
scale_x_continuous(breaks = c(0, 4, 8, 12)) +
labs(
title = "Lintasan TDS Individu dan Rerata Kelompok",
x = "Minggu ke-",
y = "TDS (mmHg)"
) +
theme_bw()
print(grafik_individu)
em_full <- emmeans::emmeans(
aov2,
~ waktu * kelompok
)
tren_interaksi <- as.data.frame(
emmeans::contrast(
em_full,
interaction = c(
waktu = "poly",
kelompok = "pairwise"
),
adjust = "none"
)
)
tren_linear <- tren_interaksi %>%
filter(waktu_poly == "linear") %>%
mutate(
p_holm = p.adjust(p.value, method = "holm")
)
print(tren_linear)
## waktu_poly kelompok_pairwise estimate SE df t.ratio p.value
## 1 linear Kontrol - DASH 18.76333 6.555447 87 2.862251 5.268628e-03
## 2 linear Kontrol - (DASH+AF) 31.45000 6.555447 87 4.797537 6.599342e-06
## 3 linear DASH - (DASH+AF) 12.68667 6.555447 87 1.935286 5.620348e-02
## p_holm
## 1 1.053726e-02
## 2 1.979803e-05
## 3 5.620348e-02
set.seed(1)
dat_miss <- dat_long
baris_hilang <- sample(
which(dat_miss$waktu != "M0"),
size = 30
)
dat_miss$tds[baris_hilang] <- NA
# Jumlah pengukuran tersedia
sum(!is.na(dat_miss$tds))
## [1] 330
# Jumlah responden yang memiliki data tidak lengkap
dplyr::n_distinct(
dat_miss$id[is.na(dat_miss$tds)]
)
## [1] 25
lmm_miss <- lmerTest::lmer(
tds ~ kelompok * waktu +
(1 + minggu_skala | id),
data = dat_miss,
REML = TRUE,
na.action = na.omit,
control = lme4::lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 200000)
)
)
# Periksa model
lmm_miss@optinfo$conv$lme4$messages
## NULL
lme4::isSingular(lmm_miss)
## [1] FALSE
# Jumlah pengukuran yang digunakan
nobs(lmm_miss)
## [1] 330
# Uji efek pada data tidak lengkap
hasil_lmm_missing <- anova(
lmm_miss,
type = 3,
ddf = "Kenward-Roger"
)
print(hasil_lmm_missing)
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## kelompok 73.58 36.79 2 86.989 1.8978 0.1560568
## waktu 3021.58 1007.19 3 168.124 51.7084 < 2.2e-16 ***
## kelompok:waktu 502.35 83.72 6 185.771 4.2934 0.0004445 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dir.create(
"hasil_analisis",
showWarnings = FALSE
)
# Simpan tabel utama sebagai CSV
tabel_simpan <- list(
deskriptif = deskriptif,
outlier = hasil_outlier,
normalitas = hasil_normalitas,
levene = hasil_levene,
box_m = hasil_boxm,
anova_campuran = tabel_anova_campuran,
posthoc_waktu = as.data.frame(posthoc_waktu),
posthoc_kelompok = as.data.frame(posthoc_kelompok),
beda_perubahan = as.data.frame(beda_perubahan),
perubahan_baseline = as.data.frame(perubahan_baseline),
efek_waktu_holm = efek_waktu_holm,
efek_kelompok_holm = efek_kelompok_holm,
tren_linear = tren_linear,
friedman = hasil_friedman,
efek_friedman = efek_friedman,
posthoc_wilcoxon = posthoc_wilcoxon
)
for (nama in names(tabel_simpan)) {
write.csv(
tabel_simpan[[nama]],
file = file.path(
"hasil_analisis",
paste0(nama, ".csv")
),
row.names = FALSE
)
}
# Simpan output model dan informasi sesi
capture.output(
print(anova_satu),
print(aov2),
print(aov1$Anova),
print(aov2$Anova),
print(tren_dash_af),
print(hasil_lmm),
print(perbandingan_model),
print(hasil_icc),
print(hasil_lmm_missing),
print(matriks_kovarians),
print(matriks_korelasi),
print(varians_selisih),
sessionInfo(),
file = "hasil_analisis/output_model.txt"
)
# Simpan grafik
ggsave(
"hasil_analisis/grafik_rerata_tds.png",
plot = grafik_tds,
width = 9, height = 6, dpi = 300
)
ggsave(
"hasil_analisis/grafik_qq.png",
plot = grafik_qq,
width = 10, height = 10, dpi = 300
)
ggsave(
"hasil_analisis/grafik_individu.png",
plot = grafik_individu,
width = 12, height = 5, dpi = 300
)
# Simpan objek analisis agar bisa dibuka kembali
save.image(
file = "hasil_analisis/analisis_tds.RData"
)
# Lihat lokasi dan daftar hasil
getwd()
## [1] "/Volumes/DATA NAZMY/MKM/Rstudio/Tugas R_anova"
list.files("hasil_analisis")
## [1] "analisis_tds.RData" "anova_campuran.csv" "beda_perubahan.csv"
## [4] "box_m.csv" "deskriptif.csv" "efek_friedman.csv"
## [7] "efek_kelompok_holm.csv" "efek_waktu_holm.csv" "friedman.csv"
## [10] "grafik_individu.png" "grafik_qq.png" "grafik_rerata_tds.png"
## [13] "levene.csv" "normalitas.csv" "outlier.csv"
## [16] "output_model.txt" "perubahan_baseline.csv" "posthoc_kelompok.csv"
## [19] "posthoc_waktu.csv" "posthoc_wilcoxon.csv" "tren_linear.csv"
install.packages(
c("rmarkdown", "knitr"),
repos = "https://cloud.r-project.org"
)
## package 'knitr' successfully unpacked and SHA256 sums checked
##
## The downloaded binary packages are in
## /var/folders/lj/l50w2ndj6wv7tfq7_9th1gg00000gp/T//Rtmp47zmG8/downloaded_packages