library(haven)
library(tidyverse)
library(estimatr)
library(stargazer)
  • ここでは差分の差分 (DID)法で使うパネルデータについて解説する　　
  • その後、差分の差分 (DID)法の仕組み、分析方法、そして実例を紹介する

1. パネルデータ

  • データは次の 3 種類に分類できる
クロスセクションデータ： 異なる個体について、同時点の観測値を集めたデータ
時系列データ： 同じ個体について、異なる時点で観測した値を集めたデータ
パネルデータ： 複数の個体について、時系列の観測値を集めたデータ

1.1 クロスセクションデータ

  • 異なる個体について、同時点の観測値を集めたデータ
  • クロスセクションデータの一例：
    ・2021年総選挙データ

1.2 時系列データ

  • 時系列データ (time-series data)
  • 同じ個体について、異なる時点で観測した値を集めたデータ
  • 時系列データの一例：
    ・高市早苗氏の総選挙結果時系列データ
  • 可視化してみる

1.3 パネルデータ

  • パネルデータ ＝ クロスセクション ＋ 時系列
  • 複数の個体について、時系列の観測値を集めたデータ
  • 観測されるデータは時点によらず同じ
  • パネルデータの一例：
    ・2021年自民党総裁選立候補者の1996年から2021年までの総選挙における得票率の推移データ

  • 可視化してみる

  • 差分の差分 (DID)法ではパネルデータを使う

2. 差分の差分法

  • 差分の差分法 (Difference-in-Differences)
  • DID, DD, dif-in-dif とも呼ばれる

DIDを使う場合

  • ある政策の効果を確かめたい場合
  • ランダム化比較試験 (RCT) が使えるならそれが望ましい
  • しかし、いつでも RCT を使えるとは限らない
  • RCT が使えない時（＝実験ができない時）
    →　調査・観察データを使って政策の効果を確かめる
  • 注意：政府や地方自治体が実行する「社会実験」と称する実験は RCT を使った因果推論における実験ではない（ラインダマイズが行われていないから）
  • このような「社会実験」では、処置を受ける企業や自治体などの個体と処置を受けない個体がランダムに決められていない = セレクションバイアスがある
  • このような「社会実験」では、国や自治体が観測単位
  • 処置群と統制群が異質（= RCT でない）
  • 全ての交絡をコントロールできない

例）：
・高速道路無料化する・・・処置群
・高速道路無料化しない・・・統制群
→　高速道路無料化しない統制群の人々から「不公平だ」という不満が寄せられるため、RCT による実験がなかなかできない

2.1 DID の仕組み

政策(処置):

  • 宮城県が、県外からの観光客のホテル・旅館の宿泊費を半額にした

結果：

  • 県外からの宿泊者数

調べたい因果効果：

  • ホテル・旅館の宿泊費半額セールが館外からの宿泊者数に与える影響
宮城県 ：ホテル・旅館の宿泊費半額サービスあり
山形県 ：ホテル・旅館の宿泊費半額サービスなし

①宿泊者数の単純比較（個体間比較）

  • 宮城県と山形県の年間観光客数の比較
都道府県名 宿泊者数半額サービス実施「後」（万人
山形県 4000
宮城県 6000
差 2000

・単純比較による因果推論結果：

・宿泊者数半額サービスが県外からの観光客数を 2000万人増やした

問題：

・宮城県と山形県は、処置「宿泊者数半額サービス」以外が同質ではない
→　そもそも宮城県の方が県外からの観光客数が多かっただけなのでは？　　

②宿泊者数の単純比較（前後比較）

  • サービス実施前の宮城県と実施後の宮城県の年間観光客数の比較
宿泊者数（万人） 宿泊者数（万人）
都道府県名 サービス実施「前」 サービス実施「後」 差（後−前）
宮城県 6500 6000 −500

・単純比較による因果推論結果：

・宿泊者数半額サービスが県外からの観光客数を 500万人減らした

問題：

  • サービス実施前の宮城県と実施後の宮城県は同質ではない
  • このままでは因果推論はできない
  • 観光客数を計測している時間が違うから
  • ここでは時間による変化がコントロールされていない：
    →　サービスを実施しなくても、観光客は減ったのでは？　　

③個体間比較と前後比較の併用　　

宿泊者数（万人） 宿泊者数（万人）
都道府県名 サービス実施「前」 サービス実施「後」 差（後−前）
宮城県 6500 6000 −500
山形県 6000 4000 -2000
差 500 2000 1500

差の差による因果推論：

  • 宿泊者数半額サービスが宿泊者を1500万人増やした
  • 個体間比較：= 2000 過大推定（宮城県の観光客はもともと多い）
  • 前後比較：= -500 過小推定（サービスを実施しないと観光客は減る）

2.2 差分の差分 (DID)

  • 差分の差分 (Difference-in-differences: DID) を利用する
  • 「個体間の差」と「時点間の差」の両方を使って因果効果を推定する方法
  • Y: 結果（宿泊者数）を計算する 2 つ方法

\(δ_{DD}\): 差分の差分による因果効果 (I)

  • 「時点間の差」と「時点間の差」の個体間（＝群間）の差を計算する

  • 宮城と山形で after の個体間を比較 = 「時点間の差」

  • 宮城と山形で before の個体間を比較 = 「時点間の差」

  • 両者（＝群間）の差を求める
    \[δ_{DD} = (Y_{宮城, after} − Y_{山形, after}) - (Y_{宮城, before} − Y_{山形, before})\\= (6000 - 6500) - (4000 - 6000)\\= -500 + 2000 = 1500\]

\(δ_{DD}\): 差分の差分による因果効果 (II)

  • 「個体間（群間）の差」と「個体間（群間）の差」の時点間の差を計算する

  • 宮城で処置の前後を比較 = 「個体間（群間）の差」

  • 山形で処置の前後を比較 = 「個体間（群間）の差」

  • 両者の時点間の差を求める

\[δ_{DD} = (Y_{宮城, after} − Y_{宮城, before}) - (Y_{山形, after} − Y_{山形, before})\\= (6000-6500) - (4000-6000)\\= -500 + 2000 = 1500\]
差分の差分による因果効果 ・宿泊者数半額サービスが宿泊者を1500万人増やした
  • 上の関係を可視化してみる

  • 図の灰色の点線は「実際にサービスを実施した宮城県が、もしサービスを実施していなかったら」を示す仮想現実 (Counter factural)を表している
    →　平行トレンド (parallel trend) の仮定

平行トレンドの仮定 ・差分の差分法で因果推論を行うために必要な仮定
・「もし処置がなかったら、処置された個体の結果（あるいは処置群の結果の平均）と統制された個体の結果（あるいは統制群の結果の平均）は平行な時間的変化を示す」

