运行一百万个回归模型的更快方法

【问题标题】：Faster way to run a million regression models运行一百万个回归模型的更快方法
【发布时间】：2018-02-15 08:09:39
【问题描述】：

我使用基因组数据，我经常需要运行一百万个或更多的回归模型。我在下面的循环可以工作，但是速度很慢，并且由于附加每条新记录的开销而持续时间越长越慢。

#### setup sample data ###
require(data.table)
data <- data.frame(
  C = rnorm(10, 5),
  D = rnorm(10, 7),
  E = rnorm(10, 9),
  A = rnorm(10, 1),
  B = rnorm(10, 3)
)
outcome <- c(rnorm(10, 5))
cov <- data.frame(cov1 = c(1, 1, 1, 2, 2, 1, 1, 1, 2, 2))
#### initialize results file ###
myresults <- data.table(NULL)
#### Run regression against same covariates and outcome for each column in data ##
for (i in 1:ncol(data)) {
  id = colnames(data)[i]
  mydata <- cbind(cov, outcome, data[, id])
  colnames(mydata)[ncol(mydata)] <- id #I can't figure out how to not have to do this
  fit <-
    glm(formula(paste0("outcome ~ as.factor(cov1) + ", id)), data = mydata)
  myresults <- rbindlist(list(
    myresults,
    data.table(
      id = id,
      estimate = signif(coef(summary(fit))[id, "Estimate"], digits = 4),
      pvalue = signif(coef(summary(fit))[id, "Pr(>|t|)"], digits = 4)
    )
  ))
}
myresults

这给出了我所需要的输出结果文件。我可以修改它以添加其他输出，在循环中运行其他模型以按协变量进行分层，然后捕获等...我的输出始终具有与初始 data 中的列相同数量的数据行数。 

   id estimate  pvalue
1:  C -0.22220 0.49230
2:  D  0.64550 0.08568
3:  E -0.06756 0.83990
4:  A  0.39750 0.54060
5:  B -0.34300 0.35410

通过切换到您在循环中看到的data.table::rbindlist，我能够获得一些改进。

我一直在尝试使用lapply(split(data, colnames(data))) 之类的东西来看看我是否可以加快速度，甚至可能使用mclapply()，但一直无法让它工作。

非常感谢您的帮助。

编辑：我支持所有回复的人，因为他们都很有帮助，我很感激所花费的时间。

以 6 倍的优势明显获胜的是 Roland 的评论。我列出了我在这里为后代所做的事情，以防它可以帮助其他人。

我合并到一个非常宽的数据集 (260 x 470,000)

require(data.table)
require(reshape2)
bigdata <- cbind(mycovs, testdata)
test <- data.table(bigdata)

然后我把它做成了一个高大的数据集：

    DT.m1 = melt(
  test,
  id.vars = c(
    "Sample_Plate",
    "BaseName",
    "Race",
    "Education",
    "mom_age_delv",
    "sex",
    "gest_age_wks",
    "MONTH_OLD",
    "DEPRESSION",
    "CD8T",
    "CD4T",
    "NK",
    "Bcell",
    "Mono",
    "Gran"
  ),
  measure.vars = c(16:ncol(test)),
  variable.name = "cpg",
  value.name = "betaval"
)

然后我运行我的完整回归模型并从系数表中取出最后一行，如下所示：

system.time(res <-
              DT.m1[, {
                fit <-
                  glm(
                    DEPRESSION ~ as.factor(Sample_Plate) + as.factor(sex) + as.factor(Education) + as.factor(Race) + MONTH_OLD + mom_age_delv + gest_age_wks + CD8T + CD4T + NK + Bcell + Mono + Gran + betaval,
                    data = .SD
                  )
                coef(summary(fit))[nrow(coef(summary(fit))), c(1, 2, 4)]
              }, by = cpg])

最后，我清理了它。

res <- cbind(res, c("beta1", "se", "pvalue"))
head(res)
final_results <- dcast(data = res, cpg ~ V2, value.var = "V1")[c(1, 2, 4, 3)]

