# Packages (installer si nécessaire)

# install.packages(c("lme4","glmmTMB","MASS","DHARMa","performance","broom.mixed","ggplot2"))

library(lme4)

library(glmmTMB)

library(MASS) # glm.nb si pas d'effet aléatoire

library(DHARMa)

library(performance)

library(broom.mixed)

library(ggplot2)

# -------------------------

# Préparation des données

# -------------------------

# On suppose que ta table s'appelle `data` et contient Danaus_per et Distance

# Si tu as un identifiant de groupe (ex: site, plot), mets son nom dans group_var.

# Sinon laisse group_var <- NA

group_var <- "site" # <-- change à NA si tu n'as pas de variable de groupe

# Faire une copie et centrer Distance (meilleure stabilité numérique)

dat <- data

dat$Distance_c <- scale(dat$Distance, center = TRUE, scale = TRUE) # centrée & réduite

# -------------------------

# 1) Choix : GLMM (si group_var) ou GLM (si pas de groupe)

# -------------------------
if(!is.na(group_var) && group_var %in% names(dat)){

# GLMM: random intercept par groupe

cat("-> On run un GLMM (random intercept):\n")

# 1.a Poisson (lme4)

m_pois <- try(glmer(Danaus_per ~ Distance_c + (1 | get(group_var)),

data = dat, family = poisson(link = "log")), silent = TRUE)

# 1.b Negative Binomial (glmmTMB)

m_nb <- try(glmmTMB(Danaus_per ~ Distance_c + (1 | get(group_var)),

data = dat, family = nbinom2(link = "log")), silent = TRUE)

# 1.c Zero-inflated NB (si beaucoup de zéros)

m_zinb <- try(glmmTMB(Danaus_per ~ Distance_c + (1 | get(group_var)),

ziformula = ~1,

data = dat, family = nbinom2(link = "log")), silent = TRUE)

# Résumés (afficher ceux qui ont convergé)

if(!inherits(m_pois, "try-error")) { cat("\nRésumé Poisson (glmer):\n"); print(summary(m_pois)) }

if(!inherits(m_nb, "try-error")) { cat("\nRésumé NB (glmmTMB):\n"); print(summary(m_nb)) }

if(!inherits(m_zinb, "try-error")) { cat("\nRésumé ZINB (glmmTMB):\n"); print(summary(m_zinb)) }

# Comparer AIC si modèles valides

mods <- list()

if(!inherits(m_pois, "try-error")) mods$pois <- m_pois

if(!inherits(m_nb, "try-error")) mods$nb <- m_nb

if(length(mods) > 1){

cat("\nAIC comparison:\n")

print(AIC(mods$pois, mods$nb, mods$zinb)[, c("df","AIC")])

# Test surdispersion (pour Poisson et NB)

if(!inherits(m_pois, "try-error")){

od_pois <- sum(residuals(m_pois, type = "pearson")^2) / df.residual(m_pois)

cat("\nOverdispersion ratio (Poisson):", round(od_pois, 3), "\n")

if(!inherits(m_nb, "try-error")){

# DHARMa test dispersion

sim_nb <- simulateResiduals(fittedModel = m_nb, plot = FALSE)

cat("DHARMa test dispersion (NB) p-value:\n"); print(testDispersion(sim_nb))

cat("DHARMa test zero inflation (NB):\n"); print(testZeroInflation(sim_nb))

# Extraire coefficients et IRR pour le meilleur modèle (ici on prend NB si disponible)

best_model <- if(!inherits(m_nb, "try-error")) m_nb else m_pois

cat("\nCoefficient summary (best model):\n")

co <- summary(best_model)$coefficients

# Pour glmmTMB, cond part accessible via $cond

if("glmmTMB" %in% class(best_model)){

co_fix <- summary(best_model)$coefficients$cond

co_fix <- co

irr <- exp(co_fix[, "Estimate"])

res_table <- data.frame(Estimate = co_fix[, "Estimate"],

SE = co_fix[, "Std. Error"],

z = co_fix[, if("z value" %in% colnames(co_fix)) "z value" else "z"],

p = co_fix[, if("Pr(>|z|)" %in% colnames(co_fix)) "Pr(>|z|)" else "Pr(>|z|)"],

IRR = irr,

row.names = rownames(co_fix))

print(res_table)

# Effets aléatoires (variance par groupe)

cat("\nVarCorr (random effects):\n"); print(VarCorr(best_model))

# Prédictions marginales (population-level) pour Distance

newd <- data.frame(Distance_c = seq(min(dat$Distance_c, na.rm=TRUE),

max(dat$Distance_c, na.rm=TRUE), length.out = 30))

pred <- predict(best_model, newdata = newd, type = "response", allow.new.levels = TRUE, se.fit =
FALSE)

newd$pred <- pred

# Remet Distance à l'échelle originale pour tracer

# (On récupère centre & scale)

center_val <- attr(dat$Distance_c, "scaled:center")

scale_val <- attr(dat$Distance_c, "scaled:scale")

newd$Distance_orig <- newd$Distance_c * scale_val + center_val

ggplot(newd, aes(x = Distance_orig, y = pred)) +

geom_line() + labs(x = "Distance", y = "Predicted count (Danaus_per)",

title = "Prédiction (modèle choisi)") + theme_minimal()

} else {

# Pas d'effet aléatoire -> GLM

cat("-> Pas de variable de groupe détectée : on run un GLM (Poisson / NB)\n")

# Poisson GLM

m_pois_glm <- glm(Danaus_per ~ Distance_c, data = dat, family = poisson(link = "log"))

summary(m_pois_glm)

# Vérifier surdispersion

od_glm <- sum(residuals(m_pois_glm, type = "pearson")^2) / df.residual(m_pois_glm)

cat("Overdispersion ratio (Poisson GLM):", round(od_glm,3), "\n")

# Si overdispersion -> glm.nb

m_nb_glm <- glm.nb(Danaus_per ~ Distance_c, data = dat)

summary(m_nb_glm)

# Zero-inflated (ps: glmmTMB peut être utilisé même sans random)

m_zinb_glm <- glmmTMB(Danaus_per ~ Distance_c,

ziformula = ~1,
 data = dat,

family = nbinom2(link = "log"))

summary(m_zinb_glm)

# Comparaison AIC

cat("\nAIC (GLM Poisson, GLM NB, ZINB):\n")

print(AIC(m_pois_glm, m_nb_glm, m_zinb_glm)[, c("df","AIC")])

# Coeffs + IRR (prendre NB si overdispersion)

best_model <- if(od_glm > 1.5) m_nb_glm else m_pois_glm

co_fix <- coef(summary(best_model))

irr <- exp(co_fix[, "Estimate"])

res_table <- data.frame(Estimate = co_fix[, "Estimate"],

SE = co_fix[, "Std. Error"],

z = co_fix[, if("z value" %in% colnames(co_fix)) "z value" else "t value"],

p = co_fix[, if("Pr(>|z|)" %in% colnames(co_fix)) "Pr(>|z|)" else NA],

IRR = irr,

row.names = rownames(co_fix))

print(res_table)

# Prévision simple

newd <- data.frame(Distance_c = seq(min(dat$Distance_c, na.rm=TRUE),

max(dat$Distance_c, na.rm=TRUE), length.out = 30))

newd$pred <- predict(best_model, newdata = newd, type = "response")

center_val <- attr(dat$Distance_c, "scaled:center")

newd$Distance_orig <- newd$Distance_c * scale_val + center_val

ggplot(newd, aes(x = Distance_orig, y = pred)) +

geom_line() + labs(x = "Distance", y = "Predicted count (Danaus_per)",

title = "Prédiction (modèle choisi)") + theme_minimal()

# -------------------------

# Fin du script

# -------------------------