Incidence rates of events occuring more than once

[Jens Moller]

Hello

I wish to calculate incidence rates with 95% confidence intervals for event occurring more than once; e.g. drug prescriptions or hospital admission within a period. In the following example I use the diet example data.

use https://www.stata-press.com/data/r18/diet.dta, clear

The fail variable contains count of CHD which occurs more than once within the time window.

* Example generated by -dataex-. For more info, type help dataex
clear
input int id byte(fail month)
  1  0  8
  2  0 12
  3  0 11
  4  0  9
  5  0  9
  6  5  3
  7  0 11
  8  0  5
  9  0  2
 10  0  7
 11  0 10
 12  0  7
 13  0  9
 14  0 12
 15  0  5
 16  3  5
 17 12  2
 18  0  2
 19 12  2
 20  0  2
end

The obvious go-to command is stptime

. stset month, id(id) failure(fail) 

Survival-time data settings

           ID variable: id
         Failure event: fail!=0 & fail<.
Observed time interval: (month[_n-1], month]
     Exit on or before: failure

--------------------------------------------------------------------------
        337  total observations
          0  exclusions
--------------------------------------------------------------------------
        337  observations remaining, representing
        337  subjects
         80  failures in single-failure-per-subject data
      2,100  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =        12

. 

. stptime

        Failure _d: fail
  Analysis time _t: month
       ID variable: id

Estimated person-time and incidence rate

    Cohort |  Person-time   Failures        Rate   [95% conf. interval]
-----------+-----------------------------------------------------------
     Total |         2100         80   .03809524   .0305987    .0474283

But the stset appears to handle the number of failures corresponding to the count of fail!=0 which lead to incidence rates with too few events.

. count if fail!=0
  80

I would rather calculate the incidence rate from the sum of fails

. tabstat fail month, statistics(sum)

   Stats |      fail     month
---------+--------------------
     Sum |       674      2100
------------------------------

. di 674/2100
.32095238

Last time I calculated the incidence rates and 95% confidence intervals manually with the formulas but I guess Stata must have a standardized solution? Preferably a command allowing for stratification and storing the estimates as scalars or matrices. If it makes any difference, I use Stata18.

[Jens Moller]

A simple work-around is simply to recode the _d variable created with the stset command.

replace _d=fail if _d==1

The following stptime command delivers incidence rate and 95% confidence intervals (almost) corresponding to manual calculations

. stptime

        Failure _d: fail
  Analysis time _t: month
       ID variable: id

Estimated person-time and incidence rate

    Cohort |  Person-time   Failures        Rate   [95% conf. interval]
-----------+-----------------------------------------------------------
     Total |         2100        674   .32095238   .2976141    .3461208

. di 674/2100
.32095238

. di `r(rate)'-1.96*sqrt(`r(rate)'/2100)
.29672164

. di `r(rate)'+1.96*sqrt(`r(rate)'/2100)
.34518312

I am not perfectly happy about this solution so I am still interested in alternative solutions.

[Carlo Lazzaro]

Jens:
have you considered -poisson, irr-?

[Jens Moller]

Thank you for your suggestion, Carlo. I believe you mean something like:

. poisson fail, exposure(month) irr

Iteration 0:  Log likelihood = -1482.3599  
Iteration 1:  Log likelihood = -1482.3599  

Poisson regression                                      Number of obs =    337
                                                        LR chi2(0)    =   0.00
                                                        Prob > chi2   =      .
Log likelihood = -1482.3599                             Pseudo R2     = 0.0000

------------------------------------------------------------------------------
        fail |  Inc. rate   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   .3209524   .0123626   -29.50   0.000     .2976141    .3461208
   ln(month) |          1  (exposure)
------------------------------------------------------------------------------

The thought should have crossed my mind as it is an obvious approach. Thank you!

[Carlo Lazzaro]

Jens:
this is exactly what I figured out.
Some more predictors, if available, can spice up your analysis, I hope.

[Jens Moller]

Yes, the Poisson regression enables both the possibility of adjusting for confounding and simple stratification:

. poisson fail ibn.job , exposure(month) irr nocons

Iteration 0:  Log likelihood = -1476.4631  
Iteration 1:  Log likelihood = -1469.3917  
Iteration 2:  Log likelihood = -1469.3761  
Iteration 3:  Log likelihood = -1469.3761  

Poisson regression                                      Number of obs =    337
                                                        Wald chi2(3)  = 868.14
Log likelihood = -1469.3761                             Prob > chi2   = 0.0000

------------------------------------------------------------------------------
        fail |        IRR   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         job |
          0  |   .3031674   .0213838   -16.92   0.000      .264024    .3481141
          1  |   .4259928   .0277298   -13.11   0.000     .3749676    .4839614
          2  |   .2684032   .0174347   -20.25   0.000     .2363177     .304845
   ln(month) |          1  (exposure)
------------------------------------------------------------------------------

stptime allows also for stratification using the by option, but the stratified estimates are not stored. Thank you once again!