1 本講義の目的

lm_robust()を用いた差分の差分法(Difference-in-Differences)の推定の基礎を学ぶ。
データの前処理やグラフ作成の復習や練習をする。

2 パッケージの読み込み

今回の分析で使うパッケージを読み込む。

pacman::p_load(tidyverse, 
               estimatr, # lm_robustによるパネルデータ分析(メイン)
               stargazer, # lmやplmのときの結果表
               texreg, # lmやplmやlm_robustのときの結果表
               plm) # plmによるパネルデータ分析(参考)

3 データの生成

パネルデータ分析の応用として、差分の差分法（difference in differences: DID）法がある。

詳しい解説は参考文献などに譲るが、ここではもっとも単純なDID（二群・二期間DID)の推定を行う。

これまでの事例と同様、能力が均等に分布した個体の所得をアウトカム変数とした分析を行う。処置変数として、職業訓練(training)を想定する。

まず、能力が1~100まで均等に分布した個体を1万人生成させる。

# 事前準備 --------------------

# 乱数の種を固定
set.seed(0)

# データの生成 ----------------
n <- 10000
# 能力は0から100まで均等に分布
ability <- runif(n, min = 0, max = 100)

# IDとabilityをデータフレームに格納する
df <- tibble(ID = 1:n, ability)

二期間DIDなので、１年目と２年目のデータ生成過程（DGP: Data Generating Process）をそれぞれ作成する。

1年目は、誰も職業訓練を受けず、所得は、定数項、能力、誤差項のみによって決まるものとする。

# 1年目

## year = 0を加える
df0 <- df %>% mutate(year = 0) 

##　職業訓練フラグ
##　一年目は誰も職業訓練をうけない
df0 <- df0 %>% mutate(training = 0)

#所得 incomeを加える
df0 <- df0 %>% mutate(income = 200 + 10*ability + 
                        rnorm(n, mean = 0, sd = 50))

2年目は、能力を部分的に反映した適正試験が 180 点以上であれば職業訓練(training)を受けられるとし、職業訓練を受けた個体は500万円、所得が上乗せされるものとする。ここで、職業訓練を受けた個体が処置群、受けなかった個体が対照群となる。

また年ショックとして、平均100、標準偏差25の正規分布\(N(\mu=100、\sigma=25)\)に従う確率変数をすべての個体に加える。

# 2年目

## year = 1　を加える
df1 <- df %>% mutate(year = 1)

## 職業訓練フラグ
## 条件：能力を部分的に反映した適正試験が 180 点以上であれば職業訓練を受けられる
## 10%くらいが該当するように恣意的に設定
## abilityとscoreの関係は非線形なものとする
## 今回もmutate()とcase_when()を使用して作成

# 適正試験の点数 score
df1 <- df1 %>% mutate(score = 30 * log10(ability) + 
                        rnorm(n, mean = 115, sd = 10))

# 職業訓練ダミー training (score 180点以上=1)
df1 <- df1 %>% mutate(training = case_when(score >= 180 ~ 1,
                                           TRUE ~ 0))

#所得 income: 年ショックとしてn(100,25)も追加
df1 <- df1 %>% mutate(income = 200 + 10*ability + 500*training + 
                        rnorm(n, mean = 100, sd = 25) + 
                        rnorm(n, mean = 0, sd = 50))

最後に1年目と2年目のデータを結合して1つのデータフレームにする。

# 1年目と2年目をくっつける。
df_binded <- bind_rows(df0, df1)
head(df_binded)

## # A tibble: 6 × 6
##      ID ability  year training income score
##   <int>   <dbl> <dbl>    <dbl>  <dbl> <dbl>
## 1     1    89.7     0        0  1116.    NA
## 2     2    26.6     0        0   488.    NA
## 3     3    37.2     0        0   511.    NA
## 4     4    57.3     0        0   717.    NA
## 5     5    90.8     0        0  1158.    NA
## 6     6    20.2     0        0   310.    NA

4 差分の差分法(DID)推定

4.1 グラフを用いた説明

上記のデータを用いてDID推定を実際に行う前に、グラフを使ってDID推定を説明しよう。

