Difference-in-differences with individual level panel data

Question

What is the correct way to specify a difference in difference model with individual level panel data?

Here is the setup: Assume that I have individual-level panel data embedded in cities for multiple years and the treatment varies on the city-year level. Formally, let $y_{ist}$ be the outcome for individual $i$ in city $s$ and year $t$ and $D_{st}$ be a dummy for whether the intervention affected city $s$ in year $t$. A typical DiD estimator such as the one outlined in Bertrand et al (2004, p. 250) is based on a simple OLS model with fixed effect terms for city and year:

$$ y_{ist} = A_{s} + B_t + cX_{ist} + \beta D_{st} + \epsilon_{ist} $$

But does that estimator ignore the individual-level panel structure (i.e. multiple observations for each individual within cities)? Does it make sense to extend this model with an individual-level fixed effect term $S_i$? Many DiD applications use repeated cross-section data without the individual-level panel data.

---

Bertrand, Marianne, Esther Duflo, and Sendhil Mullainathan. 2004. “How Much Should We Trust Differences-in-Differences Estimates?” Quarterly Journal of Economics 119(1):249–75.

Answer 1

A nice feature of difference-in-differences (DiD) is actually that you don't need panel data for it. Given that the treatment happens at some sort of level of aggregation (in your case cities), you only need to sample random individuals from the cities before and after the treatment. This allows you to estimate $$ y_{ist} = A_g + B_t + \beta D_{st} + c X_{ist} + \epsilon_{ist} $$ and get the causal effect of the treatment as the expected post-pre outcome difference for the treated minus the expected post-pre outcome difference for the control.

There is a case in which people use individual fixed effects instead of a treatment indicator and this is when we don't have a well-defined level of aggregation at which the treatment occurs. In that case you would estimate $$ y_{it} = \alpha_i + B_t + \beta D_{it} + cX_{it}+\epsilon_{it} $$ where $D_{it}$ is an indicator for the post-treatment period for individuals who received the treatment (for example, a job market program which happens all over the place). For more information on this see these lecture notes by Steve Pischke.

In your setting, adding individual fixed effects should not change anything with respect to the point estimates. The treatment indicator $A_g$ will just be absorbed by the individual fixed effects. However, these fixed effects might soak up some of the residual variance and therefore potentially reduce the standard error of your DiD coefficient.

Here is a code example which shows that this is the case. I use Stata but you can replicate this in the statistical package of your choice. The "individuals" here are actually countries but they are still grouped according to some treatment indicator.

* load the data set (requires an internet connection)
use "http://dss.princeton.edu/training/Panel101.dta"

* generate the time and treatment group indicators and their interaction
gen time = (year>=1994) & !missing(year)
gen treated = (country>4) & !missing(country)
gen did = time*treated

* do the standard DiD regression
reg y_bin time treated did

------------------------------------------------------------------------------
       y_bin |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        time |       .375   .1212795     3.09   0.003     .1328576    .6171424
     treated |   .4166667   .1434998     2.90   0.005       .13016    .7031734
         did |  -.4027778   .1852575    -2.17   0.033    -.7726563   -.0328992
       _cons |         .5   .0939427     5.32   0.000     .3124373    .6875627
------------------------------------------------------------------------------

 * now repeat the same regression but also including country fixed effects
 areg y_bin did time treated, a(country)

------------------------------------------------------------------------------
       y_bin |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        time |       .375    .120084     3.12   0.003     .1348773    .6151227
     treated |          0  (omitted)
         did |  -.4027778   .1834313    -2.20   0.032    -.7695713   -.0359843
       _cons |   .6785714    .070314     9.65   0.000       .53797    .8191729
-------------+----------------------------------------------------------------

So you see that the DiD coefficient remains the same when the individual fixed effects are included (areg is one of the available fixed effects estimation commands in Stata). The standard errors are slightly tighter and our original treatment indicator was absorbed by the individual fixed effects and therefore dropped in the regression.

In response to the comment
I mentioned the Pischke example to show when people use individual fixed effects rather than a treatment group indicator. Your setting has a well defined group structure so the way you have written your model that's perfectly fine. Standard errors should be clustered at the city level, i.e. the level of aggregation at which the treatment occurs (I haven't done this in the example code but in DiD settings the standard errors need to be corrected as demonstrated by the Bertrand et al paper).

Regarding the movers, they don't have much of a role to play here. The treatment indicator $D_{st}$ is equal to 1 for people who live in a treated city $s$ in the post-treatment period $t$. To compute the DiD coefficient, we actually just need to compute four conditional expectations, namely $$ c = \left[ E(y_{ist}|s=1,t=1) - E(y_{ist}|s=1,t=0)\right] - \left[ E(y_{ist}|s=0,t=1) - E(y_{ist}|s=0,t=0)\right] $$

So if you have 4 post-treatment periods for an individual who lives in a treated city for the first two, and then moves to a control city for the remaining two periods, the first two of those observations will be used in the computation of $E(y_{ist}|s=1,t=1)$ and the last two in $E(y_{ist}|s=0,t=1)$. To make it clear why identification comes from the group differences over time and not from the movers you can visualize this with a simple graph. Suppose the change in the outcome is truly only because of the treatment and that it has a contemporaneous effect. If we have an individual who lives in a treated city after the treatment starts but then moves to a control city, their outcome should go back to what it was before they were treated. This is shown in the stylized graph below.

You might still want to think about movers for other reasons though. For instance, if the treatment has a lasting effect (i.e. it still affects the outcome even though the individual has moved)

Answer 2

