必要なパッケージ

require(tidyverse)

実験計画法

因子と水準

　作物において、収量や抵抗性などの調査研究したい特性を形質(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品種の収量を調べる試験を行うことを考えてみる。

圃場配置 \[ \begin{matrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ \end{matrix} \]

しかし、この圃場の地力は均一ではなく、右側にいくほど収量が高くなる傾向があり、下側も収量が多少高くなっているとする。

圃場の地力（プロット1に対する相対値） \[ \begin{matrix} 0 & +1 & +2 & +3 \\ +0.2 & +1.2 & +2.2 & +3.2 \\ \end{matrix} \]

　いま、この圃場に上述したイネ4品種を2回繰り返しで栽培し、収量調査を行うことにした。ここで良くない配置の例としては、管理の容易さから品種の番号順に並べた

圃場作付け（番号順） \[ \begin{matrix} \rm{A} & \rm{B} & \rm{C} & \rm{D} \\ \rm{A} & \rm{B} & \rm{C} & \rm{D} \\ \end{matrix} \]

である。この配置では品種Aは最も地力の低い場所に栽培され、品種Dは最も地力の高い場所に栽培される。このため、品種の収量の大小と地力の大小が交絡 (confound) し、品種間の収量の差がわからなくなる。

　このような圃場でも、無作為化を導入することで系統的な地力の差を偶然誤差に変えることができる。例えば、圃場の8区画に無作為（random）に品種を割り付ける方法は、完全無作為化法 (completely randomized design) とよばれる。

　いま、\(i\)番目の品種の\(j\)番目の繰り返しの収量を \(y_{ij}\) とする。ここで、総平均を \(\mu\)、\(i\) 番目の品種の効果を \(\alpha_i\)（ただし \(\sum_i \alpha_i = 0\)）、誤差を \(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) \] とおける。以下に配置例を示す。

圃場作付け（完全無作為配置の例） \[ \begin{matrix} \rm{C} & \rm{D} & \rm{C} & \rm{B} \\ \rm{B} & \rm{D} & \rm{A} & \rm{A} \\ \end{matrix} \]

この例では、品種Aはたまたま比較的地力の高い場所に偏って配置されている。

　ここで、誤差の標準偏差を\(σ= 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 -2

地力

fert1 <- 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.2

誤差（正規乱数、標準偏差 0.5）

set.seed(1234)          # 乱数を固定するためのシード
error <- rnorm(8, sd=0.5); error    # 正規乱数誤差
## [1] -0.6035329  0.1387146  0.5422206 -1.1728489  0.2145623  0.2530279 -0.2873700
## [8] -0.2733159

　 　ここで、良くない配置例は、

x.reg <- factor(rep(1:4, 2)); x.reg     # 良くない配置例
## [1] 1 2 3 4 1 2 3 4
## Levels: 1 2 3 4

であり、完全無作為化法による配置例は、たとえば、

set.seed(123)
x.crd <- factor(sample(rep(1:4, 2))); x.crd 
## [1] 3 4 3 2 2 4 1 1
## Levels: 1 2 3 4
# 完全無作為化法による配置例. sampleが入っていることに注意

となる（上の表と同じ）。

　配置例\(x\)（x.reg, x.crd）が決まると、品種の収量は以下のように推定される。

rice <- yield[x.reg]            # 区画ごとの期待コメ収量。規則
y <- rice + land + error; y     # 観測される収量のシミュレート値
## [1] 59.39647 62.13871 60.54222 58.82715 60.41456 62.45303 59.91263 59.92668
gm <- mean(y); gm       # 総平均
## [1] 60.45143
vm <- tapply(y, x.reg, mean)        # 品種ごとの平均（yをxごとにmeanする）
est.vef <- vm - gm; est.vef     # 推定された品種効果
##          1          2          3          4 
## -0.5459175  1.8444391 -0.2240069 -1.0745146

分散分析の結果は、

summary(aov(y ~ x.reg))     # 分散分析結果
##             Df Sum Sq Mean Sq F value Pr(>F)  
## x.reg        3  9.809   3.270   9.545  0.027 *
## Residuals    4  1.370   0.343                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rice <- yield[x.crd]            # 区画ごとの期待コメ収量。規則
y <- rice + land + error; y     # 観測される収量のシミュレート値
## [1] 57.39647 58.13871 60.54222 62.82715 61.41456 58.45303 61.91263 62.92668
gm <- mean(y); gm       # 総平均
## [1] 60.45143
vm <- tapply(y, x.crd, mean)        # 品種ごとの平均（yをxごとにmeanする）
est.vef <- vm - gm; est.vef     # 推定された品種効果
##         1         2         3         4 
##  1.968225  1.669425 -1.482088 -2.155561
summary(aov(y ~ x.crd))     # 分散分析結果
##             Df Sum Sq Mean Sq F value Pr(>F)  
## x.crd        3 27.008   9.003   5.532  0.066 .
## Residuals    4  6.509   1.627                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

　推定された品種効果は以下のようになった。

品種 A B C D
真の品種効果 1 2 -1 -2
良くない配置 -0.546 1.844 -0.224 -1.075
完全無作為 1.968 1.669 -1.482 -2.156

これをみると、良くない配置では品種Aの効果は過少評価され、品種 C、D は過大評価されたことがわかる。これは、地力の差異を品種効果と見間違えたためである。一方、完全無作為法ではある程度正確な効果推定がなされたが、分散分析の結果は有意ではなかった。

1要因乱塊法

　完全無作為化法では品種Aのように、同じ品種の繰り返しが地力の高い区画に偏って配置されることがある。そこで、隣接した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） \[ \begin{matrix} ブロック1 & \rm{B} & \rm{D} & \rm{C} & \rm{A} \\ & - & - & - & - \\ ブロック2 & \rm{C} & \rm{B} & \rm{A} & \rm{D} \\ \end{matrix} \]

乱塊法配置例（パターン2） \[ \begin{matrix} 　& ブロック1 & & ブロック2 & \\ \rm{B} & \rm{D} & | & \rm{C} & \rm{B} \\ \rm{C} & \rm{A} & | & \rm{A} & \rm{D} \\ \end{matrix} \] となった。パターン1では行ごとにブロックを形成したため、ブロック内の地力は一様とは言えず、品種 A が地力の比較的高い場所に配置され、品種B は地力の比較的低い場所に配置されて偏りは解消されていない。 一方、パターン2は左右でブロックを形成したので、パターン1に比べるとブロック内の地力の差異は少なくなっている。さて、先ほどの例と同様に、上記の配置から各品種の収量がどのように推定されるかをみてみよう。

set.seed(12)
a1 <- sample(1:4); a1       # ブロック1での並べ替え
## [1] 2 4 3 1
a2 <- sample(1:4); a2       # ブロック2での並べ替え
## [1] 3 2 1 4

　パターン1の配置は、

x.rbd1 <- factor(c(a1,a2)); x.rbd1      # パターン1の配置
## [1] 2 4 3 1 3 2 1 4
## Levels: 1 2 3 4
bl1 <- factor(rep(1:2, each = 4)); bl1          # パターン1のブロック
## [1] 1 1 1 1 2 2 2 2
## Levels: 1 2

と定義され、パターン2の配置は、

x.rbd2 <- factor(c(a1[1:2],a2[1:2],a1[3:4],a2[3:4]))    ; x.rbd2    # パターン2の配置
## [1] 2 4 3 2 3 1 1 4
## Levels: 1 2 3 4
bl2 <- factor(rep(rep(1:2, each=2), 2)); bl2            # パターン2のブロック
## [1] 1 1 2 2 1 1 2 2
## Levels: 1 2

と定義される。

　配置\(x\)（x.rbd1, x.rbd2）とブロック定義\(bl\)（bl1, bl2）が決まると、品種の収量は以下のように推定される。

まず、パターン1について計算する。

rice <- yield[x.rbd1]           # 区画ごとの期待コメ収量
y <- rice + land + error; y     # 観測される収量のシミュレート値
## [1] 60.39647 58.13871 60.54222 61.82715 58.41456 62.45303 61.91263 59.92668
gm <- mean(y); gm       # 総平均
## [1] 60.45143
vm <- tapply(y, x.rbd1, mean)       # 品種ごとの平均
est.vef <- vm - gm; est.vef     # 推定された品種効果
##          1          2          3          4 
##  1.4184584  0.9733153 -0.9730408 -1.4187329
blm <- tapply(y, bl1, mean) - gm; blm   # ブロック効果
##          1          2 
## -0.2252939  0.2252939

分散分析結果は以下のようになる。

summary(aov(y ~ x.rbd1 + bl1))      # 分散分析結果
##             Df Sum Sq Mean Sq F value Pr(>F)
## x.rbd1       3 11.838   3.946   2.124  0.276
## bl1          1  0.406   0.406   0.219  0.672
## Residuals    3  5.574   1.858

次に、パターン2について計算する。

rice <- yield[x.rbd2]           # 区画ごとの期待コメ収量
y <- rice + land + error; y     # 観測される収量のシミュレート値
## [1] 60.39647 58.13871 60.54222 62.82715 58.41456 61.45303 61.91263 59.92668
gm <- mean(y); gm       # 総平均
## [1] 60.45143
vm <- tapply(y, x.rbd2, mean)       # 品種ごとの平均
est.vef <- vm - gm; est.vef     # 推定された品種効果
##          1          2          3          4 
##  1.2313967  1.1603769 -0.9730408 -1.4187329
blm <- tapply(y, bl2, mean) - gm; blm   # ブロック効果
##          1          2 
## -0.8507392  0.8507392

分散分析結果は以下のようになる。

summary(aov(y ~ x.rbd2 + bl2))      # 分散分析結果
##             Df Sum Sq Mean Sq F value Pr(>F)  
## x.rbd2       3 11.645   3.882   10.29 0.0435 *
## bl2          1  5.790   5.790   15.35 0.0296 *
## Residuals    3  1.132   0.377                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

推定された品種効果とブロック効果は以下のようになった。

品種 A B C D ブロック 1 ブロック 2
真の品種効果 1 2 -1 -2
パターン1 1.418 0.973 -0.973 -1.419 -0.225 0.225
パターン2 1.231 1.160 -0.973 -1.419 -0.851 0.851

これをみると、パターン1では品種Aは地力の高い区画に割り付けられたので過大評価されているが、パターン2では品種の効果が精度良く推定されているのがわかる。これは、パターン1のブロックは地力の差異を肩代わりしていないが、パターン2では地力の差異をブロックの差異にある程度置き換えることができたからである。

　両者を分散分析表で比較してみよう。分散分析の帰無仮説は、 \[ {\rm H_T} : \alpha_1 = \cdots = \alpha_4 = 0, \ {\rm H_B} : \sigma_B^2 = 0 \] で、品種効果は存在しないとブロック効果は無い、の２つの帰無仮説を同時に検定する。

パターン1
要因 自由度 平方和 平均平方 F値 p値
品種効果 3 11.838 3.946 2.124 0.276
ブロック効果 1 0.406 0.406 0.219 0.672
誤差 3 5.574 1.858
パターン2
要因 自由度 平方和 平均平方 F値 p値
品種効果 3 11.645 3.882 10.29 0.0435
ブロック効果 1 5.790 5.790 15.35 0.0296
誤差 3 1.132 0.377

パターン1では、地力の差がブロック効果として吸収されておらず、誤差平方和が大きいままであり、品種効果が検出されなかった。一方、パターン2では、地力の差がうまくブロック効果として吸収され、誤差平方和が小さくなり、品種効果を検出することができた。

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反復の場合の試験配置例を以下に示す。 \[ \begin{matrix} 　ブロック1 &&& ブロック2 &&& ブロック3 \\ \rm{S0} & \rm{D1} & | & \rm{S2} & \rm{D0} & | & \rm{D1} & \rm{S2} \\ \rm{S2} & \rm{D2} & | & \rm{D1} & \rm{S0} & | & \rm{S0} & \rm{S1} \\ \rm{D0} & \rm{S1} & | & \rm{S1} & \rm{D2} & | & \rm{D2} & \rm{D0} \\ \end{matrix} \]

　以下は、東京大学農学部における学生実習として行われた水稲栽培の圃場試験データを解析する。試験は、1996年から2000年( 　1998年は除く）にわたり行われ、その目的は、品種、施肥量、栽植密度が水稲の収量構成要素に与える影響を評価することであった。今回の演習では、このうち「日本晴」品種に注目して解析する。施肥量と栽植密度の2要因実験が2反復で行われた。ここでは、1999 年度のデータで対象形質を玄米千粒重 (gw) として解析してみる。

データを読み込んで確認する。

rice <- read.csv("data/ricecul.csv")  # 日本晴データ読み込み
head(rice)
##   year plot variety density fert rep    pn        gn        pr       gw
## 1 2000    2       1       1    1   1 187.2  95.57487 0.7393026 26.39719
## 2 2000   21       1       1    1   2 194.4 136.23598 0.7772188 24.98339
## 3 2000    8       1       1    2   1 216.0 118.51032 0.8497141 25.42373
## 4 2000   19       1       1    2   2 255.6 117.29109 0.8420248 25.63140
## 5 2000    6       1       1    3   1 306.0  93.27377 0.8104506 25.99653
## 6 2000   15       1       1    3   2 331.2 100.35559 0.7711600 25.59727
##         gy
## 1 307.7816
## 2 384.9432
## 3 446.6109
## 4 464.5329
## 5 565.1520
## 6 532.7308

tidyverseパッケージの関数を使ってデータを整形する。

df <- rice %>% filter(year == 1999) %>% 
  mutate(plot = factor(plot), density = factor(density), fert = factor(fert), rep = factor(rep)) %>%
  mutate(y = gw) %>%
  select(density, fert, rep, y)
head(df)
##   density fert rep        y
## 1       1    1   1 26.19850
## 2       1    1   2 25.84801
## 3       1    2   1 26.65339
## 4       1    2   2 27.12352
## 5       1    3   1 26.74632
## 6       1    3   2 26.59874
summary(aov(y ~ density * fert + rep, data = df))   # ２要因乱塊法分散分析
##              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

　なお、この例では分散分析により交互作用が検出された。これは、以下に示す密度、施肥量組合せ平均のグラフでわかるように、密植で肥料水準が少のときだけなぜか千粒重が低下したためである。

dfm <- tapply(df$y, df$density:df$fert, mean); dfm　# 密度、施肥量組合せ平均を求める
##      1:1      1:2      1:3      2:1      2:2      2:3 
## 26.02326 26.88846 26.67253 25.75824 22.38495 26.80629
plot(dfm[1:3], type="b", xaxt="n", ylim=c(20,30), xlab="fertility", ylab="gw")
lines(dfm[4:6], type="b", col="red")
axis(1, 1:3)
legend("bottomright", legend=c("sparse","dense"), pch=1, col=c("black", "red"))

分割区法(split-plot design)

　潅水区や非潅水区のような水管理は、比較的広い区画で行わないと労力、費用等がかかる。このような処理を用いた実験の場合、乱塊法でブロック内すべての組み合わせの試験区を設定することは現実的でない。水管理のように広い区画が必要な因子を1次因子、その区画内で水準を変化させる品種や施肥量などの因子を2次因子とよぶ。

先ほどの乱塊法の例で、栽植密度を潅水に変えたものとしてみてみよう。潅水区、非潅水区 (S, D) を 3 反復の分割区法に割り付けたとすると、例えば以下のような配置となる。

\[ \begin{matrix} 　ブロック1 &&& ブロック2 &&& ブロック3 \\ \rm{S0} & \rm{D1} & | & \rm{S2} & \rm{D0} & | & \rm{D1} & \rm{S2} \\ \rm{S2} & \rm{D2} & | & \rm{S0} & \rm{D1} & | & \rm{D0} & \rm{S1} \\ \rm{S1} & \rm{D0} & | & \rm{S1} & \rm{D2} & | & \rm{D2} & \rm{S0} \\ \end{matrix} \]

　1次因子（水管理）は1要因乱塊法としてブロックごとに無作為に配置され、1次因子の配置の中に2次因子が無作為に配置される。この例では、ブロック1とブロック2では右側に、ブロック3では左側にSが無作為に割り当てられており、その中に施肥条件が無索に配置されている。

　いま、\(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(df$y, df$density:df$rep, mean); drm   # 潅水（実際は密度）とブロックの組合せ平均
##      1:1      1:2      2:1      2:2 
## 26.53273 26.52343 24.07773 25.88859
yy <- drm[df$density:df$rep]; yy # yの変動のうち、組合せ平均によるものをyyとする
##      1:1      1:2      1:1      1:2      1:1      1:2      2:1      2:2 
## 26.53273 26.52343 26.53273 26.52343 26.53273 26.52343 24.07773 25.88859 
##      2:1      2:2      2:1      2:2 
## 24.07773 25.88859 24.07773 25.88859

1次因子平方和（組合せ平均を密度、ブロック、両者の交互作用の効果に分割）

df2 <- data.frame(yy, df)
summary(aov(yy ~ density + rep, data = df2))
##             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 ' ' 1

2次因子（組合せ平均を除いた変動を、施肥と、密度と施肥の交互作用に分割）

summary(aov(y - yy ~ fert + density:fert, data = df2))
##              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 ' ' 1

分割区法の分散分析を行うには、以下のようにする。

summary(aov(y ~ density * fert + Error(rep/density), data = df))    # 分割区法分散分析
## 
## 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 年次にわたるデータの分散分析を行ってみよう。まずはデータを準備する。

df <- rice %>% 
  mutate(year = factor(year), density = factor(density), fert = factor(fert), rep = factor(rep)) %>%
  mutate(y = gy) %>%
  select(year, density, fert, rep, y)
head(df)
##   year density fert rep        y
## 1 2000       1    1   1 307.7816
## 2 2000       1    1   2 384.9432
## 3 2000       1    2   1 446.6109
## 4 2000       1    2   2 464.5329
## 5 2000       1    3   1 565.1520
## 6 2000       1    3   2 532.7308

分散分析を行う。

res <- aov(y ~ density + fert + year 
           + density:fert + fert:year + density:year + year:rep, data = df)
summary(res)
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## density       1    125     125   0.033 0.85626   
## fert          2  67405   33703   9.056 0.00104 **
## year          3  54399   18133   4.872 0.00806 **
## density:fert  2   2352    1176   0.316 0.73186   
## fert:year     6  75051   12508   3.361 0.01379 * 
## density:year  3   2376     792   0.213 0.88660   
## year:rep      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(df$y, df$fert:df$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.3635
plot(1:4, x[1:4], type="b", xaxt="n", xlab="year", ylab="yield",
       ylim=range(df$y), pch=0, col="red")
axis(1, 1:4, labels=levels(df$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) という。例えば、ある果樹品種の甘さなどの形質を比較する場合、品種ごとに複数の個体（果樹）を選び、選ばれた個体から複数の果実を選び、形質を測定することが行われる。今、果樹の比較したい形質について、\(i\)番目の品種の効果を主効果 \(\alpha_i\) とし、品種ごとに得られる\(j\)番目の個体（主効果ごとの副標本）の効果を\(\beta_{ij}\)、個体ごとに収穫される\(k\)番目の果実を計測される形質値（果実の甘さなど）を\(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) \]

　ここで、注意しなければならないことは、興味の対象は主効果であり、副標本の効果についてはこのモデルでは主効果内でしか検出できないので、その効果はあくまで主効果に与える変動の大きさを測るだけである。果実の例で言えば、品種ごとに個体は異なり、品種間で共通した個体を得ることはできない。そこで、通常は、\(\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のような割り付け方もできる。

# カルシウム含量データ
y <- 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))
variety <- factor(rep(1:4, each=6))
df <- data.frame(y, leaf, variety)
head(df)
##      y leaf variety
## 1 3.28    1       1
## 2 3.09    1       1
## 3 3.52    2       1
## 4 3.48    2       1
## 5 2.88    3       1
## 6 2.80    3       1

枝分かれ分散分析を行う。

# 枝分かれ分散分析では、主標本 / 副標本　としてモデルを指定する
av <- summary(aov(y ~ variety / leaf, data = df)); 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

ここで注意しなければならないのは、上に示した分散分析では、\(F\)値が誤差平均平方との比で算出されるので、枝分かれ配置では主効果の検定が正しく行われないことである。正確に枝分かれ分散分析を行うには以下のように主効果の検定を行う必要がある。

dim(av[[1]]) # 分散分析表の数値が3 x 5の行列として保存されている
## [1] 3 5
fv <- av[[1]][1,3]/av[[1]][2,3]; fv    # F 値
## [1] 7.665167
df1 <- av[[1]][1,1]; df1     # 分子自由度
## [1] 3
df2 <- av[[1]][2,1]; df2     # 分母自由度
## [1] 8
pv <- 1 - pf(fv, df1, df2); pv   # p 値
## [1] 0.009725121

あるいは、分割区法と同様に、

summary(aov(y ~ variety + Error(variety:leaf), data = df))
## Warning in aov(y ~ variety + Error(variety:leaf), data = df): Error() model is
## singular
## 
## 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.006654

以上をまとめると、枝分かれ配置の分散分析表は以下のようになる。

枝分かれ配置の分散分析表
要因 自由度 平方和 平均平方 \(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枚の葉を平均したデータで分散分析
y2 <- tapply(y, variety:leaf, mean)
id2 <- factor(rep(1:4, each = 3))
summary(aov(y2 ~ 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\) 値は、正しい枝分かれ配置の結果と同様になり、通常の分散分析は枝分かれ配置では使えないことがわかる。主効果の検定は葉を平均したデータでも同じであるなら、枝分かれ配置で解析する必要はどこにあるのだろうか。これは、個体内の葉の変動と品種内の個体変動とが比較できることである。この例では品種内個体の効果は0.1％有意となっている。これは、同じ品種でも個体内でかなりの変動があることを示している。枝分かれ配置を行うことで、品種内の個体変動の有無が解析できるのである。

参考資料