Hello everyone,

I'm fitting a negative binomial regression model using nbreg in Stata 16.0 and then using margins to get a predicted probability of obtaining a count greater than 0 for each level of my categorical covariate. However, I am finding some of my 95% CIs contain impossible probabilities (<0 or >1). I’ve replicated an example below where I get a negative probability when prog = 1 and prog = 3.

I read that this is because the standard errors of these probability estimates are normally approximated using the delta method. I found a great package from Jeff Pitblado, transform_margins, which takes the linear prediction from the model and applies a transformation to the table produced in order to get 95% CIs within [0,1].

I have applied this to get the estimated number of events for each level of prog, but I was unsure how to use transform_margins to get the estimated probability and 95% CI of a given number of events after nbreg. Is there a way to apply it, or another method I could use?

Any advice is much appreciated, and I apologise in advance if I have made any errors in my post (this is the first one).

Code:
* Example generated by -dataex-. To install: ssc install dataex
clear
input float(daysabs prog math)
 4 2 63
 4 2 27
 2 2 20
 3 2 16
 3 2  2
13 2 71
11 2 63
 7 2  3
10 2 51
 9 3 49
 4 3 31
 5 2 22
 5 2 73
 6 1 77
 1 3 10
 0 2 89
 1 2 34
 0 2 35
 5 1 77
24 2  4
 2 2 21
 0 2 61
 1 2 41
 0 2 63
 8 2 19
 6 1 55
 7 2  6
 0 2 21
 2 2 70
 3 1 79
 0 2  1
 1 2 29
 3 2 40
 0 2 39
 0 2 68
28 2 20
 8 2 10
 5 2 84
 5 1  1
27 2 35
 5 2 32
18 2 34
19 2 40
 9 2 29
 9 2 20
 4 2 27
 2 2 74
 3 2 38
 9 2 21
20 2 81
end

// Fitting model and getting predicted probability of daysabs >=20
nbreg daysabs math i.prog
margins prog, predict(pr(20,.))

Predictive margins                              Number of obs     =         50
Model VCE    : OIM

Expression   : Pr(daysabs>=20), predict(pr(20,.))

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        prog |
          1  |   .0333275   .0539323     0.62   0.537    -.0723778    .1390328
          2  |   .0627511   .0291091     2.16   0.031     .0056982    .1198039
          3  |   .0220017   .0490054     0.45   0.653    -.0740471    .1180506
------------------------------------------------------------------------------

// Transforming margins to get estimated count of daysabs
margins prog, predict(xb)

Predictive margins                              Number of obs     =         50
Model VCE    : OIM

Expression   : Linear prediction, predict(xb)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        prog |
          1  |   1.667932   .5106264     3.27   0.001     .6671225    2.668741
          2  |   1.893794   .1681998    11.26   0.000     1.564129     2.22346
          3  |   1.539881   .6433472     2.39   0.017     .2789438    2.800819
------------------------------------------------------------------------------

transform_margins exp(@)

----------------------------------------------
             |         b         ll         ul
-------------+--------------------------------
        prog |
          1  |  5.301193   1.948622    14.4218
          2  |  6.644533    4.77851   9.239242
          3  |  4.664036   1.321733   16.45811
----------------------------------------------