地域間の従属性を考慮した頻度モデル

2020年2月1日
川上 良一（大同火災海上保険）
佐野 誠一郎（共栄火災海上保険）
1

本発表研究の背景• ASTINの論文「TERRITORIAL RISKCLASSIFICATION USING SPATIALLY
DEPENDENT FREQUENCY-SEVERITYMODELS」（Peng Shi and Kun Shi ）について研究を行っていたことがスタート• 損害保険分野では損害について頻度と損傷度に分けた分析が行われ、この中ではGLMが良く用いられる• 上記論文においては、米国マサチューセッツ州の自動車事故についてGLMに地域間の従属性を織り込んだ頻度と損傷度の考察を行っている• データサイエンスでよく用いられるR言語
にて日本の実データでも地域間の従属性を考慮した「空間モデリング」が行えないか？
2

R言語とは• R言語はオープンソース・フリーソフト
ウェアの統計解析向けのプログラミング言語及びその開発実行環境• データサイエンスに用いる言語であり、Pythonとともに人気の高い言語• 統計解析に強みがあり、アクチュアリーにも人気が高い• グラフィカルな図を描け、GISも行える• 本発表ではRにてモデリング等を行っている
3

モデルとは〇モデルとは• プラ「モデル」はプラスチック素材によるモデルであり、例えば飛行機のプラモデ
ルは質量は実物の飛行機とは大きく異なり機能も無視しているが、形や色は精工にできている• このようにモデルはあくまでも本物ではなく不必要な機能は落として、必要なエッセンスのみを取り上げたもの• 統計を用いて何らかの事象を数式で表すものを「統計モデリング」という〇モデリングの効果• 観測値の増減に対して、合理的に要因の変動要因が説明可能となる• ある観察事象はたまたま得られた結果であると考えられ、全く同一の条件でも観測結果は変わると考えられる。モデリングにより、真の値に近い値を求められるかもしれない• 例えばある地域の観測値がたまたま0であった場合に、実績値のみを用いると保険料の設定ができないが、モデリングで推計する場合は料率の設定が可能となる• 空間モデリングでは不確実性を考慮して背後に潜む空間パターンを解析することができる
4

時系列モデルと空間自己相関モデル〇時系列モデル• 時間的に変化する事象について、モデル化したもの• 時刻tの現象は時刻t-1に依存するとする右図のモデルが事例としてある（𝜌 = 1ならランダムウォーク）• 各事象は時間的に独立ではなく、時刻が近いほど同じような結果が得られやすい
似ている
𝑦𝑡 = 𝜌𝑦𝑡−1 + 𝜀〇空間自己相関モデル• 地理的に変化する事象について、モデル化したもの• 各事象は地域について独立ではなく、地域相関がある• 例えば、何らかの統計について、宜野湾市の観測値は浦添市や西原町の観測値に影響を受けている可能性が考えられるが、一方で国頭村の影響は受けていないかもしれない• 本発表では空間モデルのうち、CARモデルを取り上げる

時刻 𝑡の状況はその一つまえの𝑡 − 1に𝜌を乗じて、更にランダム項として𝜀を追加したものとして表される

5

状態空間モデル（時系列モデル・空間モデル）〇状態空間モデル• 状態空間モデルとは、実際の状態を表す変数と実際に観測できる変数が異なるような系を数式で表現したもの• 状態は実際には観測されず、実際に観測されるのは観測値のみ
ローカルレベルモデルの例
状態

𝛼𝑡−1

𝛼𝑡𝛼𝑡+1観測モデル

𝑦𝑡 = 𝛼𝑡 + 𝜀𝑡 , 𝜖𝑡 ~𝑁 0, 𝜎𝜀2

システムモデル
観測値

𝑦𝑡−1

𝑦𝑡

𝑦𝑡+1

𝛼𝑡 = 𝛼𝑡−1 + 𝜂𝑡 , 𝜂𝑡 ~𝑁 0, 𝜎𝜂2𝜀𝑡 は観測誤差、 𝜂𝑡 はシステム誤差〇状態空間モデルのメリット• 観測値は観測誤差を含んでおり、状態空間モデルでは真の値を推定することができる。
• 例えば、何らかの死亡率は実際に観察された死亡のみが観測値であり、観測されていない死亡がある。
6

