从 .lm.fit() 计算 p 值的快速方法

Question

我正在运行模拟和拟合线性模型 .lm.fit(). Although extremely fast, this function does not provide the predictors' p-values. Is there a fast way to compute them (maybe from the values returned by .lm.fit())? I am aware of this 方法来计算近似的 p 值，但我需要精确的 p 值。

更新:
Dirk Eddelbuettel 提供了拟合 lm 的最快方法，Ben Bolker 提供了计算 p 值的方法，通过结合我们得到的两个答案：

set.seed(101)
X <- cbind(1,matrix(1:10))
y <- rnorm(10)

mdl <- RcppArmadillo::fastLmPure(X, y)

pval <- 2*pt(abs(mdl$coefficients/mdl$stderr), mdl$df.residual, lower.tail=FALSE)

Answer 1

对于这个问题（获取标准错误，因此 p-values）我在包 RcppArmadillo、RcppEigen 和 RcppGSL 中编写了函数的三个不同版本 fastLm()。其中一部分当然也只是为了说明。但你可以从那里开始。确保您使用 fastLmPure() 变体获取向量和矩阵，并且 而不是 公式界面——所有时间都花在解析公式上。

在这里，只是为了好玩，是 RcppArmadillo 变体：

#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::export]]
List fastLm_impl(const arma::mat& X, const arma::colvec& y) {
    int n = X.n_rows, k = X.n_cols;

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

    // std.errors of coefficients
    double s2 = std::inner_product(res.begin(), res.end(), 
                                   res.begin(), 0.0)/(n - k);

    arma::colvec std_err = 
         arma::sqrt(s2 *
                    arma::diagvec(arma::pinv(arma::trans(X)*X)));  

    return List::create(Named("coefficients") = coef,
                        Named("stderr")       = std_err,
                        Named("df.residual")  = n - k);
}

Answer 2

Dirk 的回答会更快，但如果方便的话，这里是纯 R 中的实现（从 summary.lm 中提取您需要的位，并假设 non-full-rank 模型矩阵等没有问题.)

示例：

set.seed(101)
X <- cbind(1,matrix(1:10))
y <- rnorm(10)
m <- .lm.fit(X,y)

p-value 计算：

rss <- sum(m$residuals^2)
rdf <- length(y) - ncol(X)
resvar <- rss/rdf
R <- chol2inv(m$qr)
se <- sqrt(diag(R) * resvar)
2*pt(abs(m$coef/se),rdf,lower.tail=FALSE)

比较：

coef(summary(lm(y~X-1)))[,"Pr(>|t|)"]

从 .lm.fit() 计算 p 值的快速方法

fast method to calculate p-values from .lm.fit()

r

lm