# NAMA : DIMAS WISNU WARDHANA
# NIM : 2611018042
# MATA KULIAH : BIOSTATISTIKA INTERMEDIAT
# REPEATED MEASURE ANALYSIS DENGAN R
# Data simulasi berdasarkan:Algheshairy et al(2025).
# Dataset: 3 kelompok (CG, IG, IIG), pengukuran minggu 0, 4, 8, 12.
# Outcome: Systolic Blood Pressure (SBP; mmHg)
# 1. Paket yang diperlukan --------------------------------------------------
packages <- c("tidyverse", "rstatix", "car", "biotools")
not_installed <- packages[!vapply(packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(not_installed) > 0) {
stop(
"Paket berikut belum terpasang: ", paste(not_installed, collapse = ", "),
". Pasang sekali melalui Console dengan: install.packages(c(",
paste0('"', not_installed, '"', collapse = ", "), "))."
)
}
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(rstatix)
##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filter
# 2. Impor dan siapkan data -------------------------------------------------
data_path <- file.choose()
raw_data <- readr::read_csv(data_path, show_col_types = FALSE)
names(raw_data) <- trimws(names(raw_data))
# Terima ID maupun id.
if (!"ID" %in% names(raw_data) && "id" %in% names(raw_data)) {
raw_data <- dplyr::rename(raw_data, ID = id)
}
if (!all(c("ID", "group") %in% names(raw_data))) {
stop("CSV harus memiliki kolom ID dan group.")
}
# Jika yang dipilih file wide, ubah empat kolom SBP menjadi format long.
wide_sbp <- c("sbp_mmhg_w0", "sbp_mmhg_w4", "sbp_mmhg_w8", "sbp_mmhg_w12")
if (all(c("week", "sbp_mmhg") %in% names(raw_data))) {
data_long <- raw_data %>%
dplyr::transmute(
ID = as.character(ID),
group = as.character(group),
week_num = readr::parse_number(as.character(week)),
sbp_mmhg = as.numeric(sbp_mmhg)
)
} else if (all(wide_sbp %in% names(raw_data))) {
data_long <- raw_data %>%
dplyr::select(ID, group, dplyr::all_of(wide_sbp)) %>%
dplyr::mutate(ID = as.character(ID), group = as.character(group)) %>%
tidyr::pivot_longer(
cols = dplyr::all_of(wide_sbp),
names_to = "week_name",
names_pattern = "sbp_mmhg_w([0-9]+)",
values_to = "sbp_mmhg"
) %>%
dplyr::mutate(
week_num = as.numeric(week_name),
sbp_mmhg = as.numeric(sbp_mmhg)
) %>%
dplyr::select(ID, group, week_num, sbp_mmhg)
} else {
stop("Kolom SBP tidak dikenali. Pilih CSV long atau wide yang disiapkan.")
}
if (anyNA(data_long$week_num)) stop("Kolom week berisi nilai yang tidak dapat dibaca sebagai angka.")
if (anyNA(data_long$sbp_mmhg)) warning("Ada nilai SBP kosong/non-numerik; subjek tidak lengkap tidak masuk RM ANOVA.")
# Urutan level mengikuti kelompok pada dataset dan waktu pengukuran.
group_levels <- c("CG", "IG", "IIG")
observed_groups <- unique(data_long$group)
group_levels <- c(group_levels[group_levels %in% observed_groups],
setdiff(observed_groups, group_levels))
week_levels <- sort(unique(data_long$week_num))
data_long <- data_long %>%
dplyr::mutate(
ID = factor(ID),
group = factor(group, levels = group_levels),
week = factor(week_num, levels = week_levels,
labels = paste0("Minggu ", week_levels))
) %>%
dplyr::arrange(ID, week_num)
if (nlevels(data_long$group) < 2 || nlevels(data_long$week) < 2) {
stop("Analisis mixed ANOVA memerlukan minimal dua kelompok dan dua waktu.")
}
# Cek satu kelompok per subjek dan satu observasi per subjek-waktu.
if (data_long %>% dplyr::distinct(ID, group) %>% dplyr::count(ID) %>% dplyr::filter(n > 1) %>% nrow() > 0) {
stop("Ada ID yang tercatat pada lebih dari satu kelompok.")
}
duplicates <- data_long %>% dplyr::count(ID, week) %>% dplyr::filter(n > 1)
if (nrow(duplicates) > 0) stop("Ada lebih dari satu baris untuk kombinasi ID dan week.")
# RM ANOVA klasik memerlukan semua waktu tersedia untuk setiap subjek.
complete_ids <- data_long %>%
dplyr::group_by(ID) %>%
dplyr::summarise(
n_week = dplyr::n_distinct(week[!is.na(sbp_mmhg)]),
n_valid = sum(!is.na(sbp_mmhg)),
.groups = "drop"
) %>%
dplyr::filter(n_week == nlevels(data_long$week), n_valid == nlevels(data_long$week)) %>%
dplyr::pull(ID)
excluded_ids <- setdiff(levels(data_long$ID), as.character(complete_ids))
if (length(excluded_ids) > 0) {
message("Subjek tidak lengkap yang dikeluarkan dari RM ANOVA: ", paste(excluded_ids, collapse = ", "))
message("Jika banyak data hilang, pertimbangkan linear mixed model (LMM), bukan hanya complete-case RM ANOVA.")
}
analysis_data <- data_long %>%
dplyr::filter(ID %in% complete_ids) %>%
droplevels()
if (nrow(analysis_data) == 0) stop("Tidak ada subjek dengan pengukuran lengkap.")
cat("Jumlah subjek lengkap:", dplyr::n_distinct(analysis_data$ID), "\n")
## Jumlah subjek lengkap: 93
print(analysis_data %>% dplyr::count(group, name = "n_peserta"))
## # A tibble: 3 × 2
## group n_peserta
## <fct> <int>
## 1 CG 128
## 2 IG 124
## 3 IIG 120
print(analysis_data %>% dplyr::count(week, name = "n_pengukuran"))
## # A tibble: 4 × 2
## week n_pengukuran
## <fct> <int>
## 1 Minggu 0 93
## 2 Minggu 4 93
## 3 Minggu 8 93
## 4 Minggu 12 93
# Folder keluaran berada di samping CSV yang dipilih.
results_dir <- file.path(
dirname(data_path),
paste0(tools::file_path_sans_ext(basename(data_path)), "_hasil_R")
)
dir.create(results_dir, showWarnings = FALSE, recursive = TRUE)
# 3. Statistik deskriptif ---------------------------------------------------
descriptives <- analysis_data %>%
dplyr::group_by(group, week) %>%
dplyr::summarise(
n = sum(!is.na(sbp_mmhg)),
mean_sbp = mean(sbp_mmhg, na.rm = TRUE),
sd_sbp = sd(sbp_mmhg, na.rm = TRUE),
se_sbp = sd_sbp / sqrt(n),
.groups = "drop"
)
print(descriptives)
## # A tibble: 12 × 6
## group week n mean_sbp sd_sbp se_sbp
## <fct> <fct> <int> <dbl> <dbl> <dbl>
## 1 CG Minggu 0 32 146. 10.6 1.87
## 2 CG Minggu 4 32 146. 11.3 2.00
## 3 CG Minggu 8 32 146. 10.7 1.90
## 4 CG Minggu 12 32 142. 11.5 2.03
## 5 IG Minggu 0 31 142. 18.7 3.37
## 6 IG Minggu 4 31 138. 18.4 3.30
## 7 IG Minggu 8 31 135. 18.1 3.25
## 8 IG Minggu 12 31 131. 19.2 3.44
## 9 IIG Minggu 0 30 147. 17.7 3.24
## 10 IIG Minggu 4 30 142. 19.3 3.52
## 11 IIG Minggu 8 30 136. 19.0 3.47
## 12 IIG Minggu 12 30 131. 18.0 3.29
readr::write_csv(descriptives, file.path(results_dir, "01_deskriptif.csv"))
# 4. Pemeriksaan outlier dan normalitas -----------------------------------
# Flag outlier saja; jangan menghapus nilai tanpa memeriksa sumber/kesalahan data.
outliers <- analysis_data %>%
dplyr::group_by(group, week) %>%
rstatix::identify_outliers(sbp_mmhg)
print(outliers)
## # A tibble: 7 × 7
## group week ID week_num sbp_mmhg is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 IIG Minggu 0 P035 0 111 TRUE FALSE
## 2 IIG Minggu 0 P045 0 110. TRUE FALSE
## 3 IIG Minggu 4 P035 4 104. TRUE FALSE
## 4 IIG Minggu 4 P045 4 103. TRUE FALSE
## 5 IIG Minggu 8 P035 8 98.6 TRUE FALSE
## 6 IIG Minggu 8 P045 8 99.8 TRUE FALSE
## 7 IIG Minggu 12 P035 12 95.3 TRUE FALSE
readr::write_csv(outliers, file.path(results_dir, "02_outlier.csv"))
normality <- analysis_data %>%
dplyr::group_by(group, week) %>%
rstatix::shapiro_test(sbp_mmhg)
print(normality)
## # A tibble: 12 × 5
## group week variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 CG Minggu 0 sbp_mmhg 0.987 0.964
## 2 CG Minggu 4 sbp_mmhg 0.942 0.0859
## 3 CG Minggu 8 sbp_mmhg 0.963 0.335
## 4 CG Minggu 12 sbp_mmhg 0.975 0.649
## 5 IG Minggu 0 sbp_mmhg 0.977 0.739
## 6 IG Minggu 4 sbp_mmhg 0.977 0.712
## 7 IG Minggu 8 sbp_mmhg 0.979 0.792
## 8 IG Minggu 12 sbp_mmhg 0.979 0.777
## 9 IIG Minggu 0 sbp_mmhg 0.974 0.639
## 10 IIG Minggu 4 sbp_mmhg 0.972 0.605
## 11 IIG Minggu 8 sbp_mmhg 0.969 0.509
## 12 IIG Minggu 12 sbp_mmhg 0.965 0.403
readr::write_csv(normality, file.path(results_dir, "03_shapiro_per_kelompok_waktu.csv"))
p_box <- ggplot(analysis_data, aes(x = week, y = sbp_mmhg, fill = group)) +
geom_boxplot(alpha = 0.75, outlier.shape = 21) +
facet_wrap(~ group) +
labs(x = "Waktu pengukuran", y = "Tekanan darah sistolik (mmHg)",
title = "Sebaran tekanan sistolik per kelompok dan waktu") +
theme_minimal() + theme(legend.position = "none")
ggsave(file.path(results_dir, "01_boxplot.png"), p_box, width = 9, height = 5, dpi = 300)
p_qq <- ggplot(analysis_data, aes(sample = sbp_mmhg)) +
stat_qq() + stat_qq_line() +
facet_grid(group ~ week) +
labs(title = "Q-Q plot tekanan sistolik per kelompok dan waktu") +
theme_minimal()
ggsave(file.path(results_dir, "02_qqplot.png"), p_qq, width = 10, height = 7, dpi = 300)
# 5. Homogenitas varians dan matriks kovarians -----------------------------
# Levene berbasis median pada setiap waktu, sesuai konsep materi.
levene_results <- analysis_data %>%
dplyr::group_by(week) %>%
rstatix::levene_test(sbp_mmhg ~ group, center = median)
print(levene_results)
## # A tibble: 4 × 5
## week df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Minggu 0 2 90 3.80 0.0260
## 2 Minggu 4 2 90 2.98 0.0557
## 3 Minggu 8 2 90 3.40 0.0376
## 4 Minggu 12 2 90 2.71 0.0723
readr::write_csv(levene_results, file.path(results_dir, "04_levene_per_waktu.csv"))
# Box's M pada data wide; gunakan alpha 0,001 karena uji ini sensitif.
box_data <- analysis_data %>%
dplyr::select(ID, group, week_num, sbp_mmhg) %>%
tidyr::pivot_wider(names_from = week_num, values_from = sbp_mmhg,
names_prefix = "W") %>%
tidyr::drop_na()
box_vars <- paste0("W", week_levels)
boxm_result <- biotools::boxM(
data = as.data.frame(box_data %>% dplyr::select(dplyr::all_of(box_vars))),
grouping = box_data$group
)
print(boxm_result)
##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: as.data.frame(box_data %>% dplyr::select(dplyr::all_of(box_vars)))
## Chi-Sq (approx.) = 26.309, df = 20, p-value = 0.1558
boxm_table <- tibble::tibble(
statistic = unname(boxm_result$statistic),
df = unname(boxm_result$parameter),
p_value = boxm_result$p.value,
alpha = 0.001
)
readr::write_csv(boxm_table, file.path(results_dir, "05_box_m.csv"))
# 6. One-way repeated-measures ANOVA (efek waktu saja) ---------------------
# Analisis ini mengabaikan kelompok; hasilnya pelengkap, bukan model utama.
rm_one_way <- rstatix::anova_test(
data = analysis_data,
dv = sbp_mmhg,
wid = ID,
within = week,
type = 3,
effect.size = "pes",
ci = 0.95
)
cat("\nHASIL RM ANOVA SATU ARAH: tabel Mauchly dan epsilon disertakan\n")
##
## HASIL RM ANOVA SATU ARAH: tabel Mauchly dan epsilon disertakan
print(rm_one_way)
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes conf.low conf.high
## 1 week 3 276 150.162 1.04e-57 * 0.62 0.554 0.673
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 week 0.488 1.11e-12 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 week 0.672 2.02, 185.55 1.05e-39 * 0.687 2.06, 189.73 1.55e-40
## p[HF]<.05
## 1 *
rm_one_way_gg <- rstatix::get_anova_table(rm_one_way, correction = "GG")
print(rm_one_way_gg)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes conf.low conf.high
## 1 week 2.02 185.55 150.162 1.05e-39 * 0.62 0.554 0.673
readr::write_csv(rm_one_way_gg, file.path(results_dir, "06_rm_anova_GG.csv"))
# 7. Mixed design ANOVA: kelompok x waktu ----------------------------------
mixed_anova <- rstatix::anova_test(
data = analysis_data,
dv = sbp_mmhg,
wid = ID,
between = group,
within = week,
type = 3,
effect.size = "pes",
ci = 0.95
)
cat("\nHASIL MIXED ANOVA: tabel Mauchly dan epsilon GG/HF disertakan\n")
##
## HASIL MIXED ANOVA: tabel Mauchly dan epsilon GG/HF disertakan
print(mixed_anova)
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes conf.low conf.high
## 1 group 2 90 2.248 1.11e-01 0.048 0.000 0.146
## 2 week 3 270 323.514 4.60e-89 * 0.782 0.742 0.814
## 3 group:week 6 270 51.477 5.19e-42 * 0.534 0.453 0.594
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 week 0.925 0.224
## 2 group:week 0.925 0.224
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 week 0.95 2.85, 256.4 1.00e-84 * 0.984 2.95, 265.66 1.11e-87
## 2 group:week 0.95 5.7, 256.4 5.04e-40 * 0.984 5.9, 265.66 2.24e-41
## p[HF]<.05
## 1 *
## 2 *
# Tabel auto menerapkan GG untuk faktor within yang melanggar sphericity.
# Lihat epsilon GG pada output di atas: materi konsep menganjurkan GG bila
# epsilon < 0,75 dan HF bila epsilon >= 0,75. Tabel GG dan HF juga disimpan
# agar koreksi yang sesuai dapat dipilih untuk masing-masing efek within.
mixed_auto <- rstatix::get_anova_table(mixed_anova, correction = "auto")
mixed_gg <- rstatix::get_anova_table(mixed_anova, correction = "GG")
mixed_hf <- rstatix::get_anova_table(mixed_anova, correction = "HF")
cat("\nTabel mixed ANOVA dengan koreksi otomatis:\n")
##
## Tabel mixed ANOVA dengan koreksi otomatis:
print(mixed_auto)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes conf.low conf.high
## 1 group 2 90 2.248 1.11e-01 0.048 0.000 0.146
## 2 week 3 270 323.514 4.60e-89 * 0.782 0.742 0.814
## 3 group:week 6 270 51.477 5.19e-42 * 0.534 0.453 0.594
cat("\nTabel mixed ANOVA Greenhouse-Geisser:\n")
##
## Tabel mixed ANOVA Greenhouse-Geisser:
print(mixed_gg)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes conf.low conf.high
## 1 group 2.00 90.0 2.248 1.11e-01 0.048 0.000 0.146
## 2 week 2.85 256.4 323.514 1.00e-84 * 0.782 0.742 0.814
## 3 group:week 5.70 256.4 51.477 5.04e-40 * 0.534 0.453 0.594
cat("\nTabel mixed ANOVA Huynh-Feldt:\n")
##
## Tabel mixed ANOVA Huynh-Feldt:
print(mixed_hf)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes conf.low conf.high
## 1 group 2.00 90.00 2.248 1.11e-01 0.048 0.000 0.146
## 2 week 2.95 265.66 323.514 1.11e-87 * 0.782 0.742 0.814
## 3 group:week 5.90 265.66 51.477 2.24e-41 * 0.534 0.453 0.594
readr::write_csv(mixed_auto, file.path(results_dir, "07_mixed_anova_auto.csv"))
readr::write_csv(mixed_gg, file.path(results_dir, "08_mixed_anova_GG.csv"))
readr::write_csv(mixed_hf, file.path(results_dir, "09_mixed_anova_HF.csv"))
# Simpan ringkasan penuh termasuk Mauchly dan sphericity corrections.
capture.output(
{
cat("ONE-WAY REPEATED-MEASURES ANOVA\n")
print(rm_one_way)
cat("\nMIXED DESIGN ANOVA\n")
print(mixed_anova)
},
file = file.path(results_dir, "10_ringkasan_anova.txt")
)
# 8. Uji lanjut / efek sederhana dengan koreksi Holm -----------------------
# Pada interaksi signifikan, fokuskan interpretasi pada dua tabel simple effects:
# (a) beda kelompok pada tiap waktu dan (b) beda waktu di dalam tiap kelompok.
# Bila interaksi tidak signifikan, prioritaskan efek utama; jangan menyimpulkan
# intervensi berbeda hanya dari satu perbandingan simple effect.
posthoc_group_by_week <- analysis_data %>%
dplyr::group_by(week) %>%
rstatix::pairwise_t_test(
sbp_mmhg ~ group,
p.adjust.method = "holm",
pool.sd = FALSE
)
posthoc_week_by_group <- analysis_data %>%
dplyr::group_by(group) %>%
rstatix::pairwise_t_test(
sbp_mmhg ~ week,
paired = TRUE,
id = "ID",
p.adjust.method = "holm"
)
print(posthoc_group_by_week)
## # A tibble: 12 × 11
## week .y. group1 group2 n1 n2 statistic df p p.adj
## * <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Minggu 0 sbp_mmhg CG IG 32 31 1.03 47.1 0.310 0.673
## 2 Minggu 0 sbp_mmhg CG IIG 32 30 -0.477 46.8 0.635 0.673
## 3 Minggu 0 sbp_mmhg IG IIG 31 30 -1.23 59.0 0.224 0.673
## 4 Minggu 4 sbp_mmhg CG IG 32 31 1.98 49.6 0.0529 0.159
## 5 Minggu 4 sbp_mmhg CG IIG 32 30 0.828 46.3 0.412 0.753
## 6 Minggu 4 sbp_mmhg IG IIG 31 30 -0.891 58.6 0.376 0.753
## 7 Minggu 8 sbp_mmhg CG IG 32 31 3.11 48.4 0.00314 0.00942
## 8 Minggu 8 sbp_mmhg CG IIG 32 30 2.48 45.2 0.0168 0.0336
## 9 Minggu 8 sbp_mmhg IG IIG 31 30 -0.397 58.6 0.693 0.693
## 10 Minggu 12 sbp_mmhg CG IG 32 31 2.72 48.8 0.00906 0.0227
## 11 Minggu 12 sbp_mmhg CG IIG 32 30 2.79 48.8 0.00758 0.0227
## 12 Minggu 12 sbp_mmhg IG IIG 31 30 -0.0209 58.9 0.983 0.983
## # ℹ 1 more variable: p.adj.signif <chr>
print(posthoc_week_by_group)
## # A tibble: 18 × 11
## group .y. group1 group2 n1 n2 statistic df p p.adj
## * <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 CG sbp_mmhg Minggu 0 Minggu… 32 32 -0.326 31 7.47e- 1 7.47e- 1
## 2 CG sbp_mmhg Minggu 0 Minggu… 32 32 -1.20 31 2.38e- 1 7.14e- 1
## 3 CG sbp_mmhg Minggu 0 Minggu… 32 32 6.08 31 9.71e- 7 4.86e- 6
## 4 CG sbp_mmhg Minggu 4 Minggu… 32 32 -0.936 31 3.56e- 1 7.14e- 1
## 5 CG sbp_mmhg Minggu 4 Minggu… 32 32 6.05 31 1.05e- 6 4.86e- 6
## 6 CG sbp_mmhg Minggu 8 Minggu… 32 32 7.20 31 4.30e- 8 2.58e- 7
## 7 IG sbp_mmhg Minggu 0 Minggu… 31 31 6.26 30 6.74e- 7 1.35e- 6
## 8 IG sbp_mmhg Minggu 0 Minggu… 31 31 11.3 30 2.22e-12 8.86e-12
## 9 IG sbp_mmhg Minggu 0 Minggu… 31 31 16.7 30 1.05e-16 6.31e-16
## 10 IG sbp_mmhg Minggu 4 Minggu… 31 31 6.41 30 4.41e- 7 1.32e- 6
## 11 IG sbp_mmhg Minggu 4 Minggu… 31 31 12.5 30 1.84e-13 9.21e-13
## 12 IG sbp_mmhg Minggu 8 Minggu… 31 31 4.71 30 5.28e- 5 5.28e- 5
## 13 IIG sbp_mmhg Minggu 0 Minggu… 30 30 8.14 29 5.58e- 9 5.58e- 9
## 14 IIG sbp_mmhg Minggu 0 Minggu… 30 30 20.2 29 1.25e-18 6.27e-18
## 15 IIG sbp_mmhg Minggu 0 Minggu… 30 30 31.7 29 4.53e-24 2.72e-23
## 16 IIG sbp_mmhg Minggu 4 Minggu… 30 30 12.8 29 1.85e-13 5.56e-13
## 17 IIG sbp_mmhg Minggu 4 Minggu… 30 30 16.8 29 1.77e-16 7.08e-16
## 18 IIG sbp_mmhg Minggu 8 Minggu… 30 30 8.73 29 1.30e- 9 2.61e- 9
## # ℹ 1 more variable: p.adj.signif <chr>
readr::write_csv(posthoc_group_by_week, file.path(results_dir, "11_posthoc_kelompok_pada_waktu.csv"))
readr::write_csv(posthoc_week_by_group, file.path(results_dir, "12_posthoc_waktu_dalam_kelompok.csv"))
# 9. Plot interaksi dengan rerata +/- SE -----------------------------------
p_interaction <- ggplot(
descriptives,
aes(x = week, y = mean_sbp, group = group, color = group)
) +
geom_line(linewidth = 0.9) +
geom_point(size = 2.5) +
geom_errorbar(aes(ymin = mean_sbp - se_sbp, ymax = mean_sbp + se_sbp),
width = 0.12) +
labs(
title = "Perubahan tekanan darah sistolik menurut kelompok",
x = "Waktu pengukuran",
y = "Rerata tekanan sistolik (mmHg)",
color = "Kelompok"
) +
theme_minimal() +
theme(legend.position = "bottom")
print(p_interaction)
ggsave(file.path(results_dir, "03_plot_interaksi.png"),
p_interaction, width = 8, height = 5, dpi = 300)
cat("\nSelesai. Tabel dan grafik disimpan di:\n", normalizePath(results_dir), "\n")
##
## Selesai. Tabel dan grafik disimpan di:
## C:\Users\masdi\OneDrive\Desktop\New Biostat\Data_Biostatistika_Wide_hasil_R