aov throws "could not find function 'Error'" when used in the context of group_by

Question

I'm trying to run a repeated measures ANCOVA. The following code works fine:

tidy(aov(FA ~ sex * study + Error(PartID), data = DTI.TRACTlong))

Where FA is a continuous measure, sex and study are factors where study indicates (time 1 or time 2) and PartID is the individual id. However, I have to run this analysis for a number of regions (ROI) for two different conditions (harmonized vs. not). This seems easy enough using tidyverse with group_by (see below), but it throws Error in Error(.$PartID) : could not find function "Error". Any idea what is happening here? Why is the Error function recognized when used on its own but not when using tidyverse with group_by?

  group_by(harmonize,ROI) %>%
  do(tidy(aov(.$FA ~ .$GOBS_Gender * .$study + Error(.$PartID))))

akrun · Accepted Answer · 2021-08-25T21:48:27.627

Specify the data and remove the .$

library(broom)
library(dplyr)
...
    group_by(harmonize,ROI) %>%
    do(tidy(aov(FA ~ sex * study + Error(PartID), data = .)))

Another option may be to use to nest_by

library(tidyr)
...
  nest_by(harmonize, ROI) %>%
  mutate(out = list(tidy(aov(FA ~ sex * study + Error(PartID), data = data)))) %>%
  select(-data) %>%
  unnest(c(out))

Using a reproducible example

> data(npk)
> npk$grp <- rep(c('a', 'b'), each = 12)

-do method

> npk %>% 
     group_by(grp) %>% 
     do(tidy(aov(yield ~  N*P*K + Error(block), data = .)))
# A tibble: 18 x 8
# Groups:   grp [2]
   grp   stratum term         df   sumsq  meansq statistic  p.value
   <chr> <chr>   <chr>     <dbl>   <dbl>   <dbl>     <dbl>    <dbl>
 1 a     block   N:P:K         1  69.0    69.0      3.12    0.328  
 2 a     block   Residuals     1  22.1    22.1     NA      NA      
 3 a     Within  N             1 119.    119.       6.89    0.0786 
 4 a     Within  P             1   0.270   0.270    0.0156  0.908  
 5 a     Within  K             1  58.1    58.1      3.36    0.164  
 6 a     Within  N:P           1  12.2    12.2      0.705   0.463  
 7 a     Within  N:K           1  23.4    23.4      1.35    0.329  
 8 a     Within  P:K           1  44.6    44.6      2.58    0.207  
 9 a     Within  Residuals     3  51.8    17.3     NA      NA      
10 b     block   N:P:K         1  29.3    29.3      0.431   0.630  
11 b     block   Residuals     1  67.9    67.9     NA      NA      
12 b     Within  N             1  73.0    73.0     57.3     0.00478
13 b     Within  P             1  21.3    21.3     16.7     0.0264 
14 b     Within  K             1  38.2    38.2     29.9     0.0120 
15 b     Within  N:P           1   8.52    8.52     6.68    0.0814 
16 b     Within  N:K           1  31.5    31.5     24.7     0.0156 
17 b     Within  P:K           1  47.3    47.3     37.1     0.00888
18 b     Within  Residuals     3   3.82    1.27    NA      NA

-group_modify @andrew reece

> npk %>% 
     group_by(grp) %>% 
     group_modify(~aov(yield ~  N*P*K + Error(block), data = .) %>% 
     tidy)
# A tibble: 18 x 8
# Groups:   grp [2]
   grp   stratum term         df   sumsq  meansq statistic  p.value
   <chr> <chr>   <chr>     <dbl>   <dbl>   <dbl>     <dbl>    <dbl>
 1 a     block   N:P:K         1  69.0    69.0      3.12    0.328  
 2 a     block   Residuals     1  22.1    22.1     NA      NA      
 3 a     Within  N             1 119.    119.       6.89    0.0786 
 4 a     Within  P             1   0.270   0.270    0.0156  0.908  
 5 a     Within  K             1  58.1    58.1      3.36    0.164  
 6 a     Within  N:P           1  12.2    12.2      0.705   0.463  
 7 a     Within  N:K           1  23.4    23.4      1.35    0.329  
 8 a     Within  P:K           1  44.6    44.6      2.58    0.207  
 9 a     Within  Residuals     3  51.8    17.3     NA      NA      
