Censored Poisson regression

Let's see it work

Below we study the number of car accidents a person has during a year. The number recorded is 0, 1, 2, or 3, and 3 means 3 or more accidents. The number is right-censored.

We will model the determinants of accidents as the number of previous accidents, whether the driver is a parent, and the number of traffic tickets the driver received during the previous year.

We type

. cpoisson accidents i.past i.parent i.ntickets, ul(3) irr

Initial:      Log likelihood =  -2657.162
Rescale:      Log likelihood =  -2657.162
Iteration 0:  Log likelihood =  -2657.162
Iteration 1:  Log likelihood = -2638.9113
Iteration 2:  Log likelihood = -2638.7142
Iteration 3:  Log likelihood = -2638.6901
Iteration 4:  Log likelihood = -2638.6863
Iteration 5:  Log likelihood = -2638.6859
Iteration 6:  Log likelihood = -2638.6858
Iteration 7:  Log likelihood = -2638.6858

Censored Poisson regression                        Number of obs     =   3,000
                                                          Uncensored =   2,840
Limits: Lower = 0                                      Left-censored =       0
        Upper = 3                                     Right-censored =     160

                                                   LR chi2(8)        = 1003.81
Log likelihood = -2638.6858                        Prob > chi2       =  0.0000



   accidents 
        IRR   Std. err.      z    P>|z|     [95% conf. interval]
 

      1.past 
   2.641695   .1967576    13.04   0.000     2.282884    3.056902
    1.parent 
   .8345776   .0425139    -3.55   0.000     .7552765     .922205
             

    ntickets 

          1  
   1.994213   .1182697    11.64   0.000     1.775374    2.240027
          2  
   3.841546   .2575073    20.08   0.000      3.36859    4.380906
          3  
   6.979123   .6090534    22.26   0.000     5.881909    8.281012
          4  
   15.97291   2.579757    17.16   0.000     11.63879    21.92099
          5  
    66.9069   2547.936     0.11   0.912     2.57e-31    1.74e+34
          6  
   58.24981   4426.516     0.05   0.957     1.20e-63    2.82e+66
             

       _cons 
   .3387513    .015292   -23.98   0.000     .3100673    .3700889

Note: _cons estimates baseline incidence rate.

We interpret the model coefficients (or incidence-rate ratios) as if the censoring had not occurred. That is to say, as though we had seen all of the data, uncensored.

We find that past accidents predict more future accidents, that being a parent predicts fewer future accidents, and that the number of tickets generally predicts more future accidents, although having just 1 or 2 tickets has little significance.

Because of the censoring, we do not know which of the people coded as having 3 accidents really had exactly 3 accidents, or which had more.

We can, however, now make predictions of the expected uncensored number of accidents and the probabilities of any specified number of accidents, including values greater than 3.

We wonder, what are the chances anyone had more than 3 accidents in our data? Our data were officially top-coded, but were they practically top-coded? We can obtain each driver's probability of having four or more accidents by typing

 . predict fourplus, pr(4,.)

We now have the probability that each driver in our sample had four or more accidents. To get the expected number of drivers who had 4 or more accidents, we simply sum these probabilities

. total fourplus

Total estimation                         Number of obs = 3,000




             
 
      Total   Std. err.     [95% conf. interval]

    fourplus 
 
   74.46338   5.614501      63.45472    85.47205

We expect 74.5 drivers in our data had more than 3 accidents, and top-coding almost certainly affected our data.

Almost certainly? We have a standard error above, but the standard error and confidence interval do not account for the probabilities having themselves been estimated. If we use margins to perform the computation, it will produce the correct standard error and confidence interval

. margins, expression(predict(pr(4,.))*3000)


Predictive margins                                       Number of obs = 3,000
Model VCE: OIM

Expression: predict(pr(4,.))*3000



             
 
            Delta-method
             
 
     Margin   std. err.      z    P>|z|     [95% conf. interval]

       _cons 
 
   74.46338   5.238347    14.22   0.000     64.19641    84.73036

margins wants to report a mean, so we had to trick it into giving us a total by multiplying the probabilities by our sample size of 3000.

With such a small standard error and a lower bound of 43.6 on our confidence interval, we can definitively say, or at least as definitively as any statistician can say, that top-coding affected our data.