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 :

1. de-trend les indices/variables climatiques ayant une tendance temporelle, pour éviter toute corrélation fallacieuse entre \(y\) et la variable climatique
1. réaliser une sélection de modèle basée sur le critère AIC pour identifier les variables climatiques d’intérêt
1. 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

# CLIMATIC VARIABLES
ClimVars <- c(
  "summer_index",
  "summer_precipitation_y",
  "spring_index",
  "spring_ndvi_y",
  "winter_index",
  "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] "quadratic"
## 
## $summer_precipitation_y
## [1] "quadratic"
## 
## $spring_index
## [1] "linear"
## 
## $spring_ndvi_y
## [1] "linear"
## 
## $winter_index
## [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.

# 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_2_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. 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
    winter_index_y +
    winter_index_2_y +
    summer_index_dt +
    summer_precipitation_y_dt +
    spring_index_dt +
    spring_ndvi_y_dt +
   
   # SOCIAL x ENVIRONMENT
   winter_index_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
)

Ci-dessous la liste de l’ensemble des modèles avec un \(\Delta AICc < 2\) :

Model selection for body mass (models with ΔAICc ≤ 2)
Fixed part									Summer		Spring		Winter
\(\mu\)	year	year\(^2\)	Capture time	Inbr.	Litter size	Group size	Male sub	Sex	Su. index	Su. prec.	Sp. index	Sp. NDVI	Wi. index	Wi. index2	(Sub × index)	(Sub × index2)	AICc	df	\(\Delta\)	w
349.44	-2.15	13.26	24.27	12.11	-44.27	-21.12	+	+	-	-14.59	-12.87	9.64	-4.76	24.91	+	+	14461.74	18	0	0.37
350.97	-3.57	11.97	24.42	12	-43.54	-20.58	+	+	6.23	-	-11.96	10.21	-3.69	21.95	+	+	14463.64	18	1.89	0.15
349.64	-2.35	13.07	24.33	12.16	-44.15	-21.04	+	+	1.31	-12.21	-12.71	9.71	-4.6	24.54	+	+	14463.7	19	1.96	0.14
349.98	-2.56	12.58	23.82	11.78	-44.36	-20.89	+	+	-	-14.33	-12.53	9.91	-8.54	24.56	-	+	14463.91	17	2.16	0.13
351.45	-3.94	11.35	23.99	11.71	-43.64	-20.37	+	+	6.21	-	-11.64	10.45	-7.35	21.7	-	+	14465.48	17	3.74	0.06
350.23	-2.81	12.35	23.9	11.85	-44.21	-20.8	+	+	1.62	-11.39	-12.33	9.99	-8.32	24.11	-	+	14465.81	18	4.06	0.05
351.36	-2.5	11.9	24.3	11.13	-44.09	-20.74	+	+	-	-13.28	-12.99	10.49	-5.4	18.74	+	-	14466	17	4.25	0.04
351.63	-2.82	11.45	23.92	10.95	-44.19	-20.58	+	+	-	-13.18	-12.68	10.64	-8.54	19.04	-	-	14466.89	16	5.15	0.03
352.77	-3.8	10.71	24.41	11	-43.43	-20.24	+	+	5.52	-	-12.17	11.03	-4.42	15.98	+	-	14467.93	17	6.19	0.02
351.46	-2.6	11.8	24.33	11.16	-44.03	-20.7	+	+	0.63	-12.13	-12.91	10.53	-5.32	18.54	+	-	14468.03	18	6.29	0.02
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: 14356.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2288 -0.5090 -0.0078  0.5318  3.5846 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  mother   (Intercept) 5147     71.75   
##  Residual             3945     62.81   
## Number of obs: 1272, groups:  mother, 120
## 
## Fixed effects:
##                                     Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                          349.338      8.895  210.709  39.272  < 2e-16 ***
## annee_emergence_num                   -2.078      5.886  304.753  -0.353  0.72426    
## I(annee_emergence_num^2)              13.318      4.207  741.845   3.166  0.00161 ** 
## capture_time                          24.294      2.649 1252.232   9.172  < 2e-16 ***
## sexem                                 10.997      3.714 1149.832   2.961  0.00313 ** 
## inbreeding                            12.106      4.072  806.089   2.973  0.00304 ** 
## litter_size_naissance                -44.290      2.218 1230.839 -19.971  < 2e-16 ***
## taille_groupe_naissance              -21.145      2.934 1252.618  -7.208 9.79e-13 ***
## presence_male_sub1                    13.119      5.484 1247.091   2.392  0.01690 *  
## spring_index_dt                      -12.864      2.225 1208.059  -5.780 9.48e-09 ***
## spring_ndvi_y_dt                       9.642      2.360 1223.941   4.085 4.69e-05 ***
## summer_precipitation_y_dt            -14.606      4.496 1215.676  -3.249  0.00119 ** 
## winter_index_y                        -4.767      2.859 1214.335  -1.667  0.09572 .  
## winter_index_2_y                      24.941      3.631 1243.142   6.868 1.03e-11 ***
## presence_male_sub1:winter_index_y     -9.426      4.616 1225.929  -2.042  0.04135 *  
## presence_male_sub1:winter_index_2_y  -11.031      4.391 1228.482  -2.512  0.01213 *  
## ---
## 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_index_2_y +
    summer_index_dt +
    summer_precipitation_y_dt +
    spring_index_dt +
    spring_ndvi_y_dt +
   
   # SOCIAL x ENVIRONMENT
    winter_index_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\) :

