実験計画法
因子と水準
作物において，収量や抵抗性などの調査研究したい特性を形質 (trait) という．着目した形質に影響を与えると考えられる品種，温度，施肥量などを要因または因子 (factor) という．要因の影響を調べるためいくつかの品種を用いたり，施肥量に段階を設けたりするが，それを水準 (level) という．
要因と水準を決定したら，その反応を調べるために圃場で実際に栽培してみる必要がある．しかしながら，圃場の環境条件は場所により異なり，均一の条件で栽培管理を行うことは難しい．環境条件に多少のばらつきがあると想定される場合には，ばらつきの影響を無作為化する必要がある．これを行うのが，実験計画 (design of experiment) である．
Fisherの３原則
実験計画の創始者である R.A. Fisher が提唱した原則で，反復 (replication)，無作為化 (randomization)，局所管理 (local control) の3つからなる．これらは，処理の違いに基づく差が偶然の誤差を超えて有意なものかどうか統計的に判定するため，実験を行うにあたり必要な管理原則である．この原則に基づいて圃場での実験配置が設計される．
反復 (replication) 1 つの処理組み合わせに対して，少なくとも 2 回以上の繰り返し (repetition) が必要である．繰り返しにより，誤差の大きさが評価できる．
無作為化 (randomization) 実験順序や空間配置の無作為化．ある特定の処理条件を圃場のある特定の場所にかためてしまうと，その場所がたまたま水はけがよかったり肥沃だったりした場合，処理効果によるものか場所の効果によるものかの判別ができなくなる． このようなある特定の場所による効果を系統誤差 (systematic error) というが，無作為に処理条件を配置させることにより，系統誤差を偶然誤差 (random error) に転化させることができる．
局所管理 (local control) 実験が大規模で，実験全体を無作為化したら誤差が大きくなってしまうような場合， 実験をある程度細分化してブロックを構成し，ブロック内で処理条件を無作為化する． ブロック内のバックグラウンドはなるべく均一になるように管理できれば，系統誤差の一部がブロック間変動として除去できるので，実験誤差が小さくなる．
実験計画はなぜ必要か
完全無作為化法
ある地域ではイネ 4 品種の収量は平均的に
| 品種 A | 品種 B | 品種 C | 品種 D |
|---|---|---|---|
| 60\(kg/a\) | 61\(kg/a\) | 58\(kg/a\) | 57\(kg/a\) |
であることが知られている。ここで、下のような 8 区画の圃場でイネ 4 品種の収量を調べる試験を行うことを考えてみる。
しかしながら、この圃場の地力は均一ではなく、本当は以下のように右側にいくほど収量が高くなる傾向があり、下側も収量が多少高くなっていたが、試験者は圃場の地力について知らないとしよう。
いま、この圃場にイネ 4 品種を 2 回の繰り返しで栽培し、収量調査を行うことにした。ここで良くない配置の例としては、管理の容易さから品種の番号順に並べた
である。この配置では品種 A は最も地力の低い場所に栽培され、品種 D は最も地力の高い場所に栽培されている。このため、品種の収量の大小と地力の大小が交絡 (confound) し、品種間の収量の差がわからなくなっている。
このような場所でも、無作為化することで系統的な地力の差異を偶然誤差に変えることができる。圃場の 8 区画にランダムに品種を割り付けたとすると、これは完全無作為化法 (completely randomized design) と呼ばれる。
いま、\(i\) 番目の品種の \(j\) 番目の繰り返し（「反復」と呼ぶのはあまり適当ではない）の収量を \(y_{ij}\) とする。ここで、総平均(grand mean)を \(\mu\)、\(i\) 番目の品種効果(variety effect)を \(\alpha_i\)、ただし \(\sum_i \alpha_i = 0\)、誤差項(error term)を \(e_{ij}\)、ただし \(e_{ij} 〜 N(0, \sigma^2)\) で互いに独立である、とおくと、 \[ y_{ij} = \mu + \alpha_i + e_{ij},\ i=1,\ldots, 4,\ j=1,2, \\ \sum_i \alpha_i = 0, \ e_{ij} \sim N(0, \sigma^2) \] とおける。以下に配置例を示す。
この例では品種 B はたまたま比較的地力の低い場所に配置され、品種 C は比較的地力の高い場所に配置されて偏りがでている。
さて、誤差標準偏差を \(σ= 0.5\) と仮定して、上記の配置から各品種の収量がどのように推定されるかをみてみよう。すなわち、区画の収量データは、そこに供試された品種収量にその区画の地力と正規乱数誤差を加えたものがシミュレートされる。このデータを完全無作為化法であるとみなして品種効果を推定してみよう。
まずは品種収量と地力を定義し、固定した（毎回同じ）正規乱数誤差を設定する。
A=60; B=61; C=58; D=57 # 品種収量の定義
yield <- c(A,B,C,D)
true.vef <- yield - mean(yield); true.vef # 真の品種効果## [1] 1 2 -1 -2fert1 <- 0:3; fert2 <- fert1+0.2
land <- c(fert1, fert2); land # 地力## [1] 0.0 1.0 2.0 3.0 0.2 1.2 2.2 3.2set.seed(41) # 乱数を固定するためのシード
error <- rnorm(8, sd=0.5); error # 正規乱数誤差## [1] -0.39718417 0.09862877 0.50085216 0.64441270 0.45287671 0.24683374
## [7] 0.29964290 -0.78980350良くない配置例は、
x <- factor(rep(1:4, 2)); x # 良くない配置例## [1] 1 2 3 4 1 2 3 4
## Levels: 1 2 3 4であり、完全無作為化法による配置例は、たとえば、
set.seed(15)
a0 <- sample(1:8) # 無作為に配置を決めるための並べ替え
x <- (a0 %% 4) + 1
x <- factor(x); x # 完全無作為化法による配置例## [1] 2 4 3 3 2 1 1 4
## Levels: 1 2 3 4となる。
配置例 \(x\) が決まると、品種の収量は以下のように推定される。
rice <- yield[x] # 区画ごとの期待コメ収量
y <- rice + land + error; y # 観測される収量のシミュレート値## [1] 60.60282 58.09863 60.50085 61.64441 61.65288 61.44683 62.49964 59.41020gm <- mean(y); gm # 総平均## [1] 60.73203vm <- tapply(y, x, mean) # 品種ごとの平均（yをxごとにmeanする）
est.vef <- vm - gm; est.vef # 推定された品種効果## 1 2 3 4
## 1.2412059 0.3958139 0.3406000 -1.9776198となる。また、各種平方和や分散分析は以下のように計算される。
gmv <- rep(gm, 8) # 総平均ベクトル
vefv <- est.vef[x] # 品種効果ベクトル
errv <- y-gm-vefv # 誤差ベクトル
cbind(y, gm, vefv, errv)## y gm vefv errv
## 2 60.60282 60.73203 0.3958139 -0.5250304
## 4 58.09863 60.73203 -1.9776198 -0.6557839
## 3 60.50085 60.73203 0.3406000 -0.5717803
## 3 61.64441 60.73203 0.3406000 0.5717803
## 2 61.65288 60.73203 0.3958139 0.5250304
## 1 61.44683 60.73203 1.2412059 -0.5264046
## 1 62.49964 60.73203 1.2412059 0.5264046
## 4 59.41020 60.73203 -1.9776198 0.6557839sum((y - gmv)^2) # 総平方和## [1] 14.06799sum(vefv^2) # 品種効果（処理）平方和## [1] 11.4485sum(errv^2) # 誤差平方和## [1] 2.619488summary(aov(y ~ x)) # 分散分析結果## Df Sum Sq Mean Sq F value Pr(>F)
## x 3 11.448 3.816 5.827 0.0608 .
## Residuals 4 2.619 0.655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1推定された品種効果は以下のようになった。
| 品種 | A | B | C | D |
|---|---|---|---|---|
| 真の品種効果 | 1 | 2 | -1 | -2 |
| 良くない配置 | -0.604 | 1.541 | -0.232 | -0.705 |
| 完全無作為 | 1.241 | 0.396 | 0.341 | -1.978 |
これをみると、良くない配置では品種Aの効果は過少評価され、品種 C、D は過大評価されたことがわかる。これは、地力の差異を品種効果と見間違えたためである。一方、完全無作為法では多少正確な効果推定がなされたが、品種 B は過少評価、品種 C は過大評価されている。
両者を分散分析表で比較してみよう。分散分析の帰無仮説は、 \[ {\rm H}_0 : \alpha_1 = \cdots = \alpha_4 = 0 \] で、品種効果は存在しない、である。
良くない配置
| 要因 | 自由度 | 平方和 | 平均平方 | F値 | p値 |
|---|---|---|---|---|---|
| 品種効果 | 3 | 6.578 | 2.1928 | 6.386 | 0.0526 |
| 誤差 | 4 | 1.374 | 0.3434 | ||
完全無作為配置
| 要因 | 自由度 | 平方和 | 平均平方 | F値 | p値 |
|---|---|---|---|---|---|
| 品種効果 | 3 | 11.448 | 3.816 | 5.827 | 0.06086 |
| 誤差 | 4 | 2.619 | 0.655 | ||
前にも述べたように、データの持つ全情報量は総平方和で表現される。「良くない配置」における総平方和は、品種効果平方和の 6.578 と誤差平方和 1.374 を加えた 7.952 である。実際、品種効果は 5％有意では無く、品種間の収量差を検出できなかった。
一方、「完全無作為法」による総平方和は、11.448 + 2.619 = 14.067、となり大きく増加した。これは、「良くない配置」では、品種効果と地力効果が相殺されてしまい、データ変動全体が小さくなってしまったと考えられる。しかしながら、この方法でも品種効果は 5％有意でなく、品種間差が検出できなかった。
1要因乱塊法
完全無作為化法では品種 B は地力の低い区画に配置され、品種 C は地力の高い区画に配置されるなどアンバランスになる場合がある。そこで、隣接した 4 区画を 1 つのブロックとしてまとめて局所管理する方法が考えられる。これは Fisher の 3 原則を完全に満足する配置で、乱塊法 (randomized block design) と呼ばれる。
いま、\(j\) 番目のブロック効果を \(\rho_j\) とおく。ここで、ブロック効果は変量と考え \(\rho_j \sim N(0, \sigma_B^2)\) とおく。ブロックは系統的な地力の違いを表現するとみなせば、固定効果と考えてもよい。ブロックにはすべての試験組み合わせが入るので、これを \(j\) 番目の反復と呼ぶ。このとき、収量を \(y_{ij}\) は、 \[ y_{ij} = \mu + \alpha_i + \rho_j + e_{ij}, \\ \sum_i \alpha_i = 0, \ \rho_j \sim N(0, \sigma_B^2), \ e_{ij} \sim N(0, \sigma^2) \] とおける。これは、繰り返しの無い 2 元配置の形をしている。
ブロックの作り方としては、パターン 1 ({1, 2, 3, 4}, {5, 6, 7, 8}) とパターン 2 ({1, 2, 5, 6}, {3, 4, 7, 8}) の 2 通りが考えられる。パターン 1 とパターン 2 の配置例は、乱数によりそれぞれ、
となった。パターン 1 では行ごとにブロックを形成したため、ブロック内の地力は一様とは言えず、品種 A が地力の低い場所に配置され、品種 D は地力の比較的高い場所に配置されていて偏りは解消されていない。一方、パターン 2 は左右でブロックを形成したので、パターン 1 に比べるとブロック内の地力の差異は少なくなっている。さて、先ほどの例と同様に、上記の配置から各品種の収量がどのように推定されるかをみてみよう。
set.seed(44)
a1 <- sample(1:4); a1 # ブロック1での並べ替え## [1] 1 3 4 2a2 <- sample(1:4); a2 # ブロック2での並べ替え## [1] 2 1 3 4パターン1の配置は、
x <- factor(c(a1,a2)); x # パターン1の配置## [1] 1 3 4 2 2 1 3 4
## Levels: 1 2 3 4bl <- rep(1:2, each=4)
bl <- factor(bl); bl # パターン1のブロック## [1] 1 1 1 1 2 2 2 2
## Levels: 1 2と定義され、パターン2の配置は、
x <- factor(c(a1[1:2],a2[1:2],a1[3:4],a2[3:4])) ; x # パターン2の配置## [1] 1 3 2 1 4 2 3 4
## Levels: 1 2 3 4bl <- rep(rep(1:2, each=2), 2)
bl <- factor(bl); bl # パターン2のブロック## [1] 1 1 2 2 1 1 2 2
## Levels: 1 2と定義される。
配置 \(x\) とブロック定義 \(bl\) が決まると、品種の収量は以下のように推定される。
rice <- yield[x] # 区画ごとの期待コメ収量
y <- rice + land + error; y # 観測される収量のシミュレート値## [1] 59.60282 59.09863 63.50085 63.64441 57.65288 62.44683 60.49964 59.41020gm <- mean(y); gm # 総平均## [1] 60.73203vm <- tapply(y, x, mean) # 品種ごとの平均
est.vef <- vm - gm; est.vef # 推定された品種効果## 1 2 3 4
## 0.8915819 2.2418105 -0.9328966 -2.2004958blm <- tapply(y, bl, mean) - gm; blm # ブロック効果## 1 2
## -1.031744 1.031744となる。また、各種平方和や分散分析は以下のように計算される。
gmv <- rep(gm, 8) # 総平均ベクトル
vefv <- est.vef[x] # 品種効果ベクトル
blv <- blm[bl] # ブロック効果ベクトル
errv <- y - gmv - vefv - blv # 誤差ベクトル
cbind(y, gmv, vefv, blv, errv)## y gmv vefv blv errv
## 1 59.60282 60.73203 0.8915819 -1.031744 -0.9890548
## 3 59.09863 60.73203 -0.9328966 -1.031744 0.3312366
## 2 63.50085 60.73203 2.2418105 1.031744 -0.5047344
## 1 63.64441 60.73203 0.8915819 1.031744 0.9890548
## 4 57.65288 60.73203 -2.2004958 -1.031744 0.1530838
## 2 62.44683 60.73203 2.2418105 -1.031744 0.5047344
## 3 60.49964 60.73203 -0.9328966 1.031744 -0.3312366
## 4 59.41020 60.73203 -2.2004958 1.031744 -0.1530838sum((y-gmv)^2) # 総平方和## [1] 34.31446sum(vefv^2) # 品種効果（処理）平方和## [1] 23.06622sum(blv^2) # ブロック平方和## [1] 8.51596sum(errv^2) # 誤差平方和## [1] 2.732277summary(aov(y ~ x + bl)) # 分散分析結果## Df Sum Sq Mean Sq F value Pr(>F)
## x 3 23.066 7.689 8.442 0.0566 .
## bl 1 8.516 8.516 9.350 0.0551 .
## Residuals 3 2.732 0.911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1推定された品種効果とブロック効果は以下のようになった。
| 品種 | A | B | C | D | ブロック 1 | ブロック 2 |
|---|---|---|---|---|---|---|
| 真の品種効果 | 1 | 2 | -1 | -2 | ||
| パターン1 | -0.207 | 2.417 | -0.933 | -1.277 | -0.020 | 0.020 |
| パターン2 | 0.892 | 2.242 | -0.933 | -2.200 | -1.032 | 1.032 |
これをみると、パターン 1 では品種 A は地力の低い区画に割り付けられたので過少評価されているが、パターン 2 では品種の効果がよく推定されているのがわかる。これは、パターン 1 のブロックは地力の差異を肩代わりしていないが、パターン 2 では地力の差異をブロックの差異にある程度置き換えることができたからである。 実際、左側のブロックと右側のブロックでは平均的に 2 だけ地力の差があったが、データから推定されたブロック間差は 2.06 となり、地力の違いをうまく反映していた。
両者を分散分析表で比較してみよう。分散分析の帰無仮説は、 \[ {\rm H_T} : \alpha_1 = \cdots = \alpha_4 = 0, \ {\rm H_B} : \sigma_B^2 = 0 \] で、品種効果は存在しないとブロック効果は無い、の２つの帰無仮説を同時に検定する。
パターン1
| 要因 | 自由度 | 平方和 | 平均平方 | F値 | p値 |
|---|---|---|---|---|---|
| 品種効果 | 3 | 16.765 | 5.588 | 2.343 | 0.251 |
| ブロック効果 | 1 | 0.003 | 0.003 | 0.001 | 0.973 |
| 誤差 | 3 | 7.157 | 2.386 | ||
パターン2
| 要因 | 自由度 | 平方和 | 平均平方 | F値 | p値 |
|---|---|---|---|---|---|
| 品種効果 | 3 | 23.066 | 7.689 | 8.442 | 0.0566 |
| ブロック効果 | 1 | 8.516 | 8.516 | 9.530 | 0.0551 |
| 誤差 | 3 | 2.732 | 0.911 | ||
パターン 1 では、ブロック間に差がないので、地力の差と品種効果が交絡したままになっている。このため、誤差平方和が大きいままであり、品種効果が検出されなかった。一方、パターン 2 ではたまたま左側ブロックは地力が低く、右側ブロックが地力が高くなっているので、地力の差がうまくブロック間差に吸収された。このため誤差平方和が小さくなり、5％有意とはならなかったが、品種効果をある程度検出することができた。
2 因子要因実験
前節では実験配置の基本的考えを述べたが、取り上げた要因は品種効果の 1 因子実験であった。実際の実験では、施肥量と栽植密度、品種との関係など実験要因が複数組み合わせることがよくある。このような実験では、要因ごとの主効果 (main effect) だけでなく要因間の交互作用 (interaction) が検出できる。2 因子要因実験では 2 元配置が普通であるが、乱塊法を用いると、誤差平方和からブロック平方和が取り除けるので、有意差が検出しやすくなり、より効率的になると考えられる。
2 要因乱塊法
Fisher の 3 原則を完全に満足する実験手法で，すべての処理組み合わせの実験を1回ずつ集めたもので1 つのブロックを形成する．ブロック数が反復の数になる．いま、あるイネ品種の収量が施肥量と栽植密度の違いでどのように異なるかを調べる実験を行うとしよう。施肥量としては、無、少、多 (0, 1, 2) の 3 水準、栽植密度は疎植、密植 (S, D) の 2 水準を考えたとする。このときの処理組み合わせは 3×2 = 6 通りとなり、これが 1 ブロックとなる。2 反復の実験なら 6×2 = 12 区画、3 反復なら 6×3 = 18 区画の圃場が必要となる。
いま、\(i\) 番目の栽植密度で \(j\) 番目の肥料条件、\(k\) 番目のブロックで栽培されたイネのある形質値を \(y_{ijk}\) とすると、\(\mu\) を総平均、\(\alpha_i\) を栽植密度効果、 \(\beta_j\) を施肥量効果、\((\alpha \beta)_{ij}\) を密度と施肥量の交互作用、\(\rho_k\) をブロック効果、 \(e_{ijk}\) を誤差として、 \[ y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha \beta)_{ij} + \rho_k + e_{ijk}, \ i=1,\ldots,a, \ j=1,\ldots,b, \ k=1,\ldots,r \\ \sum_i \alpha_i = 0, \ \sum_j \beta_j = 0, \ \sum_i (\alpha \beta)_{ij} = \sum_j (\alpha \beta)_{ij} = 0, \\ \rho_k \sim N(0, \sigma_B^2), \ e_{ijk} \sim N(0, \sigma^2) \] とおける。
3 反復の場合の試験配置例を以下に示す。
東京大学農学部では 3 年生に対して，附属多摩農場で農場所属の教官や技官の指導により 水稲栽培の圃場試験実習を 行っている． 1996 年度から 2000 年度 (1998 年度を除く) までの試験目的は，品種，施肥量，栽植密度が 水稲の収量構成要素に与える影響を計測するためであった．
収量とは，単位面積 (\(m^2\)) あたりの玄米の収穫量 (\(g\))， (grain yield; gy) のことであるが，収量構成要素は，収量に影響を与える形質で収量を分解したもので，たとえば，「単位面積あたりの穂数 (panicle number; pn)」，「一穂籾数 (grain number; gn)」，「登熟歩合 (percentage of ripened grains; pr)」，「玄米千粒重 (\(g\))，(1000-grain weight; gw)」と分解できる．収量はこれらの構成要素の積で表現される．
今回の演習では「日本晴」データを用いる。このとき、年度ごとの実験では 2 反復の施肥量と栽植密度の 2 要因実験になっている。いま、1999 年度のデータで対象形質を玄米千粒重 (gw) として解析してみよう。
rice <- read.csv("data/ricecul.csv") # 日本晴データ読み込み
t <- which(rice$year==1999); rice[t,]## year plot variety density fert rep pn gn pr gw
## 13 1999 2 1 1 1 1 211.5 103.23404 0.7953421 26.19850
## 14 1999 14 1 1 1 2 202.5 101.22222 0.8349067 25.84801
## 15 1999 1 1 1 2 1 211.5 95.70213 0.8928413 26.65339
## 16 1999 20 1 1 2 2 193.5 92.79070 0.8704261 27.12352
## 17 1999 5 1 1 3 1 234.0 90.01923 0.8850673 26.74632
## 18 1999 19 1 1 3 2 252.0 99.89286 0.8246335 26.59874
## 19 1999 4 1 2 1 1 224.0 87.50000 0.8053571 25.32151
## 20 1999 13 1 2 1 2 232.4 84.69880 0.9046942 26.19497
## 21 1999 10 1 2 2 1 259.0 98.94595 0.8128927 20.83669
## 22 1999 16 1 2 2 2 266.0 95.68421 0.8894389 23.93321
## 23 1999 3 1 2 3 1 245.0 89.29143 0.8364265 26.07498
## 24 1999 24 1 2 3 2 224.0 76.56250 0.8685714 27.53759
## gy
## 13 382.1580
## 14 371.5740
## 15 404.6112
## 16 356.0760
## 17 418.8618
## 18 463.8060
## 19 335.7480
## 20 391.8432
## 21 364.6188
## 22 455.1120
## 23 400.7808
## 24 344.5680density <- factor(rice$density[t])
fert <- factor(rice$fert[t])
rep <- factor(rice$rep[t])
y <- rice[t, "gw"] # 玄米千粒重
gm <- mean(y); gm # 総平均## [1] 25.75562fm <- tapply(y, fert, mean); fm # 施肥量平均## 1 2 3
## 25.89075 24.63670 26.73941dm <- tapply(y, density, mean); dm # 密度平均## 1 2
## 26.52808 24.98316fdm <- tapply(y, density:fert, mean); fdm # 密度、施肥量組合せ平均## 1:1 1:2 1:3 2:1 2:2 2:3
## 26.02326 26.88846 26.67253 25.75824 22.38495 26.80629summary(aov(y ~ density*fert+rep)) # ２要因乱塊法分散分析## Df Sum Sq Mean Sq F value Pr(>F)
## density 1 7.160 7.160 8.964 0.0303 *
## fert 2 8.952 4.476 5.604 0.0529 .
## rep 1 2.434 2.434 3.047 0.1413
## density:fert 2 13.209 6.605 8.268 0.0260 *
## Residuals 5 3.994 0.799
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1また、この例では分散分析により交互作用が検出された。これは、下の密度、施肥量組合せ平均のグラフでわかるように、密植で肥料水準が少のときだけなぜか千粒重が低下したためである。
plot(fdm[1:3], type="b", xaxt="n", ylim=c(20,30), xlab="fertility", ylab="gw")
lines(fdm[4:6], type="b", col="red")
axis(1, 1:3, labels=0:2)
legend("bottomright", legend=c("sparse","dense"), pch=1, col=c("black", "red"))
分割区法(split-plot design)
潅水区や非潅水区のような水管理は，比較的広い区画で行わないと労力，費用等がかかる．このような処理を用いた実験の場合，乱塊法でブロック内すべての組み合わせの試験区を設定することは現実的でない．水管理のように広い区画が必要な因子を 1 次因子，その区画内で水準を変化させる品種や施肥量などの因子を 2 次因子と呼ぶ．
先ほどの乱塊法の例で、栽植密度を潅水に変えたものとしてみてみよう。潅水区、非潅水区 (S, D) を 3 反復の分割区法に割り付けたとすると、以下のような配置例となる。
1 次因子の配置は 1 要因乱塊法になっていて、この配置の中で 2 次因子がランダムに配置される。
いま、\(k\) 番目のブロックで、1次因子として \(i\) 番目の潅水条件で、2 次因子として \(j\) 番目の肥料条件で栽培されたイネの収量を \(y_{ijk}\) とするとこれは、
\[ y_{ijk} = \mu + \rho_k + \alpha_i + \epsilon^1_{ik} + \beta_j + (\alpha \beta)_{ij} + \epsilon^2_{ijk}, \\ \epsilon^1_{ik} \sim N(0, \sigma_1^2), \ \epsilon^2_{ijk} \sim N(0, \sigma_2^2) \] とおける。ここで、\(\epsilon^1_{ik}\) は 1 次誤差でありこの誤差を用いて潅水条件の検定を行う。また、\(\epsilon^2_{ijk}\) は 2 次誤差でありこの誤差で施肥量効果や潅水条件と施肥量の交互作用の検定を行う。
乱塊法の平方和分解との違いは、括弧内を平方和の自由度として、誤差平方和 (5) が 1 次誤差平方和 (1) と 2 次誤差平方和 (4) に分解されることである。このため乱塊法に比べ、分割区法では 1 次因子主効果は有意になりにくいが、2 次因子の主効果や交互作用は有意になりやすくなる（検出力が増す）と思われる。
1999 年度の玄米千粒重 (gw) のデータを 2 要因乱塊法で解析したが、密度効果を潅水効果とみなして 1 次因子とし、施肥量効果が 2 次因子である分割区法とみなして解析してみよう。
drm <- tapply(y, density:rep, mean) # 潅水、ブロック組合せ平均
drm## 1:1 1:2 2:1 2:2
## 26.53273 26.52343 24.07773 25.88859yy <- rep(0, length(y))
for(i in 1:length(drm))
yy[which(density:rep==names(drm)[i])] <- drm[i]
yy## [1] 26.53273 26.52343 26.53273 26.52343 26.53273 26.52343 24.07773 25.88859
## [9] 24.07773 25.88859 24.07773 25.88859summary(aov(yy ~ density + rep)) # 1次因子平方和## Df Sum Sq Mean Sq F value Pr(>F)
## density 1 7.160 7.160 25.935 0.000652 ***
## rep 1 2.434 2.434 8.817 0.015720 *
## Residuals 9 2.485 0.276
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(y - yy ~ fert + density:fert)) # 2次因子平方和## Df Sum Sq Mean Sq F value Pr(>F)
## fert 2 8.952 4.476 17.80 0.00300 **
## fert:density 3 13.209 4.403 17.51 0.00227 **
## Residuals 6 1.509 0.252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(y ~ density*fert+ Error(rep/density))) # 分割区法分散分析##
## Error: rep
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 2.434 2.434
##
## Error: rep:density
## Df Sum Sq Mean Sq F value Pr(>F)
## density 1 7.160 7.160 2.882 0.339
## Residuals 1 2.485 2.485
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## fert 2 8.952 4.476 11.86 0.0208 *
## density:fert 2 13.209 6.605 17.51 0.0105 *
## Residuals 4 1.509 0.377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ここで、乱塊法による分散分析表と分割区法による分散分析表を比較の意味でまとめてみよう。乱塊法と分割区法では各処理効果の平方和は同じであるが、誤差平方和が 1 次誤差平方和と 2 次誤差平方和に分解されていることを確認しよう。
乱塊法分散分析表
| 要因 | 自由度 | 平方和 | 平均平方 | F 値 | p 値 |
|---|---|---|---|---|---|
| ブロック効果 | 1 | 2.434 | 2.434 | 3.047 | 0.1413 |
| 密度効果 | 1 | 7.160 | 7.160 | 8.694 | 0.0303 |
| 施肥量効果 | 2 | 8.952 | 4.476 | 5.604 | 0.0529 |
| 密度×施肥交互作用 | 2 | 13.209 | 6.605 | 8.268 | 0.0260 |
| 誤差 | 5 | 3.994 | 0.799 | ||
分割区法分散分析表
| 要因 | 自由度 | 平方和 | 平均平方 | F 値 | p 値 |
|---|---|---|---|---|---|
| ブロック効果 | 1 | 2.434 | 2.434 | 0.979 | 0.503 |
| 密度効果 | 1 | 7.160 | 7.160 | 2.882 | 0.339 |
| １次誤差 | 1 | 2.485 | 2.485 | ||
| 施肥量効果 | 2 | 8.952 | 4.476 | 11.86 | 0.0208 |
| 密度×施肥交互作用 | 2 | 13.209 | 6.605 | 17.51 | 0.0105 |
| ２次誤差 | 4 | 1.509 | 0.377 | ||
多元配置
複数年次にわたる解析
4 年間にわたって行ったデータを用いると、気象条件などの年ごとの違いによる効果がわかる。いま、\(i\) 番目の栽植密度で \(j\) 番目の肥料条件、\(k\) 番目のブロックで栽培された \(l\) 年次におけるイネのある形質値を \(y_{ijkl}\) とする。ここで、\(\mu\) を総平均、\(\alpha_i\) を栽植密度効果、 \(\beta_j\) を施肥量効果、\(\gamma_l\) を年次効果、\((\alpha \beta)_{ij}\) を密度と施肥量の交互作用、\((\alpha \gamma)_{il}\) を密度と年次の交互作用、\((\beta \gamma)_{jl}\) を施肥量と年次の交互作用とする。ブロックは年次ごとに異なる配置になっているので、\((\rho \gamma)_{kl}\) のようにブロックと年次の交互作用として取り出すことにした。また、密度×施肥量×年次の3次交互作用は解釈が難しいので取り出さなかった。残りの \(e_{ijkl}\) を誤差とすると、形質値データ \(y_{ijkl}\) の構造モデルは、 \[ y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma_l + (\alpha \beta)_{ij} + (\alpha \gamma)_{il} + (\beta \gamma)_{jl} + (\rho \gamma)_{kl} + e_{ijkl} \] となる。形質値として収量を考えたとき、4 年次にわたるデータの分散分析を行ってみよう。
dense <- factor(rice$density) # 水準のラベル化
fert <- factor(rice$fert)
blk <- factor(rice$rep)
year <- factor(rice$year)
y <- rice$gy # 収量データ
res <- aov(y ~ dense + fert + year + dense:fert + fert:year
+ dense:year + year:blk)
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## dense 1 125 125 0.033 0.85626
## fert 2 67405 33703 9.056 0.00104 **
## year 3 54399 18133 4.872 0.00806 **
## dense:fert 2 2352 1176 0.316 0.73186
## fert:year 6 75051 12508 3.361 0.01379 *
## dense:year 3 2376 792 0.213 0.88660
## year:blk 4 21144 5286 1.420 0.25517
## Residuals 26 96762 3722
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1分散分析の結果、年次効果と施肥量効果は主効果と交互作用ともに有意 となったが、密度に関しては有意とはならなかった。どのような効果が存在するかをみるため、有意な交互作用が認められた年次と施肥量水準の 組み合わせごとの応答を、図にプロットする。
x <- tapply(y, fert:year, mean); x## 1:1996 1:1997 1:1999 1:2000 2:1996 2:1997 2:1999 2:2000
## 324.8836 449.0768 370.3308 375.6779 497.6580 518.5043 395.1045 452.4563
## 3:1996 3:1997 3:1999 3:2000
## 400.6401 440.5633 407.0041 555.3635plot(1:4, x[1:4], type="b", xaxt="n", xlab="year", ylab="yield",
ylim=range(y), pch=0, col="red")
axis(1, 1:4, labels=levels(year))
points(1:4, x[5:8], type="b", pch=2, col="blue")
points(1:4, x[9:12], type="b", pch=3, col="green")
legend("bottomright", legend=c("fert1","fert2","fert3"), pch=c(0,2,3),
col=c("red","blue","green"))
これをみると、年次ごとに収量に変動があることがわかる。また、 無肥料区 (fert1) と少肥料区 (fert2) は 年次ごとの変動が似ていて、少肥料区の方が収量が高くなっているが、 多肥料区 (fert3) は 2000 年のように高収量のときもあるが、1996 年や 1997 年のようにそれほど高い収量が得られない場合もあり、これが交互作用として検出されている。これは、多肥料区で繁茂しすぎて倒伏が起こった年があったためと考えられる。
枝分かれ配置
各標本は、いくつかの副標本から構成されており、各副標本からまた次の副標本が取られる、といったような構造をもつ標本を枝分かれ配置 (nested design) という。例えば、ある果樹品種の甘さなどの形質を比較する場合、品種ごとに複数の個体（果樹）を選び、選ばれた個体から複数の果実を選び、形質を測定することが行われる。また、オス牛の肉質に与える効果を調べるため、オス牛を複数のメス牛と交配させ、そのメス牛から生まれた子牛の肉質を調べて、オス牛の評価を行う場合などである。
果樹品種やオス牛などの比較したい形質を、主効果 \(\alpha_i\) とし、その主効果ごとの副標本（果樹やメス牛など）を \(\beta_{ij}\)、形質が計測される個体の形質値（果実の甘さや子牛の肉質など）を \(y_{ijk}\) とおくと、データは以下のようなモデルとなる。 \[ y_{ijk} = \mu + \alpha_i + \beta_{ij} + e_{ijk}, \ i=1,\ldots,a, \ j=1,\ldots,b, \ k=1,\ldots,n \\ \sum_i \alpha_i = 0, \ \beta_{ij} \sim N(0, \sigma_B^2), \ e_{ijk} \sim N(0, \sigma^2) \]
ここで、注意しなければならないことは、興味の対象は主効果であり、副標本の効果についてはこのモデルでは主効果内でしか検出できないので、その効果はあくまで主効果に与える変動の大きさを測るだけである。牛の例で言えば、メス牛の効果を正しく測るためには、比較したいオス牛とメス牛の交配をすべての組合せで行い、2 元配置モデルの実験を行えば可能である。しかし、メス牛は交配できるオスの数が限られるため、多数の牛の比較はできない。しかし、オス牛は人工交配を行えば、ほぼ無限のメス牛と交配できるので、子牛の形質からオス牛の評価は枝分かれ配置モデルで簡単に行える。このため、\(\beta_{ij}\) は固定効果でなく、変量効果にして、その分散 \(\sigma_B^2\) を仮定するのが普通である。また、場合により、主効果も変量効果にすることもある。
一般に、枝分かれ配置では、副標本の個数 (\(b\) や \(n\)) は等しくないことが多いが、取り扱いが面倒になるので、この講義では取り上げず、釣り合い型の場合のみを取り上げる。 このモデルの帰無仮説は、 \[ {\rm H_A} : \alpha_1 = \cdots = \alpha_a = 0, \ {\rm H_B} : \sigma_B^2 = 0 \] である。分散分析表は、
枝分かれ配置分散分析表
| 要因 | 自由度 | 平方和 | 平均平方 | 平均平方の期待値 | \(F\) 値 |
|---|---|---|---|---|---|
| 主効果 | \(a-1\) | \(S_A\) | \(M_A\) | \(\sigma^2+n\sigma_B^2 + bn\sum\alpha_i^2/(a-1)\) | \(M_{A}/M_{AB}\) |
| 副標本 | \(a(b-1)\) | \(S_{AB}\) | \(M_{AB}\) | \(\sigma^2+n\sigma_B^2\) | \(M_{AB}/M_E\) |
| 誤差 | \(ab(n-1)\) | \(S_E\) | \(M_E\) | \(\sigma^2\) | |
主効果の検定は、主効果平均平方を副標本平均平方で割った \(F\) 値を用いなければならない。主効果平均平方の期待値に誤差分散と副標本分散の両方が含まれるからである。これは、分割区法で、1 次因子の効果を 1 次誤差で検定したのと同じである。
かぶら菜の 4 品種から無作為に 3 個体が選ばれ、各個体から 2 枚の葉を選び、そのカルシウム含量を計測したデータを解析してみよう。副標本（葉）の割り付けには2通りが考えられるが、「葉」はすべて異なるので、leaf2のような割り付けの方が良いかも知れない。
# カルシウム含量データ
cal <- c(3.28, 3.09, 3.52, 3.48, 2.88, 2.80, 2.46, 2.44, 1.87, 1.92, 2.19, 2.19,
2.77, 2.66, 3.74, 3.44, 2.55, 2.55, 3.78, 3.87, 4.07, 4.12, 3.31, 3.31)
leaf <- factor(rep(rep(1:3, each=2), 4)) # 個体ごとの葉の番号付け
leaf2 <- factor(rep(1:12, each=2)) # 葉の番号付け、すべての葉が異なる。
variety <- factor(rep(1:4, each=6))
# leaf の割り付けでは品種ごとの葉なので、variety/leafと書く必要がある。
av <- summary(aov(cal ~ variety + variety/leaf)); av## Df Sum Sq Mean Sq F value Pr(>F)
## variety 3 7.56 2.5201 378.73 3.80e-12 ***
## variety:leaf 8 2.63 0.3288 49.41 5.09e-08 ***
## Residuals 12 0.08 0.0067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# leaf2 の割り付けは、葉がすべて異なることが分かるので普通に書いて良い。
av2 <- summary(aov(cal ~ variety + leaf2)); av2## Df Sum Sq Mean Sq F value Pr(>F)
## variety 3 7.56 2.5201 378.73 3.80e-12 ***
## leaf2 8 2.63 0.3288 49.41 5.09e-08 ***
## Residuals 12 0.08 0.0067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1R の分散分析は、\(F\)値が誤差平均平方との比で算出されるので、枝分かれ配置では主効果の検定は正しくない。正確には以下のようにする必要がある。
# 平均平方の値の取得
# av はリストなので、その1番目の行列の次元
dim(av[[1]])## [1] 3 5fv <- av[[1]][1,3]/av[[1]][2,3]; fv # F 値## [1] 7.665167df1 <- av[[1]][1,1]; df1 # 分子自由度## [1] 3df2 <- av[[1]][2,1]; df2 # 分母自由度## [1] 8pv <- 1 - pf(fv, df1, df2); pv # p 値## [1] 0.009725121これより、正確な分散分析表は以下のようになる。
枝分かれ配置の正確な分散分析表
| 要因 | 自由度 | 平方和 | 平均平方 | \(F\) 値 | \(p\) 値 |
|---|---|---|---|---|---|
| 品種 | 3 | 7.56 | 2.5201 | 7.6652 | 0.009725 ** |
| 品種内個体 | 8 | 2.63 | 0.3288 | 49.41 | 5.09e-08 *** |
| 誤差 | 12 | 0.08 | 0.0067 | ||
このようにしなければならないのは、以下の例を考えればわかる。各個体から 2 枚の葉が選ばれ、それぞれ、カルシウム含量を計測したデータから解析したが、2 枚の葉を平均したデータを個体データとして解析する場合もある。
# 2枚の葉を平均したデータで分散分析
cal2 <- tapply(cal, leaf2, mean)
id2 <- factor(rep(1:4, each=3))
summary(aov(cal2 ~ id2))## Df Sum Sq Mean Sq F value Pr(>F)
## id2 3 3.780 1.2601 7.665 0.00973 **
## Residuals 8 1.315 0.1644
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1品種効果の \(F\) 値や \(p\) 値は、正しい枝分かれ配置の結果と同様になり、R の分散分析表は枝分かれ配置では使えないことがわかる。主効果の検定は葉を平均したデータでも同じであるなら、枝分かれ配置で解析する必要はどこにあるのだろうか。これは、個体内の葉の変動と品種内個体の変動とが比較できることである。この例では品種内個体の効果は 0.1％有意となっている。これは、同じ品種でも個体内でかなりの変動がある（品種間の変動程ではないが）ことを示している。枝分かれ配置を行うことで、品種内個体の変動の有無が解析できるのである。
追加
Rの分散分析では、検定の誤差をError()で指定できるので、主効果の検定は以下のようにすれば良い。
summary(aov(cal ~ variety + Error(variety:leaf)))## Warning in aov(cal ~ variety + Error(variety:leaf)): Error() モデルは特異です##
## Error: variety:leaf
## Df Sum Sq Mean Sq F value Pr(>F)
## variety 3 7.56 2.5201 7.665 0.00973 **
## Residuals 8 2.63 0.3288
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 12 0.07985 0.006654summary(aov(cal ~ variety + Error(leaf2)))##
## Error: leaf2
## Df Sum Sq Mean Sq F value Pr(>F)
## variety 3 7.56 2.5201 7.665 0.00973 **
## Residuals 8 2.63 0.3288
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 12 0.07985 0.006654主効果の検定はこれで行えるが、副標本の検定はこれでは行えないので通常のRの分散分析で行う。