1何のための手法か

「隣の集落と似ているだけではないか」という対抗仮説を、モデルの中で扱う。

地理上に並んだ地点のデータには、ほとんどの場合、空間的な構造がある。隣り合う集落は人の行き来があり、同じ学校区に属し、似た歴史を共有する。したがって各集落の観測値が互いに独立だという前提は、そのままでは成り立たない。独立を仮定した回帰は、この重複した情報を新しい情報として数えてしまい、係数の不確実性を小さく見積もる。

対処は2通りある。ひとつは空間構造を取り除いてから回帰する方法、もうひとつは空間構造をモデルの一部として明示的に推定する方法である。ICAR は後者にあたる。各地点に空間ランダム効果 \(\phi_i\) を置き、その \(\phi\) に「隣どうしは近い値を取る」という事前分布を与える。

この研究で問いたかったのは、人口成長率と語彙消失の関連が空間的な近さで説明できてしまわないかだった。近い集落は人口動態も似ているので、空間効果を入れても係数が残るかどうかが確認の要点になる。

2起源 ―― ガンの地図から来た手法

疫学が、地域別の発生率を地図にするために作った。

Julian Besag（1945–2010、英・統計学者）が York と Mollié とともに 1991 年に発表したモデルが原型である。もともとの目的は地域別のガン発生率の地図を作ることだった。地域ごとに発生率を素朴に計算すると、人口の少ない地域では分母が小さく、たまたま数人多いだけで発生率が跳ね上がる。地図は実態ではなく標本のばらつきを描いてしまう。

そこで「隣の地域と発生率が似るのは当然だ」という知識を事前分布として持ち込み、各地域の推定を隣接地域のほうへ引き寄せる。人口の少ない地域ほど強く引き寄せられ、地図は安定する。空間的自己相関を邪魔者として消すのではなく、推定を助ける情報として使うのがこのモデルの発想である。

疫学から言語研究へ

「都市部に近い農村は都市部に似た発生率を示す」という疫学の構図は、「隣の集落の語彙変化は似ている」という言語地理学の構図にそのまま移せる。この研究が生態学の希薄化や疫学の空間モデルを使っているのは偶然ではなく、標本の大きさが地点ごとに違う集計データという問題の形が共通しているためである。

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

空間版のマルコフ性。グラフの上の確率場として書く。

ICAR の数学的な根はマルコフ連鎖と同じところにある（Markov 1906）。時間の連鎖におけるマルコフ性は「現在が決まれば過去は未来に無関係」だった。これを空間に移すと、「隣接する地点の値が決まれば、それより遠い地点は無関係」になる。この性質を持つ確率場をマルコフ確率場と呼び、Hammersley–Clifford の定理が、条件付き分布の指定から同時分布が決まる条件を与えている。

つまり ICAR は、同時分布を直接書く代わりに各地点の条件付き分布だけを書くモデルである。全体像を指定せず、局所的な関係だけを指定する点が特徴である。

4定義と仮定

条件付き自己回帰（CAR）と、その極限としての ICAR。

CAR ―― 各地点の条件付き分布で定義する
\[ \phi_i \mid \phi_{-i} \;\sim\; \mathcal{N}\!\left(\alpha\,\frac{\displaystyle\sum_{j} w_{ij}\,\phi_j}{\displaystyle\sum_{j} w_{ij}},\;\; \frac{\sigma^{2}}{\displaystyle\sum_{j} w_{ij}}\right) \]

\(w_{ij}\) … 隣接行列の成分。隣なら 1、そうでなければ 0
\(\phi_{-i}\) … 自分以外のすべての地点の空間効果
平均 … 隣接地点の平均に \(\alpha\) を掛けたもの。自分の値は隣の平均へ引き寄せられる
分散 … 隣接数が多い地点ほど小さい。情報が多いほど強く決まる

ICAR ―― \(\alpha = 1\) とした極限
\[ \phi_i \mid \phi_{-i} \;\sim\; \mathcal{N}\!\left(\bar{\phi}_{\partial i},\;\; \frac{\sigma^{2}}{n_i}\right), \qquad p(\boldsymbol{\phi}) \;\propto\; \exp\!\left[-\frac{1}{2\sigma^{2}}\sum_{i \sim j}\bigl(\phi_i - \phi_j\bigr)^{2}\right] \]

\(\bar{\phi}_{\partial i}\) … 地点 \(i\) の隣接地点における \(\phi\) の平均。\(n_i\) … 地点 \(i\) の隣接数
\(i \sim j\) … 隣接する組についての和。罰則の形は「隣どうしの差の二乗和」であり、隣接する地点の値が離れているほど確率が下がる

ICAR は事前分布として非正則である

右の式をよく見ると、\(\boldsymbol{\phi}\) に同じ定数を足しても差 \(\phi_i - \phi_j\) は変わらない。つまりこの密度は全体の水準をまったく決めていない。定数方向に平坦なので積分は発散し、確率分布として正規化できない。これが「非正則（improper）」ということである。

実害は2つある。第一に、そのままでは事前分布として不完全である。第二に、モデルに切片がある場合、切片と空間効果の水準が区別できず識別性を失う。標準的な対処は \(\sum_i \phi_i = 0\) という和ゼロ制約を課すことで、これにより水準を切片側に押しつけ、空間効果には地点間の差だけを担わせる。brms の car(type = "icar") はこの制約つきの実装である。

隣接行列は、データではなく分析者が決める

上の式に現れる \(w_{ij}\) は、どこからも自動的には出てこない。何をもって「隣」とするかは設計判断である。一定距離以内を隣とするのか、最近傍 \(k\) 個を隣とするのか、行政界を共有する地域を隣とするのか、ドロネー三角形分割で決めるのか。選び方によって隣接数の分布が変わり、引き寄せの強さが変わり、推定される空間効果が変わる。この選択はデータからは決まらない。

