之前的博文for循环与plyr包执行方差分析及其多重比较的批量分析里提到使用两种方法进行aov批量分析,但局限是方差分析模型是固定的,那么有没有办法将其改为通用型函数。
答案是当然可以!
1 首先改写通用函数
// aov batch in R
aov.batch <- function(df,mod,nf=NULL,alpha=.05) {
library(agricolae)
res<-list()
fit<-aov(mod,data=df)
aov.res<-summary(fit) # list
res1<-list()
if(is.null(nf)) {
tt<-as.character(mod[[3]])
fn<-tt[tt!="+"]
nf<-length(fn)
} else{
fn<-names(df)[nf]
nf<-length(fn)
}
for(j in 1:nf){
res1[[j]]<-duncan.test(fit,fn[j],alpha=.05)$groups
}
names(res1)<-fn
res<-list(aov.res,res1)
names(res)<-c('aov','comp1')
res
}
aovr.print<-function(res,idx){
print(res[idx])
}
请注意:暂时不能处理因子互作的多重比较。
2 运行示例
// 指定具体的方差分析模型
aov1.batch<-function(df,mod=y~Clone+Treat){
res<-aov.batch(df,mod=mod)
res
}
# 数据重构
dat1 <- tidyr::gather(dat,key=Trait,y,-c(1:3))
aov.res <- plyr::dlply(dat1,"Trait",aov1.batch)
示例数据:
> head(dat[,1:6])
Clone Treat Rep EC RWC RWD
1 1 1 1 17.25 88.43537 11.564626
2 2 1 1 16.88 92.36111 7.638889
3 3 1 1 20.05 81.14754 18.852459
4 4 1 1 23.14 62.90323 37.096774
5 5 1 1 24.06 85.32110 14.678899
6 6 1 1 15.98 65.00000 35.000000
3 查看结果
> idx<-names(aov.res)
> aovr.print(aov.res,idx)
All aov results are as following:
EC :
$`aov`
Df Sum Sq Mean Sq F value Pr(>F)
Clone 8 5269 659 55.45 <2e-16 ***
Treat 3 17594 5865 493.75 <2e-16 ***
Residuals 96 1140 12
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
$comp1
$comp1$`Clone`
y groups
5 48.6575 a
4 44.0475 b
9 40.9050 c
3 36.9675 d
1 35.2275 d
2 32.1100 e
7 29.6600 ef
8 29.1475 fg
6 26.7125 g
$comp1$Treat
y groups
4 53.53111 a
3 41.11333 b
2 29.84556 c
1 19.25889 d
...
到此,演示结束。
通用函数aov.batch编写后,不能直接使用,需要再次生产另一个专用函数aov1.batch,后续的程序才能运行。