Расщепление и запуск линейной регрессии - Использование data.table

Question

Я задал предыдущий вопрос ( Разделение данных и выполнение линейной регрессии l oop), и было предложено отличное решение с использованием Tidyverse и каналов. Я сохраняю p-значения моих параметров для 4000 «ключей» в данных, помещаю их в отдельный фрейм данных, а затем запускаю некоторые гистограммы и другие визуальные элементы, чтобы проверить значение параметра для каждого из 4000 ключей. Это работало первые несколько раз, когда я запускал его, но при запуске одного и того же кода на разных параметрах / предикторах я постоянно получаю сообщение об ошибке:

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

Если я запускаю один Для модели на этом ключе и посмотрите на сводку, pvalue определенно существует в позиции [x1,4] или [2,4], но не вернется в модель. Иногда он запускается, но затем бомбить на [3,4] или [4,4] и т. Д. c.

Кто-то предположил, что tidyverse убивает память, тем самым взрывая ее. Я знаю, что это не код, потому что он иногда будет работать, а иногда и продвигаться дальше, чем в другие времена, но это кажется странным. Итак, будет ли data.table лучшим решением для l oop по всему набору данных? Я не знаком с тем, как связать воедино data.table, поэтому как мне воссоздать следующий код, используя data.table, чтобы запустить модель для каждого из 4000 имеющихся у меня ключей с параметрами 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

Answer 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')

Answer 2

Используя разрабатываемую версию dplyr (вскоре будет выпущена как версия 1.0 в CRAN), вы можете сделать следующее:

# 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 фрейм данных с 4000 уровнями Key и не столкнулся с какими-либо проблемами с памятью на моем MacBook Pro 2018 года.

Answer 3

Я провел небольшой тест, сравнивая текущий метод OP с подходом lapply + data.table. Операции выполняются 1000 раз на 1000 строк data.table с 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 был быстрее.

Примечание. Мне не удалось протестировать версию tidyverse для разработки, которая имеет nest_by функция (проблемы с Xcode на моей установке MacOS), но было бы целесообразно включить это в тест, так как набор данных OP имеет 4000 ключей.