はじめに

Bootstrap法による回帰分析を使うかもしれないので、すこし調べてまとめました。参考にしたのは、共立出版の「Rで学ぶデータサイエンスブートストラップ入門」です。

ブートストラップ入門 (Rで学ぶデータサイエンス 4)

作者:汪金芳,桜井裕仁
発売日: 2011/12/21
メディア: 単行本

相変わらずメモレベルです。

Bootstrap法

はじめにBootstrap法についてです。

母集団から標本集団が得られた時に、標本の平均値や分散等の統計量を算出し母集団の統計量 $\theta$ を推定します。

この時の推定は点推定であり、実際には不確実性があります。また、その分布は未知の場合も多くあります。

この不確実性の程度をコンピュータ・シミュレーションによって見積もることができます。

その手法がブートストラップ法です。

具体的には、得られているデータからデータ数 $N$ のサンプルサイズで重複を許すリサンプリングを行います。そして、この得られた $N$ 個のデータを用いて統計量(ex.期待値)を算出します。

この作業を $B$ 回繰り返すことで、統計量の分布を得ることができます。

ブートストラップ法は誤差推定，信頼区間の構成，仮説検定などに用いられます。ここら辺のことについては、以下のブログを参考に...

qiita.com

あと、英語版のWikiのBootstrapping (statistics)を見ると楽しいです。

en.wikipedia.org

Bootstrap法による回帰分析

Bootstrap法よる回帰分析は、得られているデータをリサンプリングし、得られたブートストラップ標本を用いて、回帰係数の期待値の分布推定を行ったり、その有意差の検定を行ったりします。（めちゃくちゃ雑）

残差からのリサンプリングに基づく回帰分析

通常に回帰分析では、回帰分析では誤差項に正規分布を仮定しています。しかし、この仮定が満たされないことも多く、誤差項の分布をブートストラップ法により推定します。誤差 $e_i$ は互いに独立だが、 $i$ ごとに異なる分布に従う場合もあります。

このような場合に適用するために、残差を修正します。修正済み残差 $r_i$ は、平均値が0、分散が $i$ によらず $\sigma ^2$ となります。

$\displaystyle{ r_i=\frac{y_i-\hat \mu_i}{\sqrt{1-1/n-(x_i-\bar x)^2/S_{xx}} } }$

この $r_i$ を用いた場合のアルゴリズムは次にようになります。

モデル $Y_i ^\ast = \hat \beta_0 + \hat \beta_i x_i +\epsilon_i ^\ast (i=1,...,n)$ を考えます。ここで、 $\hat \beta_0$ と $\hat \beta_0$ は最小二乗推定量です。また、 $\bar r=\sum ^n _{i=1} r_i/n$ とします。

中心化した修正済み残差 $r_i-\bar r (i=1,...,n)$ から構成する経験分布関数を $F_n$ とし、 $F_n$ からの無作為復元抽出により、ブートストラップ残差 $\epsilon_1 ^\ast ,..., \epsilon ^\ast_n$ を得る
1に基づき、ブートストラップ標本 $Y_i ^\ast = \hat \beta_0 + \hat \beta_i x_i +\epsilon_i ^\ast (i=1,...,n)$ を計算する
各 $(x_1,y_1)$ を $(x_n,Y_n^\ast)$ で置き換え、ブートストラップ推定値 $( \hat \beta_0 ^\ast , \hat \beta_1 ^\ast )$ 、および $( \hat \sigma ^\ast ) ^2$ を計算する

このアルゴリズムにより得られる $\beta_1$ のブートストラップ最小二乗推定値 $\hat \beta_1^*$ は、次のように表現することができます。

$\displaystyle{ \hat \beta_1^* = \frac{\sum_{i=1}^n(x_i-\bar x)(y_i^*-\bar y^*)}{\sum_{i=1}^n (x_i-\bar x)^2} = \hat \beta_1 + \frac{\sum_{i=1}^n (x_i-\bar x)\epsilon_i^*}{S_{xx}} }$

ここで、 $\bar y ^\ast = \sum ^n _{i=1} y_i ^\ast /n$ です。

データからのサンプリングに基づく回帰分析

次のような一般的な回帰モデルを考えます。

$\displaystyle{ Y_i=\beta_0+\beta_1 x_i+\epsilon_i, i=1,...,n }$

このモデルを考える時、 $(x_i,y_i)$ の関係に興味があります。

そのためリサンプリングは、各 $i$ のデータセット $(x_i,y_i)$ に対して行います。(非常にわかりやすいものですね)

アルゴリズムは次のようになります。

n個の標本数 $(x_1,y_1),...,(x_n,y_n)$ に確率 $1/n$ ずつを付与した経験分布 $F_n$ を構成する
経験分布から、互いに独立なブートストラップ標本 $(x_1,y_1),...,(x_n,y_n)$ を無作為復元抽出する
各 $( x_1, y_1)$ を $( x_1 ^\ast , y_1 ^\ast )$ に置き換え、ブートストラップ推定値 $( \hat \beta_0 ^\ast , \hat \beta_1 ^\ast )$ を計算する

このアルゴリズムでは、誤差項の分散の均一性を仮定していません。

そのため、残差からのリサンプリングに基づく回帰分析に比べ、誤差分散の不均一性に対してロバストであるとされています。

