Ich habe einige Fragen zur Spezifikation und Interpretation von GLMMs. 3 Fragen sind definitiv statistisch und 2 beziehen sich genauer auf R. Ich poste hier, weil ich letztendlich denke, dass das Problem die Interpretation der GLMM-Ergebnisse ist.
Ich versuche gerade, ein GLMM zu montieren. Ich verwende US-Volkszählungsdaten aus der Longitudinal Tract Database . Meine Beobachtungen sind Zensusdaten. Meine abhängige Variable ist die Anzahl der freien Wohneinheiten und ich interessiere mich für die Beziehung zwischen Leerstand und sozioökonomischen Variablen. Das Beispiel hier ist einfach und verwendet nur zwei feste Effekte: Prozent der nicht weißen Bevölkerung (Rasse) und das mittlere Haushaltseinkommen (Klasse) sowie deren Interaktion. Ich möchte zwei verschachtelte zufällige Effekte einbeziehen: Traktate innerhalb von Jahrzehnten und Jahrzehnten, dh (Jahrzehnt / Traktat). Ich betrachte diese Zufälle, um die räumliche (dh zwischen Traktaten) und zeitliche (dh zwischen Jahrzehnten) Autokorrelation zu kontrollieren. Ich interessiere mich jedoch auch für das Jahrzehnt als festen Effekt, daher beziehe ich es auch als festen Faktor ein.
Da meine unabhängige Variable eine nicht negative Ganzzahl-Zählvariable ist, habe ich versucht, Poisson- und negative Binomial-GLMMs anzupassen. Ich verwende das Protokoll der gesamten Wohneinheiten als Versatz. Dies bedeutet, dass Koeffizienten als Auswirkung auf die Leerstandsquote und nicht als Gesamtzahl der leer stehenden Häuser interpretiert werden.
Ich habe derzeit Ergebnisse für ein Poisson und ein negatives binomisches GLMM, die mit glmer und glmer.nb von lme4 geschätzt wurden . Die Interpretation von Koeffizienten ist für mich aufgrund meiner Kenntnis der Daten und des Untersuchungsgebiets sinnvoll.
Wenn Sie die Daten und das Skript möchten , befinden sie sich auf meinem Github . Das Skript enthält weitere beschreibende Untersuchungen, die ich vor dem Erstellen der Modelle durchgeführt habe.
Hier sind meine Ergebnisse:
Poisson-Modell
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson ( log )
Formula: R_VAC ~ decade + P_NONWHT + a_hinc + P_NONWHT * a_hinc + offset(HU_ln) + (1 | decade/TRTID10)
Data: scaled.mydata
AIC BIC logLik deviance df.resid
34520.1 34580.6 -17250.1 34500.1 3132
Scaled residuals:
Min 1Q Median 3Q Max
-2.24211 -0.10799 -0.00722 0.06898 0.68129
Random effects:
Groups Name Variance Std.Dev.
TRTID10:decade (Intercept) 0.4635 0.6808
decade (Intercept) 0.0000 0.0000
Number of obs: 3142, groups: TRTID10:decade, 3142; decade, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.612242 0.028904 -124.98 < 2e-16 ***
decade1980 0.302868 0.040351 7.51 6.1e-14 ***
decade1990 1.088176 0.039931 27.25 < 2e-16 ***
decade2000 1.036382 0.039846 26.01 < 2e-16 ***
decade2010 1.345184 0.039485 34.07 < 2e-16 ***
P_NONWHT 0.175207 0.012982 13.50 < 2e-16 ***
a_hinc -0.235266 0.013291 -17.70 < 2e-16 ***
P_NONWHT:a_hinc 0.093417 0.009876 9.46 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) dc1980 dc1990 dc2000 dc2010 P_NONWHT a_hinc
decade1980 -0.693
decade1990 -0.727 0.501
decade2000 -0.728 0.502 0.530
decade2010 -0.714 0.511 0.517 0.518
P_NONWHT 0.016 0.007 -0.016 -0.015 0.006
a_hinc -0.023 -0.011 0.023 0.022 -0.009 0.221
P_NONWHT:_h 0.155 0.035 -0.134 -0.129 0.003 0.155 -0.233
convergence code: 0
Model failed to converge with max|grad| = 0.00181132 (tol = 0.001, component 1)
Negatives Binomialmodell
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(25181.5) ( log )
Formula: R_VAC ~ decade + P_NONWHT + a_hinc + P_NONWHT * a_hinc + offset(HU_ln) + (1 | decade/TRTID10)
Data: scaled.mydata
AIC BIC logLik deviance df.resid
34522.1 34588.7 -17250.1 34500.1 3131
Scaled residuals:
Min 1Q Median 3Q Max
-2.24213 -0.10816 -0.00724 0.06928 0.68145
Random effects:
Groups Name Variance Std.Dev.
TRTID10:decade (Intercept) 4.635e-01 6.808e-01
decade (Intercept) 1.532e-11 3.914e-06
Number of obs: 3142, groups: TRTID10:decade, 3142; decade, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.612279 0.028946 -124.79 < 2e-16 ***
decade1980 0.302897 0.040392 7.50 6.43e-14 ***
decade1990 1.088211 0.039963 27.23 < 2e-16 ***
decade2000 1.036437 0.039884 25.99 < 2e-16 ***
decade2010 1.345227 0.039518 34.04 < 2e-16 ***
P_NONWHT 0.175216 0.012985 13.49 < 2e-16 ***
a_hinc -0.235274 0.013298 -17.69 < 2e-16 ***
P_NONWHT:a_hinc 0.093417 0.009879 9.46 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) dc1980 dc1990 dc2000 dc2010 P_NONWHT a_hinc
decade1980 -0.693
decade1990 -0.728 0.501
decade2000 -0.728 0.502 0.530
decade2010 -0.715 0.512 0.517 0.518
P_NONWHT 0.016 0.007 -0.016 -0.015 0.006
a_hinc -0.023 -0.011 0.023 0.022 -0.009 0.221
P_NONWHT:_h 0.154 0.035 -0.134 -0.129 0.003 0.155 -0.233
Poisson DHARMa-Tests
One-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.044451, p-value = 8.104e-06
alternative hypothesis: two-sided
DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
data: simulationOutput
ratioObsExp = 1.3666, p-value = 0.159
alternative hypothesis: more
Negative binomiale DHARMa-Tests
One-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.04263, p-value = 2.195e-05
alternative hypothesis: two-sided
DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
data: simulationOutput2
ratioObsExp = 1.376, p-value = 0.174
alternative hypothesis: more
DHARMa-Grundstücke
Poisson
Negatives Binomial
Statistikfragen
Da ich immer noch GLMMs herausfinde, fühle ich mich hinsichtlich Spezifikation und Interpretation unsicher. Ich habe ein paar Fragen:
Es scheint, dass meine Daten die Verwendung eines Poisson-Modells nicht unterstützen, und daher bin ich mit negativem Binomial besser dran. Ich erhalte jedoch immer wieder Warnungen, dass meine negativen Binomialmodelle ihre Iterationsgrenze erreichen, selbst wenn ich die maximale Grenze erhöhe. "In theta.ml (Y, mu, Gewichte = Objekt @ bzw. $ Gewichte, Grenze = Grenze ,: Iterationsgrenze erreicht." Dies geschieht unter Verwendung einiger verschiedener Spezifikationen (dh Minimal- und Maximalmodelle für feste und zufällige Effekte). Ich habe auch versucht, Ausreißer in meinem abhängigen Bereich zu entfernen (brutto, ich weiß!), Da die oberen 1% der Werte sehr viele Ausreißer sind (untere 99% reichen von 0-1012, obere 1% von 1013-5213). Ich habe keine Auswirkungen auf Iterationen und auch nur sehr geringe Auswirkungen auf die Koeffizienten. Ich füge diese Details hier nicht hinzu. Die Koeffizienten zwischen Poisson und negativem Binom sind ebenfalls ziemlich ähnlich. Ist dieser Mangel an Konvergenz ein Problem? Passt das negative Binomialmodell gut? Ich habe auch das negative Binomialmodell mit ausgeführtAllFit und nicht alle Optimierer geben diese Warnung aus (Bobyqa, Nelder Mead und nlminbw nicht).
Die Varianz für meinen zehnjährigen festen Effekt ist durchweg sehr gering oder 0. Ich verstehe, dass dies bedeuten könnte, dass das Modell überangepasst ist. Wenn Sie ein Jahrzehnt aus den festen Effekten herausnehmen, erhöht sich die Varianz der zufälligen Effekte des Jahrzehnts auf 0,2620 und hat keinen großen Einfluss auf die Koeffizienten der festen Effekte. Ist etwas falsch daran, es zu belassen? Ich interpretiere es gut als einfach nicht nötig, um zwischen Beobachtungsvarianz zu erklären.
Zeigen diese Ergebnisse an, dass ich Modelle ohne Luftdruck ausprobieren sollte? DHARMa scheint darauf hinzudeuten, dass eine Nullinflation möglicherweise nicht das Problem ist. Wenn Sie denken, ich sollte es trotzdem versuchen, siehe unten.
R Fragen
Ich wäre bereit, Modelle mit null Inflation auszuprobieren, bin mir jedoch nicht sicher, welches Paket verschachtelte zufällige Effekte für Poisson mit null Inflation und negative binomiale GLMMs implementiert. Ich würde glmmADMB verwenden, um AIC mit null aufgeblasenen Modellen zu vergleichen, aber es ist auf einen einzelnen zufälligen Effekt beschränkt, funktioniert also nicht für dieses Modell. Ich könnte MCMCglmm ausprobieren, aber ich kenne keine Bayes'schen Statistiken, so dass dies auch nicht attraktiv ist. Irgendwelche anderen Optionen?
Kann ich potenzierte Koeffizienten innerhalb der Zusammenfassung (Modell) anzeigen oder muss ich dies außerhalb der Zusammenfassung tun, wie ich es hier getan habe?
quelle
decade
sowohl fest als auch zufällig zu haben, macht keinen Sinn. Entweder haben Sie es als fest und schließen es nur(1 | decade:TRTID10)
als zufällig ein (was der(1 | TRTID10)
Annahme entspricht, dass SieTRTID10
für verschiedene Jahrzehnte nicht die gleichen Stufen haben), oder entfernen Sie es aus den festen Effekten. Mit nur 4 Levels ist es vielleicht besser, wenn Sie es reparieren lassen: Die übliche Empfehlung ist, zufällige Effekte anzupassen, wenn man 5 Levels oder mehr hat.bobyqa
Optimierungsprogramm ausprobiert haben und es keine Warnung ausgegeben hat. Was ist dann das Problem? Einfach benutzenbobyqa
.bobyqa
konvergiert besser als das Standardoptimierungsprogramm (und ich glaube, ich habe irgendwo gelesen, dass es in zukünftigen Versionen von zum Standardoptimierer werden wirdlme4
). Ich glaube nicht, dass Sie sich über die Nichtkonvergenz mit dem Standardoptimierer Gedanken machen müssen, wenn dieser konvergiertbobyqa
.Antworten:
Ich glaube, es gibt einige wichtige Probleme, die mit Ihrer Einschätzung angegangen werden müssen.
Nach dem, was ich durch die Untersuchung Ihrer Daten gesammelt habe, sind Ihre Einheiten nicht geografisch gruppiert, dh Volkszählungsdaten innerhalb von Landkreisen. Die Verwendung von Traktaten als Gruppierungsfaktor ist daher nicht geeignet, um räumliche Heterogenität zu erfassen, da dies bedeutet, dass Sie die gleiche Anzahl von Personen wie Gruppen haben (oder anders ausgedrückt, alle Ihre Gruppen haben jeweils nur eine Beobachtung). Die Verwendung einer mehrstufigen Modellierungsstrategie ermöglicht es uns, die Varianz auf Einzelebene zu schätzen und gleichzeitig die Varianz zwischen Gruppen zu kontrollieren. Da Ihre Gruppen jeweils nur ein Individuum haben, entspricht Ihre Varianz zwischen den Gruppen Ihrer Varianz auf Einzelebene, wodurch der Zweck des Mehrebenenansatzes zunichte gemacht wird.
Andererseits kann der Gruppierungsfaktor wiederholte Messungen über die Zeit darstellen. Zum Beispiel können im Fall einer Längsschnittstudie die "Mathematik" -Ergebnisse einer Person jährlich neu bewertet werden, sodass wir für jeden Schüler einen jährlichen Wert für n Jahre haben würden (in diesem Fall ist der Gruppierungsfaktor der Schüler, wie wir n haben Anzahl der Beobachtungen, die innerhalb der Schüler "verschachtelt" sind). In Ihrem Fall haben Sie die Messungen jedes Zensus-Trakts von wiederholt
decade
. Daher können Sie IhreTRTID10
Variable als Gruppierungsfaktor verwenden, um "Varianz zwischen Jahrzehnten" zu erfassen. Dies führt zu 3142 Beobachtungen, die in 635 Traktaten verschachtelt sind, mit ungefähr 4 und 5 Beobachtungen pro Zensus-Traktat.Wie in einem Kommentar erwähnt, ist die Verwendung
decade
als Gruppierungsfaktor nicht sehr angemessen, da Sie nur etwa 5 Jahrzehnte für jeden Zensus-Trakt haben und deren Wirkung besser erfasst werden kann, wenn Sie siedecade
als Kovariate einführen .Zweitens, um zu bestimmen, ob Ihre Daten mit einem Poisson- oder einem negativen Binomialmodell (oder einem Null-Inflations-Ansatz) modelliert werden sollen. Berücksichtigen Sie das Ausmaß der Überdispersion in Ihren Daten. Das grundlegende Merkmal einer Poisson-Verteilung ist die Äquidispersion, was bedeutet, dass der Mittelwert gleich der Varianz der Verteilung ist. Wenn Sie sich Ihre Daten ansehen, ist es ziemlich klar, dass es viel Überdispersion gibt. Varianzen sind viel viel größer als die Mittel.
Um jedoch festzustellen, ob das negative Binom statistisch besser geeignet ist, besteht eine Standardmethode darin, einen Likelihood-Ratio-Test zwischen einem Poisson- und einem negativen Binomialmodell durchzuführen, was darauf hindeutet, dass das Negbin besser passt.
Nachdem dies festgestellt wurde, könnte ein nächster Test prüfen, ob der Mehrebenenansatz (gemischtes Modell) unter Verwendung eines ähnlichen Ansatzes gerechtfertigt ist, was darauf hindeutet, dass die Mehrebenenversion eine bessere Anpassung bietet. (Ein ähnlicher Test könnte verwendet werden, um eine glmer-Anpassung unter der Annahme einer Poisson-Verteilung mit dem glmer.nb-Objekt zu vergleichen, sofern die Modelle ansonsten identisch sind.)
In Bezug auf die Schätzungen der Poisson- und nb-Modelle wird erwartet, dass sie einander sehr ähnlich sind, wobei der Hauptunterschied die Standardfehler sind, dh wenn eine Überdispersion vorliegt, neigt das Poisson-Modell dazu, voreingenommene Standardfehler bereitzustellen. Nehmen Sie Ihre Daten als Beispiel:
Beachten Sie, dass die Koeffizientenschätzungen alle sehr ähnlich sind. Der Hauptunterschied besteht nur in der Signifikanz einer Ihrer Kovariaten sowie in der Varianz der Zufallseffekte, was darauf hindeutet, dass die Varianz auf Einheitsebene durch den Überdispersionsparameter in nb erfasst wird Das Modell (der
theta
Wert im Objekt glmer.nb) erfasst einen Teil der Varianz zwischen den Trakten, die durch die zufälligen Effekte erfasst wird.In Bezug auf potenzierte Koeffizienten (und zugehörige Konfidenzintervalle) können Sie Folgendes verwenden:
Letzte Gedanken zur Nullinflation. Es gibt keine mehrstufige Implementierung (zumindest ist mir bekannt) eines Poisson- oder Negbin-Modells mit Null-Inflation, mit der Sie eine Gleichung für die Null-Inflation-Komponente der Mischung angeben können. Mit dem
glmmADMB
Modell können Sie einen konstanten Null-Inflationsparameter schätzen. Ein alternativer Ansatz wäre die Verwendung derzeroinfl
Funktion impscl
Paket, obwohl dies keine Mehrebenenmodelle unterstützt. Auf diese Weise können Sie die Anpassung eines einstufigen negativen Binomials mit dem einstufigen negativen Binomial Null vergleichen. Wenn die Nullinflation für einstufige Modelle nicht signifikant ist, ist es wahrscheinlich, dass sie für die Mehrebenenspezifikation nicht signifikant ist.Nachtrag
Wenn Sie sich Sorgen über die räumliche Autokorrelation machen, können Sie dies mithilfe einer geografischen gewichteten Regression steuern (obwohl ich glaube, dass hier Punktdaten und keine Bereiche verwendet werden). Alternativ können Sie Ihre Zensusdaten nach einem zusätzlichen Gruppierungsfaktor (Bundesstaaten, Landkreise) gruppieren und diesen als zufälligen Effekt einbeziehen. Und schließlich, und ich bin nicht sicher, ob dies vollständig machbar ist, kann es möglich sein, die räumliche Abhängigkeit zu berücksichtigen, indem beispielsweise die durchschnittliche Anzahl von
R_VAC
Nachbarn erster Ordnung als Kovariate verwendet wird. In jedem Fall wäre es vor solchen Ansätzen sinnvoll festzustellen, ob tatsächlich eine räumliche Autokorrelation vorliegt (unter Verwendung von I-, LISA-Tests und ähnlichen Ansätzen von Global Moran).quelle
brms
kann null aufgeblasene negative Binomialmodelle mit zufälligen Effekten anpassen.brms
und dieses mit dem oben beschriebenen Modell glmer.nb zu vergleichen. Ich werde auch versuchen, einen von der Volkszählung definierten Ort (im Grunde Gemeinde, ~ 170 Gruppen) als Gruppierungsfaktor für die zufälligen Effekte einzubeziehen (nur 5 Landkreise in den Daten, daher werde ich das nicht verwenden). Ich werde auch die räumliche Autokorrelation von Residuen mit Global Morans I testen. Ich werde darüber berichten, wenn ich das getan habe.brms
und zu vergleichen.