Interpreting Coefficients in Poisson JWDID Models

[Frank Huang]

Dear all,

I am currently using a Poisson model to estimate a staggered DID. I first estimated the model using a Poisson two-way fixed effects (TWFE) approach (ssc install ppmlhdfe), where the estimated coefficient is typically interpreted as the percentage change in the outcome variable. However, I also use jwdid as a robustness check (ssc install jwdid), and I am somewhat uncertain about how the estimated coefficient from jwdid should be interpreted.

Below is a simple example. I estimate the effects on both the original outcome (emp) and a rescaled outcome equal to ten times the original variable (emp10) using both Poisson TWFE and jwdid. The Poisson TWFE model produces exactly the same coefficient estimates for emp and emp10, whereas the jwdid coefficient for emp10 is exactly ten times the coefficient for emp. This makes me wonder whether the jwdid coefficient should instead be interpreted as a change in the outcome level, rather than a percentage change.

*ssc install ppmlhdfe
*ssc install jwdid

ssc install frause
frause mpdta.dta, clear
gen emp = exp(lemp)
gen emp10 = 10*emp

**********************************************
* staggered poisson DID
**********************************************
* TWFE
gen post = year >= first_treat
gen treat_post = treat*post 
ppmlhdfe emp treat_post, a(countyreal year) vce(cl countyreal)
ppmlhdfe emp10 treat_post, a(countyreal year) vce(cl countyreal)

* JWDID
jwdid emp, ivar(countyreal) tvar(year) gvar(first_treat) method(ppmlhdfe)
estat simple
jwdid emp10, ivar(countyreal) tvar(year) gvar(first_treat) method(ppmlhdfe)
estat simple

[Andrew Musau]

In a Poisson DID with a log link, the coefficient on \(\text{Post} \times \text{Treated}\) is on the log scale, so it is not itself the ATET in levels.

Suppose your model is

\[
E[Y_{it}\mid X_{it}]
=
\exp\left(
\alpha
+ \beta \,\text{Post}_t
+ \gamma \,\text{Treated}_i
+ \delta (\text{Post}_t \times \text{Treated}_i)
+ X'_{it}\theta
\right)
\]

Then, \(\delta\) is the DID effect on the log conditional mean scale, \(\exp(\delta)\) is the multiplicative treatment effect
(an incidence rate ratio).

So,

\[
\exp(\delta)
\]

is the factor by which the expected outcome changes for the treated group after treatment, relative to the counterfactual trend. Therefore, in a Poisson model, the ATET in levels depends on covariates and baseline expected counts. The level effect generally is:

\[
E[Y \mid D=1,\text{Post}=1]
\times
\left(1-\exp(-\delta)\right)
\]

or equivalently the treated potential mean under treatment minus the treated counterfactual mean without treatment. If I computed this after ppmlhdfe, it will be scale dependent as with estat simple after jwdid.

*ssc install ppmlhdfe
*ssc install frause
frause mpdta.dta, clear
gen emp = exp(lemp)
gen emp10 = 10*emp

**********************************************
* staggered poisson DID
**********************************************
* TWFE
gen post = year >= first_treat
gen treat_post = treat*post
ppmlhdfe emp treat_post, a(countyreal year) vce(cl countyreal) d
predict muhat, mu
gen mu0 = muhat / exp(_b[treat_post]) if treat_post==1
gen te = muhat - mu0 if treat_post==1
sum te

ppmlhdfe emp10 treat_post, a(countyreal year) vce(cl countyreal) d
predict muhat10, mu
gen mu0_10 = muhat10/ exp(_b[treat_post]) if treat_post==1
gen te10 = muhat10 - mu0_10 if treat_post==1
sum te10

Res.:

