[R] -экстракция коэффициентов модели из вложенного списка (list-столбцы)

Question

Я следую примеру Many Models из книги R for Data Science, чтобы сгенерировать отдельные модели множественной линейной регрессии. Как показано в воспроизводимом примере ниже, мой фрейм данных состоит из трех столбцов: id, val1, val2, val3,. Я могу приспособить линейные модели, используя функцию map, как описано в книге Хэдли. Однако я изо всех сил пытался извлечь коэффициенты из каждой модели и добавить эти столбцы значений обратно во фреймы данных в пределах my.list. То, как в настоящий момент хранятся коэффициенты модели, затрудняет / усложняет вызов других частей моего кода.

Лучшее, что я до сих пор придумал, - это составить список длиной my.list и извлечь коэффициенты путем итерации по каждому фрейму данных в my.list: Name1, Name2, Name3. Это означает, что теперь у меня есть другой список в моей глобальной среде, и coef.list больше не содержит факторы Name1, Name2, Name3 из my.list; теперь они заменены на [[1]], [[2]], [[3]].

Может ли кто-нибудь предложить «более чистый» способ извлечения коэффициентов модели при работе с несколькими моделями? Мой предпочтительный вывод будет просто создать столбец для каждого коэффициента: intercept, val1, val2. Эти столбцы будут помещаться в существующие фреймы данных Name1, Name2, Name3 в пределах my.list, чтобы я мог использовать mutate непосредственно для фреймов данных:

# reproducible example
set.seed(1363)
d1 <- data.frame(id=c("Name1", "Name2", "Name3"),
             val1=c(rnorm(n=15, mean=5)), 
             val2=c(rnorm(n=15, mean=3)), 
             val3=c(rnorm(n=15, mean=8)))

# linear model function
lm.fun <- function(df){
  lm(val3 ~ val1+val2, data = df)
}

# map lm function
d1 <- d1 %>%
  group_by(id) %>%
  nest() %>%
  mutate(model = map(data, lm.fun)) %>%
  unnest(data, .drop = FALSE)

#split data frame by 'id' and convert into list
my.list <- split(d1, d1$id)

# make list of coefficients
coef.list <- list(length(my.list))
for (i in seq_along(my.list)) {
  coef.list[[i]] <- my.list[[i]][["model"]][[1]][["coefficients"]]
}

>head(coef.list, n=1)
[[1]]
(Intercept)        val1        val2 
9.03278337 -0.07096932  0.02119088 

ЖЕЛАЕМЫЙ ВЫХОД

my.list$Name1
id    val1  val2  val3  intc  coef1  coef2
Name1  1      2    3    9.03  -.070  .021   
Name1  3      1    5    9.03  -.070  .021
Name1  2      6    8    9.03  -.070  .021  

Answer 1

Используя sapply, учитывая данные d1 (содержащие столбец model):

coeff <- sapply(d1$model, function(x) return(coefficients(x)))

library(dplyr)
bind_cols(d1, data.frame(t(coeff))) %>% 
     rename_at(6:8, ~ c("intc",  "coef1",  "coef2")) %>% 
          distinct(id, .keep_all = TRUE)

Answer 2

library(tidyverse)

d1%>%
   group_by(id)%>%
   mutate(s=toString(coef(lm(val3~val1+val2))))%>%
   separate(s,c("intercept","coef1","coef2"),sep=",",convert = T)%>%
   arrange(id)
# A tibble: 15 x 7
# Groups:   id [3]
   id     val1  val2  val3 intercept   coef1   coef2
   <fct> <dbl> <dbl> <dbl>     <dbl>   <dbl>   <dbl>
 1 Name1  6.76  2.64  9.85      9.03 -0.0710  0.0212
 2 Name1  6.78  1.52  6.94      9.03 -0.0710  0.0212
 3 Name1  4.14  4.31  8.30      9.03 -0.0710  0.0212
 4 Name1  5.55  2.16  9.97      9.03 -0.0710  0.0212
 5 Name1  6.18  3.64  8.32      9.03 -0.0710  0.0212
 6 Name2  6.08  1.12  9.96      7.33  0.468  -0.488 
 7 Name2  3.54  4.71  6.43      7.33  0.468  -0.488 
 8 Name2  5.66  4.54  8.69      7.33  0.468  -0.488 
 9 Name2  6.88  4.15  7.79      7.33  0.468  -0.488 
