Negativ-Binomial-GLM vs. Log-Transformation für Zähldaten: erhöhte Typ-I-Fehlerrate

Einige von Ihnen haben vielleicht dieses schöne Papier gelesen:

O'Hara RB, Kotze DJ (2010) Zählungsdaten nicht protokollieren und transformieren. Methoden in Ökologie und Evolution 1: 118–122. klick .

In meinem Forschungsgebiet (Ökotoxikologie) beschäftigen wir uns mit schlecht replizierten Experimenten, und GLMs werden nicht häufig eingesetzt. Also habe ich eine ähnliche Simulation wie O'Hara & Kotze (2010) durchgeführt, aber ökotoxikologische Daten nachgeahmt.

Leistungssimulationen :

Ich simulierte Daten aus einem faktoriellen Design mit einer Kontrollgruppe ( ) und 5 Behandlungsgruppen ( μ 1 - 5 ). Fülle in Behandlung 1 war identisch mit der Kontrolle ( μ 1 = μ C ), in Abundanzen Behandlungen 2-5 war die Hälfte der Menge in der Kontrolle ( μ 2 - 5 = 0,5 μ c ). Für die Simulationen habe ich die Stichprobengröße (3,6,9,12) und die Häufigkeit in der Kontrollgruppe (2, 4, 8, ..., 1024) variiert. Die Abundanzen wurden aus einer negativen Binomialverteilung mit festem Dispersionsparameter ( θ = $\mu_c$$\mu_{1-5}$$\mu_1 = \mu_c$$\mu_{2-5} = 0.5 \mu_c$$\theta = 3.91$). 100 Datensätze wurden unter Verwendung einer negativen binomialen GLM und einer Gaußschen GLM + -log-transformierten Daten erzeugt und analysiert.

Das Ergebnis ist wie erwartet: Der GLM hat eine höhere Leistung, insbesondere wenn nicht viele Tiere beprobt wurden. Code ist hier.

Typ I Fehler :

Als nächstes habe ich mir den ersten Fehler angesehen. Simulationen wurden wie oben, aber alle Gruppen hatten gleiche Fülle (done ).$\mu_c = \mu_{1-5}$

Die Ergebnisse sind jedoch nicht wie erwartet: Negatives binomiales GLM zeigte einen größeren Typ-I-Fehler im Vergleich zur LM + -Transformation. Wie erwartet verschwand der Unterschied mit zunehmender Stichprobengröße. Code ist hier.

Frage:

Warum gibt es einen erhöhten Typ-I-Fehler im Vergleich zur lm + -Transformation?

Wenn wir schlechte Daten haben (kleine Stichprobengröße, geringe Häufigkeit (viele Nullen)), sollten wir dann lm + -Transformation verwenden? Kleine Probengrößen (2-4 pro Behandlung) sind typisch für solche Experimente und können nicht einfach erhöht werden.

Obwohl der neg. Behälter. GLM kann als angemessen für diese Daten begründet werden. Lm + transformation kann uns vor Fehlern vom Typ 1 bewahren.

Dies ist ein äußerst interessantes Problem. Ich habe Ihren Code überprüft und kann keinen sofort offensichtlichen Tippfehler feststellen.

$\theta$$\theta$drop1

Die meisten Tests für lineare Modelle erfordern keine Neuberechnung des Modells unter der Nullhypothese. Dies liegt daran, dass Sie die geometrische Steigung (Score-Test) und die Breite (Wald-Test) anhand von Parameterschätzungen und geschätzter Kovarianz allein unter der Alternativhypothese berechnen können.

Da das negative Binom nicht linear ist, müssen Sie meines Erachtens ein Nullmodell anpassen.

BEARBEITEN:

Ich habe den Code bearbeitet und Folgendes erhalten:

Bearbeiteter Code hier: https://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r

Die Arbeit von O'Hara und Kotze (Methoden in Ökologie und Evolution 1: 118–122) ist kein guter Ausgangspunkt für Diskussionen. Meine größte Sorge ist die Behauptung in Punkt 4 der Zusammenfassung:

