我有数据由一个响应变量(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中，该具有我正在拟合的模型中的每个因子的anovap-value列，其中每一行都是我正在运行anova的每个group。
这是我目前正在做的事情：
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,000group，所以我想知道是否有任何更快的方法来做到这一点。
我认为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
}
我不确定你为什么使用"哪个"来选择"Pr(>F("列，因为只能有一个，所以直接子集。另请注意组的基本子集，-nrow(an)删除最后一行内容。
我还尽可能多地留在循环之外，因此转换为数据框和添加组 ID 都在循环之外。rbindlapply 返回一个矩阵，使用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
真正的加速可以通过并行处理来实现，因为每个分组的方差值都是独立的 - 您有多少个 CPU 内核？在我的代码中使用parallel:mclapply将运行时间仅缩短到 4.4 秒，但您的改进可能会因数据大小和 CPU 数量而异。