R car :: Anova的更快替代方法，用于预测子集的平方叉积矩阵求和

如何解决R car :: Anova的更快替代方法，用于预测子集的平方叉积矩阵求和

我需要使用Y（n x q）和X（n x p）在多元线性模型中计算叉积平方和（实际上是该矩阵的迹线）。为此的标准R代码是：

require(MASS)
require(car)

# Example data 
q <- 10
n  <- 1000
p <- 10
Y <- mvrnorm(n,mu = rep(0,q),Sigma = diag(q))
X <- as.data.frame(mvrnorm(n,mu = rnorm(p),Sigma = diag(p)))

# Fit lm
fit <- lm( Y ~ .,data = X )

# Type I sums of squares
summary(manova(fit))$SS    

# Type III sums of squares
type = 3 # could be also 2 (II)
car::Anova(fit,type = type)$SSP

必须执行数千次，不幸的是，当预测变量的数量相对较大时，它会变慢。通常，我只对s个预测变量的子集感兴趣，因此我尝试重新实现此计算。尽管对于小样本量（n），我的实现直接转换s = 1（如下）的线性代数的速度更快，

# Hat matrix (X here stands for the actual design matrix)
H <- tcrossprod(tcrossprod(X,solve(crossprod(X))),X)

# Remove predictor of interest (e.g. 2)
X.r <- X[,-2]  
H1 <- tcrossprod(tcrossprod(X.r,solve(crossprod(X.r))),X.r) 

# Compute e.g. type III sum of squares
SS <- crossprod(Y,H - H1) %*% Y

car对于较大的n而言仍然更快：

我已经尝试过Rcpp实现，因为R中的这些矩阵产品已经使用了非常有效的代码。

关于如何更快地执行此操作的任何提示？

更新

阅读答案后，我尝试了此post中提出的解决方案，该解决方案依赖于QR / SVD / Cholesky因式分解来进行帽子矩阵计算。但是，似乎car::Anova仍然比我只计算一个（s = 1）更快，计算所有p = 30矩阵！例如n = 5000，q = 10：

Unit: milliseconds
 expr       min        lq      mean    median        uq       max neval
   ME 1137.5692 1202.9888 1257.8979 1251.6834 1318.9282 1398.9343    10
   QR 1005.9082 1031.9911 1084.5594 1037.5659 1095.7449 1364.9508    10
  SVD 1026.8815 1065.4629 1152.6631 1087.9585 1241.4977 1446.8318    10
 Chol  969.9089 1056.3093 1115.9608 1102.1169 1210.7782 1267.1274    10
  CAR  205.1665  211.8523  218.6195  214.6761  222.0973  242.4617    10

更新2

目前最好的解决方案是遍历car::Anova code（即功能car:::Anova.III.mlm，随后是car:::linearHypothesis.mlm），然后重新实现它们以解决预测变量，而不是所有预测变量。

car的相关代码如下（我跳过了检查，并简化了一点）：

B <- coef(fit)                    # Model coefficients
M <- model.matrix(fit)            # Model matrix M
V <- solve(crossprod(M))          # M'M
p <- ncol(M)                      # Number of predictors in M
I.p <- diag(p)                    # Identity (p x p)
terms <- labels(terms(fit))       # terms (add intercept)       
terms <- c("(Intercept)",terms)   
n.terms <- length(terms)
assign <- fit$assign              # assignation terms <-> p variables
  
SSP <- as.list(rep(0,n.terms))   # Initialize empty list for sums of squares cross-product matrices
names(SSP) <- terms
  
for (term in 1:n.terms){
    subs <- which(assign == term - 1)
    L <- I.p[subs,drop = FALSE]
    SSP[[term]] <- t(L %*% B) %*% solve(L %*% V %*% t(L)) %*% (L %*% B)
}

然后，只需选择条款的子集即可。

解决方法

此行及其下面类似的Rows Removed by Filter = 95 Shared Hit Blocks = 34 Plan Rows = 1 Actual Rows = 0 可能会得到改进：

H1

通常的想法是您应该很少使用H <- tcrossprod(tcrossprod(X,solve(crossprod(X))),X) ，因为它与solve(Y) %*% Z相同，但是速度较慢。我还没有完全扩展您的solve(Y,Z)调用，以查看tcrossprod和H的最佳表达式表示方式。

您也可以查看问题https://stats.stackexchange.com/questions/139969/speeding-up-hat-matrices-like-xxx-1x-projection-matrices-and-other-as，以获取通过QR分解进行此操作的说明。