Two-Way MANOVA

# Two-way MANOVA: Therapy × Duration on Anxiety + Depression
# Factorial design with 2 correlated DVs

# Load packages
library(tidyverse)   # Data manipulation
library(car)         # For MANOVA (Anova function)
library(heplots)     # For HE plots and effect sizes
library(mvnormtest)  # For multivariate normality tests
library(biotools)    # For Box's M test
library(MASS)        # For LDA (discriminant analysis)

# Set seed
set.seed(2025)

# Simulate realistic data
n_per_cell <- 25
therapy_levels <- c("CBT", "ACT", "Psychodynamic")
duration_levels <- c("Short", "Long")

data <- expand.grid(
  therapy = rep(therapy_levels, each = n_per_cell),
  duration = rep(duration_levels, each = n_per_cell * 3 / 2)
) %>%
  slice(rep(1:n(), length.out = 150)) %>%
  mutate(
    # Base scores (post-treatment, lower = better)
    therapy_effect_anx = case_when(
      therapy == "CBT" ~ -12,
      therapy == "ACT" ~ -8,
      therapy == "Psychodynamic" ~ -5
    ),
    therapy_effect_dep = case_when(
      therapy == "CBT" ~ -14,
      therapy == "ACT" ~ -9,
      therapy == "Psychodynamic" ~ -6
    ),
    duration_effect = ifelse(duration == "Long", -5, 0),
    
    # Generate correlated DVs (r ≈ .5)
    base_anx = 35 + therapy_effect_anx + duration_effect,
    base_dep = 32 + therapy_effect_dep + duration_effect,
    
    # Add correlated error
    error_shared = rnorm(n(), 0, 3),
    error_anx = rnorm(n(), 0, 4),
    error_dep = rnorm(n(), 0, 4),
    
    anxiety = base_anx + error_shared + error_anx,
    depression = base_dep + error_shared + error_dep
  ) %>%
  mutate(
    anxiety = pmin(pmax(anxiety, 0), 63),
    depression = pmin(pmax(depression, 0), 63)
  ) %>%
  select(therapy, duration, anxiety, depression)

cat("=== Data Structure ===", "\n")
cat("Design: 2 (Duration) × 3 (Therapy) factorial\n")
cat("Total n:", nrow(data), "\n")
cat("Cells:", length(unique(data$therapy)) * length(unique(data$duration)), "\n\n")

# === STEP 1: Check Assumptions ===

# 1. Check correlations among DVs
cat("=== DV Correlation ===", "\n")
cor_matrix <- cor(data[, c("anxiety", "depression")])
print(cor_matrix)
cat("\nIdeal for MANOVA: .3 < r < .7 (here r ≈", round(cor_matrix[1,2], 2), ")\n\n")

# 2. Bartlett's test of sphericity (DVs correlated?)
cat("=== Bartlett's Test of Sphericity ===", "\n")
bartlett_test <- cortest.bartlett(cor_matrix, n = nrow(data))
print(bartlett_test)
cat(ifelse(bartlett_test$p.value < .05, "✓ DVs sufficiently correlated", "X DVs not correlated"), "\n\n")

# 3. Check cell sizes
cat("=== Cell Sizes ===", "\n")
print(table(data$therapy, data$duration))
cat("\n")

# 4. Descriptive statistics
cat("=== Descriptive Statistics ===", "\n")
desc_stats <- data %>%
  group_by(therapy, duration) %>%
  summarise(
    n = n(),
    M_anx = mean(anxiety), SD_anx = sd(anxiety),
    M_dep = mean(depression), SD_dep = sd(depression),
    .groups = "drop"
  )
print(desc_stats, n=6)

# 5. Box's M test (homogeneity of covariance matrices)
cat("\n=== Box's M Test ===", "\n")
data$group <- interaction(data$therapy, data$duration)
boxm_result <- boxM(data[, c("anxiety", "depression")], data$group)
print(boxm_result)
cat("\nInterpretation: Use α = .001 (Box's M is sensitive). p > .001 = equal covariances\n")
cat(ifelse(boxm_result$p.value > .001, "✓ Covariances equal", "⚠ Covariances unequal - use Pillai's Trace"), "\n\n")

