Подходит много формул одновременно, более быстрые варианты, чем lapply?

У меня есть список формул, которые я хочу вписать в данные, вместо того, чтобы запускать цикл, я бы хотел сделать это сразу, ради производительности. Оценки все еще должны быть отдельными, я не пытаюсь оценить 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 Меня волнуют коэффициенты и некоторые стандартные ошибки.

 Hack-R27 июл. 2016 г., 02:26
Вы всегда можете попробовать использовать распараллеливание
 sheß27 июл. 2016 г., 02:17
Извините, я уточнил это. Я имею в виду, это единственные результаты, которые меня волнуют, в основном то, что я получаю, когда делаюsummary(results)$coefficients
 sheß27 июл. 2016 г., 02:24
Не совсем уверен, что я понимаю, что вы говорите
 Sumedh27 июл. 2016 г., 02:53
Вы можете посмотреть наэтот
 shayaa27 июл. 2016 г., 02:50
Всегда есть RcppEigen :: fastLm

Ответы на вопрос(2)

Как я уже сказал в своем комментарии, вам действительно нужна более эффективная, но стабильная подпрограмма, отличная отlm(), Здесь я хотел бы предоставить вам хорошо проверенный написанный мной, названныйlm.chol(), Требуетсяformula а такжеdataи возвращает:

Сводная таблица коэффициентов, как вы обычно видите вsummary(lm(...))$coef;Оценка Пирсона остаточной стандартной ошибки, полученная из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: линейная модель полного ранга

Мы берем встроенный набор данных Rtrees В качестве примера.

# 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() основанный на факторизации Холецкого с полным поворотом иlm() на основе QR-факторизации с частичным поворотом сократили различные коэффициенты до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() в 3 ~ 4 раза быстрее, чемlm(), Если вы хотите узнать причину, прочитайтемой этот ответ.

замечание

Я сосредоточился на улучшении производительности вычислительного ядра. Вы можете сделать еще один шаг, используя предложение параллельности Бена Болкера. Если мой подход дает увеличение в 3 раза, а параллельные вычисления дают увеличение в 3 раза на 4 ядрах, то в итоге вы получите увеличение в 9 раз!

Там не очень простой способ векторизовать это, ноpdredge функция отMuMIn Пакет дает вам довольно простой способ распараллелить его (предполагается, что у вас есть несколько ядер на вашем компьютере или вы можете настроить локальный кластер одним из способов, поддерживаемых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

Ваш ответ на вопрос