はじめに
千葉大学・株式会社Nospareの川久保です．
以前の記事で，線形混合モデルとその応用例をいくつか紹介しました．今回は，応用例のうち小地域推定（small area estimation）について詳しく解説していきたいと思います．
小地域推定とは
まず，小地域推定とは何を目指しているのかを，簡単に説明します．ある標本調査にもとづいて各市区町村における何らかの量（平均所得・貧困率・疾病リスクなど）を測ろうとする問題を考えます．この標本調査は，日本全国で見れば十分に大きなサンプルサイズであったとしても，各市町村においてはサンプルサイズが小さかったりゼロであったりすることがしばしばあります．サンプルサイズが小さいもとで，その地域の標本調査のデータのみにもとづいて興味のある量を推定しようとすると，推定が安定しません（推定量の分散が大きくなります）．このように，粒度の細かい情報を小さいサンプルサイズのデータから推定しようとする際，適切な統計モデルを仮定して推定精度を改善しようとする試みがなされます．こうした推定問題を，小地域推定と呼んでいます．
「仮定する適切な統計モデル」として有用であると広く認識されているのが，以前の記事でも解説した線形混合モデルと呼ばれるモデルのクラスです．線形混合モデルのうち，小地域推定で最もよく利用されている基本的な2つのモデルを次に紹介します．
小地域推定の基本的な2つのモデル
地域レベルモデル
$y_i$を地域$i$の興味のある量の直接的な推定量（direct estimator）であるとします．例えば，興味のある量が平均所得であるとき，$y_i$はある標本調査（家計調査など）において地域$i$でサンプルとして観測された所得たちの標本平均が考えられます．しかし先に述べたとおり，地域$i$におけるこの調査のサンプルサイズが小さい場合，$y_i$は推定精度がよくありません．そこで以下のような統計モデルを仮定します:
$$
y_i = x_i^\top \beta + b_i + \varepsilon_i, \quad b_i \sim \mathrm{N}(0,\tau^2), \quad \varepsilon_i \sim \mathrm{N}(0,D_i),\tag{1}
$$
ここで$x_i$は地域$i$において利用可能な補助変数ベクトル，$b_i$は変量効果として入れた地域効果，$\varepsilon_i$は標本誤差で，その分散$D_i$は既知と仮定します．$y_i$として標本平均を考えた場合は，$D_i$は地域$i$のサンプルサイズ$n_i$に反比例すると考えられます．このモデルを利用する場合は，三つ組$(y_i,x_i,D_i)$が必要です．このモデルはFay and Herriot (1979, JASA) で提案されたことから，Fay–Herriotモデルと呼ばれています．小地域推定の目的は，このモデルから$x_i^\top \beta + b_i$を予測することで，これはまさに以前の記事で解説した混合効果の予測問題で，EBLUPを用いて予測を行います．(1)式のモデルにおける$x_i^\top \beta + b_i$のEBLUPは以下です:
$$
x_i^\top\hat{\beta} + \hat{b}_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}
$$
これは，$y_i$と$x_i^\top \hat{\beta}$の加重平均になっていて，$D_i$が大きい（サンプルサイズが小さく標本誤差が大きい）地域は$x_i^\top \hat{\beta}$のウェイトが高くなり，$D_i$が小さい（サンプルサイズが大きく$y_i$の推定精度が高い）地域では直接的な推定量$y_i$のウェイトが高くなります．つまりサンプルサイズの大きい地域では安定していると考えられる$y_i$のウェイトを大きくし，サンプルサイズの小さい地域では$x_i^\top \hat{\beta}$のウェイトを大きくするよう，EBLUPはなっています．ここで$x_i^\top\hat{\beta}$は，$\hat{\beta}$を全地域のデータから推定しているため，安定していると考えられます．
(1)式のモデルは，線形混合モデルの一種ですが，以下のような階層モデルとしてとらえることも可能です:
\begin{align}
y_i \mid \theta_i &\sim \mathrm{N}(\theta_i, D_i), \\
\theta_i &\sim \mathrm{N}(x_i^\top \beta, \tau^2),
\end{align}
$\theta_i$が興味のあるパラメータ（真の平均所得）で，その周辺分布に補助変数$x_i$の情報を入れています．
個票レベルモデル
(1)式のモデルは，地域ごとに集計されたデータ$y_i$や$x_i$から構築しているのに対し，標本調査の生の個票データが利用可能である場合には，これから説明するモデルを考えることもできます．$y_{ij}$を地域$i$の$j$番目の個票における興味のある量の観測（地域$i$の$j$番目の家計の所得），$x_{ij}$を対応する補助変数ベクトルとします．このとき，
$$
y_{ij} = x_{ij}^\top \beta + b_i + \varepsilon_{ij}, \quad b_i \sim \mathrm{N}(0,\tau^2), \quad \varepsilon_{ij} \sim \mathrm{N}(0,\sigma^2), \tag{2}
$$
という線形混合モデルを仮定します．このモデルを用いた小地域推定の手法は，Battese, Harter and Fuller (1988, JASA) で提案されました．地域$i$の平均に興味がある場合，観測されていない全ての$y_{ij}$の算術平均が予測の対象となりますが，近似的に，
$$
\bar{X}_i^\top \beta + b_i
$$
を予測対象とすることが多いです．ここで$\bar{X}_i$
は，$x_{ij}$の地域$i$の母集団における算術平均です（つまり$y_{ij}$が未観測の$ij$に対しても$x_{ij}$を全部足して，地域$i$の母集団の大きさで割るということ）．
やはりこの予測問題も混合効果の予測なので，EBLUPで予測します．
(2)式のモデルは個票を用いているので，(1)式のモデルと比べて推定の効率は上がります．しかしながら，個票にアクセスできないことも多い上に，補助変数については全数または少なくとも母集団における平均$\bar{X}_i$が予測には必要なので，(2)式のモデルを使った小地域推定を行える場面は限られています．
Rによる解析例
Rではsaeというパッケージが実装されており，(1)式や(2)式のモデルの推定とEBLUPの計算を行ってくれます．(1)式のモデルはeblupFH()関数で解析できます．eblupFH()関数のhelpに掲載されているexampleを実行してみます．
library(sae)
data(milk)
attach(milk)
resultREML <- eblupFH(yi ~ as.factor(MajorArea), SD^2)
saeパッケージ内のmilkというデータセットに，(1)式のモデルを当てはめています．43地域の牛乳の消費額を推定する問題を考えていて，yiがその直接的な推定量，SDが標本誤差の標準偏差$\sqrt{D_i}$，補助変数としてmajor areaのダミー変数を考えています（同じmajor areaだと似た消費額になるという仮定）．結果はresultREMLにリストの形で格納しています．例えば各地域のEBLUPは，resultREML$eblupに43 × 1の行列として納められています．小地域推定ではあまり関心の対象とはなりませんが，回帰係数の推定値は，
> resultREML$fit$estcoef
beta std.error tvalue pvalue
(Intercept) 0.9681890 0.06936208 13.958476 2.793443e-44
as.factor(MajorArea)2 0.1327801 0.10300072 1.289119 1.973569e-01
as.factor(MajorArea)3 0.2269462 0.09232981 2.457995 1.397151e-02
as.factor(MajorArea)4 -0.2413011 0.08161707 -2.956503 3.111496e-03
で見ることができます．このデータにはmajor areaが4つあり，major area 1をreferenceとしています．
次に個票レベルモデル(2)の解析例を見てみます．
library(sae)
data(cornsoybean)
data(cornsoybeanmeans)
attach(cornsoybeanmeans)
Xmean <- data.frame(CountyIndex, MeanCornPixPerSeg, MeanSoyBeansPixPerSeg)
Popn <- data.frame(CountyIndex, PopnSegments)
resultCorn <- eblupBHF(CornHec ~ CornPix + SoyBeansPix,
dom = County, meanxpop = Xmean, popnsize = Popn,
data = cornsoybean)
cornsoybeanはsae内のデータフレームで，12地域（$i=1,\dots,12$）の37地点における，標本調査で報告されたコーンの作付面積CornHecとその補助変数たちを含んでいます．目的は，各地域における平均作付面積を推定（予測）することですが，12地域のトータルのサンプルサイズが37であることから分かるように，各地域のサンプルサイズは非常に小さくなっています．そこで(2)式のモデルを使った小地域推定を行います．cornsoybeanmeansは補助変数の各地域における平均$\bar{X}_i$のデータを含んでおり，Xmeanに必要なデータを格納しています．Popnには各地域の母集団サイズを入れています．準備が整ったところで，eblupBHF()関数を用いて解析します.12地域のEBLUPは，
> resultCorn$eblup
domain eblup sampsize
1 1 122.5825 1
2 2 123.5274 1
3 3 113.0343 1
4 4 114.9901 2
5 5 137.2660 3
6 6 108.9807 3
7 7 116.4839 3
8 8 122.7711 3
9 9 111.5648 4
10 10 124.1565 5
11 11 112.4626 5
12 12 131.2515 6
で見ることができます．
小地域推定とは（再考）
今回は小地域推定における最もシンプルなモデルを2つ紹介しましたが，これらのモデルをベースとして，
カウントデータ，二値データ，割合データ，グループデータなどの，様々なデータ形式に対応するためのそれぞれのモデル開発
空間相関や空間異質性を考慮したモデルへの拡張
時系列，時空間モデルへの拡張
など様々なモデルの拡張がなされています．例えば，カウントデータや二値データへの対応であれば一般化線形混合モデル（GLMM），空間モデルであれば線形混合モデルにおいて変量効果に空間相関の構造を入れる，などなどです．ではあらためて「小地域推定とは何か」を考えると，結局は統計学におけるこれらの一般的なモデルを，各地域の興味のある量（小地域パラメータ）を「予測」するために使っているということです．
しかし小地域推定においてより重要であると考えられているのは，EBLUPなどの予測量を使って小地域パラメータを予測した際，どの程度の誤差を持っているのかをしっかり評価することです．単に「予測が当たればいい」のではなく，予測量の平均二乗予測誤差（MSPE）や，予測区間もあわせて提示することで，小地域推定モデルから出てきた予測値をユーザーが納得して使えることを目指しています．
小地域推定に関する文献として，以下を推薦します．特に，Nospare社の菅澤先生のレビュー論文はオープンアクセスになっているため，ぜひダウンロードしてください！
まとめ
株式会社Nospareには，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．