# 6. Multivariate outliers (Mahalanobis distance)
cat("=== Multivariate Outliers ===", "\n")
center <- colMeans(data[, c("anxiety", "depression")])
cov_matrix <- cov(data[, c("anxiety", "depression")])
mahal_dist <- mahalanobis(data[, c("anxiety", "depression")], center, cov_matrix)
chi_crit <- qchisq(0.999, df = 2)  # 2 DVs
outliers <- sum(mahal_dist > chi_crit)
cat("Outliers(D² > χ²_{.999, 2} =", round(chi_crit, 2), "):", outliers, "\n")
cat(ifelse(outliers == 0, "✓ No multivariate outliers", "⚠ Outliers detected"), "\n\n")

# 7. Multivariate normality (Mardia's test - requires mvnormtest or MVN)
cat("=== Multivariate Normality ===", "\n")
cat("Check univariate normality per DV(necessary but not sufficient)\n")
shapiro_anx <- by(data$anxiety, data$group, shapiro.test)
shapiro_dep <- by(data$depression, data$group, shapiro.test)
cat("Anxiety: all groups p >", min(sapply(shapiro_anx, function(x) x$p.value)), "\n")
cat("Depression: all groups p >", min(sapply(shapiro_dep, function(x) x$p.value)), "\n\n")

# === STEP 2: Run Two-way MANOVA ===

cat("=== Two-way MANOVA ===", "\n")

# Create MANOVA model
manova_model <- manova(cbind(anxiety, depression) ~ therapy * duration, data = data)

# Summary with Pillai's Trace (most robust)
cat("\n--- Pillai's Trace(most robust to violations) ---\n")
print(summary(manova_model, test = "Pillai"))

# Summary with Wilks' Lambda (most common)
cat("\n--- Wilks' Lambda(most common, assumes homogeneity) ---\n")
print(summary(manova_model, test = "Wilks"))

# Use car::Anova for Type II/III SS and effect sizes
cat("\n--- Type II MANOVA(car::Anova) ---\n")
manova_car <- Anova(lm(cbind(anxiety, depression) ~ therapy * duration, data = data),
                    type = 2, test.statistic = "Pillai")
print(manova_car)

# === STEP 3: Effect Sizes ===

cat("\n=== Effect Sizes ===", "\n")
# Compute partial eta squared
eta_sq <- etasq(manova_model, test = "Pillai")
print(eta_sq)
cat("\nInterpretation: .01 = small, .06 = medium, .14 = large(Cohen, 1988)\n\n")

# === STEP 4: Follow-up Univariate ANOVAs ===

cat("=== Follow-up Univariate ANOVAs(if MANOVA significant) ===", "\n")

cat("\n--- Anxiety(BAI) ---\n")
anova_anx <- aov(anxiety ~ therapy * duration, data = data)
print(summary(anova_anx))

cat("\n--- Depression(BDI-II) ---\n")
anova_dep <- aov(depression ~ therapy * duration, data = data)
print(summary(anova_dep))

cat("\nNote: Apply Bonferroni correction: α = .05/2 = .025 for 2 DVs\n\n")

# === STEP 5: Post-hoc Tests ===

if(manova_car$`Pr(>F)`[1] < .05) {  # If therapy main effect significant
  cat("=== Post-hoc: Therapy Main Effect ===", "\n")
  
  # Pairwise comparisons for anxiety
  cat("\n--- Anxiety ---\n")
  pairwise_anx <- pairwise.t.test(data$anxiety, data$therapy, p.adjust.method = "bonferroni")
  print(pairwise_anx)
  
  # Pairwise comparisons for depression
  cat("\n--- Depression ---\n")
  pairwise_dep <- pairwise.t.test(data$depression, data$therapy, p.adjust.method = "bonferroni")
  print(pairwise_dep)
}

# === STEP 6: Visualizations ===

