統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

RでWebスクレイピングしたい

背景

ちょっとした用事によりリコール情報について調査する機会がありました。これまでWebスクレイピングは経験がなかったのですが、便利なライブラリ({rvest})もあることだし、挑戦してみた結果を紹介します。 内容としては、国交省のサイトにある「リコール情報検索」(こちら)からリコールデータを取得し、テキストマイニングにかけた、というものです。

分析の進め方

分析の進め方は以下の通りです:

  1. サイトのページ構成を把握
  2. 構成にマッチするようにループを組んでrvest::read_htmlで順次読み込み
  3. 取得したテキストデータをMecab形態素解析
  4. 可視化

特別なことはしておらず、サイトのページ構成に合わせて必要なデータを取得し、可視化などを行います。

1.サイトのページ構成を把握

ここは、Rではなくブラウザの機能を使いました。例えばこの辺りの記事を参考に、Google Chromeデベロッパーツールでhtmlの構成を把握しました。

2.構成にマッチするようにループを組んでrvest::read_htmlで順次読み込み

ライブラリのインストール

ここからがRによる処理となります。まずは必要なライブラリをインストールして読み込みます。今回新しくインストールしたライブラリは以下の通りで、RMecabはこちらの記事を参考にMecabのインストールから行いました。以下は、Mecabのインストールが終わっている前提です。

install.packages("rvest")
install.packages("RMeCab", repos = "http://rmecab.jp/R")
install.packages("wordcloud")

ライブラリの読み込み

インストールしたライブラリ以外に、{dplyr}や{tidyr}などの定番ライブラリ、またテキストデータを扱うので{stringr}や{stringi}なども読み込んでいます。

library(rvest)
library(dplyr)
library(tidyr)
library(stringr)
library(stringi)
library(RMeCab)
library(ggplot2)
library(wordcloud)

{RMecab}のお試し

ここで少し{RMecab}を使ってみましょう。こんな使い方ができます。

res <- RMeCabC("すもももももももものうち")
> unlist (res)
    名詞     助詞     名詞     助詞     名詞     助詞     名詞 
"すもも"     "も"   "もも"     "も"   "もも"     "の"   "うち"

{rvest}のお試し

同じく{rvest}も試してみます。read_htmlで指定したURLのページ構成をごそっと取ってきてくれます。

source_url <- "http://carinf.mlit.go.jp/jidosha/carinf/ris/search.html?selCarTp=1&lstCarNo=1060&txtMdlNm=&txtFrDat=2000/01/01&txtToDat=2017/12/31&page=1"
recall_html <- read_html(source_url, encoding = "UTF-8")

取ってきたデータの中身を確認するためには、例えば以下のようにします:

> recall_html %>% 
+    html_nodes("body") %>% # HTMLのbodyタグの情報を取り出す
+    html_text() # テキストデータを取り出す
[1] "\n\n\n  \n  \n    \n    \n      \n      \n    \n    \n    \n    \n  \n  \n    トップページ>リコール情報検索>リコール届出情報一覧\n  \n  \n  \n    \n    \n    \n      リコール届出情報一覧\n      ご利用のブラウザはJavaScriptまたはCookieが無効に設定されています。設定を確認して再度アクセスしてください。\n      \n\n\n\n\n        \n          90件のデータがヒットしました\n        \n        \n          5 / 5 ページ\n        \n        番号\n              届出番号\n              届出日\n              通称名\n            81\n              リ 国-3988-0\n              2017/02/02\n              [リバティ]\n            82\n              リ 国-3988-0\n              2017/02/02\n              [ブルーバードシルフィ]\n            83\n              リ 国-3988-0\n              2017/02/02\n              [キャラバン]\n            84\n              改 国-0513-0\n              2017/01/27\n              [デイズ]\n            85\n              改 国-0513-0\n              2017/01/27\n              [デイズ ルークス]\n            86\n              リ 国-3944-1\n              2017/01/27\n              [デイズ]\n            87\n              リ 国-3944-1\n              2017/01/27\n              [デイズ ルークス]\n            88\n              リ 国-3944-2\n              2017/01/27\n              [デイズ]\n            89\n              リ 国-3944-2\n              2017/01/27\n              [デイズ ルークス]\n            90\n              リ 国-3977-0\n              2017/01/27\n              [セレナ]\n            \n        \n\n          \n            \n          \n        \n\n\t \n\n        \n      \n\n\n\n      \n        \n        \n      \n    \n    \n    \n    \n      トップページ\n        自動車のリコール制度について\n        リコール情報検索\n        リコール届出情報一覧\n        自動車不具合情報ホットライン\n        不具合情報検索\n        事故・火災情報検索\n        よくあるお問い合わせ\n        公表資料\n        自動車を安全に使うためには\n        利用規約等\n        バナーダウンロード\n      \n        \n      \n  \n  \n\n\n  Copyright © 2001-myDate = new Date() ;myYear = myDate.getFullYear();document.write(myYear); Ministry of Land, Infrastructure, Transport and Tourism All rights reserved.\n\n\n\n\n(function(){\n    var p = ((\"https:\" == document.location.protocol) ? \"https://\" : \"http://\"), r=Math.round(Math.random() * 10000000), rf = window.top.location.href, prf = window.top.document.referrer;\n    document.write(unescape('%3C')+'img src=\"'+ p + 'acq-3pas.admatrix.jp/if/5/01/7b9f970b303989123e334299f50a384a.fs?cb=' + encodeURIComponent(r) + '&rf=' + encodeURIComponent(rf) +'&prf=' + encodeURIComponent(prf) + '\" alt=\"\"  width=\"1\" height=\"1\" '+unescape('%2F%3E'));\n})();\n\n\n\n\nnew Image(1, 1).src=\"//data-dsp.ad-m.asia/dsp/api/mark/?m=323gb&c=bcMR&cb=\" + Math.floor(new Date().getTime() / 86400);\n\n  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){\n  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),\n  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)\n  })(window,document,'script','//www.google-analytics.com/analytics.js','ga');\n\n  ga('create', 'UA-52116336-5',{'allowLinker': true});\n  ga('require','linker');\n  ga('linker:autoLink',['destination']);\n  ga('require','displayfeatures');\n  ga('send', 'pageview');\n\n"

またページの全てのtableのテキストを取り出す時はこのようになります:

> recall_html %>%
+    html_nodes(xpath="//table") %>%
+    html_text()
[1] "番号\n              届出番号\n              届出日\n              通称名\n            81\n              リ 国-3988-0\n              2017/02/02\n              [リバティ]\n            82\n              リ 国-3988-0\n              2017/02/02\n              [ブルーバードシルフィ]\n            83\n              リ 国-3988-0\n              2017/02/02\n              [キャラバン]\n            84\n              改 国-0513-0\n              2017/01/27\n              [デイズ]\n            85\n              改 国-0513-0\n              2017/01/27\n              [デイズ ルークス]\n            86\n              リ 国-3944-1\n              2017/01/27\n              [デイズ]\n            87\n              リ 国-3944-1\n              2017/01/27\n              [デイズ ルークス]\n            88\n              リ 国-3944-2\n              2017/01/27\n              [デイズ]\n            89\n              リ 国-3944-2\n              2017/01/27\n              [デイズ ルークス]\n            90\n              リ 国-3977-0\n              2017/01/27\n              [セレナ]\n            "

本番

それではここからが本番です。まずは対象となるページと、したいことが何であるかを確認しておきましょう。

f:id:ushi-goroshi:20180222163557p:plain

こちらが今回の分析対象となるページです。この検索条件として例えば車名を「ニッサン」、届出日を「2017/01/01」〜「2017/12/31」としてみましょう。

f:id:ushi-goroshi:20180222163719p:plain

このように条件に合致したリコール情報を一覧表示してくれます。例えば1個目をクリックすると、

f:id:ushi-goroshi:20180222163816p:plain

リコール情報の詳細について知ることができます。このうち、上段の表にある「車名/メーカー名」や「不具合装置」、「対象台数」などを取得したいのですが、リンクを一つずつ辿ってコピーしてくるのは大変なので、スクリプトを書いて情報を取ってきたい、というのが今回の取り組みです。

スクリプトの大まかな内容としては、

  1. 検索結果の画面から、各リコールの詳細結果画面へのURLを取得する
  2. 取得したURLに順次アクセスし、必要な情報を取り出してまとめる
  3. 次のページに移動し、繰り返し

という感じになります。

3についてですが、幸いなことに1の検索結果URLは、次ページを確認すると末尾が「page=2」となっています。ここから元のページに戻ると「page=1」となっており、数値を変更するだけで任意のページに行けそうなので、検索結果のページ数(今回は5)だけメモしておけばループで回せそうです。 また、各リコールの詳細結果画面についてはURLが「http://carinf.mlit.go.jp/jidosha/carinf/ris/detail/1141591.html」 のようになっており、末尾の「数字7桁」を変えていけば良さそうです。

というわけで、以下のように検索結果と各リコールの詳細結果画面のURLについて、変更がない部分を定義しておきます。

src_url  <- "http://carinf.mlit.go.jp/jidosha/carinf/ris/search.html?selCarTp=1&lstCarNo=1060&txtMdlNm=&txtFrDat=2017/01/01&txtToDat=2017/12/31&page="
link_url <- "http://carinf.mlit.go.jp/jidosha/carinf/ris/"

また分析に用いる項目を以下の5つとし、結果の格納用のデータフレームを準備しておきます。

