Panel Data with Count models (negative binomial and poisson) and non linear difference-in-difference

[Spirae Raïs]

Hi, I'm new here (and new to stata/statistical analysis) an I've read some previous posts but I am still confused and have some questions as I am not familiar at all with panel data analysis and nonlinear models!

I have some panel data from 1996 to 2015 and try to evaluate the impact of a particular policy occuring in 2000. My outcome variable is the number of clinical trials for a particular disease in a particular year, and I want to know if this number was impacted by the introduction of the policy. As my dataset contains many many 0 (no clinical trials, sometimes from 1996 to 2015), I understood that maybe a neg binomial or poisson panel model were more appropriate in my case, to account for 0 and overdispersion (Am I right?).
So I tried to estimate : Y_i = α + β₁Treated_it + β₂Policy_it + β₃Treated_it * Policy_it + ε
Treated is a dummy==1 if the disease in the sample is targeted by the policy and Policy is a dummy==1 after 2000.

Here are my questions:
I run the following regressions:

xtpoisson nb_trial_y treated##policy i.date, fe vce(robust)

xtnbreg nb_trial_y treated##policy i.date, fe

Does that make any sense? Because in the output the coeff for my "treated" variable is omitted and I do not understand why.
I am not sure whether to use disease fixed effects (would this simply be a dummy for either?) or year fixed effects (that would simply be a dummy for each year as I did with i.date?), and how it is redundant.
I also wanted to know if the use of "margins" was correct to conclude for the impact of the policy.

Thank you very much for your time and help!
Spirae

[Clyde Schechter]

Because in the output the coeff for my "treated" variable is omitted and I do not understand why.

This is to be expected. The variable treated is constant over time within each disease, so it is colinear with the fixed effects. To identify the model, Stata omits it. This doesn't matter anyway since the parameter of interest here is the coefficient of the Treat#Policy interaction. Even if the effect of Treat were estimable in the model it isn't really of interest.

If you look carefully at your outputs, you will also see that in addition to the usual omission of one year indicator, there is a second year that has been omitted: this is because the year indicators are colinear with the Policy variable. In this case, what Stata drops is one of the year indicators, but the principle is the same. You cannot have in the same model a constant term, all of the panel fixed effects and a variable that represents a time-invariant attribute of the panels. You cannot have in the same model a constant term, all of the time fixed effects, and a variable that is constant across all panels in each time period. Something has to get dropped to identify the model. From the perspective of estimating the policy effect, however, it doesn't matter which ones get dropped, as the interaction term, which is the key, is not colinear with either the panel or time effects, and its value is the same regardless of which variables get omitted to fix the colinearity problem.

I am not sure whether to use disease fixed effects (would this simply be a dummy for either?) or year fixed effects (that would simply be a dummy for each year as I did with i.date?), and how it is redundant.

Well, I was under the impression that when you -xtset- the data you specified disease as the panel. In that case, you have already incorporated disease fixed effects in the model by specifying the -fe- option. Were you to add disease specific indicators (dummies) you will find that Stata simply drops them anyway due to colinearity (actually, in this case, identity) with the fixed effects. The use of i.date provides your time fixed effects. For an outcome variable such as clinical trials it seems to me that using both disease and date effects is appropriate here.

I also wanted to know if the use of "margins" was correct to conclude for the impact of the policy.

Yes, although there is some controversy about exactly what margins to estimate/compare in these non-linear models. See for example the recent discussion of this issue at https://www.statalist.org/forums/for...fference-model, beginning with post #7.

[Spirae Raïs]

Clyde Schechter Thank you very much for your answer and your clear explanations.

I have another question, but maybe I should give you more info about my dataset before. I have do have more diseases in my treatment group than in my control group. My outcome variable is the number of clinical trials by disease for a given year (nb_trial_y)

. xtsum nb_trial_y if treated==1