The short answer is that using fixed effect at the unit or at the treated group level does not change estimation, only inference. Typically, using unit fixed-effect will absorb more variation, and hence you will have smaller standard errors.

Whether the units are in the same aggregated group does not change this result (it only changes how you define your treated group level, and the fact that you need panel instead of repeated cross-sections). The same applies to the time fixed effects: you can either use say year fixed effects, or just a FE for pre and post.

Note however that the equivalence only holds when there is no covariate X. As soon as you have Xs, results are different whether you use unit or group fixed effects.

The example below compares the 5 estimators. Estimators are:

1. OLS
2. FE with FE: (treated) group by (treatment) period
3. FE with FE: units by (treatment) period
4. FE with FE: (treated) group by year
5. FE with FE: units by year

                      OLS        FE period-group  FE period-unit  FE year-group  FE year-unit
---------------------------------------------------------------------------------------------
did                   -0.4028 *  -0.4028 *        -0.4028 *       -0.4028 *      -0.4028 *   
                      (0.1853)   (0.1853)         (0.1834)        (0.1758)       (0.1727)    
---------------------------------------------------------------------------------------------
Num. obs.             70         70               70              70             70          
Num. groups: time                 2                2                                         
Num. groups: treated              2                                2                         
Num. groups: country                               7                              7          
Num. groups: year                                                 10             10          
=============================================================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

And the table below show the 3 first regressions, this time with covariates:

===================================================================
                      OLS x     FE period-group X  FE period-unit X
-------------------------------------------------------------------
did                   -0.407 *  -0.407 *           -0.460 *        
                      (0.189)   (0.189)            (0.187)         
x1                     0.018     0.018              0.220          
                      (0.104)   (0.104)            (0.165)         
-------------------------------------------------------------------
Num. obs.             70        70                 70              
Num. groups: time                2                  2              
Num. groups: treated             2                                 
Num. groups: country                                7              
===================================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Code:

library(tidyverse)
library(haven)
library(texreg)
library(lfe)


dat <- read_dta("http://dss.princeton.edu/training/Panel101.dta")  %>% 
  mutate(time = (year>=1994) ,
         treated = (country>4),
         did = time*treated)

reg_ols <- lm(y_bin~ time+ treated+ did, data = dat)
reg_fe_a <- felm(y_bin~ did | time+ treated, data = dat)
reg_fe_b <- felm(y_bin~ did | time+ country, data = dat)
reg_fe_c <- felm(y_bin~ did | year+ treated, data = dat)
reg_fe_d <- felm(y_bin~ did | year+ country, data = dat)

reg_ols_x <- update(reg_ols, .~.+x1)
reg_fe_a_x <- update(reg_fe_a, .~.+x1)
reg_fe_b_x <- update(reg_fe_b, .~.+x1)


screenreg(list(reg_ols, reg_fe_a, reg_fe_b, reg_fe_c, reg_fe_d), 
          omit.coef = "time|treated|Intercept", digits=4, 
          include.rsquared = FALSE, include.adjrs = FALSE, include.rmse = FALSE,
          custom.model.names = c("OLS", "FE period-group", "FE period-unit", "FE year-group", "FE year-unit"))
#> 
#> =============================================================================================
#>                       OLS        FE period-group  FE period-unit  FE year-group  FE year-unit
#> ---------------------------------------------------------------------------------------------
#> did                   -0.4028 *  -0.4028 *        -0.4028 *       -0.4028 *      -0.4028 *   
#>                       (0.1853)   (0.1853)         (0.1834)        (0.1758)       (0.1727)    
#> ---------------------------------------------------------------------------------------------
#> Num. obs.             70         70               70              70             70          
#> Num. groups: time                 2                2                                         
#> Num. groups: treated              2                                2                         
#> Num. groups: country                               7                              7          
#> Num. groups: year                                                 10             10          
#> =============================================================================================
#> *** p < 0.001; ** p < 0.01; * p < 0.05

screenreg(list(reg_ols_x, reg_fe_a_x, reg_fe_b_x), 
          omit.coef = "time|treated|Inter", digits=3, 
          include.rsquared = FALSE, include.adjrs = FALSE, include.rmse = FALSE,
          custom.model.names = c("OLS x", "FE period-group X", "FE period-unit X"))
#> 
#> ===================================================================
#>                       OLS x     FE period-group X  FE period-unit X
#> -------------------------------------------------------------------
#> did                   -0.407 *  -0.407 *           -0.460 *        
#>                       (0.189)   (0.189)            (0.187)         
#> x1                     0.018     0.018              0.220          
#>                       (0.104)   (0.104)            (0.165)         
#> -------------------------------------------------------------------
#> Num. obs.             70        70                 70              
#> Num. groups: time                2                  2              
#> Num. groups: treated             2                                 
#> Num. groups: country                                7              
#> ===================================================================
#> *** p < 0.001; ** p < 0.01; * p < 0.05

^{Created on 2021-04-01 by the reprex package (v1.0.0)}

Answer 3

As a complement to the other answers: when you use panel data, individual are not fully randomly sampled - typically, there might be non-random attrition between waves of the survey.

It was shown by Lechner, Rodriguez-Planas and Kranz (2016) that this can lead to substantial differences between the simple DiD model, and the one using fixed-effects. This might lead to prefer the fixed-effects specifications.

Reference: Lechner, Michael, Nuria Rodriguez-Planas, and Daniel Fernández Kranz. "Difference-in-difference estimation by FE and OLS when there is panel non-response." Journal of Applied Statistics 43.11 (2016): 2044-2052. https://www.tandfonline.com/doi/abs/10.1080/02664763.2015.1126240