Mam bardzo mały zestaw danych na temat liczebności pojedynczych pszczół, które mam problemy z analizą. Są to dane zliczania i prawie wszystkie zliczenia są w jednym traktowaniu, a większość zer w drugim traktowaniu. Istnieje również kilka bardzo wysokich wartości (po jednej w dwóch z sześciu miejsc), więc rozkład zliczeń ma wyjątkowo długi ogon. Pracuję w R. Użyłem dwóch różnych pakietów: lme4 i glmmADMB.
Modele mieszane Poissona nie pasowały: modele były bardzo rozproszone, gdy nie dopasowano efektów losowych (model glm), i rozproszone, gdy dopasowano efekty losowe (model glm). Nie rozumiem, dlaczego tak jest. Projekt eksperymentalny wymaga zagnieżdżonych efektów losowych, dlatego muszę je uwzględnić. Logarytmiczny rozkład błędów Poissona nie poprawił dopasowania. Próbowałem ujemnego rozkładu błędów dwumianowych za pomocą glmer.nb i nie mogłem go dopasować - osiągnięto limit iteracji, nawet po zmianie tolerancji za pomocą glmerControl (tolPwrss = 1e-3).
Ponieważ wiele zer wynika z faktu, że po prostu nie widziałem pszczół (często są to małe czarne rzeczy), następnie wypróbowałem model z napompowaniem zerowym. ZIP nie pasował dobrze. ZINB był jak dotąd najlepszym dopasowaniem do modelu, ale nadal nie jestem zbyt zadowolony z dopasowania modelu. Nie wiem, co dalej. Wypróbowałem model przeszkodowy, ale nie mogłem dopasować skróconego rozkładu do niezerowych wyników - myślę, ponieważ tak wiele zer jest w traktowaniu kontrolnym (komunikat o błędzie brzmiał: „Błąd w modelu.frame.default (formuła = s.bee ~ tmt + lu +: różne długości różnią się (znaleziono dla „leczenia”) ”).
Ponadto myślę, że interakcja, którą podałem, robi coś dziwnego z moimi danymi, ponieważ współczynniki są nierealistycznie małe - chociaż model zawierający interakcję był najlepszy, gdy porównałem modele wykorzystujące AICctab w pakiecie bbmle.
Dołączam jakiś skrypt R, który prawie odtworzy mój zestaw danych. Zmienne są następujące:
d = data Juliana, df = data Juliana (jako czynnik), d.sq = df do kwadratu (liczba pszczół rośnie, a następnie spada przez całe lato), st = miejsce, s.bee = liczba pszczół, tmt = leczenie, lu = rodzaj użytkowania gruntów, hab = procent siedlisk półnaturalnych w otaczającym krajobrazie, ba = obszar graniczny okrągłych pól.
Wszelkie sugestie dotyczące tego, w jaki sposób mogę uzyskać dobre dopasowanie modelu (alternatywne rozkłady błędów, różne typy modeli itp.) Zostaną bardzo wdzięczne!
Dziękuję Ci.
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
