2013年生物測定基礎実験

統計解析４

東京大学大学院農学生命科学研究科　大森宏

この実験の目的

　統計解析ソフトＲを用いて，統計解析の理論と実践を学ぶ

分散分析（Analysis of Variance）

　分散分析は，ANOVA (Analysis of Variance) と略記されることもある．分散分析は，複数の処理を同時に行ったときに，処理効果を推定するための最も基本的な手法である．データ全体の持つ情報は，総平方和にまとめられているが，これを，処理の分散成分（処理平均平方）と誤差の分散成分（誤差平均平方）とに分離して，その大きさを比較することにより，処理の効果を見積もるものである．

因子と水準

　経済学では価格や成長率，工学では作業時間や故障率，農学では収量や抵抗性など，調査研究したい特性を形質(character)という．着目した形質に影響を与えると考えられるもの，例えば，収量では品種，温度，施肥量などを要因または因子(factor)という．要因の影響を調べるためいくつかの品種を用いたり，施肥量に段階を設けたりするが，それを水準(level)という．

一元配置（one-way layout）

構造モデル

　t 検定では，2 つの処理平均の比較を行ったが，この節ではこれを拡張して，複数の処理平均の比較を行う手法を考える．いま，a 水準の処理（treatment）A₁，…，A_a，があり，処理 A_i を行った n_i 個の標本， X_i1，…，X_{in_i}，が得られたとする．処理 A_i からの標本は，平均 μ_i = μ + α_i，分散 σ² の正規分布に従うと仮定する．ここで，μ を総平均（grand mean）， α_i を処理効果（treatment effect）もしくは，主効果（main effect）と言い， Σ_iα_i = 0，である．ここで，平均 0，分散 σ² を持つ誤差項（error term）e_ij を導入し，標本の構造モデル，

として表現すると，データの持つ構造が理解しやすくなる．

平方和分解

　いま，処理 A_i の標本平均を X^-_i.，標本総平均を X^-_.. とすると，これらは，

と計算される．標本全体の持つ情報は，総平方和 S_T（Total Sum of Squares）で表現される．これは，

のように誤差平方和 S_e（Error Sum of Squares）と処理平方和 S_A（Treatment Sum of Squares）とに分解される．これは，積の項が

のように 0 となるからである．
　なお，誤差平方和を群内平方和（within groups sum of squares），処理平方和を群間平方和（between groups sum of squares）と呼ぶことも多い．

平方和の期待値

　個々の標本 X_ij と処理 A_i の標本平均 X^-_i.，標本総平均 X^-_.. の構造モデルがそれぞれ，

のようになるので，誤差平方和 S_e と処理平方和 S_A の期待値は，それぞれ，

のように計算できる．

帰無仮説のもとでの平方和の比の分布

　一元配置モデルにおける帰無仮説は，すべての処理効果がない，つまり，

H₀：α_i = 0，i = 1，…，a，　

である．前節の平方和の期待値から，帰無仮説のもとで，S_e/σ² は自由度 n - a の χ² 分布に従い，S_A/σ² は自由度 a - 1 の χ² 分布に従うことがわかる．これらの χ² 分布をその自由度で割った比の F 値は，

のように自由度 a - 1，n - a の F 分布に従う．
　ここで，M_A は，処理平方和 S_A をその自由度 a - 1 で割ったもので，処理平均平方（treatment mean square）と呼ばれ，処理平均から求めた誤差分散 σ² の推定値である．一方，M_e は，誤差平方和 S_e をその自由度 n - a で割ったもので，誤差平均平方（error mean square）と呼ばれ，誤差分散の推定値である．
　帰無仮説のもとでは M_A と M_e はほぼ等しいことが期待されるので，その比 F 値は 1 に近いことが期待される．よって，F 値が大きな値をとるときは帰無仮説が正しくないと考え，帰無仮説を棄却する．F 値が大きいか小さいかの判断基準が対応する自由度の F 分布で決められる．

分散分析表と F 検定

　一元配置モデルの解析結果は，以下の分散分析表（ANOVA table）にまとめられる．

変動因	自由度（df）	平方和（S.S.）	平均平方（M.S.）	F 値
主効果	a - 1	S_A	M_A = S_A/(a - 1)	M_A/M_e
誤　差	n - a	S_e	M_e = S_e/(n - a)
全　体	n - 1	S_T

　この表から検定統計量 F 値が求められる．そして，自由度 a - 1，n - a の F 分布の 1 - γ 点（例えば 95 ％点）F（a - 1，n - a）_{1 - γ} より F 値が大きい，すなわち，

