見出し画像

Interrupted Time Series Analyses ①

Dr.KID

今回はInterrupted Time Series Analyses（ITSA：分割時系列デザイン）についてStataでの実例を解説しようと思います。元となる論文は以下：

こちらの論文を再現したものになります。

データセットのダウンロード

まずはデータセットのダウンロードです。

findit cigsales

コマンドを入力すると、以下の検索結果が出てきますので、赤字をクリックしましょう。

画像1

そこをクリックすると、以下の画面が出てきます。

画像2

"itsa"のコマンドをインストールしたことががない方は、上の赤枠でダウンロードしておきましょう。

cigsales.dtaというデータは、下の赤枠をクリックすることでダウンロードできます。データをインストールしたら、Stataでデータを開きます。

use cigsales, clear

画像3

表を開くと、図のようにデータが確認できます。

１グループでのITS: OLS + Neway-West standard error + lag(1)

次に１グループのみでのITSをしてみましょう。
まずはデータがpanel dataであり、時間 = year、ユニット = stateという点をStataに知らせます。

tsset state year

すると、以下のように表示されます。

画像4

次に、カリフォルニア州（調査では3番目の州）を治療グループ、1989年を介入開始時期とし、単一グループのITSAを設定、介入後のトレンド推定値を求めます。
以下のコマンドは、OLS + Neway-West standard error + lag(1)を設定したモデルです。

itsa cigsale, single treat(3) trperiod(1989) lag(1) posttrend figure

以下のようなアウトプットが得られます。

画像5

表に示されているように、一人当たりのたばこ販売量の開始レベルは134パックと推定され、1989年までは毎年1.78パック（P < 0.0001, CI = [-2.57,-0.99]）の減少が見られました。

介入初年度の1989年には、一人当たりのたばこ販売量が20.06パック（P <0.0001, CI = [-29.75,-10.36]）減少し、その後（介入前と比較して）一人当たりの年間販売量が1.49パック（P = 0.002, CI = [-2.39,-0.60]）と減少しました。

また、posttrendを指定して作成した"lincom"を用いた推定値から、"Proposition 99"の導入後、一人当たりのたばこ販売量は年間3.27パック減少したことがわかりました（95%CI = [-3.83,-2.72]）。

画像6

図 はこれらの結果を視覚的に示したものです。

次に、自己相関を考慮したモデルをフィットが正しかったかどうかを確認するために， actest (Baum and Schaffer 2013) を用いて自己相関の検定を行います。

actestは、インストールが必要ですので、cigsales.dtaと同じく、findit で探してインストールすると良いでしょう。

画像7

次に、自己相関の検定をします：

actest, lags(6)

画像8

表の右側のパネルに示されているように、自己相関はlag(1 )で存在していますが、それ以上のラグオーダーでは存在していません（検定はlag(6)まで）。したがって、lag(1)を指定した最初のモデルは、この自己相関を正しく説明したと判断します。

１グループでのITS: generalized least-square + AR(1)

2つ目のモデルであるpraisは、誤差がAR(1)過程に従うと仮定し、the generalized least-squares method を用いて、線形回帰モデルのパラメータを推定します。

具体的には、praisは、プールされた自己相関推定値pに基づいて元の観測値を変換し、1次誤差間の相関（すなわち、各観測期間の誤差と先行する観測期間の誤差との相関）を除去する方法です。

もう1つの方法は、praisオプションを指定してitsaを再実行することですが、これは本来AR(1)モデルをフィットするように設計されています。ここでは、残差の自己相関に基づいてpを求める"rhotype(tscorre)"を指定し、Robust standard errorを追加しています。

itsa cigsale, single treat(3) trperiod(1989) replace prais rhotype(tscorr) vce(robust)

画像9

"prais"を用いて得られた推定値は変換されているため、OLSモデルを用いて得られた"newey"の推定値と直接比較することはできません。

"prais"は，モデルが一次自己相関をどの程度補正しているかを示す指標として，Durbin-Watson d統計を提供しています。dは0から4までの値をとることができ，帰無仮説の下ではdは2に等しい．dが2より小さい値は正の自己相関（p > 0）を示唆し，dが2より大きい値は負の自己相関（p < 0）を示唆します。

複数グループでのITSA

