tidybayesパッケージで推定結果の整然化

標準化死亡比

標準化死亡比 (Standardized Mortality Ratio)は、地域間での死亡率を比較する指標としてよく用いられる指標です。

一般的には、当該地域の実際の死亡数\(d\)と、全国の年齢階級別死亡率と当該地域の死亡数かけて算出される期待死亡数\(exp_d\)の比をとって100かけた値,

\(\frac{d}{exp_d}\times 100\),

がSMRになります。

100であれば全国平均と同一、100以下であれば全国より当該地域の死亡率が低い、100以上であれば全国より当該地域の死亡率が高いことを示しています。
なお、SMRは人口規模の小さい地域では変動が大きいという欠点があり、それを補うために、階層ベイズモデルを用いた縮約統計量の利活用が推奨されています。

例えば、[関西大学の紺田先生による市町村レベルの標準化死亡比の階層モデルについての文献が勉強になります。こちらやこちら

市区町村の標準化死亡比を推定するpoisson-gammaモデル

ここでは、上記の紺田先生のモデルに基づいて標準化死亡比を階層モデルを使って推定します。

• 地域\(i\)の自殺発生数: \(d_i\)
  
• 地域\(i\)の期待死亡数: \(exp\_d_i\)
  
• 地域\(i\)の潜在的な標準化死亡比: \(\theta_i\)

possison-gamma モデル

\[d_i \sim Poisson(exp\_d_i \theta_i)\]

\(\theta\)の事前分布に、ガンマ分布指定し、

\[ \theta_i \sim Gamma (\alpha, \beta) \]

ガンマ分布の形状パラメータ(\(\alpha\))と尺度パラメータ(\(\beta\))には、指数分布とガンマ分布を指定します。

\[\alpha \sim Exp (b_\alpha = 1), \beta \sim Gamma (a = 0.1, b_\beta = 1)\]

stan コード

このpoisson-gammaモデルの推定を行うためのstanコードは、以下のようになります。

data{
  int r; // 地域数
  int d[r]; //各地域の自殺発生数
  real exp_d[r];　//各地域の期待死亡数
}
parameters{
  real <lower=0> theta[r];
  real <lower=0> alpha;
  real <lower=0> beta;
}

model {
  for(i in 1:r){
    d[i] ~ poisson(exp_d[i]*theta[i]);
    theta[i] ~ gamma(alpha,beta);
  }
  alpha ~ exponential(1);
  beta ~ gamma(0.1,1);
}

generated quantities {
  real RRmean; 
  real RRvar;
  RRmean = alpha/beta; // 市区町村の自殺相対リスクの県平均
  RRvar = alpha/(beta^2); // 市区町村の自殺相対リスクの分散
}

generated quantitiesブロックでは、gamma分布の性質利用して\(\alpha\)と\(\beta\)か県レベルの相対リスクの平均と分散を生成量として求めています。

推定の実行

library(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

modelGP<-stan_model("Sdat/HPGmodel.stan")

fit<-sampling(modelGP,data=list(d=dat$d,exp_d=dat$exp_d,r=nrow(dat)))

spread_draw関数

tidybayesパッケージでは、ほぼ全ての場合において、ベイズ推定結果が格納されているオブジェクトをこのspread_draws関数に流すことから始まります。stanオブジェクトをspread_draws関数に落とすと下記のような整然データが得られます。

library(tidyverse)
library(tidybayes)

fit %>%  spread_draws(theta[r], alpha,beta) %>% head(20)

## # A tibble: 20 x 7
## # Groups:   r [20]
##    .chain .iteration .draw     r theta alpha  beta
##     <int>      <int> <int> <int> <dbl> <dbl> <dbl>
##  1      1          1     1     1 1.04   5.74  4.71
##  2      1          1     1     2 1.29   5.74  4.71
##  3      1          1     1     3 2.20   5.74  4.71
##  4      1          1     1     4 1.09   5.74  4.71
##  5      1          1     1     5 0.941  5.74  4.71
##  6      1          1     1     6 1.01   5.74  4.71
##  7      1          1     1     7 1.58   5.74  4.71
##  8      1          1     1     8 1.05   5.74  4.71
##  9      1          1     1     9 1.10   5.74  4.71
## 10      1          1     1    10 1.28   5.74  4.71
## 11      1          1     1    11 0.828  5.74  4.71
## 12      1          1     1    12 1.32   5.74  4.71
## 13      1          1     1    13 0.875  5.74  4.71
## 14      1          1     1    14 1.17   5.74  4.71
## 15      1          1     1    15 0.480  5.74  4.71
## 16      1          1     1    16 1.26   5.74  4.71
## 17      1          1     1    17 0.574  5.74  4.71
## 18      1          1     1    18 1.95   5.74  4.71
## 19      1          1     1    19 0.988  5.74  4.71
## 20      1          1     1    20 1.46   5.74  4.71

chainごと、iterationごと、drawごと、地域(r)ごとのパラメータの推定結果が一つのtibble形式のデータフレームとして整理されます。パラメータが複数ある場合には、gather_draws関数を使って、変数群をlongフォーマットで整理することもできます(ggplotで可視化するときに便利です)。

事後分布の要約統計量系関数

tidybayesパッケーでは、spread_drawsもしくはgather_draws関数で整然化されたデータを下記の関数に流すことで、事後分布のようやく統計量の整然データを得ることができます。中央値、平均値、最頻値について、それぞれパーセンタイルに基づく方法と最高密度区間に基づく方法を選べます。

• median_hdci: 中央値で最高密度区間
• mean_hdci: 平均値で最高密度区間
• mode_hdci: 最頻値で最高密度区間
• median_qi: 中央値でパーセンタイル
• mean_qi：平均値でパーセンタイル
• mode_qi：最頻値でパーセンタイル

試しに各市町村の潜在SMR\(\theta\)の事後分布の中央値と最高密度区間を抽出して見ます。

fit %>% spread_draws(theta[r]) %>%
   median_hdci()

## # A tibble: 59 x 7
##        r theta .lower .upper .width .point .interval
##    <int> <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
##  1     1 0.778  0.556  0.999   0.95 median hdci     
##  2     2 1.00   0.668  1.41    0.95 median hdci     
##  3     3 2.22   1.71   2.78    0.95 median hdci     
##  4     4 1.01   0.784  1.26    0.95 median hdci     
##  5     5 1.28   0.746  1.85    0.95 median hdci     
##  6     6 1.14   0.699  1.62    0.95 median hdci     
##  7     7 1.06   0.546  1.64    0.95 median hdci     
##  8     8 1.29   0.650  1.99    0.95 median hdci     
##  9     9 1.10   0.622  1.65    0.95 median hdci     
## 10    10 1.46   0.824  2.18    0.95 median hdci     
## # … with 49 more rows

簡単にstanオブジェクトから整然データが得られました。地域データなので、せっかくなので地図で可視化したいと思います。

library(jpndistrict)
fukushima<-jpndistrict::jpn_pref(pref_code="7")

theta_bayes<-fit %>% spread_draws(theta[r]) %>% median_hdci()

fukushima2<-data.frame(fukushima,theta_bayes,dat)

fukushima2 %>% ggplot() +
  geom_sf(aes(geometry = geometry,fill=theta),col = "white" )+
  theme_bw()

dplyrの関数がパイプ演算子でそのまま使えるので、全ての市町村から市のデータ(r=1:13)だけを抽出したい場合に、filter関数で抽出することができます。

fit %>% spread_draws(theta[r]) %>% 
  filter(r==grep("市",dat$ct_name))%>%
   median_hdci()

## # A tibble: 13 x 7
##        r theta .lower .upper .width .point .interval
##    <int> <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
##  1     1 0.765  0.582  0.997   0.95 median hdci     
##  2     2 0.992  0.673  1.36    0.95 median hdci     
##  3     3 2.22   1.65   2.66    0.95 median hdci     
##  4     4 0.997  0.798  1.27    0.95 median hdci     
##  5     5 1.24   0.849  1.85    0.95 median hdci     
##  6     6 1.15   0.710  1.63    0.95 median hdci     
##  7     7 1.05   0.557  1.53    0.95 median hdci     
##  8     8 1.24   0.665  2.03    0.95 median hdci     
##  9     9 1.09   0.623  1.60    0.95 median hdci     
## 10    10 1.50   0.884  2.11    0.95 median hdci     
## 11    11 1.24   0.793  1.84    0.95 median hdci     
## 12    12 1.25   0.699  1.76    0.95 median hdci     
## 13    13 1.03   0.511  1.76    0.95 median hdci

事後分布のプロット関数

• geom_eyeh:事後分布を両側に出します。

• geom_harlfeyeh: 事後分布が片側だけに出ます

どちらも事後分布と共に、点推定値と確信区間が表示されます。geom_eyeh,geom_harlfeyehの引数、point_intervalで点推定値と確信区間にどれを使うか指定できます。また引数.widthで確信区間の区間幅も指定可能です。

fit %>% spread_draws(theta[r]) %>% 
  filter(r==grep("市",dat$ct_name)[1:2])%>% ggplot()+geom_eyeh(aes(x=theta,y=as.factor(r)),point_interval=median_hdci,.width=c(.66,0.96))+ylab("City")

fit %>% spread_draws(theta[r]) %>% 
  filter(r==grep("市",dat$ct_name)[1:2])%>% ggplot()+geom_halfeyeh(aes(x=theta,y=as.factor(r)),point_interval=median_hdci,.width=c(.66,0.96))+ylab("City")

モコモコしたくない人は区間推定のみのプロットも選択可能

• stat_pointintervalh:spread_drawsやspread_drawsでできた整然データをこの関数に流すことで点推定値と確信区間がプロットできます。geom_eye関数のときと同様に、引数point_intervalで点推定値と確信区間にどれを使うか指定できます。また引数.widthで確信区間の区間幅も指定可能です。

#市町村(r) 1:5までを表示
fit %>% spread_draws(theta[r]) %>% 
  filter(r==grep("市",dat$ct_name)[1:5])%>% ggplot()+stat_pointintervalh(aes(y = r, x = theta),point_interval=median_hdci,.width=c(.66,0.96))

パラメータの条件間比較

• compare_levels: この関数では,各条件の差の事後分布を求めてくれます。

ここでは、地域\(r\)(1から3)間のSMR(\(\theta\))の差を求めます。iterationごとに、条件間でのthetaの推定値の差を求める関数です。この結果をgeom_eyehやstat_intervalなどこれまで使い方を示してきた関数に流すと同様のプロットができます。美しい。

fit %>%
  spread_draws(theta[r]) %>%
  filter(r<=3) %>%
  compare_levels(theta, by = r) %>%
  ggplot(aes(y = r, x = theta)) +
  geom_halfeyeh()

用途不明なドットプロット

かわいい以外に使い道が見出せずにいます本家のページに紹介されていたので紹介しておきます。かわいいです。

fit %>%
  spread_draws(theta[r]) %>% 
  filter(r==grep("市",dat$ct_name)[1:2])%>% do(data_frame(theta = quantile(.$theta, ppoints(200)))) %>%
  ggplot(aes(x = theta)) +
  geom_dotplot(binwidth = .02) +
  facet_grid(r ~ .) +
  scale_y_continuous(breaks = NULL)

他のパッケージとの橋渡し

tidybayesパッケージでは、収束診断系の統計指標やプロットのための関数が用意されていません。それらはbayesplotで既にいい感じのんがあるからそっちにスムースにデータを送れるようにするのがいい、ということらしいです。

• unspread_draws関数は, spread.draws関数で整然化したデータarray型に戻す関数です。これを使うと、array型のデータを読み込むbayesplotの関数にデータを送れます。

例えば、一度’spread_draws’で整然化した後に、chain1だけが振る舞いがおかしかったので、その一つだけを外して自己相関などの収束結果を確認したいというとき、下記のようにできます。

# filder 関数でchain 3だけ外したあとで、unspread_drawsでアレイ型にデータを戻してbaysplotのmcmc_acf関数にデータを流す

library(bayesplot)

fit %>%
  spread_draws(theta[r]) %>% filter(.chain!=3) %>% unspread_draws(theta[r],drop_indices=T) %>% mcmc_acf(pars=c("theta[1]","theta[2]"))

ちなみにトレースプロットを書きたかったら、spread_drawsからそのままggplotで軸を指定すればOKです。

fit %>%
  spread_draws(theta[r]) %>% filter(.chain!=3,r==c(1:3)) %>% 
  ggplot(aes(x=.iteration,y=theta,group=as.factor(.chain),color=as.factor(.chain)))+geom_line(alpha=0.5)+facet_wrap(~r,ncol=1)