R: 具有多个因子值的简单方差分析

Question

我正在尝试一个具有许多因子值7000+的简单方差分析。在尝试使用aov命令时，我的计算机陷入了极大的困境（永远无法完成）。我不明白为什么会出现这种情况，因为手动编程我自己的anova似乎运行得更快。然而，挑战在于，我的方差分析在小样本量上与基本方差分析不同，我不确定确切的原因（见下文）。

我的问题是，是否有人可以选择更快的anova。此外，如果有人受到了同样的启发，我会附上我的"anova"。如果有人受到这样的启发，我将非常感谢任何反馈，看看我在写这篇文章时是否犯了错误。我的anova似乎比基础跑得快10-100倍。但话说回来，也许我犯了一些错误。

anova <- function(X, ID, na.rm=TRUE) {
  N <- sum(!is.na(X))
  K <- length(unique(ID))
  Xbari <- ave(X, ID, FUN=function(x) mean(x, na.rm=na.rm))
  Xbar  <- mean(X, na.rm=na.rm)
  ESS <- sum((Xbar-Xbari)^2, na.rm=na.rm)
  USS <- sum((X-Xbari)^2,   na.rm=na.rm)
  TDF <- (K-1)
  DFi <- (N-K)
  fstat = (ESS/TDF)/(USS/DFi)
  data.frame(p=pf(fstat, TDF, DFi),ESS=ESS,
             USS=USS,TDF=TDF, DFi=DFi, fstat=fstat)
}

Answer 1

aov基本上是lm的包装器。lm所做的第一步是创建设计矩阵，即将因子变量转换为伪变量，即将长度为n的因子转换为尺寸为n*7000的矩阵。超过7000个变量的最小二乘拟合可能会很慢。

如果anova函数中的p值计算不正确。试试这个：

anova <- function(X, ID, na.rm=TRUE) {
  stopifnot(length(X) == length(ID))
  nas <- is.na(X) | is.na(ID)
  stopifnot(!any(nas) & na.rm)
  X <- X[!nas]
  ID <- ID[!nas]
  N <- length(X)
  K <- length(unique(ID))
  Xbari <- ave(X, ID, FUN=mean)
  Xbar  <- mean(X)
  ESS <- sum((Xbar-Xbari)^2)
  USS <- sum((X-Xbari)^2)
  TDF <- (K-1)
  DFi <- (N-K)
  fstat = (ESS/TDF)/(USS/DFi)
  data.frame(p=format.pval(pf(fstat, TDF, DFi, lower.tail = FALSE)),ESS=ESS,
             USS=USS,TDF=TDF, DFi=DFi, fstat=fstat)
}
anova(npk$yield, npk$block)
#         p     ESS    USS TDF DFi    fstat
#1 0.086072 343.295 533.07   5  18 2.318386
summary(aov(yield ~ block, data = npk))
#            Df Sum Sq Mean Sq F value Pr(>F)  
#block        5  343.3   68.66   2.318 0.0861 .
#Residuals   18  533.1   29.61                 
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

通过使用data.table或dplyr来计算组平均值，可以更快地实现这一点。