この例では、複数のグループに分けて、Proposition 99がカリフォルニア州の一人当たりのたばこ販売量（パック）の減少の影響を評価します。具体的には、カリフォルニア州のデータを、他の38州のデータと比較します。

itsa cigsale, treat(3) trperiod(1989) lag(1) replace figure

画像10

表に示されているように、カリフォルニア州とその他の州の最初の平均レベルの差（_z）はほぼ同じであったが（-3.27, CI = [-13.75, 7.20]）、ベースラインの平均傾斜の差（_z_t）は（-1.23, CI = [-2.14,-0.32]）異なりました。

画像11

図を見ると，38州のたばこ販売量の平均値はカリフォルニア州よりも高くなっており，そのレベルは観察期間中ずっと高くなっています。このようなベースラインの変化パターンの違いを考慮すると、他の38州はカリフォルニア州とは比較できず、したがって_z_x_1989と_z_ x_ t1989の治療効果推定値には偏りがあると言えるかもしれません（今回のケースでは、どちらの推定値も統計的には有意ではありません）。

したがって、対照群の選択をこれら2つの変数の値が類似しているものだけに限定することで、このモデルを改善することができる。

結果変数のベースラインレベルとトレンドにおいて、カリフォルニア州と比較可能な州のみに分析を限定します。ここでの比較可能性とは、zとz tの両方でp値が0.10より大きいことと定義しています。この結果、コロラド州、アイダホ州、モンタナ州が選ばれました。

itsa cigsale, treat(3) trperiod(1989) contid(4 8 19) lag(1) replace posttrend figure

画像12

画像13

表や図を見てもわかるように、治療群と対照群はベースラインレベルとトレンドの両方で比較可能です。介入1年目（_z_x1989）には治療効果は見られませんでした。しかし、介入前のトレンドでは対照群と比較して統計的に有意な年間一人当たりのたばこ販売量の減少が見られました（-1.97, CI = [-3.26,-0.68]）。

さらに、ポストトレンドの出力から、治療群は介入後の期間に年間のたばこ販売量を3.27パック減少させ、対照群は同期間の販売量を1パックしか減少させておらず、両者の差は一人当たり年間2.28パックであることがわかりました。

これらの結果は、複数グループのITSAを実施する際に、結果変数の介入前のレベルと傾向について、治療群と対照群が比較可能であることを保証することの重要性を強調している。先に述べたように、各非治療群を治療群と個別に比較する反復プロセスを用いることができる。

画像14

上記の統計モデルにおいて、β4とβ5の両方でp値が特定の閾値よりも大きいグループは、最終モデルに含めるためのコントロールとして保持することができます。

このアプローチは、他の共変量にも容易に拡張することができますが、多くの共変量でバランスを取ることが重要な要素である場合には、ITSAの2つの代替アプローチを検討する必要があります。例えば、Abadie, Diamond, and Hainmueller (2010)が説明し、synthパッケージを用いてStataで実装された合成コントロール・アプローチ（Abadie, Diamond, and Hainmueller 2014）、またはLinden and Adams (2011)が説明した傾向スコア加重手法です。

治療期間が複数ある場合

"itsa"は、複数の治療期間の効果が注目されるデザインのバリエーションに対応できます。

例えば、研究者は、介入が導入され、撤回され、再び導入された場合や、介入の後に後の時点で別の介入が行われた場合の効果を研究することに興味があるかもしれません（他の多くのデザインの選択肢については、Barlow, Hayes, and Nelson [1984]を参照）。

以下の例では、1982年からタバコの販売データに架空の介入を加えて説明します。単一グループのITSAを、介入期間を1回追加して再推定してみましょう。

itsa cigsale, single treat(3) trperiod(1982 1989) lag(1) replace posttr figure

画像15

2 回目の介入までのすべての係数の解釈は以前と同様です。
つまり、第1回目の介入期間は介入前の期間と比較される。しかし、第2回目の介入期間に追加された係数_x1989と_x_t1989は、前（第1回目）の介入期間の係数と比較されます。

表に示されているように、また、図を目で見て確認したように、1982年から「治療効果」があったことを示す証拠があり、2回目の介入（実際には真の介入期間）の実施後には、年間売上高の追加的な減少は見られませんでした。