Variable         |      Mean   Std. Dev.       Min        Max |    Observations
-----------------+--------------------------------------------+----------------
nb_tri~y overall overall |   .033435    .471112          0         46 |     N =  195783
         between |             .2582776          0   9.619048 |     n =    9323
         within  |             .3940127  -9.585613   36.41439 |     T =      21

But in my treatment group, a large number of diseases do not have any clinical trials. Here is the summary for the diseases with at least one clinical trial during the period 1996-2015

 xtsum nb_trial_y if treated==1 & research_nul!=1

Variable         |      Mean   Std. Dev.       Min        Max |    Observations
-----------------+--------------------------------------------+----------------
nb_tri~y overall |  .3891564   1.563644          0         46 |     N =   16821
         between |             .7991886    .047619   9.619048 |     n =     801
         within  |             1.344261  -9.229891   36.77011 |     T =      21

And here is my control group

. xtsum nb_trial_y if treated==0

Variable         |      Mean   Std. Dev.       Min        Max |    Observations
-----------------+--------------------------------------------+----------------
nb_tri~y overall |  4.331633   8.732587          0         78 |     N =    1764
         between |             6.252933    .047619   26.09524 |     n =      84
         within  |             6.132085  -21.28741   56.71259 |     T =      21

Following your advices, I ran my two regressions, "treated" being my dummy variable for treatment group and "policy" a dummy variable ==1 after the introduction of the policy in 2000:

xtpoisson nb_trial_y treated##policy i.date if date<2016, fe vce(robust) 

xtnbreg nb_trial_y treated##policy i.date, fe

But when I do compute the margins using the following commands, it does not work at all:

margins treated##policy
margins treated, dydx(policy) pwcompare
margins policy, dydx(treated) pwcompare(effects)

Let me show you the results of one of the regression and the margins.

.

xtnbreg nb_trial_y treated##policy i.date , fe 
note: 8527 groups (170540 obs) dropped because of all zero outcomes
note: 2015.date omitted because of collinearity

Iteration 0:   log likelihood = -9036.9496  
Iteration 1:   log likelihood = -8897.1376  
Iteration 2:   log likelihood = -8868.5974  
Iteration 3:   log likelihood = -8865.4967  
Iteration 4:   log likelihood = -8865.4805  
Iteration 5:   log likelihood = -8865.4805  

Conditional FE negative binomial regression     Number of obs      =     17600
Group variable: disease                         Number of groups   =       880

                                                Obs per group: min =        20
                                                               avg =      20.0
                                                               max =        20

                                                Wald chi2(21)      =   2782.88
Log likelihood  = -8865.4805                    Prob > chi2        =    0.0000

--------------------------------------------------------------------------------
    nb_trial_y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
     1.treated |  -3.877878   .2585097   -15.00   0.000    -4.384548   -3.371209
      1.policy |   3.231949   .2698093    11.98   0.000     2.703132    3.760765
               |
treated#policy |
          1 1  |   1.041258   .2378041     4.38   0.000     .5751704    1.507345
               |
          date |
         1997  |   .4891134   .3353764     1.46   0.145    -.1682122    1.146439
         1998  |   .7974593   .3191896     2.50   0.012     .1718592    1.423059
         1999  |   .8576254   .3122902     2.75   0.006     .2455479    1.469703
         2000  |  -2.108035   .1392212   -15.14   0.000    -2.380903   -1.835166
         2001  |  -1.903504    .135507   -14.05   0.000    -2.169093   -1.637916
         2002  |  -1.864774   .1315328   -14.18   0.000    -2.122574   -1.606975
         2003  |  -1.646048   .1228163   -13.40   0.000    -1.886763   -1.405332
         2004  |  -.7919469   .0877606    -9.02   0.000    -.9639544   -.6199394
         2005  |  -.2644524   .0754069    -3.51   0.000    -.4122473   -.1166575
         2006  |  -.2323217   .0746455    -3.11   0.002    -.3786243   -.0860191
         2007  |   -.053566   .0714817    -0.75   0.454    -.1936676    .0865356
         2008  |   .0809973   .0698779     1.16   0.246    -.0559608    .2179555
         2009  |   .2939243   .0676548     4.34   0.000     .1613233    .4265252
         2010  |   .2833859   .0677378     4.18   0.000     .1506222    .4161496
         2011  |   .2884994     .06704     4.30   0.000     .1571034    .4198954
         2012  |   .2104043    .067859     3.10   0.002     .0774032    .3434055
         2013  |   .2102663   .0674276     3.12   0.002     .0781105     .342422
         2014  |   .1395325   .0681642     2.05   0.041     .0059331    .2731319
         2015  |          0  (omitted)
               |
         _cons |  -.9985218   .2760897    -3.62   0.000    -1.539648    -.457396