・この仮定が満たされなければ、差分の差分法で因果推論にはバイアスが生じる
  • 処置を長期にわたって観察し、平行トレンドがありそうかどうか確認する
  • 実際は、完全に平行になることは期待できない
  • 平行で「ない」ことは確認しやすいが、平行で「ある」ことを証明することはできない
  • 平行トレンドの仮定はあくまでも「仮定」であることに注意

2.3 異なる「ショック」の確認

  • 処置後に、宮城県と山形県で異なる状況が発生したなら、比較できない
  • 例えば、宿泊者数半額サービスを実施した直後に次のような「ショック」が起こった場合：

・宮城県で東京オリンピック競技が開催された
・山形県では東京オリンピック競技が開催されなかった

  • もしこのようがことがあれば、宮城県の宿泊者を増やしたのが、宿泊者数半額サービスなのか、東京オリンピックなのか区別できない
  • しかし、処置語に、宮城県と福島県が同じショックを受け、そのショックから受ける影響が同じと見なせるなら比較可能
  • 例えば、宿泊者数半額サービスを実施した直後に、新型コロナウイルスによる緊急事態宣言が発令され、旅行が制限された場合

・宮城県：緊急事態宣言の対象
・山形県：緊急事態宣言の対象

  • この場合、どちらの県も同程度のショックを受けていると考えられるので、比較可能
  • 「宿泊者数半額サービス」の実施後に、「宿泊者数」に影響を与える外生的ショックが、宮城県と山形県のどちらか一方にだけ起こっていないかどうかチェックすることが必要

2.4 DIDの利点

  • 適用できる範囲が広い
    ・処置群と統制群の両者を含むパネルデータがあれば分析可能
  • 処置群と統制群は異質でもかまわない
    ・「他の条件が等しく (ceteris paribus)」なくてもよい
    ・必要なのは「平行トレンド」の仮定だけ
  • 処置群全体の処置効果が推定できる
    ・RDDだと処置群の一部だけしか推定できないが、DID では、処置群全体の処置効果が推定できる

2.4 DIDの欠点

  • パネルデータが不可欠
    ・処置群と統制群の両者を含む長期間のパネルデータが必要
  • 平行トレンドの仮定を満たす事例を見つけるのが大変
    ・異質な個体は「異なるトレンド」を示すことが多い
    ・たとえ平行トレンドのように見えても、それが本当に平行かどうかわからない（←平行トレンドは単なる仮定に過ぎないから）

3. DID を用いた回帰分析

3.1 必要な変数

変数名 内容 データの種類
\(Y_{kt}\) ：結果変数 パネル
\(k\) ：処置が適応される単位（例：国、自治体、企業、人など）
\(t\) ：時間
\(D_k\) ：個体 \(k\) が処置群に属することを示すダミー変数 \(D = 0: 統制群、D = 1: 処置群\) クロスセクション
\(P_t\) ：時間 \(t\) が処置後 (post-treatment) であることを示すダミー変数　\(P_t = 0: 処置前、P_t = 1: 処置後\) クロスセクション

3.2 データ

  • 宮城県が、県外からの観光客のホテル・旅館の宿泊費を半額にした事例において因果効果を推定するために必要なデータは次のとおり

3.3 DID 回帰式

\[Y_{kt} = α + βD_k + γP_t + δ(D_k × P_t) + e_{kt}\]

・\(δ（デルタ）\) ・・・DID によって推定される処置効果 = 差分の差分

  • \((D_k×P_t)\) の値が 1 になるのは \(D_k=1\) でかつ \(P_k=1\) の時
    →　\(D_k=1: 処置群\) でありかつ \(P_t=1: 処置後\) の個体の時
    →　宿泊者数半額サービスを実施した 後の(\(P_t = 1\)) 宮城県(\(D_k=1\))

  • \(Y_{kt}\) ・・・時間 (\(t\)) と個体 (\(k\)) の値によって変化するパネルデータ

  • \(D_k\) ・・・ 個体 (\(k\)) に値によってのみ変化するクロスセクションデータ

  • \(P_t\) ・・・ 時間 (\(t\)) に値によってのみ変化するクロスセクションデータ

\(D_i\) と \(P_t\) の値のパターンは次の 4 つ：

\(E[Y_{kt}|D_k = 1, P_t = 1] = α + β + γ + δ\) (1)
\(E[Y_{kt}|D_k = 1, P_t = 0] = α + β\) (2)
\(E[Y_{kt}|D_k = 0, P_t = 1] = α + γ\) (3)
\(E[Y_{kt}|D_k = 0, P_t = 0] = α\) (4)
  • 処置が 1 (\(D_k = 1\)) について、処置の前 (1) から後 (2) を引く・・・前後比較

\[E[Y_{kt}|D_k = 1, P_t = 1]- E[Y_{kt}|D_k = 1, P_t = 0] ・・・(5)\]

  • 処置が 0 (\(D_k = 0\)) について、処置の前 (3) から後 (4) を引く・・・前後比較

\[E[Y_{kt}|D_k = 0, P_t = 1] - E[Y_{kt}|D_k = 0, P_t = 0]・・・(6) \]
・\((5)\) から \((6)\) を引く・・・前後比較の群間比較

\[δ = E[Y_{kt}|D_k = 1, P_t = 1]- E[Y_{kt}|D_k = 1, P_t = 0] \\ - E[Y_{kt}|D_k = 0, P_t = 1] - E[Y_{kt}|D_k = 0, P_t = 0]\]

  • DID 回帰式

\[Y_{kt} = α + βD_k + γP_t + δ(D_k × P_t) + e_{kt}\]

処置効果 \(δ\) デルタの推定 δ（\(D\) と \(P\) の交差項の係数）が差の差を推定する

3.4 実際に DID 回帰分析する上での注意

  • 差分の差分法が利用できるのは 2 個体、2 期間のパネルデータだけではない
    ・個体数は多くてもかまわない
    ・期間は長い方が良い
  • 処置は 2 値でなくても良い

DID回帰の標準誤差

  • 2 期間より長いパネルデータを分析する場合
  • 結果変数は、それぞれの個体内で強い相関をもつ
    →　通常の標準誤差では小さすぎる
    →　誤って統計的に有意だと判断してしまう
    →　観測個体ごとにクラスタ化した標準誤差 (cluster-robust standard error) を使って推定の不確実性をチェックする

DID回帰 では lm() ではなく lm_robust() を使う

DID回帰でクラスタ化した標準誤差を使う理由 ・実際の DID 分析では２期間より長い期間を分析することが多い
・２期間より長い期間のパネルデータはある程度長くなければならない

パネルデータがある程度長くなければならない理由