F ＞ F（a - 1，n - a）_{1 - γ}

であるとき，帰無仮説を棄却すると，有意水準 γ （例えば 5 ％）の検定が行える．これを，F 検定（F test）という．

多重比較（multiple comparison）

　分散分析（正確には実験計画）の文脈では，試験設計の段階で帰無仮説の設定が行われる．つまり，検定の内容が事前に決定されている．このような「先付け」のときは，検定の数がそれほど多くないなら，複数の検定を行っても有意水準についての補正を行わないのが普通だと思われる．
　しかし，データが得られた後，「後付け」でどの処理間の差が有意であるか調べたい誘惑にかられることが多い．結果として差が大きかった処理間で t 検定を繰り返して行うと，たくさんの検定を行うので，たまたま有意になる確率が名目上の有意水準（たとえば 5 ％）を超えてしまう恐れがある．これが，多重比較である．現在では，コンピュータにより多くの検定を簡単に行うことができるので，以前に比べて多重比較の問題を考慮しなければならないと考えられる．
　いま，処理平均 μ_i と μ_j の比較を行う場合を考える． 2 つの処理平均の差 μ_i - μ_j は， d_ij = X^-_i. - X^-_j. で推定される．帰無仮説（α_i = α_j = 0）のもとで，d_ij の平均と分散は，

となるので，分散 σ² をその推定量 s² で置き換えた検定統計量 t_ij は，

