多くの数式をlapplyより速く組み合わせることができますか？

Question

私は、パフォーマンスのために、これを一度にやりたいと思うループを実行するのではなく、データに合わせたい式のリストを持っています。見積もりは依然として別々でなければならない、私はSURか何かを見積もっていない。 次のコードは、私がformulae増加の長さは、本当にこれをベクトル化する方法があるよう

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) 

は、残念ながら、これはやや遅くなる欲しいものですか？

助けがあれば、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行列：今

clusterExport(clust,"df") 

がフルモデルをフィット

は（あなたがすべての変数に合わせて y~.を使用することができます）

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

：

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

クラスタは、すべてのデータセットへのアクセスを必要とするノードすべてのサブモデルに適合する（どのサブモデルが適合するかを制御する多くのオプションについては?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()と呼ばれます。普段summary(lm(...))$coefに見るように、

• 係数サマリー表;：それはformulaとdataを取り、リターン
• 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 
      } 
  

---

は、ここで私は3つの例を提供します。

---

例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! 

、部分ピボット付きQR分解に基づいて完全な旋回とlm()とコレスキー分解に基づいlm.chol()はNAに異なる係数を縮小しています。しかし、2つの推定値は等しく、適合した値と残差が同じです。

---

例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()が3〜4倍高速lm()を超えています。理由を知りたい場合は、my this answerと読んでください。

---

備考

私は、計算カーネルのパフォーマンスを向上させることに焦点を当てています。 Ben Bolkerの並列性の提案を使用すると、さらに一歩前進することができます。私のアプローチで3倍のブーストが得られ、並列コンピューティングで4コアで3倍のブーストが得られれば、9回のブーストが可能になります！