之前的博文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,后续的程序才能运行。