一次适合许多公式，比lapply更快的选择？

Question

我有一个适合数据的公式列表，而不是运行一个循环，为了性能，我希望一次执行此操作。估计应该仍然是分开的，我并不是想估计一个SUR或任何东西。 下面的代码做什么，我想

x <- matrix(rnorm(300),ncol=3) 
y <- x %*% c(1,2,3)+rnorm(100) 
formulae <-list(y~x[,1], 
       y~x[,2], 
       y~x[,1] + x[,2]) 
lapply(formulae,lm) 

可惜，这得有点慢作为的formulae长度的增加，有没有办法真正地矢量化呢？

如果有任何帮助，我只关心lm的唯一结果是系数和一些标准错误。

Answer 1

这不是一种简单的矢量化方法，但MuMIn软件包中的pdredge函数为您提供了一种非常简单的并行化方法（假定您的计算机上有多个核心，或者可以设置本地群集在由parallel包所支持的方式之一...

library(parallel) 
clust <- makeCluster(2,"PSOCK") 
library(MuMIn) 

构造数据：

set.seed(101) 
x <- matrix(rnorm(300),ncol=3) 
y <- x %*% c(1,2,3)+rnorm(100) 

这将是更容易被命名为数据帧要做到这一点，而不是一个anonymou S矩阵：

df <- setNames(data.frame(y,x),c("y",paste0("x",1:3))) 

集群节点都需要访问数据集：

clusterExport(clust,"df") 

拟合的是整个模型（你可以使用y~.适合所有变量）

full <- lm(y~x1+x2,data=df,na.action=na.fail) 

现在适合所有的子模型（见?MuMIn::dredge有更多的选项来控制安装哪些子模型）

p <- pdredge(full,cluster=clust) 
coef(p) 
## (Intercept)  x1  x2 
## 3 -0.003805107 0.7488708 2.590204 
## 2 -0.028502039  NA 2.665305 
## 1 -0.101434662 1.0490816  NA 
## 0 -0.140451160  NA  NA 

Answer 2

正如我在我的评论中所说的，你真正需要的是一个更有效但稳定的配件，而不是lm()。在这里，我会为您提供一个经过良好测试的自己写的，名为lm.chol()。它需要一个formula和data，返回：

• 系数汇总表，因为你在summary(lm(...))$coef通常看到;
• Pearson残差标准误差估计值，从summary(lm(...))$sigma得到;
• 调整-R.squared，你从summary(lm(...))$adj.r.squared得到。

---

## linear model estimation based on pivoted Cholesky factorization with Jacobi preconditioner 
lm.chol <- function(formula, data) { 
    ## stage0: get response vector and model matrix 
    ## we did not follow the normal route: match.call, model.frame, model.response, model matrix, etc 
    y <- data[[as.character(formula[[2]])]] 
    X <- model.matrix(formula, data) 
    n <- nrow(X); p <- ncol(X) 
    ## stage 1: XtX and Jacobi diagonal preconditioner 
    XtX <- crossprod(X) 
    D <- 1/sqrt(diag(XtX)) 
    ## stage 2: pivoted Cholesky factorization 
    R <- suppressWarnings(chol(t(D * t(D * XtX)), pivot = TRUE)) 
    piv <- attr(R, "pivot") 
    r <- attr(R, "rank") 
    if (r < p) { 
    warning("Model is rank-deficient!") 
    piv <- piv[1:r] 
    R <- R[1:r, 1:r] 
    } 
    ## stage 3: solve linear system for coefficients 
    D <- D[piv] 
    b <- D * crossprod(X, y)[piv] 
    z <- forwardsolve(t(R), b) 
    RSS <- sum(y * y) - sum(z * z) 
    sigma <- sqrt(RSS/(n - r)) 
    para <- D * backsolve(R, z) 
    beta.hat <- rep(NA, p) 
    beta.hat[piv] <- para 
    ## stage 4: get standard error 
    Rinv <- backsolve(R, diag(r)) 
    se <- rep(NA, p) 
    se[piv] <- D * sqrt(rowSums(Rinv * Rinv)) * sigma 
    ## stage 5: t-statistic and p-value 
    t.statistic <- beta.hat/se 
    p.value <- 2 * pt(-abs(t.statistic), df = n - r) 
    ## stage 6: construct coefficient summary matrix 
    coefficients <- matrix(c(beta.hat, se, t.statistic, p.value), ncol = 4L) 
    colnames(coefficients) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)") 
    rownames(coefficients) <- colnames(X) 
    ## stage 7: compute adjusted R.squared 
    adj.R2 <- 1 - sigma * sigma/var(y) 
    ## return model fitting results 
    attr(coefficients, "sigma") <- sigma 
    attr(coefficients, "adj.R2") <- adj.R2 
    coefficients 
    } 

