Sélection de modèles avec données climatiques
Pierre-Alexandre QUITTET
2026-01-22
DEMARCHE D’ANALYSE
L’idée de ce document est (1) d’identifier les variables climatiques qui influencent la masse et la taille des marmottons et (2) de déterminer dans quelles mesures chaque variable climatique contribue à la variation temporelle des traits.
A partir des analyses précédentes des variables climatiques,
disponible dans le script METEOROLIGCAL_DATA.R (https://rpubs.com/paquittet/1390293), j’ai
identifié 6 variables climatiques (2 par saison) d’intérêt :
Été
- Coordonnées de l’axe 1 de l’ACP (valeurs positives : étés chauds et secs)
- Précipitations estivales moyennes
Hiver
- Coordonnées de l’axe 1 de l’ACP (valeurs positives : NAO +, neige +, précipitations +)
- Température hivernale moyenne
Printemps
- Coordonnées de l’axe 1 de l’ACP (valeurs négatives : printemps chauds et secs)
- NDVI moyen du printemps
Dans ce document je vais :
- de-trend les indices/variables climatiques ayant une tendance temporelle, pour éviter toute corrélation fallacieuse entre \(y\) et la variable climatique
- réaliser une sélection de modèle basée sur le critère AIC pour identifier les variables climatiques d’intérêt
- réaliser une analyse de déviance ANODEV (Grosbois et al., 2008) qui permettra de quantifier la contribution de chaque variable climatique à la variation temporelle du trait.
INITIALISATION
# PACKAGES
library(tidyverse)
library(splines) # for function bf() to splined emergence date
# WD
dir = "C:/Users/quittet/Documents/marmot-quantitative-genetics"
# LOAD FUNCTIONS
source(paste0(dir, "/1_DATA_ANALYSIS/1_MAIN_SCRIPTS/FUNCTIONS.R"))
# LOAD DATA
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/pheno_abiotic.rds")
load("C:/Users/quittet/Documents/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/df_abiotic_short.rds")
DE-TRENDED VARIABLES
Ressources : Grosbois et al. (2008).
L’idée est d’enlever la tendance temporelle d’une variable climatique \(x\). Cela éviter les corrélations dues à des effets confondants entre l’effet temps de \(y\) et l’effet temps de \(x\).
On parlera de variable de-trended \(x_{d,i}\) à chaque date \(i\) telle que :
\[ x_{d,i} = x_i - (\alpha+\beta t_i) \]
avec \(\alpha\) et \(\beta\) les paramètres de régression de la covariable \(x\) et la tendance linéaire \(t\). La variable de-trended \(x_{d,i}\) est enfait les résidus de la régression du modèle ci-dessous. A noter qu’on peut aussi tester des effets quadratiques pour le de-trend.
Ensuite, on construit le modèle avec notre variable réponse en utilisant \(x_d\) à la place de \(x\) [@grosboisassessingimpactclimate2008] et en gardant l’effet \(t\) :
\[ y = a + b_1t+b_2x_{d} \]
Désormais, la variable de-trended s’interprète comme la variation temporelle de \(x\) autour de la tendance linéaire.
Tester les tendances temporelles
Je teste ici chaque variable climatique pour tester la présence d’une tendance temporelle. Par variable climatique, j’ai construit 3 modèles : constant, temps avec effet linéaire et temps avec effet quadratique. Chaque modèle a ensuite été comparé par AIC, et le modèle le plus parcimonieux (i.e. \(\Delta AIC\) < 2 avec le degré de liberté le plus petit) a été sélectionné.
# RENAME DF WITH CLIMATIC VARIABLES
df_clim <- df_abiotic_short
# CLIMATIC VARIABLES
ClimVars <- c(
"summer_index",
"summer_precipitation_y",
"spring_index",
"spring_ndvi_y",
"winter_index",
"winter_temperature_y"
)
# -------------------------------------------------------------------
# detrend_f
#
# Description :
# Cette fonction ajuste trois modèles pour une série temporelle :
# (1) constant : y ~ 1
# (2) linéaire : y ~ t
# (3) quadratique : y ~ t + t^2
# Elle sélectionne ensuite le modèle le plus parcimonieux selon la règle :
# - on calcule les AIC
# - on retient le modèle avec le moins de paramètres si ΔAIC < 2 du modèle le plus complexe
#
# Arguments :
# data : data.frame contenant les données
# var : nom (character) de la variable climatique à analyser
# time : nom (character) de la variable temporelle
#
# Retour :
# Liste contenant le modèle choisi, le type, les AIC, les p-values et les résidus.
# -------------------------------------------------------------------
detrend_f <- function(data, var, time) {
# Création des formules
f_const <- as.formula(paste(var, "~ 1"))
f_lin <- as.formula(paste(var, "~", time))
f_qua <- as.formula(paste(var, "~", time, "+ I(", time, "^2)"))
# Ajustement des modèles
mod_const <- lm(f_const, data = data, na.action = na.exclude)
mod_lin <- lm(f_lin, data = data, na.action = na.exclude)
mod_qua <- lm(f_qua, data = data, na.action = na.exclude)
# Calcul AIC et degrés de liberté (nombre de paramètres)
model_list <- list(constant = mod_const, linear = mod_lin, quadratic = mod_qua)
aic_list <- sapply(model_list, AIC)
df_list <- sapply(model_list, function(m) length(coef(m)))
# Identifier le modèle avec l'AIC le plus bas
min_aic <- min(aic_list)
# Règle de parcimonie : ΔAIC < 2 => choisir le modèle avec moins de paramètres
candidate_models <- names(aic_list)[(aic_list - min_aic) < 2]
best_model_name <- candidate_models[which.min(df_list[candidate_models])]
best_model <- model_list[[best_model_name]]
# Résumé et p-values
summary_model <- summary(best_model)
p_values <- summary_model$coefficients[, "Pr(>|t|)"]
# Résidus
data_res <- data.frame(
time = data[[time]],
res = residuals(best_model)
)
# Retourner résultats
return(
list(
mod_const = mod_const,
mod_lin = mod_lin,
mod_qua = mod_qua,
best_model = best_model,
best_type = best_model_name,
AIC = aic_list,
df = df_list,
summary = summary_model,
p_values = p_values,
residuals = data_res
)
)
}
Visualisation des tendances temporelles
On visualise les modèles fitted, et on regarde par variable quel modèle est le mieux expliqué la tendance (constant , effet linéaire ou effet quadratique).
## $summer_index
## [1] "quadratic"
##
## $summer_precipitation_y
## [1] "quadratic"
##
## $spring_index
## [1] "linear"
##
## $spring_ndvi_y
## [1] "linear"
##
## $winter_index
## [1] "constant"
##
## $winter_temperature_y
## [1] "constant"
Construction des variables de-trended
Je de-trend donc les variables de printemps et d’été en construisant
la variable .<variable>_dt, mais pas celles d’hiver
pour lesquelles le modèle constant est le meilleur modèle.
# VARIABLES TO DETREND
to_detrend_temp <- names(models_best)[models_best != "constant"]
# CONSTRUCT DETREND VARIABLES
for(i in to_detrend_temp) {
df_res_temp <- models_all[[i]][["residuals"]]
index_temp <- match(pheno_abiotic$annee_capture, df_res_temp$time)
pheno_abiotic[[paste0(i, "_dt")]] <- df_res_temp$res[index_temp]
}
Analyse des résidus
Si le de-trend a bien fonctionné, on est censé ne voir aucun pattern dans les résidus. Je plot les résidus de chaque variable ayant eu une trend linéaire significative :
La plupart des tendances ont bien disparues. Celle de l’indice ACP d’été est discutable, on dirait qu’on pourrait faire passer une relation quadratique en U.
MODEL SELECTION
We select the fixed part of the model with a frequentist approach.
Extract best model per trait
Data preparation
# FILTER DATA
pheno_abiotic_mod <- pheno_abiotic |>
filter(
!is.na(age_annee) &
!is.na(annee_emergence) &
!is.na(litter_size_naissance) &
!is.na(age_d) &
!is.na(presence_male_subordinates) &
age_annee == 0 &
age_d == 0 &
annee_emergence %in% 1997:2023 &
masse < 700
) |>
arrange(age_d) |>
filter(!duplicated(individual_id)) |>
mutate(
presence_male_sub = as.factor(presence_male_subordinates),
annee_emergence_num = as.numeric(as.character(annee_emergence))
)
nrow(pheno_abiotic_mod)
# REMOVE NAs
vars_model <- c(
# CONTROL VARIABLE
"annee_emergence_num",
"tibia",
"sexe",
"inbreeding",
"capture_time",
# SOCIAL
"presence_male_sub",
"litter_size_naissance",
"taille_groupe_naissance",
# CLIMATIC
"winter_snow_y",
"winter_nao_y",
"winter_temperature_y",
"winter_index_y",
"spring_index_y",
"spring_ndvi_y",
"summer_index_y",
"summer_index_dt",
"spring_index_dt",
"spring_ndvi_y_dt"
)
pheno_abiotic_mod <- pheno_abiotic_mod[complete.cases(pheno_abiotic_mod[, vars_model]), ]
# SCALE QUANTITATIVE PREDICTORS
for(i in vars_model[!vars_model %in% c("tibia", "sexe", "presence_male_sub")]) {
pheno_abiotic_mod[[i]] <- as.numeric(scale(pheno_abiotic_mod[[i]]))
}
Fixed effects selection with MuMIn::dredge()
We test the fixed part of the model with constant random effects :
mother, territory and year. We
set the argument REML = F to use ML methods.
Mass
# GET FULL MODEL
full_model_mass <- lmerTest::lmer(
masse ~
# CONTROL
inbreeding +
sexe +
capture_time +
annee_emergence_num +
I(annee_emergence_num^2) +
# SOCIAL
litter_size_naissance +
taille_groupe_naissance +
presence_male_sub +
# CLIMATIC
winter_index_y +
winter_temperature_y +
summer_index_dt +
summer_precipitation_y_dt +
spring_index_dt +
spring_ndvi_y_dt +
# SOCIAL x ENVIRONMENT
winter_index_y : presence_male_sub +
# RANDOM
(1|mother),
REML = F,
data = pheno_abiotic_mod
)
# MODEL SELECTION
options(na.action = "na.fail") # nécessaire
model_selection_masse <- MuMIn::dredge(
full_model_mass,
evaluate = TRUE,
fixed = c(
"annee_emergence_num",
"I(annee_emergence_num^2)",
"capture_time",
"inbreeding",
"sexe",
"litter_size_naissance",
"taille_groupe_naissance"
) #' ... note : fixed argument force the variables to be included in models
)
Ci-dessous la liste de l’ensemble des modèles avec un \(\Delta AICc < 2\) :
## (ntrcpt) nn_mrgnc_nm (nn_mrgnc_nm^2) cptr_tm nbrdng lttr_sz_nssnc prsnc_ml_sb sx sprng_ndx_dt sprng_ndv_dt smmr_ndx_dt
## 1 355.86 -5.17 8.46 23.49 10.54 -42.32 + + -12.45 9.37 7.3
## 2 354.06 -4.49 9.59 23.04 10.43 -42.91 + + -13.6 8.45 -
## 3 355.8 -5.57 8.33 23.15 10.56 -42.32 + + -12.25 9.42 7.52
## 4 354.92 -4.89 8.99 23.33 10.67 -42.64 + + -13 8.82 4.18
## 5 353.98 -4.89 9.48 22.66 10.42 -42.92 + + -13.42 8.51 -
## 6 354.93 -5.31 8.83 22.99 10.68 -42.62 + + -12.76 8.92 4.6## smmr_prcpttn_dt tll_grp_nssnc wntr_ndx wntr_tmprtr prsnc_ml_sb:wntr_ndx df logLik AICc delta weight
## 1 - -19.36 -5.51 12.82 + 17 -7220.99 14476.47 0 0.2
## 2 -16.21 -19.77 -6.88 15.44 + 17 -7221.09 14476.66 0.19 0.18
## 3 - -19.23 -8.5 13.29 - 16 -7222.32 14477.07 0.59 0.15
## 4 -8.55 -19.64 -6.24 14.54 + 18 -7220.52 14477.58 1.11 0.12
## 5 -16.45 -19.63 -10.1 15.92 - 16 -7222.59 14477.62 1.14 0.11
## 6 -8 -19.49 -9.25 14.91 - 17 -7221.9 14478.29 1.82 0.08
Le meilleur modèle (\(\Delta AICc < 2\) & avec le moins de paramètres) est donc :
# FORMULA OF BEST MODEL
formula_masse <- as.formula(
masse ~
# CONTROL VARIABLES
annee_emergence_num + I(annee_emergence_num^2) + capture_time + sexe + inbreeding +
# SOCIAL VARIABLE
litter_size_naissance + taille_groupe_naissance + presence_male_sub +
# CLIMATIC VARIABLES
spring_index_dt + spring_ndvi_y_dt + summer_index_dt + winter_temperature_y + winter_index_y +
# SOCIAL x ENVIRONMENT
winter_index_y:presence_male_sub +
# RANDOM EFFECTS
(1 | mother)
)
# BEST MODEL
best_model_masse <- lmerTest::lmer(formula_masse,
data = pheno_abiotic_mod)
summary(best_model_masse)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_masse
## Data: pheno_abiotic_mod
##
## REML criterion at convergence: 14379.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1974 -0.4881 -0.0111 0.5326 3.6276
##
## Random effects:
## Groups Name Variance Std.Dev.
## mother (Intercept) 4891 69.93
## Residual 4018 63.39
## Number of obs: 1272, groups: mother, 120
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 355.758 8.744 212.602 40.686 < 2e-16 ***
## annee_emergence_num -5.115 5.830 293.532 -0.877 0.38103
## I(annee_emergence_num^2) 8.525 4.191 723.685 2.034 0.04230 *
## capture_time 23.513 2.664 1249.122 8.827 < 2e-16 ***
## sexem 11.354 3.744 1152.325 3.033 0.00248 **
## inbreeding 10.531 4.055 769.606 2.597 0.00958 **
## litter_size_naissance -42.334 2.200 1232.103 -19.246 < 2e-16 ***
## taille_groupe_naissance -19.375 2.930 1254.617 -6.614 5.54e-11 ***
## presence_male_sub1 10.710 5.481 1247.553 1.954 0.05092 .
## spring_index_dt -12.440 2.165 1208.764 -5.745 1.16e-08 ***
## spring_ndvi_y_dt 9.376 2.375 1226.129 3.949 8.31e-05 ***
## summer_index_dt 7.306 2.263 1223.493 3.229 0.00128 **
## winter_temperature_y 12.825 2.547 1224.012 5.036 5.45e-07 ***
## winter_index_y -5.517 2.857 1220.309 -1.931 0.05375 .
## presence_male_sub1:winter_index_y -7.314 4.528 1229.594 -1.615 0.10646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Tibia
# GET FULL MODEL
full_model_tibia <- lmerTest::lmer(
tibia ~
# CONTROL
inbreeding +
sexe +
capture_time +
annee_emergence_num +
I(annee_emergence_num^2) +
# SOCIAL
litter_size_naissance +
taille_groupe_naissance +
presence_male_sub +
# CLIMATIC
winter_index_y +
winter_temperature_y +
summer_index_dt +
summer_precipitation_y_dt +
spring_index_dt +
spring_ndvi_y_dt +
# SOCIAL x ENVIRONMENT
winter_index_y : presence_male_sub +
# RANDOM
(1|mother),
REML = F,
data = pheno_abiotic_mod
)
# MODEL SELECTION
options(na.action = "na.fail") # nécessaire
model_selection_tibia <- MuMIn::dredge(
full_model_tibia,
fixed = c(
"annee_emergence_num",
"I(annee_emergence_num^2)",
"capture_time",
"inbreeding",
"sexe",
"litter_size_naissance",
"taille_groupe_naissance"
) #' ... note : fixed argument force the variables to be included in models
)
Ci-dessous la liste de l’ensemble des modèles avec un \(\Delta AICc < 2\) :
## (ntrcpt) nn_mrgnc_nm (nn_mrgnc_nm^2) cptr_tm nbrdng lttr_sz_nssnc prsnc_ml_sb sx sprng_ndx_dt sprng_ndv_dt smmr_ndx_dt
## 1 50.32 -0.36 1.24 0.6 0.15 -1.11 + + -0.37 1.04 -
## 2 50.27 -0.28 1.31 0.61 0.12 -1.08 + + -0.4 1.05 -
## 3 50.31 -0.31 1.27 0.6 0.12 -1.08 + + -0.38 1.07 -
## 4 50.28 -0.28 1.31 0.6 0.12 -1.08 + + -0.4 1.06 -0.05## smmr_prcpttn_dt tll_grp_nssnc wntr_ndx wntr_tmprtr prsnc_ml_sb:wntr_ndx df logLik AICc delta weight
## 1 - -1.14 -0.2 0.18 + 16 -3469.69 6971.81 0 0.21
## 2 - -1.1 -0.23 - + 15 -3470.76 6971.9 0.09 0.2
## 3 0.22 -1.11 -0.2 - + 16 -3470.21 6972.86 1.04 0.12
## 4 - -1.1 -0.23 - + 16 -3470.63 6973.7 1.89 0.08
Le meilleur modèle (\(\Delta AICc < 2\) & avec le moins de paramètres) est donc :
# FORMULA OF BEST MODEL
formula_tibia <- as.formula(
tibia ~
# CONTROL VARIABLES
annee_emergence_num + I(annee_emergence_num^2) +
capture_time + sexe + inbreeding +
# SOCIAL VARIABLE
litter_size_naissance + taille_groupe_naissance + presence_male_sub +
# CLIMATIC VARIABLES
spring_index_dt + spring_ndvi_y_dt + winter_index_y +
# SOCIAL x ENVIRONMENT
winter_index_y:presence_male_sub +
# RANDOM EFFECTS
(1 | mother)
)
# BEST MODEL
best_model_tibia <- lmerTest::lmer(formula_tibia,
data = pheno_abiotic_mod)
summary(best_model_tibia)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_tibia
## Data: pheno_abiotic_mod
##
## REML criterion at convergence: 6964.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6567 -0.5090 0.0951 0.5914 3.0427
##
## Random effects:
## Groups Name Variance Std.Dev.
## mother (Intercept) 7.157 2.675
## Residual 11.651 3.413
## Number of obs: 1272, groups: mother, 120
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 50.2676 0.3834 244.0912 131.097 < 2e-16 ***
## annee_emergence_num -0.2797 0.2545 234.1521 -1.099 0.272900
## I(annee_emergence_num^2) 1.3129 0.1958 487.6650 6.704 5.64e-11 ***
## capture_time 0.6107 0.1369 1160.6491 4.460 8.98e-06 ***
## sexem 1.1062 0.2008 1176.3016 5.509 4.44e-08 ***
## inbreeding 0.1253 0.1931 513.0596 0.649 0.516647
## litter_size_naissance -1.0828 0.1150 1258.2433 -9.418 < 2e-16 ***
## taille_groupe_naissance -1.1006 0.1491 1235.2571 -7.381 2.87e-13 ***
## presence_male_sub1 0.8336 0.2841 1258.4570 2.934 0.003411 **
## spring_index_dt -0.3994 0.1124 1238.8870 -3.553 0.000395 ***
## spring_ndvi_y_dt 1.0534 0.1249 1256.6065 8.431 < 2e-16 ***
## winter_index_y -0.2341 0.1495 1254.8722 -1.567 0.117476
## presence_male_sub1:winter_index_y -0.5895 0.2377 1257.6338 -2.480 0.013260 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MODEL DIAGNOSIS
Residuals diagnosis
Mass
We observe a slight heteroscedasticity with increased résiduals values with increased predicted mass.
Tibia
Multicollinearity
We check the collinearity of the predictors using the VIF score with
the car::vif() function.
Mass
## vif_score
## annee_emergence_num 1.320381
## I(annee_emergence_num^2) 1.287246
## capture_time 1.050813
## sexe 1.007083
## inbreeding 1.083775
## litter_size_naissance 1.046384
## taille_groupe_naissance 1.669826
## presence_male_sub 1.500745
## spring_index_dt 1.133404
## spring_ndvi_y_dt 1.121504
## summer_index_dt 1.247455
## winter_temperature_y 1.569455
## winter_index_y 1.944845
## presence_male_sub:winter_index_y 1.924264
Tibia
## vif_score
## annee_emergence_num 1.291835
## I(annee_emergence_num^2) 1.256846
## capture_time 1.042519
## sexe 1.006654
## inbreeding 1.057905
## litter_size_naissance 1.022249
## taille_groupe_naissance 1.552109
## presence_male_sub 1.445396
## spring_index_dt 1.083311
## spring_ndvi_y_dt 1.121525
## winter_index_y 1.883889
## presence_male_sub:winter_index_y 1.878828All is good !
Partial residuals
We evaluate the partial residuals of each predictors to check for potential highly non-linear variation between the response variable and the different predictors.
plot_residus_partiels <- function(modele, predicteur, ylab, xlab) {
# VERIFICATION
# Récupère les noms des coefficients fixes du modèle
variables_modele <- names(lme4::fixef(modele))
# Supprime l'intercept
variables_modele <- variables_modele[variables_modele != "(Intercept)"]
# Vérifie que le prédicteur est dans le modèle
if (!predicteur %in% variables_modele) {
stop("Le prédicteur spécifié n'est pas dans le modèle.")
}
# DONNEES DU MODELE
donnees <- model.frame(modele)
# RESIDUS
residus <- residuals(modele) # prend en compte les BLUPs
# COEFFICIENT FIXE DU PREDICTEUR
beta <- lme4::fixef(modele)[predicteur]
# VALEURS DU PREDICTEUR
x <- donnees[[predicteur]]
# RESIDUS PARTIELS
residus_partiels <- residus + beta * x
df_ggplot <- data.frame(x = x, residus_partiels = residus_partiels)
# Nom de l'axe x par défaut
if (is.null(xlab)) xlab <- predicteur
# PLOT
ggplot(data = df_ggplot, aes(x = x, y = residus_partiels)) +
geom_jitter(width = 0.1, col = "steelblue", alpha = 0.3) +
geom_smooth(method = "loess", col = "#1F1F1F", alpha = 0.5) +
labs(y = ylab, x = xlab) +
theme_classic(base_size = 20)
}
Mass
Pas de non-linéarité marquée des prédicteurs.
Tibia
Pas de non-linéarité marquée des prédicteurs.
PLOT CONDITIONAL EFFECTS OF PREDICTORS
We represent the significant interaction between snow cover and presence of male helper for both traits.
# INTERACTION EFFECTS
# Mass
marginaleffects::slopes(best_model_masse,
variables = "winter_index_y",
by = "presence_male_sub")##
## presence_male_sub Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 -5.52 2.86 -1.93 0.0535 4.2 -11.1 0.0835
## 1 -12.83 3.46 -3.71 <0.001 12.3 -19.6 -6.0597
##
## Term: winter_index_y
## Type: response
## Comparison: dY/dX# Tibia
marginaleffects::slopes(best_model_tibia,
variables = "winter_index_y",
by = "presence_male_sub")##
## presence_male_sub Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 -0.234 0.149 -1.57 0.117 3.1 -0.527 0.0587
## 1 -0.824 0.177 -4.65 <0.001 18.2 -1.170 -0.4768
##
## Term: winter_index_y
## Type: response
## Comparison: dY/dX
ANODEV
L’analyse de déviance est une démarche statistique qui permet de quantifier la contribution de variables climatiques à la variation temporelle d’un trait. L’idée est de comparer un modèle de base incluant les variables climatiques comme covariables , un modèle constant et un modèle full-time dépendant (i.e. effet année en facteur). Concrètement, on va réaliser chacun de ces 3 modèles, et extrait la vraisemblance. Les statistiques \(F\) permettent de tester l’hypothèse nulle d’absence d’effet des variables climatiques sur le trait. De plus, le \(R^2\) de l’ANODEV quantifie la proportion de la variation temporelle du trait expliquée par chaque variable climatique.
La statistique de l’ANODEV est :
\[ F = \frac{(D_0 - D_c)/(df_0 - df_c)}{(D_c - D_t)/(df_c - df_t)} \]
Avec :
- \(D_0\) : déviance du modèle constant
- \(D_c\) : déviance du modèle climatique
- \(D_t\) : déviance du modèle full time
- \(df_0, df_c, df_t\) : degrés de liberté associés
et
\[ R^2_{dev} = \frac{D_0 - D_c}{D_0 - D_t} \]
Soit en R :
ANODEV <- function(data,
response,
fixed_effects,
random_effects,
clim_vars,
time_var) {
# STOCKAGE RESULTATS
results <- data.frame(
ClimVar = character(),
R2dev = numeric(),
F_stat = numeric(),
p = numeric(),
n = numeric(),
J = numeric(),
stringsAsFactors = FALSE
)
# MODELE CONSTANT (Fcst)
form_const <- as.formula(
paste(response, "~",
paste(fixed_effects, collapse = " + "),
"+", paste(random_effects, collapse = " + "))
)
m0 <- lme4::lmer(form_const, data = data, REML = FALSE)
# ANODEV PAR VARIABLES CLIMATIQUES
for (clim in clim_vars) {
# MODELE CLIMAT (Fco)
form_clim <- as.formula(
paste(response, "~",
paste(fixed_effects, collapse = " + "),
"+", clim,
"+", paste(random_effects, collapse = " + "))
)
mc <- lme4::lmer(form_clim, data = data, REML = FALSE)
# MODELE TEMPS (Ft)
form_time <- as.formula(
paste(response, "~",
paste(fixed_effects, collapse = " + "),
"+ as.factor(", time_var, ")",
"+", paste(random_effects, collapse = " + "))
)
mt <- lme4::lmer(form_time, data = data, REML = FALSE)
# Deviances
D0 <- deviance(m0)
Dc <- deviance(mc)
Dt <- deviance(mt)
# Df
df0 <- attr(x = logLik(m0), which = "df")
dfc <- attr(x = logLik(mc), which = "df")
dft <- attr(x = logLik(mt), which = "df")
# Paramètres Skalski
J <- dfc - df0 + 1 # + 1 pour l'intercept
n <- dft - df0 + 1 # + 1 pour l'intercept
# F
F_stat <- ((D0 - Dc)/(J - 1)) /
((Dc - Dt)/(n - J))
# p-value
p_val <- pf(F_stat,
df1 = J - 1,
df2 = n - J,
lower.tail = FALSE)
# R2 dev
R2dev <- (D0 - Dc)/(D0 - Dt)
# Stockage
results <- rbind(results, data.frame(
ClimVar = clim,
R2dev = round(R2dev, 2),
F_stat = round(F_stat, 2),
p = round(p_val, 3),
n = n,
J = J
))
}
# MODELE CLIMAT TOTAL
clim_all <- paste(clim_vars, collapse = " + ")
form_clim_all <- as.formula(
paste(response, "~",
paste(fixed_effects, collapse = " + "),
"+", clim_all,
"+", paste(random_effects, collapse = " + "))
)
mc_all <- lme4::lmer(form_clim_all, data = data, REML = FALSE)
# MODELE TEMPS TOTAL
form_time_all <- as.formula(
paste(response, "~",
paste(fixed_effects, collapse = " + "),
"+ as.factor(", time_var, ")",
"+", paste(random_effects, collapse = " + "))
)
mt_all <- lme4::lmer(form_time_all, data = data, REML = FALSE)
# Deviances
D0 <- deviance(m0)
Dc <- deviance(mc_all)
Dt <- deviance(mt_all)
# Df
df0 <- attr(x = logLik(m0), which = "df")
dfc <- attr(x = logLik(mc_all), which = "df")
dft <- attr(x = logLik(mt_all), which = "df")
# Paramètres Skalski
J <- dfc - df0 + 1
n <- dft - df0 + 1
# F total
F_stat <- ((D0 - Dc)/(J - 1)) /
((Dc - Dt)/(n - J))
# p total
p_val <- pf(F_stat,
df1 = J - 1,
df2 = n - J,
lower.tail = FALSE)
# R2 total
R2dev <- (D0 - Dc)/(D0 - Dt)
# Ajout ligne TOTAL
results <- rbind(results, data.frame(
ClimVar = "TOTAL",
R2dev = round(R2dev, 2),
F_stat = round(F_stat, 2),
p = round(p_val, 3),
n = n,
J = J
))
return(results)
}
Masse
## ClimVar R2dev F_stat p n J
## 1 winter_index_y + presence_male_sub:winter_index_y 0.16 2.26 0.126 27 3
## 2 summer_index_dt 0.00 0.07 0.798 27 2
## 3 spring_index_dt 0.16 4.91 0.036 27 2
## 4 spring_ndvi_y_dt 0.10 2.80 0.106 27 2
## 5 TOTAL 0.40 2.81 0.043 27 6
Tibia
## ClimVar R2dev F_stat p n J
## 1 winter_index_y + presence_male_sub:winter_index_y 0.07 0.90 0.419 27 3
## 2 spring_index_dt 0.02 0.47 0.500 27 2
## 3 spring_ndvi_y_dt 0.22 7.22 0.013 27 2
## 4 TOTAL 0.30 2.37 0.083 27 5