Мне нужно вычислить сумму квадратов матрицы перекрестного произведения (на самом деле след этой матрицы) в многомерной линейной модели с 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, я попытался повторно реализовать этот расчет. Хотя моя реализация, напрямую переводящая линейную алгебру для s = 1 (ниже), работает быстрее для небольших размеров выборки (n),
# 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 все еще быстрее вычисляет все матрицы p = 30, чем я, вычисляющий только одну (s = 1)!! например 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 код (т. е. функции 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)
}
Тогда это просто вопрос выбора подмножества терминов.