Model selection for tibia length (models with ΔAICc ≤ 2)
Fixed part									Summer		Spring		Winter
\(\mu\)	year	year\(^2\)	Capture time	Inbr.	Litter size	Group size	Male sub	Sex	Su. index	Su. prec.	Sp. index	Sp. NDVI	Wi. index	Wi. index2	(Sub × index)	(Sub × index2)	AICc	df	\(\Delta\)	w
50.21	-0.49	1.29	0.62	0.23	-1.16	-1.22	+	+	-	-	-0.33	1.05	-0.11	0.7	+	+	6965	17	0	0.26
50.26	-0.48	1.25	0.63	0.2	-1.16	-1.21	+	+	-	-	-0.34	1.07	-0.14	0.49	+	-	6965.96	16	0.96	0.16
50.18	-0.46	1.31	0.63	0.24	-1.17	-1.23	+	+	-	-0.19	-0.34	1.03	-0.13	0.76	+	+	6966.43	18	1.43	0.13
50.2	-0.48	1.29	0.63	0.24	-1.16	-1.22	+	+	0.06	-	-0.33	1.04	-0.11	0.72	+	+	6966.79	18	1.79	0.11
50.24	-0.5	1.25	0.6	0.21	-1.17	-1.21	+	+	-	-	-0.31	1.06	-0.32	0.68	-	+	6967.59	16	2.58	0.07
50.28	-0.5	1.23	0.6	0.19	-1.16	-1.2	+	+	-	-	-0.32	1.08	-0.32	0.51	-	-	6967.65	15	2.64	0.07
50.24	-0.46	1.27	0.63	0.21	-1.16	-1.21	+	+	-	-0.14	-0.34	1.06	-0.15	0.53	+	-	6967.65	17	2.65	0.07
50.26	-0.48	1.26	0.63	0.2	-1.16	-1.21	+	+	0.03	-	-0.34	1.07	-0.14	0.5	+	-	6967.93	17	2.92	0.06
50.17	-0.46	1.32	0.63	0.24	-1.17	-1.23	+	+	-0.07	-0.31	-0.35	1.03	-0.14	0.78	+	+	6968.4	19	3.39	0.05
50.21	-0.48	1.27	0.61	0.22	-1.17	-1.21	+	+	-	-0.18	-0.32	1.05	-0.34	0.75	-	+	6969.1	17	4.1	0.03
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: 6956.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6686 -0.5176  0.0857  0.5940  2.9514 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  mother   (Intercept)  7.207   2.685   
##  Residual             11.556   3.399   
## Number of obs: 1272, groups:  mother, 120
## 
## Fixed effects:
##                                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                           50.2078     0.3847  241.2631 130.507  < 2e-16 ***
## annee_emergence_num                   -0.4860     0.2593  245.4327  -1.874 0.062072 .  
## I(annee_emergence_num^2)               1.2913     0.1993  476.6210   6.478 2.31e-10 ***
## capture_time                           0.6257     0.1367 1165.7503   4.576 5.25e-06 ***
## sexem                                  1.0808     0.2002 1173.8918   5.399 8.10e-08 ***
## inbreeding                             0.2293     0.1948  511.1105   1.177 0.239684    
## litter_size_naissance                 -1.1602     0.1169 1256.8777  -9.928  < 2e-16 ***
## taille_groupe_naissance               -1.2173     0.1522 1222.0561  -7.997 2.94e-15 ***
## presence_male_sub1                     0.9946     0.2877 1254.3739   3.457 0.000565 ***
## spring_index_dt                       -0.3298     0.1182 1242.2625  -2.790 0.005352 ** 
## spring_ndvi_y_dt                       1.0462     0.1239 1255.9042   8.445  < 2e-16 ***
## winter_index_y                        -0.1087     0.1502 1253.4218  -0.724 0.469303    
## winter_index_2_y                       0.6983     0.1713 1256.9559   4.077 4.85e-05 ***
## presence_male_sub1:winter_index_y     -0.5245     0.2445 1256.6080  -2.145 0.032107 *  
## presence_male_sub1:winter_index_2_y   -0.4035     0.2313 1255.1348  -1.744 0.081341 .  
## ---
## 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.319712
## I(annee_emergence_num^2)            1.291476
## capture_time                        1.050778
## sexe                                1.008672
## inbreeding                          1.092357
## litter_size_naissance               1.080251
## taille_groupe_naissance             1.699966
## presence_male_sub                   1.525481
## spring_index_dt                     1.194894
## spring_ndvi_y_dt                    1.150234
## summer_precipitation_y_dt           1.598148
## winter_index_y                      1.881247
## winter_index_2_y                    3.110073
## presence_male_sub:winter_index_y    1.891088
## presence_male_sub:winter_index_2_y  2.279059

