ブートストラップを用いたバイアス補正（2）

# はじめに
千葉大学・株式会社Nospareの[川久保](https://sites.google.com/site/ykawakubostat/)です．[前回の記事](https://qiita.com/ykawakubo/items/88e661b99f4465ae3445)で，ブートストラップを用いたバイアス補正法の一般論を解説しました．今回は，具体的な応用例を紹介したいと思います．


# EBLUPのMSPE推定
[線形混合モデル](https://qiita.com/ykawakubo/items/f57c38b72d2ddc7798fd)における固定効果と変量効果との混合効果は，経験最良線形不偏予測量（empirical best linear unbiased predictor, EBLUP）で予測されます．EBLUPの予測リスクは，平均二乗予測誤差（mean squared prediction error, MSPE）で評価されることが多いのですが，MSPEを推定する問題を考えたいと思います．


### モデル
例として，以下の線形混合モデルを考えます．
$$
y_i = x_i^\top \beta + b_i + \varepsilon_i, \quad b_i \overset{\mathrm{iid}}{\sim} \mathrm{N}(0,\tau^2), \quad \varepsilon_i \overset{\mathrm{indep}}{\sim} \mathrm{N}(0,D_i), \quad i = 1,\dots,n. \tag{1}
$$
このモデルは，[小地域推定](https://qiita.com/ykawakubo/items/33bd68dfa6029627e125)で用いられる地域レベルモデルとして，以前紹介しました．混合効果$\xi_i$は，平均所得などの地域ごとに関心のある量で，補助変数$x_i$と地域効果$b_i$でモデル化されます．実際のデータはサンプルから計算された$y_i$で，真の$\xi_i$に誤差$\varepsilon_i$をともない観測されます．

このモデルにおける誤差分散$D_i$は既知と仮定し，未知パラメータは$\theta = (\beta^\top, \tau^2)^\top$です．ここで，混合効果$\xi_i = x_i^\top \beta + b_i$を予測対象としたとき，そのEBLUPは，
$$
\hat{\xi}_i = {\hat{\tau}^2 \over \hat{\tau}^2 + D_i}y_i + {D_i \over \hat{\tau}^2 + D_i}x_i^\top \hat{\beta}, \tag{2}
$$
という形で与えられます．


### EBLUPのMSPE
混合効果$\xi_i$を，そのEBLUP$\hat{\xi}_i$で予測しますが，その予測誤差の大きさがどの程度であるかを見積もることが重要です．そこで，予測リスクを以下のMSPEで評価します．
$$
\mathrm{MSPE}_i = E[(\hat{\xi}_i - \xi_i)^2].
$$
ここで(2)式のEBLUPにおける$\hat{\tau}^2$を，$\tau^2$の最尤推定量としたとき，MPSEは以下のように近似できることが知られています．

$$
\mathrm{MSPE}\_i = g_{1i}(\tau^2) + g_{2i}(\tau^2) + g_{3i}(\tau^2) + o(n^{-1}),
$$
ただし，

```math
\begin{align}
g_{1i}(\tau^2) &= \frac{\tau^2 D_i}{\tau^2 + D_i} = O(1), \\
g_{2i}(\tau^2) &= \left( \frac{D_i}{\tau^2 + D_i} \right)^2 x_i^\top \left\{ \sum_{j = 1}^n (\tau^2 + D_j)^{-1} x_j x_j^\top \right\}^{-1} x_i = O(n^{-1}), \\
g_{3i}(\tau^2) &= \frac{2D_i^2}{(\tau^2 + D_i)^3} n^{-2} \sum_{j=1}^n (\tau^2 + D_j)^2 = O(n^{-1}).
\end{align}
```

すなわち，MSPE（の近似）はパラメータ$\tau^2$の関数として書けます．これを，$\tau^2$をプラグインした推定量で推定しようとすると，$g_{2i}(\hat{\tau}^2)$と$g_{3i}(\hat{\tau}^2)$についてはバイアスが$o(n^{-1})$のオーダーで十分に小さいことが知られていますが（二次漸近不偏），$g_{1i}(\hat{\tau}^2)$のバイアス$E[g_{1i}(\hat{\tau}^2)] - g_{1i}(\tau^2)$については$O(n^{-1})$のオーダーであることが知られています．このバイアスを解析的に求めることも実はできるのですが，今回はパラメトリックブートストラップ法により，以下のようにバイアス補正してみます．
$$
2g_{1i}(\hat{\tau}^2) - E_{\hat{\tau}^2}[g(\hat{\tau}^{2\ast})]. \tag{3}
$$
ここで$E_{\hat{\tau}^2}[\cdot]$は，$y^\ast = (y_1^\ast,\dots,y_n^\ast)^\top$とおいたとき，
$$
y_i^\ast = x_i^\top \hat{\beta} + b^\ast_i + \varepsilon^\ast_i, \quad b_i \overset{\mathrm{iid}}{\sim} \mathrm{N}(0,\hat{\tau}^2), \quad \varepsilon^\ast_i \overset{\mathrm{indep}}{\sim} \mathrm{N}(0,D_i)
$$
にしたがう確率分布についての期待値， $\tau^{2\ast}$は$y^\ast$から推定した変量効果の分散の最尤推定量を表します．(3)式の第2項は，十分に大きな$B$個のブートストラップサンプル$y^\ast(1),\dots,y^\ast(B)$を用いて，
$$
{1 \over B}\sum_{b=1}^B g_{1i}(\hat{\tau}^{2\ast}(b))
$$
とモンテカルロ近似します．


### Rによる数値実験
まず，関数$g_{1i}(\cdot), g_{2i}(\cdot), g_{3i}(\cdot)$と，データ生成過程を指定するための(1)式のパラメータ値と補助変数$x_i$を定義しておきます．
```R
library(sae)
set.seed(12345)

g1 <- function(t2){
  t2 * Di / (t2 + Di)
}

g2 <- function(t2){
  X <- cbind(1, x)
  
  I1 <- (Di / (t2 + Di))^2
  I2 <- diag( X %*% solve( crossprod(X, X / (t2 + Di)) ) %*% t(X) )
  
  return(I1 * I2)
}

g3 <- function(t2){
  I1 <- 2*Di^2 / (t2 + Di)^3
  I2 <- sum( (t2 + Di)^2 ) / n^2
  
  return(I1 * I2)
}

R1   <- 10000  # MSPEの真値をシミュレーションするためのreplication数
R2   <- 100    # 数値実験のreplication数
B    <- 1000   # ブートストラップサンプルサイズ
n    <- 15     # モデルのサンプルサイズ
tau2 <- 1      # 変量効果の分散
Di   <- as.vector( sapply( seq(0.7, 0.3, by = -0.1), rep, 3 ) )
               # 誤差分散
x    <- rep(c(1,0,0),5)  # 説明変数
be   <- c(1, 2)          # 回帰係数
```

`sae`というパッケージをロードしていますが，(1)式のモデルを推定するための関数が実装されています．説明変数$x_i$と，誤差分散$D_i$は数値実験のreplicationにおいて値を固定するので，関数`g1`，`g2`，`g3`の中の`Di`と`x`というオブジェクトは，グローバル環境からサーチして利用します．真のMSPEの値は分からないため，`R1`$(=10000)$個のreplicationを用いてシミュレーションし，値を`MPPE`に格納します．
```R
MSPE <- numeric(n)
for(r in 1:R1){
  bi <- rnorm(n, mean = 0, sd = sqrt(tau2))  # 変量効果
  ep <- rnorm(n, mean = 0, sd = sqrt(Di))    # 誤差項
  xi <- be[1] + x*be[2] + bi                 # 混合効果
  y <- xi + ep                               # 観測y
  
  resEBLUP <- eblupFH(y ~ x, Di, method = "ML")      # (1)式モデルの当てはめ
  MSPE <- MSPE + (xi - as.vector(resEBLUP$eblup))^2  # MSPEのシミュレーション
}
MSPE <- MSPE / R1
```

`R2`$(=1000)$個のreplicationから，MSPEの推定値を求めます．プラグイン推定値$g_{1i}(\hat{\tau}^2) + g_{2i}(\hat{\tau}^2) + g_{3i}(\hat{\tau}^2)$と，ブートストラップでバイアス補正した(3)式の推定値をそれぞれ，`mspe1`と`mspe2`に格納します．サンプルサイズが`n`(=15)なので，MSPEの推定値も15個存在し，さらにreplication数1000を列方向に伸ばして格納します．

```R
mspe1 <- mspe2 <- mspe3 <- matrix(nrow = n, ncol = R2)
for(r in 1:R2){
  bi <- rnorm(n, mean = 0, sd = sqrt(tau2))
  ep <- rnorm(n, mean = 0, sd = sqrt(Di))
  xi <- be[1] + x*be[2] + bi
  y <- xi + ep
  
  resEBLUP <- eblupFH(y ~ x, Di, method = "ML")
  tau2hat <- resEBLUP$fit$refvar  # 変量効果の分散の推定値
  g1_plug <- g1(tau2hat)
  g2_plug <- g2(tau2hat)
  g3_plug <- g3(tau2hat)
  mspe1[,r] <- g1_plug + g2_plug + g3_plug  # プラグイン推定値
  
  g1boot <- numeric(n)
  betahat <- resEBLUP$fit$estcoef$beta
  B_count <- 0

  # ブートストラップreplication
  for(b in 1:B){
    bi_boot <- rnorm(n, mean = 0, sd = sqrt(tau2hat))
    ep_boot <- rnorm(n, mean = 0, sd = sqrt(Di))
    y_boot <- betahat[1] + x*betahat[2] + bi_boot + ep_boot  # ブートストラップサンプルの発生
    
    resEBLUP <- eblupFH(y_boot ~ x, Di, method = "ML")
    tau2hat_boot <- resEBLUP$fit$refvar
    if( !is.na(tau2hat_boot) ){
      B_count <- B_count + 1
      g1boot <- g1boot + g1(tau2hat_boot)
    }
  }
  g1boot <- g1boot / B_count

  mspe2[,r] <- 2*g1_plug - g1boot + g2_plug + g3_plug  # バイアス補正推定値
```
ブートストラップreplicationのfor文中に，if文を書いていますが，ブートストラップサンプルのうちごく稀に，極端な値が発生されるケースがあり，この場合(1)式モデルの当てはめが上手くいかず，`eblupFH`が`NA`を返してしまいます．このようなブートストラップサンプルがもし出てきた場合は，それを捨てて計算するためのif文です．2つの推定量のパフォーマンスを，バイアスで評価してみます．
```R
> apply(mspe1 - MSPE, 1, mean)
 [1] -0.06812528 -0.08163627 -0.07485778 -0.05707919 -0.06407306
 [6] -0.06324323 -0.04059297 -0.05388494 -0.05379896 -0.04066726
[11] -0.04392839 -0.03911479 -0.01601116 -0.02073228 -0.02658556
> apply(mspe2 - MSPE, 1, mean)
 [1]  3.594696e-04 -1.315152e-02 -6.373026e-03  4.530277e-03
 [5] -2.463592e-03 -1.633766e-03  1.324644e-02 -4.554396e-05
 [9]  4.044568e-05  4.357612e-03  1.096486e-03  5.910087e-03
[13]  1.900246e-02  1.428134e-02  8.428058e-03
```
プラグイン推定量`mspe1`は，真のMSPEをunder-estimateしているのに対し，バイアス補正推定量`mspe2`は，ほぼ不偏に推定できていると言って良い結果になりました．


# まとめ
株式会社Nospareには，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては[株式会社Nospare](https://bit.ly/38WtGZH) までお問い合わせください．
<a href="https://bit.ly/38WtGZH">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1023506/d5c625c0-d6f4-680b-d702-3c40f11d6c39.jpeg" alt="株式会社Nospare" width=400>
</a>