运行 许多方差分析和提取某些列的快速方法

A fast way to run many anova's and extract certain columns

我的数据包含一个响应变量 (y) 和两个因子 (sextime),有几个 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"))

我想在使用 Ranova 的模型之间进行测试,对于每个 group,并将它们全部合并到一个 data.frame 中,其中有一列anova p-value 对于我正在拟合的模型中的每个因素,其中每一行都是 group 中的每一个我是 运行 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 groups 所以我想知道有没有更快的方法来做到这一点。

我觉得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

这是您原件的更简洁版本:

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"到select "Pr(>F)" 列,因为只能有一个,所以直接子集。还要注意组的基本子集化,以及 -nrow(an) 删除最后一行的东西。

我也尽可能多地离开了循环,所以转换为数据框和添加组 ID 都在循环之外。 rbind in lapply returns 一个矩阵,使用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 的代码将你的代码与我的代码进行比较:

> 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

通过并行处理可以实现真正的加速,因为每个分组的 anova 都是独立的 - 您有多少 CPU 个内核?在我的代码中使用 parallel:mclapply 得到的运行时间仅下降到 4.4 秒,但您的改进可能会因数据大小和 CPU 的数量而异。