Untitled

# 第1组 n1=10
x1_gr1 <- c(5.0, 6.0, 6.8, 7.3, 4.5, 6.2, 5.3, 6.0, 6.9, 6.6)
x2_gr1 <- c(3.0, 4.6, 6.0, 6.5, 3.2, 5.6, 5.5, 3.5, 5.2, 7.3)
# 第2组 n2=10
x1_gr2 <- c(4.6, 3.8, 5.3, 6.8, 4.9, 6.2, 7.1, 4.0, 5.0, 5.8)
x2_gr2 <- c(4.9, 5.4, 4.7, 8.0, 5.9, 6.1, 6.6, 4.1, 7.0, 7.8)

dat <- data.frame(
  group = factor(c(rep("G1",10), rep("G2",10))),
  X1 = c(x1_gr1, x1_gr2),
  X2 = c(x2_gr1, x2_gr2)
)

## 2. 自编MANOVA函数，带if‑else分支逻辑
my_manova <- function(dat, group_col, var_cols, alpha = 0.05){
  X <- as.matrix(dat[, var_cols])
  gvec <- dat[[group_col]]
  gvec <- as.factor(gvec)
  
  N <- nrow(X)
  p <- ncol(X)
  g <- nlevels(gvec)
  
  m <- g - 1        # 假设自由度 m = g‑1
  df_error <- N - g # 误差自由度
  
  # 计算组内E、组间H SSCP矩阵
  mean_total <- colMeans(X)
  E <- matrix(0, nrow = p, ncol = p)
  H <- matrix(0, nrow = p, ncol = p)
  
  for(lev in levels(gvec)){
    subX <- X[gvec == lev, , drop = FALSE]
    nk <- nrow(subX)
    mean_k <- colMeans(subX)
    E <- E + (nk - 1)*cov(subX)
    d_m <- mean_k - mean_total
    H <- H + nk * (d_m %o% d_m)
  }
  T_mat <- E + H
  lambda <- det(E) / det(T_mat)

  # Bartlett修正卡方
  r <- df_error - 0.5*(p - m + 1)
  chi_stat <- - r * log(lambda)
  df_chi <- p * m
  p_chi <- pchisq(chi_stat, df = df_chi, lower.tail = FALSE)
  
  # Box二阶修正
  s_coef <- (p*m*(p^2 + m - 5)) / 16
  cdf1 <- pchisq(chi_stat, df = df_chi)
  cdf2 <- pchisq(chi_stat, df = df_chi + 4)
  cdf_box <- cdf1 + (s_coef / r^2) * (cdf2 - cdf1)
  p_box <- max(0, min(1, 1 - cdf_box))
  
  # if‑else核心分支 
  if(m == 1){
    # m=1：两组，执行精确F检验
    df1_F <- p
    df2_F <- df_error - p + 1
    F_stat <- (df2_F / df1_F) * ((1 - lambda)/lambda)
    p_F <- pf(F_stat, df1 = df1_F, df2 = df2_F, lower.tail = FALSE)
    
    res <- data.frame(
      检验方式 = c("精确F检验"),
      统计量 = c(F_stat),
      df1 = c(df1_F),
      df2 = c(df2_F),
      p值 = c(p_F),
      结论 = ifelse(p_F < alpha, "拒绝H0", "不拒绝H0")
    )
  }else{
    # m>1多组，输出Bartlett渐近卡方、Box修正卡方
    res <- data.frame(
      检验方式 = c("Bartlett渐近卡方","Box高阶修正卡方"),
      统计量 = c(chi_stat, chi_stat),
      df1 = c(df_chi, df_chi),
      df2 = c(NA, NA),
      p值 = c(p_chi, p_box),
      结论 = ifelse(c(p_chi, p_box) < alpha, "拒绝H0", "不拒绝H0")
    )
  }
  return(res)
}

## 3. 运行自编函数
out_self <- my_manova(dat, group_col = "group", var_cols = c("X1","X2"))
cat("=====自编MANOVA输出结果=====\n")

## =====自编MANOVA输出结果=====

print(out_self)

##    检验方式   统计量 df1 df2        p值   结论
## 1 精确F检验 8.163565   2  17 0.00327366 拒绝H0

## 4. R内置manova做对照
cat("\n=====R内置man Wilks检验结果=====\n")

## 
## =====R内置man Wilks检验结果=====

fit <- manova(cbind(X1,X2) ~ group, data = dat)
print(summary(fit, test = "Wilks"))

##           Df   Wilks approx F num Df den Df   Pr(>F)   
## group      1 0.51009   8.1636      2     17 0.003274 **
## Residuals 18                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 5. 顺带一元ANOVA
cat("\n====一元方差分析结果====\n")

## 
## ====一元方差分析结果====

print(summary(aov(X1 ~ group, data = dat)))

##             Df Sum Sq Mean Sq F value Pr(>F)
## group        1  2.521   2.521   2.459  0.134
## Residuals   18 18.449   1.025

print(summary(aov(X2 ~ group, data = dat)))

##             Df Sum Sq Mean Sq F value Pr(>F)
## group        1   5.10   5.100   2.678  0.119
## Residuals   18  34.29   1.905