Conflict between Poisson confidence interval and p-value

Question

Testing whether the outcome of $x=10$ counts is compatible with a rate of $\lambda=5.22$ in R:

> poisson.test(x=10,r=5.22,alternative='two.sided')

Exact Poisson test

data:  10 time base: 1
number of events = 10, time base = 1, p-value = 0.04593
alternative hypothesis: true event rate is not equal to 5.22
95 percent confidence interval:
  4.795389 18.390356
sample estimates:
event rate 
        10 

This result leads to two contradictory conclusions:

1. The p-value is less than 0.05, which suggests that $\lambda\neq{5.22}$
2. However the 95% confidence interval is $[4.795389 < 5.22 < 18.390356]$, which keeps alive the hypothesis that $\lambda=5.22$

Thus this example violates the duality between hypothesis tests and confidence intervals. How is this possible?

Answer 1

There are several ways to define two-sided $p$-values in this case. Michael Fay lists three in his article. The following is mostly taken from his article.

Suppose you have a discrete test statistic $t$ with random variable $T$ such that larger values of $T$ imply larger values of a parameter of interest, $\theta$. Let $F_\theta(t)=\Pr[T\leq t;\theta]$ and $\bar{F}_\theta(t)=\Pr[T\geq t;\theta]$. Suppose the null value is $\theta_0$. The one-sided $p$-values are then denoted by $F_{\theta_0}(t), \bar{F}_{\theta_0}(t)$, respectively.

The three ways listed to define two-sided $p$-values are as follows:

$\textbf{central:}$ $p_{c}$ is 2 times the minimum of the one-sided $p$-values bounded above by 1: $$ p_c=\min\{1,2\times\min(F_{\theta_0}(t), \bar{F}_{\theta_0}(t))\}. $$

$\textbf{minlike:}$ $p_{m}$ is the sum of probabilities of outcomes with likelihoods less than or equal to the observed likelihood: $$ p_m=\sum_{T:f(T)\leq f(t)} f(T) $$ where $f(t) = \Pr[T=t;\theta_0]$.

$\textbf{blaker:}$ $p_b$ combines the probability of the smaller observed tail with the smallest probability of the opposite tail that does not exceed that observed probability. This may be expressed as: $$ p_b=\Pr[\gamma(T)\leq\gamma(t)] $$ where $\gamma(T)=\min\{F_{\theta_0}(T), \bar{F}_{\theta_0}(T))\}$.

If $p(\theta_0)$ is a two-sided $p$-value testing $H_0:\theta=\theta_0$, then its $100(1-\alpha)\%$ matching confidence interval is the smallest interval that contains all $\theta_0$ such that $p(\theta_{0})>\alpha$. The matching confidence limits to the $\textbf{central}$ test are $(\theta_{L},\theta_U)$ which are the solutions to: $$ \alpha/2=\bar{F}_{\theta_L}(t) $$ and $$ \alpha/2=F_{\theta_U}(t). $$

The contradiction arises because poisson.test returns $p_m$ ($\textrm{minlike}$) as the $p$-value but confidence limits that are based on the $\textrm{central}$ test!

The exactci package returns the correct matching $p$-values and confidence limits (you can set the method using the option tsmethod):

library(exactci)

poisson.exact(x=10, r=5.22, tsmethod = "central")

    Exact two-sided Poisson test (central method)

data:  10 time base: 1
number of events = 10, time base = 1, p-value = 0.08105
alternative hypothesis: true event rate is not equal to 5.22
95 percent confidence interval:
  4.795389 18.390356
sample estimates:
event rate 
        10 

Now there is no conflict between the $p$-value and the confidence intervals. In rare cases, even the exactci function will result in inconsistencies, which is mentioned in Michael Fays article.

Answer 2

The correct exact two-sided 95% confidence interval $[\lambda^{-},\lambda^{+}]$ is computed from an observation $x$ of a Poisson variable $X$ using the defining relationships

$$\Pr(X\lt x;\lambda^{-}) = \alpha/2$$

and

$$\Pr(X \gt x; \lambda^{+}) = 1 - \alpha/2.$$

We may find these limits by exploiting

$$e^{-\lambda}\sum_{i=0}^{x}\frac{\lambda^i}{i!} = F_{\text{Poisson}}(x;\lambda) = 1 - F_\Gamma(\lambda;x+1) = \frac{1}{x!}\int_\lambda^\infty t^x e^{-t}\,\mathrm{d}t$$

for natural numbers $x.$

(You can prove this inductively via repeated integrations by parts on the right hand side or you can observe that the left probability is the chance of observing $x$ or fewer points in a homogeneous, unit-rate Poisson process running for time $\lambda;$ while the right probability is the chance that its takes more than $\lambda$ time to observe the $x+1^\text{st}$ point -- which obviously is the same event.)

Thus, writing $G=F_\Gamma^{-1}$ for the Gamma quantile function, the confidence interval is

$$\left[G(\alpha/2;x), G(1-\alpha/2;x+1)\right].$$

The discreteness in the defining inequalities -- that is, the distinction between "$\lt$" and "$\le$" -- is to blame for the apparent inconsistency with the p-value. Indeed, in most circumstances replacing the lower limit by $G(\alpha/2,x+1)$ actually gives better coverage, as simulations show. Here, for instance, are simulations in R that estimate the coverages of these two procedures.

f <- function(x, alpha=0.05) qgamma(c(alpha/2, 1-alpha/2), c(x, x+1))
z <- 10
x <- matrix(rpois(2e6, f(z)), 2)
mean(x[1,] <= z & z <= x[2,])

The output, which is identical to that of poisson.test, will be close to 97.7% coverage. The altered interval is

f. <- function(x, alpha=0.05) qgamma(c(alpha/2, 1-alpha/2), x+1)
x <- matrix(rpois(2e6, f.(z)), 2)
mean(x[1,] <= z & z <= x[2,])

The output will be close to 96.3% coverage -- closer to the nominal 95% level.

The problem with this somewhat ad hoc modification is that it fails when the true rate is tiny. In the same simulation with a true rate of $1/10$ rather than $10,$ the coverage of the correct interval is around 98% but that of the modified interval is only 94.4%. If your objective is to achieve 95% or higher coverage--not going any lower--than this is unacceptable. For many applications, especially when very small values of the parameter are highly unlikely, the modified interval has much to recommend it and will produce results more consistent with the p value.

Reference

Hahn, GJ and WQ Meeker, Statistical Intervals. Wiley 1991.

Their formula (7.1), expressed in terms of quantiles of chi-squared distributions, is equivalent to the one I give in terms of Gamma distributions. (Chi-squared distributions with $2x$ degrees of freedom are scaled versions of Gamma distributions with $x$ degrees of freedom.)

Answer 3

There are two possibilities. The first, and most obvious, is that it is a bug. I looked up the documentation for poisson.test in R and, originally, it was a one-sided test. It did not support two-sided tests. The second would be that the p-value and the interval are using different loss functions, but I would suspect that is not the case. You should submit a bug report.