J'ai un très petit ensemble de données sur l'abondance des abeilles solitaires que j'ai du mal à analyser. Ce sont des données de comptage, et presque tous les comptages sont dans un traitement avec la plupart des zéros dans l'autre traitement. Il existe également quelques valeurs très élevées (une sur deux des six sites), de sorte que la distribution des dénombrements a une queue extrêmement longue. Je travaille dans R. J'ai utilisé deux packages différents: lme4 et glmmADMB.
Les modèles mixtes de Poisson ne correspondaient pas: les modèles étaient très sur-dispersés lorsque les effets aléatoires n'étaient pas ajustés (modèle glm), et sous-dispersés lorsque les effets aléatoires étaient ajustés (modèle glmer). Je ne comprends pas pourquoi c'est. La conception expérimentale nécessite des effets aléatoires imbriqués, je dois donc les inclure. Une distribution d'erreur lognormale de Poisson n'a pas amélioré l'ajustement. J'ai essayé la distribution d'erreur binomiale négative en utilisant glmer.nb et je n'ai pas réussi à l'adapter - la limite d'itération a été atteinte, même lorsque j'ai modifié la tolérance en utilisant glmerControl (tolPwrss = 1e-3).
Parce que beaucoup de zéros seront dus au fait que je n'ai tout simplement pas vu les abeilles (ce sont souvent de minuscules choses noires), j'ai ensuite essayé un modèle gonflé à zéro. Le ZIP ne convenait pas bien. Le ZINB était le meilleur modèle à ce jour, mais je ne suis toujours pas trop satisfait du modèle. Je ne sais pas quoi essayer ensuite. J'ai essayé un modèle d'obstacle mais je n'ai pas pu ajuster une distribution tronquée aux résultats non nuls - je pense parce que beaucoup de zéros sont dans le traitement de contrôle (le message d'erreur était «Erreur dans model.frame.default (formule = s.bee ~ tmt + lu +: les longueurs variables diffèrent (trouvées pour «traitement») »).
De plus, je pense que l'interaction que j'ai incluse fait quelque chose d'étrange pour mes données car les coefficients sont irréalistes - bien que le modèle contenant l'interaction était meilleur lorsque j'ai comparé des modèles utilisant AICctab dans le paquet bbmle.
J'inclus un script R qui reproduira à peu près mon ensemble de données. Les variables sont les suivantes:
d = date julienne, df = date julienne (comme facteur), d.sq = df au carré (le nombre d'abeilles augmente puis diminue tout au long de l'été), st = site, s.bee = nombre d'abeilles, tmt = traitement, lu = type d'utilisation des terres, hab = pourcentage d'habitat semi-naturel dans le paysage environnant, ba = zone périphérique autour des champs.
Toutes les suggestions sur la façon d'obtenir un bon ajustement du modèle (distributions d'erreur alternatives, différents types de modèle, etc.) seraient très appréciées!
Je vous remercie.
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
Réponses:
Ce poste a quatre ans, mais je voulais suivre ce que la société a dit dans un commentaire. Le diagnostic des résidus dans les GLMM n'est pas simple, car les graphiques résiduels standard peuvent montrer la non-normalité, l'hétéroscédasticité, etc., même si le modèle est correctement spécifié. Il existe un package R
DHARMa, spécialement adapté au diagnostic de ce type de modèles.Le package est basé sur une approche de simulation pour générer des résidus à l'échelle à partir de modèles mixtes linéaires généralisés ajustés et génère différents graphiques de diagnostic facilement interprétables. Voici un petit exemple avec les données du poste d'origine et le premier modèle ajusté (m1):
Le graphique de gauche montre un graphique QQ des résidus mis à l'échelle pour détecter les écarts par rapport à la distribution attendue, et le graphique de droite représente les résidus par rapport aux valeurs prédites tout en effectuant une régression quantile pour détecter les écarts par rapport à l'uniformité (les lignes rouges doivent être horizontales et à 0,25 , 0,50 et 0,75).
En outre, le package a également des fonctions spécifiques pour tester la dispersion sur / sous et l'inflation zéro, entre autres:
la source