适合许多glm模型：提高速度

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

fastglm 的 method 参数允许改变分解。选项 2 给出了香草 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 更不稳定。

适合许多glm模型：提高速度

Fit many glm models: improve speed

r

glm