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
version 17
poisson price i.foreign c.mpg, vce(robust) nolog
quietly margins, at((p10) mpg) at((p90) mpg) at((p50) mpg) post
return scalar p10_over_p50 = _b[1._at]/_b[3._at]
return scalar p90_over_p50 = _b[2._at]/_b[3._at]
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.