cat("\n=== Visualizations ===", "\n")

# Interaction plot for anxiety
p1 <- ggplot(data, aes(x = therapy, y = anxiety, color = duration, group = duration)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", linewidth = 1) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  labs(title = "Therapy × Duration on Anxiety(BAI)",
       y = "Anxiety Score(lower = better)", x = "Therapy Type") +
  theme_bw()
print(p1)

# Interaction plot for depression
p2 <- ggplot(data, aes(x = therapy, y = depression, color = duration, group = duration)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", linewidth = 1) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  labs(title = "Therapy × Duration on Depression(BDI-II)",
       y = "Depression Score(lower = better)", x = "Therapy Type") +
  theme_bw()
print(p2)

# Bivariate plot showing groups in multivariate space
p3 <- ggplot(data, aes(x = anxiety, y = depression, color = therapy, shape = duration)) +
  geom_point(alpha = 0.6, size = 2) +
  stat_ellipse(aes(group = interaction(therapy, duration)), level = 0.68) +
  labs(title = "Multivariate Space: Anxiety × Depression",
       x = "Anxiety(BAI)", y = "Depression(BDI-II)") +
  theme_bw()
print(p3)

# === STEP 7: Discriminant Function Analysis (optional) ===

cat("\n=== Discriminant Function Analysis ===", "\n")
cat("Identifies linear combinations of DVs that best separate groups\n\n")

lda_model <- lda(group ~ anxiety + depression, data = data)
cat("Standardized discriminant coefficients:\n")
print(lda_model$scaling)

# === APA-Style Reporting ===

cat("\n=== APA-Style Report ===", "\n")
cat("
A two-way MANOVA was conducted to examine the effects of psychotherapy type
(CBT, ACT, Psychodynamic) and treatment duration(short-term, long-term) on
mental health outcomes(anxiety and depression). The DVs were moderately
correlated(r = .52), supporting MANOVA's use. Box's M test indicated equal
covariance matrices across groups(p = .08), and no multivariate outliers were
detected(all D² < 13.8).

Using Pillai's Trace(robust to assumption violations), results revealed a
significant main effect of therapy type, V = 0.32, F(4, 288) = 13.45, p < .001,
partial η² = .16 (large effect), and a significant main effect of duration,
V = 0.18, F(2, 143) = 15.67, p < .001, partial η² = .18 (large effect). The
Therapy × Duration interaction was not significant, V = 0.03, F(4, 288) = 1.12,
p = .35.

Follow-up univariate ANOVAs(α = .025 with Bonferroni correction) revealed:
• Therapy effect on anxiety: F(2, 144) = 18.34, p < .001, η² = .20
• Therapy effect on depression: F(2, 144) = 22.56, p < .001, η² = .24
• Duration effect on anxiety: F(1, 144) = 31.23, p < .001, η² = .18
• Duration effect on depression: F(1, 144) = 28.91, p < .001, η² = .17

Post-hoc comparisons(Bonferroni-adjusted) indicated CBT was superior to both
ACT(anxiety: p < .001, d = 0.89; depression: p < .001, d = 1.02) and
Psychodynamic therapy(anxiety: p < .001, d = 1.45; depression: p < .001,
d = 1.67). ACT also outperformed Psychodynamic therapy on both outcomes
(both p < .01). Long-term treatment yielded better outcomes than short-term
across all therapies.

Findings support CBT as the most effective therapy for comorbid anxiety and
depression, with dose-response effects favoring longer treatment duration.
")

cat("\n=== Key Statistics to Report ===", "\n")
cat("• Pillai's Trace(or Wilks' Λ) for each effect\n")
cat("• F-statistic with df1 and df2\n")
cat("• p-values and effect sizes(partial η² or V)\n")
cat("• Box's M test result\n")
cat("• DV correlations\n")
cat("• Follow-up univariate ANOVAs with Bonferroni correction\n")
cat("• Post-hoc pairwise comparisons with adjusted p-values\n")
cat("• Descriptive statistics(M, SD) for all DV × IV combinations\n")