統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

ブートストラップ ロジスティック回帰分析をEZRで行う方法

kaiseki daiko banner

ブートストラップとは、小さいサンプルを復元抽出で何度もサンプリングして、疑似的に何度も試行したことにして、その結果から推定値の直接的な95%の範囲を求めたりする方法。

ブートストラップを使ったロジスティック回帰分析の方法をご紹介。

ブートストラップとは?

ブートストラップとは、小さいサンプルを復元抽出で何度もサンプリングして、疑似的に何度も試行したことにして、その結果から推定値の直接的な95%の範囲を求めたりする方法。

復元抽出とは、たまたま複数回同じデータがとられてもOKとする方法。

オリジナルのサンプルとサイズ(データの数)が同じサンプルをたくさん作る。

このほうが現実世界に近い。

こうすると、データの構成がばらつく。

オリジナルのデータから、似たようなサンプルをたくさん作って、小さいサイズのオリジナルデータでも信頼性高い、区間推定を行える方法だ。

ちなみに、ブートストラップは靴紐(ブーツのひもかな?)という意味で、アメリカのことわざで、湖に散歩に行って溺れた人が、自分の靴紐につかまって生還するというものがあり、そこから、困ったことがあっても自力で何とかするという意味だそうである。

そこから、小さいサイズのサンプルでも、そのサンプルだけで、自力で何とかするという意味で、ブートストラップ法という名前がついたのだとか。

新谷先生が、発明者のEfron先生から直々に聞いたとのことで、間違いないようだ。

24秒あたりから説明あり。

youtu.be

ググってみたところ、こんなブログ記事があった。

ameblo.jp

ブートストラップとは、ブーツのかかとについている、小さなループの部分のことを言っているみたい。

”ブートストラップ”というのは、ブーツのかかとの部分についている、小さなループの部分だ。 昔、18世紀のドイツの貴族、ミュンヒハウゼン男爵が底なし沼に落ちたが、自分の履いていたブーツの小さなブートストラップを引っぱりあげて、自力で脱出したというところから、

"Pulling himself up by his own bootstraps." 

は、”自分の力でなんとかする、小さな始まりから大きな成功をつかむ”という意味のことわざになった。(でもこれはただのお話で、本当の話ではないとされている。)

ブートストラップ ロジスティック回帰のスクリプト

以下のWebサイトでスクリプトを公開してくれていた。

感謝感謝!

Bootstrap with logistic regression | R-bloggers

使うデータとロジスティック回帰モデル

使うデータはEZR公式マニュアルのデータである。

まずlibrary(boot)で、boot パッケージを呼び出す。

Stop400が目的変数、Age, Sex, WBC45000の3つの変数を説明変数にしたロジスティック回帰を考える。

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は  1 \times 10^4 = 10000 セットの意味である。

紹介してくれている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

ただし、ブートストラップをすることで、疑似的ではあるものの、多数回の結果をもとにより再現性が高い結果が得られたと言えるかもしれない。

まとめ

小さいサイズのサンプルの時に、復元抽出を多数回行うことで、疑似的な繰り返し試行を行い、信頼性の高い解析結果を提供してくれるブートストラップ手法をロジスティック回帰分析に適用してみた。

サンプルサイズが小さい場合に、この手法を付け加えることで、小サンプルサイズの欠点を補った解析結果が得られると言えるだろう。

ぜひお試しあれ。

EZR公式マニュアル