のように自由度 n - a の t 分布に従うので d_ij = 0 の検定を行うことができる．
　有意水準 α'（たとえば 5 ％）の検定は，自由度 n - a の t 分布の 1 - α'/2（たとえば 97.5 ％) 分位点， t(n - a)_{1 - α'/2} を用いて，

が成り立つとき μ_i と μ_j の効果に違いがあると判定される．ここで， LSD（Least squared distance）は最小有意差という量で，以前は，α' = 0.05 として，処理効果のある組み合わせを見つけるためよく用いられていたが，最近は，多重比較を考慮に入れた有意水準の補正を考えるのが普通なので，単純な LSD は使用しない方が良いと思われる．
　いま，a 水準の主効果があったとすると，すべての組み合わせは r = a(a - 1)/2 通りあり，「後付け」の検定を行うときは，全体で r 回の検定を行っていると考えなければならない． R でも対比較では多重比較による有意確率の補正が簡単に行える．

なにもしない
推奨されない方法である．すなわち，補正なしの t 検定を行う．昔の LSD である．
R では，|d_ij| の p 値が出力される．
ボンフェローニ（Bonferroni）補正
いま，有意水準 α' のそれぞれ独立な検定を r 回行ったとすると，1 回の検定で正しい判断を行う確率が 1 - α' なので，r 回の検定で正しい判断を行う確率は，(1 - α')^r となる．よって，正しい判断を行わない（第 1 種の過誤の）確率は，
1 - (1 - α')^r ≒ 1 - (1 - rα') = rα'，ただし，α' ≒ 0
となる．これが，r 回の検定全体での有意水準となる．よって，検定全体での有意水準を α にするには， 1 回の検定の有意水準を α' = α/r にすればよい．これがボンフェローニ補正である．しかし，多重比較における検定は独立な検定ではないので，この補正は厳しすぎ（保守的）て，有意な組み合わせが見つからない恐れがある．
R の多重比較では，補正なしの p 値を r 倍した p 値を出力する．ただし，これが 1 を超えた場合は 1 とする．
ホルム（Holm）補正
ボンフェローニ補正を改良したものである．すべての比較組み合わせ（対比）の t 値を計算し，それを大きさの順に並べる．一番大きな t 値 t₍₁₎ の有意確率を α/r，次の大きさの t₍₂₎ の有意確率を α/(r - 1)，というように有意確率を調整する．
R ではホルム補正がデフォルトで， p 値を大きさの順に並べ最も小さな p 値を r 倍し，次に大きな p 値を r - 1 倍して出力するようである．
チューキー（Tukey）の HSD（honestly significant difference）
今までは，t 検定の有意確率を補正することにより，多重比較の問題に対処していたが，スチィーデント化された範囲の分布（Studentized range distribution）という多重比較専用の分布を用いて検定する．2 つの処理 i，j 間の比較を行うときに用いる検定統計量は，先ほどの t_ij である．

多重比較法の比較

　多重比較の方法はここで取り上げた以外の手法も知られているが，R で手軽に使える手法を解説した．ここで紹介した手法の有意水準をシミュレーションで比較してみる．標準正規乱数 50 を 1 から 5 までのグループに 10 個ずつ分ける．このデータを処理数 a = 5，処理内標本 n = 10 の一元配置データとみなす．このデータは帰無仮説が真のときで，処理平均に差がないはずである．

シミュレーションデータに対する分散分析と多重比較

a <- 5; n <- 10 # 処理水準数 a，処理内標本数 n 　 x <- NULL # 　 for(i in 1:a) x <- c(x, rep(i, n)) # グループラベル　 x <- factor(x) # ラベル化 y <- rnorm(n*a) cbind(y, x) av <- aov(y ~ x) # 分散分析 summary(av) pairwise.t.test(y, x) # 対比較ホルム補正 pairwise.t.test(y, x, p.adj = "bonf") # 対比較ボンフェローニ補正 pairwise.t.test(y, x, p.adj = "none") # 対比較補正なし TukeyHSD(av) # チューキー HSD 　

シミュレーションデータに対する分散分析と多重比較
a <- 5; n <- 10 # 処理水準数 a，処理内標本数 n 　 x <- NULL # 　 for(i in 1:a) x <- c(x, rep(i, n)) # グループラベル　 x <- factor(x) # ラベル化 y <- rnorm(n*a) cbind(y, x) av <- aov(y ~ x) # 分散分析 summary(av) pairwise.t.test(y, x) # 対比較ホルム補正 pairwise.t.test(y, x, p.adj = "bonf") # 対比較ボンフェローニ補正 pairwise.t.test(y, x, p.adj = "none") # 対比較補正なし TukeyHSD(av) # チューキー HSD

pairwise.t.test(y, x, p.adj = "none") # 対比較補正なし
だと，本来差が無いはずなのに，差があるとしてしまいゴーストを拾ってしまう．従って，多くの比較を行う場合，補正しないといけない．

品種によるコメの収量の違い

水稲の９品種をそれぞれ，６区画の水田で栽培したときのアール当たりの玄米重量は以下のようであった．このうち，A，B，D，それぞれ，同じ母本から育成された品種であり，C は標準（対照 control）品種である．
このデータは一元配置分散分析で解析できる．処理が品種で，9 水準からなっている．帰無仮説は，

H₀：収量はどの品種も同じである

である．

品種データダウンロード

品種収量の分散分析

hinsyu <- read.csv("hinsyu.csv"); hinsyu # csv データ読み込み n <- nrow(hinsyu) m <- ncol(hinsyu) boxplot(hinsyu) idname <- rep(names(hinsyu), each=n) hinsyu.vec <- as.vector(as.matrix(hinsyu)) hinsyu.fr <- data.frame(data = hinsyu.vec, id = idname) table(hinsyu.fr$id) av <- aov(hinsyu.fr$data ~ hinsyu.fr$id) # 分散分析 summary(av)

品種収量の分散分析
hinsyu <- read.csv("hinsyu.csv"); hinsyu # csv データ読み込み n <- nrow(hinsyu) m <- ncol(hinsyu) boxplot(hinsyu) idname <- rep(names(hinsyu), each=n) hinsyu.vec <- as.vector(as.matrix(hinsyu)) hinsyu.fr <- data.frame(data = hinsyu.vec, id = idname) table(hinsyu.fr$id) av <- aov(hinsyu.fr$data ~ hinsyu.fr$id) # 分散分析 summary(av)

問題１

コメデータの分散分析の結果から帰無仮説の検定を行え．
コメデータの品種効果を固定効果とみなして多重比較を行い，有意な差が認められる組を求めよ．

ピスタチオの長さ

　最初の学生実験で市販のピスタチオの長さの計測を行った．ピスタチオのサイズがコンビニにより異なっているか分散分析で調べる．
．

ピスタチオ計測データダウンロード

ピスタチオの分散分析
pis <- read.csv("pistachioSize.csv") attach(pis) # pisデータの使用を宣言 seven <- which(store=="seven") # seven ブランドのみデータを取得 boxplot(size[seven] ~ group[seven]) res <- lm(size[seven] ~ group[seven]) summary(res) anova(res) # boxplot(size ~ store) res <- lm(size ~ store) summary(res) anova(res) detach(pis) # pisデータ使用終了

ピスタチオの分散分析

pis <- read.csv("pistachioSize.csv")
attach(pis)			# pisデータの使用を宣言
seven <- which(store=="seven")			# seven ブランドのみデータを取得
boxplot(size[seven] ~ group[seven])
res <- lm(size[seven] ~ group[seven])
summary(res)
anova(res)
#
boxplot(size ~ store)
res <- lm(size ~ store)
summary(res)
anova(res)
detach(pis)			# pisデータ使用終了

問題２
ピスタチオ計測データの分析を行え．ブンランドごとの違いや，同じブランド内での袋ごとの違いなどを解析する．

多元配置

　比較したい因子が複数ある場合は，多元配置となる．このとき，どのような構造モデルにするかは，農場実習第１回で説明した実験計画による．
　1996年度から2000年度(1998年度を除く)までの試験目的は，品種，肥料水準，栽植密度が水稲の収量構成要素に与える影響を計測するためであった．その配置は２反復２要因乱塊法であった．

2000年度試験配置
 データダウンロード

　収量とは，単位面積(m2)あたりの玄米の収穫量(g)， (grain yield; gy)のことであるが，収量構成要素は，収量に影響を与える形質で収量を分解したもので，たとえば，「単位面積あたりの穂数(spikelet number; sn)」，「一穂籾数(grain number; gn)」，「登熟歩合(percentage of riped grains; pr)」，「玄米千粒重(g)，(1000-grain weight; gw)」と分解できる．収量はこれらの構成要素の積で表現される．
　今回の実験では「日本晴」データの解析を行う．

特定年次の解析

構造モデル

　ある年次における処理組み合わせでラベルをつけた収量値 X_ijk は，処理効果により，

X_ijk = μ + α_i + β_j + ρ_k + (αβ)_ij + ε_ijk, i = 1,2, j = 1,2,3, k = 1,2

と構造化（モデル化）できる．ここで，μ は総平均で α_i は栽植密度の主効果(main effect)であり，β_j は施肥量主効果， ρ_k はブロック効果，(αβ)_ij は栽植密度と施肥量の交互作用(interaction)で，各効果の和は0，たとえば，Σ α_i = 0，とする．また，ε_ijk はこれらの効果で説明できない誤差で，互いに独立に平均 0，分散 σ² の正規分布(normal distribution) に従う (ε_ijk ～ N(0, σ²)) と仮定する．各効果の推定値は，たとえば，

と表せる．

処理平方和

　効果の大きい処理に対しては各効果の水準ごとの推定値は大きく異なり，逆に各水準ごとに同じような値の推定値を与えるような処理は効果が小さいと考えられる．これは処理平方和（sum of squares; SS）で測ることができる．たとえば，密度の主効果の平方和は，　

である．実は，データがもつ全平方和（sum of squares; SS）

は，以下に示すように各処理による平方和で分解され，

となる．これと同様に，各処理平方和に対する自由度（degree of freedom）も

となる．なお自由度とは，比較の数のことで，3 水準の試験では比較は 2 通り（水準 1 対水準 2，水準 2 対水準 3）なので自由度は 2 となる．また，比較 1 つ分の平方和が平均平方（mean squares; MS）で，平方和を自由度で割った値である．たとえば，密度効果の平均平方は， MS_α = S_α/(a-1) である．また，誤差分散は，

のように誤差平均平方で推定する．各効果の大きさは，この誤差の大きさと比較することにより測ることができる．

F 検定

　各処理にまったく効果がないという帰無仮（null hypothesis），つまり，

H₀ : α_i = 0, β_j = 0

という条件のもとで各処理効果の平均平方を誤差平均平方で割った値が F 分布に従うことにより各処理効果の検定(test)が行える．たとえば，密度主効果の検定は，

を用いる．密度効果がなければ F_α は 1 に近い値をとる．大きな値を取った場合は，密度効果がないという帰無仮説が成立しておらず，密度効果が認められると考える．つまり，帰無仮説を棄却(reject)する．この大きさの基準として，慣習的に F 分布(自由度 a - 1, (ab - 1)(c - 1) )の 5 % や 1 % 点が用いられている．この基準確率を有意水準(significance level)という．

**収量に対する処理効果の分散分析表**
要因	自由度	平方和	平均平方	F 値	p 値
ブロック	1	15862	15862	3.85	0.107
密度	1	1419	1419	0.348	0.583
施肥	2	65029	32514	7.898	0.028	*
密度×施肥	2	2222	1111	0.270	0.774
誤差	5	20583	4117

* : 5 % 有意

分散分析表の見方

　各処理の効果の程度をまとめたのが分散分析（analysis of variance; ANOVA）であり，その結果は上の表にまとめられている．これをみると，施肥量の効果の平方和は3水準で S_β=65029 となり，比較は2通りなので自由度は2となる．その平均平方は平方和を2で割り，MS_β=32514 である．この値と誤差平方和の比が F 値で，7.90と大きな値になった．帰無仮説のもとで F 値より大きな値がでる確率が p 値(p - value)である．p 値が小さいことは，帰無仮説が成立する可能性が小さい，つまり，処理効果が存在する可能性が大きいことを意味する．
　この例では，統計的にみて有意水準 3 % で施肥量により収量が異なるといえるが，他の処理効果は認められなかったという結論になる．つまり，平均値から得た類推は，統計的に根拠があったということになる．

コメ収量（2000年）の分散分析

rice <- read.csv("ricecul.csv") 	# データ読み込み
rice
yield <- rice$gy 
dense <- factor(rice$density)		# 水準のラベル化
fert <- factor(rice$fert)
blk <- factor(rice$rep)  
tapply(yield[1:12], dense[1:12], mean)		# 密度水準ごとの平均
tapply(yield[1:12], fert[1:12], mean) 
tapply(yield[1:12], blk[1:12], mean) 
cm <- tapply(yield[1:12], dense[1:12]:fert[1:12], mean); cm	# 処理組み合わせごとの平均
# 施肥量水準に対する収量のグラフ
plot(1:3, cm[1:3], type="b", lwd=2, cex=1.5, xaxt="n", xlab="fertility", ylab="yield",  
ylim=c(300, 600), pch=0, cex.lab=1.0, cex.axis=1.0, col="blue") 
axis(1, 1:3, labels=c("fert1","fert2","fert3"), cex.axis=0.8) 
points(1:3, cm[4:6], type="b", lwd=2, cex=1.5, pch=2, col="red")
legend(1.2, 580, legend=c("density1","density2"), pch=c(0,2), 
col=c("blue","red"), cex=0.8)  
title(main="施肥量水準に対するコメ収量（2000年）")  
# 分散分析
ry.aov <- aov(yield[1:12] ~ blk[1:12] + dense[1:12] + fert[1:12] + dense[1:12]:fert[1:12])
summary(ry.aov)                 #分散分析表の表示

注）blk[1:12] などは主効果であり，dens[1:12]:fert[1:12] は交互作用である．

複数年次にわたる解析

構造モデル

　4年間にわたって行ったデータ表を用いると，気象条件などの年ごとの違いによる効果がわかる．年次効果を γ_k で表すとすると，

X_ijkl = μ + α_i + β_j + γ_k + (αβ)_ij + (βγ)_jk + (αγ)_ik + (αβγ)_ijk + ε_ijkl

という 3 元配置モデルで記述できる．ただし，ブロックの 1 と 2 は毎年管理者が異なるので，ブロック効果として取り出さず，たんなる繰り返しとした．

問題３

田無農場の2000年のデータで，収量構成要素を一つ選び分散分析を行い，結果を考察せよ．
田無農場の４年次にわたるデータで，収量と収量構成要素のどれか一つに対し，分散分析を行い，結果を考察せよ．

参考文献（古い順）

Introduction to the Theory of Statistics, Mood, A. M., Graubill, F. A. & Boes, D. C., 1974, McGRAW-HILL
「実験」生産環境生物学，東京大学大学院農学生命科学研究科生産・環境生物学専攻編，1999，朝倉書店
工学のためのデータサイエンス入門－フリーな統計環境Rを用いたデータ解析－，間瀬茂ら，2004，数理工学社
実践生物統計学－分子から生態まで－（第 1 章，第 2 章），東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
The R Tips データ解析環境 R の基本技・グラフィックス活用集，船尾暢男，2005，九天社
R で学ぶデータマインニング I －データ解析の視点から－，熊谷悦生・船尾暢男，2007，九天社
R で学ぶデータマインニング II －シミュレーションの視点から－，熊谷悦生・船尾暢男，2007，九天社

2013年生物測定基礎実験

統計解析４

東京大学大学院農学生命科学研究科 大森宏

この実験の目的

分散分析（Analysis of Variance）

因子と水準

一元配置（one-way layout）

構造モデル

平方和分解

平方和の期待値

帰無仮説のもとでの平方和の比の分布

分散分析表と F 検定

多重比較（multiple comparison）

多重比較法の比較

品種によるコメの収量の違い

ピスタチオの長さ

多元配置

特定年次の解析

構造モデル

処理平方和

F 検定

分散分析表の見方

複数年次にわたる解析

構造モデル

参考文献（古い順）

東京大学大学院農学生命科学研究科　大森宏