Почему `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

От RcppArmadillo Библиотека с измененным fLmTwoCasts Я также получаю такое же поведение:

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 может делать что-то вроде lm.fit быстрее, чем функции в пакете stats, заключается в том, что Armadillo использует версию QR-разложения Lapack, а пакет stats использует модифицированную версию Linpack. Следовательно, Armadillo использует код BLAS уровня 3, тогда как пакет статистики использует BLAS уровня 1. Однако Armadillo либо потерпит неудачу, либо, что еще хуже, выдаст совершенно неправильные ответы на матрицах моделей с недостаточным рангом, тогда как функции из пакета stats будут обрабатывать их должным образом из-за модифицированного кода Linpack.

Глядя на код здесь, кишки кода

 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