【R】母平均・母比率の差の検定まとめ

#*はじめに*

　古典的統計学において, 「信頼区間」という概念は主に推定(区間推定)と検定(仮説検定), 回帰分析の3つに登場する. 今回はこれらのうち「検定」を対象として, 母平均の差の検定と母比率の差の検定を確認する. 

　まず改めて統計的仮説検定とは, 母集団分布の母数に関する仮説を標本から検証する統計学的方法の1つである. `R`では`t.test()`関数などを用いることで1行のコードで検定が実行できるものの中身が*Black Box*になりがちだ. 

　そこで今回は統計量*t*や*p*値をできるだけ手計算し, 帰無仮説の分布を可視化することでより直感的な理解を目指す.

#*母平均の差の検定*

　母平均の差の検定における検定統計量*(t or z)*は下記の通り, 検証条件によって求める式が変わる. 

|  母平均の差の検定  |  標本の群数  |  標本の対応  |  母分散の等分散性  |  t値  |
|:--:|:--:|:--:|:--:|:--:|
|  One-Sample t test  |  1群  |  -  |  等分散である  |  $t=\frac{\bar{X}-\mu}{\sqrt{\frac{s^2}{n}}}$  |
|  Paired t test  |  2群  |  対応あり  |  等分散である  |  $t=\frac{\bar{X_D}-\mu}{\sqrt{\frac{s_D^2}{n}}}$  |
|  Student's test  |  2群  |  対応なし  |  等分散である  |  $t=\frac{\bar{X_a}-\bar{X_b}}{\sqrt{s_{ab}^2}\sqrt{\frac{1}{n_a}+\frac{1}{n_b}}}$  |
|  Welch test  |  2群  |  対応なし  |  等分散でない  |  $t=\frac{\bar{X_a}-\bar{X_b}}{\sqrt{\frac{s_a^2}{n_a}+\frac{s_b^2}{n_b}}}$  |

　※本記事で式中に登場する*s*は, 母分散が既知の場合は標準偏差*σ*, 母分散が未知の場合は不偏標準偏差*U*を指す
　
　以降では, 代表的なものを例題を通して確認していく. 

#*1標本のt検定(One-Sample t-test)*

　1標本の*t*検定は, ある意味区間推定とほぼ変わらない. *p*値もそうだが, 帰無仮説で差がないとする特定の数値(多くの場合は*0*)が, 設定した区間推定の上限下限に含まれているかを確認する. 

　今回は, 正規分布に従う*web*ページ*A*の滞在時間の例を用いて, 帰無仮説を以下として片側検定する.
　　

```math
H_0: \mu\geq0\\ 
H_1: \mu<0\\ 
```

　また, 1群のt検定における*t*統計量は, 以下で定義される.　　
　　

```math
t=\frac{\bar{X}-\mu}{\sqrt{\frac{s^2}{n}}}\\ 
```

　まずは, *t*値を*by hand*で計算する.

```R
#データ生成
data<-rnorm(10,30,5)

#帰無仮説よりμは0
mu<-0

#平均値
x_hat<-mean(data)

#不偏分散
uv<-var(data)

#サンプルサイズ
n<-length(data)

#自由度
df<-n-1

#t値の推計
t<-(x_hat-mu)/(sqrt(uv/n))
t
```
- *output: 36.397183465115*

　`t.test()`メソッドで, *p*値と$\bar{X}$の区間推定を確認する.

```R
t.test(before, after, paired=TRUE, alternative="less", conf.level=0.95)
```
>	One Sample t-test

>data:  data
t = 36.397, df = 9, p-value = 4.418e-11
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 28.08303 31.80520
sample estimates:
mean of x 
 29.94411 

　*p値<0.05*より, 帰無仮説を棄却する. よって母平均*μ=0*とは言えない結果となった.

#*対応のある2標本のt検定(Paired T-test)*

　「対応のある」とは, 同一サンプルから抽出された2群のデータに対する検定を指す. 対応のある2標本のt検定では, 基本的に2群の差が*0*かどうかを検定する. つまり, *前後差=0*を帰無仮説とする1標本問題として検定する．

