Ich habe einen sehr kleinen Datensatz zur Bienenfülle, den ich nur schwer analysieren kann. Es handelt sich um Zählungsdaten, und fast alle Zählungen erfolgen in einer Behandlung, wobei die meisten Nullen in der anderen Behandlung vorkommen. Es gibt auch einige sehr hohe Werte (jeweils einer an zwei der sechs Stellen), sodass die Verteilung der Zählungen einen extrem langen Schwanz hat. Ich arbeite in R. Ich habe zwei verschiedene Pakete verwendet: lme4 und glmmADMB.
Poisson-Mischmodelle passten nicht: Modelle waren sehr überstreut, wenn keine Zufallseffekte angepasst wurden (glm-Modell), und unterstreut, wenn Zufallseffekte angepasst wurden (glmer-Modell). Ich verstehe nicht, warum das so ist. Das experimentelle Design erfordert verschachtelte zufällige Effekte, daher muss ich sie einbeziehen. Eine Poisson-Lognormal-Fehlerverteilung verbesserte die Anpassung nicht. Ich habe mit glmer.nb versucht, eine negative binomiale Fehlerverteilung zu erzielen, und konnte sie nicht anpassen. Die Iterationsgrenze wurde erreicht, selbst wenn die Toleranz mit glmerControl geändert wurde (tolPwrss = 1e-3).
Da viele der Nullen auf die Tatsache zurückzuführen sind, dass ich die Bienen einfach nicht gesehen habe (es sind oft winzige schwarze Dinge), habe ich als nächstes ein Modell ohne Luftdruck ausprobiert. Die ZIP passte nicht gut. Der ZINB war bisher die beste Modellanpassung, aber ich bin immer noch nicht sehr zufrieden mit der Modellanpassung. Ich weiß nicht, was ich als nächstes versuchen soll. Ich habe ein Hürdenmodell ausprobiert, konnte aber keine abgeschnittene Verteilung auf die Nicht-Null-Ergebnisse anwenden - ich glaube, weil sich so viele der Nullen in der Kontrollbehandlung befinden (die Fehlermeldung lautete „Fehler in model.frame.default (formula = s.bee ~ tmt + lu +: variable längen unterscheiden sich (gefunden für 'behandlung') ”).
Darüber hinaus denke ich, dass die von mir eingeschlossene Interaktion meine Daten merkwürdig verändert, da die Koeffizienten unrealistisch klein sind - obwohl das Modell, das die Interaktion enthält, das Beste war, als ich Modelle mit AICctab im Paket bbmle verglichen habe.
Ich füge ein R-Skript hinzu, das meinen Datensatz ziemlich genau reproduziert. Variablen sind wie folgt:
d = julianisches Datum, df = julianisches Datum (als Faktor), d.sq = df im Quadrat (Anzahl der Bienen steigt und fällt im Sommer), st = Standort, s.bee = Anzahl der Bienen, tmt = Behandlung, lu = art der landnutzung, hab = prozentualer anteil des naturnahen lebensraums in der umliegenden landschaft, ba = grenzfläche rund um die felder.
Vorschläge, wie ich eine gute Modellanpassung erzielen kann (alternative Fehlerverteilungen, verschiedene Modelltypen usw.), werden sehr dankbar entgegengenommen!
Vielen Dank.
d <- c(80, 80, 121, 121, 180, 180, 86, 86, 116, 116, 144, 144, 74, 74, 143, 143, 163, 163, 71, 71,106, 106, 135, 135, 162, 162, 185, 185, 83, 83, 111, 111, 133, 133, 175, 175, 85, 85, 112, 112,137, 137, 168, 168, 186, 186, 64, 64, 95, 95, 127, 127, 156, 156, 175, 175, 91, 91, 119, 119,120, 120, 148, 148, 56, 56)
df <- as.factor(d)
d.sq <- d^2
st <- factor(rep(c("A", "B", "C", "D", "E", "F"), c(6,12,18,10,14,6)))
s.bee <- c(1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,4,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,3,0,0,0,0,5,0,0,2,0,50,0,10,0,4,0,47,3)
tmt <- factor(c("AF","C","C","AF","AF","C","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C",
"C","AF","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C","AF","C","AF","C",
"C","AF","AF","C","AF","C","AF","C","AF","C","C","AF","C","AF","C","AF","AF","C","AF","C",
"AF","C","AF","C","AF","C"))
lu <- factor(rep(c("p","a","p","a","p"), c(6,12,28,14,6)))
hab <- rep(c(13,14,13,14,3,4,3,4,3,4,3,4,3,4,15,35,37,35,37,35,37,35,37,0,2,1,2,1,2,1),
c(1,2,2,1,1,1,1,2,2,1,1,1,1,1,18,1,1,1,2,2,1,1,1,14,1,1,1,1,1,1))
ba <- c(480,6520,6520,480,480,6520,855,1603,855,1603,1603,855,855,12526,855,5100,855,5100,2670,7679,7679,2670,
2670,7679,2670,7679,7679,2670,2670,7679,2670,7679,2670,7679,2670,7679,1595,3000,1595,3000,3000,1595,1595,3000,1595
,3000,4860,5460,4860,5460,5460,4860,5460,4860,5460,4860,4840,5460,4840,5460,3000,1410,3000,1410,3000,1410)
data <- data.frame(st,df,d.sq,tmt,lu,hab,ba,s.bee)
with(data, table(s.bee, tmt) )
# below is a much abbreviated summary of attempted models:
library(MASS)
library(lme4)
library(glmmADMB)
library(coefplot2)
###
### POISSON MIXED MODEL
m1 <- glmer(s.bee ~ tmt + lu + hab + (1|st/df), family=poisson)
summary(m1)
resdev<-sum(resid(m1)^2)
mdf<-length(fixef(m1))
rdf<-nrow(data)-mdf
resdev/rdf
# 0.2439303
# underdispersed. ???
###
### NEGATIVE BINOMIAL MIXED MODEL
m2 <- glmer.nb(s.bee ~ tmt + lu + hab + d.sq + (1|st/df))
# iteration limit reached. Can't make a model work.
###
### ZERO-INFLATED POISSON MIXED MODEL
fit_zipoiss <- glmmadmb(s.bee~tmt + lu + hab + ba + d.sq +
tmt:lu +
(1|st/df), data=data,
zeroInflation=TRUE,
family="poisson")
# has to have lots of variables to fit
# anyway Poisson is not a good fit
###
### ZERO-INFLATED NEGATIVE BINOMIAL MIXED MODELS
## BEST FITTING MODEL SO FAR:
fit_zinb <- glmmadmb(s.bee~tmt + lu + hab +
tmt:lu +
(1|st/df),data=data,
zeroInflation=TRUE,
family="nbinom")
summary(fit_zinb)
# coefficients are tiny, something odd going on with the interaction term
# but this was best model in AICctab comparison
# model check plots
qqnorm(resid(fit_zinb))
qqline(resid(fit_zinb))
coefplot2(fit_zinb)
resid_zinb <- resid(fit_zinb , type = "pearson")
hist(resid_zinb)
fitted_zinb <- fitted (fit_zinb)
plot(resid_zinb ~ fitted_zinb)
## MODEL WITHOUT INTERACTION TERM - the coefficients are more realistic:
fit_zinb2 <- glmmadmb(s.bee~tmt + lu + hab +
(1|st/df),data=data,
zeroInflation=TRUE,
family="nbinom")
# model check plots
qqnorm(resid(fit_zinb2))
qqline(resid(fit_zinb2))
coefplot2(fit_zinb2)
resid_zinb2 <- resid(fit_zinb2 , type = "pearson")
hist(resid_zinb2)
fitted_zinb2 <- fitted (fit_zinb2)
plot(resid_zinb2 ~ fitted_zinb2)
# ZINB models are best so far
# but I'm not happy with the model check plots
Antworten:
Dieser Beitrag hat vier Jahre, aber ich wollte dem folgen, was die Gesellschaft in einem Kommentar gesagt hat. Die Diagnose von Residuen in GLMMs ist nicht einfach, da Standard-Residuendiagramme auch bei korrekt angegebenem Modell Nicht-Normalität, Heteroskedastizität usw. aufweisen können. Es gibt ein R-Paket, das
DHARMaspeziell für die Diagnose solcher Modelle geeignet ist.Das Paket basiert auf einem Simulationsansatz zur Generierung skalierter Residuen aus angepassten verallgemeinerten linearen gemischten Modellen und generiert verschiedene leicht interpretierbare diagnostische Diagramme. Hier ein kleines Beispiel mit den Daten vom Originalpfosten und dem ersten angepassten Modell (m1):
Das Diagramm auf der linken Seite zeigt ein QQ-Diagramm der skalierten Residuen, um Abweichungen von der erwarteten Verteilung zu erkennen, und das Diagramm auf der rechten Seite stellt Residuen gegenüber den vorhergesagten Werten dar, während eine Quantilregression durchgeführt wird, um Abweichungen von der Homogenität zu erkennen (rote Linien sollten horizontal und bei 0,25 liegen 0,50 und 0,75).
Darüber hinaus verfügt das Paket über spezifische Funktionen zum Testen auf Über- / Unterdispersion und Null-Inflation, unter anderem:
quelle