Быстрый способ запустить много анова и извлечь определенные столбцы - PullRequest
Новая
       34

Быстрый способ запустить много анова и извлечь определенные столбцы

0 голосов
/ 21 сентября 2018

У меня есть данные, которые состоят из переменной ответа (y) и двух факторов (sex и time), для нескольких group с: 

set.seed(1)
df <- data.frame(y = rnorm(26*18),
                 group = sort(rep(LETTERS,18)),
                 sex = rep(c(rep("F",9),rep("M",9)),26),
                 time = rep(rep(sort(rep(1:3,3)),2),26))
df$sex <- factor(df$sex, levels = c("M","F"))

Я бы хотелпроверить между моделями, используя R anova, для каждого group и объединить все в один data.frame, который имеет столбец anova p-value для каждого из факторов в модели IЯ подхожу, где каждая строка - это каждый из group s, на котором я запускаю anova.

Вот что я сейчас делаю: 

anova.df <- do.call(rbind,lapply(unique(df$group),function(i){
  an.df <- anova(lm(y ~ sex*time,data=df %>% dplyr::filter(group == i)))
  an.df <- data.frame(factor.name=rownames(an.df)[1:(nrow(an.df)-1)],p.value=an.df[1:(nrow(an.df)-1),which(colnames(an.df) == "Pr(>F)")]) %>%
    tidyr::spread(factor.name,p.value) %>%
    dplyr::mutate(group=i)
  return(an.df)
}))

Но вна самом деле у меня ~ 15 000 group с, поэтому мне интересно, есть ли более быстрый способ сделать это.

Ответы [ 2 ]

0 голосов
/ 21 сентября 2018

Вот более оригинальная версия вашего оригинала: 

br <- function(){
    andf = do.call(rbind,lapply(unique(df$group), function(g){
        an = anova(lm(y~sex*time, data=df[df$group==g,]))
        setNames(an[-nrow(an),"Pr(>F)"],rownames(an)[-nrow(an)])
    }))

    andf = data.frame(andf)
    andf$group = unique(df$group)
    andf        
}

Я не уверен, почему вы использовали «which», чтобы выбрать столбец «Pr (> F)», потому что может быть только один,так подмножество это прямо.Также обратите внимание на базовое подмножество для групп и -nrow(an) для удаления последних строк вещей.

Я также оставил как можно больше вне цикла, поэтому преобразование в фрейм данных и добавление идентификатора группывне петли.rbind in lapply возвращает матрицу, а rbind.data.frame медленнее, поэтому мне нужно преобразовать матрицу вне цикла.

Ваш код переупорядочивает столбцы: 

> head(op())
        sex    sex:time      time group
1 0.8396437 0.283887315 0.4983305     A
2 0.4137317 0.626673282 0.5004230     B
3 0.6422066 0.469439754 0.4297816     C

но мой сохраняет порядок от anova: 

> head(br())
        sex      time    sex.time group
1 0.8396437 0.4983305 0.283887315     A
2 0.4137317 0.5004230 0.626673282     B
3 0.6422066 0.4297816 0.469439754     C

вы не говорите, что порядок столбцов важен для вас или нет.

Скорость: сравнение вашего кода с моим с помощью jyjek'sкод: 

> benchmark(op(), jy(), br())
  test replications elapsed relative user.self sys.self user.child sys.child
3 br()          100   4.737    1.000     4.732    0.004          0         0
2 jy()          100   5.368    1.133     5.363    0.004          0         0
1 op()          100  12.769    2.696    12.767    0.000          0         0

Реальная скорость может быть достигнута при параллельной обработке, поскольку каждая сгруппированная анова независима - сколько у вас процессорных ядер?Использование parallel:mclapply в моем коде сократило затраченное время только до 4,4 с, но ваше улучшение может варьироваться в зависимости от размера ваших данных и количества процессоров.

0 голосов
/ 21 сентября 2018

Я думаю, purrr может вам помочь.
Возможно, это не лучшее решение, но попробуйте что-то вроде этого: 

 df%>%
   group_by(group)%>%
   nest()%>%
   mutate(fit = map(data, ~ anova(lm(y ~ sex*time, data = .x))))%>%
   select(group,data,fit)%>%
   unnest(fit)%>%
   select(group,`Pr(>F)`)%>%
   na.omit()%>%
   mutate(var=rep(c("sex","time","sex:time"),times=nrow(.)/3))%>%
   spread(var,`Pr(>F)`)
# A tibble: 26 x 4
   group   sex `sex:time`  time
   <fct> <dbl>      <dbl> <dbl>
 1 A     0.840    0.284   0.498
 2 B     0.414    0.627   0.500
 3 C     0.642    0.469   0.430
 4 D     0.423    0.569   0.567
 5 E     0.169    0.904   0.625
 6 F     0.845    0.00390 0.869
 7 G     0.937    0.318   0.473
 8 H     0.329    0.663   0.609
 9 I     0.977    0.144   0.158
10 J     0.823    0.448   0.193
# ... with 16 more rows

microbenchmark::microbenchmark(x= df%>%
                                  group_by(group)%>%
                                  nest()%>%
                                  mutate(fit = map(data, ~ anova(lm(y ~ sex*time, data = .x))))%>%
                                  select(group,data,fit)%>%
                                  unnest(fit)%>%
                                  select(group,`Pr(>F)`)%>%
                                  na.omit()%>%
                                  mutate(var=rep(c("sex","time","sex:time"),times=nrow(.)/3))%>%
                                  spread(var,`Pr(>F)`),
                                y=anova.df <- do.call(rbind,lapply(unique(df$group),function(i){
                                  an.df <- anova(lm(y ~ sex*time,data=df %>% dplyr::filter(group == i)))
                                  an.df <- data.frame(factor.name=rownames(an.df)[1:(nrow(an.df)-1)],p.value=an.df[1:(nrow(an.df)-1),which(colnames(an.df) == "Pr(>F)")]) %>%
                                    tidyr::spread(factor.name,p.value) %>%
                                    dplyr::mutate(group=i)
                                  return(an.df)
                                })),times=50)
Unit: milliseconds
 expr       min        lq     mean    median        uq      max neval cld
    x  69.98061  71.02417  74.0585  72.45625  74.08786  89.4715    50  a 
    y 166.63844 168.22296 181.6709 171.08077 184.14635 434.8872    50   b
...