--------------------------------------------------------------------------------

. 
end of do-file

. margins policy#treated

Predictive margins                                Number of obs   =      17600
Model VCE    : OIM

Expression   : Linear prediction, predict()

--------------------------------------------------------------------------------
               |            Delta-method
               |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
policy#treated |
          0 0  |          .  (not estimable)
          0 1  |          .  (not estimable)
          1 0  |          .  (not estimable)
          1 1  |          .  (not estimable)
--------------------------------------------------------------------------------

Do you have any idea why Stata cannot estimate the margins? I have the same problem with my xtpoisson regression.

Thank you very much for your help!
Spirae

[Nick Cox]

Across cultures there are many different naming conventions, but I'd surprised if Spirae Spirae was your full real name. Please follow FAQ Advice and alert the administrators to your full real name if that guess is correct.

[Spirae Raïs]

Across cultures there are many different naming conventions, but I'd surprised if Spirae Spirae was your full real name. Please follow FAQ Advice and alert the administrators to your full real name if that guess is correct.

You are right. I just read the FAQ and sent an email to the admins about my real surname. Sorry for not using it from the beginning!

[Clyde Schechter]

An imbalance between the number of entities or observations in the treatment and control groups is not a problem.

Non-estimability in the -margins- command usually arises from empty cells. Do you have observations in the estimation sample for each combination of treated and policy? (-tab treated policy if e(sample)-; bear in mind that observations with constant zero outcome are omitted from the estimation sample.) If this is your problem, then one solution is to specify the -asbalanced emptycells(reweight)- options in the -margins- command. Another solution is to specify the -noestimcheck- option. I suggest you carefully read about each of these -margins- options so that you know what they do: each of them is identifying an unidentifiable parameter, but doing so in a different way, so you need to be clear which way is better for your purposes.

[Spirae Raïs]

Thank you very much for your answer.
I have observations in the estimation sample for each combination of treated and policy

tab policy treated if e(sample)

           |        policy
   treated |         0          1 |     Total
-----------+----------------------+----------
         0 |       336      1,344 |     1,680 
         1 |     3,184     12,736 |    15,920 
-----------+----------------------+----------
     Total |     3,520     14,080 |    17,600

When I specify "asbalanced emptycells(reweight)", it doesn't work.

. margins policy#treated,  asbalanced emptycells(reweight)

Adjusted predictions                              Number of obs   =      17600
Model VCE    : Robust

Expression   : Linear prediction, predict()
Empty cells  : reweight
at           : treated          (asbalanced)
               policy           (asbalanced)
               date             (asbalanced)

--------------------------------------------------------------------------------
               |            Delta-method
               |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
policy#treated |
          0 0  |          .  (not estimable)
          0 1  |          .  (not estimable)
          1 0  |          .  (not estimable)
          1 1  |          .  (not estimable)
--------------------------------------------------------------------------------

But when I specify the option "noestimcheck", I do get some estimates :

. margins treated#policy, noestimcheck

Predictive margins                                Number of obs   =      17600
Model VCE    : Robust

Expression   : Linear prediction, predict()

--------------------------------------------------------------------------------
               |            Delta-method
               |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