誤差分散の均一性が仮定できそうな場合は、残差からのリサンプリングに基づく回帰分析の方が良いとしています。

Rでデータからのサンプリングに基づく回帰分析

data(gala, package = 'faraway')

gala_new <- gala %>% 
  mutate(Area_log = log(Area))

gala %>% 
  ggplot(aes(x=log(Area), y=Species))+
  geom_point()

f:id:saltcooky:20200224162126p:plain

グラフを見る限りポアソン回帰が良さそうなので、ポアソン回帰分析を行ってみます。

m1 <- glm(Species ~ .-Area, family = poisson(link = "log"), data = gala_new)

推定されたモデルは次のようになりました。

> summary(m1)

Call:
glm(formula = Species ~ . - Area, family = poisson(link = "log"), 
    data = gala_new)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.1377  -2.9277  -0.9179   2.0468   5.2848  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.706e-01  2.932e-01   0.923   0.3560    
Endemics     2.386e-02  1.870e-03  12.756  < 2e-16 ***
Elevation   -6.445e-04  9.520e-05  -6.770 1.29e-11 ***
Nearest      2.629e-03  1.806e-03   1.455   0.1456    
Scruz       -7.677e-04  5.888e-04  -1.304   0.1923    
Adjacent     7.057e-05  3.395e-05   2.078   0.0377 *  
Area_log     4.684e-01  5.056e-02   9.263  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3510.73  on 29  degrees of freedom
Residual deviance:  237.82  on 23  degrees of freedom
AIC: 412.65

Number of Fisher Scoring iterations: 5

それでは、データリサンプリングによるブートストラップ法で回帰係数の分布を推定してみます。

これはcarパッケージを用いると行えます。

2000回リサンプリングして回帰モデルを作成します。

library(car)
beta.boot <- Boot(m1, R = 2000, method = "case")

ここで、Boot関数の引数mathodはcaseだとデータリサンプリングによる回帰、residualsだと残差のリンサンプリングによる回帰になります。

なお、一般化線形モデルの場合、残差のリンサンプリングによる回帰は行えません。

簡単に結果を確認してみます。

> summary(beta.boot)

Number of bootstrap replications R = 2000 
               original    bootBias     bootSE     bootMed
(Intercept)  3.0809e+00 -1.1699e-02 0.22186631  3.0775e+00
Endemics     2.3858e-02  1.5258e-03 0.01065209  2.4385e-02
Elevation   -6.4452e-04 -1.4189e-04 0.00075617 -7.0778e-04
Nearest      2.6289e-03 -2.1081e-03 0.01007110  1.7596e-03
Scruz       -7.6765e-04  3.4389e-05 0.00254443 -8.7876e-04
Adjacent     7.0573e-05  1.9119e-05 0.00054507  7.8835e-05
Area_log     2.0341e-01  9.7690e-03 0.07444646  2.1234e-01

BCa法により算出された信頼区間をみてみます。

> confint(beta.boot, level=.95, type="bca")
Bootstrap bca confidence intervals

                   2.5 %       97.5 %
(Intercept)  2.661620203 3.5194048298
Endemics     0.006895978 0.0455421294
Elevation   -0.001939937 0.0011104991
Nearest     -0.016192585 0.0301755743
Scruz       -0.005363863 0.0048733422
Adjacent    -0.001132596 0.0007036367
Area_log     0.041524018 0.3335293312

推定された回帰係数の分布を確認してみます。

ggboot_histogram <- function(BD, X, bw = 1/50){
  tmp <- BD$t %>% as.data.frame() 
  
  wide <- max(tmp %>% pull(X))-min(tmp %>% pull(X))

  ggplot(tmp,aes_(x = as.name(X)))+
    geom_histogram(binwidth = wide*bw)+
    geom_density(eval(bquote(aes(y = ..count..*.(wide*bw)))), fill = 1, alpha = 0.3)+
    geom_vline(xintercept=BD$t0[X], linetype="dashed", colour = "red", size=.2) +
    ggtitle(paste0(X, " coeff by ", BN, " Bootstrap regression \n",
                   "Mean:", round(BD$t0[X],4), "   2.5%:", round(confint(beta.boot)[X,1],4), " - 97.5%:",round(confint(beta.boot)[X,2],4)))
}

ggboot_histogram(BD=beta.boot, "Area_log")

f:id:saltcooky:20200224163047p:plain

他の変数における回帰係数の分布も見てみます。

f:id:saltcooky:20200224163135p:plain

ElevationとNearestは、正規分布に従ってはなさそうです。

ブートストラップ法を用いることで、適切な分布が得られていると言えそうです。

Rで相関モデルの回帰係数の有意差検定

ブートストラップ信頼区間を用いた検定

上記で生成したBCa法により算出された信頼区間を用いて、その区間に0が含まれている場合には、統計的に有意であるといえます。実際に上記で生成したブートストラップ標本から片側5%信頼区間を見てみます。

> confint(beta.boot, level=.9, type="bca")
Bootstrap bca confidence intervals

                      5 %         95 %
(Intercept)  2.7500094380 3.4458113539
Endemics     0.0108198378 0.0410965876
Elevation   -0.0017787189 0.0005420093
Nearest     -0.0128395052 0.0187406325
Scruz       -0.0042417713 0.0042695221
Adjacent    -0.0007962402 0.0006028746
Area_log     0.0731049579 0.3059796320