manova-2 way manova

[이한숙]

# ✅ 1. 필수 패키지 설치 및 로드
install.packages(c("car", "broom", "writexl"), dependencies = TRUE)
library(car)
library(broom)
library(writexl)
library(dplyr)

# ✅ 2. 데이터 불러오기 (CSV 경로 사용자에 맞게 수정)
data1 <- read.csv("2025.04.18-1coding.csv")

# ✅ 3. 범주형 변수 변환
data1$group <- factor(data1$group, levels = c(1, 2), labels = c("정상군", "인지손상군"))
data1$condition <- factor(data1$condition, levels = c(1, 2, 3), labels = c("평지", "물컵", "숫자세기"))

# ✅ 4. 종속변수 목록 (velocity 포함)
gait_vars <- c(
  "loadingresponsetime", "loadingresponsegc", "midstancetime", "midstancegc",
  "terminalstancetime", "terminalstancegc", "preswingtime", "preswinggc",
  "swingphase", "swingsi", "cadence",  "stancephasetime", "stancesi", "cofsi", "velocity"
)

# ✅ 5. 다변량 모델 정의 및 수행
formula <- as.formula(paste0(
  "cbind(", paste(gait_vars, collapse = ", "), ") ~ group * condition"
))

manova_model <- manova(formula, data = data1)

# ✅ 6. 다변량 결과 (Pillai's Trace 포함) + Partial Eta² 계산
manova_summary <- summary(manova_model, test = "Pillai")
pillai_df <- as.data.frame(manova_summary$stats)
pillai_df$Effect <- rownames(pillai_df)
rownames(pillai_df) <- NULL
pillai_df <- pillai_df[, c("Effect", "Pillai", "approx F", "num Df", "den Df", "Pr(>F)")]
colnames(pillai_df) <- c("Effect", "Pillai", "F_approx", "num_Df", "den_Df", "p_value")

# Partial Eta² 계산
pillai_df$Partial_Eta2 <- 1 - (1 - pillai_df$Pillai)^(1 / pmin(pillai_df$num_Df, pillai_df$den_Df))

# ✅ 7. 단변량 ANOVA 요약 + Partial Eta² 계산
aov_summary <- summary.aov(manova_model)
result_list <- list()

for (i in seq_along(aov_summary)) {
  df <- as.data.frame(aov_summary[[i]])
  variable_name <- names(aov_summary)[i]
  df$Effect <- rownames(df)
  df$Variable <- variable_name

  if ("Residuals" %in% df$Effect) {
    residual_ss <- df[df$Effect == "Residuals", "Sum Sq"]
    df$Partial_Eta2 <- NA
    for (j in 1:nrow(df)) {
      if (df$Effect[j] != "Residuals") {
        df$Partial_Eta2[j] <- df$`Sum Sq`[j] / (df$`Sum Sq`[j] + residual_ss)
      }
    }
  } else {
    df$Partial_Eta2 <- NA
  }

  result_list[[i]] <- df
}

anova_df <- bind_rows(result_list) %>%
  select(Variable, Effect, Df, `Sum Sq`, `Mean Sq`, `F value`, `Pr(>F)`, Partial_Eta2)

# ✅ 8. 결과 저장
write_xlsx(
  list(
    "MANOVA_Pillai" = pillai_df,
    "Univariate_ANOVA" = anova_df
  ),
  path = "TwoWay_MANOVA_with_Pillai_Eta2.xlsx"
)

cat("✅ 분석 완료: 결과가 'TwoWay_MANOVA_with_Pillai_Eta2.xlsx'에 저장되었습니다.\n")