这导致每 1000 个模型的时间约为 10 秒。下一个关闭时间是~60 秒。

清理部分似乎应该可以在 data.table() 中完成，但我无法弄清楚。我只能为我请求的每个 coef 列重复一个高 2 列向量。

如果您对如何改进有其他想法，请告诉我，再次感谢。

【问题讨论】：

  • 您可能会研究并行性、多核和 GPU 策略等。不过，这对于 Stack Overflow 来说有点宽泛。您可能想尝试 Code Review SE。 codereview.stackexchange.com
  • 这看起来很可疑。您可能可以做一些比产生这些无用 p 值更明智的事情。无论如何，第一步：使您的数据成为 data.table 并重塑为长格式，第二步：DT[, {fit ; }，按 = 变量]。根据需要进行调整。
  • 感谢您的意见@Roland。关于有效性，阅读 GWAS SNP 测试或 CpG 甲基化标记测试，并使用调整后的 pvalues 和 qqplots 来确定显着性；如果您愿意，我们可以在某个时候进行边栏聊天。

标签： r loops split regression apply


【解决方案1】：

一个老问题，但看到我是优化的傻瓜......

OP 的更新解决方案让他们将结果嵌套在 data.table 中；基本上通过分组变量运行每个回归。这是一个很好的策略，尤其是因为data.table 会自动为您处理并行性。

但是，通过每次运行完整的 glm() 调用，您仍然要求 R 完成比它需要的更多的工作。您只需要系数（和 se + p 值）。 glm() 包括的所有其他内容——model.matrix 等——只是浪费精力和内存开销。不仅如此，我还不清楚为什么您首先要运行 generalised 模型，因为您没有调用非高斯族（二项式等）。调用标准线性模型（即lm()）也会自动减少开销。

那么，解决方案是什么？切换到矩阵数组并改用 lm.fit()。这是lm() 在底层调用的主力函数，并且更精简/更快。事实上，还有一个更精简/更快的.lm.fit() 版本，但这会使计算 SEs + p 值变得有点棘手。 （更多内容见下文。）根据 OP 的 MWE 改编的单个回归的快速基准测试：

set.seed(42)
data <- data.frame(
    C = rnorm(10, 5),
    D = rnorm(10, 7),
    E = rnorm(10, 9),
    A = rnorm(10, 1),
    B = rnorm(10, 3)
)
outcome <- c(rnorm(10, 5))
cov <- data.frame(cov1 = as.factor(c(1, 1, 1, 2, 2, 1, 1, 1, 2, 2)))

id = "A"
mydata = cbind(cov, outcome, x = data[, id])
# design matrix for *.fit versions
X = cbind(intercept = 1, cov1 = cov$cov1, x = data[, id])

microbenchmark::microbenchmark(
    .lm.fit = .lm.fit(X, outcome),
    lm.fit  = lm.fit(X, outcome),
    glm.fit = glm.fit(X, outcome),
    lm      = lm(outcome ~ cov1 + x, data = mydata),
    glm     = glm(outcome ~ cov1 + x, data = mydata),
    times   = 100
)
#> Unit: microseconds
#>     expr     min      lq      mean   median       uq      max neval  cld
#>  .lm.fit   1.684   3.479   4.92692   4.1815   4.8915   33.607   100 a   
#>   lm.fit  16.222  24.195  27.62058  27.5270  30.4225   51.694   100 a   
#>  glm.fit 135.641 149.130 195.54319 164.6615 179.8900 1920.694   100  b  
#>       lm 531.621 554.632 658.56501 632.0585 663.1655 2552.081   100   c 
#>      glm 702.767 724.999 855.46429 850.9145 863.4710 2418.207   100    d

由reprex package (v2.0.1) 于 2021 年 8 月 30 日创建

因此，您只需使用更精简的数组类型（矩阵）和回归调用 (lm.fit/.lm.fit)，即可获得 30–200 倍 (!) 的加速。请注意，您仍然可以将嵌套的 data.table 方法与矩阵对象一起使用，以获得两全其美的效果。这需要更多解释，但您可以看到一个示例 here。

