拆分和运行线性回归 - 使用 data.table答案

【问题标题】：Splitting and Running Linear Regression - Using data.table拆分和运行线性回归 - 使用 data.table
【发布时间】：2020-04-08 00:19:29
【问题描述】：

我已经问过一个先前的问题 (Splitting data and running linear regression loop)，使用 tidyverse 和管道提供了一个很好的解决方案。我正在为数据中的 4,000 个“键”保存参数的 p 值，将它们放入单独的数据框中，然后运行一些直方图和其他视觉效果来检查 4,000 个键中每个键的参数的重要性。这在我运行它的前几次有效，但是在不同的参数/预测器上运行完全相同的代码，我不断收到一条错误消息：

Error in summary(lm(y1 ~ x1 + x2 ))$coefficients['x1', : subscript out of bounds

如果我在该键上运行单个模型并查看摘要，pvalue 肯定存在于位置[x1,4] 或[2,4]，但它不会在模型中返回。有时它会运行，但随后会在 [3,4] 或 [4,4] 等上炸毁。

有人建议 tidyverse 正在杀死内存，从而将其炸毁。我知道这不是代码，因为它有时会起作用，或者有时会比其他时候进步得更远，但这似乎很奇怪。那么，data.table 会是循环整个数据集的更好解决方案吗？我不熟悉如何将data.table 链接在一起，那么我将如何使用data.table 重新创建以下代码，以在我拥有的 4,000 个键中的每一个上运行一个模型，并使用我拥有的 10 多个参数。

df
Key y1 x1 x2
A   10 1  3
A   11 2  4 
A   12 3  5
B   13 4  6 
B   14 5  7
B   15 6  8
C   16 7  9 
C   17 8  1
C   18 9  2

df %>% group_by(Key) %>%
  summarise(Intercept = lm(y1 ~ x1 + x2)$coefficients[1],
            Coeff_x1 = lm(y1 ~ x1 + x2)$coefficients[2],
            Coeff_x2 = lm(y1 ~ x1 + x2)$coefficients[3],
            R2 = summary(lm(y1 ~ x1 + x2))$r.squared,
            pvalue = summary(lm(y1 ~ x1 + x2))$coefficients["x1",4])

# A tibble: 3 x 6
  Key   Intercept Coeff_x1  Coeff_x2    R2     pvalue
  <chr>     <dbl>    <dbl>     <dbl> <dbl>      <dbl>
1 A            9.     1.00 NA            1   8.00e-16
2 B            9.     1.00 NA            1   7.00e-16
3 C            9.     1.00  7.86e-16     1 NaN

【问题讨论】：

您可以使用lapply 循环遍历key 的唯一值，并让函数返回具有所需系数的data.frame
这是上一篇文章中的建议，但我收到了一个错误Error: cannot allocate vector of size XX.x Mb，所以我认为我的机器可能太弱了……这太疯狂了，因为它非常高端18 个月前....
好吧，让我测试一下。

标签： r data.table regression

【解决方案1】：

这是一种将中间结果（例如，lm(...) 和 summary(lm(...)) 分配给临时变量）的相对直接的方式：

dt[, {LM = lm(y1 ~ x1 + x2, data = .SD)
      LM.summary = summary(LM)
    list(Intercept = LM$coefficients[1],
         Coeff_x1 = LM$coefficients[2],
         Coeff_x2 = LM$coefficients[3],
         R2 = LM.summary$r.squared,
         pvalue.x1 = LM.summary$coefficients["x1", 4],
         pvalue.x2 = LM.summary$coefficients["x2", 4])
    },
   by = Key]

{} 允许创建中间对象。然后我们只返回一个我们真正想要的列表。

数据：

library(data.table)
dt = fread('Key y1 x1 x2
A   10 1  3
A   11 2  4 
A   12 3  5
A   13 4  5
B   13 4  6 
B   14 5  7
B   15 6  8
B   15 5  9
C   16 7  9 
C   17 8  1
C   18 9  2
C   18 9  2')

【讨论】：

我收到一个错误Error in [.tbl_df(dat, , { : unused argument (by = Key) 我需要做些什么来初始化data.table吗？我假设dt 是您引用的数据框？
我使用了fread，它直接将您的表格读取为data.table。您需要将library(data.table); setDT(dt) 添加到您的代码中
哇，运行速度非常快。非常感谢你的建议。但是，当我再添加一行以获得 x2 的 pvalue 时，它告诉我它超出了范围。即使我尝试说pvalue = LM.summary$coefficients[3, 4])，它也会说它越界。但是，如果我运行pvalue = LM.summary$coefficients[2, 4])，它将正确创建我的 pvalues 的重复列（与“x1”相比）。我不明白为什么除了截距和第一个参数之外，我无法从线性模型摘要中提取 pvalues？
我认为这与每个分组只有3条记录并且契合度非常好有关。请参阅summary(lm(y1 ~ x1 + x2, data = list(y1 = 10:12, x1 = 1:3, x2 = 3:5)))$coefficients - 它也只提供两个系数。

【解决方案2】：

使用dplyr 的开发版本（不久将在 CRAN 上发布为 1.0 版），您可以执行以下操作：

# devtools::install_github("tidyverse/dplyr")

library(tidyverse)

res = df %>% 
  nest_by(Key) %>%
  mutate(model=list(lm(y1 ~ x1 + x2, data=data)))

res %>% 
  summarise(broom::tidy(model))

  Key   term        estimate  std.error statistic    p.value
  <fct> <chr>          <dbl>      <dbl>     <dbl>      <dbl>
1 A     (Intercept) 9.00e+ 0   2.71e-15   3.32e15   1.92e-16
2 A     x1          1.00e+ 0   1.26e-15   7.96e14   8.00e-16
3 B     (Intercept) 9.00e+ 0   5.57e-15   1.62e15   3.94e-16
4 B     x1          1.00e+ 0   1.10e-15   9.10e14   7.00e-16
5 C     (Intercept) 9.00e+ 0 NaN        NaN       NaN       
6 C     x1          1.00e+ 0 NaN        NaN       NaN       
7 C     x2          7.86e-16 NaN        NaN       NaN

res %>% 
  summarise(broom::glance(model))

  Key   r.squared adj.r.squared      sigma statistic    p.value    df logLik   AIC   BIC deviance df.residual
  <fct>     <dbl>         <dbl>      <dbl>     <dbl>      <dbl> <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
1 A             1             1   1.78e-15   6.34e29   8.00e-16     2   99.3 -193. -195. 3.16e-30           1
2 B             1             1   1.55e-15   8.28e29   7.00e-16     2   99.7 -193. -196. 2.42e-30           1
3 C             1           NaN NaN        NaN       NaN            3  Inf   -Inf  -Inf  0.                 0

或者，以您的问题的格式获取数据框：

library(broom)

res %>% 
  summarise(tidy(model), glance(model)) %>% 
  select(Key, term, estimate, r.squared, p.value) %>% 
  pivot_wider(names_from=term, values_from=estimate)

  Key   r.squared    p.value `(Intercept)`    x1        x2
  <fct>     <dbl>      <dbl>         <dbl> <dbl>     <dbl>
1 A             1   8.00e-16            9.  1.00 NA       
2 B             1   7.00e-16            9.  1.00 NA       
3 C             1 NaN                   9.  1.00  7.86e-16

我在 100,000 行数据帧上运行了上面的代码，其中包含 4,000 个级别的 Key，并且在我的 2018 Macbook Pro 上没有遇到任何内存问题。

【讨论】：

【解决方案3】：

我做了一个小测试，将 OP 的当前方法与 lapply + data.table 方法进行比较。这些操作在 1000 行 data.table 上执行 1000 次，具有 26 个唯一键 (keycol)：

set.seed(28)
dat <- data.table(keycol = sample(x = LETTERS, size = 1000, replace = T), 
                  x = rnorm(n = 1000, mean = 30, sd = 2), 
                  y = rnorm(n = 1000, mean = 20, sd = 2), 
                  z = rnorm(n = 1000, mean = 10, sd = 2))

speed_test <- benchmark(
  'data_table' = {
    model_list <- lapply(X = 1:26, function(z){  #X could be the unique keys or the 1:length(unique(keys))
      m <- lm(data = dat[keycol == LETTERS[z], ], formula = x ~ y + z)
      smry <- summary(m)
      ret_tbl <- data.table(intercept = smry$coefficients[1],
                            coef_y = smry$coefficients[2], 
                            coef_z = smry$coefficients[3],
                            r_squared = smry$adj.r.squared, 
                            pvale = smry$coefficients[2,4], 
                            keycol = z) 
      return(ret_tbl)

    })
    desired_tbl <- rbindlist(l = model_list, use.names = T, fill = T)
  }, 

  'tidyverse1' = {
    dat %>% group_by(keycol) %>%
      summarise(Intercept = lm(x ~ y + z)$coefficients[1],
                Coeff_y = lm(x ~ y + z)$coefficients[2],
                Coeff_z = lm(x ~ y + z)$coefficients[3],
                R2 = summary(lm(x ~ y + z))$r.squared,
                pvalue = summary(lm(x ~ y + z))$coefficients["y",4])
  }, 
  replications = 1000,
  columns = c("test", "replications", "elapsed")
)

结果

> speed_test
        test replications elapsed
1 data_table         1000  29.477
2 tidyverse1         1000  88.781

这两种方法的执行时间存在很大差距，lapply 方法在此测试中更快。

注意：我无法测试具有 nest_by 功能的 tidyverse 的开发版本（我的 MacOS 安装上的 Xcode 问题），但由于 OP 将其包含在测试中也是值得的数据集有 4000 个键。

【讨论】：