treated#policy |
          0 0  |  -.1927554   .0665191    -2.90   0.004    -.3231306   -.0623803
          0 1  |   3.325529   .2944883    11.29   0.000     2.748342    3.902715
          1 0  |  -.1927554   .0665191    -2.90   0.004    -.3231306   -.0623803
          1 1  |   5.046804   .3519192    14.34   0.000     4.357055    5.736553
--------------------------------------------------------------------------------

I do have the feeling that something's wrong but I couldn't tell what...

Thanks again for your time and help !
Spirae

[Clyde Schechter]

OK. I share your feeling that something is wrong, but I don't know in what direction to point you.

Bear in mind also that your -margins- command is giving you values of xb. For your purposes, you may want to specify a different outcome such as nu0 in the -predict()- option of -margins-.

[Spirae Raïs]

Thank you very much ! I changed the predict default option to nu0 because I have a non negative outcome.

If now I want to interpret my interaction term and conclude on the effect of the policy (see how much the policy affected the number of clinical trials in my treatment group as compared to my treatment group), I have to run:

. margins treated, dydx(policy) pwcompare predict(nu0) noestimcheck

Pairwise comparisons of average marginal effects
Model VCE    : Bootstrap

Expression   : Predicted number of events (assuming u_i=0), predict(nu0)
dy/dx w.r.t. : 1.policy

--------------------------------------------------------------
             |   Contrast Delta-method         Unadjusted
             |      dy/dx   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
1.policy     |
     treated |
     1 vs 0  |  -6.761898   2.403766     -11.47319   -2.050604
--------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the
      base level.

As I see the results, my interpretation would be : my model-predicted average difference-in-difference number of log counts is -6.76. So the policy increased the disease-specific number of clinical trials by 6.76 on average during the period 2000-2015 in my treatment group as compared to my control group? Am I right ? Is it increasing ou decreasing ? As it's written in the output treated==1 vs treated==0, I concluded the policy increased the number of trials. But maybe I am wrong.

Thank you very much again for your time and help!
Spirae

[Clyde Schechter]

It is sometimes difficult to interpret which direction 1 vs 0 means. To make it clearer, run -margins treated#policy, predict(nu0) noestimcheck- and you will be able to see directly whether the change from pre-to-post-policy is larger in the treated group or in the control group. (You actually show that in #7, but those numbers don't seem to be compatible with what you show in #9, so I'm concerned that something happened in between or that you are using a different data set or something.)

[Spirae Raïs]

You are right, I was only working on a subsample. If I take the same regression than in #7 I get large figures.

. margins treated#policy, noestimcheck

Predictive margins                                Number of obs   =      17600
Model VCE    : Robust

Expression   : Linear prediction, predict()

--------------------------------------------------------------------------------
               |            Delta-method
               |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
treated#policy |
          0 0  |  -.1927554   .0665191    -2.90   0.004    -.3231306   -.0623803
          0 1  |   3.325529   .2944883    11.29   0.000     2.748342    3.902715
          1 0  |  -.1927554   .0665191    -2.90   0.004    -.3231306   -.0623803
          1 1  |   5.046804   .3519192    14.34   0.000     4.357055    5.736553
--------------------------------------------------------------------------------

. margins treated#policy, predict(nu0) noestimcheck

Predictive margins                                Number of obs   =      17600
Model VCE    : Robust

Expression   : Predicted number of events (assuming u_i=0), predict(nu0)

--------------------------------------------------------------------------------
               |            Delta-method
               |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
treated#policy |
          0 0  |   1.200861   .1176051    10.21   0.000     .9703596    1.431363
          0 1  |   40.50086   13.69729     2.96   0.003     13.65466    67.34705
          1 0  |   1.200861   .1176051    10.21   0.000     .9703596    1.431363
          1 1  |   226.4668   89.46357     2.53   0.011     51.12141    401.8122
--------------------------------------------------------------------------------

. margins treated, dydx(policy) pwcompare predict(nu0) noestimcheck

Pairwise comparisons of average marginal effects
Model VCE    : Robust

Expression   : Predicted number of events (assuming u_i=0), predict(nu0)
dy/dx w.r.t. : 1.policy