空間隣接行列と空間重み行列〇空間隣接行列と空間重み行列• 空間隣接行列𝐶は地域の隣接状況を行列形式で表したもの• 空間重み行列𝑊は隣接状況より何らかの重みづけをしたもの〇５市町村による事例• 𝑦=𝑦1宜野湾
𝑦2中城
𝑦3浦添 のとき、空間隣接行列𝐶と空間重み行列𝑊の一例はつぎのとおり
𝑦4西原
𝑦5那覇• 𝐶=0111001• 𝑊=

2131410010100111110111133310 00 0114410 0221301200
1 ：地域が隣接していると1、隣接していないときは０
10

𝐶の1行目は宜野湾市を表しており、中城、浦添、西原に隣接し、那覇とは隣接していない00131：標準化（この例では行の合計を１とする）

𝑊は𝐶のそれぞれの地域の隣接地域数で除すことにより求める407

簡単な空間モデルの例•𝒚=𝑦1宜野湾
𝑦2中城
𝑦3浦添 , 𝑊 =𝑦4西原
𝑦5那覇

012131413130 00 0140 01412131213012001314,𝐼 =10000010000010000010000010𝒚 = 𝜌𝑊𝒚 + 𝜺
𝜺~𝑁 0, 𝜎 2 𝐼• 𝑦1宜野湾 = 𝜌1
11𝑦2中城 + 𝑦3浦添 + 𝑦4西原
333+𝜀

• 宜野湾の観測値は中城、浦添、西原のそれぞれ1/3の観測値の合計値を𝜌倍しランダム項𝜀を加えたもの• このモデルは簡単なモデルのため拡張が必要
8GLM（Generalized linear model:一般化線形モデル）とは• 線形回帰（回帰直線）を一般化したもので、線形回帰の残差（実績値と予測値の差分）は正規分布に従っているが、GLMでは残差を任意の分布としたモデル• 一般化線形モデルには線形回帰、ポアソン回帰、ロジスティック回帰などが含まれ、損害保険で基本的な予測モデリングの技法
線形回帰

GLM（ポアソン回帰）

log 𝜆 = 𝛽0 + 𝛽1 𝑥1 + 𝛽2 𝑥29階層ベイズとは• GLMを拡張した階層ベイズでは似たようなパラメータに共通の制約を与えることによってデータの個数よりもパラメータの多い場合でも統計モデルをあてはめることができる• ベイズ推計では最尤法と異なり、事後分布の形を推定することができる• 後述のCARモデルでは地域を変数とするが、階層ベイズを用いることにより、変数の数を抑えることができ、また、事前分布により条件を縛ることで滑らかさを表現する＜ベイズの定理＞
𝑃 𝑋𝜃 𝑃 𝜃P 𝜃𝑋 =∝𝑃 𝑋𝜃 𝑃 𝜃
𝑃 𝑋
事後分布＝尤度×事前分布
本論文で用いたCARBayesパッケージではMCMC（マルコフチェーン・モンテカルロ）により各パラメータを推定（上記は𝜌）10

課題〇 例えばGLMによる頻度モデルを考える
一般統計には地域を変数に持つものがよくある• 患者調査（厚生労働省）• 地域がん登録（国立がん研究センター）• 人口動態調査（厚生労働省）〇 地域間に何らかの影響が見られることもよくある• 入院発生率（地域間のベッド数の差異）• がん罹患率（生活習慣による地域差）〇 しかし、地域を変数に用いることはあまりない

患者調査より傷害による受療率（都道府県別）

11

課題
〇 なぜか
使わない場合使う場合〇 地域間の差異が明確な場合には説明力が不十分〇 各地域を独立に見るため、地域間の影響を反映しない〇 精度が低い〇 データが地域ごとに分断されることで信頼性が低下

どっちもどっち

地域間の影響を考慮したモデリングがしたい12課題〇 使用するデータ• 人口動態調査（厚生労働省）より市町村・死因・性別の死亡者数• 母数として、人口推計（厚生労働省）〇 沖縄県（うち本島の26市町村）に限定して、市町村別に見る〇 死因のうち、地域間の影響がありそうな心疾患を分析

人口動態調査より死因別死亡率トップ4
13

地域間の影響を測る〇 モデリングの前に、地域間の影響の有無を探索する〇 空間自己相関・・・近くにある地域間の観測値が持つ相関

