Negative binomial model

[Mayura Shinde]

Hi,
I am running a negative binomial regression model for count data (no. of COPD exacerbations) among patients who have exacerbation following COPD diagnosis. I used both nbreg and glm model with nbinomial family and log link to fit the model. Both models provide similar coefficient estimate and incident rate ratios can also be estimated but I wanted to know if we could get incident rates for a given independent variable. Would we have to estimate these values after running the model. Any help would be greatly appreciated.

For other glm negative binomial models used for modeling inp hospitalizations or er visits, I have used the margins var1 command to estimate the predicted number of inp hosp or predicted no. of er visits. However, I wanted to know if incident rates could be estimated for a negative binomial regression model.

Thanks,
Mayura

[Clyde Schechter]

After -nbreg-, -margins- with the -ir- option will calculate adjusted incidence rates.

[Mayura Shinde]

Hi Clyde I tried using margins var1, ir after running the nbreg model but it gace an error saying that the ir option is not allowed

[Jeff Pitblado (StataCorp)]

I think Clyde meant that predict after nbreg has an ir option for predicting incidence rates, and you can compute marginal incidence rates and their marginal effects with margins.

Here is a short example taken from the nbreg manual entry.

First, fit the model:

. webuse rod93

. nbreg deaths i.cohort, exposure(exposure)

Fitting Poisson model:

Iteration 0:   log likelihood = -2160.0542  
Iteration 1:   log likelihood =  -2159.516  
Iteration 2:   log likelihood = -2159.5158  
Iteration 3:   log likelihood = -2159.5158  

Fitting constant-only model:

Iteration 0:   log likelihood = -187.06699  
Iteration 1:   log likelihood = -151.29069  
Iteration 2:   log likelihood = -131.82867  
Iteration 3:   log likelihood = -131.58459  
Iteration 4:   log likelihood = -131.58186  
Iteration 5:   log likelihood = -131.58186  

Fitting full model:

Iteration 0:   log likelihood = -131.58186  
Iteration 1:   log likelihood = -131.38447  
Iteration 2:   log likelihood =  -131.3799  
Iteration 3:   log likelihood =  -131.3799  

Negative binomial regression                    Number of obs     =         21
                                                LR chi2(2)        =       0.40
Dispersion     = mean                           Prob > chi2       =     0.8171
Log likelihood =  -131.3799                     Pseudo R2         =     0.0015

------------------------------------------------------------------------------
      deaths |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      cohort |
  1960-1967  |  -.2676187   .7237203    -0.37   0.712    -1.686085    1.150847
  1968-1976  |  -.4573957   .7236651    -0.63   0.527    -1.875753    .9609618
             |
       _cons |  -2.086731   .5118559    -4.08   0.000     -3.08995   -1.083511
ln(exposure) |          1  (exposure)
-------------+----------------------------------------------------------------
    /lnalpha |   .5939963   .2583615                       .087617    1.100376
-------------+----------------------------------------------------------------
       alpha |   1.811212   .4679475                       1.09157    3.005294
------------------------------------------------------------------------------
LR test of alpha=0: chibar2(01) = 4056.27              Prob >= chibar2 = 0.000

Next, predict some incidence rates:

. predict irates, ir

Next, compute marginal incidence rates for each cohort:

. margins cohort, predict(ir)

Adjusted predictions                            Number of obs     =         21
Model VCE    : OIM

Expression   : Predicted incidence rate, predict(ir)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      cohort |
  1941-1949  |   .1240922   .0635173     1.95   0.051    -.0003995    .2485838
  1960-1967  |   .0949553   .0485835     1.95   0.051    -.0002667    .1901772
  1968-1976  |   .0785416   .0401794     1.95   0.051    -.0002085    .1572917
------------------------------------------------------------------------------

Finally, compare marginal incidence rates for each cohort with the base-level cohort:

. margins r.cohort, predict(ir)

Contrasts of adjusted predictions
Model VCE    : OIM

Expression   : Predicted incidence rate, predict(ir)

-------------------------------------------------------------
                          |         df        chi2     P>chi2
--------------------------+----------------------------------
                   cohort |
(1960-1967 vs 1941-1949)  |          1        0.13     0.7156
(1968-1976 vs 1941-1949)  |          1        0.37     0.5445
                   Joint  |          2        0.37     0.8309
