2群間の平均と分散の比較
Hiroyoshi Iwata
2024/11/01
1. はじめに
研究を行っていると、2つのグループ（以下、「群」とよぶ）間で観察値に違いが無いか、明らかにしたい場合が多々ある。例えば、ある処理を行った群と、行わなかった群で、何らかの特徴を計測して、群間で計測値の平均や分散に違いがあるかどうか、が知りたい場合などである。本日の講義では、2群間で平均や分散に違いがあるかどうかを統計的に検定する方法を説明する。
2. 2つの独立した正規母集団の平均と分散
最初に、2つの母集団の平均の検定を行う例を示す。
例1 エンドファイトの効果
植物の内生菌（エンドファイト）の中には、植物の病害や環境ストレスに対する抵抗性を高め、成長を旺盛にするものがある。あるエンドファイトをダイズの幼苗に接種して、その後乾燥条件下で栽培試験を行った。栽培試験後に植物体の乾物重を計測した結果、以下のような結果が得られた。エンドファイトを接種することによって、乾燥耐性は向上したのだろうか？
| 処理 | 個体数 | 乾物重の平均(g) | 乾物重の標準偏差(g) |
|---|---|---|---|
| エンドファイト接種 | 18 | 19.5 | 5.3 |
| エンドファイト非接種 | 16 | 14.8 | 4.4 |
この問題を一般化するために、2群の標本を、それぞれ、\(X_{11}, X_{12}, ..., X_{1,n_1}\)と、\(X_{21}, X_{22}, ..., X_{2,n_2}\)と表し、これらが正規分布に従う母集団\(N(\mu_1, \sigma_1^2)\)と\(N(\mu_2, \sigma_2^2)\)から得られたものとする。
上の例では、帰無仮説\(H_0: \mu_1 = \mu_2\)（2群の標本の母平均は等しい）を、対立仮説\(H_1: \mu_1 >,\ne,< \mu_2\)（対立仮説は3通り設定できる。1つは、2群の母平均が等しくない、というものであり、残りの2つは、どちらかの群の母平均が他方より大きいというものである）に対して、有意水準\(\alpha\)で検定することになる。
上の例の場合、検定される帰無仮説は、エンドファイトが乾燥耐性に影響を及ぼしているかどうか、すなわち、2つの母集団の平均が同じという仮説 \[ H_0: \mu_1 = \mu_2 \] が成り立っているか、である。
ここで、母集団は、エンドファイト処理をされた全てのダイズ個体、または、エンドファイト処理をされていない全てのダイズ個体である。なお、対立仮説\(H_1\)は、場合によって、片側、または、両側になる。両側の場合は、単純に\(H_1: \mu_1 \ne \mu_2\)、すなわち、2つの母集団平均は等しくないが、差はどちら側にも動き得ると考える。片側検定の選択は、問題設定に従って行う。この例では、エンドファイトの接種が、乾燥耐性を低下させるという片側検定の対立仮説\(H_1: \mu_1 < \mu_2\)を設定することは、エンドファイトによる乾燥耐性の向上にのみ興味がある場合は意味をなさない。その場合、目的にかなった片側仮説は、\(H_1: \mu_1 > \mu_2\)となる。
なお、母集団の平均の検定は、以下の2つのシナリオのもとで行わる。
シナリオ1：分散が未知だが等しい場合
このシナリオでは、2群の母集団の分散\(\sigma^2\)（分散は2群で等しいので1つの記号で表すことができる）を、2群の標本分散\(s_1^2\)と\(s_2^2\)から推定する。各グループの標本の大きさを\(n_1\)と\(n_2\)とすると、\(s_1^2\)と\(s_2^2\)を以下のように重み付けして平均をとる。
\[ s_p^2 = \frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}=\frac{n_1-1}{n_1+n_2-2}s_1^2+\frac{n_2-1}{n_1+n_2-2}s_2^2=ws_1^2+(1-w)s_2^2 \] これは、プールされた標本分散とよばれ、各群で推定される標本分散（\(s_1^2\)または\(s_2^2\)）に比べて、母集団の分散のより正確な推定値となる。なお、\(s_p^2\)の平方根をプールされた標本標準偏差とよび、\(s_p\)と表す。帰無仮説\(H_0\)が真のとき、すなわち、\(\mu_1 = \mu_2\)のとき、統計量 \[ t = \frac{\bar{X_1}-\bar{X_2}}{s_p\sqrt{1/n_1 + 1/n_2}} \] は、自由度\(df = n_1 + n_2 - 2\)の\(t\)分布に従う。
シナリオ2: 分散に関する仮定が無い（等しいと仮定できない）場合
2群の分散が等しいと仮定できない場合は少なくない。そのような場合には、以下の統計量に基づいて検定を行う。帰無仮説\(H_0\)が真のとき、すなわち、\(\mu_1 = \mu_2\)のとき、統計量 \[ t = \frac{\bar{X_1}-\bar{X_2}}{\sqrt{s_1^2/n_1 + s_2^2/n_2}} \] が、自由度 \[ df = \frac{(s_1^2/n_1 + s_2^2/n_2)^2}{(s_1^2/n_1)^2/(n_1-1)+(s_2^2/n_2)^2/(n_2-1)} \] をもつ\(t\)分布に従う。
両シナリオにおいて、
| 対立仮説 | \(\alpha\)水準の棄却領域 | \(p\)値 |
|---|---|---|
| \(H_1: u_1 > u_2\) | \([t_{df,1-\alpha},\infty)\) | 1 - pt(t, df) |
| \(H_1: u_1 \ne u_2\) | \((-\infty, t_{df,\alpha/2}]\cup[t_{df,1-\alpha/2},\infty)\) | 2 * pt(-abs(t), df) |
| \(H_1: u_1 < u_2\) | \((-\infty, t_{df,\alpha}]\) | pt(t, df) |
として検定する。
自由度3のときの棄却限界値を視覚化して確認してみる。
df <- 3
x <- seq(-5, 5, 0.1)
y <- dt(x, df)
plot(x, y, type = "l")
abline(h = 0)
p <- 0.05
vl <- qt(p, df)
abline(v = vl, col = "blue") # 下側5%
vu <- qt(1 - p, df)
abline(v = vu, col = "red") # 上側5%
vl <- qt(p/2, df)
abline(v = vl, col = "orange") # 両側5%
vu <- qt(1 - p/2, df)
abline(v = vu, col = "orange") # 両側5%
オレンジ色が、対立仮説\(H_1: u_1 \ne u_2\)の棄却限界値、赤色が、対立仮説\(H_1: u_1 > u_2\)の棄却限界値、青色が、対立仮説\(H_1: u_1 < u_2\)の棄却限界値である。
実験方法、標本の収集方法、その他の因子の性質に基づいて、分散が等しいとも、あるいは等しくないとも仮定できる。例えば、例1では、母分散は未知であり、等しいかどうかは分からない。検定のシナリオを決める前に、母分散が等しいかどうかを検定こともできる。
ここでは、平均の検定の説明を一旦中断して、分散の同一性の検定について説明する。
2つの正規分布の分散の同一性の検定
平均を検定する前に分散を検定することについては賛否両論がある。それは、分散の検定が平均の検定に比べて、正規分布からの逸脱に対して頑健でないためである。ただし、検定結果がボーダーラインにあるような場合を除いて、選ばれたシナリオがどちらであっても、ほとんど検定結果に影響しない。したがって、分散の同一性の仮定が疑わしい場合は、シナリオ1よりも保守的なシナリオ2を選べばよい（シナリオ2では、等分散性を仮定しなくてもよい）。
標本\(X_{11}, X_{12}, ..., X_{1,n_1}\)と、\(X_{21}, X_{22}, ..., X_{2,n_2}\)は、それぞれ、正規分布する母集団\(N(\mu_1, \sigma_1^2)\)と\(N(\mu_2, \sigma_2^2)\)から得られたとする。このとき、分散の同一性検定検定は、帰無仮説 \[ H_0: \sigma_1^2 = \sigma_2^2 \] に対して、対立仮説 \[ H_1: \sigma_1^2 \ne \sigma_2^2 \] を考える。
なお、検定のための統計量には、標本分散の比\(F=s_1^2/s_2^2\)を用いる。\(F\)は、帰無仮説\(H_0: \sigma_1^2 = \sigma_2^2\)が真のとき、自由度\(n_1-1, n_2-1\)の\(F\)分布に従う。そこで、次コマンドで計算される\(p\)値に従って、どちらのシナリオにするか決定する。
p = 2 * min(pf(F, n1-1, n2-1), 1-pf(F, n1-1, n2-1))以下は、等分散性の検定について、片側、および、両側検定についてまとめると、以下の通りとなる。
| 対立仮説 | \(\alpha\)水準の棄却領域 | \(p\)値 |
|---|---|---|
| \(H_1: \sigma_1^2 > \sigma_2^2\) | \([F_{df_1, df_2,1-\alpha},\infty)\) | 1 - pf(F, df1, df2) |
| \(H_1: \sigma_1^2 \ne \sigma_2^2\) | \([0, F_{df_1, df_2,\alpha/2}]\cup[F_{df_1,df_2,1-\alpha/2},\infty)\) | 2 * min(pf(F,df1,df2), 1-pf(F,df1,df2)) |
| \(H_1: \sigma_1^2 < \sigma_2^2\) | \([0, F_{df_1,df_2,\alpha}]\) | pf(F, df1, df2) |
例2. 例1では、等分散性を検定するための\(F\)の値は1.45であり、等分散性の仮説は有意度水準\(\alpha = 0.05\)で棄却されない（\(p > 0.473\)）。
n1 <- 18
X1bar <- 19.5
s1 <- 5.3
n2 <- 16
X2bar <- 14.8
s2 <- 4.4
Fstat <- s1^2 / s2^2
Fstat## [1] 1.45093pval <- 2 * min(pf(Fstat, n1-1, n2-1), 1-pf(Fstat, n1-1, n2-1))
pval## [1] 0.4739491ここで、なぜ、2つの確率を計算して、小さい方を選ぶのかを説明する。
pf(Fstat, n1-1, n2-1)## [1] 0.76302541-pf(Fstat, n1-1, n2-1)## [1] 0.2369746小さい方は、後者である。
次に、群1と群2を逆にして計算してみる。
# F値の計算は逆（逆数）になる。
Fstat <- s2^2 / s1^2
Fstat## [1] 0.6892132# 自由度も逆にする必要がある
n1 <- 16; n2 <- 18
# p値を確認する
pf(Fstat, n1-1, n2-1)## [1] 0.23697461-pf(Fstat, n1-1, n2-1)## [1] 0.7630254\(F\)値は1以下の値となるが、大きい方の確率と、小さい方の確率が逆になるだけである。群1と群2のどちらが分散が大きくなるかが、場合によって異なるので、2通りの確率を計算して、小さい方を選んでいる。
なお、R には、2 つの正規分布する等分散性を検定するための var.test という専用の関数もある。var.testを用いるには、2群の標本データを直接入力する。
エンドファイトデータの元データを用いて、var.testを用いた検定を行ってみる。
X1 <- c(21.4, 21.2, 21.5, 20.8, 21.0, 28.8, 15.9, 14.4, 18.2, 16.4, 15.0, 14.2, 24.7, 16.9, 11.8, 17.0, 18.7, 32.6)
X2 <- c(16.9, 5.9, 11.7, 13.6, 14.2, 15.7, 15.1, 18.0, 19.6, 16.5, 14.9, 24.2, 17.4, 8.2, 9.9, 14.3)
var.test(X1, X2)##
## F test to compare two variances
##
## data: X1 and X2
## F = 1.409, num df = 17, denom df = 15, p-value = 0.5097
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5009352 3.8368293
## sample estimates:
## ratio of variances
## 1.409028標準偏差を計算する際の丸め誤差により、若干異なる結果が得られているが、2群の分散の違いは有意ではない。
なお、var.testでは、分散の比についての信頼区間も計算される。この95%信頼区間に1が含まれていれば、両側検定の場合において等分散（分散の比が1である）という帰無仮説は5%水準で棄却されない。
2つの母集団の平均の検定
ここでは、シナリオ1、すなわち、母分散は等しいと仮定し、プールされた標準偏差に基づく\(t\)検定を行ってみる。検定統計量は、 \[ t = \frac{\bar{X_1}-\bar{X_2}}{s_p\sqrt{1/n_1 + 1/n_2}} \] で、 \[ s_p = \sqrt{\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}} \] であり、\(t\)は、先述したように自由度\(n_1+n_2 - 2\)の\(t\)分布に従う。
sp <- sqrt(((n1 - 1)*s1^2 + (n2 - 1)*s2^2) / (n1 + n2 - 2))
sp## [1] 4.842746df <- n1 + n2 - 2
tstat <- (X1bar - X2bar) / (sp * sqrt(1/n1 + 1/n2))
tstat## [1] 2.824639pval <- 1 - pt(tstat, df)
pval## [1] 0.004041259帰無仮説は、1%水準で棄却される。
\(H_0\)を棄却域を用いて検定したいとする。対立仮説は片側で、右側の裾野に対して行うため、\(\mu_1 - \mu_2 > 0\)であり、棄却域は、\([t_{n_1+n_2 - 2, 1-\alpha}, \infty)\)である。
qt(1 - 0.01, df)## [1] 2.448678棄却域を用いた場合も、観察された\(t = 2.792343\)は棄却限界値である2.448678を超えているために、帰無仮説\(H_0\)は棄却される
3. 平均の差の信頼区間
母集団平均間の差の\((1-\alpha)100\)%信頼区間は簡単に導くことができる。検定のときと同様に、母集団の分散に関する仮定（シナリオ）に依存して方法が変わる。一般にその区間は、 \[ \biggl[ \bar{X_1} - \bar{X_2} - t_{df,1-\alpha/2} s^*,\quad \bar{X_1} - \bar{X_2} + t_{df,1-\alpha/2} s^* \biggr] \] と表される。
母集団について等分散の仮定が成り立つ場合は、自由度\(df = n_1 + n_2 - 2\)で、標本標準偏差は\(s^* = s_p \sqrt{1/n_1 + 1/n_2}\)となる。また、母集団に関する仮定が無い場合は、自由度\(df\)は\(\frac{(s_1^2/n_1 + s_2^2/n_2)^2}{(s_1^2/n_1)^2/(n_1-1)+(s_2^2/n_2)^2/(n_2-1)}\)となり、標本標準偏差は\(s^* = \sqrt{s_1^2/n_1 + s_2^2/n_2}\)となる。エンドファイトの例では、\(\mu_1 - \mu_2\)の95%信頼区間は、
sp <- sqrt(((n1 - 1)*s1^2 + (n2 - 1)*s2^2)/(n1 + n2 - 2))
df <- n1 + n2 - 2
LB <- X1bar - X2bar - qt(0.975, df) * sp * sqrt(1/n1 + 1/n2)
UB <- X1bar - X2bar + qt(0.975, df) * sp * sqrt(1/n1 + 1/n2)
c(LB, UB)## [1] 1.310687 8.089313となる。
両側検定に対する2群の平均の同一性の検定は、平均の差の信頼区間をもとに行うこともできる。具体的には、両側検定の有意水準\(\alpha\)のもとで\((1-\alpha)100\)%信頼区間を計算し、この区間が0を含んでいれば、帰無仮説は棄却されない。ここでは、母平均の差の95%信頼区間が\([1.271487, 8.128513]\)である。片側検定を行わず、両側検定を行った場合も、帰無仮説\(H_0\)は棄却される。
Rでは、2標本の\(t\)検定を行うために、t.testという専用の関数が準備されている。
ここでは、エンドファイトの元データをもとに、t.testを使って検定を行ってみる。
t.test(X1, X2)##
## Welch Two Sample t-test
##
## data: X1 and X2
## t = 2.8341, df = 31.921, p-value = 0.007903
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.326220 8.105725
## sample estimates:
## mean of x mean of y
## 19.47222 14.75625平均の同一性の検定が行われるだけでなく、差の信頼区間も計算される（ここでもまた、丸め誤差によって若干異なる値が計算されている）。
4. 標本に対応がある場合の正規分布の平均の同一性の検定
2群の平均を比較する場合に、観察値が様々な因子による影響を受け、比較を難しくする場合がある。例えば、人のある疾患の治療法の効果を調べたい場合に、治療法だけでなく、年齢、性別、BMI、リスクファクターの存在などの因子が被験者の間で制御されないかたちで異なっている場合がある。こうした場合、治療効果は、治療法の違いだけでなく、上述したような他の因子の影響も強く受けている可能性がある。
このような解析の邪魔になる因子による影響を排除して、検定の検出力をあげる方法うとして、ペアリング、マッチングという方法がある。標本がペアで構成されており、群1からの観察値が、群2の観察値に一対一でペアリングされている場合を考える。すなわち、群1からの標本\(X_{11},X_{12},...,X_{1,n}\)と群2からの標本\(X_{21},X_{22},...,X_{2,n}\)について、\(i\)番目のペア\((X_{1i},X_{2i})\)が、標本として得られるような状況である。通常、ペアになった観察値が、同じ被験者（試験前-試験後の観察値、偽薬（プラセボ）-新薬投与の結果、など）あるいは、関連する被験者（きょうだい）から得られる。あるいは、同じ植物個体を2つの異なる計測法で計測した場合などもこれに当てはまる。このように、例は様々ある。
通常は、標本が異なる可能性がある平均\(\mu_1, \mu_2\)、と未知の分散\(\sigma_1^2, \sigma_2^2\)をもつ正規分布から抽出されたと仮定する。なお、正規分布の線形結合（重み付けした和）はまた正規分布であり、したがって、差\(d_i = X_{1i} - X_{2i}\)も（正規分布の線形結合なので）正規分布にしたがう。 \[ d_i \sim N(\mu_1 - \mu_2, \sigma_1^2+\sigma_2^2-2\sigma_{12}) \] \(\sigma_{12}\)は、\(X_1\)と\(X_2\)の共分散である。ただし、この正規分布の分散\(\sigma_1^2+\sigma_2^2-2\sigma_{12}\)は未知であるため、正規分布に基づく検定はできない。
そこで、次の統計量を定義する。 \[ t = \frac{\bar{d}}{s_d / \sqrt{n}} \] \(\bar{d}\)、\(s_d\)は、差\(d_i\)の平均と標本標準偏差である。ここで、標本サイズ\(n\)は、ペアの数であり、観察値の数\(2n\)ではないことに注意する。なお、ペア間の差について考える場合には、母集団の等分散性を仮定する必要はない。
検定統計量\(t\)は、帰無仮説\(H_0: \mu_1 = \mu_2\)のもとで、自由度\(n-1\)の\(t\)分布に従う。なお、この検定は、1集団のときの\(t\)検定と一致する。その場合は、標本は「全てのペア間の差」であり、帰無仮説\(H_0\)は「ペア間の差の平均」が0に等しいという仮説になる。この手法は、対応のある\(t\)検定（paired t-test）と呼ばれる。対応のある\(t\)検定は、以下のようにまとめられる。
| 対立仮説 | \(\alpha\)水準の棄却領域 | \(p\)値 |
|---|---|---|
| \(H_1: u_1 > u_2\) | \([t_{n-1,1-\alpha},\infty)\) | 1 - pt(t, n-1) |
| \(H_1: u_1 \ne u_2\) | \((-\infty, t_{n-1,\alpha/2}]\cup[t_{n-1,1-\alpha/2},\infty)\) | 2 * pt(-abs(t), n-1) |
| \(H_1: u_1 < u_2\) | \((-\infty, t_{n-1,\alpha}]\) | pt(t, n-1) |
\(\delta = \mu_1 - \mu_2\)が、母平均間の違いを表すとき、\(\delta\)の\((1-\alpha)100\)%信頼区間は、 \[ \biggl[ \bar{d} - t_{n-1, 1-\alpha/2}\frac{s_d}{\sqrt{n}}, \bar{d} + t_{n-1, 1-\alpha/2}\frac{s_d}{\sqrt{n}} \biggr] \] として求められる。
観察値をマッチングすることで、被験者間の制御不能な因子の影響を緩和できるため、対応のある\(t\)検定は、そのような実験デザインが可能な場合には、2標本の\(t\)検定よりも正確な検定ができると期待される。
例3
以下のデータは、20人の被験者について、ある治療を行った結果である。治療開始前（before）と治療後（after）に疾患スコア（低いほうが疾患の程度が低い）を得た。得られた疾患スコアから、この治療法がスコアを下げることを示す十分な証拠となるか。有意水準5%で検定せよ。
before <- c(5.9, 7.6, 12.8, 16.5, 6.1, 14.4, 6.6, 5.4, 9.6, 11.6, 11.1, 15.6, 9.6, 15.2, 21, 5.9, 10, 12.2, 20.2, 6.2)
after <- c(5.2, 12.2, 4.6, 4, 0.4, 3.8, 1.2, 3.1, 3.5, 4.9, 11.1, 8.4, 5.8, 5, 6.4, 0, 2.7, 5.1, 4.8, 4.2)
d <- before - after
n <- length(d)
dbar <- mean(d)
sdd <- sd(d)
tstat <- dbar/(sdd/sqrt(n))
tstat## [1] 5.763691critpt <- qt(1 - 0.05, n - 1)
critpt## [1] 1.729133pval <- 1 - pt(tstat, n - 1)
pval## [1] 7.439804e-06\(t\)統計量は5.764であり、棄却域は、\([1.729, \infty)\)であるため、帰無仮説は有意水準5%で棄却される。\(p\)値は、\(7.440 \times 10^{-6}\)で高度に有意である。
同じ検定をt.testを使って行ってみる。paired = Tというオプションを指定する。また、療法の効果を知りたい（beforeのほうがafterよりスコアが大きい）のでalternative = “greater”を指定する。
t.test(before, after, paired = T, alternative = "greater")##
## Paired t-test
##
## data: before and after
## t = 5.7637, df = 19, p-value = 7.44e-06
## alternative hypothesis: true mean difference is greater than 0
## 95 percent confidence interval:
## 4.448472 Inf
## sample estimates:
## mean difference
## 6.355母平均間の差について、95%信頼区間を求める。
alpha <- 0.05
LB <- dbar - qt(1 - alpha/2, n - 1)*(sdd / sqrt(n))
UB <- dbar + qt(1 - alpha/2, n - 1)*(sdd / sqrt(n))
c(LB, UB)## [1] 4.047248 8.662752信頼区間をt.test関数を用いて計算するには、両側検定のオプションを指定する（デフォルトのままでよい）。
t.test(before, after, paired = T)##
## Paired t-test
##
## data: before and after
## t = 5.7637, df = 19, p-value = 1.488e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 4.047248 8.662752
## sample estimates:
## mean difference
## 6.355# 以下も結果は同じ
t.test(before, after, paired = T, alternative = "two.sided")##
## Paired t-test
##
## data: before and after
## t = 5.7637, df = 19, p-value = 1.488e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 4.047248 8.662752
## sample estimates:
## mean difference
## 6.3555. 検出力と必要な標本数の計算
最後に、二標本の平均の差の検定を行う場合を例として、検出力と、ある検出力を得るための標本数の計算について説明する。
まず、帰無仮説\(H_0: \mu_1 = \mu_2\)が正しい場合、\(t\)値は\(t\)分布に従うが、帰無仮説が正しく無い場合は\(t\)値は非心\(t\)分布に従う。対立仮説\(|\mu_1 - \mu_2| = \Delta\)が正しいとき、\(t\)値は、自由度\(n-1\)、非心パラメーター（non-central parameter) \[ \delta = \frac{\Delta}{\sqrt{\sigma_1^2/n_1 + \sigma_2^2/n_2}} \] の非心\(t\)分布に従う。
例1のもとで、対立仮説\(H_1: \mu_1 - \mu_2 = 4.0\)に対して片側検定の場合の検出力（\(1-\beta\)）を計算してみよう。
まず、例1の標本数と標本標準偏差を代入する。
n1 <- 18
s1 <- 6.3
n2 <- 16
s2 <- 4.4次に、帰無仮説\(H_0: \mu_1 = \mu_2\)が正しい場合と、対立仮説\(H_1: \mu_1 - \mu_2 = 4.0\)が正しい場合に従う分布を図示する。
Delta <- 4
delta <- Delta / sqrt(s1^2 / n1 + s2^2 / n2)
x <- seq(-5, 7, 0.01)
y <- dt(x, n1 + n2 - 2) # 通常のt分布（帰無仮説が正しい場合）
plot(x, y, type = "l", col = "red")
y2 <- dt(x, n1 + n2 - 2, ncp = delta) # 非心t分布
lines(x, y2, col = "blue")
abline(v = qt(1- 0.025, n1 + n2 - 2), col = "orange") # 両側5%の上側の棄却限界値
abline(h = 0) # y = 0に直線を引く
青色の分布が、対立仮説\(H_1: \mu_1 - \mu_2 = 4.0\)が正しい場合に\(t\)値が従う分布（非心\(t\)分布）、赤色の分布が帰無仮説\(H_0: \mu_1 = \mu_2\)が正しい場合に\(t\)値が従う分布（通常の\(t\)分布）である。オレンジ色の垂線は、有意水準5%で両側検定を行う場合の上側の棄却限界値である。
第2種の過誤を犯す確率\(\beta\)は、非心t分布のオレンジ色の線の左側の領域の面積となる。すなわち、
beta <- pt(qt(1 - 0.025, n1+n2-2), n1+n2-2, ncp = delta)
beta## [1] 0.4445835検出率\(1 - \beta\)は、
pow <- 1 - pt(qt(1 - 0.025, n1 + n2 - 2), n1 + n2 - 2, ncp = delta)
pow## [1] 0.5554165検出力は55%であることが分かる。
ここで、同じ現象について、将来行う実験の計画を立てるとする。データが収集され解析されたときに、有意水準5%の検定により、帰無仮説\(H_1: \mu_1 - \mu_2 = 4.0\)に対して、検出力\(1-\beta = 90\)%を達成したい。どのくらいの標本サイズが必要になるだろうか。
ここでは、現在得られているデータがパイロット研究であると仮定し、さらに、母分散\(\sigma^2\)は既知でプールされた標本分散に一致すると仮定する。このとき必要な標本数は、power.t.testという関数を使って計算することができる。5%水準で両側検定を行う場合、
sp <- sqrt(((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2))
power.t.test(power = 0.9, delta = Delta, sd = sp, sig.level = 0.05, alternative = "two.sided") ##
## Two-sample t test power calculation
##
## n = 40.59824
## delta = 4
## sd = 5.491841
## sig.level = 0.05
## power = 0.9
## alternative = two.sided
##
## NOTE: n is number in *each* group40個以上の標本が必要である。
1%水準で片側検定を行う場合、
sp <- sqrt(((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2))
power.t.test(power = 0.9, delta = Delta, sd = sp, sig.level = 0.01, alternative = "one.sided") ##
## Two-sample t test power calculation
##
## n = 50.45233
## delta = 4
## sd = 5.491841
## sig.level = 0.01
## power = 0.9
## alternative = one.sided
##
## NOTE: n is number in *each* group50個以上の標本が必要である。