岩波DSを読んでから、自分の中でのベイズ統計学へ取り組む姿勢が変わりつつある。

岩波DSの良いところは、本の内容はもちろん、ソースコードやデータファイルが提供されている点やサポート記事（とそれに相当するブログ記事）が充実している点であり、自分が知る限り、これまでに岩波DSに関する書評などは以下のものがある（それぞれの記事の概要はだいたい@TJO_datasciさんの次のブログ記事にまとめられている）。

先人たちと同じことをしても仕方がない気がして、自分には何ができるのかと考えた際に、MCMCサンプリングを実行した後の結果の評価や可視化かな（@Med_KUさんがやっていた）、ということで勉強をかねてやってみる。幸いにも（？） https://sites.google.com/site/iwanamidatascience/vol1/support_tokushu で書かれているように、松浦さんの「Stan入門」の部分では収束判定の内容について割愛されているらしいので、StanによるMCMC計算後の収束判定や、推定したパラメータの扱いについてRを使ってやってみたい。知識が乏しいので、間違いや稚拙な表現に関してはご指摘いただければ幸いである。

とりあえず久保さんによる「階層ベイズ最初の一歩」の例題を取り上げてみる。

🔧 主要なRパッケージの読み込みとRStanの設定

library(ggmcmc)
library(rstan)

その他、いくつかのパッケージを利用している。

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

上記のコードはrstanの並列化のためのおまじない

🍙 給食の種類に応じて成績が変わる？を例に

例題については本書を参考に。ここでは本書の内容に沿ってRコードを実行する。まずデータを用意する必要があるので久保さんのウェブページからデータファイルをダウンロードしてくる（GitHubの岩波DSリポジトリにはないので注意）。

# http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiDs01.html からデータをダウンロード
download.file(url = "http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwnmDS1/data/data.RData",
              destfile = paste(getwd(), "iwanami_ds", "data.Rdata", sep = "/"))
load(file = "iwanami_ds/data.Rdata")
# データ構造の確認。
dplyr::glimpse(d)
## Observations: 20
## Variables: 6
## $ pref   (fctr) A, B, C, D, E, F, G, H, I, J, A, B, C, D, E, F, G, H, ...
## $ N      (int) 55, 53, 55, 53, 58, 55, 58, 59, 56, 56, 51, 49, 53, 52,...
## $ mean.Y (dbl) 151.3570, 151.5632, 152.2151, 153.0880, 153.2191, 153.3...
## $ sd.Y   (dbl) 2.938866, 3.068675, 3.202740, 2.654993, 3.071673, 3.102...
## $ Age    (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
## $ X      (int) 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0

それぞれの変数は次の通り

図１の真似

給食タイプの背景には調査した県が存在しているが、こうして見てみると確かに調査した県によって平均身長が影響を受けた、という印象がある。

本と同じく、まずは一般化線形モデルでの推定結果を見てみる。これは次のような線形予測子からなるモデルのRでのglm()関数を使った推定結果。

\(\beta_1 + \beta_2 A_i + \beta_3 A_i X_j\)

d %>% glm(mean.Y ~ 1 + Age + Age:X, data = .) %>% 
  broom::tidy() %>% 
  knitr::kable(format = "markdown", digits = 4)
term estimate std.error statistic p.value
(Intercept) 152.9048 0.3979 384.2368 0.0000
Age 5.6473 0.6893 8.1932 0.0000
Age:X -1.7196 0.7959 -2.1606 0.0453

本の内容通りだが、パラメータ\(\beta_3\)を示すXの効果は、「乙タイプの給食を食べている学校のほうが乙タイプの学校よりもおよそ1.72低い」というような結果を導く。

🔨 RStanで階層ベイズ化

というわけで、ここから給食タイプの背景にある県の効果、給食グループごとの身長増加の効果を考慮した階層ベイズモデルを構築し…というのが岩波DSの内容である。この階層ベイズモデルのStanのコードは@berobero11さんが書かれているものを写経して使わせてもらう（stan_codeというオブジェクト名で保存）。

Rでは、データをlistクラスへ格納し、{rstan}パッケージのstan()でMCMCのサンプリングを実行する際の設定を行う。Stanコードの書き方は松浦さんの「Stan入門」に詳しく書かれている。stanコードの書き方のコツや注意点についても書かれており、自分のような初学者にはとても参考になる。