10 Name2  4.89  1.27  8.72      7.33  0.468  -0.488 
11 Name3  6.41  4.38  6.22     20.1  -2.15    0.118 
12 Name3  5.06  3.28  9.42     20.1  -2.15    0.118 
13 Name3  6.25  3.16  8.15     20.1  -2.15    0.118 
14 Name3  6.03  3.63  7.78     20.1  -2.15    0.118 
15 Name3  6.51  1.46  5.90     20.1  -2.15    0.118 

Answer 3

На основании описания, если нам нужно создать несколько столбцов с коэффициентами, один из вариантов использует do с group_by. После группировки по 'id' извлеките коэффициенты как list в пределах do и rename столбцов (при необходимости)

library(tidyverse)
d1 %>% 
  group_by(id) %>% 
  do(data.frame(., as.list(coef(lm(val3 ~ val1 + val2, data = .))))) %>%
  rename_at(5:7, ~c("intc", "coef1", "coef2"))
# A tibble: 15 x 7
# Groups:   id [3]
#   id     val1  val2  val3  intc   coef1   coef2
#   <fct> <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl>
# 1 Name1  6.76  2.64  9.85  9.03 -0.0710  0.0212
# 2 Name1  6.78  1.52  6.94  9.03 -0.0710  0.0212
# 3 Name1  4.14  4.31  8.30  9.03 -0.0710  0.0212
# 4 Name1  5.55  2.16  9.97  9.03 -0.0710  0.0212
# 5 Name1  6.18  3.64  8.32  9.03 -0.0710  0.0212
# 6 Name2  6.08  1.12  9.96  7.33  0.468  -0.488 
# 7 Name2  3.54  4.71  6.43  7.33  0.468  -0.488 
# 8 Name2  5.66  4.54  8.69  7.33  0.468  -0.488 
# 9 Name2  6.88  4.15  7.79  7.33  0.468  -0.488 
#10 Name2  4.89  1.27  8.72  7.33  0.468  -0.488 
#11 Name3  6.41  4.38  6.22 20.1  -2.15    0.118 
#12 Name3  5.06  3.28  9.42 20.1  -2.15    0.118 
#13 Name3  6.25  3.16  8.15 20.1  -2.15    0.118 
#14 Name3  6.03  3.63  7.78 20.1  -2.15    0.118 
#15 Name3  6.51  1.46  5.90 20.1  -2.15    0.118 

---

Или мы можем использовать broom пакетные функции. Такие функции, как tidy, glance извлекают интересующие столбцы (и многое другое). После nest сгруппированного набора данных постройте модель, извлеките компоненты с tidy, select интересующими столбцами и unnest

library(broom)
d1 %>% 
   group_by(id) %>%
   nest %>%
   mutate(model1 = map(data, ~ 
                         lm(val3 ~ val1 + val2, data = .) %>%
                              tidy %>%
                              dplyr::select(term, estimate) %>% 
                              spread(term, estimate) %>% 
                              rename_all(~ c("intc", paste0("coef", 1:2))))) %>%  
   unnest(model1) %>%
   unnest(data)
