Unbekannte p-Wert-Berechnung

Ich habe kürzlich ein R-Skript debuggt und fand etwas sehr Seltsames. Der Autor hat seine eigene p-Wert-Funktion definiert

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

Wobei X und Y positive Werte größer als 0 sind. Der Fall <20 scheint eine Berechnung für eine Art hypergeometrische Verteilung zu sein (ähnlich dem Fisher-Test?) Und weiß jemand, was die andere Berechnung ist? Als Nebenbemerkung versuche ich, diesen Code zu optimieren, um die richtige R-Funktion zu finden, mit der dieser aufgerufen und ersetzt werden kann.

Bearbeiten: Die Papierdetailformel für die p-Wert-Berechnung finden Sie hier (klicken Sie auf pdf, um die Formeln anzuzeigen). Die Methoden beginnen auf Seite 8 des PDF, und die betreffende Formel finden Sie auf Seite 9 unter (1). Die Verteilung, die sie annehmen, ist ein Poisson.

Das zweite sieht so aus, als wäre es eine Annäherung an die Berechnung, die für den x+y < 20Fall verwendet wird, aber basierend auf der Stirling-Näherung .

Normalerweise würden Menschen, wenn es für diese Art der Approximation verwendet wird, mindestens den nächsten zusätzlichen Term verwenden (den Faktor in der Approximation für ), Was die relative Approximation für kleines erheblich verbessern würde . n$\sqrt{2\pi n}$$n!$$n$

Wenn zum Beispiel und beide 10 sind, ergibt die erste Berechnung ungefähr 0,088, während die Annäherung, wenn der Faktor in allen Begriffen enthalten ist, ungefähr 0,089 ist, was für die meisten Zwecke nahe genug ist ... aber wenn man diesen Term in der Näherung weglässt, ergibt sich 0,5 - was wirklich nicht nah genug ist! Der Autor dieser Funktion hat sich offensichtlich nicht die Mühe gemacht, die Genauigkeit seiner Annäherung im Grenzfall zu überprüfen.$x$$y$$\sqrt{2\pi n}$

Zu diesem Zweck hätte der Autor wahrscheinlich einfach die eingebaute lgammaFunktion aufrufen sollen - insbesondere, indem er diese anstelle dessen verwendet, wofür er hat log_p1:

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

was zu der Antwort führt, die er zu approximieren versucht (da lgamma(x+1)tatsächlich , genau das, was er versucht - schlecht - über die Stirling-Approximation zu approximieren).$\log(x!)$

Ebenso bin ich mir nicht sicher, warum der Autor die eingebaute chooseFunktion im ersten Teil nicht verwendet, eine Funktion, die in der Standardverteilung von R enthalten ist. In diesem Fall ist wahrscheinlich auch die relevante Verteilungsfunktion eingebaut.

Sie brauchen nicht wirklich zwei getrennte Fälle; Der lgammaeine funktioniert einwandfrei bis auf die kleinsten Werte. Auf der anderen Seite choosefunktioniert die Funktion für ziemlich große Werte (z. B. choose(1000,500)funktioniert sie einwandfrei). Die sicherere Option ist wahrscheinlich lgamma, obwohl Sie ziemlich große und haben müssten, bevor es ein Problem war.$x$$y$

Mit mehr Informationen sollte es möglich sein, die Quelle des Tests zu identifizieren. Ich vermute, der Autor hat es von irgendwoher genommen, also sollte es möglich sein, es aufzuspüren. Haben Sie einen Kontext dafür?

Wenn Sie "Optimieren" sagen, meinen Sie damit, es schneller, kürzer, wartbarer oder etwas anderes zu machen?

Bearbeiten Sie nach schnellem Lesen des Papiers:

Die Autoren scheinen in einigen Punkten falsch zu liegen. Der genaue Test von Fisher geht nicht davon aus, dass die Ränder fest sind, sondern stellt lediglich Bedingungen an sie, was überhaupt nicht dasselbe ist, wie zum Beispiel hier mit Referenzen besprochen . In der Tat scheinen sie die Debatte über die Konditionierung von Margen und warum dies getan wird, so gut wie nicht zu kennen. Die Links dort sind lesenswert.

