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

18

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 ( θ = 3,91) gezogenμcμ15μ1=μcμ25=0.5μcθ=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. Bildbeschreibung hier eingeben 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 ).μc=μ15

Die Ergebnisse sind jedoch nicht wie erwartet: Bildbeschreibung hier eingeben 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.

EDi
quelle
1
Keine Antwort auf Ihre Hauptfrage, aber etwas, das der Leser beachten sollte: Wenn Sie nicht die tatsächliche Fehlerart I für die beiden Verfahren als äquivalent festlegen, ist ein Vergleich der Leistung nicht sinnvoll. Ich kann die Leistung immer für die niedrigere erhöhen (in diesem Fall für die Take Logs und für die normale), indem ich den Fehler vom Typ I behebe. Wenn Sie andererseits die jeweilige Situation (Stichprobengröße, Häufigkeit) angeben, können Sie die Fehlerrate des Typs I ermitteln (z. B. durch Simulation) und so die zu testende Nennrate ermitteln, um die gewünschte Fehlerrate des Typs I zu erzielen , so wird ihre Leistung vergleichbar.
Glen_b -Reinstate Monica
Werden die y-Achsenwerte in Ihren Plots über die 100 Datensätze gemittelt?
Shadowtalker
Ich sollte meinen Kommentar klarstellen: In Fällen, in denen die Statistiken von Natur aus diskret sind, haben Sie keine perfekte Kontrolle über die Fehlerrate von Typ I, aber Sie können die Fehlerraten von Typ I im Allgemeinen ziemlich genau festlegen. In Situationen, in denen sie nicht nahe genug beieinander liegen, um vergleichbar zu sein, sind randomisierte Tests die einzige Möglichkeit, sie vergleichbar zu machen.
Glen_b
α
1
n

Antworten:

17

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

θθ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: Bildbeschreibung hier eingeben

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

AdamO
quelle
Aber ich denke, das drop1() passt intern wieder zum Nullmodell ...
Ben Bolker
4
glm.nbθdrop1logLikgetS3method('logLik', 'negbin'
möchte nochmal +1 geben aber ich kann nicht. Nett.
Ben Bolker
Vielen Dank! Ich habe mir nur den Code von beiden drop1()und angesehen lrtest(). Du hast recht, drop1.glmbenutzt glm.fitwas die falsche Abweichung ergibt. War nicht bewusst , dass wir nicht verwenden können drop1()mit glm.nb()!
EDi
Der typische Score und die Wald-Tests sind also im negativen Binomialmodell ungültig?
Shadowtalker
8

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.

λθλ

λ

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!

λ

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.

Robust loglinear gegen negative Binomialpassung

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"))    
John Maindonald
quelle
2
Vielen Dank für Ihren Kommentar, Herr Maindonald. In den letzten zwei Jahren gab es auch einige weitere Artikel (die sich mehr auf Hypothesentests als auf Verzerrungen konzentrierten): Ives 2015, Warton et al. 2016, Szöcs 2015.
EDi
Vielleicht ist es ein guter Ausgangspunkt für Diskussionen, auch wenn dieser Punkt problematisch ist? (Ich würde allgemeiner argumentieren, dass dies ein Grund ist, sich nicht zu sehr auf Voreingenommenheit zu konzentrieren, sondern etwas wie RMSE in Betracht zu ziehen ... [Haftungsausschluss, ich habe diese Artikel in letzter Zeit nicht erneut gelesen und ich habe nur das Abstract von gelesen die Warton-Zeitung ...]
Ben Bolker
1
Der Punkt von Warton et al. (2016), dass Dateneigenschaften die Entscheidungsgründe sein sollten, ist entscheidend. Quantil-Quantil-Diagramme sind eine gute Möglichkeit, die Details der Anpassungen zu vergleichen. Insbesondere die Anpassung an das eine oder andere oder an beide Extreme kann für einige Anwendungen wichtig sein. Modelle mit Null-Inflation oder Hürden können eine effektive Verfeinerung sein, um die richtigen Nullzählungen zu erzielen. Im oberen Extremfall kann jedes der diskutierten Modelle stark beeinträchtigt werden. Lobenswerterweise haben Warton et al ein Beispiel. Ich würde gerne Vergleiche mit einer Vielzahl von ökologischen Datensätzen sehen.
John Maindonald
Aber sind in ökologischen Datensätzen die Arten im unteren Teil (= seltene Arten) nicht interessant? Sollte nicht zu schwierig sein, einige ökologische Datensätze zusammenzustellen und zu vergleichen ...
EDi
Tatsächlich scheint das negative Binomialmodell für das untere Ende der Schadenskategorie für die Hurrikantodesdaten am wenigsten zufriedenstellend zu sein. Das Gamlss- Paket von R hat eine Funktion, die es einfach macht, Zentile der angepassten Verteilung mit Zentilen der Daten zu vergleichen:
John Maindonald,
6

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

Bob O'Hara
quelle
Willkommen zum Lebenslauf. Obwohl diese Antwort hilfreich ist, handelt es sich meistens um eine Antwort vom Typ "Nur-Link". Links ändern sich und werden entfernt. Es wäre für den Lebenslauf hilfreicher, wenn Sie die wichtigsten Punkte in jedem einzelnen erläutern würden.
Mike Hunter
Danke für die Antwort. Ich denke, die Arbeit von Warton et al. prägt den aktuellen Diskussionsstand.
EDi
Danke willkommen! Ich habe mir erlaubt, die Referenzen vollständig hinzuzufügen.
Scortchi
1
Bitte skizzieren Sie die wichtigsten Punkte, die in den neuen Referenzen gemacht werden, und beziehen Sie sich, wo es sinnvoll ist, auch auf die ursprüngliche Frage. Dies ist ein wertvoller Beitrag, der derzeit eher einem Kommentar zu einer anderen Antwort als einer Antwort auf die Frage (die beispielsweise einen Kontext für Links bieten sollte ) entspricht. Ein paar zusätzliche Kontextsätze würden dem Beitrag erheblich helfen.
Glen_b
3
Meine Kommentare beziehen sich insbesondere auf Punkt 4 in der Veröffentlichung von O'Hara und Kotze: "Wir stellten fest, dass die Transformationen schlecht abschnitten, außer ... Die Quasi-Poisson- und negativen Binomialmodelle ... [zeigten] wenig Voreingenommenheit." Die Simulationen sind ein Kommentar zum Vergleich zwischen dem erwarteten Mittelwert auf einer Skala von y (die Zählungen) und dem erwarteten Mittelwert auf einer Skala von log (y + c) für eine äußerst positive Versatzverteilung, mehr nicht. Der negative Binomialparameter Lambda ist auf der Skala von y unverzerrt, während der logarithmisch-normale Mittelwert auf einer Skala von logarithmisch (y + c) unverzerrt ist (unter Normalität auf dieser Skala).
John Maindonald