anova や glm を、for 関数で処理する方法
色々な処理区があり、各処理区の群集を特徴付ける指標がたくさん得られてきた。
そして、その処理区ごとの各指標の平均値の違いは、処理区ごとに説明できるのかを統計的に示す時に、分散分析を行うと思います。
Rプログラムの中では、anova 関数や aov 関数が分散分析を行う関数になり、この関数は、formula 型の Input ( y ~ x, data=df )を要求してきます。
この aov 関数と for 関数を使って、いろんな指標の分散分析を一括で処理する時、formula 関数が便利です。
次のような data frame (省略版) があるとします。
A~Gの処理区があって、3つの指標を測ったマトリックスになっています。各処理区にはに100個の値（１００サンプル）が、ほんとはあります。
treatment Richness Abundance Interaction 1 A 63.03062 99.183025 61.779550401 2 A 59.28171 92.629037 62.014830418 3 A 60.02215 99.844330 59.648745452 4 B 77.90848 46.798696 54.998186842 5 B 71.64935 41.451237 61.789154721 6 B 71.02360 52.430302 56.806729181 7 C 32.37181 79.320784 44.909000587 8 C 40.65879 65.529647 55.526555285 9 C 31.65834 65.384726 44.569884241 10 D 93.80565 62.985476 86.305609523 11 D 73.20372 62.672580 83.021430326 12 D 76.77968 69.194962 81.886503028 13 E 56.08483 -12.653348 4.713663180 14 E 61.87220 -5.530885 0.005947824 15 E 50.68268 5.953438 -2.610380320 16 F 72.80016 3.270578 57.132157097 17 F 76.54925 1.466419 85.587994020 18 F 70.88253 3.530791 79.221027188 19 G 97.02695 12.273005 7.767906549 20 G 86.16959 -1.815169 0.361944723 21 G 77.37816 26.637699 8.648186671
図にすると以下の通りです。
それでは、各指標、１列ごとにfor 関数を回して、分散分析を行います。
この時、paste 関数で式の形をした文字列を作り、それを formula 関数で type : formula にしてあげます。
result <- list()
for(i in colnames(summ)[-1]){
formula.tmp <- paste(i, 'treatment', sep='~') # 中身としては"Interaction~treatment"という形
formula.object <- formula(formula.tmp) # object の型をformulaにする
result[[i]] <- aov(formula.object, data=summ) # aov 関数にかける
}これで、result の中に各指標の結果を放り込むことができました。
おまけ
分散分析で有意差が出た後、Post-hoc検定、処理区間のペアワイズで何らかの検定を行うそうです。
処理区ごとに、有意差があると見せる時はgeomsignif package が便利そうです。
Rで解析：ggplot2のプロットに有意差バーを追加！「ggsignifr」パッケージ
しかし、処理区間で有意差が出まくると、線がたくさん引かれて微妙な気持ちになったので、別の見せ方のコードを提示しておきます。
線ではなく、有意差があった処理区の名前を表示するようにして見ました。
ここでのポイントは、geom_textの高さをvjust で揃えているところです。
実データではきれいにできたのですが、テストデータだと微妙でしたね...
＊TukeyHSDで２群間比較を行なっていますが、どの方法でPost-hoc検定を行うかは使うデータによって判断してください
summ <- c()
for(l in LETTERS[1:7]) {
df <- matrix(NA, nrow=3, ncol=3, dimnames=list(NULL, c('Richness', 'Abundance', 'Interaction') ))
for(i in 1:ncol(df)){
df[,i] <- rnorm( n= 3, mean=sample(1:100, 1), sd=sample(1:10,1))
}
summ <- rbind(summ, cbind(treatment=l, as.data.frame(df)))
}
library(tidyr)
library(ggplot2)
library(RColorBrewer)
lf.sum <- gather(summ, key=index, value, -1)
head(lf.sum)
ggplot(lf.sum) +
geom_boxplot(aes(x= treatment, y=value)) +
facet_wrap(~key)
result <- list()
for(i in colnames(summ)[-1]){
p.mat <- matrix(NA, ncol=7, nrow=7, dimnames=list(unique(summ$treatment), unique(summ$treatment)))
formula.tmp <- paste(i, 'treatment', sep='~')
formula.object <- formula(formula.tmp)
result[[i]] <- aov(formula.object, data=summ)
if( summary(result[[i]])[[1]][[5]][1] < 0.05 ){
test <- TukeyHSD(aov(formula.object, data=summ))$treat
for(l in 1:nrow(test)){
a <- strsplit(rownames(test), '-')[[l]]
p.mat[a, a] <- test[l,4]
}
diag(p.mat) <- NA
}
for(l in 1:nrow(p.mat)){
sig <- names(which(p.mat[l,] < 0.05))
lf.sum[which( lf.sum$treatment==rownames(p.mat)[l] & lf.sum$index==i ), 'pval'] <- paste(sig, collapse='\n')
lf.sum[which( lf.sum$treatment==rownames(p.mat)[l] & lf.sum$index==i ), 'y'] <-
quantile(lf.sum[which( lf.sum$index==i ), 'value'], probs=seq(0,1,0.1))[11]
}
}
gg <- ggplot(lf.sum)+
geom_boxplot(outlier.size = 0,aes(x=treatment, y=value, color=treatment, fill=after_scale(alpha(color, 0.5) )), width=0.5)+
geom_text( aes(x=treatment, y=y, label=pval),vjust="inward" )+
facet_wrap(~index, scales='free', ncol=2) +
scale_color_manual(values=brewer.pal(8,'Accent')[-c(6)])+
scale_shape_manual(values=c(21:25,6:8))+
theme_minimal(base_size=18) + theme(panel.border=element_rect(fill=NA, size=1),)+
scale_x_discrete(guide = guide_axis(n.dodge = 2))
