r语言 - 一次拟合多个公式，比 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

正如我在评论中所说，您真正需要的是除lm()之外更有效但更稳定的试穿程序。在这里，我会为您提供一个经过良好测试的自己写的，称为lm.chol().它需要formula和data，并返回：

• 系数汇总表，如您通常在summary(lm(...))$coef中看到的;
• 皮尔逊估计的残差标准误差，如你从summary(lm(...))$sigma得到的;
• 调整后的R.平方，正如你从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!

在这里，基于完全枢轴的 Cholesky 分解lm.chol()和基于部分枢轴的 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倍。如果你想知道原因，请阅读我的这个答案。

---

备注

我专注于提高计算内核的性能。你可以更进一步，通过使用 Ben Bolker 的并行性建议。如果我的方法提供 3 倍的提升，而并行计算在 4 个内核上提供 3 倍的提升，那么您最终会得到 9 倍的提升！

Answer 2

实际上

没有一种简单的方法可以对此进行矢量化，但是 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)

使用命名数据框而不是匿名矩阵会更容易做到这一点：

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