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