RとStanではじめる ベイズ統計モデリングによるデータ分析入門：サポートページ
『RとStanではじめる ベイズ統計モデリングによるデータ分析入門』のサポートページです。
この記事では、書籍の特徴などの紹介をしています。
本書に使用したサンプルデータとR,Stanのコードは、すべてGitHubから参照できます。
ソフトウェアのインストール方法や実行方法の補足事項なども記しています。
パッケージのバージョンが上がったことによる変更点などは『発行後の補足情報』を参照してください。
RとStanではじめる ベイズ統計モデリングによるデータ分析入門 2019年07月：初版第1刷発行 |
  |
出版社の書籍紹介ページはこちらです。
(2020年7月14日追記)正誤表は出版社の書籍紹介ページにあります。
丸善さん/ジュンク堂書店さんの在庫はこちらから確認できます
。
紀伊國屋書店さんの在庫はこちらから確認できます。
楽天さんの在庫はこちらから確認できます
いただいた書評
『六本木で働くデータサイエンティストのブログ』様
『RとStanではじめるベイズ統計モデリングによるデータ分析入門』は「みどりぼん」に取って替わる次世代の統計モデリング＋ベイジアン入門書
『主にフライフィッシングのブログ』様
ベイズ統計学とRとStanの新しいスタンダード本
『nora_goes_far』様
「RとStanではじめる ベイズ統計モデリングによるデータ分析入門」書評
目次
1．基本情報
出版社 ： 講談社
著者 ： 馬場真哉（このサイト、Logics of Blueの管理人です）
タイトル ： RとStanではじめる ベイズ統計モデリングによるデータ分析入門
発売日 ： 2019年7月10日
簡易目次 ：
第1部 【理論編】ベイズ統計モデリングの基本
第2部 【基礎編】RとStanによるデータ分析
第3部 【実践編】一般化線形モデル
第4部 【応用編】一般化線形混合モデル
第5部 【応用編】状態空間モデル
『実践Data Scienceシリーズ』というシリーズの1冊目となります。
定価は3000円で、消費税10%とすると税込み3300円となります。
2．書籍の特徴
前書きの内容を参照しつつ、書籍の特徴を記します。
以下、前書きより引用
本書は，ベイズ統計モデリングをこれから学ぼうとされる方のための入門書です．ベイズ統計モデリングによるデータ分析を体験してもらうことを目的として執筆しました．
本書では，ベイズ推論とMCMC の組合せを活用して，モデルを推定します．R とStan というともに無料で使えるソフトウェアの組合せを活用して，さまざまな分析を，実際に手を動かして実行します．ベイズ統計モデリングの良書はすでにいくつか出版されています．その中で，本書の特徴は以下のようになるでしょう．
- 数学的な議論が少ない，チュートリアル形式の入門書である
- ベイズの定理などの基本事項をしっかり復習できるようになっている
- Stan の基本だけでなく，推定結果の図示など「分析を実行する」技術も解説している
- brms やbayesplot など，優れたパッケージを積極的に活用している
- GLM からGLMM，DLM，そしてDGLM へと順にモデルを発展させていく
テーマ
完全なる主観ですが、ベイズ統計モデリングは、詳しい人とそうでない人の差が、埋まることなく開き続けているように感じています。
詳しい人は、どんどん高度な技術を推し進めていく。
そうでない人は、後でやろうと思いつつ、スタート地点からなかなか先に進めない。
このギャップを埋めたい、と願って執筆したのが『RとStanで始める ベイズ統計モデリングによるデータ分析入門』です。
本書は、基礎から順番にステップアップしていき、応用的な内容に移っていく構成になっています。
いきなり実践・応用を行うのではなく、ウォーミングアップを念入りにするわけです。
本書は5部構成となっています。
第1部の【理論編】でベイズ統計モデリングの基礎理論を学びます。
第2部の【基礎編】RとStanによるプログラミングの基礎を学びます。
第3部の【実践編】で一般化線形モデル(GLM)という実用的な統計モデルが登場します。
第4部と第5部の【応用編】で、一般化線形混合モデル(GLMM)や状態空間モデル(特に動的線形モデルDLMと動的一般化線形モデルDGLM)を学ぶ構成になっています。
階層構造を持つ複雑なモデルを自由に推定するというのが、やはり目標になるわけですが、そこに至るまでの段階を大切にした書籍となっています。
もちろん基礎理論だけでなくGLMやDLMといった個別のモデルの解説にもページ数を割いたので、RとStanを使ったベイズ統計モデリングのだいご味も味わえるはずです。
文章の書き方も工夫しました。すべての章において「テーマ・目的・概要」を予め説明してから本論に移るようにしています。
本論の流れをつかんでから詳細に移ることで、読者の方が「何が書かれているのかわからないまま読み進める」という状況に陥らないことを目指しました。
また、著者である私が「なぜこの章を書いたのか」という目的を記すことで、読者の方が「この章から何を学べばいいのか」を理解できるようになることを狙っています。
対象読者
統計学の基礎やベイズの定理などの基本事項を学んでみたものの、その有効性がピンとこない、という方を対象としています。
例えば、ベイズの定理の導出では終わらない、より実用的なベイズ統計学について学んでみたいと思った方、あるいは統計モデリングに興味が出てきた方は、本書の対象読者であるといえます。理系文系は問わず、データを分析する必要性がある人はもちろん、データサイエンスが専門でないエンジニアの方でも読めるような内容を目指して執筆しました。
逆に、R やStanを使ったベイズ統計モデリングにすでに明るいという方は、読者として想定していません。ベイズ統計学の数理的な側面もある程度省略しました。基本的な事項を、実装を何度も繰り返すことを通じて学ぶ構成になっています。
ベイズ統計モデリングの上級者をさらなる高みへと持っていく本ではないことに、留意してください。
書籍の構成
前述の通り、本書は5部構成となっています。
簡易目次を再掲します。
第1部 【理論編】ベイズ統計モデリングの基本
第2部 【基礎編】RとStanによるデータ分析
第3部 【実践編】一般化線形モデル
第4部 【応用編】一般化線形混合モデル
第5部 【応用編】状態空間モデル
各々の部の紹介を以下で行います。
第1部 【理論編】ベイズ統計モデリングの基本
本書では、確率・統計、そしてベイズ統計モデリングの基礎理論の復習から入ります。「知識の確認」以上の意味を込めて、この復習パートを、力を入れて執筆しました。
ベイズの定理からスタートする構成ではありません。さらに初歩的な、統計学の基本や確率の加法定理・乗法定理の考え方、確率分布の基本の解説からスタートする構成となっています。
そのうえで確率モデル・統計モデルの基本的な考え方を解説し、そのあとでようやくベイズ推論、そしてMCMCが登場します。
第1部ではRやStanのコードは出てきません。基礎理論を学ぶパートとなります。
数式も必要最小限にとどめました。ベイズ推論の数理的な側面やMCMC(特にNUTS)の詳細なアルゴリズムなどはある程度省略しています。
「なぜこの技術が必要になるのか」という必要性を日本語で解説することで、順を追って無理なくステップアップしていける構成を目指しました。
第2部 【基礎編】RとStanによるデータ分析
第2部では「ベイズ統計モデリングのためのプログラミングの技術」を解説しています。
R言語の復習から入り、Stanの基本的な使い方をこちらで解説します。
Stanを使ったベイズ統計モデリングは、初学者の方にとっては「コードを実行する」あるいは「コードを解読する」だけでもそれなりにハードルが高いと感じられることがあるようです。例えばR言語の使い方でさえ不安があるのに、さらにStanの文法まで覚えるのは少し大変だ、と思われる人もいるでしょう。
第1章でRの基本文法をおさえます。しかし、通り一遍の「R言語入門」だけでは終わりません。インデックスや列名・行名を使ったデータの抽出の方法、forループの仕組み、そして乱数の生成と”乱数の種”の使い方の解説などを行っています。これらを理解しているのといないのとでは、RやStanのコードを読んだ時の理解度が大きく変わると考えたからです。
第2章ではRを用いたデータの要約の方法を解説しています。この章では例えばquantile関数を使ったパーセント点(分位点)の計算の仕方や、カーネル密度推定法の直観的な解説などを載せています。そして第3章ではggplot2の使い方の解説をしています。
ベイズ統計モデリングを行うにあたって「縦横無尽に」使われる技術を、あらかじめ学んでから実践に進んでいただく構成となっています。
第4章から6章までがStanを使ったプログラミングを学ぶパートとなります。
Stanを使って「平均値と分散を推定する」という初歩的な問題を第4章で扱います。第5章では「MCMCを実行した後」に行う各種の評価の方法を解説します。それは収束の評価の方法であったり、事後予測チェックであったりします。最終的にはbayesplotパッケージをつかうことになるのですが、あえてggplot2(やggfortify)を使って”手作業で”グラフを描く方法も記しています。「Rの関数でできることが、私のできることのすべて」となってしまうのは、とてももったいないことであるからです。
第6章では、4章と5章を補足する目的で、Stanの様々な文法の解説を載せています。
第2部までを読了することで、RとStanを使った基本的な分析を行う技術が身につくはずです。
第3部 【実践編】一般化線形モデル
第3部からは、実用的な内容に移ります。一般化線形モデル(GLM)の解説をここで行います。
一般化線形モデルは、様々な統計モデルの基礎といえるモデルです。これ自体が実用的であるだけではなく、一般化線形モデルを発展させることで複雑なモデルを作っていくという「統計モデリングのスタート地点」ともいえる存在です。
ベイズ統計モデリングというと「複雑な階層構造を持つモデルでも、自由に推定できる」というのが、これを使う大きなモチベーションになっていると思います。
しかし、最初のうちは、階層構造がほとんどなくて説明変数も少ない、単純な構造を持つモデルを対象とします。初学者の方が基礎から一歩ずつ進めていけるような、ゆっくりとしたペースで進んでいく構成になっています。
第1章で一般化線形モデルのいくつかの構造を例示します。一般化線形モデルの枠組みで扱える構造を網羅することは目指していません。その代わりに、線形回帰モデルやダミー変数の使い方など、基本的な構造を、ページ数を割いて解説しました。
またデザイン行列を用いたモデルの表現方法も解説しました。入門的な書籍では行列表現を省略することが多いと思います。しかし、中・上級者向けの書籍では当たり前のように行列表現が現れます。より高度な書籍へと進むさいの足掛かりになることも狙って、行列表現の練習をするパートを用意しました。こちらでは、行列の基本的な演算の仕方から始めて、一般化線形モデルの効率的な表現の方法を学ぶ構成になっています。
第2章から5章まで、ずっと単回帰モデルを対象とした解説が続きます。
第2章では単回帰モデルをStanで推定する方法を解説し、第3章では予測の方法を、第4章でデザイン行列を使った実装の方法を解説し、第5章でbrmsという便利なパッケージを使った実装の方法を解説します。
単回帰分析をマスターしてから、応用的な構造に進んでほしいので、このような構成になっています。
第6章以降は、ダミー変数の使い方や、ポアソン回帰モデル、ロジスティック回帰モデルという、重要なモデルの構築方法を解説します。
第10章で交互作用について解説しています。交互作用は、実装は簡単なのですが、解釈を誤りやすいことで有名です。そこで「カテゴリ×カテゴリ」「カテゴリ×数量」「数量×数量」の3パターンにおいて、モデルの実装の方法と結果の解釈の方法を記しました。交互作用に関してここまで丁寧に書かれた書籍は少ないのではないかと思います。
Stanを使ったベイズ統計モデリングは、慣れてくれば、ブロックを組み合わせるようにして複雑な構造を実装していくことができます。
しかし、最初から複雑な構造を対象とすると、個別の構造が理解しにくくなるかと思います。
個別の構造を切り出して、その特徴を吟味する、というやり方で解説する方針としました。
第4部 【応用編】一般化線形混合モデル
第4部からは応用的なモデルに移ります。階層構造を持つモデルを扱うことになり、ここからは「本格的なベイズ統計モデリングだ」と思えるかもしれません。
とはいっても、第4部以降においても「個別の構造を切り出して、その特徴を吟味する」というやり方は変わりません。
第1章ではポアソン回帰モデルにおける過分散の問題を解決するための一般化線形混合モデルの構造を説明します。
第2章、第3章では「グループごとに異なったランダム効果が加わる(あるいは同じグループ内では同じランダム効果が加わっている)」モデルを解説します。
第2章ではランダム切片モデルを、第3章ではランダム係数モデルを解説します。ランダム効果の使い方における、最も初歩的でかつ頻繁に用いられる技法を、単純な事例を通して学ぶ構成になっています。
第5部 【応用編】状態空間モデル
第5部では状態空間モデルを用いた時系列分析を解説します。
基本的に線形のモデルしか対象としませんが、行列表現は避けました。個別の成分を切り出して、その成分ごとに章を分けて解説するようにしています。
第1章で状態空間モデルの基本的な用語を整理します。
第2章では、確率的に変動する水準成分を取り上げます。いわゆるローカルレベルモデルと呼ばれるモデルの解説となります。第3章ではローカルレベルモデルを対象として、予測や補間をする方法を解説しています。
第4章では、時変係数のモデルを扱います。これは回帰係数が時点によって変化することを想定したモデルです。状態空間モデルを使いたいと思う大きなモチベーションとなる構造かと思います。一般化線形モデルとの対比をしながら解説しました。
第5章では、トレンドの構造を解説します。確定的トレンド、確率的トレンド、そして2次のトレンド(平滑化トレンド)とローカル線形トレンドを解説します。「平滑化トレンドだけの解説」とか「ローカル線形トレンドだけの解説」で終わるのではなく、様々なトレンドの構造を対比させながら解説しました。
第6章では確定的周期成分、確率的周期成分を解説し、トレンドの構造と組み合わせた基本構造時系列モデルを導入します。ここまでで、線形ガウス状態空間モデル(別名は動的線形モデル：DLM)の解説はひと段落着くことになります。
第7章では補足的に、自己回帰モデルについて解説しました。状態空間モデルにおいて、自己回帰成分を組み込むことはしばしば行われるからです。構造そのものを知っておくことは有益だと思います。
第8章と第9章では、線形非ガウス状態空間モデル(別名は動的一般化線形モデル：DGLM)を扱います。
DGLMは一般化線形モデル(GLM)の拡張ともいえるし、動的線形モデル(DLM)の一般化であるともいえます。両者との対比をしながら、自然な導入ができるように工夫しました。
DGLMの構造を理解して、これを推定できるようになるのが、本書のゴールとなります。
全体を通して
本書では、線形な構造のみ、説明変数が少ないモデルのみを扱っています。その代わりに、GLM～GLMM～DLM～DGLMへと順番にモデルを発展させていって、基本的なクロスセクションデータから階層構造を持つデータそして時系列データまで、幅広いデータに対して適用できる構造を紹介しています。
応用的なモデルであっても、本書で解説した構造が使われることはしばしばあるかと思います。
より応用的な文献へ進む際の足掛かりとして、ベイズ統計モデリングへの第一歩を踏み出すために、本書を使っていただければ幸いです。
本書に載っていないこと
本書では、網羅性を満たすことよりも、個別のテーマを丁寧に記載することを目指して執筆しています。
そのため、本書でベイズ統計モデリングのすべてのテーマを扱えているわけではありません。
本書では以下のテーマを扱っていません。
- 説明変数が多く含まれるモデル
- 非線形な構造を含むモデル
- WAICなどの情報量基準を用いたモデル選択
- MCMC(特にHMC法の実装であるNUTS)の数理的な側面や詳細なアルゴリズム
- 共役な事前分布を用いて解析的に解を得る方法、など
3．本書のサポート情報
こちらでは、書籍中で使われたサンプルデータとコードを配布しています。
また、環境構築の方法や、コードを実行したときにエラーが出た場合の対応方法も紹介します。
サンプルデータとコード
本書に使用したサンプルデータとR,StanのコードはすべてGitHubから参照できます。
緑色の「Clone or download」というボタンをクリックしてから「Download ZIP」をクリックすると、すべてのファイルをZIP形式でダウンロードできます。
書籍のサンプルコードとデータ。
環境構築の手順
Rは下記のリンクからダウンロードしてください。
https://cran.ism.ac.jp
または
https://cran.r-project.org/
RStudioは下記のリンクからダウンロードしてください。無料の『RStudio Desktop Open Source License』をインストールします。
https://www.rstudio.com/products/rstudio/download/
パッケージのインストール方法は、書籍の内容も参照してください。
RStanのインストールは『RStan Getting Started (Japanese)』という外部リンクも参照してください。
以下のコードを実行すれば、本書で使用するすべてのパッケージをインストールできます。
# 第2部第3章までに使われるパッケージ
install.packages("tidyverse")
install.packages("ggfortify")
# rstanのインストール
install.packages("rstan",
repos="https://cloud.r-project.org/",
dependencies=TRUE)
# Rtoolsのインストール
pkgbuild::has_build_tools(debug = TRUE)
# Stanをより使いやすくする
install.packages("bayesplot")
install.packages("brms")
# boatデータ
install.packages("KFAS")
エラーが出た時の対応
参考までに、当方のPC上で再現したエラーと、その対応を以下に記します。
お手持ちの環境によって、対処法などが異なるかもしれません。この点を理解いただいたうえで、以下の内容を参照してください。
なお、以下の対応はすべて、OSがWindows10の64bitバージョンで確認しています。
基本的には事例2のコンパイル設定の変更だけで動くかと思いますが、この設定時にエラーが出るケースがあったので、少し詳細に手順を記しました。
rstanやRのバージョンによって、作業手順が変わる可能性があります。
あまりPCの操作に詳しくない方でも実行でき、追加のソフトウェアのインストールや仮想環境の構築を行わなくても実行できる方法を選んで載せたつもりです。他の対処方法があるかもしれません。
また、繰り返しになりますが、下記手順は、あくまでも当方の実行環境下で確認された方法となります。
PATHの追加など、PCの設定変更に関しては自己責任で実行をお願い致します。
恐れ入りますが、上記の点、ご理解ください。
事例1：パッケージのインストールの際にエラーが出る
RStudioのアイコンを右クリックして「管理者として実行」してください。
RStudioを管理者として実行した状態でパッケージのインストールを行うと、うまくいくことがあります。
事例2：Stanのコンパイルが通らない
『C++14 standard requested but CXX14 is not defined』のようなエラーがでて、コンパイルが通らないことがあります。
この場合は『RStanのインストール時に C++14 standard requested but CXX14 is not defined と怒られた時の対処｜Qiita』という外部リンクを参照して、記載されているコードを実行してください。
外部リンク先のコードを実行すると、Sys.getenv("HOME")で出力されるフォルダ上に、『.R』という名前のフォルダが作られ、そのフォルダの中に『Makevars』という名前の拡張子のないファイルが作成されます。
どんなファイルなのかが気になる場合、あるいはうっかりで2回以上、上記のコードを実行してしまった場合は、コンソール上でSys.getenv("HOME")と実行してください。フォルダ名が出力されます。
『Makevars』という名前のファイルに、以下のように記載されていればOKです。
CXX14FLAGS=-O3 -Wno-unused-variable -Wno-unused-function
CXX14 = g++ -std=c++1y
同じ設定が何度も記されていると問題です。この場合は一度ファイルを削除してから、外部リンク先のコードを再実行してください。
その他、コンパイルの設定は様々変えることができます。詳細は『RStan Getting Started (Japanese)』という外部リンクを参照してください。
事例3：事例2のコードがエラーで動かないとき
全角文字が含まれるフォルダがHOMEフォルダに指定されてしまい、事例2の外部リンクに載っているコードを実行するとエラーが出ることがあるようです。
著者の場合は「ドキュメント」という全角カタカナがフォルダ名に入っており、エラーになっていました。
様々な対策がありますが、HOMEフォルダを変えてしまうのが一つのやり方です。
Cドライブ直下に、例えば『C:\r_home』というフォルダを作ってください。
そして、以下のコードをRStudioで実行します。
Sys.setenv("HOME" = "C:/r_home")
Sys.setenv("R_USER" = "C:/r_home")
上記のコードを実行した後に、事例2の手順に従って『Makevars』ファイルを作成します。
すると、『C:\r_home』というフォルダの中に『.R』というフォルダができ、その中に『Makevars』という名前の拡張子がないファイルが作成されるはずです。
次回以降は、RStudioを起動させるたびに、以下のコードを実行します。
Sys.setenv("HOME" = "C:/r_home")
Sys.setenv("R_USER" = "C:/r_home")
本来は、設定を恒久的に変えてしまうのが良いかもしれません。
ほかの分析に悪影響を及ぼさないようにするため、ここでは、上記のやり方を記しておきます。
事例4：RtoolsのPATHが通っていない
事例2～3の時とはエラーメッセージが変わって『Compilation ERROR, function(s)/method(s) not created! sh: g++: command not found』などというエラーが出ることがあります。
この場合はRtoolsのPATHが通っていないのが理由になります。
PCによって異なる可能性がありますが、『C:\RBuildTools』などにRtoolsがインストールされているはずです。
64bitのOSをお使いの場合は、WindowsのPATHに以下の2つを追加してください。
C:\RBuildTools\3.5\bin
C:\RBuildTools\3.5\mingw_64\bin
WindowsにPATHを追加する方法は、いくつかあるでしょうが、『windows path 追加』などで検索いただければ、手順が載っているかと思います。
参考になれば幸いです。
参考資料
主にWebで閲覧できる、ベイズ統計学やStanプログラミングの参考資料のリンクを以下に記載します。
・Stan-Documentation
拙著「RとStanではじめる ベイズ統計モデリングによるデータ分析入門」で主に使用しているソフトウェアであるStanの公式ドキュメントです。英語です。
・stan-reference(日本語)
Stanのリファレンスです。有志の方々が日本語に翻訳作業を進めてくださっています。
・StatModeling Memorandum
「StanとRでベイズ統計モデリング (Wonderful R)」の著者の方のWebサイトです。
ベイズ統計モデリング全般に関して、大変有益な情報が多く公開されています。
・生態学のデータ解析 – ベイズ統計 & MCMC
「データ解析のための統計モデリング入門: 一般化線形モデル・階層ベイズモデル・MCMC」の著者の方のWebサイトです。
統計モデリングの初歩からベイズ統計モデリングまで、様々な情報が公開されています。
・作って遊ぶ機械学習
「ベイズ推論による機械学習入門」の著者の方のWebサイトです。
ベイズ推論や機械学習法など、様々なトピックが公開されています。
・六本木で働くデータサイエンティストのブログ
Stanのインストール方法から分析事例まで、様々なトピックが公開されています。
・Rで楽しむStan
「Rで楽しむベイズ統計入門」の著者の方のWebサイトです。
現在は工事中と記載がありますが、Stanのサンプルコードなどが公開されています。
・Kosugitti’s BLOG
Stanを使った分析事例など、様々なトピックが公開されています。
なお、上記リンク先では状態空間モデルを用いた欠損値補間のコードが記載されています。
欠損値補間には様々な方法がありまして、以下の記事も参考になります。
・{rstan} 欠損値あり時系列データのモデル推定
この辺りはStatModeling Memorandumの以下の記事も参考になります。
・RStanで『予測にいかす統計モデリングの基本』の売上データの分析をする
・不等間隔の状態空間モデル
なお、拙著「RとStanではじめる ベイズ統計モデリングによるデータ分析入門」では、可読性を鑑みて松浦(2016)「StanとRでベイズ統計モデリング (Wonderful R)」記載の方法を採用しています。しかし、欠損値補間の方法は1つには限らず、様々あることは覚えておくとよいかと思います。欠損値補間の使用頻度は高いと私は感じます。
・Sunny side up!
Stanを使った分析事例など、様々なトピックが公開されています。
・nora_goes_far
brmsの使用例など、様々なトピックが公開されています。
・Stan Advent Calendar 2018
Stanを用いた分析事例など、様々なトピックのリンク集となっています。
発行後の補足情報
2023年2月4日追記
※当記事にコメントをいただきました。情報提供、ありがとうございます。
Rを4.2系にあげると、RStanが実行しにくくなるようです。以下の記事(外部サイト)を参考にしてください。
・Rを4.2系にバージョンアップしたらRstanの導入でつまずいた話｜豆蔵デベロッパーサイト
・Stan & R 4.2 on Windows｜The Stan Blog
2019年9月7日追記
brmsのバージョンが2.10.0に上がったことに伴う変更点です。
ロジスティック回帰の回帰曲線を図示するp225のコードは、以下のように修正してください。
eff <- marginal_effects(glm_binom_brms,
effects = "nutrition:solar",
conditions = data.frame(size = 10))
plot(eff, points = TRUE)
2020年7月14日追記
あたらしいバージョンのbrmsを使う場合は、『marginal_effects』ではなく『conditional_effect』関数を使ってください。
例えばロジスティック回帰の回帰曲線を図示するp225のコードは、以下のようになります(brms 2.13.0で動作を確認しました)。
eff <- conditional_effects(glm_binom_brms,
effects = "nutrition:solar",
conditions = data.frame(size = 10))
plot(eff, points = TRUE)
※ 本書と同じコードで同じ結果を再現したい場合は、本書に記載の通りのバージョンに合わせてください。
2020年7月14日追記
R4.0.0以上のバージョンのRを使うとき、Stanとの連携時にエラーが出ることがあります。
管理人が実行した手順を紹介します。なお、OSはWindows10、Rは version 4.0.2で確認しました。
なお、本書と同じコードで同じ結果を再現したいならば、本書記載の通りの古いバージョンで実行するようお願いいたします。
バージョンアップは慎重に行ってください。
下記の作業は、すでにRstanを実行できている人が、バージョンアップすることを想定しています。
環境によっては、この手順の通りに実行してもうまくいかないことがあります。参考程度にご利用ください。
古いソフトウェアをアンインストールする措置があります。もとに戻すのが大変なので、自己責任で実行してください。
参考にした記事は以下の通りです。
https://discourse.mc-stan.org/t/r-4-0-rstan-and-you/15271
https://discourse.mc-stan.org/t/dealing-with-r-4-0/14586
1．古いRtoolsをアンインストールする
Rtools 3.5などがあるならば、全てアンインストールします。
Rtools 3.5のPATHも全て削除します。
※必要かどうか不明ですが、管理人は古いバージョンのR本体も全てアンインストールしました。
2．以下のURLに従ってRtools 4.0をインストールする
https://cran.r-project.org/bin/windows/Rtools/
RStudioを最新バージョンにするように指示があるのでそのようにします。
Rtools 4.0をインストールしたうえで、PATHも追加しておきます。
3．この状態でrstanパッケージをインストールして、コードが動くか確認します。
3番の対応でエラーが出た場合は、さらに以下の作業を行います。
4．古いMakevarsファイルをすべて削除、あるいは名称変更をする(Documentフォルダもすべて)
なお、管理人は、Documentフォルダに転がっていた古いMakevarsファイルのせいで、何度もコンパイル時にエラーが出ました。
それなりに手間がかかります。Rのバージョンアップは慎重に行ってください。
参考になれば幸いです。
RとStanではじめる ベイズ統計モデリングによるデータ分析入門 |
|
更新履歴
2019年06月10日：新規作成
2019年06月17日：GitHubへのリンクを追加。データ配布に関する文言を一部修正。
インストール手順や、コンパイルエラーへの対応を追記。
2019年07月01日：参考資料へのリンクを追加
2019年07月09日：出版社と紀伊国屋書店さんへのリンクを追加
2019年09月07日：brmsがバージョン2.10.0に上がったことに伴う変更点を追記
2019年10月06日：重版した旨と、いただいた書評のリンクを追記
2019年10月19日：いただいた書評のリンクを追記
2020年03月17日：重版した旨を追記(重版自体は2月)
2020年04月27日：重版した旨を追記・消費税を10%に変更
2020年07月14日：brmsとRのバージョンアップに合わせて補足事項を追記。重版情報・正誤表の情報も追記。
2023年02月04日：R4.2についての注意事項を追記
2023年07月29日：10刷の重版日付を追記
“RとStanではじめる ベイズ統計モデリングによるデータ分析入門：サポートページ” に対して71件のコメントがあります。
コメントを残す
このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください。
時系列分析の本では大変お世話になりました。
301p以降の詳しい説明を求めておりましたのでとても楽しみです！
Arai様
コメントありがとうございます。
管理人の馬場です。
興味を持っていただき、ありがとうございます。
どちらかといえば、基礎的な内容が主になりますが、こちらの新刊もお役に立てれば幸いです。
https://github.com/logics-of-blue/book-r-stan-bayesian-model-intro
よりダウンロードした圧縮ファイルを解凍すると、なぜか漢字が入っているファイル名の漢字部分がすべて文字化けしている^_^特に支障はないと思いますが。
なお、これまでに下記2冊も購入してコードを確認しながら勉強しました。いずれも大変分かりやすくて大いに役に立ちました。
Pythonで学ぶあたらしい統計学の教科書
時系列分析と状態空間モデルの基礎 RとStanで学ぶ理論と実装
申し訳ございません。
上記ファイル名にある漢字が文字化けの現象について、解凍をwindows10のexplorerを使えば、文字化け現象は起こらないことが分かりました。お騒がせしました。LhaPlusという解凍ツールでは文字化けしてしまった。おそらくファイル名の文字コードがLhaPlusのデフォルト設定と合わなかったと思われます。
李様
コメントありがとうございます。
管理人の馬場です。
内容、把握しました。
問題なかったようで良かったです。
購入して少し読みました。
行列表現をさけて初学者でも理解しやすいようにしている点が良かったです。
質問があります
ベイズ統計モデリングではpythonではなく、Rやstanが使われる例を多く見るのですが理由は何なのでしょうか。
pythonでのpyroやEdwardなど確率プログラミング言語と比べてどのような点がすぐれているのでしょうか。
gg様
コメントありがとうございます。
管理人の馬場です。
ソフトウェアの選択に関しては、用途や目的によるので一概には答えにくいです。
以下では、当方の私見を簡単に述べます。
RとStanの組み合わせが唯一の選択肢というわけではないでしょう。
pyroやPyMCそしてTensorFlow Probabilityは十分な選択肢になると思います。
RとStanの組み合わせの良いところは様々ありますが、拙著で強調しているのは有用なライブラリが多くあることです。
例えばbayesplotやbrmsというライブラリを併用することで、短くて可読性の高い分析コードとなります。
StanではHMC法の1実装であるNUTSが用いられており、収束の速さなどの速度面でも、ほかの選択肢に引けを取らない十分な性能があります。
研究者・実務者を含めて多くのユーザーがおり、情報が得られやすいのも利点の一つだと思います。
導入のハードルの低さ、性能の高さ、豊富な実績。これらを鑑みて、当方としては、ベイズ統計モデリングをこれから始めるというユーザーには、RとStanの組み合わせをお勧めしています。
P331の補足：状態空間モデルと自己回帰モデルの部分で
「自己回帰モデルは観測誤差がまったくないことを想定」との記述があるのですが、観測誤差ではなく、過程誤差ではないでしょうか？
自己回帰モデルではy_t~Normal(β_0 + B_1 * y_t-1, σ^2)のようにσ^2を観測誤差としており、状態モデルを仮定していないため過程誤差がないと思うのですが。
gg様
コメントありがとうございます。
管理人の馬場です。
誤差の解釈の問題ですね。
いろいろの解釈ができるかと思いますが、
仮に、定数項がない自己回帰モデルを想定してみます。
そして回帰係数β_1の値が「１」だったとします。
すると、これはランダムウォーク系列となります。
ローカルレベルモデルは「ランダムウォーク系列に観測誤差が加わる」という形で定式化されます。
そのため、定数項がなく、回帰係数β_1の値が「１」である自己回帰モデルは、
観測誤差のないローカルレベルモデルであるとみなすことが可能です。
このことから類推されるように
自己回帰モデルの誤差項を「過程誤差」とみなすことは、自然かと思います。
ちなみに、本書の内容からはやや外れますが、
いわゆる「和の定理」から、AR(p)にホワイトノイズが加わった系列はARMA(p,p)となります。
(ハーベイ1985『時系列モデル入門』より)
このため、自己回帰移動平均モデル(ARMA)は、観測誤差のある自己回帰モデルとしても解釈できます。
この点からみても自己回帰モデルを「観測誤差がないモデル」とみなすのが自然です。
自己回帰項は、状態方程式に組み込まれることも多いです。
この場合もやはり誤差項を「過程誤差」とみなした方が、
モデルの解釈がしやすいです。
p336のStanコードのmodelブロックにて、
// 弱情報事前分布
s_w ~ student_t(3, 0, 10);
……
というt分布に関するコードがあります。
私にとっては難しい為(汗^_^)、もう少し詳しくご教授を頂ければ有難く存じます。
(1) student_t関数の引数はそれぞれ自由度が3、平均値が0、分散(標準偏差)が10だと認識しております。
これらの引数の適用に当たって、何を根拠或いは基準にしていただいたでしょうか？
例えば、このケースでは、自由度が3としておりますが、これは、「boat」のデータのヒストグラムを事前分析して、
尖度が高くて、裾野が少々広い分布だったことが分かった為に決めて頂いたでしょうか？
また、分散(標準偏差)は10となっていますが、決め手は？
(2) Stanに関する処理プロセスについて、全くの初心者なので、この
// 弱情報事前分布
s_w ~ student_t(3, 0, 10);
コードは、ベイズの事後確率更新において、初期の1回きりでしか呼び出されないでしょうか？それとも事後確率更新の度に、
呼び出されるでしょうか？
(3) 次の質問の内容は少し本書と逸脱しますが、p336のStanコードの観測値過程は偶々モデルの為にベルヌーイ分布になっていますが、
例えばローカルレベルモデルにおいて、一様事前分布ではなく、p336の例のように、
data {
int T;
vector[T] Y;
}
parameters {
vector[T] mu;
real s_mu;
real s_Y;
}
model {
// 弱情報事前分布
s_mu ~ student_t(4, 0, 0.2);
mu[2:T] ~ normal(mu[1:(T-1)], s_mu);
Y ~ normal(mu, s_Y);
}
という弱事前分布student_t、正規分布の過程誤差、正規分布の観測誤差のような組み合わせでも、状態空間モデルの設計或いは実装上、
ローカルレベルモデルという概念からすれば、問題なるのでしょうか？或いはローカルレベルモデルにこだわらず、
ただ単に上記stanコードのように実現したい場合、理論上或いは設計上問題なるのでしょうか？
質問が多くてお手数をかけまして申し訳ございません。
李 士芳様
コメントありがとうございます。
管理人の馬場です。
返事が遅れてすみませんでした。
(1)
事前分布に関して、明確に××の分布を使うべき、とは決めにくいです。
「決め手」といわれると難しいですが、
収束を良くしつつ、それでいて制約が強く過ぎることの無いように、ということで選びました。
(2)
イメージとしては事前分布は「事前に設定しておくもの」というところでしょうか。
「途中で呼び出される」というものではないです。
事前分布を予め設定しておき、データに基づき事前分布を事後分布へと更新します。
(3)
正規分布を想定したローカルレベルモデルなどでも、
弱情報事前分布を使うことは当然あり得ます。
収束が悪かった場合などには、弱情報事前分布を設定することを検討します。
本筋とは関係ないですが、李様のサンプルコードでは、
標準偏差が負の値を取る可能性があるため、
以下のような指定があった方が良いかと思います。
real<lower=0> s_mu;
real<lower=0> s_Y;
馬場様へ
ご回答を頂きまして、ありがとうございました。
パラメータの設計について、結局は試行錯誤により、結果から良いものを選び出しますね。
(3)についてのアドバイスを有り難うございました。
新刊が出ていたのですね！ 仕事でRとstanを使っていますが、自分の知識や理解に自信がないため、入門から改めて勉強したいと思います。こちらのサイトは、カテゴリーが分かりやすく分けられていて、記述が丁寧かつ分かりやすいので、サポートページを開いていただいていることは、とても安心です。早速注文します！
TK様
コメントありがとうございます。
管理人の馬場です。
いつもありがとうございます。
拙著がお役に立てれば大変幸いです。
章毎にテーマ、目的、概要がまとめられていて、非常に学習しやすいです。
P225で回帰曲線の図示をしようとしたところ、試行回数が10で固定できていなかった様子で、
“Setting the number of trials to 1 by default if not specified otherwise.”、
というメッセージとともに下の方に偏った回帰曲線が図示されてしまいます。
試行回数を指定する方法をご教示ねがえないでしょうか。
宜しくお願い致します。
初心者様
コメントありがとうございます。
管理人の馬場です。
brms のバージョンが2.10.0に上がったことが原因である、仕様変更のようです。
marginal_effects関数において『conditions = data.frame(size = 10)』と指定すると正しく動きます。
ご連絡いただき、誠にありがとうございます。
馬場様
読了しました。今まで読んだベイズ本、統計モデリング本の中で一番分かりやすく、自分のレベルに合っていると思いました。これから座右の書として傍らに置いて、仕事で使っていきたいと思います。素晴らしいご著書を送り出していただいた感謝を込めて、僭越ながらブログに感想を書きました(https://ameblo.jp/waltham0012/entry-12527723276.html)。
本で気になった点が2つあります。
p.69 真ん中
カーネルの比Kernel(-1.27)/Kernel(-0.94)はまたしても小さな値に・・・
は
カーネルの比Kernel(-1.78)/Kernel(-0.94)はまたしても小さな値に・・・
でなくてよいでしょうか。
p.165 式(3.35)
c_ijとΣの間にイコール記号はなくてよいでしょうか。
私の理解不足でしたら申し訳ございません。
TK様
コメントありがとうございます。
管理人の馬場です。
書評を書いていただき、誠にありがとうございます。
書籍を詳細に読んでいただいたようで、著者としてとてもうれしいです。
誤植のご指摘ありがとうございます。
内容確認の上、出版社と共有いたします。
まずは記載に誤りがあったことをお詫びいたします。
また、ご指摘いただけたこと、感謝いたします。
馬場様
度々すみません。
誤植の箇所を一つ、お伝えするのを忘れていました。
p.276 i.i.dとはのところ
independend ⇒ independent
最後がdではなくt
既にご承知でしたら申し訳ございません。
TK様
コメントありがとうございます。
管理人の馬場です。
ご連絡ありがとうございます。
こちらは電子版と第2刷では修正されているかと思います。
誤植がありまして、誠に申し訳ございません。
ご指摘いただけたこと、感謝いたします。
目下、勉強させていただいております。
色分けもされて見やすく、勉強がはかどります。
グラフの描き方なので、こちらでお尋ねするべきことではないかもしれませんが、ご存じでしたらご教示いただければと思います。
２１７ページ（図３・８・３）で、weatherは晴れと曇りの2種類ですが、どちらか一方のみの状況を確認したい場合は何か方法はございますでしょうか。
凡例が３種類、４種類とある場合には個別に確認したくなる気がしました。
できれば、marginal_effects関数とplot関数を併用するという本書の方法を踏襲できればありがたいです。
ご教示いただければ幸いです。
KAZU様
コメントありがとうございます。
管理人の馬場です。
返信が遅れてしまい、失礼いたしました。
個別にグラフを描く場合は、以下のようにconditionsを設定するのが簡単です。
以下のコードで、天気がcloudyのグラフが描けます。
data.frame(weather = c(“sunny”))に変更したら、晴れの日のグラフが描けます。
conditions <- data.frame(weather = c("cloudy"))
conditions
set.seed(1)
eff_pre_condition <- marginal_effects(
glm_pois_brms,
method = "predict",
effects = "temperature",
conditions = conditions,
probs = c(0.005, 0.995), points = TRUE
)
plot(eff_pre_condition, points = TRUE)
ご丁寧にありがとうございます。
２１６ページ下段のプログラムに追加したところ、1つのみ描けました。
お世話になります。
以下の件で質問させてください。
２２３ページ上部のmodel構造の指定部分（応答変数部分）で
germination | trials(size)
となっています。
２２２ページに「size個のうちgermination個だけが発芽したことを示すため」と書いてあるのですが、もう少し詳しく教えていただけないでしょうか。
たとえば、
・「 | 」 という記号は条件付きを表しているでしょうか。そうだとしますと、なぜgerminationだけではモデルとして不適切になるのでしょうか。（図３・９・１と３・９・２の縦軸は違うのですが、プロットされている点は同じように見えます。）
・trialsというのはRに用意されている関数でしょうか。そうだとしますと、どういう働きの関数でしょうか。（?trialsとしてみたのですが、brmsやrstanパッケージ内のヘルプになかったようです。）
お手数をおかけしますが、よろしくお願いいたします。
KAZU様
コメントありがとうございます。
管理人の馬場です。
返信が遅れてしまい、失礼いたしました。
> なぜgerminationだけではモデルとして不適切になるのでしょうか
これは、試行回数が指定されていないことが問題です。
二項分布において、成功回数と試行回数の両方を明示するための、特殊なformula構文が必要で、
それが「germination | trials(size)」です。
「germination | trials(size)」という書き方は、brmsのルールのようなものですね。
縦棒(|)には条件付という意味はなく、あくまでも「成功回数｜試行回数」であることを示すためのルールになります。
「trials」に関しても、brmsのformulaにおいて、試行回数という特別な意味を持つことを表記するルールとなります。
馬場様
ベイズ統計モデリングに関する大変分かりやすい書籍を出版くださり、有り難うございます。
記載された全てのモデルを自身で再現することで、楽しみながら理解することが出来ました。
P264のランダム係数モデルを実行する際、以下のエラーメッセージが出てしまいました。
Compiling the C++ model
recompiling to avoid crashing R session
Start sampling
Warning message:
Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess
chainsを12、adapt_deltaを0.99にして再度試したものの、同じメッセージが表示されました。
rstanのバージョンは2.19.2です。
大変恐縮ですが、改善策をご教授いただければ幸いです。
何卒よろしくお願いします。
飯田真也様
コメントありがとうございます。
管理人の馬場です。
brmsのバージョンが上がったことに伴い、
推定結果が少し変わってしまったようです。
以下の設定にすればワーニングは出ないはずです。
今回の事例にかかわらず、ワーニングが出る場合は、
adapt_deltaだけではなく、iterやwarmupなどの変更も検討すると、
うまくいくことがあります。
# ランダム係数モデル
glmm_pois_brms_keisu <- brm(
formula = fish_num ~ temperature + (temperature||human),
family = poisson(),
data = fish_num_climate_4,
seed = 1,
iter = 7000,
warmup = 6000,
control = list(adapt_delta = 0.98, max_treedepth = 15)
)
お忙しい名k、立て続けにご回答いただきまして、ありがとうございます。
こちらもクリアになりました。
ベイズ統計を使えるようになりたいと独学中の者です。「RとStanではじめる ベイズ統計モデリングによるデータ分析入門」が私には一番わかり易かったのでこの本を読み進めています。最近、チラ見している数学者のツイッターの「ベイズ統計とベイズの定理は違う話です。～」
https://twitter.com/genkuroki/status/1164570497167716354
の記載内容との本のp53～は、別物なのか同じことなのか、力不足で判断ができません。
お教え頂けるとモヤモヤが解消し助かります。よろしくお願いいたします。
うどん様
コメントありがとうございます。
管理人の馬場です。
返信が遅れて失礼いたしました。
『ベイズ統計とベイズの定理は違う話です。』における
『違う話』の意味が当方にはわかりかねるので何とも言えません。
ベイズの定理そのものは、参考ツイートにある通り
『ベイズの定理は条件付き確率の定義の自明な書き直し』です。
この意見には同意できます。
逆に言えば、当然成り立つことなので、
ベイズの定理が「誤っている」ことはないはずです。
ベイズの定理をベイズ統計で「一切扱わない」こともありません。
むしろ条件付確率を扱う計算でベイズの定理が登場するのは自然でしょう。
「ベイズ統計学において、ベイズの定理を使ってはならない」ということは無いと思います。
ベイズの定理をベイズ統計で「一切扱わない」ことがないのにかかわらず
『違う話』であると主張する意図が当方にはわかりかねますが、
拙著の計算式が間違っている
(例：ベイズの定理の式が誤りである、
ベイズ推論の際にベイズの定理を用いているのは誤りである)
ということは無いのではないかという認識です。
先のツイートは
「ベイズの定理を使っているから、それがベイズ統計学なのだ」という主張は良くない、
という意図なのではないかと想像します。
これは確かにその通りで、ベイズの定理は条件付確率を扱う際にしばしば登場しますので、
ベイズの定理が現れるたびに「これはベイズ統計学なのだ」と主張するのは無理があると思います。
当方の勝手な想像で恐縮ですが、
仮にこのような意図であれば、当方も同意できますし、
拙著との矛盾もないのではないかと思います。
(拙著でも、第1部において、「計算の過程でしばしば登場するもの」としてベイズの定理を導入していますので)
馬場さま
丁寧で納得のいくご回答を頂きありがとうございました。
ご著書を頼りに、独学を進めて参ります。
ありがとうございました。m(__)m
はじめまして．
馬場さんのご出身の大学で，鳥の研究をしている大学院生です．
いつも馬場さんの本を拝読しながら，さまざまな統計を勉強しております．
現在，「RとStanではじめるベイズ統計モデリングによるデータ分析入門」を読みながら，データの解析をしております．
私は鳥につけたGPSロガーデータを扱っております（GPSポイント（サンプルサイズ）は1000個ほど）．特に，GPSに位置（緯度，経度）とともに記録される高度について，解析をしており，鳥が20m以上の高さで飛ぶのはどのような場所か，を明らかにしたいと思っています．どうもデータをプロットすると，島から近い場所のみで高くなる傾向がみられています．この関係をGAMのような非線形回帰でモデリングしたいと考えております．ただ，GPSロガーの高度のデータは，非常に誤差が大きく，実際の高度との誤差（平均±標準偏差）は0±30mほどあり，できればこの誤差についてランダム効果として変数に組み込めないかと考えております．
モデルの式としては，高度 ~ s(繁殖地からの距離）＋ランダム効果（GPSの誤差（正規分布で平均0，標準偏差30））とし，事前にランダム効果の誤差分布を指定できればうまく表現できるののではと思い，brmsを使ってできるかと思って調べてみたのですが，brmsではこのようにランダム効果の分布の指定はできないようでした．そのため，実際にstanコードを書いて．．．と思ったのですが，非線形回帰の場合のコードの書き方も含めよくわからず，四苦八苦しております．
このモデリングが正しいかも不安ではあるのですが，もしこのような場合のコードの書き方等ご存じでしたらお力を貸していただけませんでしょうか．
お忙しいところ大変恐縮ですが，どうぞよろしくお願いいたします．
JO様
コメントありがとうございます。
管理人の馬場です。
返信が遅れて失礼いたしました。
こちらは「ランダム効果の大きさに対して、事前の知識があり、それを生かしたい」ということでしょうか。
それならば、例えばランダム効果の大きさに対して、弱情報事前分布を指定すると良いかもしれません。
例えば『
normal(30, 5)』のような弱情報事前分布を指定するなどのやり方です。(標準偏差の大きさは適当です。感度分析をしても良いかもしれません)
これならbrmsで指定できます。
参考になれば幸いです。研究頑張ってくださいね。
実際に動かして、グラフを表示するまでの流れが一貫しており、非常に分かりやすいです
状態空間モデルを図示する関数を読むこむ際に
> # 状態空間モデルの図示をする関数の読み込み
> source(“plotSSM.R”, encoding=”utf-8″)
Error in source(“plotSSM.R”, encoding = “utf-8”) :
plotSSM.R:3:35: unexpected INCOMPLETE_STRING
2: state_name, graph_title, y_label,
3: date_labels = “%Y?
^
In addition: Warning message:
In readLines(file, warn = FALSE) :
invalid input found on input connection ‘plotSSM.R’
とエラーが表示されます。何が不味いでしょうか？
Ub様
コメントありがとうございます。
管理人の馬場です。
返信が遅れて申し訳ありません。
当方の環境では正常に動くようでした。
エラーメッセージを見るに、ファイルを開いたとき、
何か文字化けをしていることは無いでしょうか。
『plotSSM.R』ファイルの文字コードが、
正しくUTF-8になっているかどうか、確認していただければ幸いです。
ゴールデンウイークを使って、本書を読ませていただきました。以前、ベイズ統計やSTANに興味を持って勉強を始めましたが途中で躓いていました。 今回、本書とサポートページのコード実行で、自分でも何とか使えそうな気がしてきました。 本書は、初学者への説明が丁寧でスッキリしており、関連性の説明が分かりやすくとても助かりました。
ありがとうございます。
今回、コード実行に関して、少しトラブルがあったので報告しておきます。
当初、環境構築時にRの最新版（R4.0.0）をインストールし、これをベースに環境構築しました。その結果、まずは、スタンコード 2-4-1が動作しませんでした。
スタンコード 2-4-2は動作したので、この内容に書き換えて進めていましたが、
“3-5-brmsの使い方.R”の library(brms) でエラーになり、色々試しましたが解決できませんでした。(使用した brms のバージョンは、ver2.12.0 です。)
バージョン不整合の問題と考え、以下から、R3.5.3をインストールし環境を再構築したところ、全て動作するようになりました。
https://cran.r-project.org/bin/windows/base/old/3.5.3/
brmsのバージョンが進めば大丈夫かもしれませんが、R4.0.0で、brms(ver.2.12.0)は使用できないのかもしれません。
同じようにR4.0.0をインストールして、躓く方もいらっしゃるかと思い報告させていただきました。
SS様
コメントありがとうございます。
また、拙著をお読みいただき、ありがとうございます。
Rは、大きなアップデートがあった後の動作がしばしば不安定になりますね。
情報を共有いただき、ありがとうございます。
馬場 様
この度本を購入させていただき、ベイズ解析の勉強に大変役に立ちました。
変な質問になるかもしれませんが、実際論文を書く場合では、
頻度主義とベイズ主義の解析手法を１つの論文の中で混ざって使うのがあり得るでしょうか。例えば調査場所による光環境や地形環境の違いにはTukey HSDの検定方法を使って、光環境が植物の成長等に与えいる影響を見たい時に、また階層ベイズで分析してみるとか、こういうことは論文の中で一般的なやり方になるでしょうか？
あるいは、頻度主義かベイズ主義か、論文の中でどっちかを一貫する必要があるでしょうか？
どうぞよろしくお願い致します。
ゴ様
コメントありがとうございます。
管理人の馬場です。
拙著をお読みいただきありがとうございました。
「あり得る」かどうかでいえば、もちろんあり得ます。
ベイズモデリングとTukey HSDの検定を両方使うことが、
禁止されているわけではありませんから。
論文として公開されるかどうかも、レフェリー次第かと思います。
とはいえこのやり方が「好ましいか」や
「このようにすべきか」といわれると、
少し違うかなと思います。
ベイズモデリングをすると、自分の欲しい対象の分布を、
比較的簡単に得ることができます。
『調査場所による光環境や地形環境の違い』も
ベイズモデリングで評価することができるはずです。
Tukey HSDの検定を別途使うべき強い理由は無いように思います。
本の手順に沿って進めていたのですがstanのコードに入ってからエラーがでてプログラムがいっさい動かなくなりました。
一旦すべてアンインストールし、Rのバージョンを本のものと揃えたところ動くようになったのですが、R最新版ではstanを使うことはできないのでしょうか？
OSはwindowsです。
axio様
コメントありがとうございます。
管理人の馬場です。
拙著をお読みいただきありがとうございます。
新しいバージョンのRで動かす場合のメモをこちらに記載しました。参考になれば幸いです。
ただし、手間がかかる上に元に戻すのが大変なので、バージョンアップは慎重に行って下さい。
https://logics-of-blue.com/r-stan-bayesian-model-intro-book-support/#3-5
brmsがインストールできません。Rのバージョンは3.5.3です
axio様
コメントありがとうございます。
管理人の馬場です。
brmsの要件はR (≥ 3.5.0)なので、
R3.5.3でもbrmsはインストールできるはずです。
参考：https://cran.r-project.org/web/packages/brms/index.html
念のため当方のPCでも試してみましたが、インストールできました。
何か環境固有の原因があるかもしれません。
RStudioを「管理者として実行」するとうまくいくこともあるようです。
馬場 様
ベイズ統計モデリングの勉強のため本書を購入しました。
前提知識がない状態からでも理解しやすく、楽しみながら学んでおります。
ありがとうございます。
p159の(3.11式)と、p213の(3.50式)の部分で質問があります。
P158(3.10式)のリンク関数の対数関数(log())を、逆関数の指数関数(exp)で変換していますが、(3.11式)は以下の式が正しいのではないのでしょうか？
λi = exp(β0 + β1xi1 + β2xi2 + β3xi3)
yi ~ Poiss(λi)
現行の記載では、
(3.10式) log(λi)＝ λi (3.11式)、(3.10式) λi＝ exp(λi) (3.11式)
の関係性になるため、(3.10式)λi ≠ λi(3.11式) となる気がします。
こちらの理解不足でしたら申し訳ございません。
HZ様
コメントありがとうございます。
管理人の馬場です。
拙著をお読みいただきありがとうございます。
下記の理由で、おそらく問題ないかと思います。
HZ様の提示された式
①λi = exp(β0 + β1xi1 + β2xi2 + β3xi3)
②yi ~ Poiss(λi)
①を②に代入
③yi ~ Poiss(exp(β0 + β1xi1 + β2xi2 + β3xi3))
拙著の式
①λi = β0 + β1xi1 + β2xi2 + β3xi3
②yi ~ Poiss(exp(λi))
①を②に代入
③yi ~ Poiss(exp(β0 + β1xi1 + β2xi2 + β3xi3))
③の式では、同一になります。
馬場 様
丁寧な解説ありがとうございます。
応答変数と線形予測子をつなぐ「リンク関数」の意味を理解することができました。
馬場様
ハヤブサ本のころから、馬場様のご著書で学習しております。
RとStanによるモデリングのほか、Rで可視化する方法も大変参考になります。
次は「意思決定分析と予測の活用」も拝読させていただきます。
いま第4部・第1章を読んでおり、brmsでGLMMを実装する部分で質問がありコメントを書かせていただきました。
コーディングと実行自体はうまく出来たのですが、求まったパラメータの値だけを見ても「ランダム効果を取り入れた結果うまくモデリングできた」というイメージがしにくく、可視化をしたいと思っています。
そこで、通常のポアソン回帰の予測区間を可視化したコードを参考に、以下のようなコードで可視化をしてみました。
set.seed(1)
eff_glmm_pre <- conditional_effects(glmm_pois_brms,
method = "predict",
effects = "temperature:weather",
prob = 0.99)
plot(eff_glmm_pre, points = T)
すると可視化自体はできましたが、予測区間の範囲としてはGLMよりGLMMのポアソン回帰の方が狭く見え、かえってGLMMの方が予測区間の範囲に入らないサンプルが増えたようにも見えます（当初はすべてのサンプルが予測区間の中に入っているような絵ができると想像していました）。
また、釣獲尾数の事後期待値(事後中央値？)の太線についても、GLMよりGLMMの方が値が小さくなったように見えます。
この点をどのように考えたらよろしいでしょうか？
– GLMMの場合の予測(区間)の計算や、可視化の方法が間違っている？
– 上記コードではtemperatureとweatherで説明される範囲しか表示されておらず、過分散を表現したランダム効果は可視化されていないだけ？
– そもそも私の理解不足で、大きな勘違いをしている？
ご教示いただけますと幸いです。宜しくお願い致します。
Yoshiki様
コメントありがとうございます。
管理人の馬場です。
返信が遅れて失礼しました。
おそらく『conditional_effects』関数の引数の設定が間違っているものと思います。
p258で解説してありますが、ランダム効果をグラフに反映する場合は
『re_formula = NULL』の指定が必要です。
これを入れると、広い予測区間が描けるかと思います。
馬場様
ご返信ありがとうございます。
第4部を最後まで読み、理解いたしました。
また『re_formula = NULL』を指定することで、ご質問したグラフを作成することができました。
ご回答、ありがとうございました。
ところで再度の質問にて恐縮ですが、ベイズ統計モデリングの使い所についてご教示いただけますでしょうか？
階層ベイズモデルで扱うような、階層構造をしている事象のモデリングであれば、ベイズ統計の出番になると理解しています。
一方でGLMであれば、ベイズ推論でなくても（頻度統計の）glm関数でもモデリングはできてしまうと思います。
GLMでもベイズ推論でやるべきケースというのはあるのでしょうか？
実務で活用する際に参考にしたいと思います。どうぞよろしくお願いいたします。
Yoshiki様
コメントありがとうございます。
管理人の馬場です。
コメントを見逃しており、失礼しました。
ベイズモデリングの良さは、複雑なモデルでも統一的に推定できることと、予測分布を比較的容易に得ることができることだと思います。
ご質問のようなGLMといった単純なモデルでも、例えばポアソン回帰モデルで予測分布を得るのは、Rのglm関数（最尤法で推定する）だと難しいです。
brmsを使えば、気軽に95%予測区間を描いたり、予測分布を計算できます。
予測結果が確率で得られるというのはなかなか便利でして、人や場面によって許容できるリスクの大きさが異なる際に役立ちます。
馬場様
ご回答ありがとうございました（私もしばらく見逃しておりました）。
確かに私の仕事でも、予測結果が確率で得られると便利なことが多々ありそうです。
ちょうどいま読んでいる意思決定分析にもつながる利点だと思いました。
実務でも積極的に使っていきたいです。
いつも丁寧にご回答いただき、ありがとうございます。
馬場様
ベイズ統計の分析方法を使うため、本書を購入させていただきました。
Stanでハードルが高かったベイスをbrmsで簡単に分析ができる方法をさまざまな分析例ごとに紹介をしてくださり、本当に助けられています。
ありがとうございます。
brmsの推定結果のアウトプットについての質問です
第７章の正規線形モデルで、パラメーターの推定値の係数を、事後期待値(EAP）ではなく、事後確率最大値（MAP)にしたい場合に、どのように設定をすればよろしいのでしょうか？
結果を報告する際に、MAPが求められることが多いです。
本書の目的からずれてしまっているとは思いますが、もし可能なら教えていただけると幸いです。
よろしくお願いいたします。
SH様
返信が遅れて失礼いたしました。
コメントありがとうございます。
管理人の馬場です。
MAPを得ることはもちろんできるでしょうが、
下記のコメントでbrmsの作者が言うように、
個人的にもMAPより、事後中央値や事後期待値の方が安定していて使いやすいと思います。
https://discourse.mc-stan.org/t/obtaining-the-maximum-a-posteriori-estimate-of-jointly-credible-parameters-in-brms-models/14419
どうしてもMAPを得たい場合は、brm関数の設定ではおそらく達成できず、
MCMCサンプルを取得した結果に対して最頻値を求めることになると思います。
馬場様
ベイズ統計モデリングによるデータ分析入門で、勉強させていただいています。
本書に関して、質問がありメールさせていただきました。
p69の計算例において、「カーネル比を計算すると、0.001ほどの小さい値になります。そのため、θ2(提案）は採用されず、θ2=θ1となります」と書かれていますが、カーネル比が１より小さかった場合は、確率カーネル比でθ2（提案）を、確率（１ーカーネル比)でθ1を乱数として採用するのではないでしょうか？
ご教授いただければ幸いです。
河野理さま
コメントありがとうございます。
返信が遅れて失礼いたしました。
「カーネル比が１より小さかった場合は、確率カーネル比でθ2（提案）を、確率（１ーカーネル比)でθ1を乱数として採用するのではないでしょうか？」
はい。その通りです。
書籍中の計算事例では、カーネルの比がとても小さく、確率が0に近いため、本文のような書き方にしています。
MH法の理論としては、p68に記載の通りであり、河野理さまの理解で間違いありません。
参考になれば幸いです。
馬場様
ご回答いただきありがとうございました。
非常にわかりやすい書籍であったので、最後まで楽しく読ませていただきました。
馬場様
こちらの書籍で勉強させていただいております。
3点質問失礼します。
①P310 第5部第5章5節 「平滑化トレンドモデルの別の表現」にて
式(5.30)の下の式は、
μ_t = μ_t-1 + δ_t-1
となっていますね。しかし、素直に式変形をしていくと
μ_t = μ_t-1 + δ_t-1 + ζ_t
となるはずではないのでしょうか？どのように解釈をすればよろしいでしょうか？
②「平滑化トレンドモデル」という言葉はネットや他の書籍にあまり載っていないようですが、別名や英語での言葉、表記がありましたら教えていただけると幸いです。
③P314 第5部第5章11節 「ローカル線形トレンドモデルのためのStanファイルの実装」にて
「平滑化トレンドモデルを推定するためのStanコードは……」
は正しくは、
「ローカル線形トレンドモデルを推定するためのStanコードは……」
ではないでしょうか？
以上3点です。よろしくお願いいたします。
もん様
コメントありがとうございます。
管理人の馬場です。
返信が遅れて失礼いたしました。
① p309の下から1行目の「1階差分の大きさ」というδの定義から「δ_t = μ_t – μ_t-1」ですので、式変形して「μ_t = μ_t-1 + δ_t」で正しいです。
式(5.30)では、δ_tではなくδ_t-1となっていますが、1時点スライドするのみですので、変化のパターンとしての違いはありません。コードを書くときに見やすくするために1時点ずらしているだけです。
②「２次のトレンド成分モデル」と呼ばれることもあるようです。
例えば以下のPDFで使われています
https://www.actuaries.jp/lib/y_ronbun/pdf/2017.pdf
英語だと「Smooth trend」でしょうかね。例えばPythonの統計モデルのライブラリで使われています
https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.structural.UnobservedComponents.html
③こちら、誤植です。申し訳ございません。
誤植があったこと、お詫びいたします。
また、ご指摘いただけたこと、感謝いたします。
馬場様
当該本で学習させていただいております。
質問がありメールさせていただきました。
P.133のmcmc_aresを実行すると以下のエラーが出ます
—————-
expansion(add = c(0, 0.5 + 1/(2 * nlevels(data$parameter))), でエラー:
関数 “expansion” を見つけることができませんでした
—————-
こちら原因わかりますでしょうか。
よろしくお願いいたします。
yoshi様
コメントありがとうございます。
管理人の馬場です。
当方の環境では再現できませんでした。
推測ですが、ライブラリのバージョンの違いが原因かもしれません。
ライブラリを更新するなどして、再度試してみてはいかがでしょうか。
環境の違いが原因のようなので、サポートしにくいところがあります。ご容赦ください。
ありがとうございます！ライブラリ更新で解決いたしました。
本に書かれているバージョンのRとRStudioでStanをインストールしようとしましたが、うまくインストールできません。
エラーメッセージは以下の通りです。どうすればいいでしょうか？
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
> remove.packages(“rstan”)
Removing package from ‘C:/Users/260850/Documents/R/win-library/3.5’
(as ‘lib’ is unspecified)
Error in remove.packages : there is no package called ‘rstan’
> if (file.exists(“.RData”)) file.remove(“.RData”)
> install.packages(‘rstan’, repos=’https://cloud.r-project.org/’, dependencies=TRUE)
WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:
installation of package ‘rstan’ had non-zero exit status
chrono様
コメントありがとうございます。
返信が遅れてしまい、失礼いたしました。
当方の環境ではインストールできたので、以下は、あくまでもエラーメッセージを読んだ推測になります。
先に、下記コードを実行してRToolsをインストールしてみてはどうでしょうか。
# Rtoolsのインストール
pkgbuild::has_build_tools(debug = TRUE)
また、本記事の「環境構築の手順」や「エラーが出た時の対応」も参照してみてください。
デザイン行列を使ったポアソン回帰モデルのためのStanファイル「3-8-3-glm-pois-design-matrix.stan」(p218)を実行するためのプログラムをご教示ください。
因みに、必要なパッケージをインストールし、以下のプログラムを実行しましたが、上手くいきません。お忙しい所恐縮ですが、宜しくお願いいたします。
・・・・・・・・
mutate(fish_num_climate, weather =ifelse(weather==”sunny”, 1, 0))
formula_lm <- formula(fish_num ~ temperature + weather)
X <- model.matrix(formula_lm, fish_num_climate)
N <- nrow(fish_num_climate)
K <- 3
Y <- fish_num_climate$fish_num
data_list__design <- list (N=N, K=K, Y=Y, X=X)
mcmc_result_design <- stan(
file = "3-8-3-glm-pois-design-matrix.stan",
data = fish_num_climate,
seed = 1
)
mcmc_result_design
康青さま
コメントありがとうございます。
管理人の馬場です。
コメントを見落としており、大変失礼いたしました。
下記のGitHubから実行コードがダウンロードできますので参考になさってください。
https://github.com/logics-of-blue/book-r-stan-bayesian-model-intro/blob/master/book-data/3-8-%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E5%9B%9E%E5%B8%B0%E3%83%A2%E3%83%87%E3%83%AB.R
p.188においてas.mcmcを利用すると
as.mcmc.brmsfit is deprecated and will eventually be removed.
という警告文が出現します。こちら現在推奨の関数はわかりますでしょうか？
突然の質問失礼いたします。
そろそろベイズ統計も勉強してみたいなと思っているなかで、ほぼ初心者でこの本を引き当てました。大当たりだと思っております。
構成もよく、勉強しやすいと感じます。
ｐ249-250の内容で実行すると急に
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
0
Syntax error in ‘string’, line 11, column 27 to column 28, parsing error:
Only top-level variable declarations allowed in data and parameters blocks.
というエラーが出るようになってしまいました。
お忙しい中申し訳ありませんが解決法のアドバイスをいただけないでしょうか？
よろしくお願いいたします。
R version 4.2.1
Platform: x86_64-w64-mingw32/x64 (64-bit)
rstan version 2.26.13
西嶋さま
コメントありがとうございます。
管理人の馬場です。
拙著をお読みいただきありがとうございます。
stanファイルそのものが間違っているというエラーのようですが、当方の環境では再現できませんでした。
stanファイルはGitHubからダウンロードしたものをそのまま使っているでしょうか？
おそらくGitHubからstanファイルをダウンロードして一切変更せずに実行したら、正しく動くのではないかと思います。
あるいは、Rstanのバージョンを2.19.3に下げて実行したらうまくいくかもしれません（当方の環境がそのようになっているので）。
ご検討ください。
馬場先生
早速の返信ありがとうございます。落ち着いて見直すとstanファイルのメモで”/”を”\”で記入していました。
すみません、お恥ずかしいミスです。本書とGitHubをしっかり見て入力したつもりでしたが、メモ欄でなくコードばかり見ていたため、なんでや～とずっとなっていました。
恥かきついでにもう一つ質問をしてもよろしいでしょうか？
ｐ286のコードで図示すると結果の横軸が2010年01月でなく2009年12月スタートになってしまいます。
自身で模倣したコードだけでなく、Githubからコピペしても同様にずれてしまいます。
分析データを確認しても2010年1月1日が最初のデータになっています。
グラフの端ではなく、グラフ点自体が始まっている箇所が2009年12月となってしまっています。
申し訳ありませんが、何か解決法をご教授いただければ幸いです。
Windowsにて最新版のRとRStanをインストールすると、RStanを正常に実行できないようです。ブログに追記事項として掲載頂けませんでしょうか。
参考サイトは以下の通りです。
https://developer.mamezou-tech.com/blogs/2022/06/30/install-rstan-on-r421/
https://blog.mc-stan.org/2022/04/26/stan-r-4-2-on-windows/
コメントありがとうございます。
管理人の馬場です。
情報ありがとうございます。
上記のリンクを当ブログにも記載しておきます。