Rで行う統計検定まとめ

# はじめに
Rでつかえる統計検定(`t.test`, `var.test`, `cor.test`, `chisq.test`)の使用方法のまとめを行った。自身で作ったデータセットで検証しながら各コマンドの出力について解説する。

# この記事で用いるデータセット
## 量的変数のデータ
簡単に計算ができるデータセットが欲しかったので、[こちらの記事](https://qiita.com/sakabe/items/df7fac64e6ffbb50848c)を参考にRで実装してみる。

```R:コード
### 身長とBMIから体重を作成してデータフレーム化する
m_height <- rnorm(n=1000,mean=171.7,sd=6.6)
f_height <- rnorm(n=1000,mean=158.3,sd=5.7)

m_bmi <- rnorm(n=1000,mean=23.12,sd=4.24)
f_bmi <- rnorm(n=1000,mean=20.82,sd=3.42)

m_weight <- (m_height * 0.01)^2 * m_bmi
f_weight <- (f_height * 0.01)^2 * f_bmi

male <- data.frame(HEIGHT=m_height, WEIGHT=m_weight, SEX="male") #男性のみのデータ
female <- data.frame(HEIGHT=f_height, WEIGHT=f_weight, SEX="female") #女性のみのデータ

dat <- rbind(male, female) #結合したデータ

### できたデータをプロットしてみる
plot(0, 0, type = "n", xlim = c(min(m_weight, f_weight), max(m_weight, f_weight)), ylim = c(min(m_height, f_height), max(m_height, f_height)), xlab="weight [kg]", ylab="height [cm]" )
points(m_weight, m_height, col = "blue", pch = 1)
points(f_weight, f_height, col = "red", pch = 3)

legend("bottomright", legend = c("male", "female"), pch = c(1, 3), col = c("blue", "red"))
```
![R_WH_dataset.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/f27a3476-c624-af64-dc00-f32d3fd1eb24.png)

身長、体重で相関があるデータセットを作ることができた。

## 質的変数のデータ
質的変数も検定できるよう性別と血液型からなるデータセットを生成したい。参考として赤ちゃんの血液型の発現割合を[こちらのサイト](https://keisan.casio.jp/exec/system/1374719499)から調べて、人口比率を[こちらのサイト](https://www.stat.go.jp/data/jinsui/topics/topi1201.html)から調べた。それぞれの割合になるよう乱数を判定してデータセットを生成した。

```R:コード
myfunc.mf <- function(x) {
  a <- c()
  for(i in 1:x){
    tmp <- runif(1)
    if (tmp < 0.513986) a <- c(a, "male")
    else                a <- c(a, "female")
  }
  return(a)
}

myfunc.bt <- function(x) {
  a <- c()
  for(i in 1:x){
    tmp <- runif(1)
    if (tmp < 0.39)                    a <- c(a, "A")
    else if (0.39 <= tmp & tmp < 0.61) a <- c(a, "B")
    else if (0.61 <= tmp & tmp < 0.90) a <- c(a, "O")
    else                               a <- c(a, "AB")
  }
  return(a)
}

mf <- myfunc.mf(2000)
bt <- myfunc.bt(2000)

dat2 <- data.frame(SEX=mf, BLOOD=bt)

table(dat2)
```

```:出力
        BLOOD
SEX        A  AB   B   O
  female 381 117 197 284
  male   372 111 258 280
```

二つの変数の間に連関が無い質的変数のデータセットが完成した。


# 統計検定
## `t.test`を用いた検定。
### 平均値の検定（One Sample t-test）

正規母集団からの無作為標本である場合、標本平均 $\bar{X}$ 、標本の不偏分散の正の平方根 $\hat{\sigma}$ 、標本の個数 $n$ を用いて

```math
t = \frac{\bar{X}-\mu}{\hat{\sigma}/\sqrt{n}}
```

で計算される統計量 $t$ は帰無仮説のもとで自由度 $df=n-1$ の $t$ 分布に従う。この性質を利用して統計検定を行うRコマンドが`t.test`である。

男性の身長データについて$\mu = 170.0$ の検定を行ってみる。

```R:コード
t.test(male["HEIGHT"], mu=170.0)
```

```:出力
	One Sample t-test

data:  male["HEIGHT"]
t = 8.2108, df = 999, p-value = 6.743e-16
alternative hypothesis: true mean is not equal to 170
95 percent confidence interval:
 171.2686 172.0655
sample estimates:
mean of x 
 171.6671 
```
データ、統計量・自由度・p値、対立仮説、95%の信頼区間、標本から予測される平均値が出力される。もちろん母集団の平均は男性 171.7 cmでデータを作成しているので $p=0.05$ の水準で帰無仮説（平均値は170である）は棄却される。

### 2つの平均値の比較(Two Sample t-test)
二つの標本の母集団の平均が一致するかどうかを調べたい時も`t.test`を使うことができる。検定統計量は、

```math
t = \frac{\bar{X_1}-\bar{X_2}}{\sqrt{\frac{(n_1-1)\hat{\sigma_1}^2+(n_2-1)\hat{\sigma_2}^2}{n_1 + n_2 -2}\left(\frac{1}{n_1}+\frac{1}{n_2} \right)}}
```
となり、帰無仮説の下で自由度 $df=n_1+n_2-2$ の $t$ 分布に従う。ただしこの検定が使える条件は２つの母集団の分散が等質であることが要求される。[分散の等質性の検定](###-分散の等質性の検定)に詳細を記述する。

`t.test`を使って、男性の身長と女性の身長の母集団の平均値が等しいのか検定をしていく。

```R:コード
t.test(male["HEIGHT"], female["HEIGHT"], var.equal = TRUE)
```

```:出力

	Two Sample t-test

data:  male["HEIGHT"] and female["HEIGHT"]
t = 48.889, df = 1998, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 12.57924 13.63064
sample estimates:
mean of x mean of y 
 171.6671  158.5621 
```
データ、統計量・自由度・p値、対立仮説、95%の信頼区間、標本から予測される平均値が出力される。もちろん母集団の平均は男女で変えて行っているので $p=0.05$ の水準で帰無仮説（平均値が男女で一致する）は棄却される。

### Welch の検定
[分散の等質性の検定](###-分散の等質性の検定)で分散の等質性が得られなかった場合に、二つの平均値を比較するには Welch の検定を行う。これも `t.test` で実現できる。先ほどのコマンドのオプションを`var.equal = FLASE`とすればよいだけである。

```R:コード
t.test(male["HEIGHT"], female["HEIGHT"], var.equal = FALSE)
```

```:出力
	Welch Two Sample t-test

data:  male["HEIGHT"] and female["HEIGHT"]
t = 48.889, df = 1955.5, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 12.57924 13.63065
sample estimates:
mean of x mean of y 
 171.6671  158.5621 
```
データ、統計量・自由度・p値、対立仮説、95%の信頼区間、標本から予測される平均値が出力される。もちろん母集団の平均は男女で変えて行っているので $p=0.05$ の水準で帰無仮説（平均値が男女で一致する）は棄却される。


## `var.test`を用いた検定
### 分散の等質性の検定
通常は2群の平均値の比較を行う前に分散が等質であるかの検定を行う。そのためのRコマンドが`var.test`である。先ほど $t$ 検定で用いた男女の身長データに対して分散が等質であるかを調べてみる。`var.test`もデータフレームの入力を受け付けないので、列番号で指定する。

```R:コード
var.test(male[,1], female[,1])
```

```:出力
	F test to compare two variances

data:  male[, 1] and female[, 1]
F = 1.3458, num df = 999, denom df = 999, p-value = 2.822e-06
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 1.188720 1.523576
sample estimates:
ratio of variances 
          1.345773 
```
データ、統計量・自由度・p値、対立仮説、95%の信頼区間、標本から予測される分散の比が出力される。もちろん母集団の分散も男女で変えて行っているので $p=0.05$ の水準で帰無仮説（分散の比が１になる）は棄却される。

## `cor.test`を用いた検定
### 無相関検定（Pearson's product-moment correlation）
相関係数 $r$ を用いて計算される統計量

```math
t = \frac{r\sqrt{n-2}}{\sqrt{1-r^2}}
```
は、帰無仮説のもとで自由度 $df=n-2$ の $t$ 分布に従う。この性質を利用して相関の有無を検定するRコマンドが`cor.test`である。

`cor.test`を用いて、データの身長と体重に相関があるのか調べてみる。データはdat["WEIGHT"] などで入れるとデータフレーム扱いになってしまうので、列番号で指定した。

```R:コード
cor.test(dat[,1], dat[,2]) 
```

```:出力
	Pearson's product-moment correlation

data:  dat[, 1] and dat[, 2]
t = 36.58, df = 1998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6063250 0.6588662
sample estimates:
      cor 
0.6333248 
```
データ、統計量・自由度・p値、対立仮説、95%の信頼区間、標本から予測される相関係数が出力される。身長に基づきBMIから体重を作成しているので二つの数字に相関はある。検定では、 $p=0.05$ の水準で帰無仮説（身長と体重は無相関である）は棄却されている。

###その他の無相関検定
`cor.test`のオプションで、method="spearman", method="kendall"を選べば、それぞれ「スピアマンの順位相関係数の無相関検定」、「ケンドールの順位相関係数の無相関検定」を行うことができる。因みにデフォルトで設定されているのはタイトルの通り、method="pearson"で「ピアソンの積率相関係数の無相関検定」である。

## `chisq.test`を用いた検定
###独立性の検定
質的変数についてクロス集計表を取った場合、その観測度数 $O_i$ 、期待度数 $E_i$ から計算される統計量である

```math
\chi^2 = \sum_{i}^k \frac{\left(O_i - E_i\right)^2}{E_i}
```

は帰無仮説のもとで $\chi^2$ 分布に従う。ただし自由度はクロス集計表の「（行の数）-1」×「（列の数）-1」である。この性質を利用して独立性の検定を行うのが`chisq.test`である。

`chisq.test`を用いて性別と血液型に連関があるのか調べてみる。

```R:コード
tab2 <- table(dat2)
chisq.test(tab2)
```

```:出力
	Pearson's Chi-squared test

data:  tab2
X-squared = 7.5932, df = 3, p-value = 0.05521
```
データ、統計量・自由度・p値、が出力される。検定では、 $p=0.05$ の水準で帰無仮説（性別と血液型は連関が無い）は棄却されないので、それぞれ独立なデータであると言える。


# 参考文献

- [【機械学習】練習用データセットの作成【身長・体重】](https://qiita.com/sakabe/items/df7fac64e6ffbb50848c)
- https://keisan.casio.jp/exec/system/1374719499
- https://www.stat.go.jp/data/jinsui/topics/topi1201.html
- [Rによるやさしい統計学](https://shop.ohmsha.co.jp/shopdetail/000000001781/)