# A tibble: 15 x 7
#   id     intc   coef1   coef2  val1  val2  val3
#   <fct> <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl>
# 1 Name1  9.03 -0.0710  0.0212  6.76  2.64  9.85
# 2 Name1  9.03 -0.0710  0.0212  6.78  1.52  6.94
# 3 Name1  9.03 -0.0710  0.0212  4.14  4.31  8.30
# 4 Name1  9.03 -0.0710  0.0212  5.55  2.16  9.97
# 5 Name1  9.03 -0.0710  0.0212  6.18  3.64  8.32
# 6 Name2  7.33  0.468  -0.488   6.08  1.12  9.96
# 7 Name2  7.33  0.468  -0.488   3.54  4.71  6.43
# 8 Name2  7.33  0.468  -0.488   5.66  4.54  8.69
# 9 Name2  7.33  0.468  -0.488   6.88  4.15  7.79
#10 Name2  7.33  0.468  -0.488   4.89  1.27  8.72
#11 Name3 20.1  -2.15    0.118   6.41  4.38  6.22
#12 Name3 20.1  -2.15    0.118   5.06  3.28  9.42
#13 Name3 20.1  -2.15    0.118   6.25  3.16  8.15
#14 Name3 20.1  -2.15    0.118   6.03  3.63  7.78
#15 Name3 20.1  -2.15    0.118   6.51  1.46  5.90

---

или без использования broom

d1 %>%
    group_by(id) %>%
    nest %>%
    mutate(model = map(data, ~ lm(val3 ~ val1 + val2, data = .) %>%
                         coef %>% 
                         as.list %>%
                         as_tibble)) %>% 
    unnest(model) %>%
    unnest(data)

---

Если нам нужен только обобщенный вывод для каждого идентификатора

d1 %>% 
  group_by(id) %>% 
  nest %>%
  mutate(model1 = map(data, ~ 
                      lm(val3 ~ val1 + val2, data = .) %>%
                         tidy %>%
                         dplyr::select(term, estimate) %>%
                         spread(term, estimate) %>% 
                         rename_all(~ c("intc", paste0("coef", 1:2))))) %>%  
  dplyr::select(-data) %>% 
  unnest

---

Или с помощью data.table, мы можем сделать это более кратко, после группировки по 'id' и назначения (:=) новых столбцов, извлекая coef модели как list

library(data.table)
setDT(d1)[, c('intc', 'coef1', 'coef2') := 
           as.list(coef(lm(val3 ~ val1 + val2))), id]
d1[order(id)]
#       id     val1     val2     val3      intc       coef1       coef2
# 1: Name1 6.755964 2.642874 9.849828  9.032783 -0.07096932  0.02119088
# 2: Name1 6.776666 1.522431 6.937053  9.032783 -0.07096932  0.02119088
# 3: Name1 4.141883 4.307537 8.301940  9.032783 -0.07096932  0.02119088
# 4: Name1 5.551850 2.163882 9.971588  9.032783 -0.07096932  0.02119088
# 5: Name1 6.179506 3.635832 8.319042  9.032783 -0.07096932  0.02119088
# 6: Name2 6.083243 1.116293 9.960934  7.325156  0.46840770 -0.48806159
# 7: Name2 3.536476 4.708967 6.427627  7.325156  0.46840770 -0.48806159
# 8: Name2 5.663909 4.541081 8.691523  7.325156  0.46840770 -0.48806159
# 9: Name2 6.883746 4.150780 7.791050  7.325156  0.46840770 -0.48806159
#10: Name2 4.890291 1.269559 8.723792  7.325156  0.46840770 -0.48806159
#11: Name3 6.414915 4.383609 6.220188 20.106581 -2.14601530  0.11770877
#12: Name3 5.059774 3.276510 9.421862 20.106581 -2.14601530  0.11770877
#13: Name3 6.251416 3.157157 8.147720 20.106581 -2.14601530  0.11770877
#14: Name3 6.028100 3.630858 7.783118 20.106581 -2.14601530  0.11770877
#15: Name3 6.505153 1.460564 5.895564 20.106581 -2.14601530  0.11770877

Или сделайте join, если мы не хотим обновлять исходные данные

setDT(d1)[, as.list(coef(lm(val3 ~ val1 + val2))), id][d1, on = .(id)]

ПРИМЕЧАНИЕ: 'd1' - начальный набор данных