运行 许多方差分析和提取某些列的快速方法
A fast way to run many anova's and extract certain columns
我的数据包含一个响应变量 (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 对于我正在拟合的模型中的每个因素，其中每一行都是 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 的数量而异。
我的数据包含一个响应变量 (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 对于我正在拟合的模型中的每个因素，其中每一行都是 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 的数量而异。