Проведите тест Aov через Tibble аккуратно

Question

Я хочу запустить линейную регрессию для фрейма данных, используя ту же самую зависимую переменную.Подобный вопрос был решен здесь .Проблема в том, что функция aov для реализации ANOVA не принимает x и y в качестве аргументов (насколько я знаю).Есть ли способ провести анализ аккуратно?До сих пор я пробовал что-то вроде:

library(tidyverse)

iris %>% 
  as_tibble() %>% 
  select(Sepal.Length, Species) %>% 
  mutate(foo_a = as_factor(sample(c("a", "b", "c"), nrow(.), replace = T)),
         foo_b = as_factor(sample(c("d", "e", "f"), nrow(.), replace = T))) %>% 
  map(~aov(Sepal.Length ~ .x, data = .))

^{Создано в 2019-02-12 с помощью представительного пакета (v0.2.1)}

Желаемый результат - три анализа: Sepal.Length и Species, Sepal.Length и foo_a и последний Sepal.Length и foo_b.Это возможно или я совершенно не прав?

Answer 1

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

Чтобы сэкономить место, я работаю с iris_foo, который является вашими даннымикак вы создали его через 2 строки мутирования.Если поместить его в длинный формат, вы получите ключ к именам этих трех столбцов, которые будут использоваться в качестве независимых переменных в каждом из вызовов aov.

library(tidyverse)

iris_foo %>%
  gather(key, value, -Sepal.Length)

#> # A tibble: 450 x 3
#>    Sepal.Length key     value 
#>           <dbl> <chr>   <chr> 
#>  1          5.1 Species setosa
#>  2          4.9 Species setosa
#>  3          4.7 Species setosa
#>  4          4.6 Species setosa
#>  5          5   Species setosa
#>  6          5.4 Species setosa
#>  7          4.6 Species setosa
#>  8          5   Species setosa
#>  9          4.4 Species setosa
#> 10          4.9 Species setosa
#> # … with 440 more rows

Оттуда, гнездо по keyи создайте новый список-столбец моделей ANOVA.Это будет список aov объектов.Для простоты возврата ваших моделей вы можете удалить столбец данных.

aov_models <- iris_foo %>%
  gather(key, value, -Sepal.Length) %>%
  group_by(key) %>%
  nest() %>%
  mutate(model = map(data, ~aov(Sepal.Length ~ value, data = .))) %>%
  select(-data)

aov_models
#> # A tibble: 3 x 2
#>   key     model    
#>   <chr>   <list>   
#> 1 Species <S3: aov>
#> 2 foo_a   <S3: aov>
#> 3 foo_b   <S3: aov>

Оттуда вы можете работать с моделями так, как вам нравится.Они доступны в списке aov_models$model.Напечатано, они выглядят так, как вы ожидаете.Например, первая модель:

aov_models$model[[1]]
#> Call:
#>    aov(formula = Sepal.Length ~ value, data = .)
#> 
#> Terms:
#>                    value Residuals
#> Sum of Squares  63.21213  38.95620
#> Deg. of Freedom        2       147
#> 
#> Residual standard error: 0.5147894
#> Estimated effects may be unbalanced

Чтобы увидеть все модели, звоните aov_models$model %>% map(print).Вы также можете использовать broom функции, такие как broom::tidy или broom::glance, в зависимости от того, как вам нужно представить модели.

aov_models$model %>%
  map(broom::tidy)
#> [[1]]
#> # A tibble: 2 x 6
#>   term         df sumsq meansq statistic   p.value
#>   <chr>     <dbl> <dbl>  <dbl>     <dbl>     <dbl>
#> 1 value         2  63.2 31.6        119.  1.67e-31
#> 2 Residuals   147  39.0  0.265       NA  NA       
#> 
#> [[2]]
#> # A tibble: 2 x 6
#>   term         df   sumsq meansq statistic p.value
#>   <chr>     <dbl>   <dbl>  <dbl>     <dbl>   <dbl>
#> 1 value         2   0.281  0.141     0.203   0.817
#> 2 Residuals   147 102.     0.693    NA      NA    
#> 
#> [[3]]
#> # A tibble: 2 x 6
#>   term         df   sumsq meansq statistic p.value
#>   <chr>     <dbl>   <dbl>  <dbl>     <dbl>   <dbl>
#> 1 value         2   0.756  0.378     0.548   0.579
#> 2 Residuals   147 101.     0.690    NA      NA

Или приведение всех моделей в один фрейм данных, который содержит столбец key, вы можете сделать:

aov_models %>%
  mutate(model_tidy = map(model, broom::tidy)) %>%
  unnest(model_tidy)