RcppArmadillo fastlm results differ from R's lm, what have I done wrong?
armadillo, c++, r, rcpp
Solution
By default, R's lm fits a model with an intercept even if the design matrix you pass does not include an initial column of ones. So you'll see the following are identical:
lm(y ~ x - 1)
mylm(y, x)
If you want the 'regular' model, you need to modify your design matrix to have a first column of all 1s:
lm(y ~ x)
mylm(y, cbind(1, x))
will give identical results.
Problem
Here is the cpp file that I sourced with `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 ) ; } ``` (The code is taken from the answer of @Romain François to this question: Using `sourceCpp` to compile `fastLm`) Then 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 ```