Auflösung der Heteroskedastizität in Poisson GLMM

8

Ich habe Langzeitsammeldaten und möchte testen, ob die Anzahl der gesammelten Tiere von Wettereffekten beeinflusst wird. Mein Modell sieht wie folgt aus:

glmer(SumOfCatch ~ I(pc.act.1^2) +I(pc.act.2^2) + I(pc.may.1^2) + I(pc.may.2^2) + 
                   SampSize + as.factor(samp.prog) + (1|year/month), 
      control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=1e9,npt=5)), 
      family="poisson", data=a2)

Erklärung der verwendeten Variablen:

  • SumOfCatch: Anzahl der gesammelten Tiere
  • pc.act.1, pc.act.2: Achsen einer Hauptkomponente, die die Wetterbedingungen während der Probenahme darstellen
  • pc.may.1, pc.may.2: Achsen eines PCs, die die Wetterbedingungen im Mai darstellen
  • SampSize: Anzahl der Fallstricke oder Sammeln von Durchschnitten mit Standardlängen
  • samp.prog: Probenahmemethode
  • Jahr: Jahr der Probenahme (von 1993 bis 2002)
  • Monat: Monat der Probenahme (von August bis November)

Die Residuen des angepassten Modells zeigen eine erhebliche Inhomogenität (Heteroskedastizität?), Wenn sie gegen angepasste Werte aufgetragen werden (siehe Abb. 1):

Residuen vs. Anpassungswerte - Hauptmodell

Meine Hauptfrage lautet: Ist dies ein Problem, das die Zuverlässigkeit meines Modells in Frage stellt? Wenn ja, was kann ich tun, um das Problem zu beheben?

Bisher habe ich Folgendes ausprobiert:

  • Kontrolle der Überdispersion durch Definieren zufälliger Effekte auf Beobachtungsebene, dh Verwenden einer eindeutigen ID für jede Beobachtung und Anwenden dieser ID-Variablen als zufälliger Effekt; Obwohl meine Daten eine erhebliche Überdispersion zeigen, half dies nicht, da die Residuen noch hässlicher wurden (siehe Abb. 2).

Residuen vs. angepasste Werte - Modell mit OD-Steuerung

  • Ich habe Modelle ohne zufällige Effekte mit Quasi-Poisson glm und glm.nb angepasst. ergab auch ähnliche Residuen gegenüber angepassten Parzellen wie das ursprüngliche Modell

Soweit ich weiß, gibt es möglicherweise Möglichkeiten zur Schätzung heteroskedastisch konsistenter Standardfehler, aber ich habe keine solche Methode für Poisson (oder eine andere Art von) GLMM in R gefunden.


Als Antwort auf @FlorianHartig: Die Anzahl der Beobachtungen in meinem Datensatz beträgt N = 554, ich denke, das ist eine faire Beobachtung. Nummer für ein solches Modell, aber je mehr desto besser. Ich poste zwei Figuren, von denen die erste die DHARMa-skalierte Restkurve (von Florian vorgeschlagen) des Hauptmodells ist.

Geben Sie hier die Bildbeschreibung ein

Die zweite Figur stammt aus einem zweiten Modell, bei dem der einzige Unterschied darin besteht, dass es den zufälligen Effekt auf Beobachtungsebene enthält (die erste nicht).

Geben Sie hier die Bildbeschreibung ein

AKTUALISIEREN

Abbildung der Beziehung zwischen einer Wettervariablen (als Prädiktor, dh x-Achse) und dem Stichprobenerfolg (Antwort):

Wetter-PC und Probenahme erfolgreich

UPDATE II.

Zahlen, die Prädiktorwerte gegenüber Residuen zeigen:

Prädiktoren vs. Residuen

Z. Radai
quelle
Haben Sie darüber nachgedacht, einen nichtparametrischen Schätzer auszuführen? Oder einen Ols mit einer mittleren Regression vergleichen? Mir ist klar, dass Poisson das dominierende Modell in der Bio ist, aber GLMs sind unter Heteroskedastizität inkonsistent und OLS nicht.
Superpronker
1
Manchmal wird eine Überdispersion durch eine Nullinflation verursacht. In diesem Fall könnten Sie ein Poisson-Modell mit Null-Inflations-Parameter oder ein Hürdenmodell ausprobieren. Das glmmADMB-Paket bietet großartige Funktionen, um damit umzugehen: glmmadmb.r-forge.r-project.org/glmmADMB.html
Niek
Sehr geehrter @ Superpronker, vielen Dank für den Vorschlag. Ich habe OLS nicht überprüft. Ich wusste nicht, dass dieser Ansatz flexibel genug ist, um meine Daten zu verarbeiten. Ich werde es untersuchen
Z. Radai
Sehr geehrter @Niek, in meinen Daten gibt es keine Beobachtungen von Null - ansonsten habe ich aufgrund ihres guten Umgangs mit Überdispersion an Zeroinfl- und Hürdenmodelle (im Paket 'pscl') gedacht, aber sie können nur für Daten mit Nullen in der Antwort verwendet werden . Vor ein paar Monaten habe ich glmmADMB ausprobiert, aber es hat keine besseren Ergebnisse gebracht. Prost, ZR
Z. Radai
1
@mdewey das Grundprinzip dahinter ist, dass die Beziehung zwischen Wettereffekten und Stichprobenerfolg einem Optimum folgt: Die Wahrscheinlichkeit und das Ausmaß des Stichprobenerfolgs sind in einem bestimmten Bereich (in diesem Fall Null und seine Umgebung) des Prädiktors am höchsten. Wenn der Wert des Prädiktors weiter von diesem Optimum entfernt ist, ist der Abtasterfolg geringer und entspricht einem Suboptimum. Ich verwende den quadratischen Term, weil (1) Prädiktoren bei Null neu skaliert und neu zentriert werden und (2) dies eine bessere Annäherung an eine lineare Verbindung ergibt. Prost, ZR
Z. Radai