　今回は, 正規分布に従う*web*ページ*A*のデザイン変更前後の滞在時間の差の例を用いて, 帰無仮説を以下として片側検定する.
　　

```math
H_0: \bar{X_D}\geq\mu_D\\ 
H_1: \bar{X_D}<\mu_D\\ 
```
　　
　対応のある2標本の平均値の差の検定における*t*統計量は, 以下で定義される.　
　　

```math
t=\frac{\bar{X_D}-\mu_D}{\sqrt{\frac{s_D^2}{n}}}\\ 
\bar{X_D}=\frac{1}{n}\sum_{i=1}^n (x_{Di})\\ 
s_D^2=\frac{1}{n}\sum_{i=1}^n (x_{Di}-\bar{x_D})^2\;\;or\;\;s_D^2=\frac{1}{n-1}\sum_{i=1}^n (x_{Di}-\bar{x_D})^2\\ 
```
　　
　まずは, *t*値を*by hand*で計算する.

```R
#データ生成
before<-c(32,45,43,65,76,54)
after<-c(42,55,73,85,56,64)

#差分数列の生成
d<-before-after

#帰無仮説よりμは0
mu<-0

#差の平均
xd_hat<-mean(d)

#差の標準偏差
sd<-var(d)

#サンプルサイズ
n<-length(d)

#自由度
df<-n-1

#t値の推計
t=(xd_hat-mu)/sqrt(sd/n)
t
```
- *output: -1.4638501094228*

　次に, *p*値を計算&可視化して有意水準α(棄却域)と比較する.

```R
#棄却域の定義
t_lower<-qt(0.05, df)

#有意水準の出力
alpha<-pt(t_lower, df)
alpha

#p値
p<-pt(t,df)
p
```
- *output: 0.05* 
- *output: 0.101555331860027*

```R
#棄却域の定義
t_lower<-qt(0.05, df)

options(repr.plot.width=14, repr.plot.height=8)
curve(dt(x, df),-5,5,type="l", col="lightpink", lwd=10, main="t-distribution: df=5") 
abline(v=qt(p=0.05, df), col="salmon", lwd=4, lty=5)
abline(v=t, col="skyblue", lwd=4, lty=1)
curve(dt(x, df),-5,t,type="h", col="skyblue", lwd=4, add=T) 
curve(dt(x, df),-5,qt(p=0.05,df),type="h", col="salmon", lwd=4, add=T)
```
![pairedt.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/457537/8cf25e38-fd11-4c05-cbf1-6a9374efe275.png)

　*p値>0.05*であるようだ. `t.test()`メソッドで, *t*値と*p*値を確認する.

```R
t.test(before, after, paired=TRUE, alternative="less", conf.level=0.95)
```
> 	Paired t-test

>data:  before and after
t = -1.4639, df = 5, p-value = 0.1016
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 3.765401
sample estimates:
mean of the differences 
                    -10 

　*p値>0.05*より, 帰無仮説を採択し, 母平均*μ*は*0*とは言えない結果となった.

#*対応のないt検定(Unpaired T-test)*

#*等分散性がある場合: Studentのt検定(Student's T-test)*

　対応のない2標本の平均値の差の検定において, 2標本の母分散が等しいということが既知の場合, スタンダードな*Student*の*t*検定を用いる. その際, F検定による等分散に対する検定を行うことで判断する.

　今回は, 正規分布に従うフランス人とイタリア人の平均身長の例を用いて, 帰無仮説を以下として片側検定する.
　　

```math
H_0: \bar{X_a}\geqq\bar{X_b}\\ 
H_1: \bar{X_a}<\bar{X_b}\\ 
```

　等分散性がある2標本の平均値の差の検定における*t*統計量は, 以下で定義される.　　　
　　