負の相関

相関なし

正の相関

14

地域間の影響を測る〇 隣接関係を定義する・・・色々な隣接の型

隣接（ルーク型）

ドロネー三角網

近隣4ゾーン15

地域間の影響を測る〇 モランI統計量σ𝑖 σ𝑗 𝑤𝑖𝑗 𝑧𝑖 − 𝑧ҧ 𝑧𝑗 − 𝑧ҧ𝑛𝜌𝐼 =σ𝑖 σ𝑗 𝑤𝑖𝑗
σ𝑖 𝑧𝑖 − 𝑧ҧ 2
• 𝑧𝑖 ：地域𝑖の観測値• 𝑧：観測値の平均
ҧ• 𝑤𝑖𝑗 ：地域𝑖と地域𝑗が隣接する場合は1、そうでない場合は0• 𝑛：地域の数〇 -1（負の相関）～1（正の相関）の値を取る〇 R package「spdep」のmoran.test関数を使う
16

地域間の影響を測る〇 モランI統計量（結果）> moranMoran I test under randomisation

人口動態調査による
死因別死亡率の場合

data: tizu_data_genin$rateweights: listw
Moran I statistic standard deviate = 5.1373, p-value = 1.393e-07
alternative hypothesis: greater
sample estimates:
Moran I statisticExpectation
Variance
0.66021274-0.04000000
0.01857753

17地域間の影響を測る〇 （参考）心疾患のモラン散布図
Low-High
周辺は高いが自地域は低い

Low-Low
周辺も自地域も低い

High-High
周辺も自地域も高い

High-Low
周辺は低いが自地域は高い
18

地域間の影響を考慮したモデリング〇 条件付自己相関モデル（Conditional AutoRegressive model）
・・・地域𝑘（ 𝑘 = 1, ⋯ , 𝐾 ）の値𝑦𝑘 が次式の条件付分布に従う
𝐾𝑦𝑘 ȁ𝑦−𝑘 ~𝑁 𝜇𝑘 + ෍𝑗=1𝑐𝑘𝑗 𝑦𝑗 − 𝜇𝑗 , 𝜎𝑘 2• 𝑐𝑘𝑗 ：定数、𝑐𝑘𝑘 = 0・・・地域間の従属関係を表す• 𝑦−𝑘 ：地域𝑘以外の全ての地域〇「𝑦𝑘 は周辺地域の値𝑦𝑗 と従属性𝑐𝑘𝑗 で決まる」という事前分布を与える

19

地域間の影響を考慮したモデリング〇 𝒚 = 𝑦𝑘 の同時密度関数が一意に存在する（Brook’s lemma）
𝑦1𝒚 = ⋮𝑦𝐾~𝑁𝐾𝜇1⋮ , 𝐼𝐾 − 𝐶 −1 𝑀𝜎𝜇𝐾

• 𝐶 = 𝑐𝑘𝑗 ：𝐾 × 𝐾行列• 𝑀𝜎 ： 𝑀𝜎 𝑘𝑘 = 𝜎𝑘 2 の𝐾 × 𝐾対角行列• 𝑐𝑘𝑗 𝜎𝑘 2 = 𝑐𝑗𝑘 𝜎𝑗 2 を満たす

〇 𝒚の分散構造が（特に精度行列が𝑀𝜎 −1 𝐼𝐾 − 𝐶 として）決まっているため、計算効率がよい〇 この条件の中で色々なモデルが存在する・・・例えば
20

地域間の影響を考慮したモデリング〇 かんたんなCARモデルの例
・・・Intrinsic CARモデル（ICAR）+ 混合モデル
𝑦𝑘 = 𝑥𝑘 𝛽 + 𝜙𝑘 , 𝑘 = 1, ⋯ , 𝐾
固定効果

𝜙𝑘 ȁ𝜙−𝑘 ~𝑁

変量効果

σ𝑗 𝑤𝑘𝑗 𝜙𝑗σ𝑗 𝑤𝑘𝑗

𝜎2,σ𝑗 𝑤𝑘𝑗

• 𝜙−𝑘 ：地域𝑘以外の全ての地域に関する変量効果• 𝑤𝑘𝑗 ：地域𝑘と地域𝑗が隣接する場合は1、そうでない場合は0• 𝜎 2 ：分散パラメータ⇒「各地域は隣接する地域の加重平均」と考える
21

