生物統計学

11月18日の授業について

今日の授業用のファイルです。

イネ‘日本晴’の収穫量のデータ(エクセル形式)

イネ‘日本晴’の収穫量のデータ(CSV形式)

イネ‘日本晴’の収穫量のデータ(CSV形式、因子を楽に扱える)

‘日本晴’のデータで繰り返しのある三元配置の分散分析をaov関数を使って行う手順

処理区ごとの繰り返し数が異なる場合は線形モデルを使った分散分析を行う。(繰り返し数が同じでも問題ない)これを行うのが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%水準で有意差がある。

aov関数を使わずに、変動を分離して行う方法

繰り返しのある三元配置データで因子Aがa水準、Bがb水準、Cがc水準、繰り返しがr回の場合、一般にモデルを次のように構造化する。

Xijkl = μ+αijk+(αβ)ij+(βγ)jk+(αγ)ik+(αβγ)ijkijkl

変動は次のように分解できる。

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)=1124.4852124.4852
肥料(b-1)=267420.8733710.44
年次(c-1)=3
密度×肥料(a-1)(b-1)=22354.1031177.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)
hkokubun at faculty.chiba-u.jp