母相関係数＝０の検定をブートストラップで

「ぶーつすとらっぷ」じゃあ，ないよ。

どこぞやらで，だれかが，R を使ってブートストラップによる母相関係数の分布を推定するというブログを書いておったが，ちっとも R らしくなかったので，書いてみた。

具体例のデータは，サンプルサイズ 20 の以下のようなもの。
x = c(7.8, 25.1, 36.1, 27.1, 24.5, 11.4, 27.6, 12.3, 5.1, 4, 26.3, 21.8, 35.4, 18.5, 31.6, 7, 10.8, 24.2, 25.4, 17.9)
y = c(29.3, 13.5, 28.9, 38.8, 30.4, 8.1, 19.7, 25.9, 2.5, 16.2, 10.3, 21.2, 27.6, 12.5, 29.1, 8, 8.6, 29.7, 25.9, 13.8)

このデータから，x, y それぞれで，重複を許して 20 個のデータを抽出して相関係数を計算するというのを多数回繰り返す。

x から重複を許して 20 個のデータをリサンプリングするのは
sample(x, replace=TRUE)

y から重複を許して 20 個のデータをリサンプリングするのは
sample(x, replace=TRUE)

その相関係数を計算するのは
cor(sample(x, replace=TRUE), sample(y, replace=TRUE))

それを，たとえば 50000 回繰り返すのは
sapply(1:50000, function(i) cor(sample(x, replace=TRUE), sample(y, replace=TRUE)))

その結果を要素数 50000 のベクトル result に格納するのは
result = sapply(1:50000, function(i) cor(sample(x, replace=TRUE),
                                         sample(y, replace=TRUE)))

そのヒストグラムは
hist(result)

リサンプリングデータにおける相関係数の絶対値が，観察されたデータでの相関係数 cor(x, y) より大きい割合は，
mean(abs(result) > cor(x, y))

sum(abs(result) > cor(x, y)) / 50000 としたいかもしれないが，それは mean だ！

リサンプリングによらず，観察されたデータでの母相関係数＝0の検定（両側）の p 値は
cor.test(x, y)$p.value

sample 関数は乱数発生を伴うので，以下のプログラムでの結果は，毎回異なるが，mean(abs(result) > cor(x, y)) と cor.test(x, y)$p.value の結果はほぼ等しくなる。

> x = c(7.8, 25.1, 36.1, 27.1, 24.5, 11.4, 27.6, 12.3, 5.1, 4, 26.3, 21.8, 35.4, 18.5, 31.6, 7, 10.8, 24.2, 25.4, 17.9)
> length(sample(x, replace=TRUE))
[1] 20
> length(sample(x, replace=TRUE))
[1] 20
> x = c(7.8, 25.1, 36.1, 27.1, 24.5, 11.4, 27.6, 12.3, 5.1, 4, 26.3, 21.8, 35.4, 18.5, 31.6, 7, 10.8, 24.2, 25.4, 17.9)
> y = c(29.3, 13.5, 28.9, 38.8, 30.4, 8.1, 19.7, 25.9, 2.5, 16.2, 10.3, 21.2, 27.6, 12.5, 29.1, 8, 8.6, 29.7, 25.9, 13.8)
> result = sapply(1:50000, function(i) cor(sample(x, replace=TRUE), sample(y, replace=TRUE)))
> hist(result)
> mean(abs(result) > cor(x, y))
[1] 0.00994
> cor.test(x, y)$p.value
[1] 0.009260503