必要なパッケージ
require(lme4)## Loading required package: lme4## Loading required package: Matrix1 一元配置 (One-way layout)
1.1 構造モデル
前回の授業で 2 標本（処理）間の平均の比較、すなわち、2つの正規母集団平均間の比較、を行った。 今回の講義ではその拡張として、3 つ以上の標本（処理）間の平均の比較を行う方法を考える。
\(k\) 個の独立した母集団 \(P_1, P_2, \ldots, P_k\) の母平均 \(\mu_1, \mu_2, \ldots, \mu_k\) が同一かどうか調べたい。 \(\mu_i\) が \(i\) 番目の処理 (treatment) を表すとき処理平均 (treatment mean) とよばれる。 ここで、\(i\) 番目の母集団 \(P_i\) から、大きさ \(n_i\) の標本、\(y_{i1}, y_{i2}, \ldots, y_{in_i}\) を得たとする。また、\(N = \sum_i^k n_i\) を全標本の大きさとする。
このとき、
- 全ての母集団が正規分布に従い、
- 全ての母集団の分散が同一で、
- 標本は互いに独立である。
すなわち、 \[ y_{ij} \overset{\textrm{iid}}{\sim} N(\mu_i, \sigma^2) \] であると仮定する。なお、iid は independently identically distributed (独立で同一に分布する) という意味である。 母平均 \(\mu_i\) が互いに等しいかどうかを検定する手法が、一元配置分散分析 (One-way Analysys of Variance (ANOVA)) である。すなわち、このときの帰無仮説は、 \[ {\rm H}_0:\mu_1 = \mu_2 = \cdots = \mu_k \] であり、対立仮説は帰無仮説のあらゆる否定になり、少なくとも1つの処理が他の処理と異なっている場合を含んでいる。
母集団平均 \(\mu_i\) について、\(\mu_i = \mu + \alpha_i\) とモデル化することもできる。ここで、\(\mu\) は総平均、\(\alpha_i\) は処理効果と考えられる（先の\(\mu_i\)は処理平均だった）。このとき、\(\alpha_i\) の一意性を確保するために、\(\sum_i \alpha_i = 0\) の制約を入れる。すると、標本 \(y_{ij}\) は、\(e_{ij}\) をすべての \(i、j\) で互いに独立に、平均 0、分散 \(\sigma^2\) である正規分布に従う誤差項 (error term) を考慮し、 \[ y_{ij} = \mu + \alpha_i + e_{ij}\\ e_{ij} \overset{\textrm{iid}}{\sim} N(0, \sigma^2), \ i=1, \ldots, k, \ j = 1, \ldots, n_i \] とモデル化できる。これが、一元配置分散分析の構造モデルであり、帰無仮説は、すべての処理効果が無いという \[ {\rm H}_0:\alpha_1 = \alpha_2 = \cdots = \alpha_k = 0 \] となる。
なお、標本の大きさがいずれの母集団で等しい、 \[ n_1 = n_2 = \cdots = n_k = n \] とき、分散分析は釣り合っている (balanced) という。多くの実験計画では、釣り合い型の分散分析がデザインされる。実験途中のアクシデントなどからいくつかの観測値が欠測した不釣り合い (unbalanced) デザインになることもあるが、一般には釣り合い型のデザインの方が良い。
1.2 平方和分解とその期待値
いま、\(\bar{y}_{i.}\) を処理 \(i\) における標本平均 (sample mean)、\(\bar{y}_{..}\) を標本総平均 (average of all observations) とする。すなわち、 \[ \bar{y}_{i.} = \frac{1}{n_i} \sum_j y_{ij}\\ \bar{y}_{..} = \frac{\sum_i \sum_j y_{ij}}{\sum_i n_i}=\frac{1}{N}\sum_i\sum_j y_{ij} \] と計算される。
標本全体がもつ変動の情報は、総平方和 \(S_T\) (total sum of suquares) で表されるが、これは、 \[ S_T = \sum_i \sum_j (y_{ij} - \bar{y}_{..})^2 \\ = \sum_i \sum_j \{(y_{ij} - \bar{y}_{i.}) + (\bar{y}_{i.} - \bar{y}_{..}) \}^2 \\ = \sum_i \sum_j(y_{ij} - \bar{y}_{i.})^2 + \sum_i \sum_j (\bar{y}_{i.} - \bar{y}_{..})^2 \\ = \sum_i \sum_j(y_{ij} - \bar{y}_{i.})^2 + \sum_i n_i (\bar{y}_{i.} - \bar{y}_{..})^2 \\ = S_E + S_A \] のように誤差平方和 \(S_E\)（error sum of squares）と処理平方和 \(S_A\)（treatment sum of squares）とに分解される。これを、平方和分解 という。なお、誤差平方和を群内平方和（within groups sum of squares），処理平方和を群間平方和（between groups sum of squares）と呼ぶこともある。
個々の標本 \(y_{ij}\)、処理平均 \(\bar{y}_{i.}\)、総平均 \(\bar{y}_{..}\) を構造モデルで表すと、 \[ y_{ij} = \mu + \alpha_i + e_{ij}, \ E[e_{ij}^2] = \sigma^2 \\ \bar{y}_{i.} = \mu + \alpha_i + \bar{e}_{i.}, \ E[\bar{e}_{i.}^2] = \frac{\sigma^2}{n_i} \\ \bar{y}_{..} = \mu + \bar{e}_{..}, \ E[\bar{e}_{..}^2] = \frac{\sigma^2}{N} \] となるので、誤差平方和と処理平方和の期待値は、それぞれ、 \[ E[S_E] = E[\sum_i^k \sum_j^{n_i} (y_{ij}-\bar{y}_{i.})^2] = E[\sum_i^k \sum_j^{n_i}(e_{ij}-\bar{e}_{i.})^2]= \\ \sum_i^k E[\sum_j^{n_i} e_{ij}^2 - 2\bar{e}_{i.}\sum_j^{n_i} e_{ij} + n_i \bar{e}_{i.}^2] = \sum_i^k \{ \sum_j^{n_i} E[e_{ij}^2] - n_i E[\bar{e}_{i.}^2] \}\\ = \sum_i^k \sum_j^{n_i} \sigma^2 - \sum_i^k n_i \frac{\sigma^2}{n_i} = (N-k)\sigma^2 \] \[ E[S_A] = E[\sum_i^k n_i (\bar{y}_{i.}-\bar{y}_{..})^2] = E[\sum_i^k n_i \{\alpha_i + \bar{e}_{i.} - \bar{e}_{..} \}^2] \\ = \sum_i^k n_i \alpha_i^2 + \sum_i^k n_i E[(\bar{e}_{i.}-\bar{e}_{..})^2] = \sum_i^k n_i \alpha_i^2 + \sum_i^k n_i(E[\bar{e}_{i.}^2] - E[\bar{e}_{..}^2]) \\ = \sum_i n_i \alpha_i^2 + \sum_i n_i \frac{\sigma^2}{n_i}-\sum_i n_i \frac{\sigma^2}{N} = \sum_i n_i \alpha_i^2 + (k-1)\sigma^2 \] となる。ここで誤差平方和の期待値、\(E[S_E] = (N-k)\sigma^2\) の \(N-k\) を 平方和 \(S_E\) の自由度という。なお、平方和を自由度で割った量の 期待値は\(E[S_E/(N-k)] = \sigma^2\)である、\(S_E/(N-k)\) は\(\sigma^2\) の 不偏推定量である。また、処理平方和の期待値において、もし すべての \(\alpha_i\) が 0 であるならば、\(E[S_A/(k-1)] = \sigma^2\) と なり \(S_A/(k-1)\) も \(\sigma^2\) の不偏推定量となる。よって、\(k-1\) が処理平方和 \(S_A\) の 自由度になる。平方和をその自由度で割った量を平均平方 (mean square) という。すなわち、 \[ 平均平方 = \frac{平方和}{自由度} \] である。
総平方和 \(S_T\) が処理平方和 \(S_A\) と誤差平方和 \(S_E\) に分解されたことに応じ、総平方和の自由度 \(N-1\) も処理平方和の自由度 \(k-1\) と誤差平方和の自由度 \(N-k\) に分解されている。
1.3 \(F\) 検定
帰無仮説の元での平方和の比の分布
一元配置モデルにおける帰無仮説は、すべての処理効果がない、つまり、 \[ {\rm H}_0 : \alpha_1 = \alpha_2 = \cdots = \alpha_k = 0 \] である。帰無仮説の元では、\(S_E/\sigma^2\) は自由度 \(N-k\) のカイ二乗分布に従い、\(E[S_A]\) における \(\sum_i n_i \alpha_i^2\) の項が 0 になるので \(S_A/\sigma^2\) は自由度 \(k-1\) の カイ二乗分布に従うことがわかり、両者は独立である。これらのカイ二乗分布に従う変量をその自由度で割った比である \(F\) は、 \[ S_E/\sigma^2 \sim \chi^2(N-k)\\ S_A/\sigma^2 \sim \chi^2(k-1) \\ F = \frac{(S_A/\sigma^2)/(k-1)}{(S_E/\sigma^2)/(N-k)} = \frac{S_A/(k-1)}{S_E/(N-k)} = M_A/M_E \sim F(k-1, N-k) \] のように自由度 \(k-1, N-k\) の \(F\) 分布に従う。なお、\(k-1\) を分子自由度、\(N-k\) を分母自由度と言う。 これより、検定統計量 \(F\) を用いて帰無仮説の検定を行うことができる。なお、\(M_A\) は処理平方和 \(S_A\) をその自由度 \(k-1\) で割った量で、処理平均平方（treatment mean square）と呼ばれる。一方、\(M_E\) は誤差平方和 \(S_E\) をその自由度 \(N-k\) で割った量で、誤差平均平方（error mean square）と呼ばれ、先ほど述べたように誤差分散 \(\sigma^2\) の不偏推定量である。
帰無仮説が真のときは、処理平均平方 \(M_A\) も誤差分散の不偏推定量になるので、\(M_E\) と \(M_A\) は似た値になり、両者の比である \(F\) 値は \(1\) に近い値になると考えられる。よって、\(F\) 値が \(1\) より有意に大きな値をとるのは、\(\sum_i n_i \alpha_i^2\) が 0 より有意に大きな値を取っていた、すなわち、すべての処理効果 \(\alpha_i=0\) という帰無仮説が疑わしいと考えられる。具体的には \(F\) 分布から計算される \(p\) 値で帰無仮説の検定が行える。
分散分析表
一元配置モデルの解析結果は、以下の分散分析表（ANOVA table）にまとめられる。
| 要因 | 自由度 | 平方和 | 平均平方 | \(F\) 値 | \(p\) 値 |
|---|---|---|---|---|---|
| 処理 | \(k-1\) | \(S_A\) | \(M_A=S_A/(k-1)\) | \(F = M_A/M_E\) | \({\bf P}[F_{k-1,\ N-k} > F]\) |
| 誤差 | \(N-k\) | \(S_E\) | \(M_E=S_E/(N-k)\) | ||
この表から検定統計量 \(F\) 値が求められる．そして， 自由度 \(k - 1，N - k\) の \(F\) 分布の 1 - \(\gamma\) 点（例えば 95 ％点）\(F_{k - 1, \ N - k, \ 1 - \gamma}\) より \(F\) 値が大きい，すなわち， \[ F > F_{k - 1, \ N - k, \ 1 - \gamma} \] とき帰無仮説を棄却すれば、有意水準 \(\gamma\) (たとえば 5 ％) の検定が行える。これを、\(F\) 検定（\(F\) test）という。
例1．凝固時間
Box ら (2005) による凝固時間の例で一元配置分散分析の解析を行ってみる。24 頭のある動物が無作為に 4 つの異なる食餌 (diets) が与えられた。食餌ごとの動物の数は等しくなかった。動物ごとに血液凝固時間 (times) が計測された。食餌の違いが凝固時間に影響を与えるか調べてみる。
# ex 1
times <- c(62, 60, 63, 59, 63, 67, 71, 64, 65, 66, 68, 66, 71, 67,
68, 68, 56, 62, 60, 61, 63, 64, 63, 59)
diets <- factor(c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,
3, 4, 4, 4, 4, 4, 4, 4, 4)) # 4つのグループを示すIDを因子（factor）とする。重要な操作！
table(diets) # 食餌ごとの標本の大きさ## diets
## 1 2 3 4
## 4 6 6 8boxplot(times ~ diets)
gm <- mean(times); gm # 総平均## [1] 64N <- length(times); N # データ数## [1] 24k <- length(levels(diets)); k # 処理数## [1] 4vm <- tapply(times, diets, mean); vm # 食餌ごとの平均（times を diets ごとに mean する）## 1 2 3 4
## 61 66 68 61est.vef <- vm - gm # 推定された食餌効果
vefv <- est.vef[diets] # 食餌効果ベクトル
errv <- times - gm - vefv # 誤差ベクトル
# データの分解
cbind(times, gm, vefv, errv)## times gm vefv errv
## 1 62 64 -3 1
## 1 60 64 -3 -1
## 1 63 64 -3 2
## 1 59 64 -3 -2
## 2 63 64 2 -3
## 2 67 64 2 1
## 2 71 64 2 5
## 2 64 64 2 -2
## 2 65 64 2 -1
## 2 66 64 2 0
## 3 68 64 4 0
## 3 66 64 4 -2
## 3 71 64 4 3
## 3 67 64 4 -1
## 3 68 64 4 0
## 3 68 64 4 0
## 4 56 64 -3 -5
## 4 62 64 -3 1
## 4 60 64 -3 -1
## 4 61 64 -3 0
## 4 63 64 -3 2
## 4 64 64 -3 3
## 4 63 64 -3 2
## 4 59 64 -3 -2# 平方和の計算
SA <- sum(vefv^2); SA # 処理平方和## [1] 228SE <- sum(errv^2); SE # 誤差平方和## [1] 112MA <- SA/(k-1); MA # 処理平均平方## [1] 76ME <- SE/(N-k); ME # 誤差平均平方## [1] 5.6fvalue <- MA/ME; fvalue # F 値 ## [1] 13.57143qf(0.95, k-1, N-k) # 自由度k-1, N-k F分布95％点## [1] 3.0983911 - pf(fvalue, k-1, N-k) # p 値## [1] 4.658471e-05#
av <- aov(times ~ diets) # 一元配置分散分析
summary(av) # 結果表示## Df Sum Sq Mean Sq F value Pr(>F)
## diets 3 228 76.0 13.57 4.66e-05 ***
## Residuals 20 112 5.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.4 等分散性（homoscedasticity）の検定
一元配置分散分析では、すべての処理に対して誤差分散は一定の \(\sigma^2\) である、という仮定がおかれるが、実際に等分散性の仮定が成り立っているかは気になるところである。一般に、多くの母集団平均の比較では、比較したい処理が増えるため処理ごとの標本の大きさは小さくなる傾向にある。このため、個々の処理ごとの分散の推定値は不安定になる。誤差分散が一定という仮定を置くことにより、分散推定にかかわるデータ数を増やすことで分散推定の不安定さを軽減させている、と考えることもできる。処理ごとの標本サイズがそれほど大きくないときには、等分散性にこだわりすぎなくともよい。
等分散性の検定はいくつか知られているが、ここではバートレット (Bartlett) 検定を紹介する。
# バートレット検定
bartlett.test(times ~ diets)##
## Bartlett test of homogeneity of variances
##
## data: times by diets
## Bartlett's K-squared = 1.668, df = 3, p-value = 0.6441凝固時間の例では、バートレット検定で等分散性の仮定は棄却されなかった。等分散性の検定が有意になったときは、そのままでは分散分析は使い難いので、対数変換等の変数変換で等分散性を確保するか、等分散性を仮定しない一元配置法、もしくは一元配置ノンパラメトリック検定である Kruskal-Wallis 検定を行う。
# 等分散性を仮定しない（var=F）一元配置
oneway.test(times ~ diets, var = F)##
## One-way analysis of means (not assuming equal variances)
##
## data: times and diets
## F = 16.728, num df = 3.0000, denom df = 9.9533, p-value = 0.0003249# 一元配置ノンパラメトリック検定
kruskal.test(times ~ diets)##
## Kruskal-Wallis rank sum test
##
## data: times by diets
## Kruskal-Wallis chi-squared = 17.015, df = 3, p-value = 0.00070161.5 母数効果と変量効果
3 つ以上の母集団平均 \(\mu_j\) 間の比較を行う一元配置分散分析では、平均間の差に興味があり、 \[ \mu_i = \mu + \alpha_i + e_{ij}\\ \sum_i \alpha_i = 0\\ e_{ij} \sim N(0, \sigma^2) \] のようにモデル化して、一意性確保のため \(\sum_i \alpha_i = 0\) の制約を加えた。このとき、解析の興味は処理効果 \(\alpha_i\) の正負、大小にあり、このとき、\(\alpha_i\)を母数（固定）効果 (fixed effect) とよぶ。帰無仮説は、すべての処理効果が無い、すなわち、 \[ {\rm H}_0 : \alpha_1 = \alpha_2 = \cdots = \alpha_k = 0 \] となる。
一方、たとえば、植物育種では、ある自殖系作物の \(k\) 系統ごとにそれぞれ \(n_i\) 区画の圃場で栽培試験を行ったときに、個々の系統の優劣より、その作物の遺伝的特性による変動の大きさ (遺伝分散) と非遺伝的要因による変動の大きさ (環境分散)とを比較したい場合がある。このようなとき、\(\alpha_i\) に誤差分布と独立な分布を仮定（具体的には正規分布）する。帰無仮説は、処理間にばらつきが存在しない、すなわち、 \[ y_{ij} = \mu + \alpha_i + e_{ij}\\ \alpha_i \sim N(0, \sigma_a^2)\\ \ e_{ij} \sim N(0, \sigma^2) \\ \] のようにモデル化し、帰無仮説は、処理により生じる分散が0である、すなわち、 \[ {\rm H}_0 : \sigma_a^2 = 0 \] とする。このとき\(\alpha_i\)を変量効果 (random effect) とよぶ。先の植物育種の例では \(\sigma_a^2\) が遺伝分散であり、\(\sigma^2\) が環境分散である。遺伝分散の全分散に対する比を遺伝率 (heritability) といい、遺伝率が高い形質は、それだけ育種に期待が持てることになるが、遺伝率が低い形質では、育種より栽培管理技術や計測技術の向上を目指すべきである、などの指針が得られる。
変量モデルにおいては、標本総平均の構造モデルは先の母数効果のモデルと若干異なり、 \[ y_{..} = \mu + \bar{\alpha} + \bar{e}_{..}\\ E[\bar{\alpha}^2] = \sigma_a^2/k\\ E[\bar{e}_{..}^2] = \sigma^2/N \] となり、処理平方和の期待値は、 \[ E[S_A] = E[\sum_i n_i (\bar{y}_{i.}-\bar{y}_{..})^2] = E[\sum_i n_i \{(\alpha_i - \bar{\alpha}) + (\bar{e}_{i.} - \bar{e}_{..}) \}^2] \\ = \sum_i n_iE[ (\alpha_i - \bar{\alpha})^2] + \sum_i n_i E[(\bar{e}_{i.}-\bar{e}_{..})^2] = \sum_i n_i (E [\alpha_i^2]-E[\bar{\alpha}^2]) + \sum_i n_i(E[\bar{e}_{i.}^2] - E[\bar{e}_{..}^2]) \\ = \sum_i n_i (\sigma_a^2-\sigma_a^2/k) + \sum_i n_i \frac{\sigma^2}{n_i}-\sum_i n_i \frac{\sigma^2}{N} = \frac{N}{k}(k-1) \sigma_a^2 + (k-1)\sigma^2 \] となる。帰無仮説が真ならば \(\sigma_a^2 = 0\) なので、母数効果のモデルと同様に分散分析表が得られ、分散比である \(F\) 値が自由度 \(k-1\), \(N-k\) である \(F\) 分布に従うことを利用して検定が行える。
帰無仮説が棄却されたときは、処理効果の分散 \(\sigma_a^2\) がゼロでないと考えられ、 \[ E[M_A] = E[S_A/(k-1)] = \frac{N}{k} \sigma_a^2 + \sigma^2 \] なので、 \[ \hat{\sigma_a^2} = \frac{k}{N}(M_A - \hat{\sigma^2}) = \frac{k}{N}(M_A - M_E) \] と推定できる。
例えば、例1について、食餌(diets)が様々な食餌の中から無作為に選ばれたもので、個々の食餌の効果には興味がなく、食餌が血液凝固時間に与える変動の大きさに興味があるとする。このとき、食餌を変量効果とすると、その変量効果がしたがう正規分布の分散（分散成分, variance component）の推定値は、
sa <- k / N * (MA - ME); sa # 食餌の分散成分## [1] 11.73333である。なお、誤差の分散の推定値は、
ME## [1] 5.6である。
変量モデルや混合モデル（母数効果と変量効果の両方を含むモデル）は、lme4パッケージのlmer関数などを用いて適用できる。
model <- lme4::lmer(times ~ (1 | diets))
summary(model)## Linear mixed model fit by REML ['lmerMod']
## Formula: times ~ (1 | diets)
##
## REML criterion at convergence: 115.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.18491 -0.59921 0.09332 0.54078 2.17508
##
## Random effects:
## Groups Name Variance Std.Dev.
## diets (Intercept) 11.692 3.419
## Residual 5.599 2.366
## Number of obs: 24, groups: diets, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 64.01 1.78 35.96変量効果の分散が、上で計算した値に一致していることが確認できる。
1.6 多重比較 (\(\alpha_i =0\))
分散分析は、\(k\) 処理に対して少なくともどれか1つの処理に効果があるかを示すだけで、 どの処理間に差があるかについては、何ら情報を与えていない。 \(k\) 処理に対して、どの処理間に違いがあるかを調べるにはすべての処理間での対比較を 行う必要がある。すなわち、\({}_k {\rm C}_2 = k(k-1)/2\) 通りの検定を行うことになる。
ここで、\(k(k-1)/2\) 通りの対比較を通常の2標本\(t\)検定で行うと、同時比較に対する 有意水準を制御するのが難しい。対比較は互いに独立な検定ではないが、たとえば、有意 水準 5 %の独立な検定を 100 回行えば、効果が無いにもかかわらず5回は有意な組合せ、 すなわち、偽陽性 (false positive) が出てしまうからである。これを多重比較の問題という。 R ではこの問題に対処するためにいくつかの手法が実装されている。
ボンフェローニ (Bonferroni) 補正
いま、有意水準 \(\alpha'\) のそれぞれ独立な検定を \(r\) 回行ったとすると、1 回の検定で正しい判断を行う確率が \(1 - \alpha'\) なので、\(r\) 回の検定で正しい判断を行う確率は、\((1 - \alpha')^r\) となる。これより正しい判断を行わない（第 1 種の過誤の）確率は、 \[ 1 - (1 - \alpha')^r \approx 1 - (1 - r\alpha') = r\alpha', \ ただし, \ \alpha' \approx 0 \] となる。これが、\(r\) 回の同時検定における有意水準となる。よって、検定全体での有意水準を \(\alpha\) にするには， 1 回の検定の有意水準を \(\alpha' = \alpha/r\) にすればよいこれ。がボンフェローニ補正である。しかし、 多重比較における検定は独立な検定ではないので、この補正は厳しすぎて、有意な組み合わせを見逃す恐れがある。
なお、R の多重比較では、補正なしの \(p\) 値を \(r\) 倍した \(p\) 値を出力する。ただし、これが 1 を超えた場合は 1 とする。
pairwise.t.test(times, diets, p.adj = "bonf")##
## Pairwise comparisons using t tests with pooled SD
##
## data: times and diets
##
## 1 2 3
## 2 0.02282 - -
## 3 0.00108 0.95266 -
## 4 1.00000 0.00518 0.00014
##
## P value adjustment method: bonferroni1-2、1-3、2-4、3-4 間に有意差が認められた。有意確率は大きくなった。
ホルム (Holm) 補正
ボンフェローニ補正を改良したものである。すべての比較組み合わせの \(t\) 値を計算し、それを大きさの順に並べる。一番大きな \(t\) 値 \(t_{(1)}\) の有意確率を \(\alpha/r\)、次の大きさの \(t_{(2)}\) の有意確率を \(\alpha/(r - 1)\)、というように有意確率を調整する。R ではホルム補正がデフォルトで、\(p\) 値を大きさの順に並べ最も小さな \(p\) 値を \(r\) 倍し、次に大きな \(p\) 値を \(r - 1\) 倍して出力する。
pairwise.t.test(times, diets) # Rでのデフォルト##
## Pairwise comparisons using t tests with pooled SD
##
## data: times and diets
##
## 1 2 3
## 2 0.01141 - -
## 3 0.00090 0.31755 -
## 4 1.00000 0.00345 0.00014
##
## P value adjustment method: holm1-2、1-3、2-4、3-4 間に有意差が認められた。有意確率がボンフェローニ補正より小さくなった理由は不明。
チューキー（Tukey）の HSD（honestly significant difference）
Tukey は HSD 検定と呼ばれる対比較のために特別にデザインされた検定方法を提案した。HSD は、分散分析を行なった結果をもとに計算できる。
av <- aov(times ~ diets) # 一元配置分散分析
TukeyHSD(av)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = times ~ diets)
##
## $diets
## diff lwr upr p adj
## 2-1 5 0.7245544 9.275446 0.0183283
## 3-1 7 2.7245544 11.275446 0.0009577
## 4-1 0 -4.0560438 4.056044 1.0000000
## 3-2 2 -1.8240748 5.824075 0.4766005
## 4-2 -5 -8.5770944 -1.422906 0.0044114
## 4-3 -7 -10.5770944 -3.422906 0.00012682-1、3-1、4-2、4-3 間に有意差が認められた。
図示する
plot(TukeyHSD(av), las = 2)
あらためて、各群のとる値を確認してみよう。
boxplot(times ~ diets)
2 二元配置 (Two-way layout)
多くの実験は、2 つ以上の要因を含む。例えば、温度や濃度などである。要因の影響を詳しくみるため、要因ごとに複数の段階を設けることがある。これを水準という。ここでは2つの要因 A、B を用いた実験を考える。要因 A は \(a\) 水準、要因 B は \(b\) 水準からなっているとする。
2.1 繰り返し (repetition) の無い二元配置
構造モデル
まず最初に、要因 A、B の各水準ごとの組合せのデータが 1つしかない場合である。要因 A の第 \(i\) 水準、要因 B の第 \(j\) 水準に対する応答を \(y_{ij}\) とすると、 \[ y_{ij} = \mu + \alpha_i + \beta_j + e_{ij}\\ e_{ij} \overset{\textrm{iid}}{\sim} N(0, \sigma^2), i = 1,\ldots,a, \ j=1,\ldots,b \] とおける。ここに、\(\mu\) は総平均、\(\alpha_i\) は要因 A の第 \(i\) 水準の効果、\(\beta_j\) は要因 B の第 \(j\) 水準の効果で、\(\sum_i \alpha_i =0、\sum_j \beta_j = 0\) を満たし、\(\e_{ij}\) はすべての \(i、j\) で互いに独立な誤差項で、\(e_{ij} \sim N(0, \sigma^2)\) である。
このモデルにおける帰無仮説は、要因 A、B に対する \[ {\rm H}^\alpha_0 : \alpha_1 = \alpha_2 = \cdots \alpha_a = 0 \\ {\rm H}^\beta_0 : \beta_1 = \beta_2 = \cdots \beta_b = 0 \] の2つである。
平方和分解
\(y_{ij}\) の変動は以下のように分解される。 \[ S_T = \sum_{i=1}^a \sum_{j=1}^b (y_{ij} - \bar{y}_{..})^2 \\ = b \sum_{i=1}^a (\bar{y}_{i.} - \bar{y}_{..})^2 + a \sum_{j=1}^b (\bar{y}_{.j} - \bar{y}_{..})^2 + \sum_{i=1}^a \sum_{j=1}^b (y_{ij} - \bar{y}_{i.} - \bar{y}_{.j} + \bar{y}_{..})^2 \\ = S_A + S_B + S_E \] ここで、 \[ \bar{y}_{..} = \frac{1}{ab} \sum_{i=1}^a \sum_{j=1}^b y_{ij}\\ \bar{y}_{i.} = \frac{1}{b} \sum_{j=1}^b y_{ij}\\ \bar{y}_{.j} = \frac{1}{a} \sum_{i=1}^a y_{ij} \] である。また各平方和に対応する自由度は \[ ab-1 = (a-1) + (b-1) + (a-1)(b-1) \] になる。
分散分析表
分散分析表は、以下のようになり、この表から要因 A、B に対する帰無仮説の検定が行える。
| 要因 | 自由度 | 平方和 | 平均平方 | \(F\) 値 | \(p\) 値 |
|---|---|---|---|---|---|
| A | \(a-1\) | \(S_A\) | \(M_A=S_A/(a-1)\) | \(F_A = M_A/M_E\) | \({\bf P}[F_{a-1,\ (a-1)(b-1)} > F_A]\) |
| B | \(b-1\) | \(S_B\) | \(M_B=S_B/(b-1)\) | \(F_B = M_B/M_E\) | \({\bf P}[F_{b-1,\ (a-1)(b-1)} > F_B]\) |
| 誤差 | \((a-1)(b-1)\) | \(S_E\) | \(M_E=S_E/(a-1)(b-1)\) | ||
例2．2 品種の収量比較
イネの 2 つの品種 A、B の収量に対し、10 箇所の場所で栽培試験を行った結果が以下の表である。
| 試験地 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| 品種 A | 38.6 | 46.8 | 55.8 | 53.4 | 47.9 | 54.4 | 61.4 | 63.7 | 44.4 | 43.6 |
| 品種 B | 42.0 | 44.6 | 52.0 | 50.8 | 45.2 | 52.4 | 59.3 | 60.7 | 39.1 | 40.9 |
このデータから品種 A、B の収量に差があるかを検定してみよう。
まずはデータを準備する。
y <- c(38.6, 46.8, 55.8, 53.4, 47.9, 54.4, 61.4, 63.7, 44.4, 43.6, 42.0, 44.6,
52.0, 50.8, 45.2, 52.4, 59.3, 60.7, 39.1, 40.9) # 収量データ
variety <- factor(c(rep("A", 10), rep("B", 10))) # 品種ラベル 因子（factor）とする
place <- factor(rep(1:10, 2)) # 場所ラベル 因子（factor）とする
df <- data.frame(y, variety, place)
head(df)## y variety place
## 1 38.6 A 1
## 2 46.8 A 2
## 3 55.8 A 3
## 4 53.4 A 4
## 5 47.9 A 5
## 6 54.4 A 6aov関数を用いて、繰り返しのない2元配置分散分析を行う。
# 繰り返しのない2元配置
av <- aov(y ~ variety + place, data = df) # aov() : 分散分析 (Analysis Of Variance) を行う関数
summary(av)## Df Sum Sq Mean Sq F value Pr(>F)
## variety 1 26.4 26.45 10.63 0.00983 **
## place 9 1079.0 119.89 48.19 1.52e-06 ***
## Residuals 9 22.4 2.49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lm関数を用いて行うこともできる。
model <- lm(y ~ variety + place, data = df) # lm() : 線形モデル (Linear Model) を行う関数
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## variety 1 26.45 26.450 10.632 0.009828 **
## place 9 1079.05 119.894 48.193 1.52e-06 ***
## Residuals 9 22.39 2.488
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1結果は、品種、場所ともに有意である。
2.2 繰り返し (repetition) のある二元配置
要因 A の第 \(i\) 水準、要因 B の第 \(j\) 水準に対して \(n\) 回の繰り返し (repetition) がある場合である。繰り返しがあると、要因 A、B の主効果 (main effect) に加え、交互作用 (interaction) の検定が行える。
\(x\) 軸を施肥量、\(y\) 軸を収量とし、黒丸が品種Aの応答、赤丸が品種Bの応答であるとする。 左側の図は交互作用が無い場合で、品種 A、B とも施肥量が増加すると収量が増加し、 また、平均的にみて品種 A の収量の方が大きい。 交互作用が無い場合は主効果ごとの効果を単純に加え合わせれば良いので解釈が容易である。 一方、右図は交互作用が有る場合で、品種Aは施肥量が増えると収量も増加するが、 品種Bは施肥量が増えると収量は低下してしまった。 このように交互作用がある場合には、要因の水準ごとの組合せごとに考えなくてはならないので、 解釈が複雑になってしまう。
構造モデル
要因の水準の組合せごとに繰り返し数が異なると、平方和分解が行えず、計算や解釈が面倒になる。 ここでは、繰り返し数が \(n\) である釣り合い型の場合を考える。要因 A の第 \(i\) 水準、要因 B の第 \(j\) 水準の \(k\) 番目の繰り返しに対する応答を \(y_{ijk}\) とすると、 \[ y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + e_{ijk}\\ e_{ijk} \overset{\textrm{iid}}{\sim} N(0, \sigma^2),\ i = 1,\ldots,a, \ j=1,\ldots,b, \ k=1,\cdots,n \] とおける。ここに、\(\mu\) は総平均、\(\alpha_i\) は要因 A の第 \(i\) 水準の主効果、\(\beta_j\) は要因 B の第 \(j\) 水準の主効果で、\(\sum_i \alpha_i =0、\sum_j \beta_j = 0\) を満たし、\((\alpha\beta)_{ij}\) は交互作用で、\(\sum_i (\alpha_i\beta_j) = \sum_j (\alpha\beta)_{ij} = 0\) を満たし、\(e_{ijk}\) はすべての \(i、j、k\) で互いに独立な誤差項である。
このモデルにおける帰無仮説は、 \[ {\rm H}^\alpha_0 : \alpha_1 = \alpha_2 = \cdots \alpha_a = 0 \\ {\rm H}^\beta_0 : \beta_1 = \beta_2 = \cdots \beta_b = 0 \\ {\rm H}^{\alpha\beta}_0 : (\alpha\beta)_{11} = (\alpha\beta)_{12} = \cdots = (\alpha\beta)_{ab} = 0 \] の 3 つである。
平方和分解
\(y_{ijk}\) の変動は以下のように分解される。 \[ S_T = \sum_{i=1}^a \sum_{j=1}^b \sum_{k=1}^n (y_{ijk} - \bar{y}_{...})^2 \\ = bn \sum_i (\bar{y}_{i..} - \bar{y}_{...})^2 + an \sum_j (\bar{y}_{.j.} - \bar{y}_{...})^2 + n\sum_i \sum_j (\bar{y}_{ij.} - \bar{y}_{i..} - \bar{y}_{.j.} + \bar{y}_{...})^2 + \sum_i\sum_j\sum_k (\bar{y}_{ij.} - \bar{y}_{...})^2 \\ = S_A + S_B + S_{AB} + S_E \] ここで、 \[ \bar{y}_{...} = \frac{1}{abn} \sum_i \sum_j \sum_k y_{ijk}\\ \bar{y}_{i..} = \frac{1}{bn} \sum_j \sum_k y_{ijk}\\ \bar{y}_{.j.} = \frac{1}{an} \sum_i \sum_k y_{ijk}\\ \bar{y}_{ij.} = \frac{1}{n} \sum_k y_{ijk} \] である。また各平方和に対応する自由度は \[ abn-1 = (a-1) + (b-1) + (a-1)(b-1) + ab(n-1) \] である。
分散分析表
分散分析表は、以下のようになり、この表から要因 A、B の主効果や交互作用に対する帰無仮説の検定が行える。
| 要因 | 自由度 | 平方和 | 平均平方 | \(F\) 値 | \(p\) 値 |
|---|---|---|---|---|---|
| 主効果 A | \(a-1\) | \(S_A\) | \(M_A=S_A/(a-1)\) | \(F_A = M_A/M_E\) | \({\bf P}[F_{a-1,\ ab(n-1)} > F_A]\) |
| 主効果 B | \(b-1\) | \(S_B\) | \(M_B=S_B/(b-1)\) | \(F_B = M_B/M_E\) | \({\bf P}[F_{b-1,\ ab(n-1)} > F_B]\) |
| 交互作用 | \((a-1)(b-1)\) | \(S_{AB}\) | \(M_{AB}=S_{AB}/(a-1)(b-1)\) | \(F_{AB} = M_{AB}/M_E\) | \({\bf P}[F_{(a-1)(b-1),\ ab(n-1)} > F_{AB}]\) |
| 誤差 | \(ab(n-1)\) | \(S_E\) | \(M_E=S_E/ab(n-1)\) | ||
例3．インスリン治療
ブタによる動物実験で、インスリンの効果を検証した。供試された動物は一般的な麻酔が施された。ここで、
- LPS: 180 分間のリボ多糖注入
- HEC: 570 分間インスリン静脈内注入
の２つの処理を行った。LPS は全身性炎症を誘発し、HEC はその治療として作用する。以下の4つの実験を行った。 (1) 麻酔のみ（LPS 無し、HEC 無し）、(2) LPS のみ、(3) HEC のみ、(4) LPS および HEC。330 分経過後の腎臓中のインターロイキン -10 (IL-10) のレベルを計測した。各処理組合せにおいて \(n=8\) の釣り合い型のデータを得た。このデータを解析してみよう。
| LPS | HEC | ||||||||
|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 7.0607 | 2.6168 | 4.3489 | 3.6356 | 4.7510 | 2.9530 | 3.6137 | 5.6969 |
| 1 | 0 | 3.6911 | 4.3933 | 6.0513 | 4.2529 | 4.5554 | 3.8447 | 1.3590 | 2.1449 |
| 0 | 1 | 3.0693 | 1.6489 | 2.9160 | 2.9149 | 2.1102 | 3.1004 | 4.1170 | 3.0229 |
| 1 | 1 | 2.4159 | 3.1493 | 4.4462 | 2.8545 | 1.8944 | 3.5133 | 4.6254 | 3.8967 |
まずは、データを準備する。
# 2元配置
LPS <- as.factor(rep(c(0, 1, 0, 1), each = 8))
HEC <- as.factor(rep(c(0, 0, 1, 1), each = 8))
IL10 <- c(7.0607, 2.6168, 4.3489, 3.6356, 4.7510, 2.9530, 3.6137, 5.6969,
3.6911, 4.3933, 6.0513, 4.2559, 4.5554, 3.8447, 1.3590, 2.1449,
3.0693, 1.6489, 2.9160, 2.9149, 2.1102, 3.1004, 4.1170, 3.0229,
2.4159, 3.1493, 4.4462, 2.8545, 1.8944, 3.5133, 4.6254, 3.8967)
df <- data.frame(IL10, LPS, HEC)aov関数を用いて、交互作用を含む二元配置の分散分析を行う。
# 分散分析関数
av <- aov(IL10 ~ LPS * HEC, data = df) # 先ほどは + だったところを * にする
summary(av)## Df Sum Sq Mean Sq F value Pr(>F)
## LPS 1 0.01 0.007 0.005 0.9436
## HEC 1 7.29 7.293 5.053 0.0326 *
## LPS:HEC 1 2.14 2.141 1.483 0.2334
## Residuals 28 40.41 1.443
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lm関数を用いて、分散分析を行うこともできる。
# 線形モデル (linear model) 関数
model <- lm(IL10 ~ LPS * HEC, data = df)
anova(model)## Analysis of Variance Table
##
## Response: IL10
## Df Sum Sq Mean Sq F value Pr(>F)
## LPS 1 0.007 0.0073 0.0051 0.94363
## HEC 1 7.293 7.2932 5.0532 0.03264 *
## LPS:HEC 1 2.141 2.1409 1.4834 0.23341
## Residuals 28 40.412 1.4433
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1分散分析の結果、LPS の効果は有意でない。\(p 値= 0.9436\) である。交互作用も \(p 値 = 0.2334\) なので有意でない。HEC のみが5%有意 (\(p 値 = 0.0326\)) であった。
ここで、図を描いて効果を確認してみよう。
summary(model)##
## Call:
## lm(formula = IL10 ~ LPS * HEC, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.42795 -0.72872 0.05565 0.56202 2.72612
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3346 0.4247 10.205 6.14e-11 ***
## LPS1 -0.5476 0.6007 -0.912 0.3697
## HEC1 -1.4721 0.6007 -2.451 0.0208 *
## LPS1:HEC1 1.0346 0.8495 1.218 0.2334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.201 on 28 degrees of freedom
## Multiple R-squared: 0.1894, Adjusted R-squared: 0.1025
## F-statistic: 2.181 on 3 and 28 DF, p-value: 0.1126#
coef(model)## (Intercept) LPS1 HEC1 LPS1:HEC1
## 4.334575 -0.547625 -1.472125 1.034637#
cell00 <- coef(model)[1] #
cell01 <- coef(model)[1] + coef(model)[3] #
cell10 <- coef(model)[1] + coef(model)[2] #
cell11 <- sum(coef(model)) #
#
op <- par(mfrow = c(1, 2))
plot(c(0,1), c(cell00, cell01), ylim=c(0,5), xlab="HEC", ylab="IL10", type="b")
lines(c(0,1), c(cell10, cell11), type ="b", col="red")
legend("bottomright", c("LPS0","LPS1"), pch=1, col=c("black", "red"))
plot(c(0,1), c(cell00, cell10), ylim=c(0,5), xlab="LPS", ylab="IL10", type="b")
lines(c(0,1), c(cell01, cell11), type="b", col="red")
legend("bottomright", c("HEC0","HEC1"), pch=1, col=c("black", "red"))
par(op)上の図は、処理の平均をプロットしたものであり、交互作用プロットとよばれる。\(x\) 軸は要因の水準を表しており、この場合は無 (0) と有 (1) の 2 水準がある。\(x\) 軸が HEC のときが左図で、\(x\) 軸が LPS のときが右図である。どちらも同じ情報を異なるかたちで示している。
交互作用プロットが平行であるときは交互作用な無く、平行で無いときは交互作用がある。このグラフは平行でないので、交互作用があるように見えるが、交互作用の \(p 値 = 0.233\) なので、有意では無い。これは処理組合せごとの誤差が大きかったので、交互作用を検出するに至らなかった、と考えられる。
主効果は HEC が有意であった。これは、左図では黒 (LPS0) と赤 (LPS1) が交差しているので、LPS 間の主効果が無さそそうなことが分かる。一方、右図では黒 (HEC0) の方が赤 (HEC1) より上にあるので、有意差が有りそうなことが分かり、\(p 値 = 0.0326\) で有意差があることが確認できる。
参考：不釣り合い型分散分析
二元配置分散分析では、すべての処理組合せでの標本サイズが等しい釣り合い型であるのが普通であるが、場合によっては、実験の失敗などでデータの欠測があり、不釣り合い型になってしまう場合がある。不釣り合い型の分散分析では、平方和分解のやり方に複数の考え方があるので結果が多少異なってしまう。この講義では話が複雑になるので、詳しくは取り上げないが、どうなるかだけ、ここで示す。
# 不釣り合い型分散分析
# 5番目のデータを人為的に欠測にする。
dim(df)## [1] 32 3df2 <- df[-5,] # 5行目を除く
dim(df2)## [1] 31 3# 分散分析関数
av2 <- aov(IL10 ~ LPS * HEC, data = df2)
summary(av2)## Df Sum Sq Mean Sq F value Pr(>F)
## LPS 1 0.02 0.017 0.011 0.9163
## HEC 1 6.38 6.379 4.283 0.0482 *
## LPS:HEC 1 1.84 1.836 1.233 0.2766
## Residuals 27 40.21 1.489
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 線形モデル (linear model) 関数
model2 <- lm(IL10 ~ LPS * HEC, data = df2)
anova(model2)## Analysis of Variance Table
##
## Response: IL10
## Df Sum Sq Mean Sq F value Pr(>F)
## LPS 1 0.017 0.0168 0.0113 0.91631
## HEC 1 6.379 6.3793 4.2831 0.04819 *
## LPS:HEC 1 1.836 1.8362 1.2329 0.27664
## Residuals 27 40.214 1.4894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1参考：対比
処理平均 \(\mu_i\) のある群と他の群との違いに特に興味がある場合がある。例えば、処理数が 5 であった時、処理 1、2、3 の平均と処理 4、5 の平均との間 \[ \frac{\mu_1+\mu_2+\mu_3}{3} - \frac{\mu_4+\mu_5}{2} \] に差があるかをみたいような場合である。一般に、 \[ \sum_i c_i \mu_i, \ \sum_i c_i = 0 \] である比較を対比 (contrast) という。また線形結合、\(2\mu_1 - \mu_2 - \mu_4\)、\(\mu_3 - \mu_2\)、\(\mu_1+\mu_2+\mu_3-\mu_4-2\mu_5\) も対比である。
対比 \(C = \sum_i c_i \mu_i\) は、\(\hat{C} = \sum_i c_i \bar{X}_{i.}\) で推定されるが、帰無仮説 (\(\alpha_i =0\)) のもとでは、対比推定量の平均と分散は、 \[ E[\hat{C}] = E[\sum_i c_i \hat{X}_{i.}] = \sum_i c_i E[\bar{X}_{i.}] = \sum_i c_i \mu = \mu \sum_ic_i = 0 \\ Var[\hat{C}] = Var[\sum_i c_i \hat{X}_{i.}] = \sum_i c_i^2 Var[\hat{X}_{i.}] = \sum_i c_i^2 \frac{\sigma^2}{n_i} = \sigma^2 \sum_i \frac{c_i^2}{n_i} \] となるので、分散 \(\sigma^2\) をその推定量 \(s^2\) で置き換えた検定統計量 \(t\) は、 \[ t = \frac{\hat{C}-E[\hat{C}]}{\sqrt{\hat{Var}[\hat{C}]}} = \frac{\sum_i c_i \bar{X}_{i.}}{s \sqrt{\sum_i \frac{c_i^2}{n_i}}} \sim t_{n-k} \] のように自由度 \(n-k\) の \(t\) 分布に従うので、対比 \(C= \sum_ic_i\mu_i=0\) の検定を行うことができる。特に、各処理水準の標本サイズが \(n_i=n\) と一定で釣り合っていて、対比の係数を \(\sum c_i^2 = 1\) と標準化すると、検定統計量は、 \[ t = \frac{\sum_ic_i \bar{X}_{i.}}{s/\sqrt{n}} \sim t_{n-k} \] と簡略化される。
ところで、別の対比 \(\sum_i d_i \mu_i\) があったときに、\(\sum_i c_i d_i =0\)、となるものを直交対比という。直交対比の組は同時に検定が行える。たとえば、5 つの処理があり、処理 1、2、3 と処理 4、5 の平均間の差を考えたとき、対比係数ベクトルは、 \[ c_1 = (1/3, 1/3, 1/3, -1/2, -1/2) \] となるが、これに直交する対比としては、たとえば、 \[ c_2 = (0,0,0,1/2,-1/2), \ c_3 = (1, -1/2, -1/2,0,0) \] などが考えられる。これらの対比は互いに直交しているので、同時に検定ができる。事前に検定したい対比が決められているときは、後述する多重比較による有意性の補正を行わない。
例1の凝固時間データにおいて、\(\mu_1+\mu_2 - (\mu_3 + \mu_4) = 0\) を検定してみよう。対比係数ベクトルは、\(c = (1,1,-1,-1)\) である。
c1 <- c(1,1,-1,-1) # 対比係数ベクトル
# 対比を定義する。
contrasts(diets) <- c1
# 作成された対比の表示 (他の直交対比は勝手につくられる)
contrasts(diets)## [,1] [,2] [,3]
## 1 1 -0.5 -0.5
## 2 1 0.5 0.5
## 3 -1 0.5 -0.5
## 4 -1 -0.5 0.5fc <- lm(times ~ diets)
# 結果表示。diets1 のみに着目
summary(fc)##
## Call:
## lm(formula = times ~ diets)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.00 -1.25 0.00 1.25 5.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.0000 0.4979 128.537 < 2e-16 ***
## diets1 -0.5000 0.4979 -1.004 0.327
## diets2 6.0000 0.9958 6.025 6.85e-06 ***
## diets3 -1.0000 0.9958 -1.004 0.327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.366 on 20 degrees of freedom
## Multiple R-squared: 0.6706, Adjusted R-squared: 0.6212
## F-statistic: 13.57 on 3 and 20 DF, p-value: 4.658e-05\(\mu_1+\mu_2 - (\mu_3 + \mu_4) = 0\) は \(p 値= 0.327\) で有意では無かった。しかし、これの直交対比である \((\mu_2+\mu_3)/2 - (\mu_1 + \mu_4)/2 = 0\) は \(p \approx 0\) で強度に有意となった。しかしながら、この対比が生物学的に意味が無い場合は、「たまたま」と考えた方が良いと思われる。