适合多种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迭代次数更可取。fastglm和标准IRLS实现的另一个显著区别是，它谨慎地使用了半步，以防止发散(在许多情况下，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 

关于使用扫帚和快速的物体，我可以调查。

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