Plotting estimates obtained through predictnl with coefplot

[Maxence Morlet]

Hi all,

I have run a Poisson model using ppmlhdfe from ssc install, along the lines of

ppmlhdfe y c.x1 c.x1#c.x2, cl(unit_id) abs(fixed_effect_vectors) d(ppml1)

(x2 is not included by itself because it is perfectly collinear with the fixed effects).

This is a whole other question, but I followed Leitgoeb and McCabe et al. (2022) to obtain the interaction effect; x1 and x2 are both continuous variables. I am just interested here in how to plot estimated effects, using postestimation commands like predictnl.

In order to do this, after the estimation, I had to use predictnl:

predict fitted_value, xb
predictnl interaction_effect=((_b[x1]+_b[c.x1#c.x2]*x2)*(_b[c.x1#c.x2]*x1)+_b[c.x1#c.x2])*fitted_value, se(se_interaction_effect)

Now I would like to use coefplot from ssc install to plot the coefficient from x1 and its confidence interval, and the interaction effect and its confidence interval, generated through predictnl.

I know one normally runs est store coef_name and then plots it, however I do not know the procedure to follow when one would like to plot an estimate from predictnl.

A typical coefplot code might be, for this model, if one were to plot regression estimates:

coefplot model1, bylabel("Test"),  mlabel(cond(@pval<.05, "**", cond(@pval<.01, "***", cond(@pval<.05, "", "")))) drop(_cons) keep(x1 c.x1#c.x2) xline(0) byopts(compact cols(1))  coeflabels(x1 = "{bf:First Regressor}" c.x1#c.x2 = "{bf: Coef_on_interac}")

Please could someone help me find this code?

Below are toy data generated for illustration:

input float(id y x1 x2)
  1 2      .2488 -1.1430041
  2 1  1.3352333  1.1311661
  3 1  -.4587567 -1.4932548
  4 1   1.319296   7.923935
  5 2 -2.0332499   4.529856
  6 2   -.646085 -4.5292335
  7 1   5.043857   6.302709
  8 2   5.941895   6.645915
  9 0   3.718591  -.2793825
 10 3   5.071329   2.413525
 11 1   2.602065   7.383534
 12 4   .7288715   2.838717
 13 0    2.93099  -2.412951
 14 2   .4572504  -.7622573
 15 3 -.24726455  -5.066185
 16 2 -2.0624917   7.994617
 17 1  .04221758   3.689781
 18 3  2.0986538   3.779051
 19 3  .14305335 -1.4006286
 20 0  -.8681847   2.469832
 21 4  4.1309566   .6822363
 22 2  1.9485594  1.0957936
 23 1   7.379262   8.829028
 24 2     .07606  2.9235585
 25 5   4.741378  10.108137
 26 3    .304807   3.482954
 27 1  .11917886   4.861798
 28 3  2.0989184  1.8407828
 29 2 -1.9350915   4.587946
 30 3  2.0964124   7.992496
 31 1 -1.6835184   -.412823
 32 1 -.04724126   8.475094
 33 2   5.098127  -3.600388
 34 0  1.6679398   .8006483
 35 3 -.22119495 .025044275
 36 1  1.4548622   .4852332
 37 1  4.3716683   .8371546
 38 3   .6414328   4.800642
 39 5  -5.698874  2.7759914
 40 3   -.961277   4.419597

[Andrew Musau]

My impression is that in nonlinear models such as Poisson, especially when interactions are present, it’s better to focus on marginal effects rather than the coefficients themselves, since marginal effects are unambiguous. You can use the margins or nlcom commands, and the resulting estimates can be graphed if need be. Deriving the expression for nlcom isn't too difficult — for your model, it would be:

The model is:

$$\mathbb{E}[Y \mid X_1, X_2] = \exp(\eta), \quad \text{where } \eta = \beta_0 + \beta_1 X_1 + \beta_2 (X_1 \cdot X_2)$$

The marginal effect of \( X_1 \), is defined as:

$$\frac{\partial \mathbb{E}[Y]}{\partial X_1}= \frac{\partial}{\partial X_1} \exp(\beta_0 + \beta_1 X_1 + \beta_2 X_1 X_2)$$

Applying the chain rule:

$$= \exp(\beta_0 + \beta_1 X_1 + \beta_2 X_1 X_2) \cdot \left( \frac{\partial}{\partial X_1} (\beta_0 + \beta_1 X_1 + \beta_2 X_1 X_2) \right)$$


$$= \exp(\eta) \cdot (\beta_1 + \beta_2 X_2)$$

Therefore, the marginal effect of \( X_1 \) is:

$$\frac{\partial \mathbb{E}[Y]}{\partial X_1} = \mathbb{E}[Y \mid X] \cdot (\beta_1 + \beta_2 X_2)$$

clear
input float(id y x1 x2)
  1 2      .2488 -1.1430041
  2 1  1.3352333  1.1311661
  3 1  -.4587567 -1.4932548
  4 1   1.319296   7.923935
  5 2 -2.0332499   4.529856
  6 2   -.646085 -4.5292335
  7 1   5.043857   6.302709
  8 2   5.941895   6.645915
  9 0   3.718591  -.2793825
 10 3   5.071329   2.413525
 11 1   2.602065   7.383534
 12 4   .7288715   2.838717
 13 0    2.93099  -2.412951
 14 2   .4572504  -.7622573
 15 3 -.24726455  -5.066185
 16 2 -2.0624917   7.994617
 17 1  .04221758   3.689781
 18 3  2.0986538   3.779051
 19 3  .14305335 -1.4006286
 20 0  -.8681847   2.469832
 21 4  4.1309566   .6822363
 22 2  1.9485594  1.0957936
 23 1   7.379262   8.829028
 24 2     .07606  2.9235585
 25 5   4.741378  10.108137
 26 3    .304807   3.482954
 27 1  .11917886   4.861798
 28 3  2.0989184  1.8407828
 29 2 -1.9350915   4.587946
 30 3  2.0964124   7.992496
 31 1 -1.6835184   -.412823
 32 1 -.04724126   8.475094
 33 2   5.098127  -3.600388
 34 0  1.6679398   .8006483
 35 3 -.22119495 .025044275
 36 1  1.4548622   .4852332
 37 1  4.3716683   .8371546
 38 3   .6414328   4.800642
 39 5  -5.698874  2.7759914
 40 3   -.961277   4.419597
end

ppmlhdfe y c.x1 c.x1#c.x2, cl(id) noabsorb d(ppml1)
foreach var in x1 x2{
    qui sum `var'
    local m`var'=r(mean) 
}

*USING MARGINS
margins, dydx(x1) atmeans

*USING NLCOM 
nlcom exp(_b[_cons]+ (_b[x1]*`mx1') +(_b[c.x1#c.x2]*`mx1' *`mx2'))*(_b[x1]+(_b[c.x1#c.x2]*`mx2'))

Res.:

. *USING MARGINS
. margins, dydx(x1) atmeans

Conditional marginal effects                                Number of obs = 40
Model VCE: Robust

Expression: Predicted mean of y, predict()
dy/dx wrt:  x1
At: x1 = 1.273711 (mean)
    x2 = 2.625667 (mean)

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
          x1 |  -.1057673   .0964378    -1.10   0.273    -.2947819    .0832473
------------------------------------------------------------------------------

. 
. *USING NLCOM 
. nlcom exp(_b[_cons]+ (_b[x1]*`mx1') +(_b[c.x1#c.x2]*`mx1' *`mx2'))*(_b[x1]+(_b[c.x1#c.x2]*`mx2'))

       _nl_1: exp(_b[_cons]+ (_b[x1]*1.273710907436907) +(_b[c.x1#c.x2]*1.273710907436907 *2.625666642514989))*(_b[x1]+(_b[c.x1#c.x2]*2.6
> 25666642514989))

------------------------------------------------------------------------------
           y | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |  -.1057673   .0964378    -1.10   0.273    -.2947819    .0832473
------------------------------------------------------------------------------

To your question:

I am just interested here in how to plot estimated effects, using postestimation commands like predictnl.

Is it possible? Yes. But as you are predicting over observations, so you have a prediction and confidence interval for each observation. Do you want to graph CIs corresponding to all individuals in the dataset or you have a way to aggregate these?