Rを使っての統計検定をggplotで図示しながらやる
手を動かしながら学ぶ ビジネスに活かすデータマイニングという本を買った．
- 作者: 尾崎隆
- 出版社/メーカー: 技術評論社
- 発売日: 2014/08/22
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (1件) を見る
この本を読みながらちょっとずつ試している．
手を動かしながらできるのであれこれ考えて，試し，失敗しながら学べるスタイルがとても良いと思う．
しかし，データが図で表されていないことがあり，直感的な理解が苦しい箇所があるので，ggplot2を使い図示しながら勉強しているのでメモ程度に残しておく．
今回は3章の検定を行った．
仮説検定
2つのデータを比べて異なる結果が出たのか，それとも誤差の範囲なのかを明らかにする方法
t検定
> sleep
extra group ID
1 0.7 1 1
2 -1.6 1 2
3 -0.2 1 3
4 -1.2 1 4
5 -0.1 1 5
6 3.4 1 6
7 3.7 1 7
8 0.8 1 8
9 0.0 1 9
10 2.0 1 10
11 1.9 2 1
12 0.8 2 2
13 1.1 2 3
14 0.1 2 4
15 -0.1 2 5
16 4.4 2 6
17 5.5 2 7
18 1.6 2 8
19 4.6 2 9
20 3.4 2 10
> p <- ggplot(sleep, aes(x=group, y=extra)) + geom_point()
> p
被験者の睡眠時間，睡眠薬の有無から睡眠薬の効果があったのかを検定する
> t.test(extra ~ group, data=sleep, paired=T)
Paired t-test
data: extra by group
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.4598858 -0.7001142
sample estimates:
mean of the differences
-1.58
tが負の値 = 後者が大きい
p-value 確率 0.05より小さければ差がある
boxplotの結果
> boxplot(extra~group, data=sleep)
チルダ~でデータを区切って出力できるらしい
χ二乗検定，フィッシャーの正確確率検定
ABテストなどどちらに効果があったのかを判定する
t検定は多くのレコード（睡眠薬ならば被験者）を使って行うが，χ二乗検定では2*2の行列を使うことが一般的
> ab1
[,1] [,2]
[1,] 25 117
[2,] 16 32
> chisq.test(ab1)
Pearson's Chi-squared test with Yates' continuity correction
data: ab1
X-squared = 4.3556, df = 1, p-value = 0.03689
p値が0.05より小さければ差があると判断できる
ただし，差があるとわかるだけなのでどちらが大きい/小さいはわからず，目視しないといけない
> fisher.test(ab1)
Fisher's Exact Test for Count Data
data: ab1
p-value = 0.0267
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.1929297 0.9695458
sample estimates:
odds ratio
0.4295063
フィッシャーの正確確率検定はサンプル数が少ないときに有効である
> p <- ggplot(mdf, aes(x=AB, y=value, fill=variable))
> p + geom_bar()
この図を作りたい場合
data.frameは
> mdf
AB variable value
1 A CONV 25
2 B CONV 16
3 A NCONV 117
4 B NCONV 32
でないといけない
> df
CONV NCONV AB
1 25 117 A
2 16 32 B
をいじって図を書こうとしてもCONV とNCONVを 重ねあわせることができないので無理である
変換するにはライブラリreshapeをインストールし，meltを使えば良い
> install.packages("reshape")
> library(reshape)
> mdf <- melt(df)
Using AB as id variables
> mdf
AB variable value
1 A CONV 25
2 B CONV 16
3 A NCONV 117
4 B NCONV 32
順位和検定
順位を使って行う検定
ばらつきがあるデータの場合，t検定は使いにくい
しかし，順位であればばらつきがあっても大丈夫なので検定できる
> r
A B
1 146 157
2 162 117
3 152 116
4 149 137
5 157 123
6 132 143
7 154 122
8 156 133
9 150 113
10 143 129
11 171 134
12 157 134
13 159 119
14 167 127
15 157 121
16 135 137
17 160 126
18 148 117
19 156 155
20 155 134
21 142 121
22 163 146
23 156 128
24 152 127
25 148 124
26 154 146
27 157 133
28 152 137
29 165 121
30 157 129
> mr = melt(r)
Using as id variables
> mr
variable value
1 A 146
2 A 162
3 A 152
4 A 149
5 A 157
6 A 132
7 A 154
8 A 156
9 A 150
10 A 143
11 A 171
12 A 157
13 A 159
14 A 167
15 A 157
16 A 135
17 A 160
18 A 148
19 A 156
20 A 155
21 A 142
22 A 163
23 A 156
24 A 152
25 A 148
26 A 154
27 A 157
28 A 152
29 A 165
30 A 157
31 B 157
32 B 117
33 B 116
34 B 137
35 B 123
36 B 143
37 B 122
38 B 133
39 B 113
40 B 129
41 B 134
42 B 134
43 B 119
44 B 127
45 B 121
46 B 137
47 B 126
48 B 117
49 B 155
50 B 134
51 B 121
52 B 146
53 B 128
54 B 127
55 B 124
56 B 146
57 B 133
58 B 137
59 B 121
60 B 129
# 図示
> p = ggplot(mr, aes(x=variable, y=value))
> p + geom_point()
見た感じAが大きいように見える
これを検定で確かめる
> wilcox.test(r$A, r$B)
Wilcoxon rank sum test with continuity correction
data: r$A and r$B
W = 841.5, p-value = 7.204e-09
alternative hypothesis: true location shift is not equal to 0
警告メッセージ:
In wilcox.test.default(r$A, r$B) :
タイがあるため、正確な p 値を計算することができません
p値が充分に小さいのでABに差異があることが判明する
### 外れ値
> r2 <- r
> r2[1,1] <- 50
> r2[1,2] <- 400
> r2
A B
1 50 400
2 162 117
3 152 116
4 149 137
5 157 123
6 132 143
7 154 122
8 156 133
9 150 113
10 143 129
11 171 134
12 157 134
13 159 119
14 167 127
15 157 121
16 135 137
17 160 126
18 148 117
19 156 155
20 155 134
21 142 121
22 163 146
23 156 128
24 152 127
25 148 124
26 154 146
27 157 133
28 152 137
29 165 121
30 157 129
> mr2 = melt(r2)
Using as id variables
> p = ggplot(mr2, aes(x=variable, y=value))
> p + geom_point()
> t.test(r2$A, r2$B)
Welch Two Sample t-test
data: r2$A and r2$B
t = 1.2286, df = 38.583, p-value = 0.2267
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-7.914541 32.381208
sample estimates:
mean of x mean of y
150.5333 138.3000
> wilcox.test(r2$A, r2$B)
Wilcoxon rank sum test with continuity correction
data: r2$A and r2$B
W = 805, p-value = 1.566e-07
alternative hypothesis: true location shift is not equal to 0
警告メッセージ:
In wilcox.test.default(r2$A, r2$B) :
タイがあるため、正確な p 値を計算することができません
順位和検定ならばうまくいく事がわかる