target_column <- c("車名/メーカー名", "不具合装置", "状 況", "リコール開始日", "対象台数")
html_tbl_all  <- data_frame()

以下、リコール情報を順次取得していきます。スクリプトの流れのセクションで書いたように、検索結果の画面のページを変えつつ、各リコールの詳細結果画面へのURLを取得し、read_htmlでデータを取り出していきます。 下のスクリプトtarget_url_listsapplyすればもっと早くなるかもしれませんが、今回はパフォーマンスを求めたい訳ではないので素直にLoopを回しました。

st <- Sys.time()
## iはページ数。事前にメモしておく。今回は5
for (i in 1:5) {

   ## 検索結果の各ページのURLを指定し、データを取得
   target_page <- paste0(src_url, i)
   recall_html <- read_html(target_page, encoding = "UTF-8")

   ## 検索結果画面から、各リコール詳細結果へのURLを取得
   target_url_list <- 
      recall_html %>% 
      html_nodes(xpath = "//a") %>% # aタグに格納されている
      html_attr("href") %>% # href属性のデータを取り出す
      as_data_frame() %>% 
      filter(grepl("detail", .$value)) # 詳細結果は"detail" + 数字7桁 + .htmlで構成されている

   ## 詳細結果の数
   l <- nrow(target_url_list)

   ## ここから各詳細結果へアクセスし、データを取得する   
   for (j in 1:l) {
      ## アクセス負荷を軽減するため、少し間を置く
      Sys.sleep(2) 
      
      ## 詳細結果へのURLを指定し、データを取得
      target_url      <- paste0(link_url, target_url_list$value[j])
      recall_html_tmp <- read_html(target_url)
      html_tbl_tmp    <- html_table(recall_html_tmp)[[1]] ## 上段のテーブルのデータを取得
      
      ## 4列あるが、1・3列目に項目名が、2・4列目にデータが入っているので、2列のデータに直す
      html_tbl <- 
         html_tbl_tmp %>% 
         filter(X1 %in% target_column) %>% ## 必要な情報を抽出
         rename("Term" = X1, "Value" = X2) %>%
         select(Term, Value) %>% 
         bind_rows(
            html_tbl_tmp %>% 
            filter(X3 %in% target_column) %>% 
            rename("Term" = X3, "Value" = X4) %>%
            select(Term, Value)) %>% 
         spread(Term, Value) ## 順次追加していけるよう、wideに変換

      ## データを追加      
      html_tbl_all <- bind_rows(html_tbl_all, html_tbl)
   }
}
> Sys.time() - st
Time difference of 9.495965 mins

このスクリプトでは一年分の日産のデータを取得するのに、私の環境で約10分かかりました。結構時間がかかるので、データを保存しておきます。

save(html_tbl_all, file = "Recall_Data.Rdata")

3.取得したテキストデータをMecab形態素解析

ではこれ以降、取得したデータで分析を行います。と言ってもMecabによる形態素解析を掛けた後は集計して可視化するぐらいのものです。その前にデータを確認してみましょう。検索結果では90件と表示されていましたが、ちゃんと取れているでしょうか。

> dim(html_tbl_all)
[1] 90  5

大丈夫なようですね。データも見てみましょう。

> head(html_tbl_all)
# A tibble: 6 x 6
  リコール開始日 `車名/メーカー名` `状 況`                                        対象台数 不具合装置    Num
  <chr>          <chr>             <chr>                                           <chr>    <chr>       <dbl>
1 20171214日 いすゞ            ① 電源電圧が12Vのテールゲートリフタ装着車の後部反射器において、選定が不適切なため、反射… 1,153台  その他(保安灯火)1153
2 20171214日 いすゞ            ② 電源電圧が24Vのテールゲートリフタ装着車の後部反射器及び後退灯において、選定が不適切な… 162台    後退灯        162
3 20171215日 ニッサン          電源分配器の基板において、回路基板の製造時に不要な半田が付着した状態で防湿材がコーティングさ… 316,759… その他(電気装置)316759
4 20171215日 ニッサン          電源分配器の基板において、回路基板の製造時に不要な半田が付着した状態で防湿材がコーティングさ… 316,759… その他(電気装置)316759
5 20171215日 ニッサン          電源分配器の基板において、回路基板の製造時に不要な半田が付着した状態で防湿材がコーティングさ… 316,759… その他(電気装置)316759
6 20171201日 いすゞ            小型トラックの燃料噴射装置において、サプライポンプをエンジンに締結する取付けボルトの締付トル… 83,591台 エンジン一般…  83591

「車名」を確認すると一部にニッサン以外が含まれていますね。しかし当該の詳細結果を確認すると「いすゞ」とともに「ニッサン」がリコール対象となっており、間違いではないようです。

また「状況」を確認すると、全く同じ文言が同じ日付で出ています。これは、このリコールの届出が車種別に行われているためであり、例えば「2017年12月15日」の例では、「セレナ」「キューブ」「バネット」で同じ理由によりリコールの届出があったようです。本来この後の形態素解析では、これらのテキストを集約するべきでしょう(今回はお試しなのでやりませんが)。

分析対象となる部分を取り出す

さて、今回形態素解析の対象としたいテキストデータは「状 況」です。RMecabではテキストファイルからデータを読み込んで処理するので、テキストとして書き出しておきましょう。

load("Recall_Data.Rdata") ## 必要なら
txt_defect_situation <- html_tbl_all$`状 況`
write.csv(txt_defect_situation, file = "Situation.csv")

読み込み

書き出したテキストファイルを以下のように読み込みます。

txt_situ <- RMeCabFreq("Situation.csv")

全部で529個の単語に分割されたようです。

データ加工

テキストデータはMecabによって形態素解析され、単語ごとに分割された上で品詞を割り当てられています。このうち、単語抽出の対象となりそうなものだけを使用します。今回は名詞の頻度を確認します。

Noun_res_situ <- 
   txt_situ %>% 
   filter(Info1 == "名詞") %>% 
   filter(!Info2 %in% c("非自立", "代名詞")) %>%
   group_by(Term, Info1) %>% 
   summarise("TF" = sum(Freq)) %>% 
   ungroup() %>% 
   arrange(desc(TF)) %>% 
   mutate(Pos = factor(Term, levels = .$Term))

上のスクリプトで最後にfactorにしているのは、グラフにする時に単語の並び順ではなく頻度で表示するためです。

4.可視化

では加工済みのデータを用いて単語の頻度を可視化してみましょう。まずは棒グラフですが、単語の数が多いので上位20個に限定しています。なおMacの場合、日本語が表示されない可能性があります(私は表示されませんでした)。その場合、下記のページが参考になると思います。私は下記を全て実行したところ表示できるようになりました:

ggplot(Noun_res_situ[1:20, ], aes(x = Pos, y = TF)) +
   geom_bar(stat = "Identity") +
   theme_bw(base_family = "HiraKakuProN-W3")

f:id:ushi-goroshi:20180222164831p:plain

「"」や「","」のような変な単語(?)が混ざっていますね。これは流石に格好悪いので除いておきます。

Noun_res_situ %>% 
   filter(!Term %in% c("\"", "\",\"")) %>% 
   slice(1:20) %>% 
   ggplot(., aes(x = Pos, y = TF)) +
   geom_bar(stat = "Identity") +
   theme_classic(base_family = "HiraKakuProN-W3")

f:id:ushi-goroshi:20180222164850p:plain

なるほど、なんかそれっぽい単語が抽出されていますね〜。しかし、実際のテキストを見ていると、例えば「検査」という単語は「完成検査」という熟語として使われることが多いなど、ドメイン特有の表現があったりします。そういった特有の表現を集めた辞書がないと、こういった単語抽出はあまり効果的でなかったりします。

続いてWordcloudを作成してみます。ここでも単語の数が多いので絞ろうと思うのですが、個数ではなく出現頻度でfilterしましょう。TFが4以上の単語を抽出すると、以下のようになります。

Noun_res_situ_4 <- 
   Noun_res_situ %>% 
   filter(!Term %in% c("\"", "\",\"")) %>% 
   filter(TF >= 4)
par(family = "HiraKakuProN-W3")
wordcloud(Noun_res_situ_4$Term, Noun_res_situ_4$TF, random.color = TRUE, colors = rainbow(10))

f:id:ushi-goroshi:20180222164906p:plain

5.おまけ

以上でやりたかったことは終わりなのですが、最後におまけでリコール総台数の確認をしてみます。「対象台数」列に入っている数値を集計したいのですが、文字列として入力されているので修正します。

html_tbl_all$Num <- html_tbl_all$対象台数
html_tbl_all$Num <- str_replace_all(html_tbl_all$Num, ",", "")
html_tbl_all$Num <- str_replace_all(html_tbl_all$Num, "台", "")
html_tbl_all$Num <- as.numeric(html_tbl_all$Num)
ggplot(html_tbl_all, aes(x = Num)) +
   geom_histogram() +
   theme_classic()

f:id:ushi-goroshi:20180222164919p:plain

最後に

というわけで今回は{rvest}を用いたWebスクレイピングに挑戦しました。read_htmlで簡単にWebサイトのデータを取得することができ、html_text()html_table()で簡単に加工することが可能なため、初めての挑戦ではあったものの大きな引っかかりもなく進めることができました。今まで何となく敬遠していたのですが、積極的に使っていきたい技術ですね。

階層ベイズと状態空間モデルで広告効果を推定したい

背景

