# Re-evaluating random intercepts
summary(mnb1)$varcor

mnb1.taxless <- glm.nb(PLD ~ temp*feeding.type,
                       data=inverts)
# Here, because we are comparing a glmer with a glm, we must do something different than calling anova()
# to test importance of random intercept. We will compare the likelihood of each model:
NL1 <- -logLik(mnb1)
NL0 <- -logLik(mnb1.taxless)
devdiff <- 2*(NL0-NL1)
dfdiff <- attr(NL1,"df")-attr(NL0,"df")
pchisq(devdiff,dfdiff,lower.tail=FALSE)

# Could also look at dAIC to compare model with (mnb1) and without (mnb1.taxless) random effects
# Using AICtab() function
AICtab(mnb1,mnb1.taxless) 
# Large change in AIC if drop random intercept. Therefore, worth keeping this in.

# Diagnostic plots 
locscaleplot(mnb1)

# Plotting variance terms 
coefplot2(mnb1,ptype="vcov",intercept=TRUE,main="Random effect variance")

# Plot of fixed effects 
coefplot2(mnb1,intercept=TRUE,main="Fixed effect coefficient")

# Plotting random intercepts 
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)

# Evaluate random slopes
mnb2 <- glmer.nb(PLD ~ temp*feeding.type +
                   (PLD|taxon),
                 data=inverts)
                 
# View variance-covariance components
summary(mnb2) # option 1
attr(VarCorr(mnb2)$taxon,"correlation") # option 2
printvc(mnb2) # option 3
# Strong correlation between random effects -> not enough power to test random slopes

# Re-evaluating fixed effects 
# Note: to run the drop1 we need to speficy the theta parameter and run the NB model with 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)))
# dAIC when feeding.type x temp interaction is dropped is greater than 2, suggest to keep 
# interaction in model