Tibia

##                                    vif_score
## annee_emergence_num                 1.339880
## I(annee_emergence_num^2)            1.304284
## capture_time                        1.046082
## sexe                                1.008444
## inbreeding                          1.078719
## litter_size_naissance               1.063978
## taille_groupe_naissance             1.629416
## presence_male_sub                   1.492879
## spring_index_dt                     1.188028
## spring_ndvi_y_dt                    1.133007
## winter_index_y                      1.829132
## winter_index_2_y                    2.437391
## presence_male_sub:winter_index_y    1.872380
## presence_male_sub:winter_index_2_y  2.210807

All 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 = "lm", 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.

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

performance::r2(best_model_masse)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.694
##      Marginal R2: 0.295

Tibia

performance::r2(best_model_tibia)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.584
##      Marginal R2: 0.325

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

INTERACTION VALUE BETWEEN PRESENCE OF HELPER AND WINTER CONDITION
Trait	Climatic variable	Helper	Estimate	Std. Error	t	p
Mass	Snow cover	0	-4.767	2.859	-1.667	0.095
	Snow cover	1	-14.193	3.523	-4.029	0.000
	Winter T°	0	24.941	3.631	6.868	0.000
	Winter T°	1	13.910	3.309	4.204	0.000
Tibia	Snow cover	0	-0.109	0.150	-0.725	0.468
	Snow cover	1	-0.633	0.185	-3.421	0.001
	Winter T°	0	0.698	0.171	4.077	0.000
	Winter T°	1	0.295	0.165	1.784	0.074