これまでMarketing Mix ModelingMMM)におけるAdStock効果の推定について色々と記事を書いてきましたが、その他にも試したいと思っているモデルがいくつかあります。その一つが階層ベイズモデルと状態空間モデルを同時に取り扱うものです。

例えば「地域別の売上推移のデータ」が手元にあると考えてみましょう。地域ではなく人や商品でも構いませんが、ある要因の各水準がそれぞれ時系列データを持っている状況(いわゆるパネルデータ)で、ひとまずここでは地域とします。このようなデータはあらゆる会社で保有していることでしょう。 今、各地域についてMMMにより広告効果を推定することを考えたとき、どのようなモデリングが可能でしょうか?

シンプルに考えれば、地域ごとに一つずつモデルを作るという方法が挙げられます。例えば地域の数が2つ3つしかなかったり、モデルの作成に時間をかけることが可能であればこの方法は有効かもしれません。しかし問題もあります。地域ごとに独立してモデルを作成するということは、各地域のパラメータは互いに独立であるとの仮定を置くことになるのですが、それは事実でしょうか?

確かにTVや新聞には地方欄やローカル局があり、地域限定のコンテンツや提供される場合もあるでしょう。しかしどちらかと言えば全国で同じコンテンツが提供される割合の方が大きいでしょうし、購買行動に地域による差がそれほど認められるとは、経験的にも思えません。もし仮に地域による異質性というものがそれほど強くないのであれば、むしろそういった情報を積極的に活かした方が良いのではないでしょうか?

そのような時に、パラメータ間に緩やかな縛りをかける方法として階層ベイズモデルという手法があります。詳しくは書籍1を参考にして頂くとして、ここでは「各地域の広告効果を表すパラメータ \betaを生成する分布を仮定する」方法と説明しておきます。この分布の幅(バラつき、分散)が十分に大きければ広告効果は地域によって様々な値を取り得ますし、逆に分布の幅が極端に狭ければ地域ごとの広告効果に差はないことを意味します。

目的

この階層構造を持つモデルを、時系列データの分析でお馴染みの状態空間モデルと合わせて取り扱いたいというのが今回の試みです。具体的には、以下のようなモデルを想定しています。

 {
\begin{equation}
y_{A,t} = \mu_{A,t} + \beta_{A}X_{A,t} + \epsilon_{A,t} \
\end{equation}
}

 {
\begin{equation}
\mu_{A,t} \sim N(\mu_{A,t-1}, \sigma^{2}_{\mu}) \
\end{equation}
}

 {
\begin{equation}
\epsilon_{A,t} \sim N(0, \sigma^{2}_{\epsilon}) \
\end{equation}
}

 {
\begin{equation}
\beta_{A} \sim N(\beta_0, \sigma^{2}_{\beta})
\end{equation}
}

添字のAおよびtはそれぞれ地域(Area)と時点(time)を意味しており、 y_{A,t}は「ある地域At時点における観測値」となります。 ポイントは4行目で、地域ごとの広告効果 \beta_{A}に対して、その生成元となる分布(ここでは正規分布)を仮定しています。もし \sigma_{\beta}^{2}が大きければ \betaは様々な値を取り得ますが、値が小さければ非常に似通った値になることが理解できると思います。

ライブラリの読み込み

今回の分析では{rstan}を使用します。階層ベイズ + 状態空間モデルといった非常に複雑なモデルを表現できるのが強みです。

library(tidyverse)
library(ggplot2)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

シミュレーションデータの生成

分析用のデータを生成しますが、検証のポイントとしては以下のようになります:

  • 5地域(num_region)から4年分(num_year)の月次(num_month)データを収集し、合計240点の観測値を得たと想定
  • 回帰係数(広告効果)は平均0.5、分散0.01正規分布にしたがって発生
  • 状態変数のパラメータは各地域で共通
## シードの固定
set.seed(123)

## 観測値の数
num_region  <- 5
num_month   <- 12
num_year    <- 4
data_length <- num_month * num_year

## 回帰係数
mu_beta  <- 0.5
var_beta <- 0.01  
beta_ad  <- rnorm(num_region, mu_beta, sqrt(var_beta))
beta_ad_all <- rep(beta_ad, each = data_length)

## 分布の設定
# 状態変数の初期値
state_t0     <- 3
var_state_t0 <- 1

# 状態変数の分散
var_state <- 0.01

# 誤差項
var_error <- 0.01

## 説明変数X
scale_x  <- 5
zero_per <- 0.2
X <- ifelse(runif(num_region * data_length) > zero_per, 1, 0) * 
   rpois(num_region * data_length, scale_x)

上記の条件に従いデータを発生させます。状態変数を除けば先に示したモデルの通り簡単な線形モデルです。Area_IDは、後ほどStanでフィッティングを行うために必要となる変数で、YMは描画用です。

## 状態変数
State <- matrix(0, nrow = num_region * data_length)
for (i in 1:num_region) {
   for (j in 1:data_length) {
      if (j == 1) {
         ## 1, 49, 97, 145, 193行目が各地域の先頭となる
         State[(i-1) * data_length + j] <- rnorm(1, state_t0, sqrt(var_state_t0))
      } else {
         State[(i-1) * data_length + j] <- 
            State[(i-1) * data_length + j - 1] + rnorm(1, 0, sqrt(var_state))
      }
   }
}

## 目的変数
error <- rnorm(num_region * data_length, 0, sqrt(var_error))
Y     <- State + X * beta_ad_all + error

## 地域ごとのID
Area_ID <- 1:num_region
DF <- data.frame(
   "Area_ID" = rep(Area_ID, each = data_length), 
   "YM" = rep(1:data_length, time = num_region),
   "Y" = Y,
   "X" = X,
   "True_s" = State,
   "True_e" = error)

ここで観測値と状態変数が各地域でどのように推移しているか確認しておきます。

DF %>% 
   gather("Var", "Val", -c(Area_ID, YM)) %>% 
   filter(Var %in% c("Y", "True_s")) %>% 
   ggplot(., aes(x = YM, y = Val, color = Var)) +
   geom_line() +
   facet_wrap(~Area_ID)

f:id:ushi-goroshi:20180220113808p:plain

また、この時点で地域をまとめてlmで推定した$\beta$がどのようになるかも確認しておきましょう。

> coef(lm(Y ~ X, data = DF))

(Intercept)           X 
  3.8088184   0.5003114 

また地域ごとにデータを分割した場合も試してみます。

results <- list()
for (i in Area_ID) {
   results[[as.character(i)]] <- coef(lm(Y ~ X, data = DF, subset = Area_ID == i))
}
> print(cbind(do.call("rbind", results), beta_ad))
  (Intercept)         X   beta_ad
1    5.066906 0.4292258 0.4439524
2    2.700195 0.4837768 0.4769823
3    3.764447 0.6628493 0.6558708
4    2.953617 0.4988667 0.5070508
5    4.051693 0.5314325 0.5129288