```math
t=\frac{\bar{X_a}-\bar{X_b}}{\sqrt{s_{ab}^2}\sqrt{\frac{1}{n_a}+\frac{1}{n_b}}}\\
s_{ab}^2=\frac{\sum_{i=1}^{n_a}(X_{ai}-\bar{X_a})^2+\sum_{i=1}^{n_b}(X_{bi}-\bar{X_b})^2}{n_a+n_b-2}\\ 
```
　　
　まずは, *F*検定により等分散性を確認する.

```R
#Studentのt検定
#データ生成
france<-rnorm(8, 160, 4.5)
italy<-rnorm(10, 156, 5)

#平均値
x_hat_france<-mean(france)
x_hat_italy<-mean(italy)

#偏差平方和
vs_france<-sum((france-x_hat_france)^2)
vs_italy<-sum((italy-x_hat_italy)^2)

#サンプルサイズ
n_france<-length(france)
n_italy<-length(italy)

uv_france<-var(france)
uv_italy<-var(italy)

#不偏分散同士の比率
f_value<-uv_france/uv_italy
f_value
```
- *output: 0.37272*
![stf.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/457537/c4205f68-ff3c-2803-2951-03970dac57f9.png)

```R
#等分散性の確認
var.test(x=france, y=italy)
```
>	F test to compare two variances

>data:  france and italy
F = 0.37272, num df = 7, denom df = 9, p-value = 0.2068
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.0888051 1.7977054
sample estimates:
ratio of variances 
         0.3727192 

　*p値>0.05*より, 帰無仮説を採択し, 等分散性があるとして進める. 次に, *t*値を*by hand*で計算する.

```R
#自由度: 
df<-n_france+n_italy-2

#t値の推計
student_t<-(x_hat_france-x_hat_italy)/sqrt(((vs_france+vs_italy)/(df))*(1/n_france+1/n_italy))
student_t
```
- *output: 2.22021538587324*

　次に, *p*値を計算&可視化して有意水準α(棄却域)と比較する.

```R
#棄却域の定義
t_upper<-qt(0.95, df)

#有意水準の出力
alpha<-1-pt(qt(0.95, df), df)
alpha

#p値
p<-1-pt(student_t, df)
p
```
- *output: 0.05* 
- *output: 0.020598429743329*

```R
#棄却域の定義
t_upper<-qt(0.95, df)

options(repr.plot.width=14, repr.plot.height=8)
curve(dt(x, df),-5,5,type="l", col="lightpink", lwd=10, main="t-distribution: df=16") 
abline(v=qt(p=0.95, df), col="salmon", lwd=4, lty=5)
abline(v=student_t, col="skyblue", lwd=4, lty=1)
curve(dt(x, df),qt(p=0.95,df),5,type="h", col="salmon", lwd=4, add=T)
curve(dt(x, df),student_t,5,type="h", col="skyblue", lwd=4, add=T) 
```
![student_t.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/457537/64213511-62d8-6d9d-6ee5-61f4f97aa0b7.png)

　*p値<0.05*であるようだ. `t.test()`メソッドで, *t*値と*p*値を確認する.

```R
t.test(x=france, y=italy, var.equal=T, paired=F, alternative="greater", conf.level=0.95)
```
>	Two Sample t-test

>data:  france and italy
t = 2.2202, df = 16, p-value = 0.0206
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 0.6358715       Inf
sample estimates:
mean of x mean of y 
 160.6547  157.6784 

　*p値<0.05*より, 帰無仮説を棄却し, *2*標本の母平均に差がありそうだという結果となった.

#*等分散性がない場合: Welchのt検定(Welch's T-test)*

　一方で, 2標本の母分散は等しいと言えない場合に使われるのが*Welch*のの*t*検定である. ただし, *2*段階検定の問題から2標本のt検定を行う場合には等分散性を問わず, *Welch's T-test*を行うべきだという主張もある.

　今回は, 正規分布に従うフランス人とスペイン人の平均身長の例を用いて, 帰無仮説を以下として片側検定する.
　　

