r - 适合许多 glm 模型 : improve speed

Question

我正在编写一个适合许多 glm 模型的函数。为了给你一些关于这个功能的想法，我包含了一小部分我的代码。在几个 SO 用户的帮助下，该功能现在可以用于我的分析目的。然而，有时，特别是当样本量相对较小时，可能需要很长时间才能完成整个过程。为了减少时间，我正在考虑更改迭代最大化的一些细节，例如最大迭代次数。我还没有找到一种方法，可能是因为我仍然不熟悉 R 术语。任何这样做的建议或其他减少时间的方法都将不胜感激。

all_glm <- function(crude, xlist, data, family = "binomial", ...) {
  # md_lst include formula for many models to be fitted  
  comb_lst <- unlist(lapply(1:n, function(x) combn(xlist, x, simplify=F)), recursive=F)
  md_lst   <- lapply(comb_lst,function(x) paste(crude, "+", paste(x, collapse = "+")))
  models  <- lapply(md_lst, function(x) glm(as.formula(x), family = family, data = data))
  OR      <- unlist(lapply(models, function(x) broom::tidy(x, exponentiate = TRUE)$estimate[2]))

}    

编辑感谢 @BenBolker 指导我使用 fastglm 包，我最终得到了几个 r 包，它们可以为 glm 提供更快的替代方案。我试过 fastglm 和 speedglm。在我的机器上，两者似乎都比 glm 快。

library(fastglm)
library(speedglm)
# from 
set.seed(1)
n <- 25000
k <- 500
y <- rbinom(n, size = 1, prob = 0.5)
x <- round( matrix(rnorm(n*k),n,k),digits=3)
colnames(x) <-paste("s",1:k,sep = "")
df <- data.frame(y,x)
fo <- as.formula(paste("y~",paste(paste("s",1:k,sep=""),collapse="+")))   

# Fit three models: 
system.time(m_glm <- glm(fo, data=df, family = binomial))
system.time(m_speedglm <- speedglm(fo, data= df, family = binomial()))
system.time(m_fastglm <- fastglm(x, y, family = binomial()))

> system.time(m_glm <- glm(fo, data=df, family = binomial))
   user  system elapsed 
  56.51    0.22   58.73 
> system.time(m_speedglm <- speedglm(fo, data= df, family = binomial()))
   user  system elapsed 
  17.28    0.04   17.55 
> system.time(m_fastglm <- fastglm(x, y, family = binomial()))
   user  system elapsed 
  23.87    0.09   24.12 

Answer 1

通常用于拟合 glms 的 IRLS 算法需要在每次迭代时进行矩阵求逆/分解。 fastglm 提供了几种不同的分解选项，默认选择是一个较慢但更稳定的选项(QR with column-pivoting)。如果您只对速度感兴趣，那么两种可用的 Cholesky 型分解中的一种将显着提高速度，这比仅更改 IRLS 迭代次数更可取。 fastglm 和标准 IRLS 实现之间的另一个显着区别是它谨慎使用半步以防止发散(在许多情况下，IRLS 在实践中可能会发散)。

fastglm 的method 参数允许改变分解。选项 2 给出了 Vanilla Cholesky 分解，选项 3 给出了一个稍微更稳定的版本。在我的电脑上，您提供的示例的时间是:

> system.time(m_glm <- glm(fo, data=df, family = binomial))
   user  system elapsed 
 23.206   0.429  23.689 

> system.time(m_speedglm <- speedglm(fo, data= df, family = binomial()))
   user  system elapsed 
 15.448   0.283  15.756 

> system.time(m_fastglm <- fastglm(x, y, family = binomial(), method = 2))
   user  system elapsed 
  2.159   0.055   2.218 

> system.time(m_fastglm <- fastglm(x, y, family = binomial(), method = 3))
   user  system elapsed 
  2.247   0.065   2.337 

关于将 broom 与 fastglm 对象一起使用，我可以研究一下。

关于分解的另一个注意事项:当 fastglm 使用 QR 分解时，它直接使用设计矩阵。虽然 speedglm 在技术上提供了 QR 分解，但它的工作原理是首先计算 $X^TX$ 并对其进行分解，这在数值上比 X 上的 QR 更不稳定。