B ブートストラップ標準誤差

Efron (1979, Bootstrap Methods: Another Look at the Jackknife)は再標本化によって標準誤差や分布の分位点を推定する方法をいくつも提案し, これらをブートストラップ法と総称した. ブートストラップ法はパラメトリックなものとノンパラメトリックなものに大別されるが, 今回はノンパラメトリックなブートストラップ法を利用する. 最も単純なブートストラップ法は次の通り⁴⁵. $N$ 件のデータ $\mathcal{D}=\{z_{1},\cdots,z_{N}\}$ を取得できたとして, ここから復元抽出した $N$ 件の標本(ブートストラップ標本)を $B$ 個作成する. $\begin{aligned} D_{1}= & \{z_{1}^{(1)},\cdots,z_{N}^{(1)}\},\\ \vdots\\ D_{B}= & \{z_{1}^{(B)},\cdots,z_{N}^{(B)}\}\end{aligned}$ もし, この擬似的に作られたブートストラップ標本 $D_{b}$ が母集団から無作為抽出された標本と同等の分布に従うならば, それぞれの標本平均 $\bar{z}^{(b)}$ はすべて同一の分布にしたがう. よって, 以下のように計算した標本標準偏差を標準誤差の推定値とすることができる.

$\begin{aligned} \hat{\mathit{SE}}(\bar{z}):= & \sqrt{\frac{1}{B}\sum_{b=1}^{B}(\bar{z}^{(b)}-\bar{z}_{B})^{2}},\\ \bar{z}_{B}:= & \frac{1}{B}\sum_{b=1}^{B}\bar{z}^{(b)}\end{aligned}$

これが標本平均でなく回帰分析の係数に置き換わっても同じことが言える. つまり今回のOP法のように推定量の分布を解析的に導出するのが難しい場合にブートストラップ法を利用する場面がある.

しかし, 「ブートストラップ標本が母集団から無作為抽出された標本と同等の分布に従う」という前提を忘れてはならない. 上記の単純なケースではどの標本も等確率に抽出されており, 言い換えるなら母集団の分布がそのようになっていると暗黙に仮定していることになる. Efron (1979, Bootstrap Methods: Another Look at the Jackknife)は, 推定した回帰モデルの残差から誤差項の経験分布を構成し, これを誤差項 $\varepsilon$ の分布とみなしてサンプリングした乱数 $\hat{\varepsilon}_{i}^{\ast}$ とモデルの予測値の和 $\hat{y}_{t}^{\ast}:=\hat{y}_{t}+\hat{\varepsilon}_{t}^{\ast}$ を生成することを提案している.

標本からの計算のため, 大数の法則や中心極限定理からの類推で直観的に $B$ が大きいほど推定の精度がよくなると予想できる. しかしながら, 具体的にどれくらいなら十分であるかを示す理論は存在しないようだ. ネット上では具体的な数字を挙げる例も見られるが, 具体的な根拠は挙げられていないことが多い⁴⁶. これに関しては Rule of thumb for number of bootstrap samples - Cross Validated のAdamOの回答が参考になる. いわく, 十分な回数は問題に応じて変わる(例えばコーシー分布に従うデータならば何万, 何千万回やろうとそも存在しない平均を求めるのに十分ではない)ため, ブートストラップ計算の結果をヒストグラムに表し分布が適切な形になっているかを確認せよ, ということである. とはいえ, 彼の回答の追加コメントでも指摘されるように, 多くの場合は回数が多いほど精度が増すはずなので, 例えば先行研究と比較できるようにいったん同じ回数で計算した上でヒストグラムを確認し, もしうまくいってないなら回数を増やす (またはバグを疑う), といった方法が現実的ではないかと考える.

ブートストラップ法の計算に関するRのパッケージはsimplebootとbootパッケージの2つが代表的である⁴⁷.simplebootは名前の通り簡単なブートストラップ法の計算だけに限定したパッケージで, OP法のような特殊なモデルには対応していないからboot::boot()を使う. ただし, npパッケージによるカーネルベースの推定方法はただでさえ時間がかかるため, 今回は1段階目を多項式近似でやることにした.

boot()の使い方はパッケージのヘルプを見ればすぐ分かるが, 一応簡単に解説しておく. 第1引数dataにデータを, 第2引数statisticに計算したい統計量(のベクトル)を返す関数を与える. boot()はdataをリサンプリングしてブートストラップ標本を作り, それぞれをstatisticに与えてブートストラップ推定値を計算する. statisticに与える関数の第1引数にはdataが, 第2引数には生成ブートストラップ標本がdataのどこかを表すインデックスベクトルが与えられる. さらに関数ヘルプには補外データに対してブートストラップ法を適用する例が掲載されている. この場合, 予測値の計算にももう1つインデックスベクトルが必要になる. これが第3引数である. 第4引数以降は, ユーザーが自由に定義できる引数である. そのため, 既に書いたように回帰モデルに対してノンパラメトリックブートストラップ法を適用する場合は, 予め残差を計算して入力データに含めておく必要がある.

