Trudności ze znalezieniem odpowiedniego modelu dopasowanego do danych zliczających z mieszanymi efektami - ZINB czy coś innego?

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

Ten post ma cztery lata, ale chciałem śledzić zdanie fsociety w komentarzu. Diagnoza reszt w GLMM nie jest prosta, ponieważ standardowe wykresy resztkowe mogą wykazywać nienormalność, heteroscedastyczność itp., Nawet jeśli model jest poprawnie określony. Istnieje pakiet R DHARMa, specjalnie przystosowany do diagnozowania tego typu modeli.

Pakiet opiera się na podejściu symulacyjnym do generowania skalowanych resztek z dopasowanych uogólnionych liniowych modeli mieszanych i generuje różne łatwo interpretowalne wykresy diagnostyczne. Oto mały przykład z danymi z oryginalnego postu i pierwszego dopasowanego modelu (m1):

library(DHARMa)
sim_residuals <- simulateResiduals(m1, 1000)
plotSimulatedResiduals(sim_residuals)

Wykres po lewej stronie pokazuje wykres QQ skalowanych reszt w celu wykrycia odchyleń od oczekiwanego rozkładu, a wykres po prawej stronie przedstawia reszty w porównaniu z przewidywanymi wartościami podczas wykonywania regresji kwantowej w celu wykrycia odchyleń od jednorodności (czerwone linie powinny być poziome i przy 0,25 , 0,50 i 0,75).

Dodatkowo pakiet ma również określone funkcje do testowania dyspersji nadmiernej / niedostatecznej i zerowej inflacji, między innymi:

testOverdispersionParametric(m1)

Chisq test for overdispersion in GLMMs

data:  poisson
dispersion = 0.18926, pearSS = 11.35600, rdf = 60.00000, p-value = 1
alternative hypothesis: true dispersion greater 1

testZeroInflation(sim_residuals)

DHARMa zero-inflation test via comparison to expected zeros with 
simulation under H0 = fitted model


data:  sim_residuals
ratioObsExp = 0.98894, p-value = 0.502
alternative hypothesis: more