Подгонка моделей к группированным данным с помощью условных операторов

Question

Чтобы упростить воспроизведение, я использую набор данных goats из пакета ResourceSelection , который содержит пространственные данные для используемых (STATUS == 1) и «доступных» (STATUS == 0) GPS-местоположений горыкозы.ID для отдельных (n = 10) и ELEVATION, ..., TASP - атрибуты точек.

library(ResourceSelection)
head(goats)
  STATUS ID ELEVATION   SLOPE       ET   ASPECT       HLI      TASP
1      1  1       651 38.5216  35.3553 243.1131 0.9175926 0.9468804
2      1  1       660 39.6927  70.7107 270.0000 0.8840338 0.6986293
3      1  1       316 20.5477  50.0000 279.2110 0.7131423 0.5749115
4      1  1       334 34.0783  35.3553 266.1859 0.8643775 0.7447368
5      1  1       454 41.6187  25.0000 258.3106 0.9349181 0.8292587
6      1  1       343 28.4694 103.0776 237.0426 0.8254866 0.9756112

Я бы хотел добавить glm к каждому человеку, сгруппированному по Season (создан ниже), но структура модели варьируется в зависимости от сезона.Я использовал несколько полезных SO сообщений и других ресурсов , но все они соответствуют одной модели для каждой группы, где, как мне хотелось бы, для каждой группы Season подходят разные модели..

#Add a new `Season` field 
library(tidyverse)

goats <- goats %>% 
  mutate(Season = if_else(ID %in% 1:3, "Summer",
                          if_else(ID %in% 4:7, "Winter", "Fall")))

Ниже я создаю функцию построения модели и определяю модель, специфичную для каждого сезона, используя if else.

SeasonalMods <- function(df) {
  #Models for Archery
  if(Season == "Summer") {
    glm(STATUS ~ SLOPE + I(SLOPE^2), data = df)
  #Models for Winter  
  } else if (Season == "Winter") {
    glm(STATUS ~ SLOPE + ASPECT + TASP, data = df)
  #Models for Fall
   } else if (Season == "Fall") {
    glm(STATUS ~ ELEVATION + SLOPE + ET + ASPECT + HLI + TASP, data = df)}
  }

. Затем я пытаюсь сопоставить функцию с сгруппированными данными.и создать новый список-столбцы следующим образом.

ModelFits <- goats %>%
  group_by(Season, ID) %>% 
  nest() %>% 
  mutate(fits = map(data, SeasonalMods),
         tidied = map(fits, tidy),
         glanced = map(fits, glance),
         augmented = map(fits, augment))

, которая приводит к следующей ошибке:

Error in mutate_impl(.data, dots) : 
  Evaluation error: object 'Season' not found

Я не уверен, как правильно указать Season в функции SeasonalMods, чтобы ее можно было интерпретировать как map().

Я пытался добавить df$ перед Season в операторах if и else if, но это также приводит к ошибке.

Answer 1

После того, как я предложил использовать modelr, я пошел дальше без него.Как я уже говорил выше, ваша функция SeasonalMods не знает, что Season - это столбец во фрейме данных, который он принимает в качестве аргумента, поэтому вы получите ошибки, которые не определены.Один из способов сделать это вместо этого - добавить второй аргумент в функцию, чтобы взять сезон.Поскольку вы вкладываете свои данные, теперь легко передать данные и сезон в функцию моделирования.Я использую map2, потому что столбцы data и Season имеют одинаковую длину.

library(tidyverse)
library(broom)

Все внутренние компоненты этой функции одинаковы - я только добавил второй аргумент.

SeasonalMods <- function(df, Season) {
  ...
}

Просто чтобы проиллюстрировать используемые вами функции broom, я добавил столбец tidied и сохранил этот фрейм данных:

models <- goats %>%
  group_by(Season, ID) %>%
  nest() %>%
  mutate(fits = map2(data, Season, ~SeasonalMods(.x, .y))) %>%
  mutate(tidied = map(fits, tidy))

head(models)
#> # A tibble: 6 x 5
#>   Season    ID data                 fits      tidied          
#>   <chr>  <int> <list>               <list>    <list>          
#> 1 Summer     1 <tibble [2,106 × 7]> <S3: glm> <tibble [3 × 5]>
#> 2 Summer     2 <tibble [1,668 × 7]> <S3: glm> <tibble [3 × 5]>
#> 3 Summer     3 <tibble [1,539 × 7]> <S3: glm> <tibble [3 × 5]>
#> 4 Winter     4 <tibble [951 × 7]>   <S3: glm> <tibble [4 × 5]>
#> 5 Winter     5 <tibble [1,908 × 7]> <S3: glm> <tibble [4 × 5]>
#> 6 Winter     6 <tibble [2,184 × 7]> <S3: glm> <tibble [4 × 5]>

Просто чтобы проверить, чтомодели получили разные формулы для разных сезонов:

models$fits[[1]]
#> 
#> Call:  glm(formula = STATUS ~ SLOPE + I(SLOPE^2), data = df)
#> 
#> Coefficients:
#> (Intercept)        SLOPE   I(SLOPE^2)  
#>   -0.042618    -0.000989     0.000375  
#> 
#> Degrees of Freedom: 2105 Total (i.e. Null);  2103 Residual
#> Null Deviance:       468 
#> Residual Deviance: 337.2     AIC: 2127

models$fits[[6]]
#> 
#> Call:  glm(formula = STATUS ~ SLOPE + ASPECT + TASP, data = df)
#> 
#> Coefficients:
#> (Intercept)        SLOPE       ASPECT         TASP  
#>    0.024625     0.017838    -0.001768     0.215217  
#> 
#> Degrees of Freedom: 2183 Total (i.e. Null);  2180 Residual
#> Null Deviance:       485.3 
#> Residual Deviance: 385.7     AIC: 2421