はじめに
千葉大学・株式会社Nospareの川久保です．前回の記事で，ブートストラップを用いたバイアス補正法の一般論を解説しました．今回は，具体的な応用例を紹介したいと思います．
EBLUPのMSPE推定
線形混合モデルにおける固定効果と変量効果との混合効果は，経験最良線形不偏予測量（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}
$$
このモデルは，小地域推定で用いられる地域レベルモデルとして，以前紹介しました．混合効果$\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}),
$$
ただし，
\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$を定義しておきます．
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に格納します．
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を列方向に伸ばして格納します．
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つの推定量のパフォーマンスを，バイアスで評価してみます．
> 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 までお問い合わせください．
