Interrupted time series analysis でARIMAを使って真面目にやる

読んだ。
Interrupted time series analysis using autoregressive integrated moving average (ARIMA) models: a guide for evaluating large-scale health interventions

Interrupted time series analysis は介入の前後でlevel（増減、つまり切片）とtrend もしくはslope（単調な時系列増減、つまり傾き）を断片的に回帰するモデル（式は雑）
$Y_t=\beta_0+\beta_t t +\beta_{intervention}X_{intervention}+\beta_{level}t+\beta_{trend}t$
を推定することが多いが、時系列データなので普通は自己相関とか季節性とかそういうことを考慮しないといけないはずで、ARIMAを使ったという話で、実装上もRのforecast パッケージだけでできる、という話。

サンプルデータは2014年に処方規制がされたクエチアピン25mg錠の処方数をもとに解析していて、本文中にはstep（level のこと）とramp（trend のこと）のほかにpulse（介入が起きた直後にスパイク状に突然増える/減る効果がある場合があると思うが、これはすぐ次の時刻でもとに戻る）というものも想定している。また、介入後の増減の仕方はtransfer function という $\omega_0$ を用いてもっと複雑にできるが、これは $h$ と $r$ のパラメータを設定する必要があり、これはこれでまた複雑な話なので本文中ではさらっと流され、しかもどちらも $h=0$ 、 $r=0$ とされている。
ARIMAのモデル $(p, d, q)$ の決め方は、forecast パッケージのauto.arima関数を使ってAICもしくはBICで最適化されたものを勝手に採用してくれる。

解析結果としては、2014年の介入で処方数は3285、trendも1397/月減少したことがわかった。confint関数を使えば95%信頼区間も計算できるし、lmtest パッケージのcoeftest関数もしくはベタに計算すれば、点推定のp値も計算できる。

library(lmtest)
coeftest(model1)

z test of coefficients:

        Estimate  Std. Error  z value  Pr(>|z|)    
ar1     -0.87301     0.12396  -7.0427 1.885e-12 ***
ar2     -0.67314     0.12587  -5.3480 8.894e-08 ***
sma1    -0.60694     0.38723  -1.5674     0.117    
step -3284.77920   602.33653  -5.4534 4.942e-08 ***
ramp -1396.65226   106.63300 -13.0977 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(1-pnorm(abs(model1$coef)/sqrt(diag(model1$var.coef))))*2

         ar1          ar2         sma1         step         ramp 
1.884937e-12 8.894107e-08 1.170208e-01 4.941704e-08 0.000000e+00

ARIMAモデルを使っているので、もし介入がなかった場合の仮想現実counterfactrual も推定できる。この場合は赤線になる。