. ppmlhdfe emp treat_post, a(countyreal year) vce(cl countyreal) d
Iteration 1:   deviance = 2.3606e+05  eps = .         iters = 3    tol = 1.0e-04  min(eta) =  -2.59  P  
Iteration 2:   deviance = 3.6946e+04  eps = 5.39e+00  iters = 2    tol = 1.0e-04  min(eta) =  -3.56      
Iteration 3:   deviance = 1.6010e+04  eps = 1.31e+00  iters = 2    tol = 1.0e-04  min(eta) =  -4.47      
Iteration 4:   deviance = 1.3530e+04  eps = 1.83e-01  iters = 2    tol = 1.0e-04  min(eta) =  -5.24      
Iteration 5:   deviance = 1.3326e+04  eps = 1.53e-02  iters = 2    tol = 1.0e-04  min(eta) =  -5.74      
Iteration 6:   deviance = 1.3319e+04  eps = 5.13e-04  iters = 2    tol = 1.0e-04  min(eta) =  -5.92      
Iteration 7:   deviance = 1.3319e+04  eps = 2.39e-06  iters = 2    tol = 1.0e-04  min(eta) =  -5.93      
Iteration 8:   deviance = 1.3319e+04  eps = 1.31e-10  iters = 2    tol = 1.0e-05  min(eta) =  -5.93   S O
------------------------------------------------------------------------------------------------------------
(legend: p: exact partial-out   s: exact solver   h: step-halving   o: epsilon below tolerance)
Converged in 8 iterations and 17 HDFE sub-iterations (tol = 1.0e-08)

HDFE PPML regression                              No. of obs      =      2,500
Absorbing 2 HDFE groups                           Residual df     =        499
Statistics robust to heteroskedasticity           Wald chi2(1)    =       0.56
Deviance             =     13318.87               Prob > chi2     =     0.4558
Log pseudolikelihood = -16174.35906               Pseudo R2       =     0.9947

Number of clusters (countyreal)=       500
                           (Std. err. adjusted for 500 clusters in countyreal)
------------------------------------------------------------------------------
             |               Robust
         emp | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
  treat_post |  -.0146835   .0196899    -0.75   0.456     -.053275     .023908
       _cons |   8.109238   .0029806  2720.67   0.000     8.103396     8.11508
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
  countyreal |       500         500           0    *|
        year |         5           1           4     |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation

.
. predict muhat, mu

.
. gen mu0 = muhat / exp(_b[treat_post]) if treat_post==1
(2,209 missing values generated)

.
. gen te = muhat - mu0 if treat_post==1
(2,209 missing values generated)

.
. sum te

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
          te |        291   -20.16173    40.36058   -312.049  -.1178669

.
.
.
. ppmlhdfe emp10 treat_post, a(countyreal year) vce(cl countyreal) d
Iteration 1:   deviance = 2.3606e+06  eps = .         iters = 3    tol = 1.0e-04  min(eta) =  -2.59  P  
Iteration 2:   deviance = 3.6946e+05  eps = 5.39e+00  iters = 2    tol = 1.0e-04  min(eta) =  -3.56      
Iteration 3:   deviance = 1.6010e+05  eps = 1.31e+00  iters = 2    tol = 1.0e-04  min(eta) =  -4.47      
Iteration 4:   deviance = 1.3530e+05  eps = 1.83e-01  iters = 2    tol = 1.0e-04  min(eta) =  -5.24      
Iteration 5:   deviance = 1.3326e+05  eps = 1.53e-02  iters = 2    tol = 1.0e-04  min(eta) =  -5.74      
Iteration 6:   deviance = 1.3319e+05  eps = 5.13e-04  iters = 2    tol = 1.0e-04  min(eta) =  -5.92      
Iteration 7:   deviance = 1.3319e+05  eps = 2.39e-06  iters = 2    tol = 1.0e-04  min(eta) =  -5.93      
Iteration 8:   deviance = 1.3319e+05  eps = 1.31e-10  iters = 2    tol = 1.0e-05  min(eta) =  -5.93   S O
------------------------------------------------------------------------------------------------------------
(legend: p: exact partial-out   s: exact solver   h: step-halving   o: epsilon below tolerance)
Converged in 8 iterations and 17 HDFE sub-iterations (tol = 1.0e-08)

HDFE PPML regression                              No. of obs      =      2,500
Absorbing 2 HDFE groups                           Residual df     =        499
Statistics robust to heteroskedasticity           Wald chi2(1)    =       0.56
Deviance             =     133188.7               Prob > chi2     =     0.4558
Log pseudolikelihood = -78985.76637               Pseudo R2       =     0.9974

Number of clusters (countyreal)=       500
                           (Std. err. adjusted for 500 clusters in countyreal)
------------------------------------------------------------------------------
             |               Robust
       emp10 | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
  treat_post |  -.0146835   .0196899    -0.75   0.456     -.053275     .023908
       _cons |   10.41182   .0029806  3493.19   0.000     10.40598    10.41766
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
  countyreal |       500         500           0    *|
        year |         5           1           4     |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation

.
. predict muhat10, mu