```math
H_0: \bar{X_a}\geqq\bar{X_b}\\ 
H_1: \bar{X_a}<\bar{X_b}\\ 
```
　　
　等分散性のない2標本の差の検定における*t*統計量は, 以下で定義される.　
　　

```math
t=\frac{\bar{X_a}-\bar{X_b}}{\sqrt{\frac{s_a^2}{n_a}+\frac{s_b^2}{n_b}}}\\ 
```
　　
　まずは, *F*検定により等分散性を確認する.

```R
#データ生成
france<-rnorm(8, 160, 3)
spain<-rnorm(11, 156, 7)

#平均値
x_hat_france<-mean(france)
x_hat_spain<-mean(spain)

#不偏分散
uv_france<-var(france)
uv_spain<-var(spain)

#サンプルサイズ
n_france<-length(france)
n_spain<-length(spain)

#不偏分散同士の比率
f_value<-uv_france/uv_spain
f_value
```
- *output: 0.068597*
![f.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/457537/f38c5675-8b8b-308e-666f-b5b13d8b3efb.png)



```R
#等分散性の確認
var.test(x=france, y=spain)
```
>	F test to compare two variances

>data:  france and spain
F = 0.068597, num df = 7, denom df = 10, p-value = 0.001791
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.01736702 0.32659675
sample estimates:
ratio of variances 
        0.06859667 

　*p値<0.05*より, 帰無仮説を棄却し, 等分散性がないとして進める. 次に, *t*値を*by hand*で計算する.

```R
#自由度: Welch–Satterthwaite equationで算出(省略)
df<-11.825

#t値の推計
welch_t<-(x_hat_france-x_hat_spain)/sqrt(uv_france/n_france+uv_spain/n_spain)
welch_t
```
- *output: 0.9721899010868*

　次に, *p*値を計算&可視化して有意水準α(棄却域)と比較する.

```R
#棄却域の定義
t_upper<-qt(0.95, df)

#有意水準の出力
alpha<-1-pt(qt(0.95, df), df)
alpha

#p値
p<-1-pt(welch_t, df)
p
```
- *output: 0.05* 
- *output: 0.175211697240612*


```R
#棄却域の定義
t_lower<-qt(0.05, df)

options(repr.plot.width=14, repr.plot.height=8)
curve(dt(x, df),-5,5,type="l", col="lightpink", lwd=10, main="t-distribution: df=5") 
abline(v=qt(p=0.05, df), col="salmon", lwd=4, lty=5)
abline(v=t, col="skyblue", lwd=4, lty=1)
curve(dt(x, df),-5,t,type="h", col="skyblue", lwd=4, add=T) 
curve(dt(x, df),-5,qt(p=0.05,df),type="h", col="salmon", lwd=4, add=T)
```
![Welcht.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/457537/18ae872b-7f3f-effa-bc10-8574c4b1b064.png)

　*p値>0.05*であるようだ. `t.test()`メソッドで, *t*値と*p*値を確認する.

```R
t.test(x=france, y=spain, var.equal=F, paired=F, alternative="greater", conf.level=0.95)
```
>	Welch Two Sample t-test

>data:  france and spain
t = 0.97219, df = 11.825, p-value = 0.1752
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 -2.01141      Inf
sample estimates:
mean of x mean of y 
 158.7778  156.3704 

　*p値>0.05*より, 帰無仮説を採択し, *2*標本の母平均には差があるとは言えなさそうだという結果となった.

#*母比率の差の検定*

　母比率の差の検定では, 2つのグループのある比率が等しいかどうかを検定する. またサンプルサイズnが十分に大きいとき, 二項分布が正規分布*N(0, 1)*に近似できることと同様に, 検定統計量にも標準正規分布に従う統計量*z*を用いる. 
　　
　今回は, 正規分布に従う*web*ページ*A*の滞在時間の例を用いて, 帰無仮説を以下として検定する.
　　

