Splitting and Running Linear Regression - Using data.table

user_10380270 0
Viewed 229

I have asked a previous question (Splitting data and running linear regression loop), and a great solution was provided using tidyverse and pipes. I am saving off the p-values of my parameters for the 4,000 'keys' in the data, putting them into separate dataframe and then running some histograms and other visuals to examine the significance of the parameter for each of the 4,000 keys. This worked the first few times I was running it, but running the exact same code on different parameters/predictors, I continually recieve an error saying:

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

If i run a single model on that key and look at the summary, the pvalue definitely exists in the position [x1,4] or [2,4], but it won't return in the model. Sometimes it will run, but then bomb out on [3,4] or [4,4], etc.

Someone has suggested that tidyverse is killing the memory, thus bombing it out. I know it isn't the code, because it will sometimes work, or sometimes progress further than other times, but it just seems odd. So, would data.table be a better solution to loop over the entire dataset? I am not familiar with how to chain together data.table, so how would I recreate the following code using data.table to run a model over each of 4,000 keys that I have with the 10+ parameters that I have.

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
data.table regression r
3 Answers

This is a relatively straight forward way which assigns intermediate results (e.g., lm(...) and summary(lm(...)) to temporary variables:

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]

The {} allows the intermediate objects to be created. Then we just return a list of what we actually want.

Data:

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

Using the development version of dplyr (shortly to be released as version 1.0 on CRAN) you can do the following:

# 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

Or, to get a data frame in the format of your question:

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

I ran the code above on a 100,000-row data frame with 4,000 levels of Key and didn't run into any memory issues on my 2018 Macbook Pro.

user_4964880

I did a little test comparing OP's current method to a lapply + data.table approach. The operations are performed a 1000 times on a 1000 row data.table with 26 unique keys (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")
)

Outcome

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

There's a big gap in the execution time for these two ways with the lapply method being faster in this test.

Note: I was not able to test out the development version of tidyverse that has the nest_by function (issues with Xcode on my MacOS installation) but it would be worthwhile to include that in the test as well since OP's dataset has 4000 keys.

user_26582280
Related