今回はのように書いた. 回帰モデルのブートストラップ統計量の計算のため, dataには通常の推定に必要な変数に加え, ブートストラップ標本を生成するための $\hat{y}_{i,t},\hat{\varepsilon}_{i,t}$ の列を追加する必要がある.

今回statisticが出力すべき統計量は3つの係数パラメータ $\beta_{A},\beta_{K},\beta_{L}$ だが, 確認のため各段階で推定したモデルの当てはまりについても出力するようにしている⁴⁸.

op_params_poly <- list(beta_k = exp(fit_2nd_poly$par[2]), beta_l = unname(coef(fit_1st_poly)["l"]))

estimate_op_ordinary <- function(data, degree_1, degree_p, degree_2nd){
  fit_1st <- lm(y ~ l + poly(k, inv, degree = degree_1), data = data, na.action = na.omit)
  beta_l <- unname(coef(fit_1st)["l"])
  rmse_1st <- sqrt(mean(resid(fit_1st)^2))
  fit_exit <- glm(x ~ l + poly(k_lag, inv_lag, degree = degree_p), data = data %>% drop_na(k_lag, inv_lag), family = binomial(link = "probit"))
  std_Brier <- std_Brier(fit_exit$data$x, fitted(fit_exit))
  df_2nd <- make_2nd_sample(
    data, beta_l,
    phi = predict(fit_1st, newdata = data) - beta_l * data$l,
    p_x = predict(fit_exit, newdata = data, type = "response")
  )
  fit_2nd <- optim(
    par = c(rnorm(n = 1, sd = .5), log(ifelse(1 - beta_l <=0, .5, 1 -beta_l)) + rnorm(n = 1, sd = .3), rnorm(n = sum(1:degree_2nd + 1), sd = .5)),
    fn = obj_op_2nd, data = df_2nd, degree = degree_2nd)
  beta_a <- mean(with(df_2nd, y_bl - exp(fit_2nd$par[2]) * k), na.rm = T)
  return(c(beta_a = beta_a, beta_l = beta_l, beta_k = unname(exp(fit_2nd$par[2])), rmse_1st = rmse_1st, std_Brier = std_Brier, rmse_2nd = sqrt(fit_2nd$value)))
}

get_stat_boot <- function(data, indices, indices_, estimate, params){
  # estimate = estimate function
  # params = model hyper parameter (e.g., degrees of polynomials)
  data$y <- with(data, fitted + resid[indices])
  do.call(estimate, args = c(list(data = data), params)) }

boot_ordinary <- df_sample_2nd_poly %>%
  mutate(., fitted = if_else(is.na(y), NA_real_, do.call(get_fitted, args = c(list(data = .), op_params_poly))),
         resid = y - fitted) %>%
  boot(., get_stat_boot, R = 100, strata = .$x, estimate = estimate_op_ordinary, params = list(degree_1 = 3, degree_p = 2, degree_2nd = 2))

なお, 既に書いたようにOP法では $\beta_{A}$ を識別できないが, 他の方法との比較のため $\hat{\beta}_{A}:=(NT)^{-1}\sum_{i,t}(y_{i,t}-\hat{\beta}_{L}l_{i,t}-\hat{\beta}_{K}k_{i,t})$ として計算した. これは本来の仮定に加えて $\mathrm{E}\omega_{t}=0$ が成り立ってなければ適切な推定量とは言えないが, 今回の乱数データでは $\omega_{t}$ が平均ゼロの定常分布に収束するような方法で生成している.

引数Rはブートストラップ標本の生成数である. 本文での解説に沿って100に設定し, 各パラメータのヒストグラムが不自然でないか確認した. strataは層別抽出をするためのインデックスである. 今回は $y_{t}$ に欠損があり, OP法は欠損していないデータだけで推定しているから, ブートストラップ標本 $\hat{y}_{t}^{\ast}$ も本来のデータの欠損メカニズムを再現する必要がある. ここではstrataに元のデータの欠損判定 $x_{t}$ を与えることでこれを表現している⁴⁹. また, boot.ci()でboot()の結果から信頼区間を計算できるが, 今回は利用しない.

例えば標準誤差の計算なら100程度, 分位点の計算なら1000程度, という記述がよく見られる. この数値の根拠は Efron (1979) の論文に掲載されている実験でそのように設定されていたからではないかと思う. しかしEfronはこの数字の根拠を述べていない.↩
金明哲『Rとブートストラップ』と奥村晴彦『ブートストラップ』はいずれもこれらのRパッケージを使った計算方法も含めて言及している. 『ブートストラップのためのbootパッケージ』にもbootパッケージの言及がある.↩
サンプルプログラムでは返り値に名前を付けているが, boot()は名前の情報を保存してくれないようだ↩
欠損値を含む場合のブートストラップ法に詳しくないが, $y_{t}$ に対して $x_{t}$ は先決(predetermined)変数だから, 単純に $x_{t}$ の層化で欠損を条件付けるだけで十分だと思う.↩