Die Fastlm-Ergebnisse von RcppArmadillo weichen von denen von R ab. Was habe ich falsch gemacht?

Question

Hier ist die CPP-Datei, die ich als Quelle verwendet habe sourceCpp:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;

// [[Rcpp::export]]
List mylm(NumericVector yr, NumericMatrix Xr) {

    int n = Xr.nrow(), k = Xr.ncol();

    arma::mat X(Xr.begin(), n, k, false);       // reuses memory and avoids extra copy
    arma::colvec y(yr.begin(), yr.size(), false);

    arma::colvec coef = arma::solve(X, y);      // fit model y ~ X
    arma::colvec resid = y - X*coef;            // residuals

    double sig2 = arma::as_scalar( arma::trans(resid)*resid/(n-k) );
                                                // std.error of estimate
    arma::colvec stderrest = arma::sqrt( sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );

    return Rcpp::List::create(
        Rcpp::Named("coefficients") = coef,
        Rcpp::Named("stderr")       = stderrest
    ) ;

}

(Der Code stammt aus der Antwort von @Romain François auf diese Frage: Verwendung von „sourceCpp“ zum Kompilieren von „fastLm“.)

Dann in R:

sourceCpp('b.cpp')
set.seed(1)
x = matrix(rnorm(100), 25, 4)
y = rnorm(25)
mylm(y, x)
summary(lm(y~x))

mylm(y, x)
# $coefficients
#           [,1]
# [1,]  0.068978
# [2,]  0.117632
# [3,] -0.029917
# [4,] -0.168648
# 
# $stderr
#         [,1]
# [1,] 0.17970
# [2,] 0.24312
# [3,] 0.15138
# [4,] 0.21266

summary(lm(y~x))

# Call:
# lm(formula = y ~ x)
# 
# Residuals:
#    Min     1Q Median     3Q    Max 
# -0.869 -0.487 -0.271  0.410  1.831 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)   0.1122     0.1718    0.65     0.52
# x1            0.0499     0.1845    0.27     0.79
# x2            0.1076     0.2470    0.44     0.67
# x3           -0.0435     0.1549   -0.28     0.78
# x4           -0.1750     0.2158   -0.81     0.43
# 
# Residual standard error: 0.835 on 20 degrees of freedom
# Multiple R-squared:  0.054,     Adjusted R-squared:  -0.135 
# F-statistic: 0.285 on 4 and 20 DF,  p-value: 0.884

Answer 1

Standardmäßig passt Rs lm ein Modell mit einem Achsenabschnitt an, auch wenn die von Ihnen übergebene Entwurfsmatrix keine anfängliche Einsenspalte enthält.Sie werden also sehen, dass Folgendes identisch ist:

lm(y ~ x - 1)
mylm(y, x)

Wenn Sie das „normale“ Modell wünschen, müssen Sie Ihre Designmatrix so ändern, dass sie in der ersten Spalte nur Einsen enthält:

lm(y ~ x)
mylm(y, cbind(1, x))

wird identische Ergebnisse liefern.