画像16

この効果は、介入後のトレンドを第一介入期間と第二介入期間に分けて推定するposttrendオプションでさらに示すことができます。

ポストトレンドの結果、1982年以降のたばこ販売量の年間減少量は3.84パック/年でしたが、1989年以降の販売量の年間減少量は3.27パック/年とわずかに減少しました（差は0.569で、元の回帰表では_x_t1989と表示されています）。

この結果から、単一グループのデータを分析する際に、複数の治療期間を定義することで、介入の真の開始から離れた「中断」をテストすることができるという、さらなる有用性が明らかになりました。

この例のように、カリフォルニア州におけるたばこの販売動向の減少は、Proposition_99の実施よりも数年前に始まっていました。 この結果は、ITSAを使用する場合でも、反実仮想を表す比較可能な対照群を見つけることの重要性をさらに強調しています。

Stata code

*****0. Dataset
*find cigsales.dta
findit cigsales
*set dataset
use cigsales, clear

*****1. Single-group ITSA: OLS + Neway-West standard error + lag(1)
*declare the dataset as panel: unit = state, time = year
tsset state year
*specify
*a single-group ITSA with California (state number 3 in the study) as the treatment group 
*1989 as the start of the intervention, 
*request postintervention trend estimates, and plot the results
itsa cigsale, single treat(3) trperiod(1989) lag(1) posttrend figure

*find actest package
findit actest

*test for autocorrelation
actest, lags(6)

*****2. Single-group ITSA: generalized least-square + AR(1)
itsa cigsale, single treat(3) trperiod(1989) replace prais rhotype(tscorr) vce(robust)

*****3. Multiple-group ITSA
itsa cigsale, treat(3) trperiod(1989) lag(1) replace figure

*used 3 states
itsa cigsale, treat(3) trperiod(1989) contid(4 8 19) lag(1) replace posttrend figure

*****4. Multiple treatment periods
itsa cigsale, single treat(3) trperiod(1982 1989) lag(1) replace posttr figure

おまけ①：２グループ × ２つの介入がある場合

同じ著者らが２グループ × ２つの異なるタイミングでの介入がある場合の例を追加した論文があります。

数式的には以下の通りになります：

画像17

2つの介入を含むこのモデルでは、合計18の注目すべき測定値があります：
・治療群と対照群の介入前、介入1回目、介入2回目のトレンド: 2x3 = 6
・これらの各期間におけるトレンドのグループ間の差: 3
・治療群と対照群の各期間のトレンドの差: 2x3 = 6
（介入前対介入1回目、介入前対介入2回目、介入1回目対介入2回目）
・これらの各期間の比較におけるグループ間のコントラスト: 3
の18です。

回帰分析のアウトプットでは、これらの測定値のうち、
1) 対照群の介入前の傾向、β1
2) 治療群と対照群の介入前の傾向の差、β5
3) 対照群の介入前と介入1回目の傾向の差、β3
4) 治療群と対照群の介入前と介入1回目の傾向の差、β7
5) 対照群の介入1回目と介入2回目の傾向の変化、β9
6) 治療群と対照群の介入1回目と介入2回目の傾向の差、β11
の6つを提供しています。

12種類の指標が計算の対象となります：
1) 治療群の介入前のトレンド、β5+β1
2) 対照群の介入第1期のトレンド、β1+β3
3) 治療群の第1介入期間の傾向、β1＋β3＋β5＋β7
4) 第1介入期間における両群間の差、β5+β7
5) 対照群の第2介入期間の傾向、β1+β3+β9
6) 治療群の第2介入期間における傾向、β1＋β3＋β5＋β7＋β9＋β11
7) 第2介入期間における両群間の差、β7+β9+β11
8) 治療群の介入前と介入1期のトレンドの差、β3+β7
9) 治療群の第1介入期間と第2介入期間の間のトレンドの差、β9+β11
10) 対照群の介入前と介入第2期のトレンドの差、β3＋β9
11) 治療群の介入前と介入後のトレンドの差は、β3 +β5 +β7 +β9 +β11
12) 治療群と対照群の介入前と介入後の傾向の違いを比較したときの差、β5 + β7 + β11

まずはデータをロードして、tssetでpanel dataを指定します：

use cigsales, clear
tsset state year