10 b     block   N:P:K         1  29.3    29.3      0.431   0.630  
11 b     block   Residuals     1  67.9    67.9     NA      NA      
12 b     Within  N             1  73.0    73.0     57.3     0.00478
13 b     Within  P             1  21.3    21.3     16.7     0.0264 
14 b     Within  K             1  38.2    38.2     29.9     0.0120 
15 b     Within  N:P           1   8.52    8.52     6.68    0.0814 
16 b     Within  N:K           1  31.5    31.5     24.7     0.0156 
17 b     Within  P:K           1  47.3    47.3     37.1     0.00888
18 b     Within  Residuals     3   3.82    1.27    NA      NA

-nest_by method

> npk %>%
     nest_by(grp) %>% 
     mutate(out = list(tidy(aov(yield ~  N*P*K + Error(block), 
            data = data)))) %>%
     select(-data) %>% 
     unnest(out)
# A tibble: 18 x 8
# Groups:   grp [2]
   grp   stratum term         df   sumsq  meansq statistic  p.value
   <chr> <chr>   <chr>     <dbl>   <dbl>   <dbl>     <dbl>    <dbl>
 1 a     block   N:P:K         1  69.0    69.0      3.12    0.328  
 2 a     block   Residuals     1  22.1    22.1     NA      NA      
 3 a     Within  N             1 119.    119.       6.89    0.0786 
 4 a     Within  P             1   0.270   0.270    0.0156  0.908  
 5 a     Within  K             1  58.1    58.1      3.36    0.164  
 6 a     Within  N:P           1  12.2    12.2      0.705   0.463  
 7 a     Within  N:K           1  23.4    23.4      1.35    0.329  
 8 a     Within  P:K           1  44.6    44.6      2.58    0.207  
 9 a     Within  Residuals     3  51.8    17.3     NA      NA      
10 b     block   N:P:K         1  29.3    29.3      0.431   0.630  
11 b     block   Residuals     1  67.9    67.9     NA      NA      
12 b     Within  N             1  73.0    73.0     57.3     0.00478
13 b     Within  P             1  21.3    21.3     16.7     0.0264 
14 b     Within  K             1  38.2    38.2     29.9     0.0120 
15 b     Within  N:P           1   8.52    8.52     6.68    0.0814 
16 b     Within  N:K           1  31.5    31.5     24.7     0.0156 
17 b     Within  P:K           1  47.3    47.3     37.1     0.00888
18 b     Within  Residuals     3   3.82    1.27    NA      NA

@AmandaRodrigue thanks. I couldn't test it properly as there was no example data in the post — akrun, Aug 25 '21 at 20:32
Can also just do `... %>% aov(Petal.Width ~ grp + Error(PartID), data = .) %>% tidy`, no need for `do` — andrew_reece, Aug 25 '21 at 20:34
@andrew_reece I think that will drop the groups if I am not wrong. Can you check `npk %>% group_by(grp) %>% aov(yield ~ N*P*K + Error(block), data = .) %>% tidy` (based on the modified example I created) — akrun, Aug 25 '21 at 20:43
@akrun you are correct. i think we can still avoid `do` with `group_modify`: `group_modify(~aov(yield ~ N*P*K + Error(block), data = .) %>% tidy)` — andrew_reece, Aug 25 '21 at 21:34
@andrew_reece `group_map` is kind of experimental life cycle, so the solution may not last if they decide to deprecate in next release :-). Probably, group_split/map — akrun, Aug 25 '21 at 21:36
Good point @akrun (also I realized it makes more sense to use `group_modify` and avoid the `bind_rows`. Although it's noteworthy that the [`group_modify` docs](https://dplyr.tidyverse.org/reference/group_map.html) explicitly recommend, "group_modify() is an evolution of do()". — andrew_reece, Aug 25 '21 at 21:42

aov throws "could not find function 'Error'" when used in the context of group_by

1 Answers1