后记从 lm.fit 对象中获取标准错误和 p 值。

lm.fit 方法和 OP 所需输出的唯一轻微复杂性是我们不会自动获得 SE 和 p 值。但是，编写我们自己的函数（几乎是从summary.lm 中逐字删除）为我们执行此操作很容易。

se_pval = function(object) {
    z <- object
    p <- z$rank
    rdf <- z$df.residual
    Qr <- z$qr
    r = z$residuals
    rss <- sum(r^2)
    resvar <- rss/rdf
    p1 <- 1L:p
    R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
    se <- sqrt(diag(R) * resvar)
    est <- z$coefficients[Qr$pivot[p1]]
    tval <- est/se
    pval = 2 * pt(abs(tval), rdf, lower.tail = FALSE)
    se_pval = rbind(se, pval); names(se_pval) = colnames(pval)
    return(se_pval)
}


fit = lm.fit(X, outcome)
c(coef(fit)['x'], se_pval(fit)[, 'x'])
#>         x        se      pval 
#> 0.5784677 0.3230284 0.1164475

## Compare with glm.summary
coef(summary(glm(outcome ~ cov1 + x, data = mydata)))['x', c(1,2,4)]
#>   Estimate Std. Error   Pr(>|t|) 
#>  0.5784677  0.3230284  0.1164475

由reprex package (v2.0.1) 于 2021-08-30 创建

更新

按照 cmets 中的要求添加完全工作的迭代。

se_pval = function(object) {
    z <- object
    p <- z$rank
    rdf <- z$df.residual
    Qr <- z$qr
    r = z$residuals
    rss <- sum(r^2)
    resvar <- rss/rdf
    p1 <- 1L:p
    R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
    se <- sqrt(diag(R) * resvar)
    est <- z$coefficients[Qr$pivot[p1]]
    tval <- est/se
    pval = 2 * pt(abs(tval), rdf, lower.tail = FALSE)
    se_pval = rbind(se, pval); names(se_pval) = colnames(pval)
    return(se_pval)
}


set.seed(42)
Xlist <- list(
    C = rnorm(10, 5),
    D = rnorm(10, 7),
    E = rnorm(10, 9),
    A = rnorm(10, 1),
    B = rnorm(10, 3)
)
outcome <- c(rnorm(10, 5))
cov1 = as.factor(c(1, 1, 1, 2, 2, 1, 1, 1, 2, 2))

regs = lapply(
    names(Xlist),
    function(id) {
        X = cbind(intercept = 1, cov1 = cov1, x = Xlist[[id]])
        fit  = lm.fit(X, outcome)
        # c(id = id, coef(fit)['x'], se_pval(fit)[, 'x']) ## Could add ID here..
        c(coef(fit)['x'], se_pval(fit)[, 'x'])
        })
regs = signif(do.call("rbind", regs), digits = 4) ## numeric matrix

## Easy to add ID after the fact and coerce to e.g. a data.frame
regs = data.frame(id = names(Xlist), regs)
#>   id       x     se   pval
#> 1  C -0.4300 0.4957 0.4145
#> 2  D  0.1251 0.2590 0.6439
#> 3  E  0.4782 0.4680 0.3409
#> 4  A  0.5785 0.3230 0.1164
#> 5  B  0.1131 0.5511 0.8432

【讨论】：

  • 感谢您的跟进。我会对此进行一些测试。我提出了其他解决方案，我只需自己翻转最快的矩阵。您的解决方案虽然跳过了这样一个事实，即要优化的重要事情之一也是存储结果（字面意思是数百万）以及这样做所涉及的 IO。对吗？
  • 是的，没错……尽管存储这些结果的内存开销将是最小的，因为一切都非常精简。如果是我，我会将我所有的潜在解释变量放在一个大矩阵或列表中，然后遍历（或随机抽样）列以提取每次特定迭代的设计矩阵。我会更新我的答案来举个例子。
