grpVars <- tapply(inverts$PLD, inverts$tft, var)
summary(grpVars)

grpMeans <- tapply(inverts$PLD,inverts$tft, mean)
summary(grpMeans)

# Quasi-Poisson
lm1 <- lm(grpVars~grpMeans-1) 
phi.fit <- coef(lm1)
# Le -1 spécifie un modèle avec l'intercept fixer à zéro

# Binomial négative
lm2 <- lm(grpVars ~ I(grpMeans^2) + offset(grpMeans)-1)
k.fit <- 1/coef(lm2)
# Offset () est utilisé pour fixer la moyenne de chaque groupe à 1

# Ajustement Loess non-paramétrique
Lfit <- loess(grpVars~grpMeans)

plot(grpVars ~ grpMeans, xlab="group means", ylab="group variances" )
abline(a=0,b=1, lty=2)
text(60,200,"Poisson")
curve(phi.fit*x, col=2,add=TRUE)
# bquote()est utilisé pour remplacer des valeurs numériques dans les équations avec des symboles
text(60,800, bquote(paste("QP: ",sigma^2==.(round(phi.fit,1))*mu)),col=2)
curve(x*(1+x/k.fit),col=4,add=TRUE)
text(60,1600,paste("NB: k=",round(k.fit,1),sep=""),col=4)
mvec <- 0:120
lines(mvec,predict(Lfit,mvec),col=5)
text(50,1300,"loess",col=5)

# Poisson GLMM 
mp1 <- glmer(PLD ~ temp*feeding.type +
               (1|taxon),
             data=inverts, family="poisson")
overdisp_fun(mp1)
# rapport est significativement > 1

# NB GLMM
mnb1 <- glmer.nb(PLD ~ temp*feeding.type +
                   (1|taxon),
                 data=inverts)
overdisp_fun(mnb1)
# Semble bon!