```math
H_0: \hat{p_a}=\hat{p_b}\\ 
H_1: \hat{p_a}\neq\hat{p_b}\\ 
```
　　
　また母比率の差の検定における*t*統計量は, 以下で定義される.　なお帰無仮説が「2標本の母比率に差がない」という場合には, 分母に標本比率をプールした統合比率*(pooled proportion)*を用いることを注意したい.
　　　　

```math
z=\frac{\hat{p_a}-\hat{p_b}}{\sqrt{\hat{p}(1-\hat{p})\Bigl(\frac{1}{n_a}+\frac{1}{n_b}\Bigr)}}\\ 
\hat{p}=\frac{n_a\hat{p_a}+n_b\hat{p_b}}{n_a+n_b}
```
　　
　まずは, *z*値を*by hand*で計算する.

```R
#サンプル
new<-c(150, 10000)
old<-c(200, 12000)

#それぞれのpの期待値
p_hat_new<-new[1]/new[2]
p_hat_old<-old[1]/old[2]

#サンプルサイズ
n_new<-new[2]
n_old<-old[2]

#統合比率
p_hat_pooled<-(n_new*p_hat_new+n_old*p_hat_old)/(n_new+n_old)

#z値の推計
z<-(p_hat_new-p_hat_old)/sqrt(p_hat_pooled*(1-p_hat_pooled)*(1/n_new+1/n_old))
z
```
- *output: -0.983756784574025*

　次に, *p*値を計算&可視化して有意水準α(棄却域)と比較する.

```R
#棄却域の定義
z_lower<-qnorm(0.05, 0, 1)

#有意水準の出力
alpha<-pnorm(qnorm(p=0.05, 0, 1), 0, 1)
alpha

#p値
p<-pnorm(z, 0, 1)
p
```
- *output: 0.05* 
- *output: 0.162617556154074*


```R
#棄却域の定義
z_lower<-qnorm(0.05, 0, 1)

options(repr.plot.width=14, repr.plot.height=8)
curve(dnorm(x, 0, 1),-5,5,type="l", col="forestgreen", lwd=10, main="z-distribution") 
abline(v=qnorm(p=0.05, 0, 1), col="salmon", lwd=4, lty=5)
abline(v=z, col="skyblue", lwd=4, lty=1)
curve(dnorm(x, 0, 1),-5,z,type="h", col="skyblue", lwd=4, add=T) 
curve(dnorm(x, 0, 1),-5,qnorm(p=0.05,0, 1),type="h", col="salmon", lwd=4, add=T)
```
![z.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/457537/6420cc89-0f4d-bdb7-77e4-60eb0b33e5af.png)

　*p値>0.05*であるようだ. `prop.test()`メソッドで, *t*値と*p*値を確認する.

```R
# correct: a logical indicating whether Yates’ continuity correction should be applied where possible.
prop.test(c(new[1],old[1]), c(new[2],old[2]), p=NULL, alternative="less", correct=F, conf.level=0.95)
```
>	2-sample test for equality of proportions without continuity
	correction

>data:  c(new[1], old[1]) out of c(new[2], old[2])
X-squared = 0.96778, df = 1, p-value = 0.1626
alternative hypothesis: less
95 percent confidence interval:
 -1.000000000  0.001106871
sample estimates:
    prop 1     prop 2 
0.01500000 0.01666667 

　　*p値>0.05*より, 帰無仮説を採択し, *2*標本の母比率に差はなさそうだという結果となった. また先ほど手計算した*z*値と上記のカイ二乗値が, また*p*値が一致していることが確認できる. 

#*さいごに*

　以上で, 母平均・母比率の差の検定を終える. 今回は代表的な佐野検定だけを取り上げたが, 母分散が既知/未知などを気にすると無数に存在する. 次回はベイズ推定による差の検定をまとめる.

###◎参考文献

- https://toukeigaku-jouhou.info/2018/07/07/hypothesis-testing-list/
- https://bellcurve.jp/statistics/course/18227.html
- http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/6_t_test1.html
- https://biolab.sakura.ne.jp/z-test-proportion.html
- https://biolab.sakura.ne.jp/z-test-proportion.html