そこそこ近い値が推定されていますね(汗

Stanによるフィッティング

気を取り直して、Stanを用いたフィッティングを実行してみます。Stanに渡すデータは以下の通りです。

dat_Stan <- list(N       = data_length,
                 K       = num_region,
                 Y       = DF$Y,
                 X       = DF$X,
                 Area_ID = DF$Area_ID)

またStanのスクリプトは以下のようになります。

data {
   int N; // 地域ごとの観測値の数(data_length)
   int K; // 地域の数(num_region)
   vector[N*K] Y; // 観測値のベクトル
   vector[N*K] X; // 説明変数のベクトル
   int<lower=1, upper=K> Area_ID[N*K];
}

parameters {
   real state_t0;     // 状態変数の0期目の平均
   vector[N*K] state; // 状態変数のベクトル
   real beta_0;       // 回帰係数の事前分布の平均
   vector[K] beta;    // 地域ごとの回帰係数のベクトル
   real<lower=0> var_state_t0; // 状態変数の0期目の分散
   real<lower=0> var_state;    // 状態変数の分散
   real<lower=0> var_beta_0;   // 回帰係数の事前分布の分散
   real<lower=0> var_error;    // 誤差分散
}

model {
   // 状態変数をサンプリング
   for (k in 1:K) {
      // 1期目の値は0期目の分布からサンプルする
      state[1 + (k-1)*N] ~ normal(state_t0, var_state_t0);
      
      // 2期目以降は前期の値を平均とした分布からサンプルする
      for(i in 2:N) {
         state[i + (k-1)*N] ~ normal(state[i-1 + (k-1)*N], var_state);
      }
   }
   
   // 回帰係数をサンプリング
   for (k in 1:K) {
      beta[k] ~ normal(beta_0, var_beta_0);
   }

   // Yをサンプリング
   for(i in 1:(N*K)) {
      Y[i] ~ normal(state[i] + beta[Area_ID[i]] * X[i], var_error);
   }
}

このスクリプトHB_SSM_Sim.stanとして保存し、フィッティングを行います。

fit_01 <- stan(file = '/YourDirectory/HB_SSM_Sim.stan',
               data = dat_Stan,
               iter = 1000,
               chains = 4,
               seed = 123)

結果の確認

フィッティングが終わったので、結果を見てみましょう。fit_01からサンプルを抽出して加工します。

## サンプルを抽出する
res_01 <- rstan::extract(fit_01)

## 該当するパラメータを取り出す
ests <- summary(fit_01)$summary
t0_rows    <- rownames(ests)[grep("state_t0", rownames(ests))]
state_rows <- rownames(ests)[grep("state\\[", rownames(ests))]
b0_rows    <- rownames(ests)[grep("beta_0", rownames(ests))]
beta_rows  <- rownames(ests)[grep("beta\\[", rownames(ests))]

## 状態変数
state_par <- 
   ests %>% 
   data.frame %>% 
   select(mean) %>% 
   mutate("Par" = rownames(ests)) %>% 
   filter(Par %in% state_rows) %>% 
   mutate("Area" = DF$Area_ID)

state_cmp <- data.frame(
   True = State,
   Est = state_par$mean,
   Area = as.factor(state_par$Area),
   YM = DF$YM
)

まずは状態変数の推定結果を見てみましょう。実際の値と推定値を並べてみます。 (なお以降の作業の前にサンプルのtraceプロットとヒストグラムを確認し、収束していると判断しています)

state_cmp %>% 
   gather("Var", "Val", -c(Area, YM)) %>% 
   ggplot(., aes(x = YM, y = Val, colour = Var)) +
   geom_line() +
   facet_wrap(~Area)

f:id:ushi-goroshi:20180220114352p:plain

おぉ、かなり精度良く推定できているようですね!状態変数の初期値や分散の推定値はどうでしょうか?

> ests %>% 
   data.frame %>% 
   select(mean) %>% 
   rename("Estimated" = mean) %>% 
   mutate("Par" = rownames(ests)) %>% 
   filter(Par %in% t0_rows) %>% 
   mutate("True" = c(state_t0, var_state_t0)) %>% 
   mutate("Simulated" = c(
      mean(DF[c(1, 49, 97, 145, 193), "True_s"]),
      var(DF[c(1, 49, 97, 145, 193), "True_s"])
   )) %>% 
   select(Par, True, Simulated, Estimated)

           Par True Simulated Estimated
1     state_t0    3 3.7429556  3.731479
2 var_state_t0    1 0.8521877  1.328192

状態変数初期値の平均および分散の真の値がそれぞれ3および1であったのに対し、生成されたデータでは3.7および0.85でした。推定された値は3.7および1.3で、分散がやや過大に推定されているようです。ただし以下に示すように、推定値の95%信用区間も結構広いため、逸脱しているとまでは言えないようです。

> ests %>% 
   data.frame %>% 
   select(X2.5., X97.5.) %>% 
   rename("Lower95" = X2.5.,
          "Upeer95" = X97.5.) %>% 
   mutate("Par" = rownames(ests)) %>% 
   filter(Par %in% t0_rows) %>% 
   select(Par, everything())

           Par   Lower95  Upeer95
1     state_t0 2.3975622 5.027106
2 var_state_t0 0.5549254 3.525860

続いて回帰係数はどうでしょうか。

beta_par <- 
   ests %>% 
   data.frame %>% 
   select(mean) %>% 
   mutate("Par" = rownames(ests)) %>% 
   filter(Par %in% beta_rows)

beta_cmp <- data.frame(
   True = beta_ad,
   Est = beta_par$mean
)

ggplot(beta_cmp, aes(x = True, y = Est)) +
   geom_point() +
   coord_fixed()

f:id:ushi-goroshi:20180220114731p:plain

こちらも、事前に設定した回帰係数と推定値が似通っているようです。数値を確認しても、良い精度で推定できています。

ests %>% 
   data.frame %>% 
   select(mean) %>% 
   mutate("Par" = rownames(ests)) %>% 
   filter(Par %in% beta_rows) %>% 
   bind_cols("True" = beta_ad) %>% 
   select(Par, True, mean)
      Par      True      mean
1 beta[1] 0.4439524 0.4520343
2 beta[2] 0.4769823 0.4825661
3 beta[3] 0.6558708 0.6604340
4 beta[4] 0.5070508 0.5001494
5 beta[5] 0.5129288 0.5355829

最後に、回帰係数の事前分布についても見ておきましょう。

> ests %>% 
   data.frame %>% 
   select(mean) %>% 
   mutate("Par" = rownames(ests)) %>% 
   filter(Par %in% b0_rows) %>% 
   bind_cols("True" = c(mu_beta, var_beta)) %>% 
   select(Par, True, mean)

         Par True      mean
1     beta_0 0.50 0.5184981
2 var_beta_0 0.01 0.1385978

回帰係数の事前分布の分散はちょっと大きく推定されているようですが、平均は近しい値となっているようです。

終わりに

以上、「階層ベイズと状態空間モデルを合わせて取り扱いたい」という試みでしたが、結果としては概ね満足の行くものになったと思います。序盤に書いた通り、階層ベイズ + 状態空間モデルのような非常に複雑なモデルであっても、stanを使えば非常に簡単に推定が可能です。

これまで広告効果の推定に関していくつか記事を書いてきましたが、実はゴールとなるモデルとしてはこれを想定していました。あとはこのモデルに、これまで紹介してきたようなAdStock効果の推定を組み込めば試してみたかったモデルは一通り完了することになります。

これらについても追って紹介したいと思います。


  1. 伊庭先生のベイズモデリングの世界、久保先生の緑本、松浦さんのアヒル本など良書はたくさんあります。

スタイン推定量を確かめてみる

いきなりですが、最近購入したベイズモデリングの世界を読んでいて非常に面白い話題を発見しました。

ベイズモデリングの世界

ベイズモデリングの世界

P118から始まるスタイン推定量についてなのですが、その直感的な意味について本文から引用すると:

それぞれの y_iを、全部の{ y_i}の平均値 \frac{1}{n}\Sigma_{i=1}^{n}y_iの方向に aだけ引っ張ってやる

ことで、 y_iそのものを推定量とするよりも、パラメータの真値 \hat{\theta_{i}}との二乗誤差の期待値が小さくなるそうです。これはいわゆる縮小推定量の一種で、引っ張りの程度 aをデータから適応的に求めるのがこの推定量のキモなのだそうです。

スタイン推定量 \hat{\theta_{i}^{S}}の定義自体はそれほど難しいものでもなく、以下のような式によって表されます:

 \displaystyle \hat{\theta_{i}^{S}} = (1-a)y_{i} + \frac{n}{a}\Sigma_{i=1}^{n}y_i

ここで、

 \displaystyle a = \frac{\sigma^{2}}{s^{2}}, \
s^{2} = \frac{1}{n-3}\Sigma_{i=1}^{n}(y_i - \frac{1}{n}\Sigma_{j=1}^{n}y_i)^{2}

です。この記事では、スタイン推定量が本当に y_iよりも二乗誤差が小さくなるのか、以下のように実験してみます。

実験内容

内容としては非常に簡単なもので、平均と分散が共に既知である正規分布からいくつかのサンプルを生成します。そのサンプルをそのまま用いた場合と、スタイン推定量を用いた場合とで、真のパラメータとの二乗誤差の期待値:

 \mathbb{E}[\Sigma_{i=1}^{n}(\theta_{i}^{*}({y_i}) - \theta_i)^{2}]

にどれほど差が生じるかを確認しています。

パラメータ設定

実験では、一回の試行につき平均が{1, 2, ..., n}、分散が1である正規分布に従うサンプルをn個生成することとしました。このn個のデータから各サンプルのスタイン推定量を求め、真値である{1, 2, ..., n}との二乗誤差を計算します。同時にデータそのものを用いた場合の二乗誤差も計算しておきます。上記の実験を1000回繰り返し、最終的にどちらが小さい値となるかを比較しました。

n   <- 10
mu  <- 1:n
sig <- 1

K   <- 1000
res <- matrix(0, K, 2)
実験開始

上記のパラメータに基づき、以下のように計算を行います。

set.seed(123)
for (i in 1:K) {
   ## サンプルを発生させる
   tmp <- rnorm(n, mu, sig)
   
   ## スタイン推定値を求める
   s2  <- (1/(n-3)) * sum((tmp - mean(tmp))^2)      
   a   <- sig / s2
   s_i <- (1-a)*tmp + (a/n)*sum(tmp)

   ## 結果を保存する
   res[i, 1] <- sum((s_i - mu)^2)/n # スタイン推定値を用いた二乗誤差の期待値
   res[i, 2] <- sum((tmp - mu)^2)/n # データそのものを用いた二条誤差の期待値
}
結果

それでは結果です。

> mean(res[, 1] / res[, 2])
[1] 0.9618857

> sum(res[, 1] > res[, 2])
[1] 310

スタイン推定量を用いた場合、データそのものの場合と比較して二乗誤差が96%と小さくなっていることが確認できます。

終わりに

今回の実験は以上で終了です。ではこれの何が面白いのかと言うと、スタイン推定量と階層ベイズモデルとの関連性についてなんですね。現在、別に書こうとしている記事で階層ベイズモデルに触れているのですが、その説明として

パラメータを生成する分布を仮定することで、パラメータに緩やかな縛りをかける

といった言い回しを考えています。ここでふと「階層ベイズでは(緩やかな縛りをかけることで)推定値にバイアスを与えているわけだが、それにも関わらずモデルとして良くなるというのは、一体どういう意味なのだろう?」という疑問が浮かんできました。バイアスを与えているわけなのだから、誤差が大きくなってモデルとしては劣化しそうな気がしませんか?

その時にこのスタイン推定量のことを思い出しました。つまり、階層ベイズではパラメータに対する分布を仮定することで推定値にバイアスを与えているが、スタイン推定量の性質を以って全体的には誤差を小さくすることができているということではないかと1

書籍にはスタイン推定量が誤差を減少させる証明も付いていますが、まだ読み込めていませんので上記の理解が間違っているかもしれません。しかしこの推定量自体も非常に面白く、紹介したいと思ったので実験してみた次第です。


  1. なお「ベイズモデリングの世界」ではこれをバイアス-バリアンスのトレードオフとして解説しています。

Bayesian GLMで過適合を避ける

長らく知人に貸していた「みんなのR」が手元に戻ってきたのでパラパラと見ていたら、見逃していた面白いトピックを見つけたので紹介します。内容としてはこちらの記事とほぼ同じで、基本的には書籍の該当部分を再現しているだけですが、{purrr}のmapを使った書き方にしてみたり、少しだけ変更を加えています。

書籍はこちら。P320の「Bayesian 縮小」からです。

みんなのR -データ分析と統計解析の新しい教科書-

みんなのR -データ分析と統計解析の新しい教科書-

分析環境の準備

ライブラリの読み込み

まずは必要なライブラリの準備から。以下のライブラリを使用します。

library(MASS)
library(coefplot)
library(tidyverse)

今回使うBayesian GLMのための関数bayesglmは{arm}というライブラリにあるのですが、{coefplot}とnamespaceが衝突するためarm::bayesglmとして実行することにします。

データのロード

続いて「みんなのR」の著者が公開しているサンプルデータをこちらのリンク先から取得し、適当な場所に保存します(urlを使ってloadしたのですが、うまくいきませんでした)。

保存したデータを以下のようにloadします。

load("適当なフォルダ/ideo.rdata") # ideoというデータフレームが読み込まれる

フィッティング

まずは、(ベイズではない)通常のモデリングをした場合の問題点を明らかにするためにglmでフィッティングしてみます。ここでideo$Yearごとに同じモデルをあてはめるのですが、折角なのでサブグループごとにフィッティングを行うための方法をいくつか試してみましょう。

1. glmのsubsetで分割する

まずはそのまま、glmのオプションsubsetでYearを指定してデータを分割する方法です。

theYears          <- unique(ideo$Year)
results_01        <- vector(mode = "list", length = length(theYears))
names(results_01) <- theYears
for (i in theYears) {
   results_01[[as.character(i)]] <- 
      glm(Vote ~ Race + Income + Gender + Education, 
          data = ideo, subset = Year == i, family = binomial(link = "logit"))
}
2. purrr::mapを用いる

続いて、同じことを{purrr}ライブラリのmap関数を使って実行します。

results_02 <- 
   ideo %>% 
   split(.$Year) %>% 
   map(~glm(Vote ~ Race + Income + Gender + Education, 
            data = ., family = binomial(link = "logit")))

mapを使えばかなりスッキリ書けますね。なおsplitは{base}の関数です。

3. tidyr::nestと組み合わせる

さらに続いて、{tidyr}のnestと組合わせてみましょう。まずはmapに渡す関数を定義しておきます。

追記:この書き方は書籍「Rではじめるデータサイエンス」の20章「purrrとbroomによる多数のモデル」を参考にしています。

year_model <- function(df) {
   glm(Vote ~ Race + Income + Gender + Education, 
       data = df, family = binomial(link = "logit"))
}

次にideo$Yearでgroup_byし、nest()入れ子にします。

year_data <- 
   ideo %>% 
   group_by(Year) %>% 
   nest()

nestすると新しいデータフレームのdataにYearごとのデータフレームが格納されることになり、慣れないとかなり戸惑いがあります。

> year_data
# A tibble: 14 x 2
    Year                 data
   <dbl>               <list>
 1  1948   <tibble [374 x 7]>
 2  1952 <tibble [1,217 x 7]>
 3  1956 <tibble [1,245 x 7]>
 4  1960   <tibble [882 x 7]>
 5  1964 <tibble [1,100 x 7]>
 6  1968   <tibble [869 x 7]>
 7  1972 <tibble [1,560 x 7]>
 8  1976 <tibble [1,292 x 7]>
 9  1980   <tibble [829 x 7]>
10  1984 <tibble [1,343 x 7]>
11  1988 <tibble [1,172 x 7]>
12  1992 <tibble [1,304 x 7]>
13  1996 <tibble [1,004 x 7]>
14  2000 <tibble [1,049 x 7]>

このデータに対して先ほど定義した関数を適用します。

results_03 <- map(year_data$data, year_model)
names(results_03) <- as.character(unique(ideo$Year))

この方法では結果のlistにYearがnamesとして渡らないため後からつけているのですが、データと名前があっているのか不安になってしまいますね。

結果の確認

では結果を見てみましょう。multiplotはlistの各要素にlmglmの結果を持つオブジェクトを渡し、関心のある変数の回帰係数を比較することができます。

いずれも同じ結果となっています。

multiplot(results_01, coefficients = "Raceblack", secret.weapon = TRUE) +
   coord_flip(xlim = c(-20, 10))
multiplot(results_02, coefficients = "Raceblack", secret.weapon = TRUE) +
   coord_flip(xlim = c(-20, 10))
multiplot(results_03, coefficients = "Raceblack", secret.weapon = TRUE) +
   coord_flip(xlim = c(-20, 10))

f:id:ushi-goroshi:20180115164702p:plain
(同じグラフなので二枚は省略)

さてここで結果を見てみると、1964年の係数とその標準誤差がおかしなことになっています。multiplotでplot = FALSEとし、数値を確認してみます。

> (multi_01 <- multiplot(results_01, coefficients = "Raceblack", plot = FALSE))
          Value Coefficient   HighInner     LowInner   HighOuter    LowOuter Model
1    0.07119541   Raceblack   0.6297813   -0.4873905   1.1883673   -1.045976  1948
2   -1.68490828   Raceblack  -1.3175506   -2.0522659  -0.9501930   -2.419624  1952
3   -0.89178359   Raceblack  -0.5857195   -1.1978476  -0.2796555   -1.503912  1956
4   -1.07674848   Raceblack  -0.7099648   -1.4435322  -0.3431811   -1.810316  1960
5  -16.85751152   Raceblack 382.1171424 -415.8321655 781.0917963 -814.806819  1964
6   -3.65505395   Raceblack  -3.0580572   -4.2520507  -2.4610605   -4.849047  1968
7   -2.68154861   Raceblack  -2.4175364   -2.9455608  -2.1535242   -3.209573  1972
8   -2.94158722   Raceblack  -2.4761518   -3.4070226  -2.0107164   -3.872458  1976
9   -3.03095296   Raceblack  -2.6276580   -3.4342479  -2.2243631   -3.837543  1980
10  -2.47703741   Raceblack  -2.1907106   -2.7633642  -1.9043839   -3.049691  1984
11  -2.79340230   Raceblack  -2.4509285   -3.1358761  -2.1084548   -3.478350  1988
12  -2.82977980   Raceblack  -2.4609290   -3.1986306  -2.0920782   -3.567481  1992
13  -4.45297040   Raceblack  -3.4433048   -5.4626360  -2.4336392   -6.472302  1996
14  -2.67827784   Raceblack  -2.2777557   -3.0788000  -1.8772336   -3.479322  2000

何と1964年の回帰係数は文字通り桁が違っており、誤差も100倍といったスケールで異なっています。

では、このような結果が得られたときにどうしたら良いでしょうか?その答えの一つが、他の年の結果などから「係数はこの辺りにあるだろう」といった事前知識をモデルに取り込むことです。弱情報事前分布をモデルに取り入れることは、ベイズの観点からは縮小と見做すことができ、Bayesian Shrinkageベイズ縮小)と呼ばれます。

