僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。「母比率の検定」と書くべきところですが、理解しやすいようにここでは「比率の検定」というタイトルにしました。検定の前に二項分布について復習しておきましょう。
例題 1 「85歳以上」の転倒割合は、25.3％（約4人に1人）と報告されています（内閣府、平成17年）。この確率25.3％を転倒が起こる確率、つまり研究目的としている現象（転倒）が一般的に起こる確率 pとして進めます。S町に在住する85歳以上の高齢者10名のアンケート結果から、転倒経験ありと回答したのは6名であった（無作為抽出）。S町の85歳以上の高齢者の転倒割合は全国25.3％より高いと言えるか、有意水準5％で検証せよ！ ・帰無仮説：S町の転倒割合は、全国割合（25.3％）に等しい ・対立仮説：S町の転倒割合は、全国割合（25.3％）より高い
上側検定
ある比率より大きいことを実証した場合には、上側検定を実施します（ここ重要）
この検定をRで実行する場合は、以下のようになります
binom.test(
x=6, n=10, p=0.253,
alternative="greater"
)
上側検定なので alternative = “greater”
alternative hypothesis：対立仮説（Rは対立仮説を出力してくれるので分かりやすい）
p値は0.02となり（<0.05）、S町は全国に比べて転倒割合が有意に高いと言えます
検定の結果だけでよければ、ここまでで大丈夫です
ここからは少し詳しく勉強してみましょう
P値
ここから二項分布とP値について考えていきます
・帰無仮説：S町の転倒割合は、全国割合（25.3％）と等しい
・対立仮説：S町の転倒割合は、全国割合（25.3％）より高い
転倒回数を確率変数 $X$ と考えた場合、$X$ は二項分布に従う
\(X～B (10, 0.253)\)
この二項分布がどのようになっているのか・・・確率関数をグラフ化しましょう
面倒ですが、このような作業こそが”データサイエンスのお勉強”なのです。イメージを掴むことが大切ですね。
Rで確率分布のグラフを描いてみましょう
t <- 0:10
#確率関数の値
p <- dbinom(
x=t, size=10, prob=0.253
)
#棒グラフ（ｙ軸：確率関数の値）
barplot(
p ~ t, pch=16, xlab="", ylab=""
)
ｘ軸は転倒回数、ｙ軸は転倒回数が $x$ となる確率\(P(X=x)\)です
Rのdbinom関数は二項分布の確率関数から転倒回数ごとの確率を出力してくれます
基準となる転倒確率が25.3%の場合、10回中6回転倒する確率\(P(X=6)\)は？
dbinom(
x=6, size=10, prob=0.253
)
二項分布の確率関数に代入して確認してみましょう
サンプルサイズ $n=10$
転倒回数 $k=6$
基準となる転倒確立 $\pi=0.0253$
確率関数 \(P(X=k)={}_n C_kπ^k(1-π)^{n-k}\)
choose(10, 6)*(0.253^6)*(1-0.253)^4
当然ですが同じ値になります（chooseは組合せの関数）
0.01714845は、以下のグラフの赤色の部分になります
つまり基準となる転倒確率が25.3%の場合、「10回中6回転倒」という現象が起こる確率は約1.7%ということになります
それでは、次にｐ値を求めてみましょう（上側検定であることを忘れずに）
p値 = 6回転倒する確率 + 7回転倒する確率 +・・・+ 10回転倒する確率
この確率が0.05を下回ればS町の転倒回数が高いことが統計学的に立証できます
グラフの赤い部分の合計（p値）が0.05（有意水準）を下回れば帰無仮説が棄却されます
dbinom(x=6, size=10, prob=0.253) +
dbinom(x=7, size=10, prob=0.253) +
dbinom(x=8, size=10, prob=0.253) +
dbinom(x=9, size=10, prob=0.253) +
dbinom(x=10, size=10, prob=0.253)
Rで一気に求める場合
$\sum_{x=0}^5P(x)$ を算出して、lower.tail = FALSEとすることで $\sum_{x=6}^{10}P(x)$ となるｐ値（グラフの赤い部分）を算出します
pbinom(
q=5, size=10, prob=0.253,
lower.tail=FALSE
)
もちろん答えは同じです
#下記のように書いても同じ意味です
1 - pbinom(
q=5, size=10, prob=0.253
)
片側検定でも帰無仮説は棄却できませんでした
したがって、今回の調査ではS町の転倒割合は全国25.3％より有意に高いという結果になりました
対立仮説の勝利 true probability of success is greater than 0.253
検定ですのでサンプルサイズが変われば結果は変わります
5人中4人が転倒する確率は $P(X=4) + P(X=5)$ を求めることになります
dbinom(
x=4, size=5, prob=0.253) +
dbinom(
x=5, size=5, prob=0.253
)
Rの関数を使用します
pbinom(
3, size=5, prob=0.253,
lower.tail=F
)
やっぱり有意水準のみで検討することはよくありませんね。サンプルサイズや95％信頼区間などを考慮して結果をみんなで議論しましょう。
下側検定
ここからは例題が変わります。今度は私が想定した「10名中１名が転倒経験のあるT町」について考えます。実存しない架空の例題です。
例題 2 「85歳以上」では25.3％と4人に1人の割合で転倒することが報告されています（内閣府）。この確率25.3％を転倒が起こる確率、つまり目的としている結果が起こる確率 pとして進めます。T町に在住する85歳以上の高齢者から無作為に選んだ10名に調査した結果、1名が転倒経験ありという結果だった。T町の転倒割合は、全国割合（25.3％）より小さいと言えるか、有意水準5％で検証せよ。 ・帰無仮説：T町の転倒割合は、全国割合（25.3％）に等しい ・対立仮説：T町の転倒割合は、全国割合（25.3％）より小さい
ある比率より小さいことを実証した場合には、下側検定を実施します（ここ重要）
もうお分かりかと思いますが、以下の青色部分の検定になります
検定をしなくても確率関数のグラフからP値は0.2より大きくなることがわかります
つまりT町の転倒割合は、全国割合（25.3％）より小さいとは言えません
Rを使って求めてみます
binom.test(
x=1, n=10, p=0.253,
alternative="less"
)
下側検定なので alternative = “less”
P値=0.237で、やはり有意な結果は得られませんでした
P値
上述した検定結果のP値のみを算出する場合
pbinom(
q=1, size=10, prob=0.253,
lower.tail=T
)
転倒割合は10％のままで、サンプルサイズと転倒回数を増やしてみましょう
これは少し多めにしないとp値は小さくならないでしょう…ｎ＝40、転倒4回くらいでどうでしょうか
t <- 0:40
p <- dbinom(
x=t, size=40, prob=0.253
)
barplot(
p~t, pch=16, xlab="", ylab=""
)
これはp<0.05になっているのではないでしょうか
pbinom(
q=4, size=40, prob=0.253,
lower.tail=TRUE
)
有意になりました
サンプルサイズを増やすことが重要なのはよく分かるのですが、忙しいリハ業務のなかでサンプルサイズは簡単に増やせませんよね・・・？。だから、サンプルサイズの小さい研究では、95％信頼区間がより重要になります。
両側検定
95%信頼区間の前に両側検定について架空の例題で説明しておきます。
Rでは両側検定がデフォルト設定になっています
つまり、何も指定しなければ両側検定になります（alternative = “two.sided” ）
例題 3 U町に在住する85歳以上の高齢者から無作為に20名抽出して調査した結果、7名が転倒経験ありという結果だった。T町の転倒割合は、全国割合（25.3％）とは異なることを証明したい。 ・帰無仮説 T町の転倒割合は、全国割合（25.3％）に等しい ・対立仮説 T町の転倒割合は、全国割合（25.3％）とは異なる
Rでやってみましょう
まず、確率分布のグラフを描いておきます
t <- 0:20
p <- dbinom(
x=t, size=20, prob=0.253
)
barplot(
p~t, pch=16, xlab="", ylab=""
)
二項検定を行います
alternative = “two.sided” はデフォルトなので書く必要はありませんが、間違うこともありますので、丁寧に記載しておきましょう
binom.test(
x=8, n=20, p=0.253,
alternative="two.sided"
)
P=0.1943で帰無仮説は棄却できません
したがって全国平均とは有意に異なることが言えませんでした
このP値の求め方を勉強しましょう
両側検定のP値を求めるための3つの方法
その1
P値の２倍（上側検定または下側検定の小さい方のp値を使用する）
pbinom(
q=7, size=20, prob=0.253,
lower.tail=FALSE
)*2
$q=7$となるので注意してください
lower.tail = FALSEなのでq ≦ 7を含まない、つまりq ≧8となります （離散型）
その2
｜期待値ー観測された割合｜以上となる確率の合計
期待値：\(np = 20*0.253 = 5.06\)
\((|k – 5.06|) ≧ (|8-5.06|) =2.94 \)となる\(P(k)\)
\(k = 8, 9, …, 20\)
\(k = 2, 1, 0\)
pbinom(
q=7, size=20, prob=0.253,
lower.tail=FALSE) +
pbinom(
q=2, size=20, prob=0.253,
lower.tail=TRUE
)
lower.tail = TRUEの場合、なのでq≦3を含む、つまりq≦3となります
その3
確率が \(P(8)\) 以下となる確率の合計
Rはこの方法でP値を算しています
dbinom(
x=8, size=30, prob=0.253
)
$P(8) = 0.16$なので
$P(k | P(k) ≦ 0.16 )$ を求めます
$k = 8, 9, …, 20$
$K = 2, 1, 0$
pbinom(
q=7, size=20, prob=0.253,
lower.tail=FALSE) +
pbinom(
q=2, size=20, prob=0.253,
lower.tail=TRUE
)
今回はその2とその3は同じ値になりましたが、異なる場合もあります。
95%信頼区間
今回の調査はサンプルサイズが非常に小さかったのですが、この調査をもし100回行ったとすれば、今回の結果がどの程度の範囲に入るのでしょうか
これは確率ではありません（間違いやすいので注意が必要です）
95%信頼区間とはこの調査を100回繰り返した場合に95回が収まるであろう区間のことです
例題 「85歳以上」では25.3％と4人に1人の割合で転倒することが報告されています（内閣府）。この確率25.3％を転倒が起こる確率、つまり目的としている結果が起こる確率として進めます。S町に在住する85歳以上の高齢者から無作為に10名抽出して調査した結果、5名が転倒経験ありという結果だった。この場合の95％信頼区間を求めよ！
Rで実際に二項検定で求めてみましょう
binom.test(
5, n=10, p=0.253
)
95％信頼区間は0.187～0.813となりました
今回の割合は50％だったのですが、同じ調査を100回繰り返したとしたら、95回は18.7%～81.3%の範囲に収まるという推定結果です。
どうでしょうか？
ｐ値は有意ではなかったのですが、25.3％が信頼区間のかなり端の部分に位置していることが理解できます
有意差の有無よりも95％信頼区間を提示することが重要ですね
95％信頼区間の求め方
F分布を利用した求め方 (Clopper-Pearsonの信頼区間)
自由度
\(u1=2*(n-i+1)\)、 \(u2=2*i\)
\(l1=2*(i+1)\)、 \(l2=2*(n-i)\)
\([\frac{u2}{u2+u1F_{1-0.025}(u1, u2)}\),\(\ \ \frac{l1F_{1-0.025}(l1, l2)}{l2+l1F_{1-0.025}(l1, l2)}]\)
X～B(n, 0.253)の場合の上側P値
\(P(X≧i) = P(F_{(u1, u2)} ≧ \frac{u2(1-π)}{u1*π} )\)
\(π=\frac{u2}{u2+u1F_{1-0.025}(u1, u2)}\)
n = 10
i = 5
u1 <- 2 * (n-i+1)
u2 <- 2 * i
(π = u2 / (u2 + u1*qf(0.975, u1, u2)))
X～B(n, 0.253)の場合の下側P値
\(P(X≦i) = P(F_{(l1, l2)} ≦ \frac{l2*π}{l1*(1-π)})\)
\(π=\frac{l1F_{1-0.025}(l1, l2)}{l2+l1F_{1-0.025}(l1, l2)}\)
n = 10
i = 5
l1 <- 2 * (i+1)
l2 <- 2 * (n-i)
(π = l1 * qf(0.975, l1, l2) / (l2 + l1*qf(0.975, l1, l2)))
binom.testで出た結果と同じ値になりました
正規近似による求め方
次に正規分布に近似させる方法で求めてみましょう
近似ですのでF分布ほど厳密な方法ではありません
二項分布の期待値\(nπ\)、 分散\(nπ(1-π)\)より
\(E[\frac{i}{n}] = π\)
\(V[\frac{i}{n}] = \frac{1}{n^2}V[ i ] = \frac{1}{n^2}nπ(1-π)\)
\(X～N(π, \frac{π(1-π)}{n})\)
母比率Pの検定では帰無仮説H0=P=πとすると、この仮説のもとzは近似的に標準正規分布に従う
\(z = \frac{P-π}{\sqrt{\frac{π(1-π)}{n}}}\)
信頼区間
\(\frac{i}{n} ± z(1-0.025)\sqrt{(\frac{\frac{i}{n}(1-\frac{i}{n})}{n})}\)
Rで書くと以下のようになります
n = 10
i = 5
p1 = 0.5
p0 = 0.253
#Z統計量 N(0, 1)
z <- abs(p1-p0) / sqrt(p0*(1-p0)/n)
(p = p1 + c(-1, 1) * z * 0.975 * sqrt(p1*(1-p1)/n))
nが小さいので補正を行ってみます（Yatesの連続補正）
n = 10
il = 1 + 0.5 #i+0.5補正
iu = 1 - 0.5 #i-0.5補正
pl = 0.15
pu = 0.05
p = 0.253
#下限
z <- abs(pl-p) / sqrt(p*(1-p)/n)
pl - z * 0.975 * sqrt(pl*(1-pl)/n)
#上限
z <- abs(pu-p) / sqrt(p*(1-p)/n)
pu + z * 0.975 * sqrt(pu*(1-pu)/n)
F分布で求めた信頼区間に近づきました
サンプルサイズが小さいので信頼区間の幅もかない大きくなっています
サンプルサイズの設計
例題 「85歳以上」では25.3％と4人に1人の割合で転倒することが報告されています（内閣府）。1標本における比率の検定を行う場合、誤差を5％以内に収めるためには、どのくらいのサンプルが必要となるでしょうか？
つまり95％信頼区間が±0.05の範囲に収まるようなサンプルサイズということになります
上述した正規近似を利用して求めてみましょう
信頼区間
\(1.96\sqrt{\frac{0.253(1-0.253)}{n}}=0.05\)となるようなnを求めます
n=290程度あれば誤差の少ない信頼区間が推定できそうですが、290も無理ですよね・・・
n=290程度あれば誤差の少ない信頼区間が推定できそうですが、n=290は無理ですよね・・・
コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください