今日の授業用のファイルです。
イネ‘日本晴’の収穫量のデータ(CSV形式、因子を楽に扱える)
処理区ごとの繰り返し数が異なる場合は線形モデルを使った分散分析を行う。(繰り返し数が同じでも問題ない)これを行うのがaov関数で、内部的にlm関数を使って線形モデルの計算を行っている。イネのデータ分析手順は以下の通り。上と同様に収量(gy)について、因子Aを栽植密度(dens、水準数a=2)、Bを肥料(fert、b=3)、Cを年次効果(year、c=4)とした場合の分析法を示す。ただし、ブロックを繰り返し(r=2)として扱う。
作業ディレクトリの設定後、
data <- read.table("rice2.csv", sep=",", header=T) res <- aov(gy~dens+fert+year+dens*fert+fert*year+dens*year+dens*fert*year,data) #gyをdens, fert, year, dens*fert (densとfertの交互作用)...に分解する。 # res <- aov(gy~dens+fert+year,data) このようにすると全ての交互作用はないものとし、dens, fert, yearの3因子の変動のみを考え、残りは全て誤差にプールされる。 summary(res) Df Sum Sq Mean Sq F value Pr(>F) dens 1 124 124 0.033 0.85802 fert 2 67421 33710 8.855 0.00132 ** year 3 54399 18133 4.763 0.00961 ** dens:fert 2 2354 1177 0.309 0.73691 fert:year 6 75080 12513 3.287 0.01667 * dens:year 3 2381 794 0.209 0.88949 dens:fert:year 6 26508 4418 1.161 0.35945 Residuals 24 91362 3807 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
この出力のうち、Dfが自由度、Sum Sqが偏差平方和、Mean Sqが平均平方、F valueが分散比、Prが分散比に対する有意確率、Residualsが誤差。aovではformulaに含まれない因子や交互作用はすべて誤差に含まれる。もし、どれかの因子に有意差があった場合は、TukeyのHSD法による多重比較を行い、どの水準間に有意差があるのか確認する。
TukeyHSD(res,"調べたい因子の名前") # 例えばTukeyHSD(res,"fert") Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = gy ~ dens + fert + year, data = data) $fert diff lwr upr p adj f2-f1 85.95000 26.25310 145.64690 0.0031837 f3-f1 70.90625 11.20935 130.60315 0.0166076 f3-f2 -15.04375 -74.74065 44.65315 0.8139811
この場合、肥料の水準f1とf2の間に1%水準で、f1とf3の間に5%水準で有意差がある。
繰り返しのある三元配置データで因子Aがa水準、Bがb水準、Cがc水準、繰り返しがr回の場合、一般にモデルを次のように構造化する。
Xijkl = μ+αi+βj+γk+(αβ)ij+(βγ)jk+(αγ)ik+(αβγ)ijk+εijkl
変動は次のように分解できる。
SST = SSA+SSB+SSC+SSA×B+SSB×C+SSA×C+SSA×B×C+SSE
SSA×B = SSAB-SSA-SSB
ただし、SSABは因子AとBの各水準の組み合わせからなるa×b個の新たな水準(副次級)の変動で、各組み合わせごとの平均と全体の平均の差の平方×その組み合わせのデータ数の合計。Rでは下の手順のように計算できる。
SSA×B×C = SSABC-SSA-SSB-SSC-SSA×B-SSB×C-SSA×C
ただし、SSABCは因子A,B,Cの各水準の組み合わせからなるa×b×c個の新たな水準(副次級)の変動で、各組み合わせごとの平均と全体の平均の差の平方×繰り返し数(r)の合計。
収量(gy)について、因子Aを栽植密度(dens、水準数a=2)、Bを肥料(fert、b=3)、Cを年次効果(year、c=4)とした場合の分析法を示す。ただし、ブロックを繰り返し(r=2)として扱う。
作業ディレクトリの設定後、
rice <- read.table("rice2.csv", sep=",", header=T) m <- mean(rice$gy) ma <- by(rice$gy, rice$dens, mean) mb <- by(rice$gy, rice$fert, mean) mc <- by(rice$gy, rice$year, mean) ab <- interaction(rice$dens, rice$fert) # 因子AとBの組み合わせによる新たな水準 mab <- by(rice$gy, ab, mean) # 因子AとBの組み合わせごとの平均 bc <- interaction(rice$fert, rice$year) mbc <- by(rice$gy, bc, mean) ac <- interaction(rice$dens, rice$year) mac <- by(rice$gy, ac, mean) abc <- interaction(rice$dens, rice$fert, rice$year) mabc <- by(rice$gy, abc, mean) sst <- sum((rice$gy-m)^2) # 総変動 ssa <- sum((ma[rice$dens]-m)^2) # 因子A(密度)の級間変動 ssb <- sum((mb[rice$fert]-m)^2) # 因子B(肥料)の級間変動 ssc <- sum((mc[rice$year]-m)^2) # 因子C(年次)の級間変動 ssab <- sum((mab[ab]-m)^2) # 因子AとBの組み合わせの級間変動 ssaxb <- ssab-ssa-ssb # 因子AとBの交互作用による変動 ssbc <- sum((mbc[bc]-m)^2) # 因子BとCの組み合わせの級間変動 ssbxc <- ssbc-ssb-ssc # 因子BとCの交互作用による変動 ssac <- sum((mac[ac]-m)^2) # 因子AとCの組み合わせの級間変動 ssaxc <- ssac-ssa-ssc # 因子AとCの交互作用による変動 ssabc <- sum((mabc[abc]-m)^2) # 因子A, B, Cの組み合わせの級間変動 ssaxbxc <- ssabc-ssa-ssb-ssc-ssaxb-ssbxc-ssaxc # 因子A, B, Cの交互作用による変動 sse <- sst-ssa-ssb-ssc-ssaxb-ssbxc-ssaxc-ssaxbxc # 誤差変動 # 実は、SSABC = SSA+SSB+SSC+SSA×B+SSB×C+SSA×C+SSA×B×C なので、SSE = SST - SSABC のほうが早い。 pht <- 2*3*2*4-1 # 総変動の自由度、全データ数-1 = abcr-1 pha <- 2-1 # 因子A(密度)の級間変動の自由度、水準数-1 phb <- 3-1 # 因子B(肥料)の級間変動の自由度、水準数-1 phc <- 4-1 # 因子C(年次)の級間変動の自由度、水準数-1 phaxb <- (2-1)*(3-1) # 因子AとBの交互作用の自由度 phbxc <- (3-1)*(4-1) # 因子BとCの交互作用の自由度 phaxc <- (2-1)*(4-1) # 因子AとCの交互作用の自由度 phaxbxc <- (2-1)*(3-1)*(4-1) # 因子A, B, Cの交互作用の自由度 phe <- pht-pha-phb-phaxb-phbxc-phaxc-phaxbxc # 誤差変動の自由度 msa <- ssa/pha # 因子A(密度)の級間の平均平方 以下、略 f0a <- msa/mse # 因子A(密度)の分散比(F値) levels <- c(0.05, 0.01, 0.001) qf(levels, pha, phe, lower.tail=F) # 因子A(密度)の5, 1, 0.1%水準の棄却域、F値がこれより大きければ因子Aについての帰無仮説は棄却される。 # 得られたF値に対する有意確率は次のように求められる。 pf(f0a, pha, phe, lower.tail=F)
変動因 | 自由度 | 偏差平方和 | 平均平方 | 分散比 |
---|---|---|---|---|
密度 | (a-1)=1 | 124.4852 | 124.4852 | |
肥料 | (b-1)=2 | 67420.87 | 33710.44 | |
年次 | (c-1)=3 | |||
密度×肥料 | (a-1)(b-1)=2 | 2354.103 | 1177.051 | |
肥料×年次 | (b-1)(c-1)=6 | |||
密度×年次 | (a-1)(c-1)=3 | |||
密度×肥料×年次 | (a-1)(b-1)(c-1)=6 | |||
誤差 | 24 |
* p < 0.05
2016年11月22日作成、2016年11月22日更新
國分 尚 (Hisashi Kokubun)