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 +)
- Coordonnées de l’axe 2 de l’ACP (valeurs positives : température +)
Printemps
- Coordonnées de l’axe 1 de l’ACP (valeurs négatives : printemps chauds et secs)
- NDVI moyen du printemps
Démarche d’analyse :
- 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
library(kableExtra) # for nice table
# 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-dessus. 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\) 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
# SCALE
df_clim <- as.data.frame(apply(df_clim[-1], 2, scale))
df_clim$year <- df_abiotic_short$year
# CLIMATIC VARIABLES
ClimVars <- c(
"summer_index_1",
"summer_index_2",
"spring_index_1",
"spring_index_2",
"winter_index_1",
"winter_index_2"
)
# -------------------------------------------------------------------
# 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]]
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
## [1] "quadratic"
##
## $summer_index_2
## [1] "linear"
##
## $spring_index_1
## [1] "linear"
##
## $spring_index_2
## [1] "linear"
##
## $winter_index_1
## [1] "constant"
##
## $winter_index_2
## [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.
# GET SCALED CLIMATIC VARIABLES
to_scale_temp <- c("winter_index_1_y", "winter_index_2_y", "winter_nao_y")
index_temp <- match(pheno_abiotic$annee_capture, df_clim$year)
pheno_abiotic[to_scale_temp] <- df_clim[index_temp, c("winter_index_1", "winter_index_2", "winter_nao_y")]
# 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_index_1_y",
"winter_index_2_y",
# "spring_index_1_y",
# "spring_ndvi_y",
"spring_index_1_dt",
"spring_index_2_dt",
# "spring_ndvi_y_dt",
# "summer_index_1_y",
# "summer_precipitation_y_dt",
"summer_index_1_dt",
"summer_index_2_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", "winter_index_1_y", "winter_index_2_y", "spring_index_1_dt", "spring_index_2_dt", "summer_index_1_dt", "summer_index_2_dt")]) {
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. We set the argument REML = F to use ML
methods for model selection.
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
summer_index_1_dt +
summer_index_2_dt +
# summer_precipitation_y_dt +
winter_index_1_y +
winter_index_2_y +
spring_index_1_dt +
spring_index_2_dt +
# spring_ndvi_y_dt +
# SOCIAL x ENVIRONMENT
winter_index_1_y : presence_male_sub +
winter_index_2_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(
"I(annee_emergence_num^2)",
"annee_emergence_num",
"capture_time",
"inbreeding",
"sexe",
"litter_size_naissance",
"taille_groupe_naissance"
) #' ... note : fixed argument force the variables to be included in models
)
|
Fixed part
|
Summer
|
Spring
|
Winter
|
|||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| \(\mu\) | year | year\(^2\) | Capture time | Inbr. | Litter size | Group size | Male sub | Sex | Su. index1 | Su. index2 | Sp. index1 | Sp. index2 | Wi. index1 | Wi. index2 | (Sub × index1) | (Sub × index2) | AICc | df | \(\Delta\) | w |
| 345.12 | 1.24 | 17.13 | 24.24 | 11.12 | -44.18 | -20.31 | + | + | - | -14.17 | -11.5 | 10 | -8.78 | 27.4 | + | + | 14464.27 | 18 | 0 | 0.33 |
| 345.3 | 0.91 | 16.91 | 24.58 | 11.72 | -43.98 | -20.35 | + | + | 2.95 | -12.1 | -11.21 | 9.49 | -8.35 | 27.27 | + | + | 14464.8 | 19 | 0.53 | 0.26 |
| 345.56 | 0.91 | 16.54 | 23.81 | 10.81 | -44.26 | -20.08 | + | + | - | -14.18 | -11.11 | 10.24 | -12.65 | 27.11 | - | + | 14466.42 | 17 | 2.15 | 0.11 |
| 345.74 | 0.58 | 16.33 | 24.16 | 11.4 | -44.07 | -20.12 | + | + | 2.93 | -12.12 | -10.82 | 9.73 | -12.23 | 26.98 | - | + | 14466.97 | 18 | 2.7 | 0.09 |
| 346.41 | 1.29 | 16.26 | 24.38 | 10.38 | -44.04 | -20.03 | + | + | - | -14.39 | -11.49 | 10.74 | -9.68 | 22.4 | + | - | 14467.05 | 17 | 2.78 | 0.08 |
| 346.64 | 1 | 15.84 | 24 | 10.19 | -44.13 | -19.87 | + | + | - | -14.37 | -11.15 | 10.87 | -12.95 | 22.7 | - | - | 14468.11 | 16 | 3.84 | 0.05 |
| 346.62 | 1.06 | 16.05 | 24.63 | 10.75 | -43.9 | -20.04 | + | + | 2.04 | -12.97 | -11.29 | 10.44 | -9.44 | 21.99 | + | - | 14468.36 | 18 | 4.09 | 0.04 |
| 346.86 | 0.77 | 15.63 | 24.26 | 10.57 | -43.98 | -19.88 | + | + | 2.12 | -12.9 | -10.94 | 10.55 | -12.66 | 22.27 | - | - | 14469.36 | 17 | 5.09 | 0.03 |
| 350.92 | 1.07 | 17.07 | 24.09 | 9.28 | -43.25 | -16.01 | - | + | - | -13.52 | -11.58 | 10.77 | -12.93 | 21.07 | - | - | 14471.94 | 15 | 7.67 | 0.01 |
| 351.15 | 0.85 | 16.88 | 24.34 | 9.64 | -43.1 | -16 | - | + | 2.02 | -12.11 | -11.39 | 10.46 | -12.66 | 20.66 | - | - | 14473.26 | 16 | 8.99 | 0 |
NOTE: Italic
variables represent de-trended climatic predictors, bold correspond to
the best parsimonious model, + means interaction or factor
has been added
|
||||||||||||||||||||
Le meilleur modèle (\(\Delta AICc < 2\) & avec le moins de paramètres) est donc :
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_masse
## Data: pheno_abiotic_mod
##
## REML criterion at convergence: 14359.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1944 -0.5034 0.0052 0.5302 3.4535
##
## Random effects:
## Groups Name Variance Std.Dev.
## mother (Intercept) 5080 71.27
## Residual 3959 62.92
## Number of obs: 1272, groups: mother, 120
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 344.033 8.860 210.719 38.829 < 2e-16 ***
## annee_emergence_num 1.313 5.931 305.535 0.221 0.82497
## I(annee_emergence_num^2) 17.193 4.117 711.563 4.177 3.33e-05 ***
## capture_time 24.260 2.664 1249.804 9.106 < 2e-16 ***
## sexem 11.054 3.719 1149.828 2.972 0.00302 **
## inbreeding 11.112 4.044 809.486 2.748 0.00614 **
## litter_size_naissance -44.196 2.232 1230.504 -19.800 < 2e-16 ***
## taille_groupe_naissance -20.331 2.924 1252.900 -6.952 5.77e-12 ***
## presence_male_sub1 13.830 5.518 1247.152 2.506 0.01232 *
## summer_index_2_dt -16.007 3.578 1205.030 -4.473 8.42e-06 ***
## winter_index_1_y -9.344 3.233 1213.817 -2.890 0.00392 **
## winter_index_2_y 24.761 3.457 1227.793 7.163 1.36e-12 ***
## spring_index_1_dt -12.649 2.522 1214.683 -5.015 6.09e-07 ***
## spring_index_2_dt 12.502 3.069 1222.343 4.073 4.93e-05 ***
## presence_male_sub1:winter_index_1_y -10.021 4.913 1226.713 -2.040 0.04161 *
## presence_male_sub1:winter_index_2_y -8.680 3.946 1226.041 -2.200 0.02802 *
## ---
## 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
summer_index_1_dt +
summer_index_2_dt +
# summer_precipitation_y_dt +
winter_index_1_y +
winter_index_2_y +
spring_index_1_dt +
spring_index_2_dt +
# spring_ndvi_y_dt +
# SOCIAL x ENVIRONMENT
winter_index_1_y : presence_male_sub +
winter_index_2_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(
"I(annee_emergence_num^2)",
"annee_emergence_num",
"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\) :
|
Fixed part
|
Summer
|
Spring
|
Winter
|
|||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| \(\mu\) | year | year\(^2\) | Capture time | Inbr. | Litter size | Group size | Male sub | Sex | Su. index1 | Su. index2 | Sp. index1 | Sp. index2 | Wi. index1 | Wi. index2 | (Sub × index1) | (Sub × index2) | AICc | df | \(\Delta\) | w |
| 49.85 | -0.16 | 1.59 | 0.72 | 0.23 | -1.14 | -1.25 | + | + | - | -0.62 | - | 1.16 | -0.53 | 0.93 | + | + | 6958.49 | 17 | 0 | 0.2 |
| 49.89 | -0.14 | 1.57 | 0.72 | 0.2 | -1.14 | -1.24 | + | + | - | -0.63 | - | 1.19 | -0.56 | 0.73 | + | - | 6959.16 | 16 | 0.68 | 0.14 |
| 49.81 | -0.11 | 1.64 | 0.71 | 0.23 | -1.13 | -1.22 | + | + | - | -0.64 | -0.12 | 1.14 | -0.53 | 0.91 | + | + | 6959.49 | 18 | 1.01 | 0.12 |
| 49.84 | -0.16 | 1.6 | 0.7 | 0.21 | -1.15 | -1.25 | + | + | -0.1 | -0.69 | - | 1.18 | -0.54 | 0.94 | + | + | 6959.85 | 18 | 1.36 | 0.1 |
| 49.88 | -0.14 | 1.57 | 0.7 | 0.18 | -1.15 | -1.24 | + | + | -0.13 | -0.72 | - | 1.21 | -0.58 | 0.76 | + | - | 6960.08 | 17 | 1.6 | 0.09 |
| 49.85 | -0.09 | 1.61 | 0.72 | 0.2 | -1.13 | -1.21 | + | + | - | -0.64 | -0.12 | 1.17 | -0.56 | 0.72 | + | - | 6960.16 | 17 | 1.67 | 0.09 |
| 49.8 | -0.1 | 1.65 | 0.7 | 0.21 | -1.14 | -1.22 | + | + | -0.12 | -0.72 | -0.14 | 1.16 | -0.54 | 0.92 | + | + | 6960.65 | 19 | 2.17 | 0.07 |
| 49.9 | -0.15 | 1.55 | 0.7 | 0.19 | -1.14 | -1.23 | + | + | - | -0.63 | - | 1.2 | -0.75 | 0.75 | - | - | 6960.75 | 15 | 2.26 | 0.06 |
| 49.84 | -0.08 | 1.63 | 0.7 | 0.18 | -1.14 | -1.21 | + | + | -0.15 | -0.75 | -0.14 | 1.19 | -0.58 | 0.75 | + | - | 6960.82 | 18 | 2.33 | 0.06 |
| 49.87 | -0.16 | 1.57 | 0.69 | 0.21 | -1.15 | -1.23 | + | + | - | -0.62 | - | 1.17 | -0.74 | 0.91 | - | + | 6960.91 | 16 | 2.43 | 0.06 |
NOTE: Italic
variables represent de-trended climatic predictors, bold correspond to
the best parsimonious model, + means interaction or factor
has been added
|
||||||||||||||||||||
Le meilleur modèle (\(\Delta AICc < 2\) & avec le moins de paramètres) est donc :
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_tibia
## Data: pheno_abiotic_mod
##
## REML criterion at convergence: 6948.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4980 -0.5183 0.0776 0.5907 3.0408
##
## Random effects:
## Groups Name Variance Std.Dev.
## mother (Intercept) 7.192 2.682
## Residual 11.494 3.390
## Number of obs: 1272, groups: mother, 120
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 49.8440 0.3847 241.6745 129.570 < 2e-16 ***
## annee_emergence_num -0.1597 0.2586 239.0295 -0.618 0.537366
## I(annee_emergence_num^2) 1.5940 0.1904 448.7746 8.371 7.32e-16 ***
## capture_time 0.7172 0.1374 1159.5261 5.220 2.11e-07 ***
## sexem 1.0630 0.1996 1173.4272 5.326 1.20e-07 ***
## inbreeding 0.2282 0.1946 506.5236 1.173 0.241348
## litter_size_naissance -1.1422 0.1175 1256.7427 -9.720 < 2e-16 ***
## taille_groupe_naissance -1.2481 0.1490 1229.1495 -8.378 < 2e-16 ***
## presence_male_sub1 1.0680 0.2885 1253.2919 3.702 0.000223 ***
## summer_index_2_dt -0.6964 0.1890 1238.1943 -3.684 0.000239 ***
## spring_index_2_dt 1.4518 0.1596 1252.7137 9.099 < 2e-16 ***
## winter_index_1_y -0.5638 0.1712 1249.3392 -3.292 0.001022 **
## winter_index_2_y 0.8408 0.1818 1255.0260 4.624 4.15e-06 ***
## presence_male_sub1:winter_index_1_y -0.5451 0.2585 1256.4778 -2.109 0.035137 *
## presence_male_sub1:winter_index_2_y -0.3462 0.2084 1255.1752 -1.661 0.096997 .
## ---
## 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.347524
## I(annee_emergence_num^2) 1.239025
## capture_time 1.061332
## sexe 1.008561
## inbreeding 1.078593
## litter_size_naissance 1.091324
## taille_groupe_naissance 1.684839
## presence_male_sub 1.540098
## summer_index_2_dt 2.485242
## winter_index_1_y 2.119555
## winter_index_2_y 3.448705
## spring_index_1_dt 1.262873
## spring_index_2_dt 1.379971
## presence_male_sub:winter_index_1_y 1.893151
## presence_male_sub:winter_index_2_y 2.250988
Tibia
## vif_score
## annee_emergence_num 1.337235
## I(annee_emergence_num^2) 1.195144
## capture_time 1.061023
## sexe 1.007513
## inbreeding 1.080959
## litter_size_naissance 1.081086
## taille_groupe_naissance 1.568386
## presence_male_sub 1.508214
## summer_index_2_dt 2.449782
## spring_index_2_dt 1.327067
## winter_index_1_y 2.111976
## winter_index_2_y 3.388778
## presence_male_sub:winter_index_1_y 1.863250
## presence_male_sub:winter_index_2_y 2.213100All is good !
Goodness of fit
Avec le package performance, je compute le \(R^2\) marginal (variance expliquée par les
effets fixes) et le \(R^2\)
conditionnel (effets fixes + effets aléatoires).
Masse
## # R2 for Mixed Models
##
## Conditional R2: 0.691
## Marginal R2: 0.294
Tibia
## # R2 for Mixed Models
##
## Conditional R2: 0.589
## Marginal R2: 0.333
PLOT CONDITIONAL EFFECTS OF PREDICTORS
Je compute les valeurs d’interactions entre la présence de
helper et snow cover (axe 1 de l’ACP hiver) /
température hivernale (axe 2 de l’ACP hiver). Pour
cela, j’utilise le package marginaleffects qui calcule la
valeur d’interaction pour chaque niveau de
| Trait | Climatic variable | Helper | Estimate | Std. Error | t | p |
|---|---|---|---|---|---|---|
| Mass | Snow cover | 0 | -9.344 | 3.233 | -2.890 | 0.004 |
| 1 | -19.364 | 3.847 | -5.034 | 0.000 | ||
| Winter T° | 0 | 24.761 | 3.457 | 7.163 | 0.000 | |
| 1 | 16.081 | 3.365 | 4.778 | 0.000 | ||
| Tibia | Snow cover | 0 | -0.564 | 0.171 | -3.291 | 0.001 |
| 1 | -1.109 | 0.203 | -5.473 | 0.000 | ||
| Winter T° | 0 | 0.841 | 0.182 | 4.624 | 0.000 | |
| 1 | 0.495 | 0.178 | 2.782 | 0.005 | ||
| NOTE: Bold means \(p \leq\) 0.05, italic means \(p \leq 0.1\) |
Interaction helpers x temperature x snow
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_masse
## Data: pheno_abiotic_mod
##
## REML criterion at convergence: 14319
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0459 -0.4774 -0.0233 0.5114 3.6817
##
## Random effects:
## Groups Name Variance Std.Dev.
## mother (Intercept) 5322 72.96
## Residual 3840 61.96
## Number of obs: 1272, groups: mother, 120
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 344.5062 8.9462 206.8864 38.508 < 2e-16 ***
## annee_emergence_num -0.8062 5.9757 318.2490 -0.135 0.892759
## I(annee_emergence_num^2) 17.7860 4.1037 741.4370 4.334 1.67e-05 ***
## capture_time 24.6929 2.6351 1250.6840 9.371 < 2e-16 ***
## sexem 11.0779 3.6642 1145.7713 3.023 0.002556 **
## inbreeding 11.2565 4.0265 844.7211 2.796 0.005298 **
## litter_size_naissance -44.7599 2.2064 1225.1606 -20.286 < 2e-16 ***
## taille_groupe_naissance -16.6833 2.9929 1250.6791 -5.574 3.04e-08 ***
## presence_male_sub1 8.4939 5.5300 1242.5306 1.536 0.124803
## spring_index_1_dt -8.1931 2.6095 1207.7458 -3.140 0.001732 **
## spring_index_2_dt 11.6456 3.0317 1217.0471 3.841 0.000129 ***
## summer_index_2_dt -10.3269 3.7306 1203.7587 -2.768 0.005724 **
## winter_index_1_y -6.8152 3.2670 1208.8844 -2.086 0.037181 *
## winter_index_2_y 21.6783 3.5795 1212.6863 6.056 1.86e-09 ***
## presence_male_sub1:winter_index_1_y -10.5633 4.8481 1218.9682 -2.179 0.029533 *
## presence_male_sub1:winter_index_2_y -9.4891 3.9819 1214.7986 -2.383 0.017323 *
## winter_index_1_y:winter_index_2_y -6.1771 3.0478 1230.4961 -2.027 0.042905 *
## presence_male_sub1:winter_index_1_y:winter_index_2_y -10.2576 4.0390 1212.3681 -2.540 0.011221 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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} \]
Le code de la fonction est disponible ci-dessous, en cliquant sur
le bouton Show ->
ANODEV <- function(data,
response,
fixed_effects,
random_effects,
clim_vars,
time_var) {
# STOCKAGE RESULTATS ANODEV
results <- data.frame(
ClimVar = character(),
Dev = numeric(),
R2dev = numeric(),
F_stat = numeric(),
p = numeric(),
n = numeric(),
J = numeric(),
stringsAsFactors = FALSE
)
# STOCKAGE MODELE FULL TIME & CONSTANT
results_dev <- data.frame(
model = character(),
dev = numeric(),
k = numeric(),
stringsAsFactors = FALSE
)
# MODELE CONSTANT (Fcst)
form_const <- as.formula(
paste(response, "~",
paste(fixed_effects, collapse = " + "),
"+", paste(random_effects, collapse = " + "))
)
m0 <- lmerTest::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 <- lmerTest::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 <- lmerTest::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
n <- dft - df0
# F
F_stat <- ((D0 - Dc)/(J)) /
((Dc - Dt)/(n - J))
# p-value
p_val <- pf(F_stat,
df1 = J,
df2 = n - J,
lower.tail = FALSE)
# R2 dev
R2dev <- (D0 - Dc)/(D0 - Dt)
# Stockage
results <- rbind(results, data.frame(
ClimVar = clim,
Dev = Dc,
R2dev = round(R2dev, 2),
F_stat = round(F_stat, 2),
p = round(p_val, 3),
n = n,
J = J
))
}
# MODELE AVEC TOUTES LES VARIABLES CLIMATIQUES
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 FULL TIME DEPENDANT
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
n <- dft - df0
# F total
F_stat <- ((D0 - Dc)/(J)) /
((Dc - Dt)/(n - J))
# p total
p_val <- pf(F_stat,
df1 = J,
df2 = n - J,
lower.tail = FALSE)
# R2 total
R2dev <- (D0 - Dc)/(D0 - Dt)
# Ajout ligne TOTAL
results <- rbind(results, data.frame(
ClimVar = "TOTAL",
Dev = Dc,
R2dev = round(R2dev, 2),
F_stat = round(F_stat, 2),
p = round(p_val, 3),
n = n,
J = J
))
# MODELS
results_dev[1:3, "model"] <- c("D0", "Dc", "Dt")
results_dev[1:3, "dev"] <- c(D0, Dc, Dt)
results_dev[1:3, "k"] <- c(df0, dfc, dft)
return(list(results = results, results_dev = results_dev))
}
# ANODEV <- function(data,
# response,
# fixed_effects,
# random_effects,
# clim_vars,
# time_var) {
#
# results <- data.frame(
# ClimVar = character(),
# R2dev = numeric(),
# F_stat = numeric(),
# p = numeric(),
# n = numeric(),
# J = numeric(),
# stringsAsFactors = FALSE
# )
#
# # -----------------------------
# # MODELE CONSTANT
# # -----------------------------
# form_const <- as.formula(
# paste(response, "~",
# paste(fixed_effects, collapse = " + "),
# "+", paste(random_effects, collapse = " + "))
# )
#
# m0 <- lme4::lmer(form_const, data = data, REML = FALSE)
#
# # -----------------------------
# # MODELE CLIMAT COMPLET
# # -----------------------------
# clim_all <- paste(clim_vars, collapse = " + ")
#
# form_clim_all <- as.formula(
# paste(response, "~",
# paste(fixed_effects, collapse = " + "),
# "+", clim_all,
# "+", paste(random_effects, collapse = " + "))
# )
#
# m_all <- lme4::lmer(form_clim_all, data = data, REML = FALSE)
#
# # -----------------------------
# # MODELE TEMPS
# # -----------------------------
# 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 globales
# D0 <- deviance(m0)
# Dall <- deviance(m_all)
# Dt <- deviance(mt)
#
# df0 <- attr(logLik(m0), "df")
# dfall <- attr(logLik(m_all), "df")
# dft <- attr(logLik(mt), "df")
#
# n <- dft - df0
#
# # ============================================================
# # DROP ONE CLIMATIC VARIABLE
# # ============================================================
#
# for (clim in clim_vars) {
#
# # toutes sauf la variable testée
# remaining <- clim_vars[clim_vars != clim]
# clim_minus <- paste(remaining, collapse = " + ")
#
# form_minus <- as.formula(
# paste(response, "~",
# paste(fixed_effects, collapse = " + "),
# "+", clim_minus,
# "+", paste(random_effects, collapse = " + "))
# )
#
# m_minus <- lme4::lmer(form_minus, data = data, REML = FALSE)
#
# Dminus <- deviance(m_minus)
# dfminus <- attr(logLik(m_minus), "df")
#
# J <- dfall - dfminus
#
# # F statistic (Skalski)
# F_stat <- ((Dminus - Dall)/J) /
# ((Dall - Dt)/(n - (dfall - df0)))
#
# p_val <- pf(F_stat,
# df1 = J,
# df2 = n - (dfall - df0),
# lower.tail = FALSE)
#
# R2dev <- (Dminus - Dall)/(D0 - Dt)
#
# results <- rbind(results, data.frame(
# ClimVar = clim,
# R2dev = round(R2dev, 3),
# F_stat = round(F_stat, 3),
# p = round(p_val, 4),
# n = n,
# J = J
# ))
# }
#
# # ============================================================
# # SIGNAL CLIMATIQUE TOTAL
# # ============================================================
#
# Jtot <- dfall - df0
#
# F_tot <- ((D0 - Dall)/Jtot) /
# ((Dall - Dt)/(n - Jtot))
#
# p_tot <- pf(F_tot,
# df1 = Jtot,
# df2 = n - Jtot,
# lower.tail = FALSE)
#
# R2tot <- (D0 - Dall)/(D0 - Dt)
#
# results <- rbind(results, data.frame(
# ClimVar = "TOTAL_CLIMATE",
# R2dev = round(R2tot, 3),
# F_stat = round(F_tot, 3),
# p = round(p_tot, 4),
# n = n,
# J = Jtot
# ))
#
# return(results)
# }
ANODEV : Masse
J’inclus les tendances temporelles dans la partie fixe du modèle en plus des variables climatiques comme recommandé par Grosbois et al. (2008).
| Variables climatiques (+trend) | Dev | \(R^2\) | \(F_{cst/co/t}\) | p | n | J |
|---|---|---|---|---|---|---|
| (sub × winter snow cover) | 14570.91 | 0.06 | 0.75 | 0.483 | 25 | 2 |
| (sub × winter T°) | 14503.14 | 0.30 | 4.91 | 0.017 | 25 | 2 |
| summer NDVI | 14560.31 | 0.10 | 2.62 | 0.118 | 25 | 1 |
| spring aridity | 14545.81 | 0.15 | 4.22 | 0.051 | 25 | 1 |
| spring NDVI | 14560.27 | 0.10 | 2.63 | 0.118 | 25 | 1 |
| TOTAL | 14445.04 | 0.50 | 2.61 | 0.048 | 25 | 7 |
| NOTE: Bold means \(p \leq\) 0.05, italic means \(p \leq 0.1\) |
ANODEV : Tibia
| Variables climatiques (+trend) | Dev | \(R^2\) | \(F_{cst/co/t}\) | p | n | J |
|---|---|---|---|---|---|---|
| summer NDVI | 7041.758 | 0.03 | 0.77 | 0.388 | 24 | 1 |
| (sub × winter snow cover) | 7038.882 | 0.04 | 0.46 | 0.638 | 24 | 2 |
| (sub × winter T°) | 7018.393 | 0.09 | 1.12 | 0.343 | 24 | 2 |
| spring aridity | 7042.064 | 0.03 | 0.76 | 0.394 | 24 | 1 |
| spring NDVI | 6983.149 | 0.18 | 5.15 | 0.033 | 24 | 1 |
| TOTAL | 6922.948 | 0.34 | 1.24 | 0.336 | 24 | 7 |
| NOTE: Bold means \(p \leq\) 0.05, italic means \(p \leq 0.1\) |
MODELE NAO
On regarde par trait l’effet du NAO d’hiver pour justifier du rôle du changmeent climatique dans les variations des traits. Je compare ici l’AIC de 3 modèles par traits : (1) Sans NAO, (2) Effet principal NAO et (3) Interaction NAO x helpers.
Mass
| \(\mu\) | year | year\(^2\) | Capture time | Inbr. | Litter size | Group size | HELP | Sex | NAO | (HELP × NAO) | AICc | df | \(\Delta\) | w |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 349.56 | -8.45 | 14.88 | 20.53 | 8.39 | -42.98 | -16.93 | + | + | 16.31 | + | 14575.51 | 13 | 0 | 0.84 |
| 350.8 | -9.15 | 14.05 | 20.19 | 7.9 | -42.82 | -16.59 | + | + | 11.06 | - | 14578.83 | 12 | 3.32 | 0.16 |
| 353.42 | -8.34 | 13.23 | 19.32 | 6.03 | -40.87 | -15.46 | + | + | - | - | 14600.08 | 11 | 24.57 | 0 |
NOTE: Italic
variables represent de-trended climatic predictors, bold correspond to
the best parsimonious model, + means interaction or factor
has been added
|
TIBIA
| \(\mu\) | year | year\(^2\) | Capture time | Inbr. | Litter size | Group size | HELP | Sex | NAO | (HELP × NAO) | AICc | df | \(\Delta\) | w |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 49.82 | -0.58 | 1.73 | 0.47 | 0.15 | -1.13 | -0.98 | + | + | 0.73 | + | 7061.59 | 13 | 0 | 0.98 |
| 49.9 | -0.61 | 1.67 | 0.45 | 0.11 | -1.12 | -0.96 | + | + | 0.37 | - | 7069.02 | 12 | 7.43 | 0.02 |
| 49.98 | -0.57 | 1.65 | 0.42 | 0.06 | -1.05 | -0.92 | + | + | - | - | 7076.66 | 11 | 15.07 | 0 |
NOTE: Italic
variables represent de-trended climatic predictors, bold correspond to
the best parsimonious model, + means interaction or factor
has been added
|
| Trait | Climatic variable | Helper | Estimate | Std. Error | t | p |
|---|---|---|---|---|---|---|
| Mass | Winter NAO | 0 | 16.315 | 3.206 | 5.089 | 0.000 |
| 1 | 5.824 | 3.208 | 1.815 | 0.069 | ||
| Tibia | 0 | 0.733 | 0.167 | 4.392 | 0.000 | |
| 1 | 0.001 | 0.170 | 0.004 | 0.997 | ||
| NOTE: Bold means \(p \leq\) 0.05, italic means \(p \leq 0.1\) |
