> codeiqのRの問題の解説編なんだが，え？ あの問題って引っかけじゃなくて，本当に検定の多重性を無視していたんだ，と驚いた。うーん，「東京大学で医療分野を中心に統計学を研究」しているなら，疫学・生物統計学教室と無関係ではないだろう経歴の人なのに，これってまずいんでは？
ちょっと，的外れでは？
私も別の記事でちょっと触れた記事なんだけど，まあ，色々問題はあるかも知れないが，「検定」についてシミュレーションをしてαエラー，βエラーについて確認してみたというだけのことではないでしょうか。中澤さんも，2種の条件下の2つの検定，あるいはそれぞれの状況下で企図された1000回の標本調査の多重性を考えているわけではないと思いますけどね？
英語版を見たら robust test だったよ。
統計学で robust といえば，「頑健」のこと。ロバスト検定といってもいいし，頑健検定といってもよいけど「耐久」とはひどいね。
締め切りにより，解答案公表も解禁でしょうから
http://gihyo.jp/dev/serial/01/codeiq/0007
「第7回 t検定による問題解決，Rで実践できますか？～データサイエンティストの統計学─倉橋一成からの問題」について
// 10万人の集団A，B，Cの群があります。これらの集団の身長の「平均値」が等しいか等しくないか，全員の身長を測定せずに判断したいと思います。それぞれの集団から100人ずつサンプリングしてt検定したらどうなるか，Rで計算してみましょう。
//
// 準備として，以下のコードを実行し，10万人分の身長データを3群作成してください。それぞれの群は以下のように想定しています。
//
// A群とB群⇒ 平均値が170cmの集団
// C群⇒ 平均値が175cmの集団
//
// # A群のデータ生成
// set.seed(1)
// heightA <- 170 + 10*rnorm(100000)
// # B群のデータ生成
// set.seed(2)
// heightB <- 170 + 10*rnorm(100000)
// # C群のデータ生成
// set.seed(3)
// heightC <- 175 + 10*rnorm(100000)
まずこの段階では，rnorm 関数が mean, sd 引数を取ることを指摘しておこう。つまり，
heightA <- rnorm(100000, mean=170, sd=10)
ということである（引き数名は省略可能だが，説明のためには書いておく方が無難）。
また，このような一連の手続きをプログラムにするにあたっては，何回も出てくる定数は一箇所で定義しておき，二回目以降はその変数名を使うというのは定石である。このような部分は求める解答ではないが，出題者も留意する方が善いと思われる。もっとも，回答者がこのようなコメントをつけるかどうかを見るというなら別であるが。
ということで，この部分までについては，以下のように書く方がスマート。（なお，英数字と日本語は半角空白で区切った方が読みやすいと思うので）。
n <- 100000 ################ 母集団のサイズ
# A 群のデータ生成
set.seed(1)
heightA <- rnorm(n, mean=170, sd=10)
# B 群のデータ生成
set.seed(2)
heightB <- rnorm(n, mean=170, sd=10)
# C 群のデータ生成
set.seed(3)
heightC <- rnorm(n, mean=175, sd=10)
// 準備ができたら，下記の4つの問題に解答してください。
// 【問1】 A，B，C群からそれぞれ100人ずつサンプリングしたいです。下記の空欄(1)～(3)にあてはまるコードを埋めてください。
//
// # A群のサンプリング
// set.seed(11)
// heightASmple <- (1)
// # B群のサンプリング
// set.seed(12)
// heightBSmple <- (2)
// # C群のサンプリング
// set.seed(13)
// heightCSmple <- (3)
heightASmple というのは，変数名として紛らわしい。Smple というのは綴り間違いではなく出題者がわざと縮めたのだろうが，プログラムする側は Sample と綴りがち。実行してみればエラーが出るのでわかるけど，いらいらしないためにも変な変数名にしない方がよい。また，複数の単語からなる変数名は "." などで区切るのがよいと思われる。ということで，以下では heightA.Sample などのように変更する。
解答
m <- 100 ################ サンプルサイズ（標本の大きさ）
# A 群のサンプリング
set.seed(11)
heightA.Sample <- sample(heightA, m)
# B 群のサンプリング
set.seed(12)
heightB.Sample <- sample(heightB, m)
# C 群のサンプリング
set.seed(13)
heightC.Sample <- sample(heightC, m)
なお，母集団の 100000 人から 100 人をサンプリングするというのは，現実のシミュレーションとしては妥当である。しかし，以降の問題では本質的なものではなく，また，このような条件でサンプリングするとシミュレーションが面倒になる（つまり，結果は同じようになるが，無限母集団として扱った方が簡単ということ）。つまり，
set.seed(11)
heightA.Sample <- rnorm(m, 170, 10)
heightB.Sample <- rnorm(m, 170, 10)
heightC.Sample <- rnorm(m, 175, 10)
で差し支えないということ。なお，3 群のサンプリングそれぞれの前に set.seed する必要はない。気になるならば，A 群からサンプリングする前に 1 回行えば，それで十分。
// 【問2】A群 vs B群，A群 vc C群の身長の平均値をt検定し，5％有意であったか計算してください。計算に使ったRのコードと計算結果を提出してください。
母標準偏差が 10 ということがどの程度強制力を持つのか分からないが，とりあえず var.equal=TRUE としておく。
解答
t.test(heightA.Sample, heightB.Sample, var.equal=TRUE)$p.value
t.test(heightA.Sample, heightC.Sample, var.equal=TRUE)$p.value
// 【問3】問1と問2の操作を1000回繰り返したときのP値を求めてください。下記の空欄(1)と(2)にあてはまるコードを埋めてください。
//
// PvalAB <- NA
// PvalAC <- NA
//
// # 問1のP値を算出
// for( i in 1:1000 ){
//
// (1)
//
// }
//
// # 問2のP値を算出
// for( i in 1:1000 ){
//
// (2)
//
// }
出題文が間違えている。問 1 では P 値は求めていない。
ということで，問 1 と問 2 を別のループにできるわけもない。for ループを使いなさいという強制と理解しておく。
A 群と B 群の比較と A 群と C 群の比較を別々の for ループでやりなさいということかもしれないが，heightA.Sample は別々に 2 回やる必要はないし，無駄なので，for ループは 1 つでよい。
そのほかに，いくつか不適切な点がある。
(a) PvalAB <- NA; PvalAC <- NA のように初期化して，for ループの中で結果を数珠つなぎにしていく方法がだめな方法であるということはよく指摘される。最初から要素数が分かっているときには（今の場合は 1000，なおかつ，その数値を使うのではなく，その数値が入っている変数名を使って），loop <- 1000; PvalAB <- numeric(loop) のようにすべしと。
(b) PvalAB とか PvalAC も，Pval.AB とか Pval.AC のようにわかりやすい（？）名前にする。
(c) 定数を明示的に書かないことにすると，i in 1:1000 は i in 1:loop となるので，もし loop <- 0 などと変なことをするといけないので，i in seq_len(loop) とせよというのも，よくされる指摘である。
解答
loop <- 1000 ################ シミュレーションの試行回数
Pval.AB <- Pval.AC <- numeric(loop)
for (i in seq_len(loop)) {
heightA.Sample <- sample(heightA, m)
heightB.Sample <- sample(heightB, m)
heightC.Sample <- sample(heightC, m)
Pval.AB[i] <- t.test(heightA.Sample, heightB.Sample, var.equal=TRUE)$p.value
Pval.AC[i] <- t.test(heightA.Sample, heightC.Sample, var.equal=TRUE)$p.value
}
// 【問4】問3で作成した1000回分のP値を使って，αエラーとβエラーに関してどのようなことが言えるか記述してください。計算に使ったRのコードと計算結果を提出してください。
解答
################ 群 A と 群 B の平均値の差の検定で P 値 <= 0.05 となる割合
mean(Pval.AB <= 0.05)
################ 群 A と 群 c の平均値の差の検定で P 値 <= 0.05 となる割合
mean(Pval.AC <= 0.05)
念のため，mean(Pval.AC <= 0.05) は sum(Pval.AC <= 0.05) / loop ということなので，P 値の平均ということではない。
全部まとめて解答
シミュレーションというのは，条件を変えて行うところに意味があるので，今までのプログラムを関数化して，変化させる条件を引数で指定できるようにする。
sim <- function(n=100000, m=100, loop=1000, mean.AB=170, mean.C=175, SD=10, seed=1) {
set.seed(seed)
heightA <- rnorm(n, mean.AB, SD)
heightB <- rnorm(n, mean.AB, SD)
heightC <- rnorm(n, mean.C, SD)
# set.seed(11)
heightA.Sample <- sample(heightA, m)
# set.seed(12)
heightB.Sample <- sample(heightB, m)
# set.seed(13)
heightC.Sample <- sample(heightC, m)
t.test(heightA.Sample, heightB.Sample, var.equal=TRUE)$p.value
t.test(heightA.Sample, heightC.Sample, var.equal=TRUE)$p.value
Pval.AB <- Pval.AC <- numeric(loop)
for (i in seq_len(loop)) {
heightA.Sample <- sample(heightA, m)
heightB.Sample <- sample(heightB, m)
heightC.Sample <- sample(heightC, m)
Pval.AB[i] <- t.test(heightA.Sample, heightB.Sample, var.equal=TRUE)$p.value
Pval.AC[i] <- t.test(heightA.Sample, heightC.Sample, var.equal=TRUE)$p.value
}
return(c(mean(Pval.AB <= 0.05), mean(Pval.AC <= 0.05)))
}
実行例
> sim()
[1] 0.052 0.933
> sim(loop=10000)
[1] 0.0465 0.9448
> sim(n=1000000, loop=10000)
[1] 0.0473 0.9429
解釈
A 群と B 群の母平均は等しいので，そこから取り出される標本の標本平均も等しいと期待される。しかし，偶然により標本平均が等しいという帰無仮説が棄却される確率がαエラー（今の場合 0.05）である。loop が大きくなれば，loop 回の検定で p.value <= 0.05 になる割合は 0.05 から大きく外れるようなことは少なくなる。微妙な言い回しをしたのは，これらは確率的なものなので，loop が大きくなれば一様に 0.05 に近づいていくということではないからである。αエラーというのは，帰無仮説が正しい（A 群と B 群の母平均は等しい）にもかかわらず，帰無仮説を棄却する誤りである。有意水準αで検定すると，この誤りは確率的にαに等しい。
A 群と C 群の母平均は等しくないので，そこから取り出される標本の標本平均も等しくないと期待される。その確率は検出力（1 - β）である。βエラーとは，帰無仮説が誤っているのに帰無仮説を棄却できない確率である。どれくらいの確率でそのようなことが起きるかは，母平均の差と母標準偏差によって異なる。今の場合（平均値の差が 5，母標準偏差が 10 ）は，その確率が 1-0.05 = 0.95 のように見えるかもしれないが，いつも 0.95 になるのではないということに注意が必要である。
二群の，サンプルサイズ，平均値，標準偏差（計 6 個の統計量）と有意水準αが決まれば，そのときの検出力 power =（1 - β）は計算（power.t.test）により求めることができる。βエラー = 1-検出力
> sim()
[1] 0.052 0.933
> power.t.test(n=100, delta=5, sd=10, sig.level=0.05)
Two-sample t test power calculation
n = 100
delta = 5
sd = 10
sig.level = 0.05
power = 0.9404272
alternative = two.sided
NOTE: n is number in *each* group
========================================================
> sim(mean.C=180)
[1] 0.052 1.000
> power.t.test(n=100, delta=10, sd=10, sig.level=0.05)
Two-sample t test power calculation
n = 100
delta = 10
sd = 10
sig.level = 0.05
power = 0.9999998
alternative = two.sided
NOTE: n is number in *each* group
> sim(SD=8)
[1] 0.052 0.992
========================================================
> power.t.test(n=100, delta=5, sd=8, sig.level=0.05)
Two-sample t test power calculation
n = 100
delta = 5
sd = 8
sig.level = 0.05
power = 0.992614
alternative = two.sided
NOTE: n is number in *each* group
========================================================
★★ 出題の枠にとらわれない解答 ★★
sim <- function(n=100000, m=100, loop=1000, mean.AB=170, mean.C=175, SD=10) {
heightA <- rnorm(n, mean.AB, SD)
heightB <- rnorm(n, mean.AB, SD)
heightC <- rnorm(n, mean.C, SD)
heightA.Sample <- sample(heightA, m)
heightB.Sample <- sample(heightB, m)
heightC.Sample <- sample(heightC, m)
t.test(heightA.Sample, heightB.Sample, var.equal=TRUE)$p.value
t.test(heightA.Sample, heightC.Sample, var.equal=TRUE)$p.value
Pval.AB <- Pval.AC <- numeric(loop)
for (i in seq_len(loop)) {
heightA.Sample <- sample(heightA, m)
heightB.Sample <- sample(heightB, m)
heightC.Sample <- sample(heightC, m)
Pval.AB[i] <- t.test(heightA.Sample, heightB.Sample, var.equal=TRUE)$p.value
Pval.AC[i] <- t.test(heightA.Sample, heightC.Sample, var.equal=TRUE)$p.value
}
return(c(mean(Pval.AB <= 0.05), mean(Pval.AC <= 0.05)))
}
replicate を使おう。
sim2 <- function(m=100, loop=1000, mean.AB=170, mean.C=175, SD=10) {
Pval <- replicate(loop, {
A <- rnorm(m, mean.AB, SD)
c(t.test(A, rnorm(m, mean.AB, SD), var.equal=TRUE)$p.value,
t.test(A, rnorm(m, mean.C, SD), var.equal=TRUE)$p.value)
})
return(rowMeans(Pval <= 0.05))
}
> system.time(print(sim()))
[1] 0.056 0.929
ユーザ システム 経過
0.998 0.014 1.009
> system.time(print(sim2()))
[1] 0.049 0.957
ユーザ システム 経過
0.608 0.005 0.610
ベンチマークで見てみる
> library(rbenchmark)
> benchmark(sim(), sim2(), replications=10)
test replications elapsed relative user.self sys.self user.child sys.child
1 sim() 10 10.063 1.638 10.035 0.078 0 0
2 sim2() 10 6.145 1.000 6.308 0.060 0 0
やはり，sim2 の方が速い。
