シンプソンのパラドックスと階層ベイズ

丸の内オアゾに入っている書店をぶらついているときに『データ分析をマスターする12のレッスン』という本を見つけました。最近出た学部3〜4年生向けの本のようです。パラパラと立ち読みしたところ，商学と政策の先生が書いた本で自分にちょうどよさそうなレベルだったのでお買い上げ。

データ分析をマスターする12のレッスン (有斐閣アルマBasic)

作者: 畑農鋭矢,水落正明
出版社/メーカー: 有斐閣
発売日: 2017/10/12
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

　
その中で，おもしろい例題がありました。各都道府県の合計特殊出生率と女性労働力率の関係です。グーグルによると，それぞれ次のとおりです。

合計特殊出生率（tfr）：一人の女性が一生に産む子供の平均数
女性労働力率（flp）： 15歳以上の女性の人口に占める、実際に働いている、もしくは求職中の女性の割合。

　
横軸を女性労働力率，縦軸を合計特殊出生率でプロットすると，次のようになります。

f:id:hikaru1122:20171023220827p:plain

　
散布図を見る限り，負の相関があります。実際に相関係数は -0.19。しかし，このデータは1980年と2000年のデータがいっしょになっていて，年代を色分けし，回帰直線を引いてみると次のようになります。青い回帰直線は年代別にしなかった場合のもの。

f:id:hikaru1122:20171023220845p:plain

　
全体と層別にしたときに結果が異なるシンプソンのパラドックスが起きています。層別にすると，女性労働力率が大きいほど合計特殊出生率が大きくなって「女性が働くほど子どもが生まれやすい」なんて説明がつきそうです。

　
『データ分析をマスターする12のレッスン』ではトレンド項というのを入れたりして，うまく解決していくのですが，年代でネストされていて，都道府県の個別事情（効果）も踏まえたマルチレベル分析でもやれるのではないかと思い，練習がてら，やってみることにしました。

ベイズ推定でやる！ちょうど犬4匹本の輪読会も開催するし

kobe-r-annex.connpass.com

　
各都道府県の年代別で線を引いてみると，そんなに傾きは変わらないようなので，ランダム切片だけに絞ります。

f:id:hikaru1122:20171023220950p:plain

　
このレベルの Stan コードはまだ自前で書けないので，brms パッケージ頼み。

prior = c(set_prior("cauchy(0, 3)", class = "sd", group = "prefID"),
           set_prior("cauchy(0, 3)", class = "sd", group = "yearID"))

library(brms)
kekka.1 = brm(tfr ~ 1 + flp + (1 | prefID) + (1 | yearID), data = d, iter = 5000, prior = prior,
                 control = list(adapt_delta = 0.9, max_treedepth = 12))

　
上段はランダム切片の事前分布です。この事前分布にしたのは，stan レファレンスマニュアル（日本語版）を参考にして設定しました。もっといい設定があれば教えてください。下段がベイズ推定をしているところ。control の部分はオプションです。実行すると次のようなエラーが出ます。アドバイスに従って，最大値まで上げても消えなかったので，まあまあの落としどころとして設定しました。

警告メッセージ: 
There were ●●●● divergent transitions after warmup. Increasing adapt_delta above 0.9 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

推定結果

　
次のとおりで，Rhat だけで判断すると収束しています（有効サンプルサイズは気にしなくていいかな？）。

> kekka.1
 Family: gaussian(identity) 
Formula: tfr ~ 1 + flp + (1 | prefID) + (1 | yearID) 
   Data: d (Number of observations: 94) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
    ICs: LOO = NA; WAIC = NA; R2 = NA
 
Group-Level Effects: 
~prefID (Number of levels: 47) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.13      0.02     0.10     0.17       1232 1.00

~yearID (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     1.22      1.21     0.15     4.48        953 1.01

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     1.69      1.09    -0.77     4.10        634 1.01
flp          -0.00      0.00    -0.01     0.01       1095 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.05      0.01     0.04     0.06       1590 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
 警告メッセージ: 
There were 402 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup