7
$\begingroup$

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.

$\endgroup$
2

2 Answers 2

4
$\begingroup$

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

Figure

"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")
$\endgroup$
8
  • $\begingroup$ Great answer. Would be possible to examine the R implementation and learn from it? Many thanks. $\endgroup$
    – Maximilian
    Sep 19, 2019 at 19:21
  • 2
    $\begingroup$ @Max I added the code at the end. Cheers! $\endgroup$
    – whuber
    Sep 19, 2019 at 19:24
  • $\begingroup$ @whuber, I am not understanding, because if I change your code for the ci function from c(alpha/2, 1-alpha/2) to c(0, 1-alpha), to make it one-sided, I still get coverage that never goes below 95%. $\endgroup$
    – Jdub
    Sep 19, 2019 at 21:00
  • 1
    $\begingroup$ I just checked an ad hoc modification of ci whereby when $k=0$ the CI is replaced by $[0, u]$ where $u$ is the one-sided UCL of level $\alpha.$ For $\alpha=0.05$ and values of $\lambda$ between $3.0$ and $3.44,$ the coverage is lower than $1-\alpha$. However, because it's not that much lower--it only drops to $93.8\%$--this is a procedure worthy of consideration. With more sophisticated adjustments (for other small values of $k$) you could do even better if you're willing to sacrifice the guarantee of $1-\alpha$ coverage. $\endgroup$
    – whuber
    Sep 19, 2019 at 21:14
  • 1
    $\begingroup$ Related question: stats.stackexchange.com/questions/428048 $\endgroup$
    – Jdub
    Sep 20, 2019 at 15:12
1
$\begingroup$

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$$

$\endgroup$
3
  • 1
    $\begingroup$ What about poisson.test(0, 50) in R? This gives a 95%-CI of $(0.000, 0.074)$. My intuition is that the lower limit should be $0$ in the case of $0$ counts. $\endgroup$ Sep 19, 2019 at 15:21
  • 1
    $\begingroup$ @COOLSerdash, that is an interesting observation. That is actually giving you a 97.5% confidence interval: -log(0.025)/50 = 0.07377759. The code does not realize it is up against a zero boundary where a two-sided confidence interval does not make sense. $\endgroup$
    – Jdub
    Sep 19, 2019 at 15:28
  • 2
    $\begingroup$ I think you need to be a little careful--that reasoning about a zero boundary doesn't hold up. $\endgroup$
    – whuber
    Sep 19, 2019 at 18:50

Your Answer

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

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