# データをリストに格納。
# Stanコードのdataブロックで宣言した変数。データの数が一致するように注意する
d_list <- list(N_r = 2,
               N_pref = length(unique(d$pref)), # データにおける都道府県の数
               Mean_Y = t(matrix(d$mean.Y, length(unique(d$pref)), 2)),
               Sigma_Y = t(matrix(d$sd.Y/sqrt(d$N), length(unique(d$pref)), 2)), # dim() 2 10... N_r  N_pref
               X = t(matrix(d$X, length(unique(d$pref)), 2)),
               Age = t(matrix(d$Age, length(unique(d$pref)), 2)))

# stanコードのコンパイル
#   少ないサンプル数などで試し射ちをしたあとなど、
#   Stanコードに修正を加えない段階になれば実行するのが良い？
stan_code_c <- stan_model(model_code = stan_code)

# MCMCサンプリングの実行
fit_sample <- sampling(object = stan_code_c, 
                       data = d_list,
                       iter = 30000, 
                       chain = 3,
                       warmup = 1000, 
                       seed = 71,
                       sample_file = "iwanami_ds/kubo_model2.csv",
                         # kubo_model2_X.csv という名称の
                         # mcmcサンプリング実行結果を出力したファイルが保存される
                         # （Xの部分には1からchain引数に指定した数が入る）
                       diagnostic_file = "iwanami_ds/diagnostic.txt")

🔬 推定結果の確認と収束判定

sampe_file引数で実行結果を保存しておくと（それに加えseed引数で乱数の固定をすると再現性のある結果になる）、一度Rセッションを終了してからサンプリング結果をよびだせる。その際は次のように、sprintf()関数などを使ってcsvファイルの一式をread_stan_csv()関数に渡すと良い。

(csvfiles <- sprintf("iwanami_ds/kubo_model2_%s.csv", 1:3))

[1] “iwanami_ds/kubo_model2_1.csv” “iwanami_ds/kubo_model2_2.csv” [3] “iwanami_ds/kubo_model2_3.csv”

fit_sample <- read_stan_csv(csvfiles = csvfiles)
fit_sample %>% summary() %>% .$summary %>% head(3) %>% knitr::kable(digits = 4)
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta[1] 152.9105 0.0038 0.4313 152.0460 152.6480 152.9120 153.1730 153.7760 12751.02 1.0001
beta[2] 5.1791 0.0041 0.5795 4.0189 4.8301 5.1798 5.5268 6.3402 19946.56 1.0002
beta[3] -0.8148 0.0058 0.8159 -2.4585 -1.3080 -0.8153 -0.3167 0.8044 20064.20 1.0001

Stanコードで指定したパラメータについての推定値の要約表（EAP推定値 mean、収束判定指標 Rhatなど）を見てみと今回求めたすべてのパラメータでRhatが1に近い値となっている。これはサンプリングされた値が定常分布へ収束していることを示唆している（Rhatは1.2よりも小さいことが収束の判定基準として採用されることが多く（Gelman-Rubinの統計量）、岩波DSサポートページでの松浦さん曰く、Chain数３以上で実行した際にRhatが1.1以下であればだいたい問題ないという）。

Rhat以外にも、MCMCサンプリングの収束を確認する方法が幾つか存在する。その中でも、視覚的に収束を確認できる方法として以下のものがある。traceplot()関数の引数parsには、パラメータを指定することができる（しないと多すぎてアレなことになりかねない）。

「trace plotは見ろ」… 松浦さんもそう言っている。

traceplot(fit_sample, pars = "beta")

実行した３回の連鎖の軌跡が、ほぼ同一線状に並んでいることからも、先ほどのMCMC計算は収束している、という結果で良いだろう。

{rstan}でのMCMC実行結果はstanfitクラスとして保存されるオブジェクトであり、これはS4クラスに属している。そのため要素へのアクセスには@スロットを利用する。そもそものスロット名はstr()なり、slotName()で確認できる。