.
. gen mu0_10 = muhat10/ exp(_b[treat_post]) if treat_post==1
(2,209 missing values generated)

.
. gen te10 = muhat10 - mu0_10 if treat_post==1
(2,209 missing values generated)

.
. sum te10

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
        te10 |        291   -201.6173    403.6057   -3120.49  -1.178671

. 

[Jeff Wooldridge]

This is a bit subtle. The jwdid command reports the coefficients in the original output as the difference in the log of the mean. These coefficients are obtained using ppml, and so you should be able to reproduce them "by hand." But when you aggregate them by exposure time using estat even or obtain a single weighted effect using estat simple, these are average partial effects on the level. So they will not be comparable to using ppml directly.

I actually wish the command allowed the user to choose whether the effect is on the linear index -- so, as a proportionate effect -- or on the level. This is true for logit, flogit, poisson. So the choices would be scale(index) and scale(response), or something like that.

But in the example you gave, the only reason you're getting a difference is because you use estat simple, correct?

[Frank Huang]

Dear Andrew Musau and Jeff Wooldridge :

Thank you both very much for the clear and helpful clarification regarding the interpretation of the jwdid results.
The original Poisson coefficient estimates reported by jwdid are actually identical between emp and emp10. The difference only comes from the post-estimation results from commands such as estat simple (or estat event).

[George Ford]

a few changes to jwdid_estat makes this possible, I believe. worth checking.

frause mpdta, clear
gen emp = exp(lemp)
jwdid emp, ivar(countyreal) tvar(year) gvar(first_treat) method(ppmlhdfe)
estat simple                    // level ATT
estat simple , scale(response)  // level ATT
estat simple , scale(index)     // log-scale ATT


jwdid emp, ivar(countyreal) tvar(year) gvar(first_treat) method(ppmlhdfe)


estat simple                    // level ATT
------------------------------------------------------------------------------
             |            Delta-method
             | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      simple |  -28.91383    20.5633    -1.41   0.160    -69.21716     11.3895
------------------------------------------------------------------------------

estat simple , scale(response)  // level ATT
------------------------------------------------------------------------------
             |            Delta-method
             | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      simple |  -28.91383    20.5633    -1.41   0.160    -69.21716     11.3895
------------------------------------------------------------------------------

estat simple , scale(index)     // log-scale ATT
------------------------------------------------------------------------------
             |            Delta-method
             | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      simple |  -.0286905   .0130213    -2.20   0.028    -.0542117   -.0031694
------------------------------------------------------------------------------


gen emp10 = 10*emp
jwdid emp10, ivar(countyreal) tvar(year) gvar(first_treat) method(ppmlhdfe)
estat simple                    // 10x the emp level result (scale-dependent)
estat simple, scale(index)      // same as emp result (scale-invariant)



estat simple                    // 10x the emp level result (scale-dependent)
------------------------------------------------------------------------------
             |            Delta-method
             | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      simple |  -289.1383    205.633    -1.41   0.160    -692.1716     113.895
------------------------------------------------------------------------------

estat simple, scale(index)      // same as emp result (scale-invariant)
------------------------------------------------------------------------------
             |            Delta-method
             | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      simple |  -.0286905   .0130213    -2.20   0.028    -.0542117   -.0031694
------------------------------------------------------------------------------

[Jeff Wooldridge]

This is great stuff, George! I used the new .ado and it replicated my "by hand" calculations exactly. This is a big help to empirical researchers. Has it been incorporated at ssc yet?

[George Ford]

Fernando is going to post soon.

And confirmed to work with all estat options (group, calendar, etc...).

[Jeff Wooldridge]

Excellent! I'm hoping that soon a "method(fraclogit)" option, and even "method(fracprobit)," will be allowed. This is actually a straightforward addition because everything would be just like logit/probit but no data check. Either fracreg or glm can be used. But I would mess it up, so I'll leave it to others.

[Frank Huang]

Dear George Ford,

Thank you very much for the new ado file.

After rechecking my results, I found that the estimates from jwdid (index) are very consistent with the main results in my study.

Just to confirm: the estimates from jwdid (index) can still be interpreted as proportional / percentage changes in the outcome, similar to standard Poisson DID coefficients. Is that correct?

Thank you again for your help.

[Jeff Wooldridge]

Frank: Your interpretation is correct. The quantities reported when you use estat event, scale(index) are weighted averages of the coefficients from the unrestricted Poisson DiD estimation.