# Re-évaluer les intercepts aléatoires
summary(mnb1)$varcor

mnb1.taxless <- glm.nb(PLD ~ temp*feeding.type, data=inverts)

# Ici, parce que nous comparons un glmer avec un glm, nous devons faire quelque chose de différent que 
# anova(). Pour tester l'importance de l’intercept aléatoire, nous allons comparer la vraisemblance de 
# chaque modèle :
NL1 <- -logLik(mnb1)
NL0 <- -logLik(mnb1.taxless)
devdiff <- 2*(NL0-NL1)
dfdiff <- attr(NL1,"df")-attr(NL0,"df")
pchisq(devdiff,dfdiff,lower.tail=FALSE)

# Nous pourrions aussi comparer l'AIC du modèle avec (mnb1) et sans (mnb1.taxless) effets aléatoires avec 
# la fonction AICtab() 
AICtab(mnb1,mnb1.taxless) 
# Changement important du AIC si nous supprimons l’intercept aléatoire. Donc, ça vaut la peine de garder 
# cet effet.

# Graphique diagnostic 
locscaleplot(mnb1)

# Graphique des paramètres de variance
coefplot2(mnb1,ptype="vcov",intercept=TRUE, main="Random effect variance")

# Graphique des effets fixes
coefplot2(mnb1,intercept=TRUE,main="Fixed effect coefficient")

# Graphique des intercepts aléatoires 
pp <- list(layout.widths=list(left.padding=0, right.padding=0))
r2 <- ranef(mnb1,condVar=TRUE)
d2 <- dotplot(r2, par.settings=pp)
grid.arrange(d2$taxon,nrow=1)

# Évaluer pentes aléatoires
mnb2 <- glmer.nb(PLD ~ temp*feeding.type +
                   (PLD|taxon),
                 data=inverts)
                 
# Examiner composant de variance-covariance
summary(mnb2) # option 1
attr(VarCorr(mnb2)$taxon,"correlation") # option 2
printvc(mnb2) # option 3
# Forte corrélation entre les effets aléatoires -> pas assez de puissance pour tester pentes aléatoires

# Re-évaluer les effets fixes
# Remarque : pour utiliser la fonction de drop1 nous devons spécifier le paramètre thêta et exécuter 
# le modèle NB avec glmer :
theta.mnb1 <- theta.md(inverts$PLD, fitted(mnb1), dfr = df.residual(mnb1))
mnb1 <- glmer(PLD ~ temp*feeding.type +
                (1|taxon),
              data=inverts, family=negative.binomial(theta=theta.mnb1))

(dd_LRT <- drop1(mnb1,test="Chisq"))
(dd_AIC <- dfun(drop1(mnb1)))
# Lorsque l’interaction feeding.type x température est supprimée, dAIC change de plus de 2 unité → suggère
# de garder l'interaction dans le modèle