1何のための手法か

「人口が増えた集落ほど方言語彙が失われたか」を、係数の事後分布として問う。

ベイズ推定の一般的な定式化は7つの観点を式で読むの観点3にある。このページはその繰り返しではなく、実際に階層モデルを推定するときに何を決めなければならないかを扱う。決めるべきものは、尤度（応答の分布とリンク関数）、事前分布、サンプラーの設定、そして収束と適合をどう確かめるかの4つである。

この研究で欲しかったのは、標準化した人口成長率の係数 \(\beta\) について「\(\beta > 0\) である事後確率」という形の答えだった。応答は語彙の消失率で、厳密に \((0,1)\) に収まる割合である。ガウス分布を当てるのは分布の誤指定になるため、ベータ分布に logit リンクを置いた（Ferrari & Cribari-Neto 2004）。集落は中学校区に入れ子になっているので、区ごとの切片をランダム効果とした。

標本は25集落しかない。小標本では事前分布が結果に効く。だから事前分布を既定値まかせにせず明示し、その選択が何を意味するのかを言えるようにしておく必要がある。

2起源 ―― 原理は18世紀、道具は2017年

定理が使えるようになるまでに250年かかっている。

Stan と brms が変えたのは原理ではなく参入障壁である

Stan 以前、ベイズ階層モデルを使うにはサンプリングのコードを自分で書く必要があり、事実上、統計の専門家しか使えなかった。brms は lm() とほぼ同じ書き方でモデルを指定できるようにした。このページで扱う実装も brm(loss ~ ... , family = Beta(link="logit")) という一行で書かれている。応用研究者が直接この手法を使えるようになったのは、ごく最近のことである。

3数学的にどこから来ているか

確率論の公理、マルコフ連鎖、そして19世紀の解析力学。

ベイズ推定そのものは Kolmogorov (1933) の確率の公理化の上に乗っている。「\(P(\beta > 0) = 0.99\)」と書けるのは、確率が集合に割り当てられた測度として厳密に定義されているからである。この基盤は頻度論の手続きとも共有されており、ベイズだけのものではない。

計算のほうは別の系譜から来ている。事後分布の分母は積分であり、階層モデルでは解析的に解けない。そこで分布そのものを求めず、そこから標本を取る方針に切り替える。標本を取る仕組みがマルコフ連鎖である（Markov 1906）。現在の状態が決まれば過去は未来に無関係、という性質を使って、事後分布を定常分布とする連鎖を作り、そこを歩かせる。

求めたいものと、実際に計算するもの
\[ p(\theta \mid D) \;=\; \frac{p(D \mid \theta)\,p(\theta)}{\displaystyle\int p(D\mid\theta)\,p(\theta)\,d\theta} \qquad\longrightarrow\qquad \theta^{(1)}, \theta^{(2)}, \dots, \theta^{(S)} \;\sim\; p(\theta \mid D) \]

左辺の分母は階層モデルでは解けない。右のように事後分布からの標本列が得られれば、事後平均も区間も \(P(\beta>0)\) も、標本から数え上げるだけで出る。MCMC は積分を標本で置き換える手続きである。

Stan が使うのはハミルトン力学である

素朴な MCMC は行き当たりばったりのランダムウォークで、次元が上がると効率が落ちる。Stan が採るハミルトニアン・モンテカルロは、事後分布の対数を位置エネルギーの地形と見なし、Hamilton (1833) の運動方程式に従ってその上でボールを転がす。勾配の情報を使うぶん、高次元でも遠くまで一度に移動できる。刻み幅と経路長を自動で決める No-U-Turn Sampler（Hoffman & Gelman 2014）が既定である。19世紀の解析力学が、21世紀のベイズ計算の中核にいる。

4定義と仮定 ―― 何を決めているのか

尤度・事前分布・サンプラー設定の3つを、すべて分析者が決める。

この研究のモデル
\[ \mathrm{loss}_i \;\sim\; \mathrm{Beta}\bigl(\mu_i \phi,\; (1-\mu_i)\phi\bigr), \qquad \mathrm{logit}(\mu_i) \;=\; \alpha + u_{d(i)} + \beta\, x_i + \gamma_1 z_{1i} + \gamma_2 z_{2i} \]

\(\mu_i\) … 集落 \(i\) の平均消失率。logit リンクで実数直線に写す
\(\phi\) … 精度パラメータ。ベータ分布の散らばりを決める
\(u_{d(i)}\) … 集落 \(i\) が属する中学校区 \(d\) のランダム切片（5区）
\(x_i\) … 標準化した人口成長率。\(z_1, z_2\) … 標準化した対数標本サイズ（中学生・高年層）

