Estimating risk ratios in Poisson model in Stata

Question

This is a coding question. Basically, in a Poisson model, I would like to calculate risk ratio between 10th percentile and 50th percentile (i.e., predicted value at 10th percentile divided by predicted value at 50th percentile) and 90th percentile over 50th percentile.

The current codes that I know are:

poisson wt_num s1 c.s1#c.s1 c.s1#c.s1#c.s1, vce(robust) irr
_pctile s1, p(10 50 90)
margins, at(s1 = (`r(r1)' `r(r2)' `r(r3)')) contrast(atcontrast(rb2)) 

The last 2 line of codes _pctile s1, p(10 50 90) and margins, at(s1 = ('r(r1)' 'r(r2)' 'r(r3)')) contrast(atcontrast(rb2)) allow me to get predicted probability at these 3 values, but using contrast function is incorrect as they take predicted prob(10th %ile) - predicted prob(50th %ile) instead of dividing to get ratios. Any suggestion would be appreciated!

Answer 1

The next step is to take the predicted margins, which are themselves random since they are means, and calculate their ratios. Since a ratio of two random variables is a non-linear transformation, to get the standard error of that ratio, you need to either apply the delta method or else bootstrap.

Personally, I would report both or favor the bootstrap version unless I had a huge sample and relatively few parameters.

Here's an example of both:

. sysuse auto, clear
(1978 automobile data)

. poisson price i.foreign c.mpg, vce(robust) nolog

Poisson regression                                      Number of obs =     74
                                                        Wald chi2(2)  =  33.91
                                                        Prob > chi2   = 0.0000
Log pseudolikelihood = -28478.503                       Pseudo R2     = 0.3526

------------------------------------------------------------------------------
             |               Robust
       price | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
     foreign |
    Foreign  |   .2849739   .0876098     3.25   0.001     .1132618     .456686
         mpg |  -.0524904   .0094258    -5.57   0.000    -.0709647   -.0340162
       _cons |   9.723688   .1967522    49.42   0.000     9.338061    10.10932
------------------------------------------------------------------------------

. margins, at((p10) mpg) at((p90) mpg) at((p50) mpg) post // coefl

Predictive margins                                          Number of obs = 74
Model VCE: Robust

Expression: Predicted number of events, predict()
1._at: mpg = 14 (p10)
2._at: mpg = 29 (p90)
3._at: mpg = 20 (p50)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |   8798.503   688.4214    12.78   0.000     7449.221    10147.78
          2  |   4003.724   350.1841    11.43   0.000     3317.376    4690.072
          3  |   6421.418   282.5083    22.73   0.000     5867.712    6975.124
------------------------------------------------------------------------------

. 
. /* Delta Method SEs */
. nlcom ///
> (p10_over_p50:_b[1._at]/_b[3._at]) ///
> (p90_over_p50:_b[2._at]/_b[3._at])

p10_over_p50: _b[1._at]/_b[3._at]
p90_over_p50: _b[2._at]/_b[3._at]

------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
p10_over_p50 |   1.370181   .0774903    17.68   0.000     1.218302    1.522059
p90_over_p50 |   .6234954   .0528925    11.79   0.000      .519828    .7271628
------------------------------------------------------------------------------

. 
. /* Bootstrap SEs */
. capture program drop myratios

. program myratios, rclass
  1.         version 17
  2.         poisson price i.foreign c.mpg, vce(robust) nolog
  3.         quietly margins, at((p10) mpg) at((p90) mpg) at((p50) mpg) post
  4.         return scalar p10_over_p50 = _b[1._at]/_b[3._at]
  5.         return scalar p90_over_p50 = _b[2._at]/_b[3._at]
  6. end

. bs ///
> p10_over_p50 = r(p10_over_p50) ///
> p90_over_p50 = r(p90_over_p50) ///
> , seed(18092023) nodots reps(1000): myratios

Bootstrap results                                        Number of obs =    74
                                                         Replications  = 1,000

      Command: myratios
 p10_over_p50: r(p10_over_p50)
 p90_over_p50: r(p90_over_p50)

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
p10_over_p50 |   1.370181   .1111034    12.33   0.000     1.152422    1.587939
p90_over_p50 |   .6234954   .0704366     8.85   0.000     .4854422    .7615485
------------------------------------------------------------------------------

As expected, the bootstrapped SEs are a bit wider here.

To explain the code a bit, the post option in margins replaces the Poisson coefficients with the average predictions and their VCE as estimation results. This allows you to use nlcom on them to get the ratios of the averages and their SEs.

This gets you $$\frac{\frac{1}{N}\sum_i^N \mathbf E[price \vert mpg = X, foreign_i]}{\frac{1}{N}\sum_i^N \mathbf E[price \vert mpg = 20, foreign_i]}$$

where X is either the 90th or the 10th percentile of mpg. That can be simplified even further to

$$\frac{1}{N}\sum_i^N \frac{\mathbf E[price \vert mpg = X, foreign_i]}{\mathbf E[price \vert mpg = 20, foreign_i]} =\exp(\beta_{mpg}\cdot(X-20))$$

Your example is more complicated with all the interactions, so no such simplification may be possible.