Sunday, March 30, 2014

Computing Confidence Interval for Poisson Mean

For Poisson distribution, there are many different ways for calculating the confidence interval. The paper by Patil and Kulkarni discusses 19 different ways to calculate a confidence interval for the mean of a Poisson distribution.

The most commonly used method is the normal approximation (for large sample size) and the exact method (for small sample size)

Normal Approximation:

For Poisson, the mean and the variance are both lambda (λ).
The standard error is calculated as: sqrt(
λ /n) where λ is Poisson mean and n is sample size or total exposure (total person years, total time observed,…)
The confidence interval can be calculated as:
                  λ ±z(α/2)*sqrt(λ/n).
The 95-percent confidence interval is calculated as: λ ±1.96*sqrt(λ/n).
The 99-percent confidence interval is calculated as: λ ±2.58*sqrt(λ/n).
  
EXACT method:

Refer to the following paper for the description of this method: 

The confidence interval for event X is calculated as: 

           (qchisq(α/2, 2*x)/2, qchisq(1-α/2, 2*(x+1))/2 )

Where x is the number of events occurred under Poisson distribution.

In order to calculate the exact confidence interval for Poisson mean, the obtained confidence interval for the number of events need to be converted to the confidence interval for Poisson mean.

Here are two examples from the internet:

Example 1:
Would like to know how confident I can be in my λ. Anyone know of a way to set upper and lower confidence levels for a Poisson distribution?
  • Observations (n) = 88
  • Sample mean (λ) = 47.18182
what would the 95% confidence look like for this?
With Normal Approximation, the 95% confidence interval is calculated as:

                47.18182 +/- 1.96* sqrt(47.18182/88)

This gives 45.7467, 48.617

With Exact method, we first need to calculate x (# of events):
                X = n * λ = 88 * 48.18182 = 4152

The compute the 95% confidence interval for X = 4152. This will give the 95% confidence interval for X as (4026.66, 4280.25)

The 95% confidence interval for mean (λ) is therefore:
lower bound = 4026.66 / 88 = 45.7575
upper bound = 4280.25 /88 = 48.6392

  
Say that 14 events are observed in 200 people studied for 1 year and 100 people studies for 2 years. Calculate the 95% confidence interval for Poisson mean

In this example, the number of events (X) is given, the Poisson rate (λ) or mean needs to be calculated.

First step is to calculate the person year:
The person time at risk is 200 + 100 x 2 = 400 person years
The poisson rate / poisson mean (λ) is :
  • Events observed = 14
  • Time at risk of event = 400
  • Poisson (e.g. incidence) rate estimate = 14/400 = 0.035

Normal Approximation: 95% confidence interval is calculated as: 
                0.035 +/- 1.96* sqrt(0.035/400)
This will give the 95% confidence interval of (0.0167, 0.0533)

Exact approach:  Calculate the 95% confidence interval for the number of events (X) using: 
(qchisq(0.025, 2*x)/2, qchisq(0.975, 2*(x+1))/2 )

and the result is: [7.65, 23.49]
Exact 95% confidence interval for Poisson mean is:
         Lower bound = 7.65 / 400 =0.019135 for lower bound and
         Upper bound = 23.49 / 400 = 0.058724 for upper bound


We will then say the Poisson mean is 0.035 with 95% confidence interval of (0.019, 0.059).

The following SAS programs can illustrate the calculations above: 

data normal;
  input lambda n ;
  lower = lambda - probit(0.975)*sqrt(lambda/n);
  upper = lambda + probit(0.975)*sqrt(lambda/n);
  datalines;
  47.18182 88
  0.035       400
;
proc print data=normal;
  title 'Normal Approximation for 95% confidence interval for Poisson mean';
run;

data exact;
  input X;
  lower = quantile('CHISQ',.025,2*x)/2;
  upper = quantile('CHISQ',.975,2*(x+1))/2;
  datalines;
 4152
 14
;
proc print data=exact;
  title 'Exact method for 95% confidence interval for Poisson mean';
run;

5 comments:

  1. Is there a way we can use proc genmod to generate 95% CI for incidence rate for total cohort? Thank you!

    ReplyDelete
  2. Yes, you can use proc genmod to calculate 95% confidence interval. The sample codes will be something like:

    proc genmod ;
    model eventcount = / dist=poisson link=log offset=logduration lrci alpha=0.05 ;
    run ;

    ReplyDelete
  3. Anonymous3:01 PM

    Thanks for your explanation! How do I calculate the Standard Error from the estimates of the confidence intervals obtained by the exact method?

    ReplyDelete
  4. Hi Sir,

    proc genmod ;
    model eventcount = / dist=poisson link=log offset=logduration lrci alpha=0.05 ;
    run ;

    For the above code, we need to exponentiate the output estimates right? Also event count variable contains the total event for a subject in the respective duration, where the log of this duration is passed as offset so that we get the rates using the property of the log function

    ReplyDelete