0
$\begingroup$

I am struggling with calculating a confidence interval for my data using a Poisson distribution. The context is manufacturing process quality control, and I think that a Poisson distribution is the correct model for my process.

My confusion comes from the fact that my sample observations are not all the same size, and I don't understand how to account for this. Here is a simplified example.

My factory grinds flour. During the grinding, sometimes rocks fall in the flour. This occurs at random time intervals.

I have 500 samples randomly taken from that process, and I found a total of 7 rocks. I would like to calculate a confidence interval for the number of rocks per 1000 kg of flour.

If all the samples were 1 kg, I would just use the formula for an exact Poisson confidence interval. The rate is low, so my understanding is that I have to use the exact form rather than a normal approximation. For 7 observations out of 500 samples, I would calculate the mean as $7/500=0.0140$ with a 95% confidence interval of $(0.0056, 0.0288)$. Converted to a per 1000 kg basis, this would be 5.6 to 28.8 rocks per 1000 kg.

However, in my case, some of the samples are 1 kg, and some of them are 2 kg. Here's a summary of my data:

\begin{matrix} \text{Sample Number} & \text{Size} & \text{Particle Count} \\ 1 & 1 \text{ kg} & 2 \\ 2 & 1 \text{ kg} & 1 \\ 3 & 1 \text{ kg} & 1 \\ 4 & 2 \text{ kg} & 1 \\ 5 & 2 \text{ kg} & 1 \\ 6 & 2 \text{ kg} & 1 \\ 7 \text{ to } 300 & 1 \text{ kg} & 0 \\ 301 \text{ to } 500 & 2 \text{ kg} & 0 \\ \end{matrix}

So, the observed occurrence rate for rocks is 7 rocks per 703 kg of flour, or 9.96 rocks per 1000 kg. How do I calculate a 95% confidence interval for this mean? I considered just treating each 2 kg sample as two separate 1 kg samples, but I am not sure if this is really valid.

I am somewhat familiar with the concept of "offsets" in Poisson regression, but I don't understand how to apply that for calculating a confidence interval.

Improve this question
$\endgroup$
1
  • $\begingroup$ Can you post a link to your data? then we can see. $\endgroup$
    – kjetil b halvorsen ♦
    Commented Mar 17, 2018 at 1:40

1 Answer 1

Reset to default
1
$\begingroup$

If you use a log(kg of flower) offset in Poisson regression, then you will assume that $$\log E Y_i = \log \mu + \log \text{weight}_i $$ or $$\log E(Y_i / \text{weight}_i) = \log \mu.$$ I.e. the regression coefficient for the intercept of a model with an intercept ($ \log \mu$) and a log weight offset is giving you the log-mean rate per kg. If you want the mean rate of events per kg, you exponentiate the estimated coefficient of the intercept.

Most software will also automatically give you a standard error and confidence interval (for a 95% CI approximately estimated coefficient $\pm$ 1.96 $\times$ SE) for the intercept. To get this on the mean rate scale you can just exponentiate the confidence interval limits.

In fact, if we assume as above that there's no difference between the batches you can do the whole calculation by hand. The estimated mean rate per kg is $$\hat{\mu}=(\sum_i y_i)/(\sum_i \text{weight}_i)$$ with 95% CI limits $$\exp ( \log(\sum_i y_i) - \log(\sum_i \text{weight}_i) \pm 1.96 \times 1/\sqrt{\sum_i y_i} ) .$$

Improve this answer
$\endgroup$
6
  • $\begingroup$ True, but at least that approximation is a lot better on the log scale. Doing a regression model and getting a profile likelihood CI avoids that. $\endgroup$
    – Björn
    Commented Mar 19, 2018 at 16:16
  • $\begingroup$ Using $1.96 \times 1 / \sqrt{\sum{y_I}}$ is the normal approximation. I don't think it is valid for small values of $\bar{y}$. I am trying use the exact confidence interval (Garwood, 1936). $\endgroup$
    – James Duvall
    Commented Mar 19, 2018 at 16:39
  • $\begingroup$ You can also easily take a Bayesian approach, if your prior information can be expressed in terms of a Gamma(prior shape, prior rate) distribution. The posterior (assuming a single rate) is a Gamma(prior shape+total particle count, prior rate+total kg). Almost any software can calculate the percentiles (e.g. 2.5th and 97.5th) of a gamma distribution for you. You might arrive at a prior by assuming one equivalent to 0.01 particles in 1 kg (i.e. corresponding to very little information). If a single Gamma is not enough to express your prior, mixtures of gamma also allow such conjugate updating.. $\endgroup$
    – Björn
    Commented Mar 19, 2018 at 16:40
  • $\begingroup$ To give a R example with your example data: the lower bound of a percentile credible interval is qgamma(0.025, shape=0.01+7, rate=1+(3+6+294*1+200*2)) = 0.0040 and the upper bound is qgamma(0.975, shape=0.01+7, rate=1+(3+6+294*1+200*2)) = 0.0186. $\endgroup$
    – Björn
    Commented Mar 19, 2018 at 16:46
  • $\begingroup$ Thanks, I will have to do some reading on that approach to understand it I think. $\endgroup$
    – James Duvall
    Commented Mar 19, 2018 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.