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")

Answer 2

I answered my own question after some research. Please comment if something is wrong.

Frequentist confidence interval

An exact confidence interval can be derived based on the probability mass function for a Poisson distribution,

$$ P(X \le k) = \frac{\lambda^ke^{-\lambda}}{k!} $$ in which $k$ is some possible count and $\lambda$ is the mean and variance.

In our case, with a count of zero, $k = 0$, so $$ P(X = 0) = \frac{\lambda^0e^{-\lambda}}{0!}=e^{-\lambda} $$ For a given level of confidence, $1-\alpha$, we can use this to solve for an upper bound on the estimate of $\lambda$: $$ \alpha = e^{-\lambda}$$ $$ log(\alpha) = log(e^{-\lambda})$$ $$ log(\alpha) = -\lambda$$ $$ \hat{\lambda}_{upper} = -log(\alpha)$$ For a 95% confidence interval the upper confidence limit for a single observation of a Poisson random variable with a count of zero is $-log(.05)=2.995732$, and for a 99% confidence interval the upper limit is $-log(.01)=4.60517$.

The sum of multiple Poisson random variables is also a Poisson random variable, with mean $n\lambda$, so to convert these confidence limits for our case with 50 observations, we can simply divide by 50.

$$UCL_{95}=2.995732/50=0.05991464$$ $$UCL_{99}=4.60517/50=0.0921034$$

Bayesian credible interval

The likelihood is from a Poisson distribution: $$L(\lambda|x)=\prod_{i=1}^n\frac{e^{-\lambda}\lambda^{x_i}}{x_i!}=\frac{e^{-n\lambda}\lambda^{\sum x_i}}{\prod_{i=1}^n(x_i!)} $$

If you make the prior a $gamma(\alpha,\beta)$: $$p(\lambda)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}\lambda^{\alpha-1}e^{-\beta\lambda} $$

then the posterior is a $gamma(\sum x_i + \alpha, n+\beta) $:

$$p(\lambda|x)=\frac{p(x|\lambda)p(\lambda)}{p(x)}\propto p(x|\lambda)p(\lambda) $$

$$p(\lambda|x) \propto e^{-n\lambda}\lambda^{\sum x_i}\lambda^{\alpha-1}e^{-\beta\lambda} = \lambda^{\sum x_i + \alpha - 1} e^{-(n+\beta)\lambda} $$

If you use $\alpha = 1$ and $\beta = 0$ in the prior then the upper credible limits are the same as the upper confidence limits.

For this particular case, we can use the 95th and 99th percentiles of a $gamma(1, 50)$ distribution to get the upper limits of 95% and 99% credible intervals (remember the $\sum x_i = 0$ for this particular case).

In R you could use:

> qgamma(.95, shape = 1, rate = 50)
[1] 0.05991465
> qgamma(.99, shape = 1, rate = 50)
[1] 0.0921034

Difference in approaches

The two approaches result in the same interval bounds, but they have different estimates.

For the frequentist approach, the point estimate for the Poisson mean is the maximum likelihood estimate, which for a Poisson distribution is just the average of the sample: $$\hat{\lambda}=\frac{0}{50}= 0$$

For the Bayesian approach, the point estimate for the Poisson mean is the mean of the posterior distribution, which is a $gamma(1, 50)$: $$\hat{\lambda}=\frac{\alpha}{\beta}=\frac{1}{50}= 0.02$$