#####################################################################################################################################################
############################################################ Functions for mixed models ############################################################# 
#####################################################################################################################################################

# Intellectual property of Maryse Marchand (NRCAN, CFL)

## Fonction pour vérifier les suppositions du modèles mixte linéaire de façon graphique
# Homogénéité de la variance
# Normalité des résidus
# Normalité des effets aléatoires
# Évaluation de l'ajustement des données prédites-observées

## Arguments

# mod: objet de classe "lme", représentant un modèle mixte linéaire
# response: nom de la variable réponse du modèle "mod", entre guillemets
# level: niveau de groupement à utiliser pour obtenir les prédictions du modèle (level zéro = prédictions pour la population)
# ** fonctionne présentement pour les modèles avec 1 ou 2 niveau d'effets aléatoires **
# w_width et w_height: paramètres pour ajuster les dimensions de la fenêtre graphique

diag_plots <- function(mod, response, level, w_width = 10, w_height = 10) {
  
  # Vérifier classe de l'objet
  if(class(mod) != "lme")
    stop("mod must be of class lme")
  
  # 1 niveau d'effets aléatoire
  if (level == 1){
    
    # Calculer résidus
    res <- residuals(mod, level = level, type = "pearson")
    
    x11(w = w_width, h = w_height)
    par(mfrow = c(2,2))
    
    # Homogénéité de la variance
    plot(fitted(mod, level = level), res,
         xlab = "Fitted values", ylab = "Standardized residuals", main = "Homoscedasticity")
    abline(h=0, lty=2)
    abline(h=-3, lty=2, col="red")
    abline(h=3, lty=2, col="red")
    # Normalité des résidus
    qqnorm(res, main="Normality of residuals")
    qqline(res)
    # Normalité des effets aléatoires (sites)
    qqnorm(ranef(mod, level = level)$"(Intercept)", main="Normality of random effects")
    qqline(ranef(mod, level = level)$"(Intercept)")
    # Données prédites vs observées
    plot(fitted(mod), mod$data[response][,1], main = "Fitted vs observed values", xlab = "Fitted values", ylab = "Observed values"); abline(0,1)
    
  }
  
  # 2 niveaux d'effets aléatoires
  if (level == 2){
    
    # Calculer résidus
    res <- residuals(mod, level = level, type = "pearson")
    
    x11(w = w_width, h = w_height)
    par(mfrow = c(2,2))
    
    # Homogénéité de la variance
    plot(fitted(mod, level = 2), res,
         xlab = "Fitted values", ylab = "Standardized residuals", main = "Homoscedasticity")
    abline(h=0, lty=2)
    abline(h=-3, lty=2, col="red")
    abline(h=3, lty=2, col="red")
    # Normalité des résidus
    qqnorm(res, main="Normality of residuals")
    qqline(res)
    # Normalité des effets aléatoires 
    qqnorm(ranef(mod, level = 2)$"(Intercept)", main="Normality of random effects (level 2)")
    qqline(ranef(mod, level = 2)$"(Intercept)")
    qqnorm(ranef(mod, level = 1)$"(Intercept)", main="Normality of random effects (level 1)")
    qqline(ranef(mod, level = 1)$"(Intercept)")
    
  }
  
  if (!level %in% c(1,2)) stop("Only for 1 or 2 levels of random effects")
  
  
}