# =============================================================================
# Nama : [Bayu Rosandy]
# NIM : [2611018001]
# =============================================================================
#
# REPEATED MEASURE ANALYSIS DENGAN R
# Contoh terapan: Program pelatihan keselamatan kerja dan kepatuhan SOP
# pada operator produksi (DATA SIMULASI)
#
# Isi:
# 0. Paket & pengaturan
# 1. Simulasi data (format panjang & lebar)
# 2. Eksplorasi data: statistik deskriptif, profile plot, spaghetti plot
# 3. Repeated Measure ANOVA satu arah (within-subject: Waktu)
# 3a. Uji asumsi: outlier, normalitas, sfierisitas (Mauchly)
# 3b. ANOVA + koreksi Greenhouse-Geisser / Huynh-Feldt
# 3c. Pendekatan multivariat (MANOVA) sebagai pembanding
# 3d. Post hoc berpasangan & kontras polinomial (tren)
# 3e. Alternatif nonparametrik: uji Friedman
# 4. Mixed Design ANOVA (between: Kelompok x within: Waktu)
# 4a. Uji asumsi: outlier, normalitas, Levene, Box's M, Mauchly
# 4b. ANOVA campuran + ukuran efek
# 4c. Analisis efek sederhana (simple effects) & post hoc
# 4d. Kontras interaksi (perubahan dari baseline antarkelompok)
# 5. Pembanding: Linear Mixed Model (LMM)
# 6. Menyimpan data & ringkasan hasil
# =============================================================================
# 0. PAKET & PENGATURAN
# install.packages(c("tidyverse","afex","emmeans","rstatix","car",
# "effectsize","lme4","lmerTest","performance","ggpubr","Hmisc"))
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. SIMULASI DATA
# Skenario: 90 operator produksi diacak ke tiga kelompok:
# - Kontrol : briefing keselamatan kerja standar
# - Video : pelatihan keselamatan kerja berbasis video
# - Video+Simulasi : video + simulasi praktik prosedur keselamatan
# Skor kepatuhan SOP (0–100) diukur pada minggu 0, 4, 8, 12.
#
# Struktur korelasi dibuat dari intersep acak + kemiringan acak per subjek
# agar data realistis dan cenderung melanggar asumsi sfierisitas.
set.seed(2027)
n_per <- 30
kel_lab <- c("Kontrol","Video","Video+Simulasi")
minggu <- c(0,4,8,12)
mu <- rbind(
"Kontrol" = c(68,68.5,69,69.2),
"Video" = c(68,74,78,80),
"Video+Simulasi" = c(68,76,82,86)
)
sd_int <- 6
sd_slope <- 0.35
sd_eps <- 3.5
dat_wide <- lapply(seq_along(kel_lab), function(g) {
id <- (g-1)*n_per + seq_len(n_per)
b0 <- rnorm(n_per,0,sd_int)
b1 <- rnorm(n_per,0,sd_slope)
y <- sapply(seq_along(minggu), function(t)
pmin(pmax(mu[g,t] + b0 + b1*minggu[t] + rnorm(n_per,0,sd_eps),0),100))
colnames(y) <- paste0("SOP_M",minggu)
data.frame(
id = sprintf("K%02d%02d",g,seq_len(n_per)),
kelompok = kel_lab[g],
usia = sample(20:55,n_per,replace=TRUE),
jk = sample(c("L","P"),n_per,replace=TRUE,prob=c(.65,.35)),
round(y,1)
)
}) |> bind_rows()
dat_wide$kelompok <- factor(dat_wide$kelompok,levels=kel_lab)
dat_wide$id <- factor(dat_wide$id)
dat_long <- dat_wide |>
pivot_longer(starts_with("SOP_M"),names_to="waktu",values_to="skor") |>
mutate(
waktu=factor(waktu,levels=paste0("SOP_M",minggu),labels=paste0("M",minggu)),
minggu=as.numeric(sub("M","",waktu))
)
head(dat_wide)
## id kelompok usia jk SOP_M0 SOP_M4 SOP_M8 SOP_M12
## 1 K0101 Kontrol 51 L 64.0 66.6 69.5 69.6
## 2 K0102 Kontrol 26 L 62.3 63.8 66.2 67.5
## 3 K0103 Kontrol 40 L 60.9 66.8 61.7 71.3
## 4 K0104 Kontrol 39 P 67.7 71.4 79.7 75.7
## 5 K0105 Kontrol 23 L 77.0 77.3 81.0 86.7
## 6 K0106 Kontrol 21 L 72.4 72.8 72.0 72.3
head(dat_long)
## # A tibble: 6 × 7
## id kelompok usia jk waktu skor minggu
## <fct> <fct> <int> <chr> <fct> <dbl> <dbl>
## 1 K0101 Kontrol 51 L M0 64 0
## 2 K0101 Kontrol 51 L M4 66.6 4
## 3 K0101 Kontrol 51 L M8 69.5 8
## 4 K0101 Kontrol 51 L M12 69.6 12
## 5 K0102 Kontrol 26 L M0 62.3 0
## 6 K0102 Kontrol 26 L M4 63.8 4
str(dat_long)
## tibble [360 × 7] (S3: tbl_df/tbl/data.frame)
## $ id : Factor w/ 90 levels "K0101","K0102",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ kelompok: Factor w/ 3 levels "Kontrol","Video",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ usia : int [1:360] 51 51 51 51 26 26 26 26 40 40 ...
## $ jk : chr [1:360] "L" "L" "L" "L" ...
## $ waktu : Factor w/ 4 levels "M0","M4","M8",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ skor : num [1:360] 64 66.6 69.5 69.6 62.3 63.8 66.2 67.5 60.9 66.8 ...
## $ minggu : num [1:360] 0 4 8 12 0 4 8 12 0 4 ...
# 2. EKSPLORASI DATA
desk <- dat_long |>
group_by(kelompok,waktu) |>
get_summary_stats(skor,type="mean_sd")
desk
## # A tibble: 12 × 6
## kelompok waktu variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Kontrol M0 skor 30 68.3 6.99
## 2 Kontrol M4 skor 30 68.9 7.58
## 3 Kontrol M8 skor 30 69.8 7.10
## 4 Kontrol M12 skor 30 69.9 7.65
## 5 Video M0 skor 30 67.8 5.95
## 6 Video M4 skor 30 75.0 7.72
## 7 Video M8 skor 30 77.4 7.80
## 8 Video M12 skor 30 80.0 7.51
## 9 Video+Simulasi M0 skor 30 69.2 7.02
## 10 Video+Simulasi M4 skor 30 77.0 6.81
## 11 Video+Simulasi M8 skor 30 83.3 6.88
## 12 Video+Simulasi M12 skor 30 86.3 10.0
S <- cov(dat_wide[,paste0("SOP_M",minggu)])
R <- cor(dat_wide[,paste0("SOP_M",minggu)])
round(S,1); round(R,2)
## SOP_M0 SOP_M4 SOP_M8 SOP_M12
## SOP_M0 43.9 34.0 33.1 35.7
## SOP_M4 34.0 65.2 59.2 69.3
## SOP_M8 33.1 59.2 82.5 87.7
## SOP_M12 35.7 69.3 87.7 116.6
## SOP_M0 SOP_M4 SOP_M8 SOP_M12
## SOP_M0 1.00 0.64 0.55 0.50
## SOP_M4 0.64 1.00 0.81 0.79
## SOP_M8 0.55 0.81 1.00 0.89
## SOP_M12 0.50 0.79 0.89 1.00
pasangan <- combn(paste0("SOP_M",minggu),2)
var_selisih <- apply(pasangan,2,function(p)
var(dat_wide[[p[1]]] - dat_wide[[p[2]]]))
names(var_selisih) <- apply(pasangan,2,paste,collapse=" - ")
round(var_selisih,1)
## SOP_M0 - SOP_M4 SOP_M0 - SOP_M8 SOP_M0 - SOP_M12 SOP_M4 - SOP_M8
## 41.1 60.2 89.2 29.3
## SOP_M4 - SOP_M12 SOP_M8 - SOP_M12
## 43.2 23.5
p_profil <- ggplot(dat_long,aes(minggu,skor,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=.6) +
scale_x_continuous(breaks=minggu) +
labs(x="Minggu ke-",y="Skor kepatuhan SOP",colour="Kelompok",
title="Profil rerata skor SOP (±95% CI)") +
theme(legend.position="bottom")
p_profil
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `fun.data()`:
## ! The package "Hmisc" is required.
p_spag <- ggplot(dat_long,aes(minggu,skor,group=id)) +
geom_line(alpha=.3) +
stat_summary(aes(group=kelompok),fun=mean,geom="line",
linewidth=1.2) +
facet_wrap(~kelompok) +
scale_x_continuous(breaks=minggu) +
labs(x="Minggu ke-",y="Skor SOP",title="Lintasan individu dan rerata kelompok")
p_spag
# 3. REPEATED MEASURE ANOVA SATU ARAH
# Pertanyaan: apakah skor kepatuhan SOP berubah selama 12 minggu
# pada kelompok Video+Simulasi?
d1 <- droplevels(filter(dat_long,kelompok=="Video+Simulasi"))
d1w <- filter(dat_wide,kelompok=="Video+Simulasi")
# 3a. Uji asumsi
d1 |> group_by(waktu) |> identify_outliers(skor)
## [1] waktu id kelompok usia jk skor minggu
## [8] is.outlier is.extreme
## <0 rows> (or 0-length row.names)
d1 |> group_by(waktu) |> shapiro_test(skor)
## # A tibble: 4 × 4
## waktu variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M0 skor 0.964 0.395
## 2 M4 skor 0.980 0.832
## 3 M8 skor 0.980 0.831
## 4 M12 skor 0.927 0.0409
ggpubr::ggqqplot(d1,"skor",facet.by="waktu")
aov1_rs <- anova_test(data=d1,dv=skor,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 3 87 110.982 1.22e-29 * 0.793
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 waktu 0.637 0.029 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF]
## 1 waktu 0.796 2.39, 69.27 4.35e-24 * 0.873 2.62, 75.91 3.61e-26
## p[HF]<.05
## 1 *
get_anova_table(aov1_rs,correction="auto")
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 waktu 2.39 69.27 110.982 4.35e-24 * 0.793
# 3b. ANOVA + koreksi GG/HF
aov1 <- aov_ez(id="id",dv="skor",data=d1,within="waktu",
anova_table=list(es=c("ges","pes"),correction="GG"))
aov1
## Anova Table (Type 3 tests)
##
## Response: skor
## Effect df MSE F ges pes p.value
## 1 waktu 2.39, 69.27 19.48 110.98 *** .423 .793 <.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) 747972 1 5708.1 29 3800.07 < 2.2e-16 ***
## waktu 5164 3 1349.4 87 110.98 < 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.63727 0.028747
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.79622 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8725611 3.611093e-26
eta_squared(aov1,partial=TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## waktu | 0.79 | [0.73, 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.41 | [0.27, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
# 3c. Pendekatan multivariat
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.99243 3800.1 1 29 < 2.2e-16 ***
## waktu 1 0.92859 117.0 3 27 1.382e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3d. Post hoc & kontras tren
em1 <- emmeans(aov1,~waktu)
em1
## waktu emmean SE df lower.CL upper.CL
## M0 69.2 1.28 29 66.6 71.8
## M4 77.0 1.24 29 74.5 79.6
## M8 83.3 1.26 29 80.7 85.9
## M12 86.3 1.83 29 82.6 90.1
##
## Confidence level used: 0.95
pairs(em1,adjust="bonferroni")
## contrast estimate SE df t.ratio p.value
## M0 - M4 -7.84 0.951 29 -8.249 <0.0001
## M0 - M8 -14.11 0.732 29 -19.279 <0.0001
## M0 - M12 -17.13 1.260 29 -13.639 <0.0001
## M4 - M8 -6.26 0.907 29 -6.905 <0.0001
## M4 - M12 -9.29 1.180 29 -7.870 <0.0001
## M8 - M12 -3.02 0.986 29 -3.066 0.0280
##
## P value adjustment: bonferroni method for 6 tests
contrast(em1,"trt.vs.ctrl",ref=1,adjust="holm")
## contrast estimate SE df t.ratio p.value
## M4 - M0 7.84 0.951 29 8.249 <0.0001
## M8 - M0 14.11 0.732 29 19.279 <0.0001
## M12 - M0 17.13 1.260 29 13.639 <0.0001
##
## P value adjustment: holm method for 3 tests
contrast(em1,"poly")
## contrast estimate SE df t.ratio p.value
## linear 57.65 3.90 29 14.801 <0.0001
## quadratic -4.82 1.18 29 -4.068 0.0003
## cubic -1.66 2.97 29 -0.559 0.5807
# 3e. Alternatif nonparametrik
friedman_test(d1,skor~waktu|id)
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 skor 30 69.6 3 5.10e-15 Friedman test
friedman_effsize(d1,skor~waktu|id)
## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 skor 30 0.774 Kendall W large
d1 |> wilcox_test(skor~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 skor M0 M4 30 30 9 0.0000000615 3.69e-7 ****
## 2 skor M0 M8 30 30 0 0.00000000186 1.12e-8 ****
## 3 skor M0 M12 30 30 0 0.00000000186 1.12e-8 ****
## 4 skor M4 M8 30 30 15.5 0.000000264 1.59e-6 ****
## 5 skor M4 M12 30 30 16.5 0.000000322 1.93e-6 ****
## 6 skor M8 M12 30 30 95 0.00374 2.24e-2 *
# 4. MIXED DESIGN ANOVA
# Pertanyaan: apakah pola perubahan skor SOP berbeda antarkelompok?
# 4a. Uji asumsi
dat_long |> group_by(kelompok,waktu) |> identify_outliers(skor)
## # A tibble: 2 × 9
## kelompok waktu id usia jk skor minggu is.outlier is.extreme
## <fct> <fct> <fct> <int> <chr> <dbl> <dbl> <lgl> <lgl>
## 1 Kontrol M8 K0120 31 L 86.7 8 TRUE FALSE
## 2 Video M0 K0229 25 L 86 0 TRUE FALSE
dat_long |> group_by(kelompok,waktu) |> shapiro_test(skor)
## # A tibble: 12 × 5
## kelompok waktu variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Kontrol M0 skor 0.934 0.0646
## 2 Kontrol M4 skor 0.981 0.855
## 3 Kontrol M8 skor 0.984 0.914
## 4 Kontrol M12 skor 0.976 0.713
## 5 Video M0 skor 0.939 0.0853
## 6 Video M4 skor 0.983 0.907
## 7 Video M8 skor 0.952 0.190
## 8 Video M12 skor 0.984 0.928
## 9 Video+Simulasi M0 skor 0.964 0.395
## 10 Video+Simulasi M4 skor 0.980 0.832
## 11 Video+Simulasi M8 skor 0.980 0.831
## 12 Video+Simulasi M12 skor 0.927 0.0409
ggpubr::ggqqplot(dat_long,"skor",ggtheme=theme_bw()) +
facet_grid(waktu~kelompok)
dat_long |> group_by(waktu) |> levene_test(skor~kelompok)
## # A tibble: 4 × 5
## waktu df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 M0 2 87 1.71 0.187
## 2 M4 2 87 0.287 0.751
## 3 M8 2 87 0.141 0.869
## 4 M12 2 87 2.98 0.0562
box_m(dat_wide[,paste0("SOP_M",minggu)],dat_wide$kelompok)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 27.1 0.132 20 Box's M-test for Homogeneity of Covariance Matric…
# 4b. ANOVA campuran
aov2 <- aov_ez(id="id",dv="skor",data=dat_long,
between="kelompok",within="waktu",
anova_table=list(es=c("ges","pes"),correction="GG"))
aov2
## Anova Table (Type 3 tests)
##
## Response: skor
## Effect df MSE F ges pes p.value
## 1 kelompok 2, 87 175.89 16.35 *** .228 .273 <.001
## 2 waktu 2.59, 225.42 18.42 114.90 *** .220 .569 <.001
## 3 kelompok:waktu 5.18, 225.42 18.42 23.26 *** .102 .348 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
summary(aov2)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1992596 1 15302.3 87 11328.770 < 2.2e-16 ***
## kelompok 5752 2 15302.3 87 16.351 9.37e-07 ***
## waktu 5485 3 4153.2 261 114.897 < 2.2e-16 ***
## kelompok:waktu 2220 6 4153.2 261 23.257 < 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.79822 0.0016768
## kelompok:waktu 0.79822 0.0016768
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## waktu 0.86368 < 2.2e-16 ***
## kelompok:waktu 0.86368 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## waktu 0.8925231 1.213461e-42
## kelompok:waktu 0.8925231 7.444534e-20
eta_squared(aov2,partial=TRUE)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## kelompok | 0.27 | [0.14, 1.00]
## waktu | 0.57 | [0.51, 1.00]
## kelompok:waktu | 0.35 | [0.27, 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.25 | [0.13, 1.00]
## waktu | 0.22 | [0.14, 1.00]
## kelompok:waktu | 0.10 | [0.03, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
afex_plot(aov2,x="waktu",trace="kelompok",error="within",
mapping=c("colour","shape","linetype")) +
labs(y="Skor SOP",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"
# 4c. Efek sederhana & post hoc
em2 <- emmeans(aov2,~waktu|kelompok)
joint_tests(aov2,by="kelompok")
## Warning in pf(conf$F.ratio, conf$df1, conf$df2, lower.tail = FALSE): NaNs
## produced
## kelompok = Kontrol:
## model term df1 df2 F.ratio p.value
## waktu 3 87 0.767 0.5156
##
## kelompok = Video:
## model term df1 df2 F.ratio p.value
## waktu 3 87 34.597 <0.0001
##
## kelompok = Video+Simulasi:
## model term df1 df2 F.ratio p.value
## waktu 3 87 71.908 <0.0001
joint_tests(aov2,by="waktu")
## waktu = M0:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 0.350 0.7059
##
## waktu = M4:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 9.797 0.0001
##
## waktu = M8:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 25.937 <0.0001
##
## waktu = M12:
## model term df1 df2 F.ratio p.value
## kelompok 2 87 28.793 <0.0001
contrast(em2,"trt.vs.ctrl",ref=1,adjust="holm")
## kelompok = Kontrol:
## contrast estimate SE df t.ratio p.value
## M4 - M0 0.597 1.02 87 0.588 0.5583
## M8 - M0 1.497 1.05 87 1.420 0.4774
## M12 - M0 1.543 1.26 87 1.227 0.4774
##
## kelompok = Video:
## contrast estimate SE df t.ratio p.value
## M4 - M0 7.203 1.02 87 7.094 <0.0001
## M8 - M0 9.687 1.05 87 9.192 <0.0001
## M12 - M0 12.203 1.26 87 9.704 <0.0001
##
## kelompok = Video+Simulasi:
## contrast estimate SE df t.ratio p.value
## M4 - M0 7.843 1.02 87 7.724 <0.0001
## M8 - M0 14.107 1.05 87 13.386 <0.0001
## M12 - M0 17.130 1.26 87 13.621 <0.0001
##
## P value adjustment: holm method for 3 tests
em2b <- emmeans(aov2,~kelompok|waktu)
pairs(em2b,adjust="tukey")
## waktu = M0:
## contrast estimate SE df t.ratio p.value
## Kontrol - Video 0.56 1.72 87 0.325 0.9435
## Kontrol - (Video+Simulasi) -0.87 1.72 87 -0.505 0.8692
## Video - (Video+Simulasi) -1.43 1.72 87 -0.830 0.6856
##
## waktu = M4:
## contrast estimate SE df t.ratio p.value
## Kontrol - Video -6.05 1.91 87 -3.173 0.0059
## Kontrol - (Video+Simulasi) -8.12 1.91 87 -4.259 0.0002
## Video - (Video+Simulasi) -2.07 1.91 87 -1.086 0.5251
##
## waktu = M8:
## contrast estimate SE df t.ratio p.value
## Kontrol - Video -7.63 1.88 87 -4.065 0.0003
## Kontrol - (Video+Simulasi) -13.48 1.88 87 -7.182 <0.0001
## Video - (Video+Simulasi) -5.85 1.88 87 -3.117 0.0069
##
## waktu = M12:
## contrast estimate SE df t.ratio p.value
## Kontrol - Video -10.10 2.19 87 -4.618 <0.0001
## Kontrol - (Video+Simulasi) -16.46 2.19 87 -7.524 <0.0001
## Video - (Video+Simulasi) -6.36 2.19 87 -2.906 0.0128
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# 4d. Kontras interaksi: perubahan M12-M0 antarkelompok
em_full <- emmeans(aov2,~waktu*kelompok)
contrast(
em_full,
interaction=list(waktu=list("M12-M0"=c(-1,0,0,1)),
kelompok="pairwise"),
adjust="holm"
)
## waktu_custom kelompok_pairwise estimate SE df t.ratio p.value
## M12-M0 Kontrol - Video -10.66 1.78 87 -5.994 <0.0001
## M12-M0 Kontrol - (Video+Simulasi) -15.59 1.78 87 -8.764 <0.0001
## M12-M0 Video - (Video+Simulasi) -4.93 1.78 87 -2.770 0.0068
##
## P value adjustment: holm method for 3 tests
contrast(em2,"poly")[c(1,4,7)]
## contrast kelompok estimate SE df t.ratio p.value
## linear Kontrol 5.53 4.04 87 1.370 0.1741
## linear Video 39.09 4.04 87 9.688 <0.0001
## linear Video+Simulasi 57.65 4.04 87 14.288 <0.0001
tren_int <- summary(
contrast(em_full,interaction=c(waktu="poly",kelompok="pairwise"),
adjust="none")
)
tren_lin <- subset(tren_int,waktu_poly=="linear")
tren_lin$p.holm <- p.adjust(tren_lin$p.value,"holm")
tren_lin
## waktu_poly kelompok_pairwise estimate SE df t.ratio
## 1 linear Kontrol - Video -33.56333 5.706644 87 -5.881448
## 4 linear Kontrol - (Video+Simulasi) -52.12333 5.706644 87 -9.133797
## 7 linear Video - (Video+Simulasi) -18.56000 5.706644 87 -3.252349
## p.value p.holm
## 1 7.382172e-08 1.476434e-07
## 4 2.389727e-14 7.169180e-14
## 7 1.629722e-03 1.629722e-03
# 5. PEMBANDING: LINEAR MIXED MODEL (LMM)
lmm1 <- lmer(skor~kelompok*waktu+(1|id),data=dat_long,REML=TRUE)
lmm2 <- lmer(skor~kelompok*waktu+(1+minggu|id),data=dat_long,REML=TRUE)
anova(lmm1,lmm2,refit=FALSE)
## Data: dat_long
## Models:
## lmm1: skor ~ kelompok * waktu + (1 | id)
## lmm2: skor ~ kelompok * waktu + (1 + minggu | id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lmm1 14 2245.5 2299.9 -1108.8 2217.5
## lmm2 16 2224.2 2286.3 -1096.1 2192.2 25.322 2 3.172e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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 381.19 190.60 2 87.00 16.351 9.370e-07 ***
## waktu 2750.82 916.94 3 185.37 78.301 < 2.2e-16 ***
## kelompok:waktu 1125.87 187.64 6 206.40 16.006 4.337e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(lmm1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.715
## Unadjusted ICC: 0.428
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))
# Simulasi 30 nilai hilang (MCAR) untuk menunjukkan keunggulan LMM
set.seed(1)
dat_miss <- dat_long
dat_miss$skor[sample(which(dat_miss$waktu!="M0"),30)] <- NA
lmm_miss <- lmer(skor~kelompok*waktu+(1+minggu|id),data=dat_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 353.43 176.72 2 86.969 16.140 1.093e-06 ***
## waktu 2391.94 797.31 3 168.298 72.476 < 2.2e-16 ***
## kelompok:waktu 935.45 155.91 6 185.950 14.156 2.960e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
n_distinct(dat_miss$id[is.na(dat_miss$skor)])
## [1] 25
# 6. MENYIMPAN DATA & RINGKASAN
write.csv(dat_wide,"data_kasus_baru_wide.csv",row.names=FALSE)
write.csv(dat_long,"data_kasus_baru_long.csv",row.names=FALSE)
saveRDS(list(wide=dat_wide,long=dat_long,anova=aov2,lmm=lmm2),
"hasil_repeated_measure_kasus_baru.rds")