はじめに

小地域推定とスタイン推定量について興味を持ったので、簡単にメモをまとめます。

参考になったのは、これらの資料やページです。

小地域推定についてのノート [理論編]

小地域推定のためのノート [実習編]

推定における縮小法の展開—高次元解析と小地域推定—

多変量線形混合モデルと小地域推定への応用

縮小推定のはなし.pdf - Speaker Deck

読書日記: 覚え書き：スタイン推定量とはなんぞや

Small area estimation II

Small area estimation with mixed models: a review

スタイン推定量について

一般的な推定では推定量の最小二乗誤差が小さくなるようにパラメータ $\theta$ を推定します。

$\displaystyle{ E \left [\sum_{i=1}^n (\theta_i^* ({y_i})-\theta_i)^2 \right ] }$

一方、スタイン推定量は次のように定義される。

$\displaystyle{ \hat \theta^{s}_i =(1-\alpha ) y_i+\frac{\alpha}{n} \sum_{i=1}^n y_i }$

$\displaystyle{ \alpha=\frac{\sigma^2}{s^2} }$

$\displaystyle{ s^2 =\frac{1}{n-3}\sum_{i=1}^n \left (y_i-\frac{1}{n}\sum_{j=1}^ny_i \right )^2 }$

これは、各 $y_i$ を平均値の方向に $\alpha$ だけ引っ張ることを意味している

$\alpha$ は推定量 $\hat \theta _i ^s$ に全体平均をどの程度混ぜるかを意味している。

この $\alpha$ は、期待値のばらつき $s ^2$ が大きくなると０に近づくことになり、推定値の縮小が大きくなる。

一方で、期待値のばらつきが小さい場合は、 $\alpha$ が大きくなり縮小効果が大きくなる。

これにより、推定値の二乗誤差の期待値が小さくなること期待され証明もされている。

この時の誤差が小さくなるとは個別の推定精度が上がるというわけではなく、サンプル全体における誤差が小さくなるということに気をつけなければならない。

スタイン推定量とベイズ推定の関係

このスタイン推定量は経験ベイズ推定量に一致することが知られている。

$\theta_i$ の事前分布は次のように仮定する。

$\displaystyle{ p(\theta_i|\delta^2, \theta_0) = \frac{1}{\sqrt{ 2 \pi \delta^2}} exp \left( -\frac{(\theta_i-\theta_0)^2}{2 \delta^2} \right) }$

ベイズの公式から $\theta_i$ の事後分布は次のようになる。

$\displaystyle{ p(\theta_i|y_i, \delta^2, \theta_0) = C exp \left [ -\frac{(y_i-\theta_0)^2}{2 \sigma^2} -\frac{(\theta_i-\theta_0)^2}{2 \delta^2} \right ] }$

ここで $C$ は正則化定数ある。

ここから $\theta_i$ のMAP推定量は次のようになる。

$\displaystyle{ \theta_i^{MAP}=(1-b)y_i+b \theta_0,　b=\frac{\sigma^2}{\delta^2+\sigma^2} }$

事前分布の平均値 $\theta_0$ と分散 $\sigma ^2$ は未知であるため、経験ベイズ法による近似をして推定を行う。

その結果、平均値と分散の推定量は次のようになる。

$\displaystyle{ \hat \theta_0=\frac{1}{n} \sum_{i=1}^n yi }$

$\displaystyle{ \hat \delta^2 + \sigma^2 =\frac{1}{n}\sum_{i=1}^n \left (y_i-\frac{1}{n}\sum_{j=1}^n y_i \right )^2 }$

これを上記のMAP推定量に代入すると次のようになる。

$\displaystyle{ \hat \theta_i^{EB} =(1-c)y_i+\frac{c}{n}\sum_{i=1}^n \left (y_i-\frac{1}{n}\sum_{j=1}^n y_i \right )^2 }$

$\displaystyle{ c=\frac{\sigma^2}{\hat s^2} }$

$\displaystyle{ \hat s^2 =\frac{1}{n}\sum_{i=1}^n \left (y_i-\frac{1}{n}\sum_{j=1}^ny_i \right )^2 }$

この推定量は、ほぼスタイン推定量と一緒であり、 $\hat s ^2$ の分数の分母が異なるのみであるが、等しい推定量であることが証明されている。

Rで縮小推定

Rで二項分布の場合の縮小推定を行うためにはebbr packageがあります。

github.com

> trial_df
# A tibble: 10 x 6
# Rowwise: 
    prob samplesize     n  rate  lower upper
   <dbl>      <int> <int> <dbl>  <dbl> <dbl>
 1  0.2         172    35 0.203 0.146  0.271
 2  0.23        117    23 0.197 0.129  0.280
 3  0.25         31    11 0.355 0.192  0.546
 4  0.22        229    52 0.227 0.174  0.287
 5  0.3         300    94 0.313 0.261  0.369
 6  0.22         98    26 0.265 0.181  0.364
 7  0.26        287    70 0.244 0.195  0.298
 8  0.19         81    17 0.210 0.127  0.315
 9  0.12        232    31 0.134 0.0926 0.184