Wir stellten fest, dass die Transformationen außer schlecht abschnitten. . .. Die Quasi-Poisson- und Negativ-Binomial-Modelle ... [zeigten] wenig Voreingenommenheit.

$\lambda$$\theta$$\lambda$

$\lambda$

Der folgende R-Code veranschaulicht den Punkt:

x <- rnbinom(10000, 0.5, mu=2)  
## NB: Above, this 'mu' was our lambda. Confusing, is'nt it?
log(mean(x+1))
[1] 1.09631
log(2+1)  ## Check that this is about right
[1] 1.098612

mean(log(x+1))
[1] 0.7317908

Oder Versuche

log(mean(x+.5))
[1] 0.9135269
mean(log(x+.5))
[1] 0.3270837

Die Skala, auf der die Parameter geschätzt werden, ist sehr wichtig!

$\lambda$

Beachten Sie, dass die Standarddiagnose auf einer Protokollskala (x + c) besser funktioniert. Die Wahl von c mag nicht allzu wichtig sein; oft machen 0,5 oder 1,0 sinn. Es ist auch ein besserer Ausgangspunkt für die Untersuchung von Box-Cox-Transformationen oder der Yeo-Johnson-Variante von Box-Cox. [Yeo, I. und Johnson, R. (2000)]. Weitere Informationen finden Sie auf der Hilfeseite zu powerTransform () im Fahrzeugpaket von R. Das Gamlss-Paket von R ermöglicht es, negative Binomialtypen I (die gebräuchliche Sorte) oder II oder andere Verteilungen, die sowohl die Streuung als auch den Mittelwert modellieren, mit Leistungstransformationsverbindungen von 0 (= log, dh log link) oder mehr zu kombinieren . Passungen konvergieren möglicherweise nicht immer.

Beispiel: Todesfälle versus Grundschaden Daten beziehen sich auf Hurrikane mit Namen im Atlantik, die das US-amerikanische Festland erreichten. Daten sind aus einer aktuellen Version des DAAG-Pakets für R verfügbar (Name hurricNamed ). Die Hilfeseite für die Daten enthält Details.

Das Diagramm vergleicht eine angepasste Linie, die unter Verwendung einer robusten linearen Modellanpassung erhalten wurde, mit der Kurve, die durch Transformieren einer negativen Binomialanpassung mit logarithmischer Verknüpfung auf die logarithmische Skala (Anzahl + 1) erhalten wurde, die für die y-Achse im Diagramm verwendet wurde. (Beachten Sie, dass Sie eine logarithmische Skala (Anzahl + c) mit positivem c verwenden müssen, um die Punkte und die angepasste "Linie" aus der negativen Binomialanpassung in demselben Diagramm anzuzeigen.) Beachten Sie die große Verzerrung offensichtlich für die negative Binomialanpassung auf der logarithmischen Skala. Die robuste lineare Modellanpassung ist auf dieser Skala viel weniger voreingenommen, wenn man eine negative Binomialverteilung für die Zählungen annimmt. Eine lineare Modellanpassung wäre unter den Annahmen der klassischen Normaltheorie unverzerrt. Ich fand die Voreingenommenheit erstaunlich, als ich das erste Mal das obige Diagramm erstellte! Eine Kurve würde besser zu den Daten passen. Der Unterschied liegt jedoch innerhalb der Grenzen der üblichen Standards der statistischen Variabilität. Die robuste lineare Modellanpassung leistet bei Zählungen am unteren Ende der Skala schlechte Arbeit.

Anmerkung --- Studien mit RNA-Seq-Daten: Der Vergleich der beiden Modellstile war für die Analyse von Zähldaten aus Genexpressionsexperimenten von Interesse. Die folgende Arbeit vergleicht die Verwendung eines robusten linearen Modells, das mit log (count + 1) arbeitet, mit der Verwendung negativer Binomialanpassungen (wie im Bioconductor-Package edgeR ). Die meisten Zählungen in der RNA-Seq-Anwendung, die in erster Linie im Auge behalten wird, sind groß genug, dass geeignet gewogene log-lineare Modellanpassungen sehr gut funktionieren.