fit_sample %>% {
  print(class(.))
  print(slotNames(.))
  print(.@model_pars)
  str(., max.level = 2)
}
## [1] "stanfit"
## attr(,"package")
## [1] "rstan"
##  [1] "model_name" "model_pars" "par_dims"   "mode"       "sim"       
##  [6] "inits"      "stan_args"  "stanmodel"  "date"       ".MISC"     
## [1] "beta" "s_r"  "r"    "mu"   "lp__"
## Formal class 'stanfit' [package "rstan"] with 10 slots
##   ..@ model_name: chr "kubo_model2"
##   ..@ model_pars: chr [1:5] "beta" "s_r" "r" "mu" ...
##   ..@ par_dims  :List of 5
##   ..@ mode      : int 0
##   ..@ sim       :List of 12
##   ..@ inits     : list()
##   ..@ stan_args :List of 3
##   ..@ stanmodel :Formal class 'stanmodel' [package "rstan"] with 4 slots
##   ..@ date      : chr "Sat Nov 07 18:49:47 2015"
##   ..@ .MISC     :<environment: 0x7fb791983318>

💹 {ggmcmc}パッケージを使った推定値の可視化

{ggmcmc}はMCMC計算を実行した後のパラメータの扱いや、各種の図を描く際に便利である。しかしパラメータ数が多いと扱いが面倒なので注意する（ってberobero11さんが言っていた）。そのため、必要なパラメータ名だけを抽出しておくと良い（データフレームオブジェクトなので操作が楽）。

df_param <- ggs(fit_sample)
df_param %>% {
  print(class(.))
  dplyr::tbl_df(.)
}
## [1] "tbl_df"     "tbl"        "data.frame"
## Source: local data frame [3,915,000 x 4]
## 
##    Iteration Chain Parameter   value
##        (int) (int)    (fctr)   (dbl)
## 1          1     1    beta.1 153.544
## 2          2     1    beta.1 153.644
## 3          3     1    beta.1 153.772
## 4          4     1    beta.1 153.506
## 5          5     1    beta.1 153.148
## 6          6     1    beta.1 153.269
## 7          7     1    beta.1 153.517
## 8          8     1    beta.1 153.185
## 9          9     1    beta.1 152.449
## 10        10     1    beta.1 152.795
## ..       ...   ...       ...     ...
# betaとつくパラメータだけを抽出
fit_param_beta1 <- df_param %>% dplyr::filter(grepl("^beta", Parameter))

ggs_traceplot()はrstan::traceplot()と同様、MCMC計算における各連鎖の軌道をプロットする。

p1 <- ggs_histogram(fit_param_beta1) # 事後分布のヒストグラム
p2 <- ggs_autocorrelation(fit_param_beta1) # 事後分布の自己相関
p3 <- ggs_geweke(fit_param_beta1) # Gewekeの診断方法
p4 <- ggs_Rhat(fit_param_beta1) + xlab("R_hat") # 各パラメータのRhatのばらつき

# １真衣の図にまとめる
gridExtra::grid.arrange(p1, p2, p3, p4,
                        ncol = 2)

# ggmcmc()を実行すると、plot引数で指定した種類の図をまとめたpdfが作成される
# ggmcmc(fit_param_beta1, plot=c("density", "running", "caterpillar"))

また、ci()関数は、ggs()で処理したパラメータの信用区間をテーブル上に出力する。

ci(fit_param_beta1) %>% knitr::kable(digits = 4)
Parameter low Low median High high
beta.1 152.0460 152.2130 152.9120 153.6031 153.7760
beta.2 4.0189 4.2544 5.1798 6.1090 6.3402
beta.3 -2.4585 -2.1246 -0.8153 0.4930 0.8044 
# fit_sample %>% summary() %>% .$summary %>% head(3) %>% knitr::kable()
# を実行した時と同じような要約表を出力する

♨️ {rstan}によるサンプリングの実行関数

おまけ。

{rstan}ではMCMCサンプリングを実施するための類似した名称の関数が存在する。それぞれの違いが理解できていなかったので整理してみる。

外部に保存した.stan拡張子のファイルであればfile引数にファイルのパスを、Rのオブジェクトとしてstanコードを作成した場合はmodel_code引数にオブジェクト名を渡す。

関数名 コンパイル MCMCサンプリング 返り値 その他
stan() 必要（どちらでも良い） accept stanfitクラス StanコードのコンパイルとMCMCサンプリングを同時に行える
stan_model() （必要） しない stanmodelクラス StanコードのC++コードとしてのコンパイル。この関数だけでMCMCサンプリングを実行するわけではない
sampling() しない ok stanfitクラス stan_model()であらかじめコンパイルされたstanmodelクラスオブジェクトを引数objectに渡してサンプリングを実施

📚 参考

Enjyoto 😄