这里我要提出三个例子。

---

实施例1：满秩线性模型

我们采取的r内置数据集trees作为一个例子。

# using `lm()` 
summary(lm(Height ~ Girth + Volume, trees)) 
#Coefficients: 
#   Estimate Std. Error t value Pr(>|t|)  
#(Intercept) 83.2958  9.0866 9.167 6.33e-10 *** 
#Girth  -1.8615  1.1567 -1.609 0.1188  
#Volume  0.5756  0.2208 2.607 0.0145 * 
#--- 
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

#Residual standard error: 5.056 on 28 degrees of freedom 
#Multiple R-squared: 0.4123, Adjusted R-squared: 0.3703 
#F-statistic: 9.82 on 2 and 28 DF, p-value: 0.0005868 

## using `lm.chol()` 
lm.chol(Height ~ Girth + Volume, trees) 
#    Estimate Std. Error t value  Pr(>|t|) 
#(Intercept) 83.2957705 9.0865753 9.166905 6.333488e-10 
#Girth  -1.8615109 1.1566879 -1.609346 1.187591e-01 
#Volume  0.5755946 0.2208225 2.606594 1.449097e-02 
#attr(,"sigma") 
#[1] 5.056318 
#attr(,"adj.R2") 
#[1] 0.3702869 

结果完全一样！

---

实施例2：秩亏线性模型

## toy data 
set.seed(0) 
dat <- data.frame(y = rnorm(100), x1 = runif(100), x2 = rbeta(100,3,5)) 
dat$x3 <- with(dat, (x1 + x2)/2) 

## using `lm()` 
summary(lm(y ~ x1 + x2 + x3, dat)) 
#Coefficients: (1 not defined because of singularities) 
#   Estimate Std. Error t value Pr(>|t|) 
#(Intercept) 0.2164  0.2530 0.856 0.394 
#x1   -0.1526  0.3252 -0.469 0.640 
#x2   -0.3534  0.5707 -0.619 0.537 
#x3    NA   NA  NA  NA 

#Residual standard error: 0.8886 on 97 degrees of freedom 
#Multiple R-squared: 0.0069, Adjusted R-squared: -0.01358 
#F-statistic: 0.337 on 2 and 97 DF, p-value: 0.7147 

## using `lm.chol()` 
lm.chol(y ~ x1 + x2 + x3, dat) 
#    Estimate Std. Error t value Pr(>|t|) 
#(Intercept) 0.2164455 0.2529576 0.8556595 0.3942949 
#x1     NA   NA   NA  NA 
#x2   -0.2007894 0.6866871 -0.2924030 0.7706030 
#x3   -0.3051760 0.6504256 -0.4691944 0.6399836 
#attr(,"sigma") 
#[1] 0.8886214 
#attr(,"adj.R2") 
#[1] -0.01357594 
#Warning message: 
#In lm.chol(y ~ x1 + x2 + x3, dat) : Model is rank-deficient! 

这里，lm.chol()基于Cholesky分解基于QR分解具有部分枢转完整的枢转和lm()已缩减不同系数NA。但是两个估计是等价的，具有相同的拟合值和残差。

---

实施例3：用于大型线性模型

n <- 10000; p <- 300 
set.seed(0) 
dat <- as.data.frame(setNames(replicate(p, rnorm(n), simplify = FALSE), paste0("x",1:p))) 
dat$y <- rnorm(n) 

## using `lm()` 
system.time(lm(y ~ ., dat)) 
# user system elapsed 
# 3.212 0.096 3.315 

## using `lm.chol()` 
system.time(lm.chol(y ~ ., dat)) 
# user system elapsed 
# 1.024 0.028 1.056 

lm.chol()性能比lm()快3〜4倍。如果你想知道原因，请阅读my this answer。

---

备注

我的重点是提高对计算内核的性能。通过使用本博尔克的平行主义建议，你可以进一步进一步。如果我的方法提供3倍提升，并行计算提供4倍内核提升3倍，则最终提升9倍！