弱情報事前分布を選ぶ理由

事前分布には3つの選択肢がある。強い事前分布は結論を先に決めてしまう。無情報に近い平坦な事前分布は、小標本ではサンプラーが不安定になりやすく、logit スケールでは現実にはありえない大きさの係数にまで確率を配ってしまう。弱情報事前分布はその中間で、「おおよそこの範囲」という薄いガイドだけを与える。

実際に置いた事前分布
\[ \beta,\gamma \sim \mathcal{N}(0,\,1), \qquad \alpha \sim \mathcal{N}(0,\,2), \qquad \phi \sim \mathrm{Gamma}(0.1,\,0.1) \]

係数に \(\mathcal{N}(0,1)\) を置くことは、logit スケールで「効果はおおむね \(-3\) から \(+3\) の範囲だろう」と述べるのに等しい。標準化した予測変数に対しては十分に広く、かつ発散を防ぐ程度には狭い。ゼロを中心に置いているので、この事前分布は効果を大きく見せる方向には働かない。

事前分布は「影響がない」のではなく「影響が小さいことを確かめる」ものである

この研究では、データが事前分布を上回る情報量を持つため最終結果への影響はほとんど出ない。ただしそれは結果として確かめられたことであって、事前分布を置けば自動的にそうなるわけではない。n=25 という標本では、より狭い事前分布を置けば推定は確実にそちらへ引かれる。事前分布を明示して報告するのは、読み手がこの点を検証できるようにするためである。

5この研究での使い方

4本の連鎖を回し、収束を確かめ、事後分布から確率を数える。

  1. データ　消失率が得られた25自治会。中学校区は5区。応答 loss は厳密に \((0,1)\) 内にあり、端点変換は不要だった。
  2. モデル指定　brm(loss ~ z_pop_growth + z_log_ss_ms + z_log_ss_gp + (1|district), family = Beta(link="logit"))。標本努力の交絡を統制するため、両世代の対数標本サイズを共変量に入れている。
  3. サンプラー設定　4連鎖・各4,000反復・うち2,000をウォームアップ。ウォームアップ後の標本は合計8,000。adapt_delta = 0.97、max_treedepth = 12、seed = 1。バックエンドは cmdstanr。
  4. 収束の確認　\(\hat{R}\) と有効標本サイズ、発散遷移の数を見る。報告したモデルは \(\hat{R} = 1.00\)、bulk-ESS > 3,000、発散は1以下だった。
  5. 事後確率を数える　8,000本の標本のうち \(\beta > 0\) であるものの割合を数える。それが \(P(\beta > 0)\) である。
  6. 頑健性の確認　共変量を足したモデル、空間ランダム効果を入れたモデル、集落内の不確実性を伝播させたカウントレベルモデル、1集落ずつ除く感度分析、単位を統合した再解析、検出力シミュレーションを併走させた。
モデル\(\beta\)95% 確信区間\(P(\beta>0)\)
base（主モデル）+0.231[+0.06, +0.40]0.993
cov（面積・標高・建物用地を追加）+0.229[+0.03, +0.42]0.987
sp（ICAR 空間ランダム効果を追加）+0.228[+0.04, +0.40]0.992

標準化した人口成長率の係数。数値は R/12_brms_loss.R の実行ログからの引用。主モデルは \(\hat{R} = 1.00\)、Bulk-ESS 5,184、Tail-ESS 5,006。空間モデルの詳細は ICAR 空間ランダム効果 のページで扱う。

確信区間と信頼区間は、確率をどこに置くかが違う

上の区間は事後分布の 2.5 パーセンタイルと 97.5 パーセンタイルである。「この区間に真の値がある確率が95%」と読んでよい。確率が置かれているのはパラメータのほうである。頻度論の信頼区間は逆で、パラメータは固定した値と見なし、確率は手続きのほうに置く。「同じ手続きを繰り返すと、作られる区間の95%が真の値を含む」という主張であり、いま手元にある一つの区間について確率を述べてはいない。

この違いは推定原理の違いではない。尤度はベイズでも使う。上の式の \(p(D\mid\theta)\) がそれで、最尤法はその尤度を最大にする点を選ぶ原理のひとつにすぎない。ベイズと頻度論を分けているのは、尤度を使うかどうかではなく、パラメータに分布を認めるかどうかである。

系統別の部分プーリング