10  0.24        136    27 0.199 0.135  0.276

スタイン推定量と通常の平均値のMSEを見てみます。

スタイン推定量の方がMSEが小さいです。

> mean(sae_trial_df$stain_se)
[1] 0.000524215
> mean(sae_trial_df$obs_se)
[1] 0.001696015

各データの結果を可視化してみます。

# 描画
sae_trial_df %>%
  ggplot(aes(.fitted, rank(.fitted))) +
  geom_point() +
  geom_point(aes(x = prob), color = "red") +
  geom_point(aes(x = rate), color = "blue") +
  geom_errorbarh(aes(xmin = .low, xmax = .high, height = .4)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper, height = .4), col = "blue") +
  labs(x = "推定量",
       y = "対象",
       title = "縮小推定結果",
       subtitle = "黒:推定量, 青:観測値, 赤:真値")+
  theme_gray (base_family = "HiraKakuPro-W3")

f:id:saltcooky:20210227162831p:plain

一番上振れしたデータがあり、平均値への縮小がうまく聞いているようです。

小地域推定

人口統計などを考えるとき、狭い地域や人口が粗な地域に対しては十分なサンプルサイズが得られない。

このような地域では、その地域のデータだけでは統計量の十分な推測ができず、統計量の分散が大きくなってしまう。

この時の推定量は、得られている情報のみを利用しているため直接推定量と呼ばれている。

このような状況での推定問題を小地域推定とよび、周辺地域のデータを組み込んで推定精度を高めることができる。

代表的なモデルとして、Fay-Herriotモデルや枝分かれ誤差回帰モデル (nested error regression model, NERM) が挙げられる。

Fay-Herriotモデル

$i$ 番目の小地域の観測データを $y_i$ , 共変量を $x_{i} = ( 1, x_{i} ^1, \dots , x_{i} ^{p−1} ) ^T$ とすると、Fay-Herriotモデルは次のように定義される。

$\displaystyle{ y_{i} = x^T_{i} \beta + v_i + \epsilon_{i} ,　 i = 1, \dots , m }$

ここで、 $v_i$ は地域の差異を表す変量効果で $v_i \sim N (0, \tau^2 )$ に従い，それと独立に誤差項 $\epsilon_i$ が $\epsilon_{i} \sim N(0, \sigma^2 )$ に従う。この時の $\sigma^2$ は未知である。

各地域の特性値が特徴量 $x_i$ を用いて $\xi_i = x^T_i \beta + v_i$ と表されるとき，この値を予測する問題を考える．

このとき，小地域のデータ $y_i$ を与えたときの $x_i$ の条件付き期待値を求めると、次のようにかける。

$\displaystyle{ \hat \xi_i^B = E [ \xi_i | y_i ] = x_i^T \beta+ E[v_i|y_i] = x_i^T \beta + \frac{\tau^2}{\tau^2 + \sigma^2}( y_i - x_i^T \beta) }$

これはベイズ推定量に対応し、 $x_i ^T \beta$ の方向へ縮小する作用が得られることになる。

$\rho = \tau ^2 / \sigma ^2$ 遠くと $\beta$ の一般化最小二乗推定量（GLS）は次のようになる。

$\displaystyle{ \tilde{\mathbf{\beta}}_{GLS} = \left[ \sum_i^m \frac{x_i y_i}{\tau^2 + \sigma^2} \right] \left[ \sum_i^m \frac{x_i x_i^t}{\tau^2 + \sigma^2} \right] ^{-1} }$

これを $\hat \xi _i ^B$ に代入すると次のようにある。

$\displaystyle{ \hat \xi_i^{BLUP}(\rho) = c_i^T \tilde{\mathbf{\beta}}_{GLS} + \frac{n_i \rho}{1 + n_i \rho}( y_i - x_i^T \tilde{\mathbf{\beta}}_{GLS}) }$

最良線形不偏予測量 (best linear unbiased predictor, BLUP) とも呼ばれている。

枝分かれ誤差回帰モデル

このモデルでは、一地域から複数の観測値が得られている場合のモデルである。

$i$ 番目の小地域の $j$ 番目のデータを $y_{ij}$ , 共変量を $x_{ij} = (1, x_{ij}^1, . . . , x_{ij}^{p−1})^T$ とし，全体で $m$ 個の小地域から構成され，各小地域から $n_i$ 個のデータが観測されているとする。

この時の枝分かれ誤差回帰モデルは次のようになる。

$\displaystyle{ y_{ij} = x^T_{ij} β + v_i + \epsilon_{ij} ,　 i = 1, \dots , m, j = 1, \dots , n_i }$

