Solved – Confidence interval for mean of Poisson with only zero counts

Question

I have 50 sample units in group A all with counts of zero. I have other groups that have some non-zero counts.

I would like to make a 95% confidence interval for the mean count per sample unit within group A.

Is there a confidence interval procedure for a Poisson variable that can create a confidence interval with an upper bound that is greater than zero?

In other words, if you haven't observed anything 50 times, how big could the Poisson mean be so that seeing 50 zeroes would be reasonable?

I am assuming that with further sampling some observation would eventually result in a non-zero count in the group of interest.

Answer 1

The standard procedure (Hahn & Meeker, section 7.2.2) exploits the basic relationship between Poisson and Chi-squared variates; namely, when $F_{\lambda}$ is the Poisson PDF of parameter $\lambda$ and $G_{\nu}$ is the Chi-squared PDF of parameter $\nu,$ then for any $k\in\{0,1,2,\ldots\},$

$$1-F_\lambda(k) = G_{2k+2}(2\lambda).\tag{1}$$

An upper confidence limit of size $1-\alpha$ for $\lambda$ based on observing a Poisson variable $K_\lambda$ is, by definition, a function $u$ for which

$$1-\alpha = \inf_{\lambda\in\mathbb{R}^+}\Pr(\lambda \le u(K_\lambda)).$$

If we choose a suitable inverse of $u$ and write $k=K_\lambda$ for the observed value, we may exploit $(1)$ to re-express this criterion as

$$\eqalign{ 1-\alpha &= \inf_{\lambda\in\mathbb{R}^+}\Pr(u^{-1}(\lambda) \le K_\lambda) \\ &= \inf_{\lambda\in\mathbb{R}^+}1-F_\lambda(k)) \\ &= \inf_{\lambda\in\mathbb{R}^+}G_{2k+2}(2\lambda), }$$

with unique solution

$$\lambda_+(\alpha) = \frac{1}{2} G^{-1}_{2k+2}(1-\alpha).$$

Similar reasoning arrives at a lower $1-\alpha$ confidence limit

$$\lambda_{-}(\alpha) = \frac{1}{2} G^{-1}_{2k}(\alpha).$$

One of the many possible two-sided confidence interval procedures splits the risk between the upper and lower endpoints by using $[\lambda_{-}(\alpha/2), \lambda_{+}(\alpha/2)].$

---

When $k=0,$ the function $G_{0},$ or the distribution of a "chi-squared variate with zero degrees of freedom," has to be understood as the distribution of the constant zero, whence "$G^{-1}_0(\alpha)$" is always zero no matter what $\alpha\gt 0$ may be. In this case $G_{2k+2} = G_2$ is the Exponential distribution with scale factor $2,$ entailing

$$\lambda_{+}(\alpha/2) = G^{-1}_2(1-\alpha/2) = -2\log(\alpha/2).$$

For instance, with $\alpha=5\%$ this UCL is $7.38,$ whereas the one-sided upper confidence limit for the same $\alpha$ is only $3.00.$ If you are tempted to use the latter because it produces a shorter confidence interval, consider these simulation results for a large range of $\lambda$ (from $0.1$ to $1,000,$ after which a Normal approximation will work well):

"Coverage" is the proportion of samples for which the confidence interval, nominally set at $1-\alpha = 95\%,$ includes $\lambda.$ Each red point in this plot summarizes 400,000 independently simulated samples. The gray graph is the calculated coverage based on Poisson probabilities only.

The discreteness of the Poisson distributions causes the actual coverage to oscillate, but a trend is clear: coverage really is close to the nominal value for large $\lambda,$ but can be substantially greater for small $\lambda.$

Some of the conclusions we may draw are

1. The foregoing analysis produces confidence intervals with the correct coverage.

2. The coverage tends to be higher than intended (greater than $1-\alpha$) when $\lambda$ is smaller than $10$ or so, approaching $100\%$ in the limit as $\lambda\to 0.$

In retrospect this behavior is obvious: because the confidence limits depend only on $k,$ the limits for $k=0$ have to be fairly large to allow for the possibility that $\lambda$ is fairly large. Consequently, when $\lambda$ actually is small, the coverage must be greater than the nominal coverage.

If you know (or assume) $\lambda$ is small at the outset, you could modify this procedure accordingly to produce confidence intervals that tend to be shorter.

---

Reference

G. J. Hahn and W. Q. Meeker (1991), Statistical Intervals. A Guide for Practitioners. J. Wiley & Sons.

---

Code

#
# Poisson confidence intervals (symmetric, two-sided).
# `k` may be a vector of observations.
#
ci <- function(k, alpha=0.05) {
  matrix(qchisq(c(alpha/2, 1-alpha/2), rbind(2*k, 2*k+2))/2, 2)
}
#
# Simulation study of coverage.
# Takes a few seconds with n=4e5.
#
n <- 4e5
lambda <- 10^seq(-1, 3, length.out=21)
set.seed(17)
coverage <- sapply(lambda, function(lambda) {
  mean((function(x) x[1,] <= lambda & lambda <= x[2,])(ci(rpois(n, lambda))))
})
#
# Calculation of coverage.
#
lambda.calc <- 10^seq(-1, 3, length.out=4021)
x <- max(lambda.calc)
CI <- ci(k <- 0:(x + 8*sqrt(x)))
coverage.calc <- sapply(lambda.calc, function(l) {
  covers <- CI[1,] <= l & l <= CI[2,]
  sum(dpois(k, l)[covers])
})
#
# Plot of results.
#
library(ggplot2)
ggplot(data.frame(lambda=lambda, Coverage=coverage), 
       aes(lambda, Coverage)) + 
  geom_line(data=data.frame(lambda=lambda.calc, Coverage=coverage.calc), col="#a0a0a0") + 
  geom_point(color="Red") + 
  scale_x_log10() + 
  coord_cartesian(ylim=c(min(0.9499, min(coverage.calc)), 1), expand=FALSE) + 
  geom_hline(yintercept=0.95) + 
  xlab(expression(lambda)) + 
  ggtitle("Simulated Coverage Rates of 95% Two-Sided Poisson Confidence Intervals")