Подгонка модели с использованием каждого столбца предиктора индивидуально сохраняет результаты во фрейме данных

У меня есть фреймворк с одним столбцом переменной ответа и несколькими столбцами переменных-предикторов. Я хочу подобрать модели для переменной ответа, используя каждую из переменных-предикторов отдельно, и, наконец, создать фрейм данных, содержащий коэффициенты модели. Раньше я бы сделал так:

data(iris)

iris_vars <- c("Sepal.Width", "Petal.Length", "Petal.Width")
fits.iris <- lapply(iris_vars, function(x) {lm(substitute(Sepal.Length ~ i, list(i = as.name(x))), data = iris)})

# extract model coeffs, so forth and so on, eventually combining into a result dataframe
iris.p <- as.data.frame(lapply(fits.iris, function(f) summary(f)$coefficients[,4]))
iris.r <- as.data.frame(lapply(fits.iris, function(f) summary(f)$r.squared))

Однако теперь, когда я начал использовать dplyr, broom и т. Д., Это кажется немного громоздким. Используя purrr::map, я могу более или менее воссоздать этот список моделей:

# using purrr, still uses the Response variable "Sepal.Length" as a predictor of itself
iris %>%
select(1:4) %>%
# names(select(., 2:4)) %>% this does not work 
names() %>%
paste('Sepal.Length ~', .) %>% 
map(~lm(as.formula(.x), data = iris))

Однако я не уверен, как преобразовать этот список в соответствующую форму для использования с broom::tidy. Если бы я использовал сгруппированные строки, а не столбцы, я бы сохранил подгонки модели и использовал broom::tidy, чтобы сделать что-то вроде этого:

iris.fits <- group_by(Species) %>% do(modfit1 = lm(Sepal.Length~Sepal.Width,data=.))
tidy(iris.fits, modfit1)

Конечно, я делаю не это, но я надеялся, что при использовании столбцов данных будет аналогичная процедура. Есть ли способ, возможно, использовать purrr::nest или что-то подобное для создания желаемого результата?


r dplyr purrr broom
person nofunsally    schedule 20.12.2016    source источник
comment
Как должен выглядеть окончательный результат? Выход map_df(tidy), добавленный в вашу трубу, близок к тому, что вы хотите?   -  person aosmith    schedule 20.12.2016
comment
да, map_df(tidy) близко. Хотя, map_df(glance) не работает.   -  person nofunsally    schedule 20.12.2016
comment
Вы можете setNames перед использованием glance, если вы определите вектор переменных, используемых в качестве независимых переменных. Например, если вы просто работаете с iris_vars, вы можете сделать что-то вроде ... %>% setNames(., iris_vars) %>% map_df(glance, .id = "variable")   -  person aosmith    schedule 20.12.2016
comment
код setNames приводит к Error in data.frame(r.squared = r.squared, adj.r.squared = adj.r.squared, : object 'fstatistic' not found   -  person nofunsally    schedule 20.12.2016
comment
Эта ошибка связана с glance и тем фактом, что одна из ваших моделей имеет проблемы с использованием одной и той же переменной в левой и правой частях формулы lm. Почувствуйте разницу с map(iris_vars, ~lm(as.formula(paste("Sepal.Length ~", .x)), data = iris)) %>% setNames(., iris_vars) %>% map_df(glance, .id = "variable")   -  person aosmith    schedule 20.12.2016

Ответы (4)


arrow_upward
3
arrow_downward

1) Это дает выходные данные glance и tidy для соответствия модели:

library(broom)

make_model <- function(nm) lm(iris[c("Sepal.Length", nm)])
fits <- Map(make_model, iris_vars)

glance_tidy <- function(x) c(unlist(glance(x)), unlist(tidy(x)[, -1]))
out <- sapply(fits, glance_tidy)

1a) или как магистральный трубопровод:

library(magrittr)
out <- iris_vars %>% Map(f = make_model) %>% sapply(glance_tidy)

Любой из них дает следующую матрицу:

> out
                Sepal.Width   Petal.Length    Petal.Width
r.squared      1.382265e-02   7.599546e-01   6.690277e-01
adj.r.squared  7.159294e-03   7.583327e-01   6.667914e-01
sigma          8.250966e-01   4.070745e-01   4.779948e-01
statistic      2.074427e+00   4.685502e+02   2.991673e+02
p.value        1.518983e-01   1.038667e-47   2.325498e-37
df             2.000000e+00   2.000000e+00   2.000000e+00
logLik        -1.829958e+02  -7.702021e+01  -1.011107e+02
AIC            3.719917e+02   1.600404e+02   2.082215e+02
BIC            3.810236e+02   1.690723e+02   2.172534e+02
deviance       1.007561e+02   2.452503e+01   3.381489e+01
df.residual    1.480000e+02   1.480000e+02   1.480000e+02
estimate1      6.526223e+00   4.306603e+00   4.777629e+00
estimate2     -2.233611e-01   4.089223e-01   8.885803e-01
std.error1     4.788963e-01   7.838896e-02   7.293476e-02
std.error2     1.550809e-01   1.889134e-02   5.137355e-02
statistic1     1.362763e+01   5.493890e+01   6.550552e+01
statistic2    -1.440287e+00   2.164602e+01   1.729645e+01
p.value1       6.469702e-28  2.426713e-100  3.340431e-111
p.value2       1.518983e-01   1.038667e-47   2.325498e-37