ここで、 $v_i$ は地域の差異を表す変量効果で $v_i \sim N ( 0, \tau ^2 )$ に従う確率変数とし、それと独立に誤差項 $\epsilon _{ij}$ が $\epsilon _{ij} \sim N(0, \sigma^2)$ に従う。

各地域の特性値を p次元のベクトル $c_i$ を用いて $\xi _i = c^T_i \beta + v_i$ と表される。

例えば $c_i = \bar x_i = n^{−1}_i \sum ^{n_i}_{ j=1} x_{ij}$ とおくと $\xi_i = \bar x ^T_i \beta + v_i$ は平均に関する量を推定することになる。

このときの小地域のデータの組 $y_i = (y_{i1}, \dots , y_{in_i})^T$ を与えたときの $x_i$ の条件付き期待値を求めると次のようにかける。

$\displaystyle{ \hat \xi_i^B = E [ \xi_i | y_i ] = c_i^T \beta+ E[v_i|y_i] = c_i^T \beta + \frac{\tau^2}{\tau^2 + \sigma^2/n_i}(\bar y_i - \bar x_i^T \beta) }$

これはベイズ推定量に対応し、 $\bar y_i$ と $\bar x_i ^T$ の方向へ縮小する作用が得られることになる。

$\rho = \tau^2 / \sigma^2$ 遠くと $\beta$ の一般化最小二乗推定量（GLS）は次のようになる。

$\displaystyle{ \hat \beta (\rho) = \left ( \sum_{i=1}^m \frac{n_i \bar x_i \bar x_i ^T}{1+n_i\rho} + \sum_{i=1}^m \sum_{j=1}^n x_{ij}x_{ij}^T \right)^{-1} \left (\sum_{i=1}^m \frac{n_i \bar x_i \bar y_i ^T}{1+n_i\rho} + \sum_{i=1}^m \sum_{j=1}^n x_{ij}y _{ij}^T \right ) }$

そして、最良線形不偏予測量 (best linear unbiased predictor, BLUP) は次のようになる。

$\displaystyle{ \hat \xi_i^{BLUP}(\rho) = c_i^T \hat \beta (\rho) + \frac{n_i \rho}{1 + n_i \rho}(\bar y_i - \bar x_i^T \hat \beta(\rho)) }$

実際には、(一般化)線形混合モデルを用いて推定することになる。

RでFay-Herriotモデル

いくつかパッケージがありますが、BayesSAEパッケージを用いました。

以下のモデルに従うサンプルデータを作成しました。

$\displaystyle{ \theta_i = 0.23x_1 + 0.32x_2 + 0.15x_3 + v_i + \epsilon_i }$

　クソコードですね。

set.seed(1234)
N <- 20
beta <- c(0.23, 0.32, 0.15, 1)

model_df <- 
  data.frame(x1 = rgamma(N, 1.2),
             x2 = rgamma(N, 2.1),
             x3 = rbernoulli(N),
             v1 = rnorm(N, 0, 0.1)) %>% 
  mutate(theta = beta[1]*x1 + beta[2]*x2 +beta[3]*x3 + beta[3]*v1) 

sample_df <- 
  model_df %>% 
  mutate(sample_size = 2*floor(rnorm(n=N,30,25))) %>% 
  tibble() %>% 
  mutate(theta_mu = map2(.x = sample_size, .y = theta, ~ rnorm(.x, .y, 0.2))) %>% 
  mutate(Y = map(sample_df$theta_mu, mean) %>% unlist(),
         sd = map(sample_df$theta_mu, sd) %>% unlist(),
         Var = sd^2)

平均値と補助変数を用いた場合のMSEを見てみました。

今回の場合では、補助変数1つの場合が最も小さいという結果になりました。

> # 平均値のMSE
> mean((sample_df$theta - sample_df$Y)^2)
[1] 0.9035524
> # 補助変数1つの場合の縮小推定のMSE
> mean((sample_df$theta - summary(bsae_res)$theta)^2) 
[1] 0.8464218
> # 補助変数3つの場合の縮小推定のMSE
> mean((sample_df$theta - summary(bsae_val_res)$theta)^2)
[1] 0.8622283

補助変数3つの場合、推定結果をグラフ化してみます。

# 描画
sample_df %>% 
  mutate(fitted = summary(bsae_res)$theta) %>% 
  ggplot(aes(fitted, rank(theta))) +
  geom_point(aes(x = theta), color = "red") +
  geom_point(aes(x = Y), color = "black", alpha = 0.5) +
  geom_point(aes(x = fitted), color = "blue") +
  geom_errorbarh(aes(xmin = Y-1.96*sd, xmax = Y+1.96*sd, height = .4), color = "black", alpha = 0.5) +
  labs(x = "推定量",
       y = "対象",
       title = "小地域推定結果",
       subtitle = "黒:観測平均値, 青:推定量, 赤:真値")+
  theme_gray (base_family = "HiraKakuPro-W3")