Kasus: membandingkan SBP dan LDL antara Perokok vs Bukan Perokok

Saya ingin mengetahui apakah perokok dan bukan perokok berbeda secara bersama-sama pada dua indikator kesehatan kardiovaskular: tekanan darah sistolik (SBP) dan kolesterol LDL (LDL).

set.seed(1024)

# --- Data contoh ---
# Kolom: SBP (mmHg), LDL (mg/dL)
perokok <- data.frame(
  SBP = c(142, 135, 150, 138, 148, 144, 139, 151, 145, 140),
  LDL = c(132, 125, 140, 128, 135, 130, 129, 142, 136, 131)
)

bukan_perokok <- data.frame(
  SBP = c(130, 128, 135, 124, 129, 133, 131, 127, 136, 125),
  LDL = c(115, 118, 120, 112, 117, 119, 116, 114, 121, 110)
)
# Ukuran sampel dan p
n1 <- nrow(perokok)
n2 <- nrow(bukan_perokok)
p  <- ncol(perokok)
# Mean vectors
mean1 <- colMeans(perokok)        # rata-rata perokok
mean2 <- colMeans(bukan_perokok) # rata-rata bukan perokok
diff_mean <- as.matrix(mean1 - mean2)
# Sample covariance matrices
S1 <- cov(perokok)
S2 <- cov(bukan_perokok)
# Pooled covariance matrix
Spooled <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n1 + n2 - 2)
# Hitung Hotelling T^2
T2 <- as.numeric( (n1 * n2) / (n1 + n2) * t(diff_mean) %*% solve(Spooled) %*% diff_mean )
# --- Konversi ke statistik F ---
df1 <- p
df2 <- n1 + n2 - p - 1
Fstat <- ( (n1 + n2 - p - 1) * T2 ) / ( (n1 + n2 - 2) * p )
p_value <- 1 - pf(Fstat, df1, df2)
# Output hasil
cat("n1 (perokok) =", n1, "\n")
## n1 (perokok) = 10
cat("n2 (bukan)   =", n2, "\n\n")
## n2 (bukan)   = 10
cat("Mean Perokok:      ", round(mean1,2), "\n")
## Mean Perokok:       143.2 132.8
cat("Mean Bukan Perokok:", round(mean2,2), "\n\n")
## Mean Bukan Perokok: 129.8 116.2
cat("Pooled covariance matrix:\n")
## Pooled covariance matrix:
print(round(Spooled,2)); cat("\n")
##       SBP   LDL
## SBP 22.40 20.04
## LDL 20.04 20.73
cat("Hotelling T2 statistic =", round(T2,4), "\n")
## Hotelling T2 statistic = 78.0604
cat("F-statistic =", round(Fstat,4), " with df1 =", df1, " df2 =", df2, "\n")
## F-statistic = 36.8619  with df1 = 2  df2 = 17
cat("p-value =", signif(p_value,4), "\n\n")
## p-value = 6.579e-07
alpha <- 0.05
if (p_value < alpha) {
  cat("Keputusan: Tolak H0 (ada perbedaan vektor rata-rata antara perokok dan bukan perokok) \n")
} else {
  cat("Keputusan: Gagal tolak H0 (tidak ada bukti perbedaan vektor rata-rata pada alpha =", alpha, ")\n")
}
## Keputusan: Tolak H0 (ada perbedaan vektor rata-rata antara perokok dan bukan perokok)