まず、処置前と処置後、処置群と対照群について、アウトカム変数(income)の４つの平均値を作成する。

今回の場合、処置前と処置後を区別する時間ダミー変数(year)はすでにあるが、処置群と対照群を区別する処置群ダミー変数はないため、作成する必要がある。

なお、trainingは処置群ダミーではなく処置ダミーであり、\(t=0\)のときには全員training\(=0\)であることに注意。

#処置群ダミーの作成

##　処置群の抽出
df_treated <- df_binded %>% filter(training == 1) %>% dplyr::select(ID)

##　処置群ダミーの作成
df_treated <- df_treated %>% mutate(treated = 1)

## df_panelにdf_treatedを結合: dplyr::left_join(a, b, by = "x1") bをaに対応付け、結合する(aが全て残る)
df_panel <- left_join(df_binded, df_treated, by ="ID")

# treated(処置群ダミー）のNAを0をreplace関数を使って置き換える
df_panel <- df_panel %>% mutate(treated = replace(treated, which(is.na(treated)), 0))

#以下のような方法など、他の置き換える方法がある。
# df_panel <- df_panel  %>% mutate_at(vars(treated), funs(ifelse(is.na(.),0,.)))
# df_panel$treated[is.na(df_panel$treated)] <- 0

次に、treatedとyearのグループごとに、アウトカム変数(income)の４つの平均値を作成する。

# summarize(.by)でグループごとに平均値を計算
df_panel_mean <- df_panel %>% summarize(
  income = mean(income),
  .by = c("treated", "year")  # グループごとにmean(income)を計算
)

（なお、別の書き方としてgroup_by()してからsummarize()する方法もある）

# （参考） group_byでグループ分けし、summariseでグループごとに平均値を計算
df_panel_mean <- df_panel %>% 
  group_by(treated, year) %>% # グループわけ
  summarize(income = mean(income), .groups = "keep") %>%  # グループごとに平均値を計算
  ungroup() # グループ分けを解除

この４つの平均値を図示すると、以下のようになる。

# 処置群と対照群のラベルづけ
df_panel_mean <- df_panel_mean %>% mutate(treated_label =
                      case_when(treated == 1 ~ "Treated", 
                                treated == 0 ~ "Control"))

# A, B, C, Dのラベルづけ（図に挿入）
df_panel_mean <- df_panel_mean %>% mutate(point_label =
                      case_when(treated == 1 & year == 1 ~ "A", 
                                treated == 1 & year == 0 ~ "B", 
                                treated == 0 & year == 1 ~ "C",
                                treated == 0 & year == 0 ~ "D"))
# 図の作成
g <- df_panel_mean %>% ggplot(aes(x = year, 
                                  y = income, 
                                  color = treated_label)) +
  geom_point() + # 散布図
  geom_line() +  # 直線
  labs(title = "Income by year and group") + # タイトル
  scale_x_continuous(breaks = c(0,1)) +　# X軸は0と1のみにする
  geom_text(aes(label = point_label), vjust = -0.45) #A,B,C,Dの文字を点の0.45上に表示

# 図の表示  
g

ここで、職業訓練(training)の平均効果は、

\[ \hat{\beta} = (Y^A-Y^B) - (Y^C-Y^D) \]

と推定できる。つまり、処置群(Treated)におけるincomeの平均値の伸び（=職業訓練効果+年ショック）から、対照群(Control)におけるincomeの平均値の伸び（=年ショック）を引いたものが、職業訓練の平均効果と推定される。

4.2 回帰モデルを用いたDID推定

4.2.1 回帰モデル

上記のような二群(処置群と対照群）・二期間(\(t=0,1\))のDID推定は、以下の回帰式を用いて推定できる。

\[ Y_{it} = \mu + \theta D_i + \pi T_t + \beta D_iT_t + \varepsilon_{it} \tag{1} \] ここで、\(D_i\)は処置群ダミー、\(T_t\)は時間ダミー(\(t=0\)だと0、\(t=1\)だと1)であり、DIDパラメータ（処置群に対する平均処置効果:ATET)は\(\beta\)である。

また、この式の二年間の差分（一階の階差）をとると、以下のように表すこともできる。