Bayesian GLMによるフィッティング

では実際にやってみましょう。関数は{arm}ライブラリのbayesglmで、glmと同じように使えますが、事前分布を指定することが可能です。まずは先ほどのモデルに、回帰係数の事前分布として尺度パラメータ2.5のコーシー分布(=自由度1のt分布)を用いて当てはめてみます。

results_04 <- 
   ideo %>% 
   split(.$Year) %>% 
   # {arm}はlibrary(arm)とせず、arm::bayesglmとして直接呼び出す
   map(~arm::bayesglm(Vote ~ Race + Income + Gender + Education, 
                      data = ., family = binomial(link = "logit"),
                      prior.scale = 2.5, prior.df = 1))
結果の確認

1964年の異常値が解消されたことを確認してみましょう。

multiplot(results_04, coefficients = "Raceblack", secret.weapon = TRUE) +
   coord_flip()

f:id:ushi-goroshi:20180115165807p:plain

係数とその誤差がかなり改善されていますね!

事前分布の仮定が他の年の係数に影響を与えてはいないでしょうか?それも確認してみましょう。

multi_04 <- multiplot(results_04, coefficients = "Raceblack", plot = FALSE)
> cbind(sprintf("%.3f", multi_01$Value), sprintf("%.3f", multi_04$Value))
      [,1]      [,2]    
 [1,] "0.071"   "0.073" 
 [2,] "-1.685"  "-1.635"
 [3,] "-0.892"  "-0.867"
 [4,] "-1.077"  "-1.050"
 [5,] "-16.858" "-5.100"
 [6,] "-3.655"  "-3.524"
 [7,] "-2.682"  "-2.654"
 [8,] "-2.942"  "-2.864"
 [9,] "-3.031"  "-2.971"
[10,] "-2.477"  "-2.446"
[11,] "-2.793"  "-2.744"
[12,] "-2.830"  "-2.779"
[13,] "-4.453"  "-4.151"
[14,] "-2.678"  "-2.620"