Recht, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: Präzisionsgewichte ermöglichen lineare Modellanalysewerkzeuge für RNA-Seq-Lesezahlen. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

NB auch das kürzlich erschienene Papier:

Schurch NJ, Schofield P., Gierlinski M., Cole C., Sherstnev A., Singh V., Wrobel N., Gharbi K., Simpson G. G., Owen-Hughes T., Blaxter M., Barton G. J. (2016). Wie viele biologische Replikate werden in einem RNA-seq-Experiment benötigt und welches Differential-Expressions-Tool sollten Sie verwenden? RNA http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

Es ist interessant, dass das lineare Modell unter Verwendung des Limma- Pakets (wie edgeR von der WEHI-Gruppe) im Vergleich zu Ergebnissen mit vielen Replikaten, wie die Anzahl der Replikate ist, sehr gut abschneidet (im Sinne einer geringen Abweichung ) reduziert.

R-Code für die obige Grafik:

library(latticeExtra, quietly=TRUE)
hurricNamed <- DAAG::hurricNamed
ytxt <- c(0, 1, 3, 10, 30, 100, 300, 1000)
xtxt <- c(1,10, 100, 1000, 10000, 100000, 1000000 )
funy <- function(y)log(y+1)
gph <- xyplot(funy(deaths) ~ log(BaseDam2014), groups= mf, data=hurricNamed,
   scales=list(y=list(at=funy(ytxt), labels=paste(ytxt)),
           x=list(at=log(xtxt), labels=paste(xtxt))),
   xlab = "Base Damage (millions of 2014 US$); log transformed scale",
   ylab="Deaths; log transformed; offset=1",
   auto.key=list(columns=2),
   par.settings=simpleTheme(col=c("red","blue"), pch=16))
gph2 <- gph + layer(panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Name"], pos=3,
           col="gray30", cex=0.8),
        panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Year"], pos=1, 
           col="gray30", cex=0.8))
ab <- coef(MASS::rlm(funy(deaths) ~ log(BaseDam2014), data=hurricNamed))

gph3 <- gph2+layer(panel.abline(ab[1], b=ab[2], col="gray30", alpha=0.4))
## 100 points that are evenly spread on a log(BaseDam2014) scale
x <- with(hurricNamed, pretty(log(BaseDam2014),100))
df <- data.frame(BaseDam2014=exp(x[x>0])) 
hurr.nb <- MASS::glm.nb(deaths~log(BaseDam2014), data=hurricNamed[-c(13,84),])
df[,'hatnb'] <- funy(predict(hurr.nb, newdata=df, type='response'))
gph3 + latticeExtra::layer(data=df,
       panel.lines(log(BaseDam2014), hatnb, lwd=2, lty=2, 
           alpha=0.5, col="gray30"))    

Der Originalbeitrag spiegelt Tony Ives 'Artikel: Ives (2015) wider . Es ist klar, dass Signifikanztests unterschiedliche Ergebnisse für die Parameterschätzung liefern.

John Maindonald erklärt, warum die Schätzungen voreingenommen sind, aber seine Unkenntnis des Hintergrunds ist ärgerlich - er kritisiert uns dafür, dass eine Methode, der wir alle zustimmen, fehlerhaft ist, fehlerhaft ist. Viele Ökologen führen blindlings Transformationen durch, und wir haben versucht, auf die Probleme hinzuweisen, die damit verbunden sind.

Hier gibt es eine differenziertere Diskussion: Warton (2016)

Ives, AR (2015), Um die Signifikanz von Regressionskoeffizienten zu testen, gehen Sie vor und führen Sie eine Protokolltransformation durch. Methods Ecol Evol, 6: 828–835. doi: 10.1111 / 2041-210X.12386

Warton, DI, Lyons, M., Stoklosa, J. und Ives, AR (2016), Drei Punkte, die bei der Auswahl eines LM- oder GLM-Tests für die Zähldaten zu berücksichtigen sind. Methods Ecol Evol. doi: 10.1111 / 2041-210X.12552