или транспонировано:

> t(out)
              r.squared adj.r.squared     sigma  statistic      p.value df
Sepal.Width  0.01382265   0.007159294 0.8250966   2.074427 1.518983e-01  2
Petal.Length 0.75995465   0.758332718 0.4070745 468.550154 1.038667e-47  2
Petal.Width  0.66902769   0.666791387 0.4779948 299.167312 2.325498e-37  2
                 logLik      AIC      BIC  deviance df.residual estimate1
Sepal.Width  -182.99584 371.9917 381.0236 100.75610         148  6.526223
Petal.Length  -77.02021 160.0404 169.0723  24.52503         148  4.306603
Petal.Width  -101.11073 208.2215 217.2534  33.81489         148  4.777629
              estimate2 std.error1 std.error2 statistic1 statistic2
Sepal.Width  -0.2233611 0.47889634 0.15508093   13.62763  -1.440287
Petal.Length  0.4089223 0.07838896 0.01889134   54.93890  21.646019
Petal.Width   0.8885803 0.07293476 0.05137355   65.50552  17.296454
                  p.value1     p.value2
Sepal.Width   6.469702e-28 1.518983e-01
Petal.Length 2.426713e-100 1.038667e-47
Petal.Width  3.340431e-111 2.325498e-37

2) Если мы удалим первый unlist из определения функции glance_tidy, мы получим 2-й список (а не 2-ю числовую матрицу):

glance_tidy_l <- function(x) c(glance(x), unlist(tidy(x)[, -1]))
iris_vars %>% Map(f = make_model) %>% sapply(glance_tidy_l)

              Sepal.Width  Petal.Length  Petal.Width  
r.squared     0.01382265   0.7599546     0.6690277    
adj.r.squared 0.007159294  0.7583327     0.6667914    
sigma         0.8250966    0.4070745     0.4779948    
statistic     2.074427     468.5502      299.1673     
p.value       0.1518983    1.038667e-47  2.325498e-37 
df            2            2             2            
logLik        -182.9958    -77.02021     -101.1107    
AIC           371.9917     160.0404      208.2215     
BIC           381.0236     169.0723      217.2534     
deviance      100.7561     24.52503      33.81489     
df.residual   148          148           148          
estimate1     6.526223     4.306603      4.777629     
estimate2     -0.2233611   0.4089223     0.8885803    
std.error1    0.4788963    0.07838896    0.07293476   
std.error2    0.1550809    0.01889134    0.05137355   
statistic1    13.62763     54.9389       65.50552     
statistic2    -1.440287    21.64602      17.29645     
p.value1      6.469702e-28 2.426713e-100 3.340431e-111
p.value2      0.1518983    1.038667e-47  2.325498e-37
person G. Grothendieck    schedule 20.12.2016

arrow_upward
2
arrow_downward

Мой ответ по духу похож на ответ Джулии Силге и wysiwyg, но я хочу избежать ручного ввода имен переменных и сохранить имена для ответа и предиктора в формуле объекта модели:

require(tibble)
require(dplyr)
require(tidyr)
require(purrr)
require(broom)

df <- iris
response_var <- "Sepal.Length"

vars <- tibble(response=response_var,
               predictor=setdiff(names(df), response_var))

compose_formula <- function(x, y)
  as.formula(paste0("~lm(", y, "~", x, ", data=.)"))

models <- tibble(data=list(df)) %>%
           crossing(vars) %>%
           mutate(fmla = map2(predictor, response, compose_formula),
                  model = map2(data, fmla, ~at_depth(.x, 0, .y)))

models %>% unnest(map(model, tidy))
models %>% unnest(map(model, glance), .drop=T)

Выход:

# A tibble: 9 x 7
      response    predictor              term   estimate  std.error statistic
         <chr>        <chr>             <chr>      <dbl>      <dbl>     <dbl>