1964年は5行目ですが、その他の年の係数に大きな変動はないようです。また異常値が解消されたことで係数のバラつきはかなり小さくなりました。

plot(cbind(multi_01$Value, multi_04$Value))

f:id:ushi-goroshi:20180115170220p:plain

> apply(cbind(multi_01$Value, multi_04$Value), 2, sd)
[1] 4.037096 1.335844

また異常値となっていた5行目を除外すると、相関係数はほぼ1となります。

> correlate(cbind(multi_01$Value[-5], multi_04$Value[-5])) %>% fashion()
  rowname   V1   V2
1      V1      1.00
2      V2 1.00  
plot(cbind(multi_01$Value[-5], multi_04$Value[-5]))

f:id:ushi-goroshi:20180115170528p:plain

おまけ

bayesglmのオプションにあるprior.scaleとprior.dfを変更することで様々な事前分布を仮定することができます。ここをどちらもInfとすると正規分布となり、glmの結果と一致するそうです。

ということで試してみたのですが、Yearごとのデータで当てはめるとエラーとなってしまったので、全体のデータに対して当てはめた結果を比較してみます。

results_05 <- glm(Vote ~ Race + Income + Gender + Education, 
                  data = ideo, family = binomial(link = "logit"))
results_06 <- arm::bayesglm(Vote ~ Race + Income + Gender + Education, 
                            data = ideo, family = binomial(link = "logit"),
                            prior.scale = Inf, prior.df = Inf)

# 小数点以下7桁ぐらいまでは一致する
> cbind(sprintf("%.7f", coef(results_05)), sprintf("%.7f", coef(results_06)))
      [,1]         [,2]        
 [1,] "0.1592985"  "0.1592984" 
 [2,] "-0.2833290" "-0.2833290"
 [3,] "-2.4806024" "-2.4806025"
 [4,] "-0.8534029" "-0.8534029"
 [5,] "-0.3919568" "-0.3919568"
 [6,] "-0.5816898" "-0.5816898"
 [7,] "-0.3790876" "-0.3790876"
 [8,] "-0.0601803" "-0.0601803"
 [9,] "0.0230734"  "0.0230734" 
[10,] "0.1551473"  "0.1551473" 
[11,] "0.7373030"  "0.7373030" 
[12,] "0.1753482"  "0.1753482" 
[13,] "0.0954255"  "0.0954255" 
[14,] "0.3491240"  "0.3491240" 
[15,] "-0.3277760" "-0.3277760"
[16,] "-0.1064615" "-0.1064615"
[17,] "0.1802366"  "0.1802366" 
[18,] "-0.1227574" "-0.1227574"
> correlate(cbind(coef(results_05), coef(results_06))) %>% fashion()
  rowname   V1   V2
1      V1      1.00
2      V2 1.00
plot(cbind(coef(results_05), coef(results_05)))

f:id:ushi-goroshi:20180115170829p:plain

ほとんど一致していますね!

終わりに

この記事では、bayesglm関数を使ったBayesian GLMの当てはめについて紹介しました。また事前分布を変更することでbayesglmの結果がglmと一致することを確認しました。

ベイズでは事前分布のバラつきの程度によってパラメータに対する分析者の「確信度」が表されます。この確信度とデータを観測した時の影響について、追ってブログで紹介したいと思っています。

ベイジアンネットワークで変数間の関連性を見る

実務でデータ分析を行うときにしばしば変数間の因果関係の有無およびその強さについて問われることがあるのですが、その回答としてベイジアンネットワークを使いたいと思っていたので調べてみました。
この記事ではベイジアンネットワークを扱うためのパッケージ ({bnlearn}および{BNSL})を用いた構造の推定およびそのTipsを示します。

分析環境の準備

ライブラリの読み込み

まずは必要なライブラリを読み込みます。今回はRにおけるベイジアンネットワークのメジャーなライブラリである{bnlearn}に加え、{BNSL}を用います。

library(bnlearn)
library(BNSL)
サンプルデータの作成

続いてこちらの記事を参考にサンプルデータを作成します。一部の係数を修正しましたが、記事同様にHeightとSBPがBMIに刺さり、BMIがFBSに刺さる構造となっています。

set.seed(123)
norm <- rnorm(4000)
Data <- matrix(norm, nrow = 1000, ncol = 4, byrow = T)
colnames(Data) <- c("Height", "BMI", "SBP", "FBS")

Data2 <- Data <- as.data.frame(Data)
Data2$Height <- 170 + Data$Height * 10
Data2$SBP    <- 120 + Data$SBP * 10
Data2$BMI    <- 22  + Data$Height * 5 + Data$SBP * 5
Data2$FBS    <- 90  + (Data2$BMI - 22) * 5 + Data$FBS * 10

以下のようなデータが得られます。

> head(Data2)
    Height      BMI      SBP       FBS
1 164.3952 26.99116 135.5871 115.66090
2 171.2929 24.95102 124.6092  92.10449
3 163.1315 24.68614 132.2408 107.02886
4 174.0077 21.22465 114.4416 103.99239
5 174.9785 27.99603 127.0136 115.25225
6 159.3218 11.53086 109.7400  30.36538

分析実行その1

構造学習

それではベイジアンネットワークによる関連性の推定に移ります。ベイジアンネットワークはデータから変数間の関連性を推定する構造学習と,推定された構造を元にパラメータを推定するフィッティングの二段階で実施されます。
今回、構造学習のアルゴリズムにはGrow-Shrinkを指定しました。Webで調べた限りだとHill-Climbingが使われる事も多い様子ですが、先ほど参照した記事によるとGrow-Shrinkグラフ理論による裏付けがあるとのこと。関数はgsです(Hill-Climbingはhc)。

str_01_gs <- gs(Data2)
str_01_hc <- hc(Data2)
par(mfrow = c(1, 2))
plot(str_01_gs)
plot(str_01_hc)

f:id:ushi-goroshi:20180111224629p:plain

gsでは設定した通りの構造が推定されましたが、hcではBMIとSBPの関係性が逆転しています。これ以降、基本的にgsを使うことにします。

フィッティング

ではこの構造を用いてパラメータを推定しましょう。なおここで得られるパラメータ推定値は個別にlmをかけたものと同じとなり、例えばBMIであればlm(BMI ~ Height + SBP)と等しい値となります。

fit_01 <- bn.fit(str_01_gs, Data2)
coef_01 <- coefficients(fit_01)

> coef_01
$Height
(Intercept) 
    170.091 

$BMI
(Intercept)      Height         SBP 
     -123.0         0.5         0.5 

$SBP
(Intercept) 
   120.1408 

$FBS
(Intercept)         BMI 
 -20.240828    5.019595 
> coef(lm(BMI ~ Height + SBP, Data2))
(Intercept)      Height         SBP 
     -123.0         0.5         0.5 

同じ値となっていますね。

分析実行その2

構造学習

次にサンプルデータのgaussian.testを用いてパラメータを推定します。元の構造はわかりませんが、gsによって推定した形は以下のようです。

str_02 <- gs(gaussian.test)
plot(str_02)

f:id:ushi-goroshi:20180111225340p:plain

フィッティング

この構造を用いてパラメータの推定を行おうとするとエラーが出ます。どうやらグラフの構造が悪く、フィッティングできないようです。

> fit_02 <- bn.fit(str_02, gaussian.test)
 bn.fit(str_02, gaussian.test) でエラー: 
  the graph is only partially directed.

そこでstr_02の中身を確認すると、modelが[partially directed graph]となっていることがわかります。またarcsのundirected arcsも1となっています。

> str_02

  Bayesian network learned via Constraint-based methods

  model:
    [partially directed graph]
  nodes:                                 7 
  arcs:                                  7 
    undirected arcs:                     1 
    directed arcs:                       6 
  average markov blanket size:           4.00 
  average neighbourhood size:            2.00 
  average branching factor:              0.86 

  learning algorithm:                    Grow-Shrink 
  conditional independence test:         Pearson's Correlation 
  alpha threshold:                       0.05 
  tests used in the learning procedure:  144 
  optimized:                             FALSE 

再度str_02の構造を確認するとノードBD間のエッジが無向になっていることがわかります。

f:id:ushi-goroshi:20180111225340p:plain


そこでこの構造にB -> Dのエッジを設定してみましょう。ベイジアンネットワークではこのように、試行錯誤の結果や分析者の知見を取り込みながら、データ全体の構造を推定していきます。
この修正によりstr_03のundirected arcsが0となり、modelがpartially directed graphではなくノード間の関連性を表したものとなっていることを確認しましょう。

> (str_03 <- set.arc(str_02, "B", "D"))

  Bayesian network learned via Constraint-based methods

  model:
   [A][B][E][G][C|A:B][D|B][F|A:D:E:G] 
  nodes:                                 7 
  arcs:                                  7 
    undirected arcs:                     0 
    directed arcs:                       7 
  average markov blanket size:           4.00 
  average neighbourhood size:            2.00 
  average branching factor:              1.00 

  learning algorithm:                    Grow-Shrink 
  conditional independence test:         Pearson's Correlation 
  alpha threshold:                       0.05 
  tests used in the learning procedure:  144 
  optimized:                             FALSE 

ちゃんとB -> D間の矢印が有向になっています。

> plot(str_03)

f:id:ushi-goroshi:20180111230018p:plain

では、この構造を用いて再度フィッティングしてみましょう。