--------------------------------------------------------------
             |   Contrast Delta-method         Unadjusted
             |      dy/dx   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
1.policy     |
     treated |
     1 vs 0  |   185.9659   81.09337      27.02584     344.906
--------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the
      base level.

So here I can reasonably say that the policy affected positively the number of clinical trials in my treatment group ? The average DID number of clinical trials is 186? Can I go much further than that ? Please correct me if I am totally wrong.

I also did the same analysis with my regression "xnbreg", and I saw that the predicted number of events was negative. How is it possible that poisson and neg binomial models contradicts each other? And how I am suppose to choose between the two models ? Lrtest?

Thanks a lot for your help and time !
Spirae

[Clyde Schechter]

Yes, according to this model, the average DID is 186. The confidence interval gives you a sense of how much uncertainty there is in that estimate. It's a very wide confidence interval. Still, even the lower limit, about 27, looks like a large effect.

As for choosing between models giving contradictory results, there are a couple of things to think about:

1. Poisson models are pretty robust when it comes to estimating means, even when the error distribution does not look very Poisson-like at all. The same cannot be said of negative binomial models. So, all else equal, I'd assign greater credibility to the effect estimates from a Poisson model (though not so much to the standard errors, confidence intervals, and p-values.)

2. You should check the fit of both models. Calculate predicted values and see how those line up with the observed values in the two models. If both models fit pretty well, then based on #1 I would go with Poisson. But if one of them clearly fits the data better than the other, that would override #1.

[Spirae Raïs]

Thank you very much for your answers! After all, Poisson model seems to be more appropriate as it fits the data better than Neg binomial! Thanks again !

Today I also wanted to explore two different questions :
- the common trend assumption
- try to see how my treatment effect changes with time (my guess is that the impact of the policy was not immediate but took several years, but I may be wrong).

