# ============================================================
# TUGAS BIOSTATISTIKA INTERMEDIATE
# NAMA : LUKMAN HAKIM
# NIM : 2611018052
# REPEATED MEASURE ANALYSIS DENGAN R
# Studi: Movingcall - physical activity promotion
# Dataset: SIMULASI berbasis Fischer et al. (2019)
# Outcome: objective MVPA (menit/minggu)
# 3 kelompok x 3 waktu
# ============================================================
# 0. PAKET ---------------------------------------------------
# Paket yang dibutuhkan. Paket yang belum ada akan dipasang otomatis.
required_pkgs <- c(
"dplyr", "tidyr", "ggplot2", "afex", "emmeans", "rstatix",
"car", "effectsize", "lme4", "lmerTest", "performance", "ggpubr"
)
missing_pkgs <- required_pkgs[
!vapply(required_pkgs, requireNamespace, FUN.VALUE = logical(1), quietly = TRUE)
]
## Registered S3 method overwritten by 'lme4':
## method from
## na.action.merMod car
if (length(missing_pkgs) > 0) {
message("Menginstal paket yang belum tersedia: ", paste(missing_pkgs, collapse = ", "))
install.packages(missing_pkgs, dependencies = TRUE)
}
suppressPackageStartupMessages({
invisible(lapply(required_pkgs, library, character.only = TRUE))
})
options(contrasts = c("contr.sum", "contr.poly"))
theme_set(theme_bw(base_size = 12))
# 1. DATA ----------------------------------------------------
# Skrip lama langsung mencari CSV di working directory.
# Versi ini akan memakai file tersebut bila tersedia; jika tidak, RStudio
# akan meminta Anda memilih file CSV secara manual.
default_file <- "Data_Simulasi_MVPA_Movingcall_3x3.csv"
if (file.exists(default_file)) {
data_file <- default_file
} else {
message(
"File '", default_file, "' tidak ditemukan di working directory: ",
getwd(),
"\nSilakan pilih file CSV data MVPA pada jendela yang muncul."
)
data_file <- file.choose()
}
# Deteksi pemisah CSV: koma (,) atau titik koma (;).
first_line <- readLines(data_file, n = 1, warn = FALSE)
sep_detected <- if (grepl(";", first_line, fixed = TRUE)) ";" else ","
dec_detected <- if (sep_detected == ";") "," else "."
dat_wide <- read.table(
data_file,
header = TRUE,
sep = sep_detected,
dec = dec_detected,
check.names = FALSE,
stringsAsFactors = FALSE,
strip.white = TRUE
)
# Bersihkan nama kolom dari spasi yang tidak sengaja.
names(dat_wide) <- trimws(names(dat_wide))
required_cols <- c("id", "kelompok", "Baseline", "M6", "M12")
missing_cols <- setdiff(required_cols, names(dat_wide))
if (length(missing_cols) > 0) {
stop(
"Kolom wajib tidak ditemukan: ", paste(missing_cols, collapse = ", "),
"\nKolom yang terbaca: ", paste(names(dat_wide), collapse = ", "),
"\nFormat yang dibutuhkan: id, kelompok, Baseline, M6, M12"
)
}
# Bersihkan isi variabel kelompok.
dat_wide$kelompok <- trimws(as.character(dat_wide$kelompok))
valid_groups <- c("Control", "Coaching", "Coaching+SMS")
invalid_groups <- setdiff(unique(dat_wide$kelompok), valid_groups)
if (length(invalid_groups) > 0) {
stop(
"Nama kelompok tidak sesuai. Nilai yang tidak dikenali: ",
paste(invalid_groups, collapse = ", "),
"\nGunakan tepat: Control, Coaching, Coaching+SMS"
)
}
# Pastikan outcome numerik. Jika angka tersimpan sebagai karakter dengan koma
# desimal, ubah terlebih dahulu menjadi titik sebelum dikonversi ke numeric.
for (v in c("Baseline", "M6", "M12")) {
if (!is.numeric(dat_wide[[v]])) {
dat_wide[[v]] <- as.numeric(gsub(",", ".", trimws(as.character(dat_wide[[v]])), fixed = TRUE))
}
}
if (anyNA(dat_wide[, c("Baseline", "M6", "M12")])) {
bad_rows <- which(!complete.cases(dat_wide[, c("Baseline", "M6", "M12")]))
stop(
"Ada nilai kosong/tidak numerik pada Baseline, M6, atau M12. Baris bermasalah: ",
paste(bad_rows, collapse = ", ")
)
}
if (anyDuplicated(dat_wide$id)) {
stop("Kolom 'id' harus unik: satu baris per peserta pada format wide.")
}
dat_wide$kelompok <- factor(dat_wide$kelompok, levels = valid_groups)
dat_wide$id <- factor(dat_wide$id)
dat_long <- dat_wide |>
pivot_longer(
cols = c(Baseline, M6, M12),
names_to = "waktu",
values_to = "mvpa"
) |>
mutate(
waktu = factor(waktu, levels = c("Baseline", "M6", "M12")),
bulan = case_when(
waktu == "Baseline" ~ 0,
waktu == "M6" ~ 6,
waktu == "M12" ~ 12
)
)
message("Data berhasil dibaca: ", nrow(dat_wide), " peserta; ", nrow(dat_long), " observasi longitudinal.")
## Data berhasil dibaca: 90 peserta; 270 observasi longitudinal.
print(head(dat_wide))
## id kelompok Baseline M6 M12
## 1 P001 Control 86.0 59.3 54.3
## 2 P002 Control 116.2 122.3 68.0
## 3 P003 Control 82.6 97.1 53.3
## 4 P004 Control 88.2 83.1 111.1
## 5 P005 Control 73.6 57.5 52.5
## 6 P006 Control 95.2 107.8 55.6
print(head(dat_long))
## # A tibble: 6 × 5
## id kelompok waktu mvpa bulan
## <fct> <fct> <fct> <dbl> <dbl>
## 1 P001 Control Baseline 86 0
## 2 P001 Control M6 59.3 6
## 3 P001 Control M12 54.3 12
## 4 P002 Control Baseline 116. 0
## 5 P002 Control M6 122. 6
## 6 P002 Control M12 68 12
str(dat_long)
## tibble [270 × 5] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 90 levels "P001","P002",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ kelompok: Factor w/ 3 levels "Control","Coaching",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ waktu : Factor w/ 3 levels "Baseline","M6",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ mvpa : num [1:270] 86 59.3 54.3 116.2 122.3 ...
## $ bulan : num [1:270] 0 6 12 0 6 12 0 6 12 0 ...
# 2. EKSPLORASI ----------------------------------------------
desk <- dat_long |>
group_by(kelompok, waktu) |>
get_summary_stats(mvpa, type = "mean_sd")
desk
## # A tibble: 9 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Control Baseline mvpa 30 90.0 27.0
## 2 Control M6 mvpa 30 84.9 26.4
## 3 Control M12 mvpa 30 63.9 26.1
## 4 Coaching Baseline mvpa 30 90.0 26.6
## 5 Coaching M6 mvpa 30 116. 20.6
## 6 Coaching M12 mvpa 30 96.9 20.5
## 7 Coaching+SMS Baseline mvpa 30 90 23.3
## 8 Coaching+SMS M6 mvpa 30 118. 23.9
## 9 Coaching+SMS M12 mvpa 30 106. 28.2
p_profil <- ggplot(dat_long, aes(waktu, mvpa, group = kelompok, colour = kelompok)) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", size = 2.5) +
labs(x = "Waktu", y = "MVPA objektif (menit/minggu)", colour = "Kelompok",
title = "Profil rerata MVPA") +
theme(legend.position = "bottom")
p_profil
p_spag <- ggplot(dat_long, aes(waktu, mvpa, group = id)) +
geom_line(alpha = .25) +
stat_summary(aes(group = 1), fun = mean, geom = "line", linewidth = 1.2) +
facet_wrap(~kelompok) +
labs(x = "Waktu", y = "MVPA (menit/minggu)",
title = "Lintasan individu dan rerata kelompok")
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH -----------------------
# Contoh fokus: kelompok Coaching+SMS
d1 <- droplevels(filter(dat_long, kelompok == "Coaching+SMS"))
# 3a. Asumsi
# Outlier
d1 |> group_by(waktu) |> identify_outliers(mvpa)
## [1] waktu id kelompok mvpa bulan is.outlier is.extreme
## <0 rows> (or 0-length row.names)
# Normalitas
d1 |> group_by(waktu) |> shapiro_test(mvpa)
## # A tibble: 3 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Baseline mvpa 0.938 0.0807
## 2 M6 mvpa 0.968 0.479
## 3 M12 mvpa 0.990 0.993
ggpubr::ggqqplot(d1, "mvpa", facet.by = "waktu")
# Sphericity + ANOVA
aov1_rs <- anova_test(data = d1, dv = mvpa, wid = id, within = waktu,
effect.size = "pes")
aov1_rs
## ANOVA Table (type III tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2 58 22.843 4.82e-08 * 0.441
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.998 0.975
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05
## 1 waktu 0.998 2, 57.9 4.94e-08 * 1.072 2.14, 62.17 4.82e-08 *
get_anova_table(aov1_rs, correction = "auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2 58 22.843 4.82e-08 * 0.441
# 3b. afex + effect size
aov1 <- aov_ez(id = "id", dv = "mvpa", data = d1, within = "waktu",
anova_table = list(es = c("ges","pes"), correction = "GG"))
aov1
## Anova Table (Type 3 tests)
##
## Response: mvpa
## Effect df MSE F ges pes p.value
## 1 waktu 2.00, 57.90 267.96 22.84 *** .181 .441 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov1)
## Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): HF eps > 1
## treated as 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 986630 1 39833 29 718.300 < 2.2e-16 ***
## waktu 12221 2 15514 58 22.843 4.824e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.99823 0.97548
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.99823 4.942e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 1.071969 4.824455e-08
eta_squared(aov1, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.44 | [0.28, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 3c. Post hoc
em1 <- emmeans(aov1, ~ waktu)
pairs(em1, adjust = "holm")
## contrast estimate SE df t.ratio p.value
## Baseline - M6 -28.5 4.16 29 -6.856 <0.0001
## Baseline - M12 -15.6 4.20 29 -3.713 0.0017
## M6 - M12 12.9 4.31 29 2.994 0.0056
##
## P value adjustment: holm method for 3 tests
contrast(em1, "trt.vs.ctrl", ref = 1, adjust = "holm")
## contrast estimate SE df t.ratio p.value
## M6 - Baseline 28.5 4.16 29 6.856 <0.0001
## M12 - Baseline 15.6 4.20 29 3.713 0.0009
##
## P value adjustment: holm method for 2 tests
# 3d. Non-parametrik pembanding
friedman_test(d1, mvpa ~ waktu | id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 mvpa 30 18.2 2 0.000112 Friedman test
friedman_effsize(d1, mvpa ~ waktu | id)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 mvpa 30 0.303 Kendall W moderate
d1 |> pairwise_wilcox_test(mvpa ~ waktu, paired = TRUE, p.adjust.method = "holm")
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 mvpa Baseline M6 30 30 8 0.0000000466 1.40e-7 ****
## 2 mvpa Baseline M12 30 30 91 0.00277 5.53e-3 **
## 3 mvpa M6 M12 30 30 352. 0.0122 1.22e-2 *
# 4. MIXED DESIGN ANOVA: KELOMPOK x WAKTU -------------------
# 4a. Asumsi
# Outlier per sel
dat_long |> group_by(kelompok, waktu) |> identify_outliers(mvpa)
## # A tibble: 1 × 7
## kelompok waktu id mvpa bulan is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 Control Baseline P011 24.1 0 TRUE FALSE
# Normalitas per sel
dat_long |> group_by(kelompok, waktu) |> shapiro_test(mvpa)
## # A tibble: 9 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Control Baseline mvpa 0.971 0.580
## 2 Control M6 mvpa 0.987 0.968
## 3 Control M12 mvpa 0.961 0.336
## 4 Coaching Baseline mvpa 0.989 0.988
## 5 Coaching M6 mvpa 0.969 0.510
## 6 Coaching M12 mvpa 0.949 0.164
## 7 Coaching+SMS Baseline mvpa 0.938 0.0807
## 8 Coaching+SMS M6 mvpa 0.968 0.479
## 9 Coaching+SMS M12 mvpa 0.990 0.993
# Homogenitas varians antar kelompok pada tiap waktu
dat_long |> group_by(waktu) |> levene_test(mvpa ~ kelompok)
## # A tibble: 3 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 Baseline 2 87 0.0595 0.942
## 2 M6 2 87 0.783 0.460
## 3 M12 2 87 1.16 0.320
# Box's M untuk homogenitas matriks kovarians
tryCatch(
box_m(dat_wide[, c("Baseline", "M6", "M12")], dat_wide$kelompok),
error = function(e) message("Box's M tidak dapat dihitung: ", e$message)
)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 6.59 0.883 12 Box's M-test for Homogeneity of Covariance Matric…
# 4b. Mixed ANOVA
aov2 <- aov_ez(id = "id", dv = "mvpa", data = dat_long,
between = "kelompok", within = "waktu",
anova_table = list(es = c("ges","pes"), correction = "GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: mvpa
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 1313.67 12.63 *** .171 .225 <.001
## 2 waktu 2.00, 173.64 271.59 33.00 *** .100 .275 <.001
## 3 kelompok:waktu 3.99, 173.64 271.59 15.82 *** .096 .267 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2)
## Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): HF eps > 1
## treated as 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 2444128 1 114289 87 1860.534 < 2.2e-16 ***
## kelompok 33196 2 114289 87 12.635 1.523e-05 ***
## waktu 17889 2 47159 174 33.003 7.056e-13 ***
## kelompok:waktu 17156 4 47159 174 15.825 4.580e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## waktu 0.99794 0.91509
## kelompok:waktu 0.99794 0.91509
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.99794 7.416e-13 ***
## kelompok:waktu 0.99794 4.777e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 1.02135 7.056007e-13
## kelompok:waktu 1.02135 4.579800e-11
eta_squared(aov2, partial = TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.23 | [0.10, 1.00]
## waktu | 0.28 | [0.18, 1.00]
## kelompok:waktu | 0.27 | [0.17, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# Versi rstatix
aov2_rs <- anova_test(data = dat_long, dv = mvpa, wid = id,
between = kelompok, within = waktu,
effect.size = "pes", type = 3)
get_anova_table(aov2_rs, correction = "auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 kelompok 2 87 12.635 1.52e-05 * 0.225
## 2 waktu 2 174 33.003 7.06e-13 * 0.275
## 3 kelompok:waktu 4 174 15.825 4.58e-11 * 0.267
# 4c. Simple effects dan post hoc
joint_tests(aov2, by = "kelompok") # efek waktu pada tiap kelompok
## kelompok = Control:
## model term df1 df2 F.ratio p.value
## waktu 2 87 22.058 <0.0001
##
## kelompok = Coaching:
## model term df1 df2 F.ratio p.value
## waktu 2 87 20.078 <0.0001
##
## kelompok = Coaching+SMS:
## model term df1 df2 F.ratio p.value
## waktu 2 87 22.035 <0.0001
joint_tests(aov2, by = "waktu") # efek kelompok pada tiap waktu
## waktu = Baseline:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.000 1.0000
##
## waktu = M6:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 18.936 <0.0001
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 22.938 <0.0001
em_time <- emmeans(aov2, ~ waktu | kelompok)
contrast(em_time, "trt.vs.ctrl", ref = 1, adjust = "holm")
## kelompok = Control:
## contrast estimate SE df t.ratio p.value
## M6 - Baseline -5.09 4.31 87 -1.181 0.2410
## M12 - Baseline -26.09 4.15 87 -6.281 <0.0001
##
## kelompok = Coaching:
## contrast estimate SE df t.ratio p.value
## M6 - Baseline 26.51 4.31 87 6.148 <0.0001
## M12 - Baseline 6.90 4.15 87 1.662 0.1002
##
## kelompok = Coaching+SMS:
## contrast estimate SE df t.ratio p.value
## M6 - Baseline 28.50 4.31 87 6.611 <0.0001
## M12 - Baseline 15.61 4.15 87 3.757 0.0003
##
## P value adjustment: holm method for 2 tests
em_group <- emmeans(aov2, ~ kelompok | waktu)
pairs(em_group, adjust = "tukey")
## waktu = Baseline:
## contrast estimate SE df t.ratio p.value
## Control - Coaching 0.00000 6.63 87 0.000 1.0000
## Control - (Coaching+SMS) -0.00667 6.63 87 -0.001 1.0000
## Coaching - (Coaching+SMS) -0.00667 6.63 87 -0.001 1.0000
##
## waktu = M6:
## contrast estimate SE df t.ratio p.value
## Control - Coaching -31.59667 6.12 87 -5.159 <0.0001
## Control - (Coaching+SMS) -33.59667 6.12 87 -5.485 <0.0001
## Coaching - (Coaching+SMS) -2.00000 6.12 87 -0.327 0.9430
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Control - Coaching -32.99667 6.50 87 -5.079 <0.0001
## Control - (Coaching+SMS) -41.70667 6.50 87 -6.420 <0.0001
## Coaching - (Coaching+SMS) -8.71000 6.50 87 -1.341 0.3767
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4d. Kontras interaksi: apakah perubahan Baseline -> M12 berbeda antar kelompok
# Hitung perubahan M12 - Baseline pada tiap kelompok, lalu bandingkan perubahan tersebut.
em_change <- contrast(
em_time,
method = list("M12 - Baseline" = c(-1, 0, 1))
)
pairs(em_change, adjust = "holm")
## kelompok = Control:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## kelompok = Coaching:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## kelompok = Coaching+SMS:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
# 5. PEMBANDING: LINEAR MIXED MODEL -------------------------
lmm1 <- lmer(mvpa ~ kelompok * waktu + (1 | id), data = dat_long, REML = TRUE)
lmm2 <- lmer(mvpa ~ kelompok * waktu + (1 + bulan | id), data = dat_long, REML = TRUE,
control = lmerControl(optimizer = "bobyqa"))
## boundary (singular) fit: see help('isSingular')
anova(lmm1, lmm2, refit = FALSE)
## Data: dat_long
## Models:
## lmm1: mvpa ~ kelompok * waktu + (1 | id)
## lmm2: mvpa ~ kelompok * waktu + (1 + bulan | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 11 2406.0 2445.6 -1192 2384.0
## lmm2 13 2409.9 2456.7 -1192 2383.9 0.0552 2 0.9728
anova(lmm2, 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 6846 3423.0 2 87.00 12.635 1.523e-05 ***
## waktu 17889 8944.6 2 115.30 32.826 5.199e-12 ***
## kelompok:waktu 17145 4286.1 4 129.03 15.700 1.712e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.562
## Unadjusted ICC: 0.398
tryCatch(
performance::check_model(lmm2),
error = function(e) message("Plot diagnostik LMM tidak dapat dibuat: ", e$message)
)
# 6. SIMPAN OUTPUT ------------------------------------------
output_dir <- file.path(getwd(), "output_repeated_measure")
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)
write.csv(
desk,
file.path(output_dir, "Ringkasan_Deskriptif_MVPA.csv"),
row.names = FALSE
)
ggsave(
file.path(output_dir, "Profil_MVPA.png"),
p_profil, width = 7, height = 4.5, dpi = 300
)
ggsave(
file.path(output_dir, "Spaghetti_MVPA.png"),
p_spag, width = 9, height = 4.5, dpi = 300
)
message("Selesai. Output tersimpan di: ", normalizePath(output_dir, winslash = "/", mustWork = FALSE))
## Selesai. Output tersimpan di: C:/Users/kesba/Downloads/FERI/ARS/BIOSTATISTIKA/KIM HAKIM/output_repeated_measure