こんにちは，株式会社Nospareの小林です．
Spatial Statistics に掲載されたMacNab (2022)によるベイズ疾病マッピングに関するレビュー論文``Bayesian disease mapping: Past, present, and future''の紹介を２回くらいに渡って行います（オープンアクセスです）．疾病マッピング（disease mapping）は日本語では他にも疾病地図（なんかダサいですがこちらのほうが一般的なようです）などとも呼ばれ，病気の発生や死亡率を地図上に示すことで，それらの地理上の分布を可視化することができる統計学的な手法のひとつです（下図は2017-2019年の35歳以上の成人男性の郡ごとの心臓病の死亡率，出典：CDC Division for Heart Disease and Stroke Prevention）．本記事では，疾病マッピングとはなにか，1980-1990年代での疾病マッピングの流れ，多次元疾病マッピング，条件付き分布に基づいたガウス・マルコフ確率場モデルについて紹介します．
はじめに
例えば，ある地域において観測された病気やがんなどの発生・死亡率を地図上に示すことで病気と環境との関連性を説明する，という試みが早く（19世紀後半）から行われていました．ところが，がんといった比較的まれな病気については，観測された死亡率や相対リスクをそのまま用いると病気の地理上の分布としては安定性に欠けるものでした．病気の発生率や死亡率が高かったり低かったりする地域（郡など）を統計的に把握しようとする試みから，格子データ（lattice data）に対する空間統計の枠組みとして疾病マッピングが生まれました．疾病マッピングでは統計モデルに基づいて，リスク地図作成の安定化を目指します．疾病マッピングの初期の研究では，経験ベイズ（EB）法による推定によって死亡率や相対リスクの地図の安定化を目指していましたが，相対リスクに関する統計的推論という点が欠けていました．その後，ベイズの枠組みで疾病マッピングに関する研究が進んでいき，初期の１変量の疾病マッピングのための階層モデルの開発から，空間・時間・多変量・配列などのデータための多次元の疾病マッピングへと発展することになりました．
1980年代〜1990年代
1980年代に入って，がん登録（cancer registry，日本でもがん登録に関する法律があります）によって，行政区域レベルなどでのデータが入手できるようになったことから，統計学の分野では疾病マッピングが注目されるようになりました．上述のように，研究の目標は，安定したがんのリスクマップを作成するためのモデルベースの手法を開発することでした．これにより，様々ながんの環境的決定要因や病因などに関して仮説を立て，さらなる調査を行うことができると考えられました．
モデルの説明を行うために簡単な設定のもと表記を導入します．$O=(O_1,\dots,O_n)$を地域レベルの観測された死亡数（実死亡数）とします．$n$は地域数で，「地域」とは郡，区などの行政区域や地方などに当たります（研究によって異なります）．一方で，$E=(E_1,\dots,E_n)$を期待死亡数とします．期待死亡数は多くの場合，地域別・年齢別の暴露人口とそれに対応する年齢別の参照率から計算されます．相対リスク$RR_i=O_i/E_i$は標準化死亡率（SMR）で定義されることが多いです．
がんの死亡などといったまれな事象を取り扱う際，実死亡数はポアソン分布に従っていると見ることができます：
O_i\sim Poisson(\phi_i),\quad \phi_i=E_i\varphi_i
ここで$\phi_i$はポアソン分布の平均で，$\varphi_i$は未知の相対リスクです．死亡数は地理上で独立なポアソン事象だと仮定すると，相対リスクの最尤推定量は$\hat{\varphi}_i=O_i/E_i$，標準誤差は$\text{se}(\hat{\varphi}_i)=\sqrt{\hat{\varphi}_i)}/E_i$となります．この最尤推定量は小さな死亡数のもとでは安定性と効率性に欠けることが知られています．また空間上の相関から，独立性の仮定は制約的と言えます．さらに，ポアソン分布の期待値と分散が等しいという性質も，疾病マッピングでよく用いられるデータでは現実的ではないことが多いです（過分散：分散のほうが期待値よりも大きい）．
経験ベイズ（EB）のアプローチでは，$\varphi$に事前分布$f(\varphi|\theta)$を仮定し，ハイパーパラメータ$\theta$はデータから推定します．1980年代はたくさんのEBアプローチに関する研究が行われました．
リスクは非負なので，例えば$\varphi$の事前分布としてガンマ分布$\varphi_i\sim \text{gamma}(u,v)$を仮定することができます．これは空間事前分布ではない単純な事前分布なので，$O_i$の周辺分布は負の二項分布になり，過分散の分布として知られています．$\varphi_i$のEB推定量は事後平均に当たる$E[\varphi_i|O_i]=\frac{O_i+v}{E_i+u}$ですが，負の二項分布に基づいてハイパーパラメータ$u$と$v$を最尤推定することによって，
\tilde{\varphi}_i=\frac{O_i+\hat{v}}{E_i+\hat{u}}
とします．ここで，$\tilde{\varphi}_i$は，標準化死亡率$O_i/E_i$と$\varphi_i$の平均の推定値$\hat{v}/\hat{u}$の折衷となっていることがわかります（グローバルスムージングと呼びます）．
グローバルスムージングはデータに空間相関がある場合には望ましくはありません．空間相関を考慮し空間スムージングを行うために，条件付き分布によって定式化されたガウス・マルコフ確率場（Gauss-Markov random field; GMRF）に基づく事前分布（conditional autoregressive prior; CAR）を$\psi=(\psi_1,\dots,\psi_n)=(\log\varphi_1,\dots,\log\varphi_n)$に仮定することが提案されました：
\psi\sim GMRF(\mu,\Omega),\quad \Omega=\sigma^{-2}(I_n-\lambda W)
ここで$\Omega$は精度行列，$\sigma$は尺度パラメータ，$\lambda$は空間依存の度合いを決定するパラメータ，$W$は隣接行列です．CARモデルによって，対称的な空間相関，空間スムージング，リスクの予測や推測において情報を借り合うことができるようになります．このCARモデルの平均と精度はそれぞれ
E[\psi_i|\psi_{-i}]=\lambda\sum_{j\sim i}\psi_j,\quad \text{Prec}(\psi_i|\psi_{-i})=\sigma^{-2}
で与えられ，$\psi_{-i}$は$\psi$から$\psi_i$を除いたもの，$j\sim i$は地域$j$と$i$が隣接していることを意味します．地域別の対数リスクの条件付き期待値は隣接する地域の対数リスクの和に比例し，$\lambda$によってその度合が決定されます．ポアソン-ガンマモデルと違って，このモデルは$O$の解析的な分布を求めることができないのでEMアルゴリズムなどによって$\psi$とハイパーパラメータの推定を行います．
（ナイーブな）EBによる$\varphi$の短所としては，点推定値はほぼ不偏なのに対し，標準誤差は推定による不確実性を過小評価してしまうことにあります．これは，EB推定量にハイパーパラメータの推定値を代入することによって，推定の不確実性が十分に考慮されないということに由来します．1990年代，2000年代ではEB推定に対する改善が行われてきました．1990年代には，マルコフ連鎖モンテカルロ法（MCMC）に関する研究と実用化が急速に進んだおかげで，疾病マッピングの推定をフルベイズで行えるようになり，現在ではフルベイズがメインストリームとなっています．RではCARBayesやR-INLAなどといったパッケージがあり，これを使って空間モデルに対するベイズ推測を行うことができます．
多次元の疾病マッピング
多次元の疾病マッピングは非常に幅広いトピックで，時空間，多変量，配列などに対するモデリングが含まれます．
ここでは対数相対リスクの行列を$\psi=(\psi_{ij})$として表記することにします．時空間疾病マッピングの場合は$j=1,\dots,T$というように$T$時点を考え，多変量の場合は$j=1,\dots,J$というように$J$種類の病気を考えることになります．このようなモデルは地域間，時間上で情報を借り合ったりリスクのスムージングを行うことが可能になり，空間リスクモデルの拡張と見ることができます．
時間と空間の交差効果を考えることも可能で，一般的な特定化として次のようなモデルが挙げられます：
\psi_{it}=\alpha_t + \varphi_t + \phi_i + \eta_i + \delta_{it},\quad i=1,\dots,n, \quad t=1,\dots,T.
ここで$\delta_{it}$は時空間の交差効果を表します（添字$t$と$i$がついている変数が２つずつ入っているのは，それぞれ片方が相関のある空間あるいは時間の効果を表す項，もう片方が独立な空間・時間の効果を表す項とするのが多いです）．もっと単純に
\psi_{it}=\alpha_t + \theta_i + \delta_{it}
とすることもできます．これらの空間・時間効果は自己回帰（autoregressive; AR）モデルやCARモデルなどでモデリングすることができます．
また，多次元配列の疾病マッピングのためにこのモデルを拡張することも可能です．例えば，空間・時間効果に加えて，年齢効果を加え，この年齢効果は例えばスプラインなどでモデリングできます．
多変量疾病マッピングでは共通のリスク要因を持つような複数の病気を同時にモデリングすることで，情報を効果的に借り合うことができ空間スムージングを改善させることができます．多変量モデルでは，例えば
\psi_{ij}=\theta_i+\delta_{ij}
というようなリスクを２つの項の和としてモデリングしたりします．もちろん，より柔軟なアプローチとして多変量GMRF事前分布を$\text{vec}(\psi)$に仮定することもできます．
地域レベルの共変量を用いる
疾病マッピングにおいて，共変量を加えて回帰モデルを考えることは実務上重要であることが多いです．空間回帰モデルでは病気リスクの変動を，リスク要因によって説明できる部分と説明できないに分けて定量化することができます．このような回帰モデルはecological regression（ER）と呼ばれ，環境暴露，社会経済的な貧困，その他地域レベルの特徴がどのように病気のリスクに影響を与えるかを特定することができます．
一方でERはいくつかの課題が挙げられます．例えば，リスク暴露と病気の関連の度合いは地域レベルと個人レベルでは異なるため，バイアスが生じる可能性があります．また空間GLMMにおける固定効果と変動効果の共線性や，観測されない交絡因子，共変量の地理上の観測単位のずれなどの問題もあります．
条件付き分布で構成されるガウス・マルコフ確率場
条件付き分布で構成されるGMRFは疾病マッピングの文脈で重要な役割を果たしてきました．こういったモデルは，空間相関や共分散関数のモデリングのために定式化されてきました．CARモデルのように条件付き分布からGMRFを構成することはコンセプト的にも，数値計算的にも（ギブスサンプラーが使えるので）望ましいことでした．ここではいくつかのCARモデルを紹介します．
proper CAR (pCAR)： pCARの条件付き期待値と分散は次のとおりです
E[\psi_i|\psi_{-i}]=\frac{\lambda\sum_{j\sim i} \psi_j}{w_{i+}},\quad Var(\psi_i|\psi_{-i})=\frac{\sigma^2}{w_{i+}}
ここで$\lambda\in(0,1)$は空間スムージングをコントロールするパラメータで，$w_{i+}=\sum_{j\sim i} w_{ij}$は隣接地域の数を表します．pCARの条件付き期待値は隣接地域のリスクの平均に比例します．また隣接地域数が多いと，条件付き分散が小さくなり，より多くの情報を借り合うことになります．
互いに隣接している２つの地域$i$, $j$について（$i\sim j$，それぞれの隣接地域数が異なるとする：$w_{i+}\neq w_{j+}$），$\psi_i$に対して$\psi_j$が影響を与える度合いは$\lambda/w_{i+}$なのに対し，$\psi_j$に対して$\psi_i$が影響を与える度合いは$\lambda/w_{j+}$です．これは空間上の依存性が非対称であることを表していて，別の言い方をすると，隣接地域数が多い地域は，その隣接地域のうち隣接地域数が少ないものの影響を受ける度合いが低くなるということを意味しています．疾病マッピングにおいては，リスクの予測精度が高い地域は，リスクの予測精度が低い地域の影響を受けにくいはずですので，このモデルは自然なものであると考えられます．
LerouxのCAR (LCAR)： LCARの条件付き期待値と分散は次のとおりです
E[\psi_i|\psi_{-i}]=\frac{\lambda\sum_{j\sim i} \psi_j}{1-\lambda+\lambda w_{i+}},\quad Var(\psi_i|\psi_{-i})=\frac{\sigma^2}{1-\lambda + \lambda w_{i+}}
LCARでは隣接地域からの影響は，$\lambda$の非線形な増加関数になっています．$\lambda$は空間スムージングのパラメータで，iCAR（次に紹介）とiid正規分布の精度行列を混合し，局所的なスムージングとグローバルスムージングの折衷を行い，また条件付き分散のコントロールも行います．
intrinsic CAR (iCAR)： iCARの条件付き期待値と分散は次のとおりです
E[\psi_i|\psi_{-i}]=\frac{\sum_{j\sim i} \psi_j}{w_{i+}},\quad Var(\psi_i|\psi_{-i})=\frac{\sigma^2}{w_{i+}}
iCARの条件付き期待値は隣接地域のリスクの平均となっています．iCARでは精度行列のランクが$n-1$となっており，pCARとLCARの$\lambda\to1$の極限となっています．
BYMとMBYM: iCARは疾病マッピングではよく用いられており，対数相対リスク$\psi$は空間相関のある効果$\psi^s\sim \text{iCAR}(\sigma^2_s)$と空間相関のない効果$\psi^h\sim N(0,\sigma^2_h I_n)$の和$\psi=\psi^s+\psi^h$でモデル化することが多いです（本論文中においてBYMというモデルに該当）．またBYMを改善したものとして$\psi=\sqrt{\lambda}\psi^s+\sqrt{1-\lambda}\psi^h$とするものも提案されています（MBYM）．
これら5つの異なる空間構造を仮定するモデルの中でどれが一番よいかというのは問題に依存するようで，モデル比較・選択を通して手元の問題に最も適したものを使用するのが通常のようです．
おわりに
病気の発生や死亡がどのように空間上に分布しているかを明らかにする疾病マッピングは，がんだけではなく最近のCOVID-19を含めいろいろな病気に対して応用がされています．本記事では，疾病マッピングモデルの基本的な形や空間上の相関構造をどのようにモデリングするかについて紹介しました．次回はモデルのさらなる拡張について紹介したいと思います．
一緒にお仕事しませんか！
株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリーや統計学に関する研修につきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！