【解决方案2】：

你应该试试map和map_df

library(tidyverse)
myfun <- function(data, outcome, cov) {
    require(tidyverse)
      numcol <- ncol(data)
    newdata <- data %>%
                mutate(outcome = outcome, cov = cov$cov1)

    fmla <- map(names(newdata[,1:numcol]), ~glm(formula(paste0("outcome ~ as.factor(cov) + ", .x)), data=newdata))
    ans <- map_df(fmla, ~as_tibble(matrix(coef(summary(.x))[2, c(1,4)], ncol=2, byrow=TRUE)), .id="id") %>%
               rename(estimate=V1, pvalue=V2)
}

基准

拥有更大的数据

biggerdata <- as_tibble(matrix(rnorm(2000), nrow=10))
library(microbenchmark)
microbenchmark(myfun(biggerdata,outcome,cov), OP(biggerdata,outcome,cov))

Unit: milliseconds
                            expr      min         lq       mean     median
 myfun(biggerdata, outcome, cov)   71.534   72.98252   77.82994   76.31598
    OP(biggerdata, outcome, cov) 1936.986 1994.03518 2048.96934 2018.33299
         uq       max neval
   79.97554  106.9852   100
 2085.44655 2297.3878   100

OP函数

OP <- function(data, outcome, cov) {
    myresults <- data.table(NULL)
    #### Run regression against same covariates and outcome for each column in data ##
    for (i in 1:ncol(data)) {
        id = colnames(data)[i]
        mydata <- cbind(cov, outcome, data[, id])
        colnames(mydata)[ncol(mydata)] <- id #I can't figure out how to not have to do this
        fit <- glm(formula(paste0("outcome ~ as.factor(cov1) + ", id)), data = mydata)
        myresults <- rbindlist(list(
                        myresults,
                        data.table(
                              id = id,
                              estimate = signif(coef(summary(fit))[id, "Estimate"], digits = 4),
                              pvalue = signif(coef(summary(fit))[id, "Pr(>|t|)"], digits = 4)
                        )
        ))
    }
    myresults
}

处理 cov 中多个协变量的新函数

set.seed(20)
newcov <- data.frame(cov1 = sample(c(1,2), 10, replace=TRUE),
            cov2 = sample(c(1,2), 10, replace=TRUE),
            cov3 = sample(c(1,2), 10, replace=TRUE))

mynewfun <- function(data, outcome, cov) {
                require(tidyverse)
                numcol <- ncol(data)
                newdata <- data %>%
                        mutate(outcome = outcome) %>%
                          cbind(cov)

                covname <- names(cov)
                fmla <- map(names(newdata[,1:numcol]), ~glm(formula(paste0("outcome ~ ", paste0(covname, collapse=" + "), " + ", .x)), data=newdata))
                ans <- map_df(fmla, ~as_tibble(matrix(coef(summary(.x))[2, c(1,4)], ncol=2, byrow=TRUE)), .id="id") %>%
                       rename(estimate=V1, pvalue=V2)
                return(ans)
         }

mynewfun(data,outcome,newcov)

【讨论】：

  • data 的每一列都应该是~ 的 RHS 上的一个术语（通过迭代）。我正在研究语法，但看起来data 列上的迭代被建模为自变量？我想要outcome ~ allcovariates + data[,i] 或者我在阅读你的代码时还差得远？
  • 我对此有点困惑。 allcovariates 是什么意思？在cov 中有一个单独的列cov1。我以为你想为每个data[,i] 制作以下公式outcome ~ cov$cov1 + data[,i]。
  • 我只使用了一个协变量作为示例。我的模型如下所示：fit &lt;- glm(formula(paste0(outcome, " ~ as.factor(cov1) + as.factor(cov2) + cov3", data[,i])), data=data) 其中协变量是静态的，结果也是如此。可变部分是名为data 的大数据表的每一列。实际上，我需要一次容纳 467,000 列数据。事实上，结果变量和所有协变量都在一个数据集中，因为它们是静态的，而大数据集是独立的。
  • 你有时间去聊天室吗？
  • 查看更新答案（见底部功能）是否符合您的预期
