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

method的参数允许fastglm人们改变分解。选项 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 更不稳定。