以整洁的方式通过 tibble 运行 aov 测试答案

【问题标题】：Run a aov test through a tibble in a tidy way以整洁的方式通过 tibble 运行 aov 测试
【发布时间】：2019-02-12 18:55:16
【问题描述】：

我想使用相同的因变量对数据框进行线性回归。 here 解决了一个类似的问题。问题是实现方差分析的aov 函数不接受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 = .))

^{由reprex package (v0.2.1) 于 2019 年 2 月 12 日创建}

所需的输出是三个分析：Sepal.Length 和 Species、Sepal.Length 和 foo_a 以及最后一个 Sepal.Length 和 foo_b。有可能还是我完全错了？

【问题讨论】：

看看咕噜声。现在无法测试。
您可以重塑为长数据，嵌套在组合的Species/foo_a/foo_b 列上，并使用many models 方法
@camille 你能解释一下你的答案吗？我正在考虑使用嵌套数据，但找不到正确的方法...
@TitoSanz。如果您对以下解决方案不感兴趣，请告诉我。我可以删除它
@akrun 我对您的解决方案感兴趣，但提出的解决方案暗示我将所有因子变量作为列表提供给map2 函数。我的真实数据集有 25 个预测变量数字和因子，所以我试图弄清楚如何给出类似 is.factor 的东西，或者只选择因子变量加上依赖变量并将所有变量作为独立变量传递......

标签： r dplyr

【解决方案1】：

一种方法是将其制成长形数据框，按感兴趣的自变量分组，并使用"many models" 方法。我通常更喜欢这样的事情，而不是尝试跨多个列进行 tidyeval - 它只是让我更清楚地了解正在发生的事情。

为了节省空间，我正在使用iris_foo，这是您通过 2 条 mutate 行创建的数据。将其转换为长格式会为您提供这三列名称的键，这些名称将在每个 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)

【讨论】：

这完全有道理！谢谢