Unknown p-value calculation

Question

I was recently debugging an R script and I found something very weird, the author defined their own p-value function

pval <- function(x, y){
    if (x+y<20) { # x + y is small, requires R.basic
        p1<- nChooseK(x+y,x) * 2^-(x+y+1);
        p2<- nChooseK(x+y,y) * 2^-(x+y+1);
        pvalue = max(p1, p2)
    }
    else { # if x+y is large, use approximation
        log_p1 <- (x+y)*log(x+y) - x*log(x) - y*log(y) - (x+y+1)*log(2);
        pvalue<-exp(log_p1);
    }
    return(pvalue)
}

Where X and Y are values positive values greater than 0. The <20 case seems to be a calculation for some kind of hypergeometric distribution (something similar to Fisher test?) and does anyone know what the other calculation is? As a sidenote, I am trying to optimize this code so trying to figure out the proper R function to call and replace this with.

Edit: Paper detailing formula for p-value calculation can be found here (need to click pdf to see the formulas) Methods start on page 8 of the pdf and the formula in question can be found on page 9 under (1). The distribution they assume is a Poisson.

Answer 1

The second thing looks like it is an approximation to the calculation being used for the x+y < 20 case, but based off the Stirling approximation.

Normally when it's being used for this sort of approximation, people would use at least the next additional term (the factor of $\sqrt{2\pi n}$ in the approximation for $n!$), which would improve the relative approximation substantially for small $n$.

For example, if $x$ and $y$ are both 10, the first calculation gives about 0.088 while the approximation when the factor of $\sqrt{2\pi n}$ is included in all the terms is about 0.089, close enough for most purposes... but omitting that term in the approximation gives 0.5 - which is really not close enough! The author of that function clearly hasn't bothered to check the accuracy of his approximation at the boundary case.

For this purpose, the author should probably simply have called the built in lgamma function - specifically, by using this instead of what he has for log_p1:

log_p1 <- lgamma(x+y+1)-lgamma(x+1)-lgamma(y+1)-(x+y+1)*log(2)

which results in the answer he's trying to approximate (since lgamma(x+1) actually returns $\log(x!)$, the very thing he's trying to approximate - poorly - via Stirling approximation).

Similarly, I'm not sure why the author doesn't use the built in choose function in the first part, a function which comes in the standard distribution of R. For that matter, the relevant distribution function is probably built-in too.

You don't really need two separate cases; the lgamma one works just fine right down to the smallest values. On the other hand, the choose function works for pretty big values (e.g. choose(1000,500) works just fine). The safer option is probably lgamma, though you'd need to have quite large $x$ and $y$ before it was an issue.

With more information it should be possible to identify the source of the test. My guess is the writer has taken it from somewhere, so it should be possible to track it down. Do you have some context for this?

When you say 'optimize' do you mean make it faster, shorter, more maintainable, or something else?

---

Edit after quickly reading over the paper:

The authors appear to be wrong on a number of points. Fisher's exact test doesn't assume the margins are fixed, it simply conditions on them, which is not the same thing at all, as discussed, for example, here, with references. Indeed, they seem pretty much completely ignorant of the debate over conditioning on margins and why it is done. The links there are worth reading.

The paper's authors do at least seem to understand that the probabilities they give must be cumulated to give p-values; for example near the middle of the first column of page 5 (emphasis mine):

The statistical significance according to Fisher’s exact test for such a result is 4.6% (two-tail P-value, i.e., the probability for such a table to occur in the hypothesis that actin EST frequencies are independent of the cDNA libraries). In comparison, the P-value computed from the cumulative form (Equation 9, see Methods) of Equation 2 (i.e., for the relative frequency of actin ESTs to be the same in both libraries, given that at least 11 cognate ESTs are observed in the liver library after two were observed in the brain library) is 1.6%.

(though I am not sure I agree with their calculation of the value there; I'd have to check carefully to see what they're actually doing with the other tail.)

I don't think the program does that.

Beware, however, that their analysis is not a standard binomial test; they use a Bayesian argument to derive a p-value in an otherwise frequentist test. They also seem - somewhat oddly, to my mind - to condition on $x$, rather than $x+y$. This means that they must end up with something like a negative binomial rather than a binomial, but I find the paper really badly organized and terribly badly explained (and I'm used to working out what's going on in statistics papers), so I am not going to be certain unless I go through carefully.

I'm not even convinced that the sum of their probabilities are 1 at this point.

There's a lot more to be said here, but the question isn't about the paper, it's about the implementation in the program.

--

Anyway, the upshot is, at least the paper correctly identifies that p-values consist of a sum of probabilities like those in equation 2, but the program doesn't. (See eqn 9a and 9b in the Methods section of the paper.)

The code is simply wrong on that.

[You could use pbinom, as @whuber's comment would imply, to work out the individual probabilities (but not the tail, since it's not a binomial test as they structure it) but then there's an extra factor of 1/2 in their equation 2 so if you want to replicate the results in the paper, you need to alter them.]

You can obtain it, with some fiddling, from pnbinom -

The usual forms of the negative binomial are either the number of trials to the $k^\rm{th}$ success or the number of failures to the $k^\rm{th}$ success. The two are equivalent; Wikipedia gives the second form here. The probability function is:

$$ {k+r-1 \choose k}\cdot (1-p)^r p^k,\! $$

The equation 2 on p4 (and so also eqn 1 on p3) is a negative binomial, but shifted by 1. Let $p = N_1/(N_1+N_2)$, $k=x$ and $r = y+1$.

This makes me concerned that since the limits on $y$ have not been similarly shifted, that their probabilities may not even add to 1.

That would be bad.