-------------------------------------------------------------

---------------------------------------------------------------------------
                          |            Delta-method
                          |   Contrast   Std. Err.     [95% Conf. Interval]
--------------------------+------------------------------------------------
                   cohort |
(1960-1967 vs 1941-1949)  |  -.0291369    .079967     -.1858693    .1275955
(1968-1976 vs 1941-1949)  |  -.0455505   .0751583      -.192858    .1017569
---------------------------------------------------------------------------

[Mayura Shinde]

Thanks Jeff for sharing the example. I do understand the logic now but it still does not completely provide information that I need to present from the analyses. I need to provide the exponentiated values of incident rates and 95% CI for spirometry below. Could you please help how they could be estimated.
Here is the output from my analysis:

. nbreg episode_ct i.spirometry i.age_cat age i.plan_type_cat i.commercial i.regcol_cat i.index_cat pre_ECI_SCO i.pre_smoke_flg pr
> e_all_INP_cnt, offset(log_risk_years) robust nolog irr

Negative binomial regression Number of obs = 140987
Dispersion = mean Wald chi2(15) = 8149.98
Log pseudolikelihood = -141231.85 Prob > chi2 = 0.0000


Robust
episode_ct IRR Std. Err. z P>z [95% Conf. Interval]

1.spirometry 1.378698 .0176037 25.15 0.000 1.344623 1.413636

age_cat
2 1.084357 .0220537 3.98 0.000 1.041982 1.128454
3 .9312842 .0303648 -2.18 0.029 .8736322 .9927408

age .9943475 .0009077 -6.21 0.000 .9925701 .9961282

plan_type_cat
2 .966913 .0133982 -2.43 0.015 .9410064 .9935329
3 .9301332 .0241905 -2.78 0.005 .8839088 .9787748

1.commercial .763994 .0110013 -18.69 0.000 .7427332 .7858634

regcol_cat
Midwest 1.018662 .016341 1.15 0.249 .9871322 1.051198
South 1.178552 .0196236 9.87 0.000 1.140711 1.217648
West 1.111342 .0212351 5.52 0.000 1.070491 1.153751

index_cat
ER 3.668882 .084592 56.38 0.000 3.506775 3.838483
INP 4.256831 .0999632 61.68 0.000 4.065347 4.457334

pre_ECI_SCO 1.003722 .0027709 1.35 0.178 .9983055 1.009167
1.pre_smoke_flg 1.110758 .0151819 7.69 0.000 1.081397 1.140916
pre_all_INP_cnt .9720787 .0109344 -2.52 0.012 .9508821 .9937477
_cons .5200266 .0247271 -13.75 0.000 .4737522 .570821
log_risk_years 1 (offset)

/lnalpha .1402715 .0160082 .108896 .1716469

alpha 1.150586 .0184188 1.115046 1.187259



predict irates, ir
(1414 missing values generated)

. margins spirometry, predict(ir)

Predictive margins Number of obs = 140987
Model VCE : Robust

Expression : Predicted incidence rate, predict(ir)

------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
spirometry |
0 | .3434874 .0021252 161.63 0.000 .3393222 .3476527
1 | .4735654 .005397 87.75 0.000 .4629875 .4841433
------------------------------------------------------------------------------

[Clyde Schechter]

Sorry, yes. Jeff is correct about what I meant to say. I was referring to -margins, predict(ir)-.

I need to provide the exponentiated values of incident rates and 95% CI for spirometry below

Really? Why do you want that. It sounds beyond bizarre.

The incidence rates themselves are exponentiated values of the linear predictor derived from the negative binomial regression. The adjusted incidence rates for those with and without spirometry are shown in the very last table in the output you posted. These are bread-and-butter statistics for this kind of setting. Why are they not what you need?

By the way, to enhance readability, in the future please post Stata output by pasting it in between code delimiters, as Jeff did. (See FAQ #12 if you are not familiar with code delimiters.)

[Mayura Shinde]

Thank you Clyde for the explanation. I didn't know the predicted incident rates estimated in the last table were exponentiated values. Thanks again for your help!

[Hassen Ali]

Dear Clyde,Jeff,Thank you very much!! I have learned a lot from your daily posts and replies.
Respectfully,Hassen