① ある程度長くなければ平行トレンドの仮定を確認できないから
→ ある程度長いデータがないと前後のトレンドを比較できない
② 処置のタイミングが異なることがあり得るから
・極端に早期に処置した個体と、極端に遅く処置した個体がある場合
→ 全体の期間を含むデータが必要なので、２期間より長いデータが必要

DID回帰でクラスタ化した標準誤差を使う理由
  • 時系列データでは結果変数の値が、同一個体内（= 同じ州内）での観測値が類似するという系列相関 (serial correlation)があり相互に強い相関をもつため

☆ 例）若者の死者数が多い州では、その前年も翌年も死者数も多い

・前の年のデータが分かれば次の年のデータも推測ができる
→ 一つ一つの観測値が独立に得られたとはみなせない
・時系列データは相関が強い
→ バラバラに集めたデータと比較すると、時系列データは情報量が少ない
→ 標本サイズを割り引いて考えないと、標準誤差の計算を間違える
→ 標準誤差が小さくなりすぎて、誤って統計的に有意だと判断してしまう
・ある固体内で、結果変数の値をたくさん観測しているとしても、たくさんある観測値のひとつひとつが独立だとは見なせない
・結果変数の値は固体の中でクラスタを形成している
→ それらを割り引いて標準誤差を計算する必要がある
・時系列データでは似ている個体がデータにたくさん含まれていることを想定する必要がある
・似ている個体データ（＝クラスタ）が集まっている単位を指定し、標準誤差を計算する

→ クラスタ化した標準誤差を使う必要がある
→ lm() ではなく lm_robust() を使う

4. DID の実例（最低賃金の影響）

  • ここでは、Card and Krueger (1994) が分析した、最低賃金の上昇が失業に与える影響の研究を例に考える

  • 1992年4月にアメリカのニュージャージー州 (NJ)で最低賃金が引き上げられた

  • 標準的な経済理論は、最低賃金が引き上げられれば、フルタイム労働者はレイオフされると想定

ここでのリサーチ・クエスチョン

