ブートストラップ入門

```

1 Limits of Parametric Statistics 既知の確率分布を使った推測の限界

平均値の区間推定やクロス表の独立性の検定など、社会学でのオーソドックスな分析の多くは、t 値やカイ二乗値といった統計量が、一定の条件下で t分布やカイ二乗分布のような確率分布に近似することを利用した分析法である。こういった分析法は比較的計算が簡単で多くの研究者に知られているという大きなメリットがあるが、万能ではない。例えば、中央値やジニ係数の区間推定をしたい場合には、既知の確率分布を使った推定法は知られていない（と思う）。

このような問題はしばしば生じる。例えば、

• 期待値が非常に小さい場合のクロス表の独立性の検定、
• 変数が正規分布していない場合の相関係数の区間推定、
• 残差分散が不均一な場合の分散分析、
• サンプルサイズが小さく、変数が正規分布していない場合の平均値の差の検定

といった例が考えられる。このような問題に対する対処法はいくつかあるが、ここではブートストラップという方法を紹介する。以下の議論は統計学的には曖昧な（もしかしたら多少不正確な）記述も含まれているが、社会学者がその基本的な考え方を習得する上では、これぐらいがちょうどいいと判断している。

2 Bootstrapping ブートストラップ

ブートストラップとは、得られたサンプルから何回もリサンプリングを行い、それらのリサンプルから統計量を求めることで、推測を行うような方法のことである。

2.1 Resampling リサンプリング

ふつう統計的なデータ分析は、母集団から抽出されたサンプルを使って行われる。この得られたサンプルの中からさらに抽出（サンプリング）を行うことをリサンプリング、リサンプリングによって得られたサンプルをリサンプルと呼ぶ。なお以下では表現を具体化しわかりやすくするために、研究対象は人間で人間の集団が母集団だという前提で話を進めるが、研究対象は何でも構わない。また特に断りのない限り、抽出（サンプリング）は無作為抽出のことである。

例えば、ある大学の中から 20 人の学生を無作為に抽出し、先月のアルバイト収入をたずねたところ、以下の様なデータが得られたとしよう。

x1 <- c(0, 0, 0, 0, 2000, 4000, 5000, 9000, 10000, 10500, 11000, 11500, 12000, 14000, 17000, 20000, 20500, 21000, 30000, 50000)
x1

##  [1]     0     0     0     0  2000  4000  5000  9000 10000 10500 11000
## [12] 11500 12000 14000 17000 20000 20500 21000 30000 50000

これが標本（サンプル）である。この標本の中からさらに 5人を無作為に選んだら、以下の様になった。

set.seed(1986) # 乱数を発生させるが、結果に再現性を持たせるため。引数の数字は何でもいい
x2 <- sample(x1, 5)  # sample(ベクトル, n) で、ベクトルの中から n 人だけ無作為に非復元抽出する 
x2

## [1] 10500     0 11000  4000 50000

このような 2段階目のサンプリングがリサンプリングであり、標本の標本をリサンプルと呼ぶ。

2.2 Sampling with/without Replacement 復元抽出と非復元抽出

標本抽出には、復元抽出と非復元抽出の二種類がある。非復元抽出とは、一度サンプルに選ばれた人は母集団から除外され、二度と抽出されないように行うサンプリングであるが、復元抽出とは、一度サンプルに選ばれた人も母集団から除外されず、複数回同じ人が選ばれる可能性があるようなサンプリングのことである。以下は上の 5人のリサンプルの中からさらに 4人非復元抽出、そして復元抽出した実例。

# x2 の5人の中から 4人を抽出
set.seed(1986)
x3.1 <- sample(x2, 4) # 非復元抽出。同じ人が複数回選ばれることはない
x3.1

## [1] 11000 10500     0  4000

set.seed(1986)
x3.2 <- sample(x2, 4, replace = TRUE) # 復元抽出。同じ人が複数回選ばれる可能性あり
x3.2

## [1] 11000     0  4000     0

上の非復元抽出では

• 一人目の 11000 円の人が抽出された後、二人目は、残りの 4人 (10500, 0, 4000, 50000) の中から抽出される。
• そして二人目に 10500 円の人が抽出された後、三人目は、残りの 3人 (0, 4000, 50000) の中から抽出される。
• 以下同様にして4人目まで抽出、

というプロセスが繰り返されている。いっぽう復元抽出の場合、

• 一人目の 11000 円の人が抽出された後、二人目は、もう一度、5人全体の中から抽出される。
• 三人目以降も、5人全体の中から抽出されるので、ぐうぜん同じ人が2回以上選ばれることもある。この例では 0円の人が2回選ばれている、

というわけである。

2.3 Steps of Bootstrapping ブートストラップの計算

ブートストラップ区間推定の最も簡単なやり方は、以下のとおりである。

1. サンプルと同じサイズのリサンプルを 復元 抽出する。
2. リサンプルから区間推定や検定をしたい統計量（例えば中央値や相関係数）を計算する。
3. 上の 2つのプロセスを何回も（例えば 1000回）繰り返し、統計量をすべて記録する。このベクトルを統計量の分布と呼んでおく。
4. 統計量の分布の \(\alpha/2\)パーセンタイルと \(100-\alpha/2\)パーセンタイルが \(100-\alpha\) % 区間推定の下限値と上限値になる。例えば \(\alpha=5\)、つまり 95% 信頼区間を求めるには、得られた統計量の分布の 2.5 パーセンタイルと、97.5 パーセンタイルを求めればよい。

2.3.1 Calculation: 計算例

上で使った x1 と名付けたサンプルからブートストラップで中央値の区間推定をしてみよう。ブートストラップの前にまず、このサンプルから単純に中央値を計算すると、

x1

##  [1]     0     0     0     0  2000  4000  5000  9000 10000 10500 11000
## [12] 11500 12000 14000 17000 20000 20500 21000 30000 50000

median(x1)

## [1] 10750

である。次に 1000 回繰り返しのブートストラップしてみよう。まず 1 回めは以下のとおり。

median1000 <- rep(NA, 1000) # 1000回分の中央値を記録するための空のベクトルを用意
x4.1 <- sample(x1, 20, replace = TRUE)
sort(x4.1)  # 中央値がわかりやすいように昇順で並べ替え。 sort(ベクトル) はベクトルを昇順で並べ替える 

##  [1]     0     0     0     0  2000  2000  2000  4000 10000 10000 11500
## [12] 11500 11500 14000 14000 17000 20000 20000 21000 30000

(median1000[1] <- median(x4.1)) # 計算した中央値を median1000 の 1番めの要素として記録

## [1] 10750

2回めの計算は以下のとおり。

x4.2 <- sample(x1, 20, replace = TRUE)
sort(x4.2) 

##  [1]     0     0  2000  2000  4000  4000 10500 11500 11500 11500 14000
## [12] 17000 20000 20000 20000 20000 20500 30000 50000 50000

(median1000[2] <- median(x4.2))

## [1] 12750

以下同様にあと 998 回同じプロセスを繰り返せばよい。しかし、上のようなスクリプトを 998回も書くのは面倒なので、以下のように 1000回分まとめてループを回すとよい。

for(i in 1:1000){
  x4 <- sample(x1, 20, replace = TRUE)
  median1000[i] <- median(x4)
}
median1000[1:30] # 1000個中央値を表示させるのは煩雑なので、最初の30回分だけ表示

##  [1] 10500 11500  5000 10750 10500 10500 11500 10750 10000  9500  8000
## [12]  9750 12750 11500 11250 10250 12000 10500 10750 11250  9500 12500
## [23] 10500 10500 17000 10250 10500  2000 10500  6500

hist(median1000) # 得られた統計量（1000個の中央値）の分布を確認

最後に得られた中央値の分布の 2.5パーセンタイルと97.5パーセンタイルを求めると、

q1 <- quantile(median1000, c(0.025, 0.975)) # quantile(データベクトル, パーセンタイル位置を指定するベクトル) で指定したパーセンタイルを返す
q1

##  2.5% 97.5% 
##  4500 15500

となり、この大学の学生のアルバイト収入の中央値は 95% の確率で 4500～15500 のあいだの値をとることがわかる。

2.3.2 apply() を使った繰り返し計算

余談だが、for() を使うよりも、apply() を使ったほうが計算が早い。上の計算は、以下のように書いてもよい。

x5 <- sample(x1, 20 * 1000, replace=TRUE) # 20人 x 1000回分のデータを一度に作る
x5[1:40]  # 最初の 40個だけ（つまり２回分のリサンプルの）中身を見る

##  [1]  5000  9000     0     0 11000     0  5000 21000  5000     0 11500
## [12]     0  2000 30000 20500  9000  5000     0  2000  9000 50000     0
## [23] 21000     0  9000     0     0 20500 30000 12000 12000 10500     0
## [34]  9000     0 20500  9000     0 11000 20500

x5 <- matrix(x5, 20, 1000) # 上のベクトルを 20行 x 1000列の行列に変換（各列が１回分のリサンプル）
x5[,1:4] # 最初の４列分だけ（つまり４回分のリサンプルの）中身を見る

##        [,1]  [,2]  [,3]  [,4]
##  [1,]  5000 50000 17000 20500
##  [2,]  9000     0 20000  2000
##  [3,]     0 21000 17000  4000
##  [4,]     0     0  5000 17000
##  [5,] 11000  9000     0  4000
##  [6,]     0     0  5000 17000
##  [7,]  5000     0  9000 20000
##  [8,] 21000 20500 20000  9000
##  [9,]  5000 30000     0  9000
## [10,]     0 12000  9000  5000
## [11,] 11500 12000  9000 30000
## [12,]     0 10500  2000 14000
## [13,]  2000     0 20000  4000
## [14,] 30000  9000 21000     0
## [15,] 20500     0     0     0
## [16,]  9000 20500 50000 12000
## [17,]  5000  9000     0 20500
## [18,]     0     0 20500 21000
## [19,]  2000 11000     0  9000
## [20,]  9000 20500  4000 11500

median1000 <- apply(x5, 2, median) # 列ごとに中央値を計算
median1000[1:10]  # 最初の 10個だけ（つまり 10回分のリサンプルの）中央値を見る

##  [1]  5000  9750  9000 10250  6500  4500  9000 12000 11000 13000

quantile(median1000, c(0.025, 0.975)) # 1000個の中央値の 2.5パーセンタイルと 97.5パーセンタイルを求める

##     2.5%    97.5% 
##  4500.00 15506.25

2.3.3 simpleboot: 専用の関数で計算

上のやり方でまったく問題ないが、R にはブートストラップを行う専用の関数がいくつかある。ここでは simpleboot パッケージの one.boot() 関数と boot.ci() 関数を使った計算例を紹介しておこう。

one.boot(サンプルのベクトル, 統計量を計算する関数, 繰り返し数)

boot.ci(one.boot()の結果, type=“perc”) # type=“perc” を省略すると、何通りかの区間推定結果を返す

library(simpleboot)

## Loading required package: boot
## Simple Bootstrap Routines (1.1-3 2008-04-30)

boot.x1 <- one.boot(x1, median, 2000)
boot.ci(boot.x1, type="perc")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot.x1, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 4500, 15500 )  
## Calculations and Intervals on Original Scale

2.4 Warnings: 以下の点に注意

• このようなブートストラップの結果は無作為のリサンプリングにもとづいているので、毎回異なる結果が得られることもある。
• 繰り返し数は多いほど推測の精度は高まる（少なくとも低まらない）が、たくさん繰り返すほど時間がかかる。95% 信頼区間を求めるならば、1000回以上は繰り返したい。
• サンプルが母集団からかけ離れたものであればブートストラップの結果も歪むので、サンプルは母集団に近似している方がよい（無作為抽出でサンプルサイズが大きいほうがいい）。

2.5 Multiple Variables: 複数の変数を扱う場合の計算例

扱う変数が複数でも、考え方は中央値のような一変数だけを扱う統計量と同じだが、プログラミングに多少の工夫が必要である。以下のような例で計算してみよう。

学生のアルバイト収入がどの程度時間的に変化するのか知りたいとしよう。そこで、先ほどの 20人の大学生に、三ヶ月後に再び一ヶ月間のアルバイト収入をたずねた。その結果が以下のとおりである。

y1 <- c(0, 10000, 0, 20000, 11000, 0, 5000, 0, 9500, 5000,
        8500, 11500, 12000, 0, 15000, 21000, 13000, 19000, 60000, 30000)
d1 <- data.frame(x1, y1)
d1

##       x1    y1
## 1      0     0
## 2      0 10000
## 3      0     0
## 4      0 20000
## 5   2000 11000
## 6   4000     0
## 7   5000  5000
## 8   9000     0
## 9  10000  9500
## 10 10500  5000
## 11 11000  8500
## 12 11500 11500
## 13 12000 12000
## 14 14000     0
## 15 17000 15000
## 16 20000 21000
## 17 20500 13000
## 18 21000 19000
## 19 30000 60000
## 20 50000 30000

cor(d1)   # cor(データフレーム) でデータフレーム内の変数間の相関行列を返す

##          x1       y1
## x1 1.000000 0.654781
## y1 0.654781 1.000000

アルバイト収入は正規分布しているとは思えないので（下の図1を参照）、ふつうの相関係数の検定は適切ではない。

par(mfrow=c(1,2))
hist(y1,col="red")
hist(x1, col = "blue")

図1：2時点のアルバイト収入のヒストグラム

そこで、この 2 時点のアルバイト収入の相関係数の 95% 信頼区間をブートストラップで求めてみよう。まずこの 20 人のサンプルの中から復元抽出でリサンプルを得るが、データが 2 変数からなるので、このデータフレームの何行目の人をリサンプルするかまず決めるという手間が余分にかかる。

set.seed(1986)
no.resample <- sample(1:20, 20, replace=TRUE) # 1～20の整数の中から20回復元抽出してリサンプルするケースを決定
no.resample

##  [1] 10  5 13  7 12 14  4 15 14 19  5  6 16 12  5  9  2 18  3  2

これでリサンプルすべき人が d1 の何行目の人か決まった。これに従ってリサンプリングすると以下のようになる。

d2 <- d1[no.resample, ]
d2

##         x1    y1
## 10   10500  5000
## 5     2000 11000
## 13   12000 12000
## 7     5000  5000
## 12   11500 11500
## 14   14000     0
## 4        0 20000
## 15   17000 15000
## 14.1 14000     0
## 19   30000 60000
## 5.1   2000 11000
## 6     4000     0
## 16   20000 21000
## 12.1 11500 11500
## 5.2   2000 11000
## 9    10000  9500
## 2        0 10000
## 18   21000 19000
## 3        0     0
## 2.1      0 10000

各行の名前に “14.1”, “5.1” といった小数付きものもあるが、これはもとのサンプル (d1) の 14行目や 5行目が 複数回選ばれていて、2回めに選ばれた時にそのような名前がつけられている。これは行や列の名前はすべて互いに異なっていなければならない（同じ名前は2回使えない）からである。

このリサンプルで x1 と y1 の相関係数を計算すればよい。

cor(d2)

##           x1        y1
## x1 1.0000000 0.6051769
## y1 0.6051769 1.0000000

この作業を 1000回続ければいい。ここから先は 1変数のブートストラップと同じである。

cor1000 <- rep(NA, 1000)
for(i in 1:1000){
     no.resample <- sample(1:20, 20, replace=TRUE)
     d2 <- d1[no.resample, ]
     cor1000[i] <- cor(d2)[1,2]
}
quantile(cor1000, c(0.025, 0.975))

##      2.5%     97.5% 
## 0.2975676 0.8740514

2.5.1 pairs.boot()

simpleboot パッケージを使って相関係数のブートストラップをするならば、pairs.boot() 関数が使える。one.boot() ではベクトルしかデータとして引数に指定することが出来ないので、これで相関係数の計算はできないのである。

pairs.boot(変数1, 変数2, 統計量を計算する関数, 繰り返し数)

boot.x1y1 <- pairs.boot(x1, y1, cor, 1000)
boot.ci(boot.x1y1, type="perc")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot.x1y1, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 0.2518,  0.8837 )  
## Calculations and Intervals on Original Scale

ただし、pairs.boot() は 2変数しか使わない統計量にしか使えないので、もっと多くの変数を扱う場合には役に立たない。OLS で重回帰分析する場合には、lm.boot(), 一般化線形モデルなら car パッケージの Boot() 関数が使えるし、その他にも様々なブートストラップ対応の関数があるが、専用の関数がない場合は、上で紹介したように自分でループを回すか、boot パッケージの boot() 関数を使う必要がある。boot() は非常に汎用性が高い優れた関数だが、ややプログラミングが難しい。

2.5.2 by() で相関係数のブートストラップ

さらに余談だが、apply() 関数族である、by() を使っても相関係数のブートストラップは可能である。

cor1000 <- rep(NA, 1000)
no.resample <- sample(1:20, 20 * 1000, replace=TRUE) # 2万回ぶん復元抽出する行を一度に決める
no.resample[1:30] # 最初の30回分の結果だけ示す

##  [1] 13  7 20 14  9 11 18  2  9 11 10 14 17 14 17  1  8 16 11  6 20 20 10
## [24]  8  1  5 16 17 14 19

d2 <- d1[no.resample, ] # ２万回分のリサンプルをひとまとめにする
head(d2, 30) # 最初の30行だけ表示

##         x1    y1
## 13   12000 12000
## 7     5000  5000
## 20   50000 30000
## 14   14000     0
## 9    10000  9500
## 11   11000  8500
## 18   21000 19000
## 2        0 10000
## 9.1  10000  9500
## 11.1 11000  8500
## 10   10500  5000
## 14.1 14000     0
## 17   20500 13000
## 14.2 14000     0
## 17.1 20500 13000
## 1        0     0
## 8     9000     0
## 16   20000 21000
## 11.2 11000  8500
## 6     4000     0
## 20.1 50000 30000
## 20.2 50000 30000
## 10.1 10500  5000
## 8.1   9000     0
## 1.1      0     0
## 5     2000 11000
## 16.1 20000 21000
## 17.2 20500 13000
## 14.3 14000     0
## 19   30000 60000

sample.group <- gl(1000, 20) # 2万ケースのりサンプルを20人ずつの1000グループに分けるためのグループの名前（番号）を作成
head(cbind(d2, sample.group), 30) # 最初の30行だけ表示

##         x1    y1 sample.group
## 13   12000 12000            1
## 7     5000  5000            1
## 20   50000 30000            1
## 14   14000     0            1
## 9    10000  9500            1
## 11   11000  8500            1
## 18   21000 19000            1
## 2        0 10000            1
## 9.1  10000  9500            1
## 11.1 11000  8500            1
## 10   10500  5000            1
## 14.1 14000     0            1
## 17   20500 13000            1
## 14.2 14000     0            1
## 17.1 20500 13000            1
## 1        0     0            1
## 8     9000     0            1
## 16   20000 21000            1
## 11.2 11000  8500            1
## 6     4000     0            1
## 20.1 50000 30000            2
## 20.2 50000 30000            2
## 10.1 10500  5000            2
## 8.1   9000     0            2
## 1.1      0     0            2
## 5     2000 11000            2
## 16.1 20000 21000            2
## 17.2 20500 13000            2
## 14.3 14000     0            2
## 19   30000 60000            2

corr <- function(x) cor(x)[1,2] # この分析では数値を一つだけ返す関数のほうが後のデータ処理が簡単なので、相関行列の1行2列目だけを返す関数を作成
cor1000 <- by(d2, sample.group, corr) # d2 を sample.group のグループで分類して corr() をすべてのグループに関して計算
quantile(cor1000, c(0.025, 0.975)) 

##      2.5%     97.5% 
## 0.2652717 0.8730850

3 Test: 検定したい場合

例えば、帰無仮説として

\(H_0\): 相関係数 \(=0\)

のような仮説をおき、これを棄却できるか検定したい場合は、上で計算した信頼区間に 0 が含まれているかどうかチェックすればよい。含まれていれば帰無仮説は棄却できないし、含まれていなければ、帰無仮説は棄却できる。上の例では、95% 信頼区間を計算し、信頼区間の中に 0 は含まれていないので、5% 水準で帰無仮説は棄却できる。

ただし、ブートストラップ検定というと、違うやり方の方が一般的な印象があるので、注意。しかし、そのやり方は回帰分析のようなモデルで使うのが普通なので、社会調査士 D 科目の範疇を超えている気がするので、割愛する。

4 Variety of Bootstrapping: ブートストラップの多様性

ブートストラップにもいろいろな方法がある。ここで紹介した方法は完全ブートストラップ (full bootstrapping) と呼ばれるものだが、セミパラメトリック・ブートストラップとか、パラメトリック・ブートストラップといった方法もある。

セミパラメトリック・ブートストラップとは、回帰分析で残差が正規分布していないとかんがえられる場合、分析結果から得られた残差ベクトルをリサンプリングして、回帰モデルから従属変数のリサンプルを作るような方法である。パラメトリック・ブートストラップとは、残差をリサンプリングするかわりに、推定された残差分散の正規分布に従う乱数を発生させることで、リサンプルの従属変数を作る方法である。

また、今回は信頼区間を計算するときに単純にパーセンタイルをとったが、他にもやり方は何通りかある。