library(tidyverse)
library(car)
library(rstatix)
library(biotools)
library(psych)
df <- read.csv("hsb2.csv")
glimpse(df)
## Rows: 200
## Columns: 11
## $ id <int> 70, 121, 86, 141, 172, 113, 50, 11, 84, 48, 75, 60, 95, 104, 3…
## $ female <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ race <int> 4, 4, 4, 4, 4, 4, 3, 1, 4, 3, 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4,…
## $ ses <int> 1, 2, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 1, 1, 3, 2, 3, 2, 2,…
## $ schtyp <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1,…
## $ prog <int> 1, 3, 1, 3, 2, 2, 1, 2, 1, 2, 3, 2, 2, 2, 2, 1, 2, 1, 2, 1, 1,…
## $ read <int> 57, 68, 44, 63, 47, 44, 50, 34, 63, 57, 60, 57, 73, 54, 45, 42…
## $ write <int> 52, 59, 33, 44, 52, 52, 59, 46, 57, 55, 46, 65, 60, 63, 57, 49…
## $ math <int> 41, 53, 54, 47, 57, 51, 42, 45, 54, 52, 51, 51, 71, 57, 50, 43…
## $ science <int> 47, 63, 58, 53, 53, 63, 53, 39, 58, 50, 53, 63, 61, 55, 31, 50…
## $ socst <int> 57, 61, 31, 56, 61, 61, 61, 36, 51, 51, 61, 61, 71, 46, 56, 56…
names(df)
## [1] "id" "female" "race" "ses" "schtyp" "prog" "read"
## [8] "write" "math" "science" "socst"
memastikan apakah dalam data tersebut terdapat missing value
colSums(is.na(df))
## id female race ses schtyp prog read write math science
## 0 0 0 0 0 0 0 0 0 0
## socst
## 0
dv <- c("read", "write", "math")
iv <- "prog"
covar <- "science"
X <- na.omit(df[, dv])
mardia_result <- psych::mardia(X)
mardia_result
## Call: psych::mardia(x = X)
##
## Mardia tests of multivariate skew and kurtosis
## Use describe(x) the to get univariate tests
## n.obs = 200 num.vars = 3
## b1p = 0.52 skew = 17.25 with probability <= 0.069
## small sample skew = 17.64 with probability <= 0.061
## b2p = 13.22 kurtosis = -2.3 with probability <= 0.021
pvalue skewness= 0.061 dengan p walue kurtosis = 0.021 (salah satunya memenuhi)
boxM(df[, dv], df$prog)
##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: df[, dv]
## Chi-Sq (approx.) = 16.574, df = 12, p-value = 0.1663
p-value = 0.1663 > 0.05 (memenuhi asusmsi)
R <- cor(X)
bartlett_result <- cortest.bartlett(
R,
n = nrow(X)
)
bartlett_result
## $chisq
## [1] 229.6047
##
## $p.value
## [1] 1.683784e-49
##
## $df
## [1] 3
p-value (1.683784e-49) < 0.05 artinya memenuhi asumsi
par(mfrow=c(1,3))
for (y in dv) {
plot(
df[[covar]],
df[[y]],
xlab = "science",
ylab = y,
main = paste("science vs", y)
)
abline(
lm(df[[y]] ~ df[[covar]]),
col = "red"
)
}
manova_model <- manova(cbind(read, write, math) ~ prog, data = df)
summary(manova_model, test = "Pillai")
## Df Pillai approx F num Df den Df Pr(>F)
## prog 1 0.035319 2.392 3 196 0.06984 .
## Residuals 198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P-value = 2.304e-10 < 0.05 yang artinya terdapat perbedaan signifikan pada nilai read, write, dan math berdasarkan prog.
summary.aov(manova_model)
## Response read :
## Df Sum Sq Mean Sq F value Pr(>F)
## prog 1 381.1 381.10 3.674 0.05671 .
## Residuals 198 20538.3 103.73
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response write :
## Df Sum Sq Mean Sq F value Pr(>F)
## prog 1 586.4 586.42 6.7146 0.01027 *
## Residuals 198 17292.5 87.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response math :
## Df Sum Sq Mean Sq F value Pr(>F)
## prog 1 393.5 393.53 4.564 0.03388 *
## Residuals 198 17072.3 86.22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Read p = 4.283e-09 < 0.05 (signifikan)
Write p = 4.31e-09 < 0.05, (signifikan)
Math p = 7.364e-12 < 0.05, (signifikan)
Hasil uji univariat menunjukkan bahwa seluruh variabel dependen, yaitu read, write, dan math, memiliki perbedaan yang signifikan berdasarkan kelompok prog.
slope_model <- manova(cbind(read, write, math) ~ prog * science, data = df)
summary(slope_model, test = "Pillai")
## Df Pillai approx F num Df den Df Pr(>F)
## prog 1 0.06205 4.278 3 194 0.005967 **
## science 1 0.48618 61.188 3 194 < 2.2e-16 ***
## prog:science 1 0.00739 0.482 3 194 0.695366
## Residuals 196
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
prog : science = p = 0.3466 > 0.05 (tidak signifikan) Tidak terdapat interaksi yang signifikan antara prog dan science terhadap variabel dependen.
mancova_model <- manova(cbind(read, write, math) ~ prog + science, data = df)
summary(mancova_model, test = "Pillai")
## Df Pillai approx F num Df den Df Pr(>F)
## prog 1 0.06198 4.295 3 195 0.005828 **
## science 1 0.48602 61.463 3 195 < 2.2e-16 ***
## Residuals 197
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mancova_model, test = “Pillai”)
summary.aov(mancova_model)
## Response read :
## Df Sum Sq Mean Sq F value Pr(>F)
## prog 1 381.1 381.1 5.9551 0.01556 *
## science 1 7931.1 7931.1 123.9309 < 2e-16 ***
## Residuals 197 12607.2 64.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response write :
## Df Sum Sq Mean Sq F value Pr(>F)
## prog 1 586.4 586.4 9.6572 0.002165 **
## science 1 5329.8 5329.8 87.7703 < 2.2e-16 ***
## Residuals 197 11962.7 60.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response math :
## Df Sum Sq Mean Sq F value Pr(>F)
## prog 1 393.5 393.5 7.3828 0.007172 **
## science 1 6571.5 6571.5 123.2850 < 2.2e-16 ***
## Residuals 197 10500.8 53.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
read Variabel read dipengaruhi secara signifikan oleh prog (F = 33.195; p = 3.84e-13) dan science (F = 111.267; p < 2.2e-16).
write Variabel write dipengaruhi secara signifikan oleh prog (F = 29.327; p = 7.209e-12) dan science (F = 75.564; p = 1.410e-15).
math Variabel math dipengaruhi secara signifikan oleh prog (F = 46.009; p < 2.2e-16) dan science (F = 113.565; p < 2.2e-16).
# ANOVA (tanpa covariate)
anova_read <- aov(read ~ prog, data = df)
anova_write <- aov(write ~ prog, data = df)
anova_math <- aov(math ~ prog, data = df)
# ANCOVA (dengan covariate science)
ancova_read <- aov(read ~ science + prog, data = df)
ancova_write <- aov(write ~ science + prog, data = df)
ancova_math <- aov(math ~ science + prog, data = df)
hasil_aa <- data.frame(
DV = c("read", "write", "math"),
F_ANOVA = c(
summary(anova_read)[[1]][1, "F value"],
summary(anova_write)[[1]][1, "F value"],
summary(anova_math)[[1]][1, "F value"]
),
p_ANOVA = c(
summary(anova_read)[[1]][1, "Pr(>F)"],
summary(anova_write)[[1]][1, "Pr(>F)"],
summary(anova_math)[[1]][1, "Pr(>F)"]
),
F_ANCOVA = c(
summary(ancova_read)[[1]][2, "F value"],
summary(ancova_write)[[1]][2, "F value"],
summary(ancova_math)[[1]][2, "F value"]
),
p_ANCOVA = c(
summary(ancova_read)[[1]][2, "Pr(>F)"],
summary(ancova_write)[[1]][2, "Pr(>F)"],
summary(ancova_math)[[1]][2, "Pr(>F)"]
)
)
hasil_aa$F_naik <- round(hasil_aa$F_ANCOVA - hasil_aa$F_ANOVA, 2)
print(hasil_aa)
## DV F_ANOVA p_ANOVA F_ANCOVA p_ANCOVA F_naik
## 1 read 3.674017 0.05670663 0.07990165 0.7777275 -3.59
## 2 write 6.714608 0.01027459 1.61966902 0.2046379 -5.09
## 3 math 4.564033 0.03387766 0.31357142 0.5761330 -4.25
Hasil analisis menunjukkan bahwa pada uji ANOVA, variabel write dan math menunjukkan pengaruh signifikan (p < 0.05), sedangkan variabel read tidak signifikan.
Namun, setelah dilakukan analisis ANCOVA dengan memasukkan kovariat science, seluruh variabel dependen menunjukkan nilai p > 0.05, yang berarti tidak terdapat pengaruh signifikan setelah pengaruh kovariat dikontrol.
Nilai F pada ANCOVA juga lebih kecil dibandingkan ANOVA (F_naik negatif), yang menunjukkan bahwa penambahan kovariat science tidak meningkatkan kekuatan model.
library(ggplot2)
plot_aa <- data.frame(
DV = rep(c("read","write","math"), 2),
Metode = rep(c("ANOVA","ANCOVA"), each = 3),
F_val = c(hasil_aa$F_ANOVA, hasil_aa$F_ANCOVA)
)
plot_aa$Metode <- factor(plot_aa$Metode, levels = c("ANOVA","ANCOVA"))
ggplot(plot_aa, aes(x = DV, y = F_val, fill = Metode)) +
geom_bar(stat = "identity", position = "dodge", width = 0.6, color = "white") +
geom_text(aes(label = round(F_val, 1)),
position = position_dodge(width = 0.6),
vjust = -0.5, size = 4, fontface = "bold") +
scale_fill_manual(values = c("#4472C4","#70AD47")) +
labs(
title = "ANOVA vs ANCOVA: F-value efek 'prog'",
subtitle = "ANCOVA mengontrol 'science' --> F-value prog meningkat",
x = "Variabel Dependen", y = "F-value", fill = "Metode"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
legend.position = "top")
manova_aov <- summary.aov(manova_model)
mancova_aov <- summary.aov(mancova_model)
# Ekstraksi manual dengan indeks kolom untuk keamanan
# Kolom 4: F value, Kolom 5: Pr(>F)
hasil_mv <- data.frame(
DV = c("read", "write", "math"),
F_MANOVA = c(
manova_aov[[1]][1, 4],
manova_aov[[2]][1, 4],
manova_aov[[3]][1, 4]
),
p_MANOVA = c(
manova_aov[[1]][1, 5],
manova_aov[[2]][1, 5],
manova_aov[[3]][1, 5]
),
F_MANCOVA = c(
mancova_aov[[1]][2, 4],
mancova_aov[[2]][2, 4],
mancova_aov[[3]][2, 4]
),
p_MANCOVA = c(
mancova_aov[[1]][2, 5],
mancova_aov[[2]][2, 5],
mancova_aov[[3]][2, 5]
)
)
hasil_mv$F_naik <- round(hasil_mv$F_MANCOVA - hasil_mv$F_MANOVA, 2)
print(hasil_mv)
## DV F_MANOVA p_MANOVA F_MANCOVA p_MANCOVA F_naik
## 1 read 3.674017 0.05670663 123.93090 1.204359e-22 120.26
## 2 write 6.714608 0.01027459 87.77026 1.746477e-17 81.06
## 3 math 4.564033 0.03387766 123.28500 1.471013e-22 118.72
Hasil analisis menunjukkan bahwa pada uji MANOVA, variabel write dan math menunjukkan pengaruh signifikan (p < 0.05), sedangkan variabel read tidak signifikan (p > 0.05).
Namun setelah dilakukan analisis MANCOVA dengan memasukkan kovariat science, seluruh variabel dependen menunjukkan hasil yang sangat signifikan (p < 0.05).
Nilai F pada MANCOVA meningkat secara signifikan dibandingkan MANOVA, yang menunjukkan bahwa kovariat science memberikan kontribusi besar dalam meningkatkan kekuatan model.
plot_mv <- data.frame(
DV = rep(c("read","write","math"), 2),
Metode = rep(c("MANOVA","MANCOVA"), each = 3),
F_val = c(hasil_mv$F_MANOVA, hasil_mv$F_MANCOVA)
)
plot_mv$Metode <- factor(plot_mv$Metode, levels = c("MANOVA","MANCOVA"))
ggplot(plot_mv, aes(x = DV, y = F_val, fill = Metode)) +
geom_bar(stat = "identity", position = "dodge", width = 0.6, color = "white") +
geom_text(aes(label = round(F_val, 1)),
position = position_dodge(width = 0.6),
vjust = -0.5, size = 4, fontface = "bold") +
scale_fill_manual(values = c("#ED7D31","#C00000")) +
labs(
title = "MANOVA vs MANCOVA: F-value follow-up univariat efek 'prog'",
subtitle = "MANCOVA mengontrol 'science' → F-value prog meningkat",
x = "Variabel Dependen", y = "F-value", fill = "Metode"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
legend.position = "top")