> fit_03 <- bn.fit(str_03, gaussian.test)
> (coef_03 <- coefficients(fit_03))
$A
(Intercept) 
   1.007493 

$B
(Intercept) 
   2.039499 

$C
(Intercept)           A           B 
   2.001083    1.995901    1.999108 

$D
(Intercept)           B 
   5.995036    1.498395 

$E
(Intercept) 
   3.493906 

$F
 (Intercept)            A            D            E            G 
-0.006047321  1.994853041  1.005636909  1.002577002  1.494373265 

$G
(Intercept) 
   5.028076 

今度は上手くいったようですね!

分析実行その3

続いて直接効果と間接効果の分離を試みます。再度gsで構造を推定し、B -> Dのエッジと共にB -> Fのエッジも追加しましょう。これでBはFに対して直接刺さるエッジと、Dを経由してFに刺さるエッジを持つことになります。

str_04 <- set.arc(str_03, "B", "F")
par(mfrow = c(1, 3))
plot(str_02) # gsで推定した構造
plot(str_03) # B -> Dに修正
plot(str_04) # B -> Fを追加

f:id:ushi-goroshi:20180111230351p:plain

フィッティング

続いてフィッティングです。B -> Fのエッジの有無で係数がどのように変化するかを比較してみましょう。必要な部分だけ抜粋します。

fit_04 <- bn.fit(str_04, gaussian.test)

> coefficients(fit_04)
~ 省略 ~
$F
(Intercept)           A           B           D           E           G 
 -0.4123157   1.9947919  -0.1025991   1.0737532   1.0023721   1.4943286 
~ 省略 ~
> coefficients(fit_03)
~ 省略 ~
$F
 (Intercept)            A            D            E            G 
-0.006047321  1.994853041  1.005636909  1.002577002  1.494373265 
~ 省略 ~

切片に多少の変化があっただけで、係数としては大きな変化はないようです。どうやらBの変動による寄与はほとんどDを経由してFに影響しており、BのFへの直接の効果はなさそうです。

分析実行その4

構造学習

今度は構造の推定に{BNSL}を試してみましょう。

str_05 <- bnsl(gaussian.test)
plot(str_05)

f:id:ushi-goroshi:20180111231012p:plain

どうもgsを用いた時と比べて構造がかなり異なるようですね。比較してみましょう。

par(mfrow = c(1, 2))
plot(str_02)
plot(str_05)

f:id:ushi-goroshi:20180111231120p:plain

"A -> C"や"G -> F"など、方向性が逆転する関係が見られる上、Eも独立してしまいました。やはりデータのみから因果関係を推定するのは不安定ということになるのでしょうか。

フィッティング

この構造でフィッティングしてみます。フィッティング自体はbnlearn::bn.fitを用いています(結果は省略)。

fit_05 <- bn.fit(str_05, gaussian.test)

分析実行その5

さらに続いて、データを標準化した場合の結果も確認してみます。scaleでデータを標準化します。

my_gaussian <- as.data.frame(scale(gaussian.test))
str_06 <- bnsl(my_gaussian)
par(mfrow = c(1, 2))
plot(str_04)
plot(str_06)

f:id:ushi-goroshi:20180111231456p:plain

何と、標準化すると構造が変わってしまいました。。。フィッティングは省略します。

分析実行その6

最後に、ライブラリ{MASS}に格納されているサンプルデータbirthwtを用いて推定してみます。なおbnlearnではintegerは使えないので、各列をnumericに変換し、簡単のために変数も一部に限定します。

構造学習

最後なので、構造学習のアルゴリズムには色々なものを指定してみます。

vars <- c("low", "race", "smoke", "ht", "ui")
my_birthwt <- as.data.frame((do.call("cbind", lapply(birthwt[, vars],
                                                     as.numeric))))
str_07_01 <- hc(my_birthwt)
str_07_02 <- gs(my_birthwt)
str_07_03 <- iamb(my_birthwt)
str_07_04 <- fast.iamb(my_birthwt)
str_07_05 <- inter.iamb(my_birthwt)
str_07_06 <- tabu(my_birthwt)
str_07_07 <- mmhc(my_birthwt)
str_07_08 <- rsmax2(my_birthwt)
str_07_09 <- mmpc(my_birthwt)
str_07_10 <- si.hiton.pc(my_birthwt)
str_07_11 <- chow.liu(my_birthwt)
str_07_12 <- aracne(my_birthwt)
par(mfrow = c(1, 3))
plot(str_07_01)
plot(str_07_02)
plot(str_07_03)
par(mfrow = c(1, 3))
plot(str_07_04)
plot(str_07_05)
plot(str_07_06)
par(mfrow = c(1, 3))
plot(str_07_07)
plot(str_07_08)
plot(str_07_09)
par(mfrow = c(1, 3))
plot(str_07_10)
plot(str_07_11)
plot(str_07_12)

f:id:ushi-goroshi:20180111231826p:plain
f:id:ushi-goroshi:20180111231837p:plain
f:id:ushi-goroshi:20180111231930p:plain
f:id:ushi-goroshi:20180111231942p:plain

アルゴリズムによって少しずつ推定された構造が違うのですが、low -> smokeおよびrace -> smokeのエッジは概ね共通して見られるようです。しかし残念ながらこのデータは「出生児の低体重の有無」と「母体環境」であるため、エッジの向きは逆ですね。。。

終わりに

この記事ではベイジアンネットワークを扱うためのパッケージ およびその使い方を示しました。結果的にデータのみから因果関係を推定するのは困難であり、良い構造推定のためにはドメイン知識を持ったメンバーとの協業が必要となりそうですが、こういったデータの構造を確認しながらのディスカッションは非常に楽しそうですね。

WAICを計算してみる

ある統計モデルの予測精度を表す指標としてAICAkaike Information Criterion)というものがよく使われますが、stanを用いてモデリングした場合にはAICを計算することができません。しかし2009年に渡辺澄夫先生が発表したWAICWidely Applicable Information Criterion、広く使える情報量規準)はAICをより一般化したものとして定義され、MCMCによる事後分布から計算することができます。

現在のところ{rstan}では直接WAICを計算することができませんが、{loo}というパッケージを用いることでWAICを得ることができます。この記事では、stanの結果からWAICを計算により求め、それがlooの結果と一致することを確認します。

WAICの定義

まず始めにWAICの定義を確認しましょう*1。WAICは以下のように定義されます:

WAIC = -2(lppd - p_{WAIC})

ここで

\displaystyle lppd = \sum_{i = 1}^{N}\log{Pr(y_{i})}

であり、また

\displaystyle p_{WAIC} = \sum_{i = 1}^{N}V(y_{i})

です。Pr(y_{i})およびV(y_{i})は、それぞれ観測値iの事後分布から得られる対数尤度の平均および対数尤度の分散を表してます。なおlppdはlog-pointwise-predictive-densityです。

シミュレーション

簡単な数値例を発生させ、そのWAICを計算してみましょう。パラメータの推定には当然ながら{rstan}を利用します。

# シミュレーションデータの発生
set.seed(123)
N <- 100 # サンプルサイズ
b <- 1.2 # 回帰係数
X <- rnorm(N, 0, 1) # 説明変数
E <- rnorm(N, 0, 2) # 誤差項
Y <- b * X + E
D <- data.frame(Y, X) # データフレーム

回帰係数1.2、誤差標準偏差2として発生させたデータに対し、{rstan}を用いて線形回帰モデルを当てはめてみましょう。いつものように{rstan}をロードします。

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

stanに渡すデータを以下のように定義します。今回の例では単純な線形モデルを用いているので、渡すパラメータもほとんどありません。

dat_stan <- list(
   N = N,
   X = D$X,
   Y = D$Y
)

定義したdat_stanをstanに渡し、フィッティングします。

fit_01 <- stan(file = './WAIC.stan', 
               data = dat_stan, 
               iter = 3000,
               chains = 4,
               seed = 1234)

ここでWAIC.stanの内容は以下の通りです。簡単です。

data {
   int<lower=1> N;
   vector[N] X;
   vector[N] Y;
}

parameters {
   real b0;
   real b1;
   real<lower=0> sigma;
}

model {
   Y ~ normal(b0 + b1 * X, sigma);
}

フィッティングが終わると、以下のような結果が得られます:

> summary(fit_01)$summary
             mean     se_mean        sd         2.5%          25%
b0      -0.201859 0.002533264 0.1962258   -0.5791061   -0.3360031
b1       1.092886 0.002755139 0.2134122    0.6655672    0.9504853
sigma    1.964687 0.001920081 0.1421800    1.7080842    1.8637401
lp__  -116.169685 0.021008195 1.2014656 -119.2035323 -116.7161937
               50%           75%        97.5%    n_eff     Rhat
b0      -0.2018323   -0.06931008    0.1877817 6000.000 1.000041
b1       1.0919604    1.23311123    1.5082127 6000.000 1.000551
sigma    1.9570918    2.05945046    2.2612886 5483.250 1.000177
lp__  -115.8524219 -115.28245822 -114.7994731 3270.734 1.000839

色々と書いてあって情報が多いので、必要なところを抜き出しましょう。真の回帰係数および誤差標準偏差はそれぞれ1.2と2でした。

> summary(fit_01)$summary[c("b1", "sigma"), c("mean", "50%")]
          mean      50%
b1    1.092886 1.091960
sigma 1.964687 1.957092

回帰係数が少し低いようですが、まぁ良しとしましょう。

WAICの計算