I understood that it was possible to answer both questions using leads and lags ( SLIDE 25 on https://www2.warwick.ac.uk/fac/soc/e...re_3_-_did.pdf )

So I tried to apply this using the following commands:

xtpoisson nb_trial_y treated##i.date i.date, fe vce(robust)

Please correct me if I wronly interpreted this method.

note: 8529 groups (162051 obs) dropped because of all zero outcomes

Iteration 0:   log pseudolikelihood =  -11931.45  
Iteration 1:   log pseudolikelihood = -8952.2168  
Iteration 2:   log pseudolikelihood = -8439.5994  
Iteration 3:   log pseudolikelihood = -8382.7699  
Iteration 4:   log pseudolikelihood = -8379.0648  
Iteration 5:   log pseudolikelihood = -8379.0087  
Iteration 6:   log pseudolikelihood = -8379.0086  

Conditional fixed-effects Poisson regression    Number of obs      =     16074
Group variable: disease                         Number of groups   =       846

                                                Obs per group: min =        19
                                                               avg =      19.0
                                                               max =        19

                                                Wald chi2(36)      =   1738.59
Log pseudolikelihood  = -8379.0086              Prob > chi2        =    0.0000

                                (Std. Err. adjusted for clustering on disease)
------------------------------------------------------------------------------
             |               Robust
  nb_trial_y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   1.treated |          0  (omitted)
             |
        date |
       1998  |  -.2231436    .466369    -0.48   0.632     -1.13721    .6909228
       1999  |   .2623643    .524094     0.50   0.617    -.7648412     1.28957
       2000  |   .0953102   .4433158     0.21   0.830    -.7735729    .9641932
       2001  |   1.193922   .4478497     2.67   0.008     .3161532    2.071692
       2002  |   .9162907   .4009988     2.29   0.022     .1303476    1.702234
       2003  |   1.458615   .3279152     4.45   0.000      .815913    2.101317
       2004  |   2.054124   .3534185     5.81   0.000     1.361436    2.746811
       2005  |   2.151762   .3517091     6.12   0.000     1.462425    2.841099
       2006  |   2.292535   .3131281     7.32   0.000     1.678815    2.906254
       2007  |   2.617396   .3254541     8.04   0.000     1.979517    3.255274
       2008  |   2.701361   .3405391     7.93   0.000     2.033917    3.368805
       2009  |   2.766319   .3416073     8.10   0.000     2.096781    3.435857
       2010  |   2.617396   .3507009     7.46   0.000     1.930035    3.304757
       2011  |   2.772589   .3415269     8.12   0.000     2.103208    3.441969
       2012  |   2.791165   .3594391     7.77   0.000     2.086678    3.495653
       2013  |   2.809403   .3743131     7.51   0.000     2.075762    3.543043
       2014  |    2.97553   .3534238     8.42   0.000     2.282832    3.668227
       2015  |   2.772589   .3479853     7.97   0.000      2.09055    3.454627
             |
treated#date |
     1 1998  |   .0408264   .7643067     0.05   0.957    -1.457187     1.53884
     1 1999  |   .1431005   .7596997     0.19   0.851    -1.345883    1.632085
     1 2000  |   1.578666    .635776     2.48   0.013     .3325679    2.824764
     1 2001  |     .01005   .6761922     0.01   0.988    -1.315262    1.335362
     1 2002  |    .624154   .6335978     0.99   0.325    -.6176749    1.865983
     1 2003  |   -.072321   .5862076    -0.12   0.902    -1.221267    1.076625
     1 2004  |   1.235521   .5540551     2.23   0.026     .1495929    2.321449
     1 2005  |   2.035364   .5449883     3.73   0.000     .9672064    3.103521
     1 2006  |   1.847951   .5250434     3.52   0.000     .8188845    2.877017
     1 2007  |   1.728571   .5288043     3.27   0.001      .692134    2.765009
     1 2008  |   1.966783   .5438115     3.62   0.000     .9009325    3.032634
     1 2009  |   2.271715   .5395353     4.21   0.000     1.214245    3.329185
     1 2010  |   2.300149   .5464446     4.21   0.000     1.229137    3.371161
     1 2011  |    1.99243   .5393955     3.69   0.000     .9352341    3.049626
     1 2012  |   1.963861   .5504633     3.57   0.000     .8849723    3.042749
     1 2013  |   1.754945   .5591928     3.14   0.002     .6589475    2.850943
     1 2014  |   1.127113   .5510712     2.05   0.041     .0470338    2.207193
     1 2015  |    .969041   .5484857     1.77   0.077    -.1059713    2.044053

So if I correctly understand, my lag variables are not significant, meaning that my common parallel assumption before the policy adoption cannot be rejected ?
Then I seems that it took 4 years for the policy to significatively affect the number of clinical trials in my treatment group.

I tried to have a look at the margins for each year but stata is taking forever to give the results of:

margins treated#i.date, predict(nu0) noestimcheck

Last question, is there any way I can graph the change of my treatment effect with time (I show you the graph I am thinking about in picture)?

Thank you very much again for your great advices and help!
Spirae

[Clyde Schechter]

Neither the command you show, nor the output, has any mention of lags or leads in it. Your code redundantly specifies i.year twice, and is other wise the basic DID model with year variables. If you think that the treatment effect will only be seen after a four year delay, then the way to model that is with

xtpoisson nb_trial_y treated##L4.policy, fe

If you want to model a gradual onset of effect over 4 years, it would be

xtpoisson nb_trial_y treated##L(1/4).policy, fe

[Spirae Raïs]

Thanks I didn't know how to use time series operators!

So if I want to model this :

I would have to do :

xtpoisson nb_trial_y  treated##F1.policy  treated##policy i.date, fe vce(robust)

xtpoisson nb_trial_y  treated##F2.policy  treated##policy i.date, fe vce(robust)

xtpoisson nb_trial_y  treated##F2.policy  treated##policy treated##L.policy treated##L2.policy  treated##L3.policy i.date, fe vce(robust)

Is it correct?

Thank you very much for your time and help!
Spirae