ブートストラップとは、小さいサンプルを復元抽出で何度もサンプリングして、疑似的に何度も試行したことにして、その結果から推定値の直接的な95％の範囲を求めたりする方法。
ブートストラップを使ったロジスティック回帰分析の方法をご紹介。
- ブートストラップとは？
- ブートストラップ ロジスティック回帰のスクリプト
- 使うデータとロジスティック回帰モデル
- ブートストラップする
- ブートストラップ95％信頼区間を出力する
- まとめ
- EZR公式マニュアル
ブートストラップとは？
ブートストラップとは、小さいサンプルを復元抽出で何度もサンプリングして、疑似的に何度も試行したことにして、その結果から推定値の直接的な95％の範囲を求めたりする方法。
復元抽出とは、たまたま複数回同じデータがとられてもOKとする方法。
オリジナルのサンプルとサイズ（データの数）が同じサンプルをたくさん作る。
このほうが現実世界に近い。
こうすると、データの構成がばらつく。
オリジナルのデータから、似たようなサンプルをたくさん作って、小さいサイズのオリジナルデータでも信頼性高い、区間推定を行える方法だ。
ちなみに、ブートストラップは靴紐（ブーツのひもかな？）という意味で、アメリカのことわざで、湖に散歩に行って溺れた人が、自分の靴紐につかまって生還するというものがあり、そこから、困ったことがあっても自力で何とかするという意味だそうである。
そこから、小さいサイズのサンプルでも、そのサンプルだけで、自力で何とかするという意味で、ブートストラップ法という名前がついたのだとか。
新谷先生が、発明者のEfron先生から直々に聞いたとのことで、間違いないようだ。
24秒あたりから説明あり。
ググってみたところ、こんなブログ記事があった。
ブートストラップとは、ブーツのかかとについている、小さなループの部分のことを言っているみたい。
”ブートストラップ”というのは、ブーツのかかとの部分についている、小さなループの部分だ。 昔、１８世紀のドイツの貴族、ミュンヒハウゼン男爵が底なし沼に落ちたが、自分の履いていたブーツの小さなブートストラップを引っぱりあげて、自力で脱出したというところから、
"Pulling himself up by his own bootstraps."
は、”自分の力でなんとかする、小さな始まりから大きな成功をつかむ”という意味のことわざになった。（でもこれはただのお話で、本当の話ではないとされている。）
ブートストラップ ロジスティック回帰のスクリプト
以下のWebサイトでスクリプトを公開してくれていた。
感謝感謝！
Bootstrap with logistic regression | R-bloggers
使うデータとロジスティック回帰モデル
使うデータはEZR公式マニュアルのデータである。
まずlibrary(boot)で、boot パッケージを呼び出す。
Stop400が目的変数、Age, Sex, WBC45000の３つの変数を説明変数にしたロジスティック回帰を考える。
fit <- のところをご自分のモデルに置き換えればOK。
library(boot) logit_test <- function (d,indices) { d <- d[indices,] fit <- glm(Stop400 ~ Age + Sex + WBC45000, family=binomial(logit), data=d) return(coef(fit)) }
dは使うデータセット名、indicesは、データの行で、一部を使うなら、その部分を指定することができる。
logit_test(d=Imatinib_Stop, indices=1:100)
ブートストラップする
dataでデータセットを、statisticで上記で設定したロジスティックモデルを指定する。
Rはブートストラップサンプリングの数で、1e4は セットの意味である。
紹介してくれているWebサイトのオリジナルでは1e5で、10万セットになっていた。
boot_fit <- boot( data = Imatinib_Stop, statistic = logit_test, R = 1e4 )
ブートストラップ95％信頼区間を出力する
1万回ブートストラップ ロジスティック回帰分析を行った結果から、偏回帰係数の95％信頼区間を出力してみる。
パラメータとしては、切片と、Age, Sex, WBC45000の3つの合計4つが得られる。
tとあるのがパラメータのブートストラップで得られた値で、切片、Age, Sex, WBC45000の順に1~4行目となっている。
Rは10000回のブートストラップを意味している。
boot.conf.int <- function(m, R){ dat0 <- data.frame() dat1 <- data.frame() for (i in 1:m){ ll <- sort(boot_fit$t[,i])[R*0.025] ul <- sort(boot_fit$t[,i])[R*0.975] name <- paste("t",i,sep="") dat0 <- data.frame(name, ll, ul) dat1 <- rbind(dat1, dat0) } print(dat1) } boot_fit boot.conf.int(m=4, R=10000)
出力結果は以下のとおりである。
boot_fitの結果の、originalというところが、ブートストラップ前の単にロジスティック回帰分析をしたときのパラメータの推定値である。
その95％範囲が、llとulで表現されている。
ゼロをまたがない、t2(Age)とt3(Sex)が、統計学的有意と判断できる。
> boot_fit ORDINARY NONPARAMETRIC BOOTSTRAP Call: boot(data = Imatinib_Stop, statistic = logit_test, R = 10000) Bootstrap Statistics : original bias std. error t1* -1.06315563 -0.068497961 1.05253432 t2* 0.04600895 0.003603028 0.01824306 t3* -1.84683615 -0.148550788 0.77479533 t4* -0.08261246 -0.004526698 0.52341161 > boot.conf.int(m=4, R=10000) name ll ul 1 t1 -3.10373706 0.75911213 2 t2 0.01628948 0.08893698 3 t3 -3.48208222 -0.82632481 4 t4 -1.13685233 0.93418265
ブートストラップをする前の95％信頼区間も同様の結果であったので、顕著な違いはなかった。
X1が点推定値、X2が95％信頼区間下限値、X3が95％信頼区間上限値。
> odds X1 X2 X3 (Intercept) -1.06315563 -2.89982986 0.77351860 Age 0.04600895 0.01453602 0.07748187 Sex[T.M] -1.84683615 -2.94329176 -0.75038054 WBC45000 -0.08261246 -1.02916292 0.86393800
ただし、ブートストラップをすることで、疑似的ではあるものの、多数回の結果をもとにより再現性が高い結果が得られたと言えるかもしれない。
まとめ
小さいサイズのサンプルの時に、復元抽出を多数回行うことで、疑似的な繰り返し試行を行い、信頼性の高い解析結果を提供してくれるブートストラップ手法をロジスティック回帰分析に適用してみた。
サンプルサイズが小さい場合に、この手法を付け加えることで、小サンプルサイズの欠点を補った解析結果が得られると言えるだろう。
ぜひお試しあれ。