そして、trperiodに１期および２期の開始時期を指定して、itsaコマンドを走らせます：

itsa cigsale, treat(3) trperiod(1982 1989) lag(1) posttrend figure

Coefficient（β）のデータは以下の通りです：

画像18

第１相と第２相のトレンドと、その比較は以下の通りです：

画像19

第２相にトレンドが治療群で大きく減少しているのが分かります。

画像20

表に示された推定値を解釈するには、この図を視覚的に確認することが役立ちます。例えば、表では、対照州の一人当たりの年間たばこ販売量の介入前の傾向は増加しているように見えるが、カリフォルニア州の販売量の介入前の傾向は減少しているように見える。
表によると、対照州では1970年から1982年の間に一人当たりの年間たばこ販売量が1.593(95% CI: [0.424,2.762])パック増加しているのに対し、カリフォルニア州では同期間に一人当たりの年間たばこ販売量が-3.390(95% CI: [-5.802,-0.978])パック減少していることが確認されました。
カリフォルニア州と対照州の介入前のトレンドの差は-1.797(95% CI: [-3.107,-0.486])でした。

*****5. Multiple treatment periods and groups
use cigsales, clear
tsset state year

*fit ITSA model
itsa cigsale, treat(3) trperiod(1982 1989) lag(1) posttrend figure

その他、知りたい内容はlincomを使用して推定可能です：

画像25

おまけ②：ITSA + マッチングを使用する場合

前回と同じく、まずデータセットをアップロードし、tssetを行います：

use cigsales, clear
tsset state year

介入前のアウトカム変数「cigsale」を用いて、治療を受けたユニット(3)と一致するコントロールを見つけます。P値 < 0.20 をカットオフ値をとし、lag(1)での自己相関を指定し、保存されるローカルマクロの名前を「controls」と指定します。

itsamatch cigsale, treatid(3) trperiod(1989) pr(0.20) lag(1) local(controls)

画像21

コントロールユニット4、8、9の3つのマッチが見つかりました。保存されているlocal(controls)を使って、itsaモデルのcontid()にこれらの値を入力します。

 itsa cigsale, treatid(3) trperiod(1989) contid(`controls') lag(1) replace figure(xlabel(1970(5)2000)) posttrend

画像22

画像23

＊上の２つのコードは同時に走らせる必要があります。

このように、治療群と対照群の介入前のトレンドはほぼ同じで、介入後のトレンドがそのまま介入効果として推定できそうです。

さらに共変数をマッチングの因子として使用する場合は、itsamatchのモデルに組み入れれば行えます：

*We now add the covariate "retprice" to varlist and rerun itsamatch.
itsamatch cigsale retprice, treatid(3) trperiod(1989) pr(0.20) lag(1) local(controls)
*  The same three matches were found, control units 4, 8, and 9. We now verify (numerically and visually) the balance on "retprice" by specifying "retprice" as the dependent variable in the itsa model.
itsa retprice, treatid(3) trperiod(1989) contid(`controls') lag(1) replace figure(xlabel(1970(5)2000)) posttrend

画像24

retail priceに関しては、２つのグループのトレンドでほとんど差がないのが分かります。

画像26

***itsa match
use cigsales, clear
tsset state year

*find controls that match the treated unit
itsamatch cigsale, treatid(3) trperiod(1989) pr(0.20) lag(1) local(controls)

*Three matches were found, control units 4, 8, and 9. We now plug those values in contid() in an itsa model using the stored local(controls)
itsa cigsale, treatid(3) trperiod(1989) contid(`controls') lag(1) replace figure(xlabel(1970(5)2000)) posttrend

*We now add the covariate "retprice" to varlist and rerun itsamatch.
itsamatch cigsale retprice, treatid(3) trperiod(1989) pr(0.20) lag(1) local(controls)
*  The same three matches were found, control units 4, 8, and 9. We now verify (numerically and visually) the balance on "retprice" by specifying "retprice" as the dependent variable in the itsa model.
itsa retprice, treatid(3) trperiod(1989) contid(`controls') lag(1) replace figure(xlabel(1970(5)2000)) posttrend


良質な記事をお届けするため、サポートしていただけると嬉しいです！英語論文などの資料取り寄せの費用として使用します。