# REPEATED MEASURE ANALYSIS DENGAN R
# TUGAS BIOSTATISTIKA INTERMEDIATE
# NAMA : FERI FEBRIAN ARISTIA
# NIM : 2611018048
# Rpubs : ferifebrian
# Data simulasi berdasarkan Ding et al. (2024), PLOS Global Public Health
# Desain: 3 treatment x 3 waktu; outcome = berat badan (kg)
# 0. PAKET & PENGATURAN ------------------------------------------------------
# Paket yang belum terpasang akan dipasang otomatis (butuh internet saat pertama kali)
paket <- c("dplyr", "tidyr", "ggplot2", "afex", "emmeans", "rstatix", "car",
"effectsize", "lme4", "lmerTest", "pbkrtest", "performance", "ggpubr")
baru <- paket[!paket %in% rownames(installed.packages())]
if (length(baru) > 0) install.packages(baru)
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))
# 1. DATA -------------------------------------------------------------------
# Skrip membaca file CSV di working directory (cek dengan getwd()).
# Kolom yang dibutuhkan: id, treatment, weight_baseline, weight_6m, weight_2y
# Nilai treatment: "Control", "Basic MCP", "Full MCP".
# Jika file CSV tidak ditemukan, data simulasi dibuat otomatis (seed tetap) dan
# disimpan sebagai CSV tersebut, sehingga seluruh skrip tetap bisa dijalankan.
file_data <- "Data_Simulasi_Berat_Badan_MCP_3x3.csv"
if (file.exists(file_data)) {
dat_wide <- read.csv(file_data, stringsAsFactors = FALSE)
message("Data dibaca dari: ", file_data)
} else {
message("File ", file_data, " tidak ditemukan -> membuat data simulasi.")
set.seed(2024)
n_per <- 40
grp <- c("Control", "Basic MCP", "Full MCP")
# Perubahan rerata berat badan (kg) terhadap baseline: 6 bulan, 2 tahun
# (nilai ASUMSI untuk simulasi; ganti dengan data/parameter Anda bila perlu)
delta <- rbind("Control" = c(0.2, 0.5),
"Basic MCP" = c(-1.5, -1.2),
"Full MCP" = c(-3.0, -2.5))
dat_wide <- bind_rows(lapply(seq_along(grp), function(g) {
b0 <- rnorm(n_per, 0, 8) # perbedaan berat badan dasar antarpasien
b1 <- rnorm(n_per, 0, 1.2) # perbedaan respons antarpasien
data.frame(
id = sprintf("S%03d", (g - 1) * n_per + seq_len(n_per)),
treatment = grp[g],
weight_baseline = round(78 + b0 + rnorm(n_per, 0, 1.0), 1),
weight_6m = round(78 + delta[g, 1] + b0 + b1 + rnorm(n_per, 0, 1.0), 1),
weight_2y = round(78 + delta[g, 2] + b0 + 1.5 * b1 + rnorm(n_per, 0, 1.2), 1)
)
}))
write.csv(dat_wide, file_data, row.names = FALSE)
}
## Data dibaca dari: Data_Simulasi_Berat_Badan_MCP_3x3.csv
# Pemeriksaan struktur data
kol_perlu <- c("id", "treatment", "weight_baseline", "weight_6m", "weight_2y")
stopifnot(all(kol_perlu %in% names(dat_wide)))
print(colSums(is.na(dat_wide[, kol_perlu])))
## id treatment weight_baseline weight_6m weight_2y
## 0 0 0 0 0
dat_wide$treatment <- factor(dat_wide$treatment,
levels = c("Control", "Basic MCP", "Full MCP"))
stopifnot(!any(is.na(dat_wide$treatment))) # berhenti jika ada label treatment yang tidak cocok
dat_wide$id <- factor(dat_wide$id)
dat_long <- dat_wide |>
pivot_longer(c(weight_baseline, weight_6m, weight_2y),
names_to = "time", values_to = "weight") |>
mutate(time = factor(time,
levels = c("weight_baseline", "weight_6m", "weight_2y"),
labels = c("Baseline", "6 months", "2 years")))
head(dat_wide); head(dat_long)
## id treatment weight_baseline weight_6m weight_2y
## 1 S001 Control 84.6 84.7 86.9
## 2 S002 Control 82.7 81.6 82.2
## 3 S003 Control 78.0 76.6 78.5
## 4 S004 Control 74.5 78.7 82.2
## 5 S005 Control 87.4 85.3 85.2
## 6 S006 Control 89.4 88.7 88.3
## # A tibble: 6 × 4
## id treatment time weight
## <fct> <fct> <fct> <dbl>
## 1 S001 Control Baseline 84.6
## 2 S001 Control 6 months 84.7
## 3 S001 Control 2 years 86.9
## 4 S002 Control Baseline 82.7
## 5 S002 Control 6 months 81.6
## 6 S002 Control 2 years 82.2
# 2. EKSPLORASI DATA ---------------------------------------------------------
desk <- dat_long |> group_by(treatment, time) |>
get_summary_stats(weight, type = "mean_sd")
print(desk)
## # A tibble: 9 × 6
## treatment time variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Control Baseline weight 40 76.4 8.58
## 2 Control 6 months weight 40 76.7 8.30
## 3 Control 2 years weight 40 77.5 8.52
## 4 Basic MCP Baseline weight 40 78.7 6.40
## 5 Basic MCP 6 months weight 40 76.8 6.30
## 6 Basic MCP 2 years weight 40 77.2 6.28
## 7 Full MCP Baseline weight 40 79.4 8.50
## 8 Full MCP 6 months weight 40 76.5 8.83
## 9 Full MCP 2 years weight 40 77.1 8.85
p_profile <- ggplot(dat_long, aes(time, weight, group = treatment, colour = treatment)) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", size = 2.5) +
labs(x = "Time", y = "Body weight (kg)", colour = "Treatment",
title = "Mean body weight profile")
print(p_profile)
p_spag <- ggplot(dat_long, aes(time, weight, group = id)) +
geom_line(alpha = .25) +
stat_summary(aes(group = treatment), fun = mean, geom = "line",
colour = "firebrick", linewidth = 1.2) +
facet_wrap(~ treatment) + labs(x = "Time", y = "Body weight (kg)")
print(p_spag)
# 3. REPEATED MEASURE ANOVA SATU ARAH: FULL MCP ----------------------------
d1 <- droplevels(filter(dat_long, treatment == "Full MCP"))
print(d1 |> group_by(time) |> identify_outliers(weight))
## [1] time id treatment weight is.outlier is.extreme
## <0 rows> (or 0-length row.names)
print(d1 |> group_by(time) |> shapiro_test(weight))
## # A tibble: 3 × 4
## time variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline weight 0.967 0.290
## 2 6 months weight 0.960 0.171
## 3 2 years weight 0.981 0.726
print(ggpubr::ggqqplot(d1, "weight", facet.by = "time"))
aov1 <- aov_ez(id = "id", dv = "weight", data = d1, within = "time",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
print(aov1); print(summary(aov1))
## Anova Table (Type 3 tests)
##
## Response: weight
## Effect df MSE F ges pes p.value
## 1 time 1.61, 62.66 2.31 52.38 *** .021 .573 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 724071 1 8769.3 39 3220.168 < 2.2e-16 ***
## time 194 2 144.7 78 52.382 3.785e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.75527 0.0048295
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.80338 1.257e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.8326476 5.29292e-13
em1 <- emmeans(aov1, ~ time)
print(pairs(em1, adjust = "holm"))
## contrast estimate SE df t.ratio p.value
## Baseline - X6.months 2.93 0.277 39 10.575 <0.0001
## Baseline - X2.years 2.38 0.371 39 6.414 <0.0001
## X6.months - X2.years -0.55 0.252 39 -2.184 0.0350
##
## P value adjustment: holm method for 3 tests
print(friedman_test(d1, weight ~ time | id))
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 weight 40 45.0 2 1.66e-10 Friedman test
# 4. MIXED DESIGN ANOVA: TREATMENT x TIME ----------------------------------
# 4a. Asumsi
print(dat_long |> group_by(treatment, time) |> shapiro_test(weight))
## # A tibble: 9 × 5
## treatment time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Control Baseline weight 0.958 0.141
## 2 Control 6 months weight 0.960 0.168
## 3 Control 2 years weight 0.948 0.0629
## 4 Basic MCP Baseline weight 0.935 0.0242
## 5 Basic MCP 6 months weight 0.933 0.0209
## 6 Basic MCP 2 years weight 0.931 0.0168
## 7 Full MCP Baseline weight 0.967 0.290
## 8 Full MCP 6 months weight 0.960 0.171
## 9 Full MCP 2 years weight 0.981 0.726
print(dat_long |> group_by(time) |> levene_test(weight ~ treatment))
## # A tibble: 3 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 117 2.73 0.0693
## 2 6 months 2 117 2.55 0.0824
## 3 2 years 2 117 2.84 0.0624
print(box_m(dat_wide[, c("weight_baseline", "weight_6m", "weight_2y")], dat_wide$treatment))
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 8.82 0.718 12 Box's M-test for Homogeneity of Covariance Matric…
# 4b. Mixed RM ANOVA
aov2 <- aov_ez(id = "id", dv = "weight", data = dat_long,
between = "treatment", within = "time",
anova_table = list(es = c("ges", "pes"), correction = "GG"))
print(aov2); print(summary(aov2))
## Anova Table (Type 3 tests)
##
## Response: weight
## Effect df MSE F ges pes p.value
## 1 treatment 2, 117 184.24 0.13 .002 .002 .875
## 2 time 1.61, 188.54 2.26 37.40 *** .006 .242 <.001
## 3 treatment:time 3.22, 188.54 2.26 22.21 *** .007 .275 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 2154631 1 21555.8 117 11694.8496 < 2.2e-16 ***
## treatment 49 2 21555.8 117 0.1337 0.875
## time 136 2 426.8 234 37.4044 8.027e-15 ***
## treatment:time 162 4 426.8 234 22.2104 1.471e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.75886 1.1206e-07
## treatment:time 0.75886 1.1206e-07
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.80571 2.053e-12 ***
## treatment:time 0.80571 5.573e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.8152804 1.562204e-12
## treatment:time 0.8152804 4.158068e-13
print(eta_squared(aov2, partial = TRUE))
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## treatment | 2.28e-03 | [0.00, 1.00]
## time | 0.24 | [0.17, 1.00]
## treatment:time | 0.28 | [0.19, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 4c. Simple effects dan post hoc
print(joint_tests(aov2, by = "treatment"))
## treatment = Control:
## model term df1 df2 F.ratio p.value
## time 2 117 5.967 0.0034
##
## treatment = Basic MCP:
## model term df1 df2 F.ratio p.value
## time 2 117 25.292 <0.0001
##
## treatment = Full MCP:
## model term df1 df2 F.ratio p.value
## time 2 117 60.238 <0.0001
print(joint_tests(aov2, by = "time"))
## time = Baseline:
## model term df1 df2 F.ratio p.value
## treatment 2 117 1.656 0.1954
##
## time = X6.months:
## model term df1 df2 F.ratio p.value
## treatment 2 117 0.012 0.9880
##
## time = X2.years:
## model term df1 df2 F.ratio p.value
## treatment 2 117 0.029 0.9713
em2 <- emmeans(aov2, ~ time | treatment)
print(pairs(em2, adjust = "holm"))
## treatment = Control:
## contrast estimate SE df t.ratio p.value
## Baseline - X6.months -0.312 0.272 117 -1.150 0.2524
## Baseline - X2.years -1.133 0.368 117 -3.075 0.0052
## X6.months - X2.years -0.820 0.253 117 -3.237 0.0047
##
## treatment = Basic MCP:
## contrast estimate SE df t.ratio p.value
## Baseline - X6.months 1.877 0.272 117 6.911 <0.0001
## Baseline - X2.years 1.423 0.368 117 3.863 0.0004
## X6.months - X2.years -0.455 0.253 117 -1.796 0.0750
##
## treatment = Full MCP:
## contrast estimate SE df t.ratio p.value
## Baseline - X6.months 2.933 0.272 117 10.794 <0.0001
## Baseline - X2.years 2.382 0.368 117 6.469 <0.0001
## X6.months - X2.years -0.550 0.253 117 -2.171 0.0319
##
## P value adjustment: holm method for 3 tests
# 4d. Kontras perubahan dari baseline antarkelompok
dat_change <- dat_wide |>
mutate(change_6m = weight_6m - weight_baseline,
change_2y = weight_2y - weight_baseline)
print(anova_test(data = dat_change, dv = change_6m, between = treatment))
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 treatment 2 117 37.124 3.28e-13 * 0.388
print(anova_test(data = dat_change, dv = change_2y, between = treatment))
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 treatment 2 117 24.337 1.45e-09 * 0.294
print(pairwise_t_test(dat_change, change_6m ~ treatment, p.adjust.method = "holm"))
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <chr> <dbl> <chr>
## 1 change_6m Control Basic… 40 40 9.11e- 8 **** 1.82e- 7 ****
## 2 change_6m Control Full … 40 40 9.48e-14 **** 2.84e-13 ****
## 3 change_6m Basic MCP Full … 40 40 6.99e- 3 ** 6.99e- 3 **
print(pairwise_t_test(dat_change, change_2y ~ treatment, p.adjust.method = "holm"))
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <chr> <dbl> <chr>
## 1 change_2y Control Basic … 40 40 3.04e- 6 **** 6.07e-6 ****
## 2 change_2y Control Full M… 40 40 6.04e-10 **** 1.81e-9 ****
## 3 change_2y Basic MCP Full M… 40 40 6.78e- 2 ns 6.78e-2 ns
# 5. PEMBANDING: LINEAR MIXED MODEL -----------------------------------------
lmm <- lmer(weight ~ treatment * time + (1 | id), data = dat_long, REML = TRUE)
print(anova(lmm, ddf = "Kenward-Roger"))
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 0.488 0.244 2 117 0.1337 0.875
## time 136.446 68.223 2 234 37.4044 8.027e-15 ***
## treatment:time 162.041 40.510 4 234 22.2104 1.471e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(lmm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ treatment * time + (1 | id)
## Data: dat_long
##
## REML criterion at convergence: 1793.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.56858 -0.45134 0.01732 0.48765 2.37786
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 60.805 7.798
## Residual 1.824 1.351
## Number of obs: 360, groups: id, 120
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 77.3633 0.7154 117.0000 108.143 < 2e-16 ***
## treatment1 -0.5192 1.0117 117.0000 -0.513 0.608806
## treatment2 0.2042 1.0117 117.0000 0.202 0.840420
## time1 0.7967 0.1007 234.0000 7.914 9.91e-14 ***
## time2 -0.7025 0.1007 234.0000 -6.979 3.04e-11 ***
## treatment1:time1 -1.2783 0.1424 234.0000 -8.980 < 2e-16 ***
## treatment2:time1 0.3033 0.1424 234.0000 2.131 0.034151 *
## treatment1:time2 0.5333 0.1424 234.0000 3.746 0.000226 ***
## treatment2:time2 -0.0750 0.1424 234.0000 -0.527 0.598804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn1 trtmn2 time1 time2 trt1:1 trt2:1 trt1:2
## treatment1 0.000
## treatment2 0.000 -0.500
## time1 0.000 0.000 0.000
## time2 0.000 0.000 0.000 -0.500
## trtmnt1:tm1 0.000 0.000 0.000 0.000 0.000
## trtmnt2:tm1 0.000 0.000 0.000 0.000 0.000 -0.500
## trtmnt1:tm2 0.000 0.000 0.000 0.000 0.000 -0.500 0.250
## trtmnt2:tm2 0.000 0.000 0.000 0.000 0.000 0.250 -0.500 -0.500
print(performance::icc(lmm))
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.971
## Unadjusted ICC: 0.956
# 6. SIMPAN RINGKASAN --------------------------------------------------------
write.csv(desk, "Ringkasan_Deskriptif_MCP.csv", row.names = FALSE)
message("Selesai. Ringkasan disimpan di: ", file.path(getwd(), "Ringkasan_Deskriptif_MCP.csv"))
## Selesai. Ringkasan disimpan di: C:/Users/kesba/Downloads/FERI/ARS/BIOSTATISTIKA/Feri Febrian Arisrtia/Feri Febrian Arisrtia/Ringkasan_Deskriptif_MCP.csv