« 読了: Wall (2004) SARモデル・CARモデルは実は君が思っているような空間モデルではないかもしれないよ？ | メイン | 読了：久保川(2006) 線形混合モデルで小地域推定 »
2018年7月20日 (金)
Chandra, H., Kumar, S., Aditya, K. (2016) Small area estimation of proportions with different levels of auziliary data. Biometrical Journal, 60(2), 395-415.
二値データの小地域推定について、推定方法をシミュレーションで比較した実験。仕事の都合で、まさにこういうのを探していたんです。著者らはインド・ニューデリーの人。ありがとう、インドの先生。
ちょっと細かめにメモを取った。
1. イントロダクション
二値データの小地域推定、つまり小地域における割合の推定においては、ロジット・リンクの一般化線形混合モデル(GLMM)が用いられることが多い。ユニットレベル共変量がある場合、小地域における割合の推定には経験的プラグイン予測量(EPP)が使われることが多い。実はMSEを最小化する推定量は他にあるのだが(経験的最良予測量, EBP)、閉形式でなく数値近似が必要なのであまり使われていない。
一般に小地域モデルは地域レベルとユニットレベルに大別される。地域レベルではFay-HerriotモデルによるEBLUP推定量が知られている。地域レベルモデルもユニットレベルモデルも線形混合モデルの特殊ケースである。
二値の場合、地域レベルでもユニットレベルでも、結局はロジット・リンクのGLMMが使われているわけだが、ユニットレベルモデルのEPPと地域レベルモデルのEPPは異なる。最近では地域レベルポアソン混合モデルのEBPというのも開発されている。
二値データの場合、どういう共変量があったらどういう推定量を使うべきか。本論文では、ロジット・リンクのGLMM(すなわちロジスティック正規混合モデル)のもとでのいろんな推定量について述べる。
2. 小地域推定
母集団$U$(サイズ$N$)が、小地域$U_1, \ldots, U_D$(サイズ$N_1, \ldots, N_D$)に重なりなく完全に分割されているとする。母集団から標本$s$(サイズ$n$)を取ってきた。小地域別にみると$s_1, \ldots, s_D$(サイズ$s_1, \ldots, s_D$)である。抽出されなかった残りの部分を$r_1, \ldots, r_D$とする。
小地域$i$のユニット$j$の関心ある変数の値を$y_{ij}$とする。本論文ではこれを二値とする。
いま推測したいのは$P_i = (1/N_i) \sum_{j \in U_i} y_{ij}$である($U_i$は$s_i$と$r_i$からなるという点に注意)。
デザインベース直接推定量(DIR)は
$\hat{P}^{DIR}_i = \sum_{j \in s_i} w_{ij} y_{ij}$
と書ける。ここで$w_{ij}$は調査ウェイトで、$\sum_{j \in s_i} w_{ij} = 1$と基準化されている。
$\hat{P}^{DIR}$のデザインベース分散は
$var(\hat{P}^{DIR}_i) \approx \sum_{j \in s_j} w_{ij}(w_{ij}-1)(y_{ij} - \hat{P}^{DIR}_i)^2$
と近似できる。
SRSならば、
$\hat{P}^{DIR}_i = (1/n_i) \sum_{j \in s_i} y_{ij}$
$var(\hat{P}^{DIR}_i) \approx (1/n_i) \hat{P}^{DIR}_i (1-\hat{P}^{DIR}_i)$
である。
では小地域モデルについて。
ユニットレベル共変量を$\mathbf{x}_{ij}$(長さ$p$)とする[書いてないけどどうやら切片項もコミである]。$\pi_{ij} = Prob(y_{ij} = 1)$として、ロジットリンクGLMM
$logit(\pi_{ij}) = \mathbf{x}^T_{ij} \mathbf{\beta} + u_i, \ \ u_i \sim N(0, \sigma^2)$
を考える。$u_i$は互いに独立。[←ここで$j$が$N_i$まで動く点に注意。つまりこれは母集団モデルである]
このモデルによる経験的プラグイン予測量(EPP)は
$\hat{P}^{EPP}_i = (1/N_i) (\sum_{j \in s_i} y_{ij} + \sum_{j \in r_i} \hat{\mu}_{ij})$
ただし$\hat{\mu}_{ij} = logit^{-1} (\mathbf{x}^T_{ij} \hat{\beta} + \hat{u}_i)$である。
[原文では、$\hat{\mu}_{ij} = \exp(\mathbf{x}^T_{ij} \hat{\beta} + \hat{u}_i) / (1+\exp(\mathbf{x}^T_{ij} \hat{\beta} + \hat{u}_i))$とあるけれど、書くのが面倒なので$logit^{-1}(\cdot)$と略記する。以下同様]
ここで$\hat{\beta}$は固定効果パラメータの推定値、$\hat{u}_i$はランダム効果パラメータの予測値である。
GLMMの推定手続きとしてもっともよく使われているのはPQLである。この方法では、反応変数の非正規な分布を線形近似し、線形化された従属変数が近似的に正規分布に従うと仮定する。場合によっては一致性・不偏性が失われるものの、経験的には上手くいく。
ここからは、共変量の入手可能性に応じた、３つのケースについて考える。
ケース1、標本ユニットでのみ共変量が入手可能であるケース。
この場合、モデルのあてはめは可能だが、EPPは手に入らない。そこで次の推定量(EPP1)を提案する。
$\hat{P}^{EPP1}_i = \sum_{j \in s_i} w_{ij} \hat{\mu}_{ij}$
EPP1はDIRより効率的だが、モデルの誤指定によるバイアスはDIRより大きい。EPP1は、標本抽出ウェイトに対応する回数だけ共変量の値を反復し、共変量のセンサスを再作成し、このセンサスを使って、ロジット混合モデルで確率を予測していることになる[←ややこしい...]
EPP1の欠点は、標本ユニットのアウトカムを使っていない点である。
[↑そりゃそうだよなあ、これは変だ。もっと良い推定量があるんじゃなかろうか...]
ケース2, 標本ユニットで共変量が入手可能で、かつそれらの共変量の母集団での集計値が入手可能であるケース。
$\bar{\mathbf{X}}_i = (1/N_i) \sum_{j \in U_i} \mathbf{x}_{ij}$, $\bar{\mathbf{x}}_i = (1/N_i) \sum_{j \in s_i} \mathbf{x}_{ij}$とすると、ここから非標本ユニットにおける平均$\bar{\mathbf{X}}_{r_i}$がわかる[面倒なので式は省略するが、そりゃわかるよな]。
まずモデルをあてはめて$\hat{\beta}$, $\hat{u}_i$を得る。で、$\hat{\mu}^{EPP2}_{ij} = logit^{-1} (\bar{\mathbf{X}}^T_{r_i} \hat{\beta}+ \hat{u}_i)$を求める。推定量は
$\hat{P}^{EPP2}_i = (1/N_i) (\sum_{j \in s_i} y_{ij} + \sum_{j \in r_i} \hat{\mu}^{EPP2}_{ij} )$
標本平均$p_i = (1/n_i) \sum_{j \in s_i} y_{ij}$, 標本抽出割合$f_i = n_i / N_i$を使って書き直すと
$\hat{P}^{EPP2}_i = f_i p_i + (1-f_i) \hat{\mu}^{EPP2}_{ij}$
$f_i$が小さく、かつ$\bar{\mathbf{X}}_i \approx \bar{\mathbf{X}}_{r_i}$なら、$\bar{\mathbf{X}}_{r_i}$の代わりに$\bar{\mathbf{X}}_i$を使って
$\hat{P}^{EPP2.1} = f_i p_i + logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\beta} + \hat{\mathbf{u}}_i)$
でもよい。$f_i$が無視できる大きさなら、
$\hat{P}^{EPP2.2} = logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\beta} + \hat{\mathbf{u}}_i)$
でもよい。
ケース3, 共変量の母集団集計値だけが入手可能であるケース。
この場合はさきほどのモデルのあてはめ自体が無理で、
$logit(\pi_{ij}) = \bar{\mathbf{X}}^T_i \mathbf{\alpha} + v_i, \ \ v_i \sim N(0, \sigma_v^2)$
をあてはめるしかない。で、$\hat{\mu}^{EPP3}_{ij} = logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\alpha}+ \hat{u}_i) $ を求める。推定量は
$\hat{P}^{EPP3}_i = (1/N_i) (\sum_{j \in s_i} y_{ij} + \sum_{j \in r_i} \hat{\mu}^{EPP3}_{ij} ) = f_i p_i + (1-f_i) \hat{\mu}^{EPP3}_{ij}$
となる。標本抽出割合が無視できる大きさなら、
$\hat{P}^{EPP3.1}_i = \hat{\mu}^{EPP3}_{ij}$
でもよい。これはEPP2.2と同じになるけど、元のモデルは異なる。[ああそうか。どちらも横着して非標本ユニットだけに注目している点では同じだが、EPP2.2は非標本ユニットにおける共変量の母平均を考えるのが面倒だから全体の母平均を使っているのに対して、EPP3.1は最初から全体の母平均しかわかんないわけだ]
...という路線のほかに、ふたつ代替路線がある。両方とも、ケース2と3で使える。
代替路線その１，地域レベルでモデルを組む。
$logit(\pi_i) = \bar{\mathbf{X}}^T_i \lambda + a_i, \ \ a_i \sim N(0, \sigma_a^2)$
とする。推定量は
$\hat{P}^{EPP.A}_i = logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\lambda}+ \hat{a}_i)$
となる。まず$\sigma_a^2$をREML推定し、次に$\lambda$と$a_1, \ldots, a_D$を
PQL推定する。詳細はChandra et al. (2011 J.App.Stat.), Johnson et al. (2010 J.Official Stat.)をみよ。
代替路線その２，Fay-Herriotモデルを組んじゃってEBLUP推定量を使う(EBLUP.FH)。これは推定値が$[0,1]$の外側に出ちゃうかもしれない。[Fay-Herriotモデルの場合、小地域ごとの標本抽出誤差の分散が既知でないといけないと思うんだけど、そこはどうやったんだろう？ $n_i p_i (1-p_1)$としたのかな？]
[話が長くなってしまったので、EBLUP.FHを除くすべての推定量を同一の記法で書き直しておく。なお、
$logit(\pi_{ij}) = \mathbf{x}^T_{ij} \mathbf{\beta} + u_i, \ \ u_i \sim N(0, \sigma^2)$
$logit(\pi_{ij}) = \bar{\mathbf{X}}^T_i \mathbf{\alpha} + v_i, \ \ v_i \sim N(0, \sigma_v^2)$
である。うーん、こうしてみると、標本ウェイトがEPP1でしか使われていないという点も不思議だなあ...]
- $\hat{P}^{EPP}_i = f_i p_i + (1-f_i) logit^{-1} (\mathbf{x}^T_{ij} \hat{\beta} + \hat{u}_i)$
- $\hat{P}^{EPP1}_i = \sum_{j \in s_i} w_{ij} logit^{-1} (\mathbf{x}^T_{ij} \hat{\beta} + \hat{u}_i)$
- $\hat{P}^{EPP2}_i = f_i p_i + (1-f_i) logit^{-1} (\bar{\mathbf{X}}^T_{r_i} \hat{\beta}+ \hat{u}_i)$
- $\hat{P}^{EPP2.1}_i = f_i p_i + logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\beta} + \hat{\mathbf{u}}_i)$
- $\hat{P}^{EPP2.2}_i = logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\beta} + \hat{\mathbf{u}}_i)$
- $\hat{P}^{EPP3}_i = f_i p_i + (1-f_i) logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\alpha}+ \hat{u}_i)$
- $\hat{P}^{EPP3.1}_i = logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\alpha}+ \hat{u}_i)$
- $\hat{P}^{EPP.A}_i = logit^{-1} (\bar{\mathbf{X}}^T_i \hat{\lambda}+ \hat{a}_i)$
3. MSE推定[略]
4. シミュレーション
というわけでシミュレーションをふたつやります。(1)モデル・ベース・シミュレーション。データ生成モデルをつくって母集団を都度生成する。(2)デザイン・ベース・シミュレーション。実データを母集団にみたてて標本抽出する。
まず評価基準を定義しておく。
$k$回目の反復における小地域$i$の割合を$p_{ik}$とする[直接推定量ではなくて母平均のことらしい。よってデザインベースの場合は$k$を通じて同一]。その推定値を$\hat{p}_{ik}$, MSE推定値を$mse_{ik}$とする。$p_{ik}$の全$K$回を通じた平均を$P_i$とする[デザインベースの場合は$p_i$と等しいってことでしょうね]。
推定量のMSEの推定値を$mse_{ik}$, 真値を$MSE_i = mean_k ( (\hat{p}_{ik} - p_{ik})^2 )$とする。
[以下、簡略のため $(1/K)\sum_k^K(\cdot)$を$mean_k (\cdot)$と略記する]
- 予測量の相対バイアス：
$\displaystyle RB(P)= mean_i \left( \frac{mean_k(\hat{p}_{ik})}{P_i} - 1 \right) \times 100$
[よくわからない。この式だと、小地域予測量$\hat{p}_{ik}$の反復を通じた平均を、母集団の反復を通じた平均$P_i$をベースラインに評価していることになる。デザイン・ベースの場合は分かるけど、モデル・ベースの場合、母集団生成モデルがどうであろうが各反復における真値は$p_{ik}$なんだから、ベースラインは$p_{ik}$であるべきなのでは？ つまり$mean_i (mean_k (\hat{p}_{ik} / p_{ik}) -1) \times 100$のほうが良いのでは...？] - 予測量の相対RRMSE：
$\displaystyle RRMSE(P) = mean_i \left( \frac{\sqrt{(mean_k((\hat{p}_{ik} - p_{ik})^2)}}{P_i} \right) \times 100$
[これもよくわからん... 真値を既知としたRMSEの地域間平均$mean_i \left( \sqrt{mean_k ((\hat{p}_{ik} - p_{ik})^2)} \right) \times 100$じゃだめなの？$P_i$で割っているのは異なる地域を通じて平均するためのスケーリングっていうことかなあ...] - 予測量のMSEの相対バイアス：
$\displaystyle RB(M) = mean_i \left( \frac{mean_k(mse_{ik} - MSE_i)}{MSE_i} \right) \times 100$
[あー、ここで$MSE_i$で割ってるのはスケーリングだわ。上で$P_i$で割ってたのもスケーリングであったか。でも $RB(P)$の式は依然として謎だなあ] - 予測量の信頼区間のカヴァレッジレート：
$CR(M) = mean_i \left( mean_k ( I(|\hat{p}_{ik} - p_{ik}| \leq 2\sqrt{mse_{ik}}) ) \right)$
[これはわかりやすい。予測誤差が正規とは限らないので、95%になるはずとはいえないけど...] - 予測量のMSEの相対RMSEのパーセント平均：
$\displaystyle RRMSE(M) = mean_i \left( \sqrt{mean_k \left( (\frac{mse_{ik} - MSE_i}{MSE_i})^2 \right)}\right)$
シミュレーション１，モデルベース・シミュレーション。
$D=30, N_i=500$とする。母集団生成モデル(後述)に従って母集団を生成し、地域を層にした層別無作為抽出で標本を抽出し、いろんな推定量で推定する。以上を$K=1000$回繰り返す。
分析モデルは：
$y_{ij} \sim Binomial(1, \mu_{ij})$
$\mu_{ij} = logit^{-1} (1+x_{ij} + u_i)$
$u_i \sim N(0, \sigma^2)$ (iid)
$x_{ij} \sim \chi^2(1)$ (iid)
要因は：
- $\sigma^2$: $\{0.05, 0.10, 0.20, 0.25\}$
- $n_i$: $\{5, 10, 20, 50\}$
- 母集団生成モデル: {(1)分析モデルと同じ, (2)実は$\mu_{ij} = logit^{-1} (1+x_{ij}+x^2_{ij} + u_i)$, (3)実は$u_i \sim N(1, \sigma^2)$, (4)実は(2)(3)の両方}。
[他にMSE評価のために別のシミュレーションをやっているけど、省略]
結果。
[ここでは考察はなく、表のみ示されている。表を眺めていて結構びっくりしたので、メモしておく。
推定量のパフォーマンス(Table 1), $n_i=5, \sigma^2=0.25$をみると、まずEBLUP.FHはバイアスRB(P)が法外に大きいので除外してよい。他の推定量について分散RRMSE(P)を比べると、DIRが20.3, EPPが9.2。肝心の提案推定量は、EPP1, EPP2系, EPP3系がいずれも10から12なのに対して、なんとEPP.Aの成績が良い(9.56)。かつ、EPP2系はなぜかバイアスRB(P)が大きくなるのに対して、EPP.AはバイアスにおいてもDIRと遜色ない小ささを誇る。つまり、EPP.A (地域レベルGLMM)の成績は、EPP(ユニットレベルGLMM)に匹敵するのである。
まじか。それって、仮に共変量の母平均$\bar{\mathbf{X}}_i$が手に入っている場合(たとえば小地域が国勢調査小地域、共変量が性年代だというような場合)、非標本の共変量がわからないならば、たとえ標本の共変量をユニットレベルで測っていたとしても(ケース2)、それを捨てて地域レベルGLMMを組んだほうがよい、ということだ。言い換えると、$y_{ij}$と$x_{ij}$の地域内での共分散は捨てて良いということだ。不思議だ...
このメモは論文を一通り読み終えてから書き直しているので、ここで感想をメモしておくと、上記の知見がどこまで一般化可能なのかという点について非常にもやもやした気持ちが残る。頭が整理できていないんだけど、生成モデルにおいて共変量が地域と直交しているという点がどうもひっかかる。このとき、$y_{ij}$と$x_{ij}$の地域内での共分散は存在するいっぽう、($N_i$が十分に大きければ)小地域特性の予測に共変量は効かないわけだ。だからどうだと論理的に説明できないんだけど、これってすごく特殊な状況なのではないかと...]
シミュレーション２，デザインベース・シミュレーション。
インドのNational Sample Survey Office (NSSO)による、ビハール州のdebt investment についての標本調査データを使う。38地域, 母集団サイズは52～140。ここから20件づつ層別無作為抽出して分析する。
目的変数は「世帯が貸付残高を持っているか」。共変量は所有している土地の広さ。
[面倒になってきたのでここからのメモは省略するが、動かすのは標本サイズだけである]
結果の考察。
推定量についてみると、
- DIRはバイアスが小さくRMSEが大きい。
- 間接予測のなかではEPPのRMSEがいちばん小さい。
- ケース1で使えるのはDIRとEPP1だが、EPP1のRMSEはDIRより小さいし、バイアスもEPPに匹敵する。
- EPP2系はどれもだいたい同じ。EPP3系もどれもだいたい同じ。
- ケース2で使えるのはDIR, EPP1, EPP2系, EPP.A, EBLUP.FHだが、バイアスもRMSEも良いのはEPP.A。
- EPP2系よりEPP3系のほうが良い。
- EBLUP.FHは小標本のときひどいが、$n_i=50$あたりでEPP2系を追い越しEPP3に追いつく。
- ここまでをまとめていうと、ケース1ではEPP1が、ケース2と3ではEPP.Aがよい。
- モデルを誤指定した場合でも上記の知見は変わらない。
MSEについてみると、
- EPP1のMSE推定は基本的に過小。
- 解析的MSE推定のほうがブートストラップMSE推定より安定している。
5. 実データへの適用[略]
6. 結論
提案手法はEPPのかわりとして使えることがわかった。EPP1のMSEも使えることがわかった。
今後の課題: 空間モデルへの拡張、よりよいMSE推定の開発。
... やれやれ、疲れた。
論文の内容に対する疑問をメモしておくと：(1)ベイジアンならどうなるのか。これは自分で試せって話でしょうね。(2)上にメモしたように、なぜEPP.Aの成績がよいのかという点に興味を惹かれる。
実のところ、この論文を読んだ最大の理由は、すべてのシミュレーションのRコードが公開されているという点にある。そういうの、ほんとに助かります。まだちゃんとみてないけど、lme4パッケージを使っているようだ。
論文：データ解析(2018-) - 読了：Chandra, Kumar, Aditya (2018) 二値データの小地域推定