語彙を関西系・美濃系・東京系の3系統に分け、系統ごとの成長 slope を部分プーリングで推定した。関西 \(\beta = 0.233\) [0.073, 0.394]、美濃 \(\beta = 0.240\) [0.083, 0.403]、東京 \(\beta = 0.219\) [0.022, 0.402]。事後確率はそれぞれ 0.997、0.996、0.983 である。slope がほぼ同一なのは強い縮約が効いているためで、系統間で効果が似ていることを反映している。

6限界 ―― この手法で言えないこと

事後確率は、モデルが正しいという条件のもとでの確率である。

収束していない連鎖から出た数値は、いかなる意味も持たない

\(\hat{R}\) は複数の連鎖が同じ分布を探索できているかを測る診断量で、収束していれば 1 に近づく。Vehtari ほか (2021) は従来の \(\hat{R}\) が裾の重い分布や連鎖間で分散が違う場合に収束failureを見落とすことを示し、順位正規化した改良版を提案している。現在の brms が報告するのはこの改良版である。

同じ研究の多様性モデルでは、これが実際に問題になった。非空間モデルの一部で \(\hat{R} = 1.24\) に達し、未収束と判断して結果を採用していない。空間モデルでは max_treedepth に毎回到達し、空間分散が0付近に貼りついてファネル状の事後分布になっていた。いずれも「数値が出た」ことと「推定できた」ことが別である例である。事後分布の要約を読む前に、まず診断を見る必要がある。

「事後確率0.99」は「頑健」を意味しない

この研究の主結果は \(P(\beta>0) = 0.99\) だが、著者はこれを頑健な結論としては提示していない。集落内の標本不確実性を伝播させるカウントレベルモデルでは、対応する交互作用は \(\beta = -0.17\)、95%確信区間 [\(-0.50\), \(0.18\)] と0をまたぐ。1集落ずつ除く感度分析では、高年層が1名しかいない単一集落を除いたときだけ \(\beta = 0.12\) まで落ちる。高年層が5人以上になるまで集落を統合して再解析すると効果は消える。

検出力シミュレーションがこの食い違いを説明する。実際の標本設計のもとで、真の交互作用効果 \(-0.15\) と \(-0.30\) を検出できる確率はそれぞれ 12% と 39% にすぎない（\(-0.50\) でようやく82%）。集落レベルの解析が強い結論を返すのは、各集落の消失率を誤差のない精密な量として扱っているからである。事後確率が高いことは、モデルが置いた仮定のもとでの確からしさを示すにとどまる。

検定と同じ誤用は、ベイズでも起きる

「区間が0を含まない」を有意性の代わりに使えば、名前が変わっただけで同じことをしている。この研究では、効果がないことを主張したい箇所について、区間に加えて \(\pm 0.05\)・\(\pm 0.10\) の実質的等価領域に事後確率がどれだけ入るかを報告している。「効果がない」と断定するのではなく、どの程度の精度で言えているのかを述べるための措置である。

7出典とデータ

対象は滋賀県米原市から岐阜県関ケ原町にかけての集落で行われた3世代方言語彙調査である。原調査は2015年にはびろネットが実施し、本研究はその匿名化データの二次解析にあたる。解析は R/01 から R/25 までのスクリプトに分けて記述され、原データファイルは改変せず作業コピー経由でのみ読む構成になっている。論文は準備中で、本稿執筆時点では公開されていない。

手法の文献

  • Bürkner, P.-C. (2017) brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1–28. doi:10.18637/jss.v080.i01
  • Carpenter, B., Gelman, A., Hoffman, M. D. ほか (2017) Stan: A probabilistic programming language. Journal of Statistical Software, 76(1), 1–32. doi:10.18637/jss.v076.i01
  • Hoffman, M. D. & Gelman, A. (2014) The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15, 1593–1623.
  • Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. & Bürkner, P.-C. (2021) Rank-normalization, folding, and localization: an improved \(\hat{R}\) for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718. doi:10.1214/20-BA1221
  • Ferrari, S. & Cribari-Neto, F. (2004) Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815.

この手法がどの観点に属するか

ベイズ推定は、パラメータのほうに分布を認め、観測前の知識とデータを合わせて事後分布に更新する。返すのは点推定ではなく分布であり、「効果が正である確率」という形の答えを直接与える。7つの観点を式で読むの観点3（ベイズ）にそのまま重なる。同じ研究で使っているICAR 空間ランダム効果も、空間構造を事前分布として書くという意味で同じ観点に属する。

このページは、小野原が実際の研究で使っている手法を、同じ分野の読み手に向けて整理したものです。数値はすべて解析スクリプトの実行ログと出力ファイルから引用し、推定値や概算は記載していません。
手法の位置づけの全体像は 7つの観点を式で読む、どの授業回がどの観点に当たるかは 手法と数学の地図 にあります。