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

为什么fastLm()当我用一项观察进行回归时返回结果吗？

下面为什么不lm() and 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

来自RcppArmadillo修改过的库电影双播 https://github.com/RcppCore/RcppArmadillo/blob/master/inst/examples/fastLm.r我也得到同样的行为：

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

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

From ?fastLm:

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

看代码here https://github.com/RcppCore/Rcpp/blob/82d8a3a922fc158f2403b233b2ac6861c538d514/inst/examples/Attributes/Depends.cpp，代码的核心是

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

这会进行 QR 分解。我们可以匹配lmFast结果与qr()来自基本 R （这里我不只使用基本 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