最低賃金の上昇は、フルタイム労働者を減らすのか？

  • 1992年4月にニュージャージーの最低賃金が$4.25から$5.05に上げられた
  • 処置前の1992年2月と処置後の11月のニュージャージーとペンシルベニアにおけるファストフード産業における雇用を比較
  • 処置 (1992年4月）の前後でのニュージャージーのみの雇用の変化を単純比較すると、天候やマクロ経済学的要因などの変数をコントロールし損なう
  • ペンシルベニアを差分の差分法におけるコントロールとして含めると、ニュージャージーとペンシルベニアで共通の変数からもたらされるあらゆるバイアスが、たとえそれらの変数が観測できないとしても、コントロールされる
  • ニュージャージーとペンシルベニアは時間を通じて並行なトレンドを持つと仮定
  • ペンシルベニアでの雇用の変化は最低賃金の上昇がなかった場合のニュージャージーで起こるはずだった雇用の変化 (Counter factural) として考えることが出来る
  • ニュージャージーにおける最低賃金の上昇は、標準的な経済理論が示唆するような、失業の増加はもたらさず、フルタイム労働者の増加をもたらした
  • 下の表はカードとクルーガーが推定した雇用の処置効果を図示

処置前後のフルタイム労働者の割合 (1992年)

処置前 処置後 差
ニュージャージー 29.7 32 -2.3
ペンシルバニア 31 27.2 3.8
差 -1.3 4.8 6.1
  • ニュージャージーにおける $0.80の最低賃金の上昇がフルタイム労働者の雇用を 6.1 % ポイント増加させたと推定

4.1 データの準備

  • Card and Krueger (1994) が分析で使っているオリジナルデータを加工した card.csv をダウンロードして RProject フォルダの中にある data フォルダに入れる
card <- read_csv("data/card.csv", na = ".")
  • 変数を確認する
names(card)
 [1] "SHEET"    "CHAIN"    "CO_OWNED" "STATE"    "SOUTHJ"   "CENTRALJ"
 [7] "NORTHJ"   "PA1"      "PA2"      "SHORE"    "NCALLS"   "EMPFT"   
[13] "EMPPT"    "NMGRS"    "WAGE_ST"  "INCTIME"  "FIRSTINC" "BONUS"   
[19] "PCTAFF"   "MEALS"    "OPEN"     "HRSOPEN"  "PSODA"    "PFRY"    
[25] "PENTREE"  "NREGS"    "NREGS11"  "TYPE2"    "STATUS2"  "DATE2"   
[31] "NCALLS2"  "EMPFT2"   "EMPPT2"   "NMGRS2"   "WAGE_ST2" "INCTIME2"
[37] "FIRSTIN2" "SPECIAL2" "MEALS2"   "OPEN2R"   "HRSOPEN2" "PSODA2"  
[43] "PFRY2"    "PENTREE2" "NREGS2"   "NREGS112"
  • 加工対象となる 6 つの変数の型を数値 (numeric) に変換する
card$EMPFT = as.numeric(card$EMPFT)
card$EMPPT = as.numeric(card$EMPPT)
card$WAGE_ST = as.numeric(card$WAGE_ST)
card$EMPFT2 = as.numeric(card$EMPFT2)
card$EMPPT2 = as.numeric(card$EMPPT2)
card$WAGE_ST2 = as.numeric(card$WAGE_ST2)
  • コードブック (codebook) を読み、分析に使う変数を分かりやすい名前に変えて抜き出す
  • 変数の変換は mutate() で、変数名の変換はrename() でできるが、変換した変数だけをデータフレームに残したいので、transmute() を使う
minwage <- card %>% 
  transmute(
    state = ifelse(STATE == 1, "NJ", "PA"),  # NJ なら state = 1、PA なら state = 0
    fulltime_before = EMPFT,    # 最低時給上昇前のフルタイム労働者の数
    parttime_before = EMPPT,    # 最低時給上昇前のパートタイム労働者の数
    wage_before     = WAGE_ST,  # 最低時給上昇前の賃金
    fulltime_after  = EMPFT2,   # 最低時給上昇後のフルタイム労働者の数
    parttime_after  = EMPPT2,   # 最低時給上昇後のパートタイム労働者の数
    wage_after      = WAGE_ST2, # 最低時給上昇後の賃金
    full_prop_before = fulltime_before / (fulltime_before + parttime_before),
    full_prop_after  = fulltime_after  / (fulltime_after + parttime_after)
  ) %>% 
  na.omit()                     # 説明のため、完全ケース分析にする
DT::datatable(minwage)%>%
    DT::formatRound(columns=c('full_prop_before', 'full_prop_after'),
                    digits=2)
  • これで分析に必要なデータフレームが完成した

Card and Krueger (1994)のオリジナルデータ ・Card and Krueger (1994)が分析で使っているオリジナルデータ (Stata 用にフォーマットされている）は下のコマンドでダウンロードできる
・下のコマンドを実行する前にパソコンの RProject フォルダ内にあらかじめ data フォルダを作っておくこと
・データ (public.dat) は data フォルダ内にダウンロードされる

dir.create("tmp")
download.file(url = "http://davidcard.berkeley.edu/data_sets/njmin.zip",
              destfile = "tmp/njmin.zip")
unzip("tmp/njmin.zip", exdir = "data")
card_org <- read_table("data/public.dat", 
                  col_names = FALSE, na = ".")

注意：public.dat は .dta (Stata形式) ではない
・このデータには変数名がついておらず、1行目から値が保存されているので、col_names = FALSE と指定する
・欠測が 「.」で表されているので、na = "." を指定

・変数を確認する

names(card_org)
 [1] "X1"  "X2"  "X3"  "X4"  "X5"  "X6"  "X7"  "X8"  "X9"  "X10" "X11" "X12"
[13] "X13" "X14" "X15" "X16" "X17" "X18" "X19" "X20" "X21" "X22" "X23" "X24"
[25] "X25" "X26" "X27" "X28" "X29" "X30" "X31" "X32" "X33" "X34" "X35" "X36"
[37] "X37" "X38" "X39" "X40" "X41" "X42" "X43" "X44" "X45" "X46" "X47"

・実際に分析する前に、このオリジナルデータを分析者にとって使いやすくクリーニングする必要がある
・このセクション 4. DID の実例（最低賃金の影響）では コードブック (codebook) を読み、分析に使う変数を分かりやすい名前に変えたデータセット (minwage) を使っている

4.2 最低時給の確認

  • 1992年4月にニュージャージーで実施された最低時給の上昇が、実際に最低時給を引き上げたかどうかを確認する
  • 最低時給上昇「前」の賃金 (wage_before) と最低時給上昇「後」の賃金 (wage_after) が $5.05未満のファーストフード店の割合を計算する
minwage %>% 
  group_by(state) %>% 
  summarize(before = mean(wage_before < 5.05),
            after = mean(wage_after < 5.05),
            .groups = "drop") %>% 
  knitr::kable(digits = 4)
state before after
NJ 0.9107 0.0034
PA 0.9403 0.9552
  • NJ での最低時給引き上げ前 (before) には、どちらの州でも 90%以上の労働者の時給が $5.05未満 (91%: NJ, 94%: PA)
  • NJ での最低時給引き上げ後 (after)には
    ・処置があった NJ・・・時給が $5.05未満の Fast Food 店の割合 → 0.34% (0.0034)
    ・処置がなかった PA・・・時給が $5.05未満の Fast Food 店の割合は 95.5% (0.9952)
  • NJ における法律上の最低時給の引き上げは、実際に NJ における最低時給を引き上げた

4.3 単純比較 I： 個体間比較

政策(処置):

  • 1992年4月に NJ の最低賃金が $4.25から$5.05に上げられた

結果：

  • フルタイム労働者の割合

調べたい因果効果：

  • NJ の最低賃金の上昇がフルタイム労働者の割合に与える影響
ニュージャージー (NJ) ：最低賃金の引き上げあり
ペンシルベニア (PA) ：最低賃金の引き上げなし
  • 最低時給の引き上げが雇用にどのような影響を与えたか、単純比較で推定してみる
  • 処置後の NJ と PA のフルタイム労働者の割合(full_prop_after) を単純に比較してみる
minwage %>% 
  group_by(state) %>% 
  summarize(fulltime = mean(full_prop_after),
            .groups = "drop") %>% 
  knitr::kable(digits = 4)
state fulltime
NJ 0.3204
PA 0.2723

単純比較による因果推論結果：

  • 処置後、NJ でのフルタイム労働者の割合は 32% (0.3204)
  • 処置後、PA でのフルタイム労働者の割合は 27% (0.2723)
    ・NJ の最低賃金の上昇がフルタイム労働者の割合を 4.8%ポイント (27.2 → 32) 増やした

問題：

・NJ と PA は、処置「最低賃金の上昇」以外が同質ではない
- このままでは因果推論はできない
→　そもそも NJ の方がもともとフルタイム労働者の割合が高かったのでは？

4.4 単純比較 II：前後比較

  • 処置の前後でフルタイム労働者の割合を単純に比較 する
  • 比較するのは次の 2 つの変数
    ・full_prop_before: 処置の「前」フルタイム労働者の割合
    ・full_prop_after: 処置「後」のフルタイム労働者の割合
minwage %>% 
  filter(state == "NJ") %>% 
  summarize(across(.cols = starts_with("full_"), mean)) %>% 
  knitr::kable(digits = 4)
full_prop_before full_prop_after
0.2965 0.3204

・単純比較による因果推論結果：

  • 処置前、NJ でのフルタイム労働者の割合は 29.65% (0.2965)
  • 処置後、NJ でのフルタイム労働者の割合は 32% (0.3204)
    ・NJ の最低賃金の上昇がフルタイム労働者の割合を 2.4%ポイント (29.65 → 32) 増やした

問題：

  • サービス実施前の NJ と実施後の NJ は同質ではない
  • このままでは因果推論はできない
  • フルタイム労働者の割合を計測している時間が違うから
  • ここでは時間による変化がコントロールされていない：
    →　サービスを実施しなくても、フルタイム労働者の割合は増えたのでは？　　

4.5 差分の差分 (DID)

  • 個体間と処置前後のフルタイム労働者の割合を表にする
  • ここでは次の 2 つの変数を使う
    ・full_prop_before: 処置の「前」フルタイム労働者の割合
    ・full_prop_after: 処置「後」のフルタイム労働者の割合
d_full <- minwage %>% 
  group_by(state) %>%
  summarize(across(.cols = starts_with("full_"), mean),
            .groups = "drop")
knitr::kable(d_full, digits = 3)
state full_prop_before full_prop_after
NJ 0.297 0.320
PA 0.310 0.272
  • ここで示した4つの数字の関係を可視化できる
p_did <- d_full %>% 
  pivot_longer(cols = starts_with("full"),
               names_to = "time",
               names_prefix = "full_prop_",
               values_to = "prop") %>% 
  mutate(time = factor(time, levels = c("before", "after"),
                       labels = c("処置前", "処置後"))) %>% 
  ggplot(aes(x = time, y = prop, group = state)) +
  geom_line(aes(color = state)) +
  geom_point(aes(shape = state)) +
  ylim(0.25, 0.35) +
  labs(x = "", y = "フルタイム労働者の平均割合") +
  scale_color_discrete(name = "") +
  scale_shape_discrete(name = "")
plot(p_did)
  • NJ と PA のフルタイム労働者の割合に平行トレンドが仮定できるなら、差分の差分によって、最低時給上昇の処置効果 \(δ\)（デルタ）を推定することができる。

  • 平行トレンドがあると仮定して、差の差を計算してみる
minwage %>% 
  group_by(state) %>%
  summarize(across(.cols = starts_with("full_"), mean),
            .groups = "drop") %>% 
  mutate(dif_ba = full_prop_after - full_prop_before) %>% 
  with(dif_ba[state == "NJ"] - dif_ba[state == "PA"])
[1] 0.06155831
  • 差分の差分 (DID) を使うと、最低賃金の引き上げは、フルタイム労働者の割合を 6.1%ポイント上昇させる
  • この推定値は、個体間の単純比較 (4.8%ポイント) や前後比較による推定値 (2.4%ポイント) よりも大きい

処置前後のフルタイム労働者の割合 (1992年)

処置前 処置後 差
ニュージャージー 29.7 32 -2.4
ペンシルバニア 31 27.2 3.8
差 -1.3 4.8 6.1

4.6 回帰分析による DID の推定

  • DID による推定値を、回帰分析によって得る方法を考えてみる
  • DID 回帰のために必要なのは、次の変数である
変数名 内容
\(prop\) ：フルタイム労働者の割合（結果変数）
\(D\) ：処置群を表すダミー変数 \(D = 0\): 統制群、\(D = 1\): 処置群
\(P\) ：処置後を表すダミー変数 \(P_t = 0\): 処置前、\(P_t = 1\): 処置後
\(DP\) ：\(D\) と \(P\) の交差項
  • ここまで使ってきたデータ (minwage) は横長 (wide)で 処置前の結果変数 (full_prop_before) と処置後の結果変数 (full_prop_after) の値が異なる列にある形式
  • この 2 つの結果変数は一つの変数 (prop) に統一する必要がある　　
DT::datatable(minwage)%>%
    DT::formatRound(columns=c('full_prop_before', 'full_prop_after'),
                    digits=2)
  • pivot_longer()を使って縦長 (long) に変換し、必要な変数を作る
minwage_long <- minwage %>% 
  dplyr::select(state, starts_with("full_")) %>% 
  pivot_longer(cols = starts_with("full"),
               names_to = "time",
               names_prefix = "full_prop_",
               values_to = "prop") %>% 
  mutate(D = ifelse(state == "NJ", 1, 0), # NJ なら D = 1, PA なら D = 0
         P = ifelse(time == "after", 1, 0)) # 処置後なら P = 1, 処置前なら P = 0
  • 作成したデータフレームを確認
DT::datatable(minwage_long) %>%
    DT::formatRound(columns=c('prop'),
                    digits=2)
  • このデータフレームを使って回帰分析を行う
did_fit00 <- lm(prop ~ D * P, data = minwage_long)
tidy(did_fit00) %>% 
  select(term, estimate) %>% 
  knitr::kable(digits = 2)
term estimate
(Intercept) 0.31
D -0.01
P -0.04
D:P 0.06
  • D:P（DとPの交差項）の係数 = 0.06 が、DIDによる推定値（上で計算した差の差の値と同じ）であることがわかる
結論 ・1992年4月にアメリカのニュージャージー州 (NJ)で最低賃金が引き上げられたため、フルタイム労働者の割合は6.1%ポイント増えた

5. DID の実例（法定飲酒年齢）

  • 法定飲酒年齢が死者数に与える影響: Angrist and Pischke(2015)

Research Question:

5.1 DID 回帰式

  • ここでは差分の差分を利用した回帰分析によって、MLDA の平均処置効果を推定する
  • MLDA の平均処置効果は、上で作成した legal を使うことで、legal の平均処置効果として代替する
  • ここでは次の DID 回帰式を推定する

\[Y_{st} = α + δLEGAL_{st} +{\sum_{k=1}^{50}}{β_kSTATE_{ks}} + {\sum_{j=1971}^{1983}}{γ_jYEAR_{jt}} + e_{st}\]

変数の説明

\(s = {[1, 2, 3,..., 51]}\) : 州（\(s = 51\)は D.C.）
\(t = {[1970, 1972, ... , 1983]}\) : 年（1970年は参照カテゴリ）
\(Y_{st}\) : \(t\) 年の \(s\) 州での 18-20歳の死者数（10万人あたり）
\(STATE_{ks}\) : \(k = s\) の時に 1 になる変数（D.C.を除く50州のダミー）例）CAなら \(STATE_{CA,s} = 1\)、その他の州は 0
\(YEAR_{jt}\) : \(j = t\) の時に 1 になる変数（1970年を除く13個の年ダミー）例）1971年なら \(YEAR_{1971,s} = 1\)、その他の年は 0
\(LEGAL_{st}\) :処置効果（18歳から20歳までの若者のうち、合法的に飲酒できる人の割合）

5.2 データの準備

download.file(url = "http://masteringmetrics.com/wp-content/uploads/2015/01/deaths.dta",
              destfile = "data/deaths.dta")
  • Stata 形式のデータなので、haven::read_dta() で読み込む
MLDA <- read_dta("data/deaths.dta")
  • データフレーム MLDA が含む変数の確認
names(MLDA)
 [1] "year"         "state"        "legal1820"    "dtype"        "agegr"       
 [6] "count"        "pop"          "age"          "legal"        "beertaxa"    
[11] "beerpercap"   "winepercap"   "spiritpercap" "totpercap"    "mrate"
DID 回帰の結果変数: mrate
  • 結果変数 \(Y_{st}\) は mrate
  • mrate は次の 4 種類の死に関するデータを含んでいる
  • 4 種類の死は dtype によって分類されている
すべての死 All deaths dtype = 1
自動車事故による死 Motor vehicle accidents dtype = 2
自殺 Suicide dtype = 3
内蔵疾患による死 All internal causes dtype = 6
  • ここでは「すべての死」(dtype = 1) のみを扱う
    注意：（7. Exercise では他の 3 種類の死に関しても取り扱う）
  • 分析に必要な変数だけを残す
df5 <- MLDA %>% 
  filter(year <= 1983,
         agegr == 2,      # 18-21 years old
         dtype == 1) %>%  # all deaths      
  mutate(state = factor(state),
         year_fct = factor(year))
  • このデータは、州 (state) と年 (year) の組み合わせが1つひとつの行を構成する 州-年パネル (state-year panel)
  • 州の数と観測期間の数を確認してみる
df5 %>% 
  summarize(across(.cols = c(state, year), 
                   .fns = n_distinct))
# A tibble: 1 × 2
  state  year
  <int> <int>
1    51    14
  • 分析の対象となる州は51（50州とワシントンD.C.）
  • 分析の対象となる年は14年（1970-1983）: 1970年は参照カテゴリ
  • 分析に使う結果変数は mrate
mrate： mortality rate: 死亡率（10万人あたりの死者数）
  • mrate の分布を確認する
hist_mrate <- ggplot(df5, aes(x = mrate, y = after_stat(density))) +
  geom_histogram(color = "black", binwidth = 10) +
  labs(x = "死者数（10万人あたり）", y = "確率密度")
plot(hist_mrate)

  • 死亡率の時系列変化を可視化する
ts_mrate <- ggplot(df5, aes(x = year, y = mrate, group = state)) +
  geom_line() +
  labs(x = "年", y = "死者数（10万人あたり）") +
  theme(legend.position = "none")
plot(ts_mrate)

DID 回帰に必要な説明変数
D： 個体の差（処置群と統制群の区別）を表す state: 州
P： 時間（処置・施策が実施される前後の区別）を表す yeat_fct: 年
D×P： 処置を受けた後の観測値であることを示す state:year_fct
  • 個体の差には州を表す state を使う
  • 時間には年を表す year_fct を使う
    （year をそのまま使うと数値として扱われてしまうので、factor 型に変換した year_fct を使う）
  • しかし、処置後の観測値であるかどうかは、これら 2 つの変数を使った交差項 (state:year_fct) では表現できない
その理由：
  • この分析で考える処置（介入）は MLDA の変更であるが、MLDA は18歳, 19歳, 20歳, 21歳のいずれかであるため
  • MLDA が変更されるタイミングは、州によって異なるため
  • 理論的には、複数回の MLDA 変更もあり得るため
解決策：
  • 処置を表す変数 (legal) を作る
    legal: 18歳から20歳までの若者のうち、合法的に飲酒できる人の割合 (0 ≦ legal ≦ 1)

  • 法定飲酒年齢 (MLDA) の引き下げは、legal の値を大きくする
  • 年の途中で MLDA が変更された場合には、その期間に応じて値が調整される
  • legal の値が大きいほど、若者（18-20歳）がアルコールにアクセスしやすい
  • MLDA を引き下げるという処置（施策）を実行すると、legal の値が大きくなる

5.3 DID 回帰

  • ここでは次の DID 回帰式を推定する

\[Y_{st} = α + δLEGAL_{st} +{\sum_{k=1}^{50}}{β_kSTATE_{ks}} + {\sum_{j=1971}^{1983}}{γ_jYEAR_{jt}} + e_{st}\]

  • 上で説明したように、この回帰式を実際に実行するために必要な変数は次のとおりである
\(Y_{st}\)： 結果変数（死亡率：10万人あたりの死者数） mrate
\(LEGAL_{st}\)： 18歳から20歳までの若者のうち、合法的に飲酒できる人の割合 legal
\(STATE_{ks}\)： 個体の差（処置群と統制群の区別）を表す state
\(YEAR_{jt}\)： 時間（処置・施策が実施される前後の区別）を表す year_fct
  • これらの変数の記述統計を示す
summ.stat <- df5 %>% 
  select(mrate, legal, state, year)
 stargazer(as.data.frame(summ.stat), type = "html")
Statistic N Mean St. Dev. Min Max
mrate 714 136.442 39.118 52.474 331.887
legal 714 0.557 0.439 0.000 1.000
year 714 1,976.500 4.034 1,970 1,983 
DT::datatable(summ.stat) %>%
    DT::formatRound(columns=c('mrate', 'legal'),
                    digits=2)
  • DID 回帰してみる
fit_ad0 <- lm(mrate ~ 0 
              + legal 
              + state 
              + year_fct, # year を factor 型に変換した変数 
              data = df5)
tidy(fit_ad0) %>% 
  select(term, estimate, std.error, p.value) %>% 
  filter(term == "legal") %>% 
  knitr::kable(digits = 2)
term estimate std.error p.value
legal 10.8 3.14 0

  • 変数 legal の効果が 10.80 →　 これが平均処置効果

  • 18-20歳の人々がアルコールを摂取できない状況から摂取できる状況に変化すると(legal = 0 → legal = 1)、10万人あたりの死者数が約11人増える

  • つまり、MLDA を下げると、酒を飲めるようになった年齢の人たちのなかで死者が増える

  • Angrist and Pischke (2009) (p.196)の表5.2 の “All deaths” の行の (1) の列にある数値と同じ

  • Angrist and Pischke (2009) (p.196)の表5.2の標準誤差・・・ 4.59
  • ここで得られた標準誤差 (3.14) のほうが小さい
その理由:
  • この分析では時系列データを使う場合、同一州内の観測値が類似するという系列相関 (serial correlation) を考慮していないため
    →　標準誤差を実際よりも過小評価している
    →　不確実性を過小評価している
解決策:

  • 州というクラスタを考慮に入れた標準誤差 (cluster-robust standard error) を計算する
    →　stimatr パッケージの estimatr::lm_robust() を使う　　

  • lm_robust() でクラスタ標準誤差を得るためには、次の 2 つを指定

  • cluster には、クラスタを表す変数を指定 (cluster = state)

  • se_type には、標準誤差を計算する方法を指定

  • 通常は se_type = "CR2" と指定するが、ここでは Angrist and Pischke (2009) と同じ結果を得るために se_type = "stata" と 指定

fit_ad1 <- lm_robust(mrate ~ 0 
                     + legal 
                     + state 
                     + year_fct, 
                     data = df5,
                     clusters = state, 
                     se_type = "stata")
tidy(fit_ad1) %>% 
  select(term, estimate, std.error, p.value) %>% 
  filter(term == "legal") %>% 
  knitr::kable(digits = 2)
term estimate std.error p.value
legal 10.8 4.59 0.02
  • legal の係数の推定値は先ほどとまったく同じ (10.80)
  • 標準誤差が先ほどより大きくなった (4.59)
  • クラスタ化を考慮しない分析では、標準誤差が過小評価される
  • この係数の推定値と標準誤差の値は、Angrist and Pischke (2009) (p.196)の表5.2 の “All deaths” の行の (1) の列にある数値に一致

5.4 平行トレンドの仮定

  • DID 分析では、平行トレンドが仮定されているが、分析対象となる個体と期間の両者が多い場合には、仮定を緩めた分析をすることが可能

  • 「傾きが同じ」という平行トレンドの仮定を「傾きが異なる」に緩める

  • ここでは次の DID 回帰式を推定する

\[Y_{st} = α + δLEGAL_{st} +{\sum_{k=1}^{50}}{β_kSTATE_{ks}} + {\sum_{j=1971}^{1983}}{γ_jYEAR_{jt}} + {\sum_{k=1}^{50}}{θ_k(STATE_{ks}×t)} + e_{st}\]

  • 個体ごとの時間トレンドの違いを説明する変数 \((STATE_{ks}×t)\) を回帰分析に加える
    →　トレンドが完全に平行でなくても、平均置効果を推定することができるようになる
  • ここでは、それぞれの州の時間トレンドは線形（つまり、処置がない場合の死亡者数は、変化なし、単調増加、単調減少のいずれか）であることを仮定
  • Angrist and Pischke (2009) (p.196) の表5.2 の (2) 列では、州ごとの時間トレンドが考慮されている
  • 州ごとの時間トレンドは、state × year（factor 型でなく、数値）で表す
  • コマンドを入力する際の注意：
    state*year と入力すると、sate, year, state:year の 3 つの項が全てモデルに含まれてしまうので、state:year と入力する
fit_ad2 <- lm_robust(mrate ~ 0 
                     + legal 
                     + state  
                     + year_fct 
                     + state:year, 
                     data = df5,
                     clusters = state, 
                     se_type = "stata")
tidy(fit_ad2) %>% 
  select(term, estimate, std.error, p.value) %>% 
  filter(term == "legal") %>% 
  knitr::kable(digits = 2)
term estimate std.error p.value
legal 8.47 5.1 0.1
  • 表5.2 の “All deaths” の行の (2) の列と同じ結果が得られた
  • State trends の項目に Yes と書かれていることが、平行トレンド (parallel trend) を考慮したという意味

5.5 重み付け (weight) の考慮

  • 表5.2 の (3) の列では、州の人口を考慮した重み付き回帰が行われている
  • 人口の多い Texas や California などと、人口の少ない Wyoming や Vermont は異なる
  • lm() 同様、lm_robust() でも weights で重みが指定可能
  • これを加重最小自乗法 (weighted least square: WLS) という
  • 人口は、pop という変数に記録されているので、weights = pop と指定する
fit_ad3 <- lm_robust(mrate ~ 0 
                     + legal 
                     + state 
                     + year_fct, 
                     data = df5, 
                     weights = pop,
                     clusters = state, 
                     se_type = "stata")
tidy(fit_ad3) %>% 
  select(term, estimate, std.error, p.value) %>% 
  filter(term == "legal") %>% 
  knitr::kable(digits = 2)
term estimate std.error p.value
legal 12.41 4.6 0.01

  • 表5.2 の “All deaths” の行の (3) の列と同じ結果が得られた

5.6 平行トレンドの仮定 + 重み付け

  • 平行トレンドの仮定と重み付けを両方同時に考慮した回帰分析も可能である
fit_ad4 <- lm_robust(mrate ~ 0 
                     + legal 
                     + state 
                     + year_fct 
                     + state:year,
                     data = df5, 
                     weights = pop,
                     clusters = state, 
                     se_type = "stata")
tidy(fit_ad4) %>% 
  select(term, estimate, std.error, p.value) %>% 
  filter(term == "legal") %>% 
  knitr::kable(digits = 2)
term estimate std.error p.value
legal 9.65 4.64 0.04
  • 表5.2 の “All deaths” の行の (4) の列と同じ結果が得られた

・4 つの結果のうち、州ごとのトレンドと人口の違いを考慮した最後のモデル (4) から、次のような結論が得られる

結論 法定飲酒年齢を21歳から18歳に引き下げると、18-20歳の死者は10万人あたり約10人増える
・この因果効果は有意水準0.04（4%）で統計的に有意
・10万人あたり10人の死者数増加は、それなりに大きな効果といえそう

6. 平行トレンドのシミュレーション

  • ここでは平行トレンドの仮定が成り立たない場合について、シミュレーションを使って考えてみる

6.1 シミュレーションの設定

  • 50の個体（state, 州：1から50までの整数で区別される）について、1990年から2009年までの値をもつパネルデータを作る
set.seed(2020-07-01)
sim_df <- tibble(year = 1990:2009)
init <- rnorm(50, mean = 100, sd = 30)
trend <- -3 + 1:50 * 0.2 + rnorm(50, mean = 0, sd = 1)
outcome <- sapply(1:50,
                  function(i) init[i] + trend[i] * 0:19)
colnames(outcome) <- str_c("state_", 1:50)
sim_df <- outcome %>% 
  as_tibble() %>% 
  mutate(year = 1990:2009) %>% 
  pivot_longer(cols = starts_with("state"),
               names_to = "state",
               names_prefix = "state_",
               values_to = "y0")
DT::datatable(sim_df) %>%
    DT::formatRound(columns=c('y0'),
                    digits=2)
  • ここで作成したデータセットには次の 3 つの変数が含まれる
state： 州を表す変数: (1,2, ..., 50)
year： 年を表す変数: (1990, 1991, ..., 2009)
y0： 処置がない場合の結果変数
  • 処置がまったくない場合に結果変数 y0 がどのような変化示すか、図示してみる
ggplot(sim_df, aes(x = year, y = y0, group = state)) +
  geom_line()

  • この図から分かるとおり、トレンドが平行ではない
  • ただし、どのトレンドも線形（直線）である
  • 処置を受ける県を決める
  • state の型を調べる
str(sim_df$state)
 chr [1:1000] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" ...
  • state は文字型 (character)
  • 1 から 25 までが処置を受け（処置群 = 1）
  • 26 から 50 までは処置を受けない（統制群 = 0）
sim_df <- sim_df %>% 
  mutate(D = ifelse(state %in% as.character(1:25), 1, 0))
  • 次に、処置を受けるタイミングを決める
  • 話を単純にするため、処置を受けるタイミングは、どの個体も同じ 2000年とする
  • 2000年に処置を受けるので、効果は2001年（2000年と2001年の差）に現れることになる
sim_df <- sim_df %>% 
  mutate(P = ifelse(year > 2000, 1, 0))
  • 処置群の個体に、平均 delta の処置効果を与える
  • ここでは、delta = 10 とする
delta <- 10
effect <- rep(rnorm(50, mean = delta, sd = 2), each = 20)
sim_df <- sim_df %>%
  arrange(state, year) %>% 
  mutate(outcome = case_when(
    D == 0 ~       y0,
    P == 0 ~       y0,
    TRUE         ~ y0 + effect 
  ))
  • これで1回分のシミュレーションデータができた
DT::datatable(sim_df) %>%
    DT::formatRound(columns=c('y0', 'outcome'),
                    digits=2)
  • 結果変数 (outcome) を可視化してみる
  • 処置群と統制群を異なる色で示す
ggplot(sim_df, aes(x = year, y = outcome, group = state)) +
  geom_line(aes(color = as.factor(D))) +
  scale_color_discrete(name = "", labels = c("統制群", "処置群")) +
  labs(x = "年", y = "結果変数") +
  geom_vline(xintercept = 2001, linetype = "dotted") +
  theme_set(theme_classic(base_size = 10,
                        base_family = "HiraginoSans-W3"))

  • トレンドが平行ではなく、統制群のトレンド（トマト色）と比較すると、処置群のトレンド（青色）の方が負の傾きになりやすいことが分かる
  • 2000年に処置を受けるので、効果は2001年（2000年と2001年の差）に現れていることが分かる（黒色の縦の点線）
  • このデータを使って、トレンドをコントロールせずに DID 回帰を実行してみる
fit_s1 <- lm_robust(outcome ~ D * P, 
                    data = sim_df,
                    clusters = state, 
                    se_type = "stata")
tidy(fit_s1) %>% 
  filter(term == "D:P") %>% 
  select(term, estimate, std.error, p.value) %>% 
  knitr::kable(digits = 2)
term estimate std.error p.value
D:P -42.94 4.99 0
  • 当初、設定した平均処置効果は 10 であった
  • しかし、ここで得られた平均処置効果は -42.94
  • 平行トレンドの仮定が成り立たないときに、差の差によって平均処置効果を推定しようとしてもうまくいかない
  • トレンドをコントロールしてみる
  • 県によって異なる線形の時間トレンド state:year をコントロールして、回帰分析を行ってみる
fit_s2 <- lm_robust(outcome ~ D * P 
                    + state:year, 
                    data = sim_df,
                    clusters = state, 
                    se_type = "stata")
tidy(fit_s2) %>% 
  filter(term == "D:P") %>% 
  select(term, estimate, std.error, p.value) %>% 
  knitr::kable(digits = 2)
term estimate std.error p.value
D:P 10.07 0.45 0
  • ここで得られた平均処置効果は 10.07
  • 県ごとのトレンドをコントロールすることにより、設定した値 (10) に近い平均処置効果が推定できた

6.2 シミュレーションの繰り返し実行

  • 上の結果は偶然かもしれない
    →　複数回、シミュレーションを繰り返すとどうなるか確かめてみる
  • 1 回だけのシミュレーションを行う関数を作る
did_trend_sim <- function(delta = 10, 
                          n_units = 50, 
                          n_periods = 20,
                          n_treat = 25, 
                          treat_timing = 10) {
  
  if (n_treat > n_units) stop("n_treat must be smaller than n_units.")
  if (treat_timing >= n_periods | treat_timing < 2) stop("treat_timing must be within time periods.")
  
  sim_df <- tibble(year = 1:n_periods)
  init <- rnorm(n_units, mean = 100, sd = 20)
  trend <-  -n_units/20 + 1:n_units * 0.2 + rnorm(50, mean = 0, sd = 1)
  outcome <- sapply(1:n_units, function(i) init[i] + trend[i] * ((1:n_periods) - 1))
  colnames(outcome) <- str_c("unit_", 1:n_units)
  sim_df <- outcome %>% 
    as_tibble() %>% 
    mutate(time = 1:n_periods) %>% 
    pivot_longer(cols = starts_with("unit"),
                 names_to = "unit",
                 names_prefix = "unit_",
                 values_to = "y0") %>% 
    mutate(D = ifelse(unit %in% as.character(1:n_treat), 1, 0)) %>% 
    mutate(P = ifelse(time > treat_timing, 1, 0)) %>% 
    mutate(outcome = case_when(
      D == 0 ~ y0,
      P == 0 ~ y0,
      TRUE   ~ y0 + rnorm(n(), mean = delta, sd = 2)
    ))
  
  ##  標準誤差は記録しないので、lm() で推定する
  f1 <- lm(outcome ~ D * P, data = sim_df)
  f2 <- lm(outcome ~ D * P + unit:time, data = sim_df)
  return(c(coef(f1)[4], coef(f2)[4]))
}
  • 上の関数を使い、シミュレーションを1000回繰り返す
res <- replicate(1000, did_trend_sim(delta = 10))
  • 結果を図示してみる
res_df <- res %>% 
  t() %>% 
  as_tibble(.name_repair = c("unique"))
names(res_df) <- c("no_trend", "with_trend")

hist_res1 <- ggplot(res_df, aes(x = no_trend, y = after_stat(density))) +
  geom_histogram(color = "lightgrey") +
  geom_vline(xintercept = 10, color = "tomato") +
  labs(x = "推定値", y = "確率密度", title = "トレンドのコントロールなし")

hist_res1

  • トレンドをコントロールしないと、設定値 (10) からかけ離れた値を推定してしまう
hist_res2 <- ggplot(res_df, aes(x = with_trend, y = after_stat(density))) +
  geom_histogram(color = "lightgrey") +
  geom_vline(xintercept = 10, color = "tomato") +
  labs(x = "推定値", y = "確率密度", title = "トレンドのコントロールあり")

hist_res2

  • トレンドをコントロールすると、当初設定した値 (10) のあたりを正しく推定できる

シミュレーションの結論 ・平行トレンドの仮定が成り立たない場合、トレンドを考慮に入れないと正しい推定ができない
・トレンドをコントロールすることによって、平均処置効果の推定が可能になる

・ただし、これは作成したデータのトレンドが線形だった（そして私たちはそれを事前に知っていた）からうまくいっただけ
・非線形のトレンドがある場合、線形トレンドをコントロールしても、推定はうまくいかないことが予想される
・パネルデータの分析では、どのようなトレンドがあるのか見極めることが重要

7. Excercise

Angrist and Pischke (2009) の例について、以下を実行しなさい

Q1:「自動車事故による死」(Motor vehicle accidents: dtype = 2) を使って DID 回帰を実行し、上の表5.2と同じ結果が再現できることを確認しなさい　　

Q2:「自殺」(Suicide: dtype = 3) を使って DID 回帰を実行し、上の表5.2と同じ結果が再現できることを確認しなさい　　

Q3:「内蔵疾患による死」(All internal causes: dtype = 6) を使って DID 回帰を実行し、上の表5.2と同じ結果が再現できることを確認しなさい　　

Reference　　
  • Angrist, Joshua D., and Jörn-Steffen Pischke. 2009. Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton: Princeton University Press.
  • Card, David, and Alan B. Krueger. 1994. “Minimum Wages and Employment: A Case Study of the Fast-Food Industry in New Jersey and Pennsylvania.” American Economic Review 84 (4): 772–93.
  • 矢内勇生「KUT計量経済学応用」
  • 宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」
  • 浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年
  • Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017