NOTE: Bold means \(p \leq\) 0.05, italic means \(p \leq 0.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

\[ 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).

ANODEV MASS WITH TREND COVARIATES-SPECIFIC
Variables climatiques (+trend)	Dev	\(R^2\)	\(F_{cst/co/t}\)	p	n	J
(sub × winter snow cover)	14558.96	0.07	0.82	0.456	24	2
(sub × winter T°)	14496.25	0.30	4.66	0.020	24	2
summer precipitation	14569.88	0.03	0.69	0.415	24	1
spring index	14526.42	0.19	5.31	0.031	24	1
spring NDVI	14557.32	0.07	1.86	0.185	24	1
TOTAL	14425.20	0.56	3.05	0.028	24	7
NOTE: Bold means \(p \leq\) 0.05, italic means \(p \leq 0.1\)

ANODEV : Tibia

ANODEV TIBIA WITH TREND COVARIATES-SPECIFIC
Variables climatiques (+trend)	Dev	\(R^2\)	\(F_{cst/co/t}\)	p	n	J
(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 index	7042.064	0.03	0.76	0.394	24	1
spring NDVI	6982.125	0.19	5.25	0.031	24	1
TOTAL	6930.516	0.32	1.40	0.268	24	6
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

NAO model selection for body mass (models with ΔAICc ≤ 2)
\(\mu\)	year	year\(^2\)	Capture time	Inbr.	Litter size	Group size	HELP	Sex	NAO	(HELP × NAO)	AICc	df	\(\Delta\)	w
350.33	-8.45	14.88	20.53	8.39	-42.98	-16.93	+	+	16.71	+	14575.51	13	0	0.84
351.32	-9.15	14.05	20.19	7.9	-42.82	-16.59	+	+	11.33	-	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

NAO model selection for body mass (models with ΔAICc ≤ 2)
\(\mu\)	year	year\(^2\)	Capture time	Inbr.	Litter size	Group size	HELP	Sex	NAO	(HELP × NAO)	AICc	df	\(\Delta\)	w
49.86	-0.58	1.73	0.47	0.15	-1.13	-0.98	+	+	0.75	+	7061.59	13	0	0.98
49.92	-0.61	1.67	0.45	0.11	-1.12	-0.96	+	+	0.38	-	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

Sélection de modèles avec données climatiques

Pierre-Alexandre QUITTET

2026-01-22

DEMARCHE D’ANALYSE

INITIALISATION

DE-TRENDED VARIABLES

Tester les tendances temporelles

Visualisation des tendances temporelles

Construction des variables de-trended

Analyse des résidus

MODEL SELECTION

Extract best model per trait

Data preparation

Fixed effects selection with `MuMIn::dredge()`

Mass

Tibia

MODEL DIAGNOSIS

Residuals diagnosis

Mass

Tibia

Multicollinearity

Mass

Tibia

Partial residuals

Mass

Tibia

Goodness of fit

PLOT CONDITIONAL EFFECTS OF PREDICTORS

ANODEV

ANODEV : Masse

ANODEV : Tibia

MODELE NAO

Mass

TIBIA

Sélection de modèles avec données climatiques

Pierre-Alexandre QUITTET

2026-01-22

DEMARCHE D’ANALYSE

INITIALISATION

DE-TRENDED VARIABLES

Tester les tendances temporelles

Visualisation des tendances temporelles

Construction des variables de-trended

Analyse des résidus

MODEL SELECTION

Extract best model per trait

Data preparation

Fixed effects selection with MuMIn::dredge()

Mass

Tibia

MODEL DIAGNOSIS

Residuals diagnosis

Mass

Tibia

Multicollinearity

Mass

Tibia

Partial residuals

Mass

Tibia

Goodness of fit

PLOT CONDITIONAL EFFECTS OF PREDICTORS

ANODEV

ANODEV : Masse

ANODEV : Tibia

MODELE NAO

Mass

TIBIA

Fixed effects selection with `MuMIn::dredge()`