【解决方案3】：

在我的电脑上，这大约快 40%：

timestart <- Sys.time()

mydata <- cbind(cov, outcome, data)
my.glm <- function (mycol) {
  fit <- glm(eval(parse(text = paste("outcome ~ cov1 +", mycol))), data = mydata)
    data.table(
    id = mycol,
    estimate = signif(coef(summary(fit))[mycol, "Estimate"], digits = 4),
    pvalue = signif(coef(summary(fit))[mycol, "Pr(>|t|)"], digits = 4)
  )
}
(res.l <- do.call(rbind, lapply(colnames(data), my.glm)))

Sys.time() - timestart

【讨论】：

    【解决方案4】：

    我不知道加速是否有意义，但我对您的代码做了一些简化。
    首先，在循环外只调用一次factor 和colnames(data)。

    cov2 <- data.frame(cov1 = factor(c(1, 1, 1, 2, 2, 1, 1, 1, 2, 2)))
#
cnames <- colnames(data)
mydata2 <- data.frame(cov2, outcome, other = NA)

    现在，定义一个供lapply 使用的函数。请注意，此函数使用全局环境中存在的多个数据对象，这通常是一种不好的做法。

    fun <- function(i){
    id <- cnames[i]
    mydata2[, 3] <- data[, id]
    names(mydata2)[3] <- id
    fit <- glm(formula(paste0("outcome ~ cov1 + ", id)), data = mydata2)
    data.table(
        id = id,
        estimate = signif(coef(summary(fit))[id, "Estimate"], digits = 4),
        pvalue = signif(coef(summary(fit))[id, "Pr(>|t|)"], digits = 4)
    )
}

myresults2 <- rbindlist(lapply(seq_len(ncol(data)), fun))

identical(myresults, myresults2)
[1] TRUE

    对象myresults 是您的代码获得的对象。如您所见，结果是相同的。

    【讨论】：

      【解决方案5】：

      每次调用myresults &lt;- rbindlist(list(myresults, ...)) 时，都会复制整个myresults，修改副本，然后让名称指向副本。 R 中循环效率低下的最常见原因是“增长对象”。您知道结果的确切尺寸（ncol(data) 乘以 3），所以就让它开始吧。然后用data.table引用赋值（不可复制）。

      看看这是否有助于提高效率：

      #### initialize results file ###
myresults <- data.table(
  id       = character(length(data)),
  estimate = numeric(  length(data)),
  pvalue   = numeric(  length(data))
)

#### Run regression against same covariates and outcome for each column in data ##
for (i in seq_along(data)) {
  id = colnames(data)[i]
  mydata <- cbind(cov, outcome, data[, id])
  colnames(mydata)[ncol(mydata)] <- id #I can't figure out how to not have to do this
  fit <-
    glm(formula(paste0("outcome ~ as.factor(cov1) + ", id)), data = mydata)
  set(
    myresults,
    i     = i,
    j     = c("id", "estimate", "pvalue"),
    value = list(
      id       = id,
      estimate = signif(coef(summary(fit))[id, "Estimate"], digits = 4),
      pvalue   = signif(coef(summary(fit))[id, "Pr(>|t|)"], digits = 4)
    )
  )
}

      我还将for (i in 1:ncol(data)) 替换为for (i in seq_along(data))，因为当data 没有列时，第一种方式的行为很糟糕。你可能认为它永远不会发生，但这样写循环是个坏习惯。

      【讨论】：

      • 这非常有帮助。在开始运行时每 1000 个模型花费的时间要多几秒钟，但速度始终保持不变，并且大大提高了整体时间。我发布的 data.table 方法在用我的所有协变量制作了一个高数据集之后最终使对象太大。这种方法最终允许我在不到 3 小时的时间内运行约 500k 模型。再次感谢。
      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2020-08-12
      • 2023-03-22
      • 2020-10-16
      相关资源
      最近更新 更多
      热门标签
      Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode