Interpret margins for menbreg with an exposure variable.

[Brad Anderson]

I asked a somewhat similar question a couple of days ago and apologize for asking again. I just don't understand what the margins command is returning.

risk2 = total # of injections (persons can report multiple injections / day)
blrisk = total # of injections reported in the 30 days prior to baseline (can be > 30)
age, gender, wrace latinx and homenite are planned covariates
month = month of assessment (1, 3, 6, 9, 12)
dayexp = days persons were at risk
cond = intervention arm
sid = subject identifier

.
local covars ib12.month blrisk age i.gender i.wrace i.latinx homenite

. menbreg risk2 `covars'  i.cond if month>0 ///
>         , exp(dayexp) || sid:, vce(robust) irr nolog

Mixed-effects nbinomial regression              Number of obs     =        668
Overdispersion:            mean
Group variable:             sid                 Number of groups  =        213

                                                Obs per group:
                                                              min =          1
                                                              avg =        3.1
                                                              max =          5

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(11)     =      82.88
Log pseudolikelihood = -2677.3591               Prob > chi2       =     0.0000
                                  (Std. Err. adjusted for 213 clusters in sid)
------------------------------------------------------------------------------
             |               Robust
       risk2 |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       month |
       Mo 1  |   2.237339   .6230537     2.89   0.004     1.296252    3.861661
       Mo 3  |    2.88821   .7463729     4.10   0.000     1.740447    4.792882
       Mo 6  |   1.560984   .3566892     1.95   0.051      .997462     2.44287
       Mo 9  |   1.461995   .3357811     1.65   0.098     .9320677    2.293211
             |
      blrisk |   1.001802   .0005316     3.39   0.001     1.000761    1.002845
         age |   .9261572    .013625    -5.21   0.000      .899834    .9532503
             |
      gender |
       Male  |   3.264972   1.067272     3.62   0.000     1.720418    6.196192
             |
       wrace |
      White  |    1.34672   .4738177     0.85   0.398     .6757694    2.683838
             |
      latinx |
     Latinx  |   .9710289   .4488155    -0.06   0.949     .3924656    2.402496
    homenite |   .9997214   .0001893    -1.47   0.141     .9993505    1.000092
             |
        cond |
     Active  |   .3426471   .0971267    -3.78   0.000     .1965919    .5972121
       _cons |   1.354801   .8692848     0.47   0.636     .3852271    4.764687
  ln(dayexp) |          1  (exposure)
-------------+----------------------------------------------------------------
    /lnalpha |   1.153256   .1081295                      .9413257    1.365185
-------------+----------------------------------------------------------------
sid          |
   var(_cons)|   8.682843   1.156382                      6.688039    11.27262
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.
Note: _cons estimates baseline incidence rate (conditional on zero random effects).

. margins cond

Predictive margins                              Number of obs     =        668
Model VCE    : Robust

Expression   : Marginal predicted mean, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        cond |
    Control  |    8673.77   7431.901     1.17   0.243    -5892.488    23240.03
     Active  |   2972.042   2491.355     1.19   0.233    -1910.923    7855.007
------------------------------------------------------------------------------

. margins cond, exp(predict(mu)/dayexp)

Predictive margins                              Number of obs     =        668
Model VCE    : Robust

Expression   : predict(mu)/dayexp

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        cond |
    Control  |   91.84315   73.67205     1.25   0.213    -52.55142    236.2377
     Active  |   31.46979   24.68277     1.27   0.202    -16.90755    79.84713
------------------------------------------------------------------------------

. margins cond, exp((predict(mu)/dayexp)*30)

Predictive margins                              Number of obs     =        668
Model VCE    : Robust

Expression   : (predict(mu)/dayexp)*30

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        cond |
    Control  |   2755.294   2210.162     1.25   0.213    -1576.543    7087.132
     Active  |   944.0937   740.4831     1.27   0.202    -507.2264    2395.414
------------------------------------------------------------------------------

The crude 30-day rate during follow-up is about 40 events / 30 days. "I would like to be able to say something like: The estimated 30-day marginal mean rate is about XX in the active intervention arm and about YY in the control arm." I don't understand what the marginal estimates reported above mean with respect to expected events / some unit of time, and would like to express the outcome as estimated events / 30 days. Thanks.

[Clyde Schechter]

The output you see in -margins- is showing you the result of predicting the number of events for each observation in your data set, with each observation treated as a Control, dividing that by days of exposure, and multiplying by 30. Then do the same thing with each observation treated as Active. Those are the numbers you are seeing.

They are probably not what you want. In particular, they are not event rates for the Control and Active groups. They are average individual rates for those groups. If you have even a small number of entities where exposure is brief, the individual rate for that person will be estimated in a very unstable way, and some of them may turn out to be huge numbers which then drag up the mean.

What I suspect you want is a different order of operations. You want to have -margins- give you the average number of predicted events in each group, and then divide that by the average days of exposure in each group, and multiply by 30. -margins- can't do that directly. You have to give it a little help. In the code below, I assume that cond is coded with Active = 1 and Control = 0

forvalues c = 0/1 {
    summ expdays if cond == `c', meanonly
    local mean_exposure_`c' `r(mean)'
    margins, at(cond = `c') exp(30*predict(mu)/`mean_exposure_`c'')

These results, because they are adjusted for your covariates, will not exactly replicate the crude rates in each group, but they probably will be in the same ballpark. (No guarantees: a small number of extreme outliers can wreak havoc here as well.)

[Brad Anderson]

Thanks so much Clyde! I would've never figured that out on my own. I knew they wouldn't replicate the crude rates, but I really had no idea what they were estimating.