地域間の影響を考慮したモデリング〇 色々なCARモデル・・・R package「CARBayes」より
条件付分布種類

1

特徴
2

𝜙𝑘 = 𝜙𝑘 + 𝜙𝑘
BYM

2

𝜏
σ𝑗 𝑤𝑘𝑗

ICARに無構造の分散𝜙𝑘 を加える
ことで、周辺地域に影響されない誤差を反映したモデル

𝜆 σ𝑗 𝑤𝑘𝑗 𝜙𝑗𝜏2𝜙𝑘 ȁ𝜙−𝑘 ~𝑁,σ𝑗 𝑤𝑘𝑗
σ𝑗 𝑤𝑘𝑗

ICARに対して𝜆を置くことにより分散共分散行列を正則にしたモデル

𝜌 σ𝑗 𝑤𝑘𝑗 𝜙𝑗
𝜏2
𝜙𝑘 ȁ𝜙−𝑘 ~𝑁,𝜌 σ𝑗 𝑤𝑘𝑗 + 1 − 𝜌 𝜌 σ𝑗 𝑤𝑘𝑗 + 1 − 𝜌

𝜌により変量効果の中で従属性の強さをコントロールできるようにしたモデル
𝜌 = 1のときはICARに等しく、
𝜌 = 0のときは全て独立となる

11

𝜙𝑘 ቚ 𝜙−𝑘 ~𝑁

σ𝑗 𝑤𝑘𝑗 𝜙𝑗 1σ𝑗 𝑤𝑘𝑗

2,

2

𝜙𝑘 ~𝑁 0, 𝜎 2ProperLeroux22

地域間の影響を考慮したモデリング〇 分析例（Lerouxモデル+GLM）・・・死亡者数𝑦𝑖 を予測する

𝑦𝑖 ~𝑃𝑜𝑖𝑠𝑠𝑜𝑛 𝜆𝑖log 𝜆𝑖 = log𝐸𝑖 + 𝑥𝑖 𝛽𝑔𝑒𝑛𝑑𝑒𝑟 + 𝜙𝑘 𝑖
• 地域の変量効果𝜙𝑘 𝑖 にはLerouxモデルを適用する• 𝐸𝑖 ：人口（オフセット項）• 𝛽𝑔𝑒𝑛𝑑𝑒𝑟 ：性別〇 比較としてGLM（地域変数（固定効果）の有無別）も見る• log 𝜆𝑖 = log𝐸𝑖 + 𝑥𝑖 𝛽𝑔𝑒𝑛𝑑𝑒𝑟 （地域変数なし）• log 𝜆𝑖 = log𝐸𝑖 + 𝑥1𝑖 𝛽𝑔𝑒𝑛𝑑𝑒𝑟 + 𝑥2𝑖 𝛽𝑟𝑒𝑔𝑖𝑜𝑛 （地域変数あり）
23

Lerouxモデル+GLMの解説〇 𝑦𝑖 ~𝑃𝑜𝑖𝑠𝑠𝑜𝑛 𝜆𝑖• 心疾患による死亡者数はポアソン分布に従う。ただし観察値の属性（性別や市町村）により、ポアソン分布のパラメータは異なる〇 log 𝜆𝑖 = log𝐸𝑖 + 𝑥𝑖 𝛽𝑔𝑒𝑛𝑑𝑒𝑟 + 𝜙𝑘 𝑖
• 式変形すると𝜆𝑖 = 𝐸𝑖 ∙ exp 𝑥𝑖 𝛽𝑔𝑒𝑛𝑑𝑒𝑟 ∙ exp 𝜙𝑘 𝑖
• リンク関数として対数をとることで、線形予測子としてポアソン分布のパラメータ𝜆𝑖 を表すことができる〇 𝜙𝑘 𝑖 ห𝜙−𝑘 𝑖 ~𝑁

𝜌 σ𝑗 𝑤𝑘 𝑖 𝑗 𝜙𝑗,

𝜏2𝜌 σ𝑗 𝑤𝑘 𝑖 𝑗 +1−𝜌 𝜌 σ𝑗 𝑤𝑘 𝑖 𝑗 +1−𝜌