[Sie gehen von "Fischers Test ist immer konservativer als unser" zu der Behauptung über, dass Fischers Test zu konservativ ist ... was nicht unbedingt folgt, es sei denn, es ist falsch zu konditionieren . Sie müssten das feststellen, aber angesichts der Tatsache, dass Statistiker seit etwa 80 Jahren darüber streiten und diese Autoren nicht zu wissen scheinen, warum Konditionierung durchgeführt wird, glaube ich nicht, dass diese Leute diesem Problem auf den Grund gegangen sind .]

Die Autoren des Papiers scheinen zumindest zu verstehen, dass die von ihnen angegebenen Wahrscheinlichkeiten kumuliert werden müssen, um p-Werte zu erhalten. Zum Beispiel in der Mitte der ersten Spalte von Seite 5 (Hervorhebung von mir):

Die statistische Signifikanz gemäß dem exakten Fisher-Test für ein solches Ergebnis beträgt 4,6% (zweiseitiger P-Wert, dh die Wahrscheinlichkeit, dass eine solche Tabelle in der Hypothese auftritt, dass die Aktin-EST-Frequenzen unabhängig von den cDNA-Bibliotheken sind). Im Vergleich dazu ist der aus der kumulativen Form (Gleichung 9, siehe Methoden) von Gleichung 2 berechnete P-Wert (dh, dass die relative Häufigkeit von Aktin-ESTs in beiden Bibliotheken gleich ist, vorausgesetzt, dass in mindestens 11 verwandte ESTs beobachtet werden Die Leberbibliothek, nachdem zwei in der Gehirnbibliothek beobachtet wurden, beträgt 1,6%.

(obwohl ich nicht sicher bin, ob ich mit der Berechnung des Werts dort einverstanden bin; ich müsste sorgfältig prüfen, was sie tatsächlich mit dem anderen Schwanz machen.)

Ich glaube nicht, dass das Programm das tut.

Beachten Sie jedoch, dass ihre Analyse kein Standard-Binomialtest ist. Sie verwenden ein Bayes'sches Argument, um einen p-Wert in einem ansonsten häufig auftretenden Test abzuleiten. Sie scheinen auch - meiner Meinung nach etwas seltsam - eher auf als auf . Das bedeutet, dass sie eher ein negatives Binom als ein Binom haben müssen, aber ich finde das Papier wirklich schlecht organisiert und schrecklich schlecht erklärt (und ich bin es gewohnt, herauszufinden, was in Statistikpapieren vor sich geht), also bin ich es Ich werde nicht sicher sein, wenn ich nicht sorgfältig durchgehe.$x$$x+y$

Ich bin nicht einmal davon überzeugt, dass die Summe ihrer Wahrscheinlichkeiten zu diesem Zeitpunkt 1 ist.

Hier gibt es noch viel mehr zu sagen, aber die Frage bezieht sich nicht auf das Papier, sondern auf die Implementierung im Programm.

- -

Das Ergebnis ist jedenfalls, dass zumindest das Papier korrekt identifiziert, dass p-Werte aus einer Summe von Wahrscheinlichkeiten wie denen in Gleichung 2 bestehen, das Programm jedoch nicht . (Siehe Gleichung 9a und 9b im Abschnitt Methoden des Papiers.)

Der Code ist einfach falsch.

[Sie könnten verwenden pbinom, wie @ whubers Kommentar implizieren würde, um die einzelnen Wahrscheinlichkeiten zu berechnen (aber nicht den Schwanz, da es sich nicht um einen Binomialtest handelt, da sie ihn strukturieren), aber dann gibt es in ihrer Gleichung 2 einen zusätzlichen Faktor von 1/2 Wenn Sie die Ergebnisse im Papier replizieren möchten, müssen Sie sie ändern.]

Sie können es mit einigem Fummeln erhalten von pnbinom-

Die üblichen Formen des negativen Binomials sind entweder die Anzahl der Versuche zum Erfolg von oder die Anzahl der Fehler zum Erfolg von . Die beiden sind gleichwertig; Wikipedia gibt das zweite Formular hier . Die Wahrscheinlichkeitsfunktion ist: k t $k^\rm{th}$$k^\rm{th}$

Die Gleichung 2 auf p4 (und damit auch Gleichung 1 auf p3) ist ein negatives Binom, aber um 1 verschoben . Sei , und .$p = N_1/(N_1+N_2)$$k=x$$r = y+1$

Dies macht mir Sorgen, dass, da die Grenzen für nicht ähnlich verschoben wurden, ihre Wahrscheinlichkeiten möglicherweise nicht einmal zu 1 addieren.$y$

Das wäre schlecht