1. はじめに
近年の計量経済学と機械学習の融合の成果として CausalImpact パッケージが統計的因果推論の枠組みで用いられている。
Brodersen et al. (2015) が考案した causal impact フレームワークは, Rubin 流の因果推論である Difference in Differences (DID, DD) 推定法DIなどの考え方を踏まえて, ベイズ構造時系列モデルの特性を取り入れているのが特徴である。詳しくは、参考のブログ記事を参照。
2.野生イノシシへの豚熱感染による分布の減少
さて、ここで本題の統計的因果推論による検証 したいデータを紹介する。まずは下記のニュースを見てください。
2018年9月以降、各県で野生イノシシにおける豚熱の感染死亡個体が確認され、全国に拡大している。
豚熱というのは、以前は豚コレラと呼ばれていて、豚とイノシシにとって高熱性の伝染病であり、強い伝染力と致死率が特徴であるが、人に感染することはなく、仮にその肉を人が食べても問題はない。しかし、養豚農家にとっては一大事になる。
その豚熱の感染拡大によって、野生イノシシの分布に影響を与えていることが予想されており、上記のニュースでもふれられているが、このことを統計的に検証したデータ解析はネットでは見当たらない。そこで、全国の野生イノシシの捕獲数データから検証してみたい。
2．検証するデータ
データは、各県のHPから公開されている平成30年（2018年）から令和３年（2021年）までの野生イノシシの捕獲数のデータを用いた。
（ただし、全ての県のデータを用いていない）
Aグループ
2018年に野生イノシシの死亡個体に感染が確認された岐阜、愛知、三重県の３県
Bグループ
2019年に野生イノシシの死亡個体に感染が確認された福井、滋賀、栃木県、富山県の４県
Cグループ
2020年に野生イノシシの死亡個体に感染が確認された、奈良県、大阪府、京都府、新潟県、神奈川県、千葉県、茨城県、福島県、山形県
の９県
Dグループ
2021年までは野生イノシシの死亡個体に感染が確認されていない、徳島県、香川県、愛媛県、高知県、鹿児島県、熊本県 、長崎県、佐賀県、山口県、岡山県、島根県、鳥取県、宮城県、岩手県の14県
3．捕獲頭数の推移
捕獲頭数の推移を各グループごとに見てみる。点線のラインが野生イノシシの死亡個体に感染が確認された年になる
ここで使用するRのライブラリ
library(tidyverse)
library(CausalImpact)
library(patchwork)
データをインポートする
d = read.csv("https://raw.githubusercontent.com/kenkenvw/data/main/wild_boar.csv",fileEncoding = "shift-jis")
各グループごとの捕獲頭数推移を見る
d1 = d %>% select(1,3:5) %>%
pivot_longer(2:4, names_to = "pref",values_to = "numbers")
colnames(d1)[1]="year"
g1 = ggplot(d1,aes(x=year ,y=numbers )) + geom_line(aes(color=pref)) +geom_vline(xintercept = 2018,linetype=2) + theme_bw()+xlab("A-group")
d2 = d %>% select(1,6:9) %>%
pivot_longer(2:5, names_to = "pref",values_to = "numbers")
colnames(d2)[1]="year"
g2 = ggplot(d2,aes(x=year ,y=numbers )) + geom_line(aes(color=pref)) +geom_vline(xintercept = 2019,linetype=2) + theme_bw() +xlab("B-group")
d3 = d %>% select(1,10:18) %>%
pivot_longer(2:10, names_to = "pref",values_to = "numbers")
colnames(d3)[1]="year"
g3 = ggplot(d3,aes(x=year ,y=numbers )) + geom_line(aes(color=pref)) +geom_vline(xintercept = 2020,linetype=2) + theme_bw()+xlab("C-group")
d4 = d %>% select(1,19:32) %>%
pivot_longer(2:15, names_to = "pref",values_to = "numbers")
colnames(d4)[1]="year"
g4 = ggplot(d4,aes(x=year ,y=numbers )) + geom_line(aes(color=pref)) + theme_bw()+xlab("D-group")
(g1|g2)/(g3|g4)
確かに、各県に野生イノシシの死亡個体に感染が確認された年の後には捕獲頭数が減っている傾向が見られる。
4．CausalImpact による豚熱感染の影響の検証
このままのデータでは、各県の捕獲頭数そのものの違い（スケール）が大きいので（例：岩手県と長崎県）、平均で引いて標準偏差で割る処理、いわゆるデータの正規化を行うことでスケールを統一し、それぞれグループごとにデータフレームに分割した。
df = map_df(d[,-c(1,2)],scale) #西暦と和暦を除去する
df_A = df %>% select(1:3)
df_B = df %>% select(4:7)
df_C = df %>% select(8:16)
df_D = df %>% select(17:30)
1）Aグループの検証（2018年の豚熱感染の影響）
Aグループを処理群、Dグループをコントロール群とする。
test_capture = df_A
ctr_capture = df_D
豚熱感染確認の事前と事後にわける。Aグループは2018年のため、時系列の順番で言えば、３になる
pre_period = c(1,3)
post_period = c(4,6)
y = rowMeans(test_capture)
感染確認後はNAを入れる。
post_period_response = y[post_period[1] : post_period[2]]
y[post_period[1] : post_period[2]] = NA
y
>[1] 0.1969665 0.7024283 0.6373450 NA NA NA
ここでbstsモデルを走らせる
※bsts: Bayesian Structural Time Series（ベイズ構造時系列）の時系列分析を行うパッケージ。
状態空間モデルを扱うことができる。ローカル線形トレンドモデルとした。
ss = AddLocalLinearTrend(list(), y)
bsts_model = bsts(y ~ as.matrix(ctr_capture), ss, niter = 4000)
予測モデルを見る
plot(bsts_model,"comp")
時点４からは予測モデルとなっている。trendは上昇傾向が見られる。
それではCausalImpactを用いて影響を見てみる。
ci = CausalImpact(bsts.model = bsts_model,
post.period.response = post_period_response)
plot(ci)
青い部分が95%信用区間となっている。一番下の累積による影響を見てみると、統計的に有意に減少傾向が確認された。
さらに自動作成されたレポートを見る（DeepLで翻訳）
summary(ci,"report")
Analysis report {CausalImpact}
During the post-intervention period, the response variable had an average value of approx. -0.51. By contrast, in the absence of an intervention, we would have expected an average response of 1.05. The 95% interval of this counterfactual prediction is [0.062, 2.00]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is -1.56 with a 95% interval of [-2.51, -0.57]. For a discussion of the significance of this effect, see below.
Summing up the individual data points during the post-intervention period (which can only sometimes be meaningfully interpreted), the response variable had an overall value of -1.54. By contrast, had the intervention not taken place, we would have expected a sum of 3.16. The 95% interval of this prediction is [0.19, 6.00].
The above results are given in terms of absolute numbers. In relative terms, the response variable showed a decrease of -157%. The 95% interval of this percentage is [-326%, -121%].
This means that the negative effect observed during the intervention period is statistically significant. If the experimenter had expected a positive effect, it is recommended to double-check whether anomalies in the control variables may have caused an overly optimistic expectation of what should have happened in the response variable in the absence of the intervention.
The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = 0.004). This means the causal effect can be considered statistically significant.
分析報告{CausalImpact}。
介入後の期間中、反応変数の平均値は約-0.51であった。対照的に、介入がない場合、平均1.05の反応が予想された。この反実仮想予測の95%区間は[0.062, 2.00]である。観察された反応からこの予測を引くと、介入が反応変数に及ぼした因果効果の推定値が得られる。この効果は-1.56で、95%区間は[-2.51, -0.57]である。この効果の有意性については、下記を参照。
介入後の期間の個々のデータポイントを合計すると（意味のある解釈ができるのは時々だけである）、反応変数の全体的な値は-1.54であった。対照的に、介入が行われなかったとしたら、合計値は3.16になると予想された。この予測の95%区間は[0.19, 6.00]である。
上記の結果は絶対数で示されている。相対値では、応答変数は-157%の減少を示した。このパーセンテージの95%区間は[-326%, -121%]である。
これは、介入期間中に観察された負の効果が統計的に有意であることを意味する。実験者がプラスの効果を期待していた場合、コントロール変数の異常が、介入がない場合に応答変数で起こるはずのことを過度に楽観的に期待する原因となっていないかどうかを再チェックすることが推奨される。
偶然にこの効果が得られる確率は、非常に小さい（ベイズ片側テールエリア確率 p = 0.004）。これは、因果効果が統計的に有意であると考えられることを意味する。
*** Translated with www.DeepL.com/Translator (free version) ***
豚熱による影響の効果は「この効果は-1.56で、95%区間は[-2.51, -0.57]である」と報告され、信用区間が０をまたがないため、統計的に有意である。
「偶然にこの効果が得られる確率は、非常に小さい（ベイズ片側テールエリア確率 p = 0.004）。これは、因果効果が統計的に有意であると考えられることを意味する」
ここで述べられている効果の単位は、正規化されたスケールであることに留意。単位１は１標準偏差のスケールになる。
2）Bグループの検証（2019年の豚熱感染の影響）
引き続きBグループを処理群、Dグループをコントロール群とする。
test_capture = df_B
ctr_capture = df_C
pre_period = c(1,4)
post_period = c(5,6)
y = rowMeans(test_capture)
感染確認後はNAを入れる。
post_period_response = y[post_period[1] : post_period[2]]
y[post_period[1] : post_period[2]] = NA
y
>[1] 0.3508710 0.3074623 0.6279529 0.9792244 NA NA
bstsモデルを走らせる。
ss = AddLocalLinearTrend(list(), y)
bsts_model = bsts(y ~ as.matrix(ctr_capture), ss, niter = 4000)
予測モデルを見る。（時点５からが予測）
plot(bsts_model,"comp")
予測は先ほどよりも上昇トレンドが強く見られた。それではCausalImpactを用いて影響を見てみる。
ci = CausalImpact(bsts.model = bsts_model,
post.period.response = post_period_response)
plot(ci)
Aグループよりも影響は、はっきりしてその効果も大きそうだ。自動レポートを見る。
summary(ci,"report")
分析報告{CausalImpact}。
介入後の期間中、反応変数の平均値は約-1.13であった。対照的に、介入がない場合、平均1.29の反応が予想された。この反実仮想予測の95%区間は[0.58, 2.07]である。観察された反応からこの予測を引くと、介入が反応変数に及ぼした因果効果の推定値が得られる。この効果は-2.43で、95%区間は[-3.20, -1.71]です。この効果の有意性については、下記を参照。
介入後の期間の個々のデータポイントを合計すると（意味のある解釈ができるのは時々だけである）、反応変数の全体的な値は-2.27であった。対照的に、介入が行われなかったとしたら、合計値は2.59になると予想された。この予測の95%区間は[1.16, 4.14]である。
上記の結果は絶対数で示されている。相対値では、応答変数は-196%の減少を示した。このパーセンテージの95%区間は[-294%, -155%]である。
これは、介入期間中に観察された負の効果が統計的に有意であることを意味する。実験者がプラスの効果を期待していた場合、コントロール変数の異常が、介入がない場合に応答変数で起こるはずのことを過度に楽観的に期待する原因となっていないかどうかを再チェックすることが推奨されます。
偶然にこの効果が得られる確率は、非常に小さい（ベイズ片側テールエリア確率 p = 0）。これは，因果効果が統計的に有意であると考えられることを意味する．
*** Translated with www.DeepL.com/Translator (free version) ***
「この効果は-2.43で、95%区間は[-3.20, -1.71]」となり、効果は統計的に有意となっている。Ａグループの-1.56よりも影響が大きい。
3）Ｃグループの検証（2020年の豚熱感染の影響）
残念ながら、このグループは事後の年のデータが１年しかなく、分散（標準偏差）がでないので、計算ができない。R4年のデータが出そろうと計算可能になる。
5．まとめ
2018年、2019年に豚熱の感染んが確認された各県と豚熱の感染が確認されていない県のデータを比較した結果、両方のグループともに統計的に有意に捕獲数が減り、2018年よりも2019年の感染確認の方が影響が大きい結果となった（理由は不明）2020年に感染確認された県のデータの解析ができれば、もしかすると、何か他に与えている影響があることがわかるかもしれない。
（上記のことは、あくまでも個人的見解です）
6．参考
[R] CausalImpact でできること, できないこと
蔓延防止等重点措置（まん防）の効果検証を「あえて」DID＋TSclustによる時系列クラスタリング＋CausalImpactでやってみた