• 地域の変量効果を、次の階層構造を持つモデルでベイズ推計する• 𝜏 2 ~Inverse−Gamma 1,0.01 、𝜌~Uniform 0,124

地域間の影響を考慮したモデリング〇 分析例（Lerouxモデル+GLM）• 2012年度～2016年度の5年分のデータを使用• 2017年度を評価に用いる

〇 結果

実績

CARGLM25

地域間の影響を考慮したモデリング〇 結果（地域変量効果の予測値（MAP推定量））
市町村

exp 𝜙

市町村

exp 𝜙

市町村

exp 𝜙

市町村

exp 𝜙

那覇市
宜野湾市
浦添市
名護市
糸満市
沖縄市

0.90
0.72
0.63
1.16
0.92
0.85

豊見城市
うるま市
南城市
国頭村
大宜味村
東村

0.71
0.99
1.15
1.88
2.22
1.38

今帰仁村
本部町
恩納村
宜野座村
金武町
読谷村
嘉手納町

1.44
1.61
1.30
1.27
1.19
0.88
1.01北谷町
北中城村
中城村
西原町
与那原町
南風原町
八重瀬町

0.72
0.84
0.68
0.66
0.94
0.66
0.83• 数値は死亡率に対する相対的な倍率を表す。最大は大宜味村、次に国頭村、本部町と続き、北部地域が上位を占めている• 最少は浦添市、その次に西原町、南風原町と中南部寄りの地域が多い• 例えば地域の変量効果が死亡率に与える影響として、大宜味村と浦添市の格差は約３.5倍（2.22Τ0.63）ある26

地域間の影響を考慮したモデリング〇 結果（q-qプロット）
CAR

RMSE
0.0002375022

GLM（地域変数なし）

RMSE
0.0003095988

GLM（地域変数あり）

RMSE
0.0002472108

RMSEが小さいほど精度が高い
27

まとめ〇 通常のGLMでは難しい、地域間の従属性という多面的な構造を反映したモデルとして、CARモデルを提案した〇 今回のケースでは、CARモデルは、地域相関を持つデータに対して通常のGLMよりも（少し）高い予測精度を示した〇 アクチュアリー業務のうち、結果より説明力が求められるケースが多い料率算定において、 CARモデルは、一定の条件に対して納得案のある事前情報を持つ特性から、使い勝手のよい手法であると考えられる

28

地域間の影響を考慮したモデリング（その先）〇 地域間の従属性について、時間による変化を考慮する
・・・CARモデル×時系列モデル• 変量効果に時系列を加える• R Package「CARBayesST」が対応
国立感染症研究所統計より
手足口病の発生率（都道府県・週別）

29

参考文献
• Shi, P., & Shi, K. (2017). TERRITORIAL RISK CLASSIFICATION USING SPATIALLY DEPENDENTFREQUENCY-SEVERITY MODELS. ASTIN Bulletin, 47(2), 437-465. doi:10.1017/asb.2017.7• BRECHMANN, E. and CZADO, C. (2014) Spatial modeling. In Predictive Modeling Applicationsin Actuarial Science: Volume I, Predictive Modeling Techniques (eds. E.W. Frees, G. Meyers andR.A. Derrig), pp. 260–279. Cambridge: Cambridge University Press.
• C. Brunsdon, L. Comber, 湯谷啓明(訳), 工藤和奏(訳), 市川太祐 (訳) (2018) Rによる地理空間デー
タ解析入門. 共立出版.
• Bivand R, Wong DWS (2018). “Comparing implementations of global and local indicators ofspatial association.” TEST, 27(3), 716–748. https://doi.org/10.1007/s11749-018-0599-x.
• Duncan Lee (2013). CARBayes: An R Package for Bayesian Spatial Modeling with Conditional
Autoregressive Priors. Journal of Statistical Software, 55(13), 1-24. URL
http://www.jstatsoft.org/v55/i13/.
• Lee D, Rushworth A, Napier G (2018). “Spatio-Temporal Areal Unit Modeling in R withConditional Autoregressive Priors Using the CARBayesST Package.” _Journal of Statistical
Software_, *84*(9), 1-39. doi: 10.18637/jss.v084.i09 (URL:
https://doi.org/10.18637/jss.v084.i09).

• 金明哲(2010). Rで学ぶデータサイエンス７ 地理空間データ分析.共立出版.

30