20
$\begingroup$

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?

$\endgroup$
19
$\begingroup$

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.

$\endgroup$
  • 4
    $\begingroup$ This is a useful and informative explanation, thank you (+1). The analysis in my answer suggests the "correct matching p-values" are too high, though, and that employing them will (unnecessarily) lose power. The consequences are greater than they might appear: when a p-value changes from, say, 0.05 to 0.08 (roughly the magnitude in this example), it indicates you need about a third more data to achieve the same power -- and that could blow your sampling budget. $\endgroup$ – whuber Sep 23 '20 at 19:40
  • 2
    $\begingroup$ Yes. The difficult question is judging whether, with the null hypothesis $\lambda=5.22$, whether an outcome of $1$ is more or less extreme than an outcome of $10$ ($0$ clearly is, as is $11$ or anything greater). The probability of observing $1$ is greater than the probability of $10$, but the probability of $\{0,1\}$ is lower than the probability of $\{10,11,12,\ldots\}$ $\endgroup$ – Henry Sep 24 '20 at 0:13
9
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Why does the correct confidence interval have the same mass is each tail? Hahn and Meeker have 2 paragraphs in their section 2.7 discussing it, but it's light on details. $\endgroup$ – user257566 Sep 24 '20 at 18:27
  • 1
    $\begingroup$ @user257566 In the sense that there exist many valid confidence interval procedures in this setting, there is no single best one. But my use of the word "correct" was intended in the sense of this is the very one the OP requested by their use of poisson.test, that's all. $\endgroup$ – whuber Sep 24 '20 at 18:30
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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