Antworten:

9

Es ist schwierig, die Anpassung des Poisson (oder eines anderen GLM mit ganzzahligem Wert) an Pearson- oder Abweichungsreste zu bewerten, da auch ein perfekt passendes Poisson-GLMM inhomogene Abweichungsreste aufweist.

Dies ist insbesondere dann der Fall, wenn Sie GLMMs mit REs auf Beobachtungsebene durchführen, da die durch OL-REs erzeugte Dispersion von den Pearson-Residuen nicht berücksichtigt wird.

Um das Problem zu demonstrieren, erstellt der folgende Code überdisperse Poisson-Daten, die dann mit einem perfekten Modell ausgestattet werden. Die Pearson-Residuen sehen Ihrem Plot sehr ähnlich - daher kann es sein, dass es überhaupt kein Problem gibt.

Dieses Problem wird durch das DHARMa R-Paket gelöst , das aus dem angepassten Modell simuliert, um die Residuen eines GL (M) M in einen standardisierten Raum umzuwandeln. Sobald dies erledigt ist, können Sie Restprobleme wie Abweichungen von der Verteilung, Restabhängigkeit von einem Prädiktor, Heteroskedastizität oder Autokorrelation auf normale Weise visuell bewerten / testen. In der Paketvignette finden Sie Beispiele. Sie können im unteren Diagramm sehen, dass dasselbe Modell jetzt gut aussieht, wie es sollte.

Wenn Sie nach dem Plotten mit DHARMa immer noch Heteroskedastizität feststellen, müssen Sie die Dispersion als Funktion von etwas modellieren, was kein großes Problem darstellt, aber wahrscheinlich erfordern würde, dass Sie zu JAGs oder einer anderen Bayes'schen Software wechseln.

library(DHARMa)
library(lme4)

testData = createData(sampleSize = 200, overdispersion = 1, randomEffectVariance = 1, family = poisson())

fittedModel <- glmer(observedResponse ~ Environment1 + (1|group) + (1|ID), family = "poisson", data = testData, control=glmerControl(optCtrl=list(maxfun=20000) ))

# standard Pearson residuals
plot(fittedModel, resid(., type = "pearson") ~ fitted(.) , abline = 0)

# DHARMa residuals
plot(simulateResiduals(fittedModel))

Geben Sie hier die Bildbeschreibung ein

Geben Sie hier die Bildbeschreibung ein

Florian Hartig
quelle
Lieber @FlorianHartig! Vielen Dank für Ihre Erkenntnisse. Ich habe versucht, mit DHARMa zu planen. Basierend auf der Darstellung gibt es immer noch etwas, was dazu führt, dass das untere Quantil eher wie eine reziproke Kurve als wie eine gerade Linie geformt wird. Sie haben erwähnt, dass in diesem Fall eine Lösung darin bestehen könnte, die Dispersion als Funktion von etwas zu modellieren. Können Sie genau helfen, wie ich eine solche Funktion bewerten kann? Prost, ZR
Z. Radai
Kannst du mir die Handlung schicken oder posten? Eine geringfügige Variabilität wird erwartet, wenn Ihr N klein ist
Florian Hartig
Lieber @FlorianHartig, die Frage wurde bearbeitet und zeigt nun auch die DHARMa-Diagramme!
Z. Radai
@ Z.Radai - Was ich in den Darstellungen sehe, ist, dass Ihre Residuen systematisch zu hoch für niedrige Modellvorhersagen sind. Dies sieht eher nach einem Problem der Modellstruktur (fehlende Prädiktoren?) Als nach einem Verteilungsproblem aus. Ich würde versuchen, die Residuen gegen mögliche und möglicherweise fehlende Prädiktoren zu zeichnen.
Florian Hartig
1
Ich würde mir keine Sorgen um die Heteroskedastizität machen, in Ihrem Fall ist sie moderat und die Auswirkung auf die Inferenz sollte gering sein - das einzige Problem, das ich sehe, ist die systematische Unterprognose für kleine Werte, die nicht durch Modellierung der Varianz gelöst werden kann. Aber wenn Sie wissen müssen, sehen Sie hier stats.stackexchange.com/questions/247183/…
Florian Hartig