それではこの結果を用いて、WAICを計算してみましょう。まずはlppdからですが、lppdは前述の通り以下によって計算されます:

\displaystyle lppd = \sum_{i = 1}^{N}\log{Pr(y_{i})}

まずはstanの結果からposterior sampleを得ます。stanを実行するときにwarmupを指定していないためsampleサイズ(ite)の半分が切り捨てられた結果、残る6,000サンプル(1,500 * 4 chain)がposterior sampleとして保存されています。

post_samples <- extract(fit_01)
> str(post_samples)
List of 4
 $ b0   : num [1:6000(1d)] -0.146 -0.13 -0.195 -0.535 -0.162 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ iterations: NULL
 $ b1   : num [1:6000(1d)] 1.036 0.928 1.004 1.309 1.005 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ iterations: NULL
 $ sigma: num [1:6000(1d)] 2.2 2.01 1.64 2.03 1.86 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ iterations: NULL
 $ lp__ : num [1:6000(1d)] -116 -115 -118 -117 -115 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ iterations: NULL

この6,000サンプルのそれぞれをパラメータとして与えた時の各 y_{i}の対数尤度を求めると:

Des <- cbind(1, X) # 計画行列(Design Matrix)
B   <- cbind(post_samples$b0, post_samples$b1) # パラメータ(Beta)
tmp <- matrix(NA, length(post_samples$b0), N) # 6000行、 100列の行列
for (i in 1:N) {
   tmp[, i] <- dnorm(Y[i], mean = B %*% Des[i, ], 
                     sd = post_samples$sigma)
}

となります。ここからlppdは:

lppd <- sum(log(colMeans(tmp)))

で得られます。
同じように p_{waic}を求めてみましょう。 p_{waic}の定義は以下のようでした:

\displaystyle p_{WAIC} = \sum_{i = 1}^{N}V(y_{i})

ここから、

pwaic <- sum(apply(tmp, 2, var))

として求められます。
最後に、これらの値を用いてWAICを求めます:

WAIC <- -2 * (lppd - pwaic)

これでようやくWAICを得ることがことができましたので、これらの値を出力してみると:

> cat("WAIC:", WAIC, "\n",
+     "lppd:", lppd, "\n",
+     "pwaic:", pwaic, "\n")
WAIC: 420.8022 
lppd: -207.2305 
pwaic: 3.170572 

となります。ちなみに同じモデルについて、lmを用いてフィッティングした時のAICを求めると:

fit_02 <- lm(Y ~ X, D)
> AIC(fit_02)
[1] 420.4523

であり、非常に近い値となります。

パッケージ{loo}を使ってWAICを求める。

では、今度は上記の行程で求めた数値をパッケージで得た値と比較してみましょう。WAICの計算には{loo}を用います。

library(loo)

looによりWAICを計算させるためには、stanのスクリプトにgenerated quantitiesブロックを追加しておく必要があります。以下をmodelブロックの後に追加します。

generated quantities {
   vector[N] log_lik;
   
   for (n in 1:N)
      log_lik[n] = normal_lpdf(Y[n] | b0 + b1 * X[n], sigma);
}

normal_lpdfというのはRで言うところのdnormですね。
このスクリプトをWAIC_02.stanとして保存し、stanを実行します。

fit_03 <- stan(file = './WAIC_02.stan', 
               data = dat_stan, 
               iter = 3000,
               chains = 4,
               seed = 1234)

結果は以下のようになります。前の結果と同じであることを確認してください。

> summary(fit_03)$summary[c("b1", "sigma"), c("mean", "50%")]
          mean      50%
b1    1.092886 1.091960
sigma 1.964687 1.957092

では、今度は{loo}のextract_log_likを使ってposterior sampleを抽出し、loo::waicによりWAICを計算します:

tmp02   <- extract_log_lik(fit_03)
WAIC_02 <- waic(tmp02)

> cat("WAIC by loo:", WAIC_02$waic, "\n",
+     "WAIC by manual:", WAIC, "\n")
WAIC by loo: 420.8022 
WAIC by manual: 420.8022 

ちゃんと同じ値が得られていますね!

終わりに

分析内容は以上となります。この記事ではWAICの定義を確認した上で、その計算のための方法を示し、既存のパッケージによる数値と一致することを確かめました。WAICは階層モデルのような複雑なものであっても適用可能で、stanを使っていく上では抑えておきたい指標なので今後も勉強を続けていきます。

*1:Statistical RethinkingのChapter 6から

Stanで疑似相関を見抜く

ベイズ統計をちゃんと勉強しようと思い、最近Statistical Rethinkingという本を読んでいます。

Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

その中でタイトルの通り面白い話があったので紹介がてら実行してみます。以下は5章1節からの引用ですが、x_spurとx_realはそれぞれyと相関の「ない」変数と相関の「ある」変数を表しています:

But when you include both x variables in a linear regression predicting y, the posterior mean for the association between y and x_spur will be close to zero, while the comparable mean for x_real will be closer to 1.

この記事の内容を端的に言えば「yと相関のあるx1と、x1と相関のあるx2(つまりyとは直接相関がない)が説明変数として与えられている時、stanを用いることで不要な変数(この場合x2)を除外できる」というものになります。
またこちらの本では全体を通して{rethinking}パッケージによってパラメータの推定が行われていますが、ここではせっかくなので{rstan}でやってみます。

前準備

まずは{rstan}を読み込みます。

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

続いてy、x1、x2をそれぞれ以下のように定義します。

## 乱数のシードの設定
set.seed(15)
N  <- 100 # サンプルサイズ
x1 <- rnorm(N) # yと相関のある(真の)説明変数
x2 <- rnorm(N, x_real) # x1と相関することでyと相関する(偽の)説明変数
y  <- rnorm(N, x_real) # x1と相関する目的変数

## データフレームにする
dat <- data.frame(y, x_real, x_spur)

定義からx2はyと相関しないため、本来であれば統計モデルから除外されるべきものです。しかし通常はデータを取得した段階で上記のような情報は得られないため、モデリングした結果を見て判断することになります。ちょっとやってみましょう。

## 回帰分析の実行(一部省略)
> summary(res.lm <- lm(y ~ ., data = dat))

Call:
lm(formula = y ~ ., data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1262 -0.5336  0.0269  0.6975  1.8557 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.10912    0.10485  -1.041    0.301    
x_real       1.27236    0.14407   8.832 4.45e-14 ***
x_spur      -0.19652    0.09488  -2.071    0.041 *

このデータではx1、x2ともにyに対して有意に(P < 0.05)影響していることとなり*1、本来であれば不要であるx2をモデルから落とすことができません。当然かもしれませんが、両者ともに有意であるためstep関数などを用いた変数選択も意味がありません。

## stepによる変数選択の実行
> step(res.lm)

Start:  AIC=10.25
y ~ x_real + x_spur

         Df Sum of Sq    RSS    AIC
<none>                104.34 10.248
- x_spur  1     4.615 108.95 12.576
- x_real  1    83.901 188.24 67.255

Call:
lm(formula = y ~ x_real + x_spur, data = dat)

Coefficients:
(Intercept)       x_real       x_spur  
    -0.1091       1.2724      -0.1965 

Stanによる推定

それでは上記のデータセットを用いてstanで推定を実行した場合はどのようになるか試してみましょう。まずはstanに渡すデータを定義します。

## stanのdataブロックに渡すデータの定義
dat_stan <- list(
   N = nrow(dat),
   D = 2,
   X = dat[, 2:3],
   Y = dat[, 1]
)

stanのスクリプトは、以下のような内容になります:

data {
   int<lower=0> N;
   int<lower=0> D;
   matrix[N, D] X;
   vector[N] Y;
}

parameters {
   real a;
   vector[D] beta;
   real<lower=0> sig;
}

model {
   for (d in 1:D) {
      beta[d] ~ normal(0, 1);
   }
   Y ~ normal(a + X * beta, sig);
}

特に工夫した点などもなく、単に重回帰を推定するためのスクリプトとなっていますが、このスクリプトを用いてstan()を実行してみましょう*2

## フィッティング
fit_01 <- stan(file = './StanModel/Detect_Confound_Factor.stan', 
               data = dat_stan, 
               iter = 3000,
               chains = 4,
               seed = 1234)

結果はこちらとなります*3

> summary(fit_01)$summary[, c("mean", "2.5%", "50%", "97.5%", "Rhat")]

               mean        2.5%         50%         97.5%      Rhat
a        -0.1055297  -0.3129441  -0.1052853   0.102064821 0.9998115
beta[1]   1.2470319   0.9587653   1.2454155   1.530978156 1.0001602
beta[2]  -0.1827657  -0.3682030  -0.1828429   0.003176939 1.0005802
sig       1.0497539   0.9158756   1.0444427   1.210491329 0.9995408
lp__    -54.9435610 -58.5631720 -54.6322411 -53.160356164 1.0010219

ここでbeta[2]に注目すると95%ベイズ信用区間に0が含まれており、モデルから変数を除外することを検討できるようになりました。

おわりに

分析内容は以上となります。冒頭で紹介したStatistical Rethinkingには他にも多重共線性が生じている時の挙動など、関心を引く話題が多く散りばめられており、非常に読み応えのある本なので、しばらくこれで勉強しようと思います。

*1:というかそうなるようにseedを色々と試しました

*2:ファイル名やフォルダは適当です

*3:ちなみに、最近私のMacではXcodeが更新された際に使用許諾を同意していなかったためコンパイルが通らず、しばらく悩みました