\[ \Delta Y_i = \pi + \beta D_i + \Delta \varepsilon_{it} \tag{2} \]

(2)式の\(\Delta Y_i\)はincomeの一階の階差変数であり、以下のように作成できる。

# incomeのラグ変数および一階の階差変数の作成
# dplyr::lagとplm::lagの二つがあるためか、dplyr::をつけないとうまくいかなない
df_panel <- df_panel %>% 
  group_by(ID) %>% # IDでグループ分け
  mutate(
    lag_income = dplyr::lag(income, order_by = year), # lag変数の作成
    diff_income = income - lag_income # 一階の階差変数の作成
  ) %>% 
  ungroup() # グループ分けの解除

4.2.2 lmを用いた推定

lmを用いて、(1)式および(2)式に基づいてDID推定を行う。ただし標準誤差は均一分散を仮定して算出されている点に注意する。

# DID regression

# pooled OLS:(1)式
DID_pooled <- lm(income ~ treated + year + treated:year, 
                 data = df_panel)

# first-differenced (一階の階差) OLS: (2)式
DID_FD <- lm(diff_income ~ treated, 
             data = df_panel)

# stargazer
stargazer(DID_pooled, 
          DID_FD, 
          type = "text", 
          digits = 3, 
          df = FALSE,
          column.labels = c("DID (Pooled data)","DID (FD data)"),
          model.names = F, 
          model.numbers = F)

## 
## ===================================================
##                           Dependent variable:      
##                     -------------------------------
##                          income        diff_income 
##                     DID (Pooled data) DID (FD data)
## ---------------------------------------------------
## treated                298.569***      503.117***  
##                          (9.016)         (2.419)   
##                                                    
## year                    98.896***                  
##                          (4.206)                   
##                                                    
## treated:year           503.117***                  
##                         (12.750)                   
##                                                    
## Constant               668.077***       98.896***  
##                          (2.974)         (0.798)   
##                                                    
## ---------------------------------------------------
## Observations             20,000          10,000    
## R2                        0.344           0.812    
## Adjusted R2               0.344           0.812    
## Residual Std. Error      280.734         75.334    
## F Statistic           3,500.464***    43,247.750***
## ===================================================
## Note:                   *p<0.1; **p<0.05; ***p<0.01

DID推定値は概ね500近傍となっており、(1)式と(2)式で完全に一致する。一方、DID推定量の標準誤差が(1)式と(2)式で大きく異なるのは、(1)式では処置群全体の固定効果\(\theta\)を除去しているのに対し、(2)式では、事実上、個体ごとの固定効果\(\theta_i\)を除去しているためである。

4.2.3 lm_robustを用いた推定（メイン）

パネルデータであれば、(1)式の処置群ダミー\(\theta D_i\)の代わりに、個体ダミーを用いてもよい。（ただし、繰り返しクロスセクションデータを用いたDID推定では個体ダミーは用いることはできない点は留意が必要である。）

この場合、個体ダミーを\(n-1\)個用いて推定してもよいが、lm_robust関数で個体固定効果モデルを推定してもいい。

# DID with a two-by-two DID model (オーソドックスなDIDモデル)
DID_robust1 <- lm_robust(income ~ treated + year + treated:year,
                         data = df_panel,
                         clusters = ID,
                         se_type = "stata")

# DID with two-way FE model (年ダミー変数は活用し、個体固定効果はlm_robustで指定） 
DID_robust2 <- lm_robust(income ~ treated:year + year,
                         fixed_effects = ~ ID, #個体固定効果
                         data = df_panel,
                         clusters = ID,
                         se_type = "stata")

# DID with two-way FE model（個体固定効果と年固定効果をlm_robustで指定）
DID_robust3 <- lm_robust(income ~ treated:year,
                         fixed_effects = ~ ID + year, #個体固定効果と年固定効果
                         data = df_panel,
                         clusters = ID,
                         se_type = "stata")

#結果表の出力
screenreg(list(DID_robust1, 
               DID_robust2, 
               DID_robust3), 
          custom.model.names = c("Pooled", 
                                 "Oneway", 
                                 "Twoway"),
          include.ci = FALSE, 
          digits = 3)