1 Sepal.Length  Sepal.Width       (Intercept)  6.5262226 0.47889634 13.627631
2 Sepal.Length  Sepal.Width       Sepal.Width -0.2233611 0.15508093 -1.440287
3 Sepal.Length Petal.Length       (Intercept)  4.3066034 0.07838896 54.938900
4 Sepal.Length Petal.Length      Petal.Length  0.4089223 0.01889134 21.646019
5 Sepal.Length  Petal.Width       (Intercept)  4.7776294 0.07293476 65.505517
6 Sepal.Length  Petal.Width       Petal.Width  0.8885803 0.05137355 17.296454
7 Sepal.Length      Species       (Intercept)  5.0060000 0.07280222 68.761639
8 Sepal.Length      Species Speciesversicolor  0.9300000 0.10295789  9.032819
9 Sepal.Length      Species  Speciesvirginica  1.5820000 0.10295789 15.365506
# ... with 1 more variables: p.value <dbl>

а также

# A tibble: 4 x 13
      response    predictor  r.squared adj.r.squared     sigma  statistic
         <chr>        <chr>      <dbl>         <dbl>     <dbl>      <dbl>
1 Sepal.Length  Sepal.Width 0.01382265   0.007159294 0.8250966   2.074427
2 Sepal.Length Petal.Length 0.75995465   0.758332718 0.4070745 468.550154
3 Sepal.Length  Petal.Width 0.66902769   0.666791387 0.4779948 299.167312
4 Sepal.Length      Species 0.61870573   0.613518054 0.5147894 119.264502
# ... with 7 more variables: p.value <dbl>, df <int>, logLik <dbl>, AIC <dbl>,
#   BIC <dbl>, deviance <dbl>, df.residual <int>
person LmW.    schedule 14.05.2017

arrow_upward
1
arrow_downward

Если вы готовы немного поработать по настройке квази-вложенного фрейма данных со столбцами-списками, чтобы начать работу, шаг map / model / unnest / tidy выпадает довольно хорошо.

Сначала настройте фрейм данных:

> library(dplyr)
> 
> nested_df <- data_frame(data = list(iris %>% 
                                          select(response = Sepal.Length, 
                                                 predictor = Sepal.Width), 
                                      iris %>% 
                                          select(response = Sepal.Length, 
                                                 predictor = Petal.Length),
                                      iris %>% 
                                          select(response = Sepal.Length, 
                                                 predictor = Petal.Width)))
> 
> nested_df
# A tibble: 3 × 1
                    data
                  <list>
1 <data.frame [150 × 2]>
2 <data.frame [150 × 2]>
3 <data.frame [150 × 2]>

Теперь используйте purrr, tidyr и broom, чтобы получить результаты моделирования.

> library(tidyr)
> library(purrr)
> library(broom)
> 
> nested_df %>%
      mutate(models = map(data, ~ lm(response ~ predictor, .))) %>%
      unnest(map(models, tidy))
# A tibble: 6 × 5
         term   estimate  std.error statistic       p.value
        <chr>      <dbl>      <dbl>     <dbl>         <dbl>
1 (Intercept)  6.5262226 0.47889634 13.627631  6.469702e-28
2   predictor -0.2233611 0.15508093 -1.440287  1.518983e-01
3 (Intercept)  4.3066034 0.07838896 54.938900 2.426713e-100
4   predictor  0.4089223 0.01889134 21.646019  1.038667e-47
5 (Intercept)  4.7776294 0.07293476 65.505517 3.340431e-111
6   predictor  0.8885803 0.05137355 17.296454  2.325498e-37

Вы можете использовать filter, чтобы извлечь только наклоны (term == "predictor"), или вы можете использовать glance вместо tidy в последней строке кода, чтобы получить эти результаты.

person Julia Silge    schedule 29.12.2016

arrow_upward
1
arrow_downward

Другой вариант для квазивложенного фрейма данных, который бы быстро работал с большим количеством переменных-предикторов, - это использование purrr :: map_df. Это также может дать вам предиктор в виде столбца во включающем фрейме данных. Затем следуйте примеру Джулии при подгонке моделей.

> library(dplyr)
> library(purrr)
>
> def_nested_df <- function(x) {
    data_frame("covariate" = x,
               "data" = list(iris %>% tbl_df %>%
                               select_("response" = "Sepal.Length", 
                                       "predictor" = x)))
  }    
> 
> nested_df <- 
    c("Sepal.Width", "Petal.Length", "Petal.Width") %>%
    map_df(def_nested_df)
>
> nested_df
# A tibble: 3 × 2
     covariate               data
         <chr>             <list>
1  Sepal.Width <tibble [150 × 2]>
2 Petal.Length <tibble [150 × 2]>
3  Petal.Width <tibble [150 × 2]>
>
> nested_df[[1, "data"]]
# A tibble: 150 × 2
   response predictor
      <dbl>     <dbl>
1       5.1       3.5
2       4.9       3.0
3       4.7       3.2
4       4.6       3.1
5       5.0       3.6
6       5.4       3.9
7       4.6       3.4
8       5.0       3.4
9       4.4       2.9
10      4.9       3.1
# ... with 140 more rows
person Matt Mulvahill    schedule 06.02.2017