当我对一个观察值进行回归时，为什么“fastLm()”会返回结果？

Question

当我对一个观察值运行回归时，为什么 fastLm() 返回结果？

在下面，为什么lm() 和fastLm() 结果不相等？

library(Rcpp)
library(RcppArmadillo)
library(data.table)
set.seed(1)
DT <- data.table(y = rnorm(5), x1 = rnorm(5), x2 = rnorm(5), my.key = 1:5)
#             y         x1         x2 my.key
# 1: -0.6264538 -0.8204684  1.5117812      1
# 2:  0.1836433  0.4874291  0.3898432      2
# 3: -0.8356286  0.7383247 -0.6212406      3
# 4:  1.5952808  0.5757814 -2.2146999      4
# 5:  0.3295078 -0.3053884  1.1249309      5

lm(y ~ 1 + x1 + x2, data = DT[my.key == 1])
# Coefficients:
# (Intercept)           x1           x2  
#     -0.6265           NA           NA

fastLm(X = model.matrix(y ~ 1 + x1 + x2, data = DT[my.key == 1]), y = DT[my.key == 1]$y)
# Coefficients:
# (Intercept)          x1          x2 
#    -0.15825     0.12984    -0.23924 

model.matrix(y ~ 1 + x1 + x2, data = DT[my.key == 1])
#   (Intercept)        x1       x2
#             1 -0.8204684 1.511781
# attr(,"assign")
# [1] 0 1 2

DT[my.key == 1]$y
# [1] -0.6264538

当我使用整个 DT 时，结果是相等的：

 all.equal(fastLm(X = model.matrix(y ~ 1 + x1 + x2, data = DT), y = DT$y)$coef, 
           lm.fit(x = model.matrix(y ~ 1 + x1 + x2, data = DT), y = DT$y)$coef)
# [1] TRUE

从带有修改的fLmTwoCasts 的RcppArmadillo 库中，我也得到了相同的行为：

src <- '
Rcpp::List fLmTwoCastsOnlyCoefficients(Rcpp::NumericMatrix Xr, Rcpp::NumericVector yr) {
    int n = Xr.nrow(), k = Xr.ncol();
    arma::mat X(Xr.begin(), n, k, false);
    arma::colvec y(yr.begin(), yr.size(), false);
    arma::colvec coef = arma::solve(X, y);
    return Rcpp::List::create(Rcpp::Named("coefficients")=trans(coef));
}
'
cppFunction(code=src, depends="RcppArmadillo")

XX <- model.matrix(y ~ 1 + x1 + x2, data = DT[my.key == 1])
YY <- DT[my.key == 1]$y
fLmTwoCastsOnlyCoefficients(XX, YY)$coef
#            [,1]      [,2]       [,3]
# [1,] -0.1582493 0.1298386 -0.2392384

使用整个DT，系数应该是相同的：

lm(y ~ 1 + x1 + x2, data = DT)$coef
# (Intercept)          x1          x2 
#   0.2587933  -0.7709158  -0.6648270

XX.whole <- model.matrix(y ~ 1 + x1 + x2, data = DT)
YY.whole <- DT$y
fLmTwoCastsOnlyCoefficients(XX.whole, YY.whole)
#           [,1]       [,2]      [,3]
# [1,] 0.2587933 -0.7709158 -0.664827

Answer 1

因为fastLm 不担心排名不足；这是您为速度付出的代价的一部分。

来自?fastLm：

... Armadillo 可以比 stats 包中的函数更快地执行类似 lm.fit 的操作的原因是，Armadillo 使用 QR 分解的 Lapack 版本，而 stats 包使用修改后的 Linpack 版本。因此，Armadillo 使用 3 级 BLAS 代码，而 stats 包使用 1 级 BLAS。然而，犰狳会失败，或者更糟的是，在秩不足的模型矩阵上产生完全错误的答案，而由于修改了 Linpack 代码，stats 包中的函数将正确处理它们。

看代码here，代码的胆子是

 arma::colvec coef = arma::solve(X, y);

这会进行 QR 分解。我们可以将lmFast 结果与来自base R 的qr() 匹配（这里我不只使用base R 构造，而是依赖data.table）：

set.seed(1)
dd <- data.frame(y = rnorm(5), 
      x1 = rnorm(5), x2 = rnorm(5), my.key = 1:5)

X <- model.matrix(~1+x1+x2, data=subset(dd,my.key==1))
qr(X,dd$y)
## $qr
##   (Intercept)         x1       x2
## 1           1 -0.8204684 1.511781

你可以看lm.fit()的代码，看看R在拟合线性模型时对秩亏做了什么；它调用的底层 BLAS 算法通过旋转来执行 QR ...

如果你想标记这些情况，我认为Matrix::rankMatrix() 可以解决问题：

library(Matrix)
rankMatrix(X) < ncol(X)  ## TRUE
X1 <- model.matrix(~1+x1+x2, data=dd)
rankMatrix(X1) < ncol(X1)  ## FALSE