Calcul de la valeur de p inconnu

Je déboguais récemment un script R et j'ai trouvé quelque chose de très étrange, l'auteur a défini sa propre fonction de valeur p

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

Où X et Y sont des valeurs positives supérieures à 0. Le cas <20 semble être un calcul pour une sorte de distribution hypergéométrique (quelque chose de similaire au test de Fisher?) Et quelqu'un sait-il quel est l'autre calcul? En tant que sidenote, j'essaie d'optimiser ce code afin d'essayer de comprendre la fonction R appropriée pour appeler et remplacer cela par.

Edit: La formule détaillant le papier pour le calcul de la valeur p peut être trouvée ici (il faut cliquer sur pdf pour voir les formules) Les méthodes commencent à la page 8 du pdf et la formule en question se trouve à la page 9 sous (1). La distribution qu'ils supposent est un Poisson.

La deuxième chose semble être une approximation du calcul utilisé pour le x+y < 20cas, mais basée sur l' approximation de Stirling .

Normalement, lorsqu'il est utilisé pour ce type d'approximation, les gens utiliseraient au moins le terme supplémentaire suivant (le facteur de dans l'approximation den! ), ce qui améliorerait sensiblement l'approximation relative pour les petitsn.$\sqrt{2\pi n}$$n!$$n$

Par exemple, si et y sont tous deux 10, le premier calcul donne environ 0,088 tandis que l'approximation lorsque le facteur de $x$$y$ est inclus dans tous les termes soit environ 0,089, assez proche pour la plupart des fins ... mais en omettant ce terme dans l'approximation donne 0,5 - ce qui n'est vraiment pas assez proche! L'auteur de cette fonction n'a clairement pas pris la peine de vérifier l'exactitude de son approximation au cas limite.$\sqrt{2\pi n}$

À cet effet, l'auteur aurait probablement dû simplement appeler la lgammafonction intégrée - en particulier, en l'utilisant à la place de ce qu'il avait pour log_p1:

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

lgamma(x+1)$\log(x!)$

De même, je ne sais pas pourquoi l'auteur n'utilise pas la choosefonction intégrée dans la première partie, une fonction qui vient dans la distribution standard de R. Pour cette question, la fonction de distribution pertinente est probablement également intégrée.

lgammachoosechoose(1000,500)lgamma$x$$y$

Avec plus d'informations, il devrait être possible d'identifier la source du test. Je suppose que l'écrivain l'a pris quelque part, il devrait donc être possible de le retrouver. Avez-vous un contexte pour cela?

Lorsque vous dites «optimiser», voulez-vous dire le rendre plus rapide, plus court, plus facile à entretenir ou autre chose?

Modifier après avoir lu rapidement sur le papier:

Les auteurs semblent se tromper sur un certain nombre de points. Le test exact de Fisher ne suppose pas que les marges sont fixes, il les conditionne simplement , ce qui n'est pas du tout la même chose, comme discuté, par exemple, ici , avec des références. En effet, ils semblent à peu près complètement ignorants du débat sur le conditionnement sur les marges et pourquoi cela se fait. Les liens qui s'y trouvent méritent d'être lus.

[Ils passent du «test de Fisher est toujours plus conservateur que le nôtre» à l'affirmation selon laquelle le test de Fisher est trop conservateur ... ce qui ne suit pas nécessairement à moins qu'il ne soit mauvais de conditionner . Ils devraient établir cela, mais étant donné que c'est quelque chose dont les statisticiens se disputent depuis environ 80 ans, et ces auteurs ne semblent pas savoir pourquoi le conditionnement est fait, je ne pense pas que ces gars-là aient tout à fait au fond de ce problème .]

Les auteurs de l'article semblent au moins comprendre que les probabilités qu'ils donnent doivent être cumulées pour donner des valeurs de p; par exemple vers le milieu de la première colonne de la page 5 (c'est moi qui souligne):

La signification statistique selon le test exact de Fisher pour un tel résultat est de 4,6% (valeur P bilatérale, c'est-à-dire la probabilité qu'un tel tableau se produise dans l'hypothèse que les fréquences d'actine EST sont indépendantes des bibliothèques d'ADNc). En comparaison, la valeur P calculée à partir de la forme cumulative (équation 9, voir Méthodes) de l'équation 2 (c'est-à-dire, pour que la fréquence relative des EST d'actine soit la même dans les deux bibliothèques, étant donné qu'au moins 11 EST apparentés sont observés dans la bibliothèque hépatique après deux observations dans la bibliothèque cérébrale) est de 1,6%.

(même si je ne suis pas sûr d'être d'accord avec leur calcul de la valeur là-bas; je devrais vérifier attentivement pour voir ce qu'ils font réellement avec l'autre queue.)

Je ne pense pas que le programme fasse cela.

$x$$x+y$

Je ne suis même pas convaincu que la somme de leurs probabilités soit 1 à ce stade.

Il y a beaucoup plus à dire ici, mais la question ne concerne pas le document, c'est la mise en œuvre dans le programme.

-

Quoi qu'il en soit, le résultat est qu'au moins l'article identifie correctement que les valeurs de p consistent en une somme de probabilités comme celles de l'équation 2, mais le programme ne le fait pas . (Voir les eqn 9a et 9b dans la section Méthodes du document.)

Le code est tout simplement faux à ce sujet.

[Vous pouvez utiliser pbinom, comme le laisse entendre le commentaire de @ whuber, pour calculer les probabilités individuelles (mais pas la queue, car ce n'est pas un test binomial tel qu'il le structure), mais il y a un facteur supplémentaire de 1/2 dans leur équation 2, donc si vous souhaitez reproduire les résultats dans le document, vous devez les modifier.]

Vous pouvez l'obtenir, avec quelques bidouilles, auprès de pnbinom-

$k^\rm{th}$$k^\rm{th}$

$p = N_1/(N_1+N_2)$$k=x$$r = y+1$

$y$

Ce serait mauvais.