## 
## =========================================================
##               Pooled         Oneway         Twoway       
## ---------------------------------------------------------
## (Intercept)     668.077 ***                              
##                  (3.062)                                 
## treated         298.569 ***                              
##                  (6.512)                                 
## year             98.896 ***     98.896 ***               
##                  (0.797)        (1.128)                  
## treated:year    503.117 ***    503.117 ***    503.117 ***
##                  (2.435)        (3.443)        (3.443)   
## ---------------------------------------------------------
## R^2               0.344          0.988          0.988    
## Adj. R^2          0.344          0.976          0.976    
## Num. obs.     20000          20000          20000        
## RMSE            280.734         53.269         53.269    
## N Clusters    10000          10000          10000        
## =========================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05

4.2.4 plmを用いた推定(参考)

plm関数で個体固定効果モデルを推定してもよいが、標準誤差は均一分散を仮定したものである点は注意が必要。

# df_panelのデータをpanelとして認識
df_panel <- pdata.frame(df_panel, 
                        index = c("ID", "year")) # 個体=firm、時間=yearと認識させたdf

# DID regression

# DID with two-way FE model (年ダミーは活用） 
DID_plm1 <- plm(income ~ treated:year + year, 
                data = df_panel, 
                model = "within")

# DID with two-way FE model（年ダミーも入れない）
DID_plm2 <- plm(income ~ treated:year, 
                data = df_panel, 
                model = "within", 
                effect = "twoways")

# stargazer
stargazer(DID_plm1, 
          DID_plm2, 
          type = "text", 
          digits = 3, 
          df = FALSE,
          column.labels = c("DID (onewayFE)","DID (twowayFE)"),
          model.names = F, 
          model.numbers = F)

## 
## ===========================================
##                    Dependent variable:     
##               -----------------------------
##                          income            
##               DID (onewayFE) DID (twowayFE)
## -------------------------------------------
## year1           98.896***                  
##                  (0.798)                   
##                                            
## treated:year0  -503.117***    -503.117***  
##                  (2.419)        (2.419)    
##                                            
## -------------------------------------------
## Observations      20,000         20,000    
## R2                0.895          0.812     
## Adjusted R2       0.789          0.624     
## F Statistic   42,419.510***  43,247.750*** 
## ===========================================
## Note:           *p<0.1; **p<0.05; ***p<0.01

5 参考文献

安藤道人「計量経済学２」あるいは「計量経済特論２」の「差分の差分法（差の差法）」の講義資料
計量経済学応用差分の差分法 (矢内勇生氏の講義ノート)
Chapter 18 Difference-in-Differences (Nick Huntington-Klein (2021)“The Effect: An Introuction to Research Design and Causality”)
Chapter 8 Difference-in-Differences (Scott Cunningham (2021) “Causal Inference: The Mixtape)

Rで学ぶ計量経済学と機械学習 9

計量経済学5：差分の差分(DID)法

安藤道人（立教大学）　三田匡能 (株式会社 GA technologies)

2023-12-09

1 本講義の目的

2 パッケージの読み込み

3 データの生成

4 差分の差分法(DID)推定

4.1 グラフを用いた説明

4.2 回帰モデルを用いたDID推定

4.2.1 回帰モデル

4.2.2 lmを用いた推定

4.2.3 lm_robustを用いた推定（メイン）

4.2.4 plmを用いた推定(参考)

5 参考文献

Rで学ぶ計量経済学と機械学習 9 計量経済学5：差分の差分(DID)法

安藤道人（立教大学） 三田匡能 (株式会社 GA technologies)

2023-12-09

1 本講義の目的

2 パッケージの読み込み

3 データの生成

4 差分の差分法(DID)推定

4.1 グラフを用いた説明

4.2 回帰モデルを用いたDID推定

4.2.1 回帰モデル

4.2.2 lmを用いた推定

4.2.3 lm_robustを用いた推定（メイン）

4.2.4 plmを用いた推定(参考)

5 参考文献

Rで学ぶ計量経済学と機械学習 9

計量経済学5：差分の差分(DID)法

安藤道人（立教大学）　三田匡能 (株式会社 GA technologies)