5この研究での使い方

3km 以内を隣と定義し、空間項の有無を比較した。

  1. 距離の計算　各集落の緯度経度から大円距離（ハヴァサイン公式、地球半径 6371 km）を計算する。
  2. 隣接行列の構成　距離が 0 より大きく 3 km 以下の組を隣接（1）、それ以外を 0 とする二値行列 \(W\) を作る。しきい値の 3 km は、先行する2016年の解析に倣った値である。
  3. 隣接構造の確認　得られた \(W\) は 25×25、平均隣接数 4.80、孤立点 0。孤立点があると ICAR は不安定になるので、この確認は必須である。
  4. モデルへの投入　car(W, gr = settlement, type = "icar") を階層ベータ回帰に加える。推定の枠組みは ベイズ推定 のページと同じで、4連鎖・8,000標本、adapt_delta = 0.97。
  5. 空間項の要否を判定　空間ありとなしを LOO 近似交差検証で比較し、空間項がデータによって支持されるかを確かめる。
比較elpd_looelpd_diffse_diff判定
消失率・空間なし（base）33.80.00.0こちらが選ばれる
消失率・ICAR あり（sp）33.1−0.70.3空間項は適合を改善しない

数値は R/12_brms_loss.R の実行ログからの引用。空間項を入れても人口成長率の係数はほぼ動かず、\(\beta = +0.228\)、95%確信区間 [+0.04, +0.40]、\(P(\beta>0) = 0.992\) だった（空間なしは \(+0.231\)）。

空間分散はほぼゼロに縮退した

同じ研究の多様性モデル（51自治会）では空間分散 sdcar そのものを読むことができる。関西系で 0.04、95%確信区間 [0.00, 0.16]。総形質数で 0.02、[0.00, 0.09]。いずれも下限が 0 に接しており、空間信号が検出されていない。この場合、事後分布は空間分散が 0 に近づくほど個々の \(\phi_i\) が縮む漏斗状の形になり、サンプラーが苦しむ。実際このモデルでは max_treedepth に毎回到達し、有効標本サイズも 179 や 57 と低く、\(\hat{R}\) は 1.02 と 1.07 だった。

結論：主モデルには入れず、頑健性確認として使う

以上から、この研究では空間項を主モデルに入れていない。ICAR は「空間構造を補正した」ためではなく、「空間構造では説明できないことを確かめる」ために使われている。25地点という標本に25個の空間ランダム効果を載せるのは過剰であり、LOO もそれを支持していない。

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

隣接行列の定義を変えると、結論が変わりうる

これは一般論ではなく、この研究群で実際に起きたことである。2016年の解析では、距離しきい値 1/2/3 km で隣接を定義し、Moran's I と Geary's C によって得失率にも回帰残差にも有意な空間自己相関は検出されなかった。ところが2017–2018年の解析では、ドロネー三角形分割や \(k\) 近傍など複数の隣接行列を定義し直した結果、空間自己相関があることを示す指標も得られている。

作業記録はこの食い違いを「要再確認の重要点」として残しており、どちらが正しいかはまだ決着していない。ここで述べられるのは、同じデータに対して隣接の定義を変えると空間自己相関の有無について異なる指標が出た、という事実そのものである。したがって本ページの 3 km という値も、ひとつの設計判断として読む必要がある。しきい値 1 km・2 km での感度分析は、作業記録に残された未了の課題である。

空間効果が「ない」ことの証明にはならない

LOO が空間なしモデルを選び、sdcar が 0 に貼りついたとしても、それはこの隣接定義・この標本サイズのもとで空間信号を検出できなかったということである。25地点という標本では検出力そのものが低い。2016年の解析でも「Moran は『なし』だが小標本で検出力が低い」という留保が明記されており、ベイズ CAR による再検証はもともとその留保に応えるために導入されたものだった。

加えて、空間モデルの LOO 計算では Pareto \(k\) が 0.7 を超える観測が1つ（4.0%）現れており、相対有効標本サイズも 0.1–0.3 と低い。比較の土台自体が完全には安定していない。

空間的な近さは、近さの一種にすぎない

ICAR が表現するのは地理的な隣接だけである。方言の伝播は道路網・通学区・通婚圏・地形の障壁に沿って進むのであって、直線距離だけで決まるわけではない。3 km の円は、そうした経路を距離で代理させたものにすぎない。この点は、同じ問題を別の手法で扱ったMantel 検定のページで「空間距離が測定値ではなく設計物である」と述べたことと同じ構図である。

7出典とデータ

対象は滋賀県米原市から岐阜県関ケ原町にかけての集落で行われた3世代方言語彙調査である。原調査は2015年にはびろネットが実施し、本研究はその匿名化データの二次解析にあたる。隣接行列の構成と空間モデルの推定は R/03_brms_car_template.R、R/09_brms_diversity.R、R/12_brms_loss.R に記述されている。論文は準備中で、本稿執筆時点では公開されていない。

手法の文献

  • Besag, J., York, J. & Mollié, A. (1991) Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43(1), 1–20. doi:10.1007/BF00116466
  • 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
  • Ferrari, S. & Cribari-Neto, F. (2004) Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815.

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

ICAR は、空間的な近さについての知識を事前分布の形でモデルに書き込む。隣どうしの差を罰する項がそのまま事前分布の対数密度であり、推定はデータと事前分布を合わせた事後分布から行う。7つの観点を式で読むの観点3（ベイズ）に属する。同じ観点にあるベイズ推定のページが推定の枠組み全体を扱い、このページはその中の空間構造の書き方を扱っている。

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