How to predict effect after poisson regression

[Karl Yesler]

Hello Statalist,

I am trying to manually predict the effect of a 1 unit increase in a predictor variable after a poisson regression.

To illustrate what I mean, I've done it (correctly) for a probit regression as follows:

use http://www.ats.ucla.edu/stat/stata/dae/binary.dta, clear

probit admit gre gpa i.rank

margins, dydx(gre)  ///gives same results as manual dydx prediction below

predict xb, xb

gen yhat = (normal(xb))

replace gre = gre +1

predict xb_2, xb

gen yhat_2 = (normal(xb_2))

gen dydx = yhat_2 - yhat

collapse dydx

///gives same results as margins dydx prediction above
list dydx

For the poisson regression, I attempt to repeat as follows:

I assume there is something wrong with my understanding of the distribution function for poisson? If so, what is the appropriate code?

use http://www.ats.ucla.edu/stat/stata/dae/poisson_sim, clear

poisson num_awards math i.prog 

margins, dydx(math) 

predict xb , xb 

gen yhat = (lnfactorial(xb))

replace math = math + 1

predict xb_2 , xb 

gen yhat_2 = (lnfactorial(xb_2))

gen dydx = yhat_2 - yhat

collapse dydx

list dydx

Thanks!

[William Lisowski]

"lnfactorial". You keep using that function. I do not think it means what you think it means. (Apologies if you're not acquainted with the movie "The Princess Bride".)

I think you need to read up on the use of predict, for instances from the output of help predict. Instead of predicting xb and transforming it into yhat, predict can be used (without the eb option) to predict the "natural" yhat for the model that was fit. And I'm reasonably certain the transformation of eb into yhat does not involve computing the log of the factorial of eb.

[Karl Yesler]

William,

Thanks for your suggestion but I do not think the "natural" yhat is quite the same as what I am trying to get with the margins dydx command. Compare the results from margins dydx with "natural" predict below:

 use http://www.ats.ucla.edu/stat/stata/dae/poisson_sim, clear
(highschool and beyond (200 cases))

. 
. poisson num_awards math i.prog 

Iteration 0:   log likelihood = -182.75759  
Iteration 1:   log likelihood = -182.75225  
Iteration 2:   log likelihood = -182.75225  

Poisson regression                                Number of obs   =        200
                                                  LR chi2(3)      =      98.22
                                                  Prob > chi2     =     0.0000
Log likelihood = -182.75225                       Pseudo R2       =     0.2118

------------------------------------------------------------------------------
  num_awards |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        math |   .0701524   .0105992     6.62   0.000     .0493783    .0909265
             |
        prog |
   academic  |   1.083859    .358253     3.03   0.002     .3816962    1.786022
   vocation  |   .3698092   .4410703     0.84   0.402    -.4946727    1.234291
             |
       _cons |  -5.247124   .6584531    -7.97   0.000    -6.537669    -3.95658
------------------------------------------------------------------------------

. 
. margins, dydx(math) 

Average marginal effects                          Number of obs   =        200
Model VCE    : OIM

Expression   : Predicted number of events, predict()
dy/dx w.r.t. : math

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        math |    .044196   .0077519     5.70   0.000     .0290027    .0593894
------------------------------------------------------------------------------

. 
. predict yhat 
(option n assumed; predicted number of events)

. 
. replace math = math + 1
(200 real changes made)

. 
. predict yhat_2
(option n assumed; predicted number of events)

. 
. gen dydx = yhat_2 - yhat

. collapse dydx

. 
. list dydx

     +----------+
     |     dydx |
     |----------|
  1. | .0457831 |
     +----------+

[Maarten Buis]

That is because with the way you used predict you got the effect of a discrete change in your x (math) rather than a marginal effect.

[William Lisowski]

Try the code below for two runs of your probit analysis. The first predicts yhat by transforming the output of predict, xb as you did in #1 above. The second predicts yhat directly with the output of predict. The results are identical. I suggest, without going to the effort to confirm it, that if you used the correct formula to transform the prediction of xb into yhat for the poisson model (rather than lnfactorial(xb), which yielded 200 missing values for the 200 observations) those results would also be identical to the direct prediction of yhat you showed in #3 above.

We can conclude that the reason for the difference between your two approaches for the poisson model lies somewhere other than in the technique used for predicting yhat.

In particular, it lies in your understanding of margins, dydx(math). That computes the derivative of the prediction relative to the specified variable. Your calculation of (yhat(math+1) - yhat(math))/((math+1)-math) is at best an approximation to the derivative of the prediction (or conversely, the derivative of the prediction is at best an approximation to the results of a one-unit change). Apparently for your probit model that was a good approximation; for your poisson model, not quiet as good, although not bad, as your results in #3 above show.

use http://www.ats.ucla.edu/stat/stata/dae/binary.dta, clear
quietly probit admit gre gpa i.rank
margins, dydx(gre)
predict xb, xb
gen yhat = (normal(xb))
replace gre = gre +1
predict xb_2, xb
gen yhat_2 = (normal(xb_2))
gen dydx = yhat_2 - yhat
collapse dydx
list dydx

use http://www.ats.ucla.edu/stat/stata/dae/binary.dta, clear
quietly probit admit gre gpa i.rank
margins, dydx(gre)
predict yhat
replace gre = gre +1
predict yhat_2
gen dydx = yhat_2 - yhat
collapse dydx
list dydx

[Karl Yesler]

William and Martin,

Thank you for your explanations. That makes sense to me now.