iteration = 1000 # nombre d'iteration pour les bootstraps 
options(digits=2)  # precision of decimal numbers
set.seed(123)  # set random

1 VARIATION TEMPORELLE DES VARIABLES SOCIALES

On trie le jeu de données de façon à éviter la pseudoréplication i.e. deux observations d’individus sur un même territoire et une même année sont identiques en terme de structure sociale. On ne conserve donc que les territoire x année unique.

##   annee_capture territoire taille_groupe litter_size sex_ratio_subordinates
## 1          1990          b             8           3                   1.00
## 2          1991          b            10           1                   0.80
## 4          1990          a             7           0                   1.00
## 5          1993          j             9           5                   0.33
## 7          2006          n             4           3                     NA
## 8          2006          d            NA          NA                     NA
##   nb_male_subordinates presence_male_subordinates
## 1                    1                          1
## 2                    4                          1
## 4                    3                          1
## 5                    1                          1
## 7                    0                          0
## 8                   NA                         NA

Effectif par année

Le nombre de territoires échantillonnés augmentent au cours du temps, avec de faibles / très faibles effectifs dans les premières années.

Méthode de mesure de la variation temporelle

On mesure maintenant la variation temporelle des variables sociales à partir de ce jeu de données. On calcule la moyenne par année $\pm$ IC95%.

Les variations temporelles des variables sociales suivantes ont été explorées :

Taille du groupe [continue]
Taille de portée [continue]
Nombre de helpers [continue]
Présence de helpers [binaire]

Données manquantes par variables sociales

## GROUP SIZE
##    2    3    4    5    6    7    8    9   10   11   12   13   14   15   17 <NA> 
##  152   93   93   74   69   50   45   28   18    4    5    2    2    1    1  144

## LITTER SIZE
##    0    1    2    3    4    5    6    7    8 <NA> 
##  213   30   43   98  149   80   17    6    1  144

## HELPERS
##    0    1    2    3    4    5    6 <NA> 
##  350  144   74   36   15    4    4  154

1.1 Taille de groupe et taille de portée

NOTE : La taille de groupe n’inclut pas les marmottons.

Au dessus de la moyenne En dessous de la moyenne

1.2 Helpers

Le nombre d’helpers à $t$ au territoire $i$ a été calculé comme suit : \[ \text{HELP}_{i,t} = (n_{mâle} -1)_{i,t} \] Avec $n_{mâle}$ le nombre d’individu mâle d’âge $\geq 2$ ans (n’inclut donc pas les yearlings) (Allainé et al., 2000). Dans la formule, le $-1$ sert à retirer le mâle dominant du calcul.

Au dessus de la moyenne En dessous de la moyenne

1.3 Corrélation entre les variables sociales

##             Group Litter Nb Helpers Helpers bin
## Group                 ns        ***         ***
## Litter       0.04                ns          ns
## Nb Helpers   0.61  -0.06                    ***
## Helpers bin  0.56  -0.04       0.95

La taille de groupe est corrélée positivement au nombre de helpers. Contrairement à ce que j’avais imaginé, la taille de portée n’est pas du tout corrélée au nombre de helpers, ni à la taille du groupe. J’aurais en effet imaginé que les grands groupes avec beaucoup de helpers auraient des tailles de portée plus grandes, les mères ayant des contextes sociaux plus favorables.

2 PROXY DE FITNESS

2.1 Proxy de fitness et scaling

Proxy de fitness

On découpe la fitness totale en 5 épisodes de sélection :

Survie au premier hiver ($W_1$)
Survie jusqu’à 2 ans ($W_2$)
Survie jusqu’à 3 ans ($W_3$)
Accès à la dominance ($W_4$)
Durée de la dominance ($W_5$)
LRS ($W$)

Ces paramètres (à part le LRS) ont été estimé par CMR avec l’algorithme VITERBI.

Scaling des traits
Dans les parties suivantes, la masse et la taille des marmottons ont été centré réduit par classe d’âge en jour (de 0 à 15 jours). Ensuite, une seconde normalisation a centré réduit la masse et la taille sur l’ensemble des données (i.e. population), de façon à pouvoir à comparer les deux traits entre eux.

2.2 Distribution de la fitness absolue

2.3 Compute de la fitness relative

On compute la fitness individuelle (reconstruite avec VITERBI) qu’on va diviser par la moyenne de fitness absolue pour chaque année pour obtenir la fitness relative individuelle. La fitness relative $w$ est donc définie pour l’individu $i$ au $k^{th}$ épisode de sélection pour l’année $t$ tel que :

\[ w_{i,k,t} = \frac{W_{i,k, t}}{\bar{W}_{k,t}} \]

Le problème ici c’est qu’avec l’algorithme Viterbi, les proxy de fitness ne sont plus binaires. C’est donc moins direct pour calculer les moyennes conditionnellement aux épisodes de sélection précédent. En effet, dans le cas d’un proxy binaire (e.g. survie ‘classique’), il suffit pour l’épisode de sélection 2, de ne considérer que le survivants de l’épisode de sélection 1. Or ici, les estimations de fitness $W$ sont continues. On doit donc calculer la moyenne de fitness par années à l’aide de la formule de A. Arnold (2023).

\[ \bar{W}(z) = \int p(z) \ W(z) \ dz \]

Ici $p(z)$ est la distribution des fréquences du trait pré-sélection. On doit prendre en compte le fait que cette distribution va changer après chaque épisode $k$ de sélection pour devenir $p^*(z)$. La distribution des fréquences phénotypique $p^*(z)$ après l’épisode de sélection $k_1$ peut être obtenue en multipliant $p(z)$ par la fitness relative à l’épisode $k-1$ : $w_{k=1}$ (Arnold, 2023).

On peut obtenir la distribution post-sélection $p^*(z)$ en pondèrant $z$ par la fitness relative $w(z)$ associée, soit :

\[ p^*(z) = \int \frac{W(z) \cdot p(z)}{\bar{W}} = p(z) \cdot w(z) \]

Cette pondération est parfaitement équivalente à calculer la moyenne pondérée de la fitness où les poids de chaque phénotype correspond à leur fitness relative.

# RELATIVE FITNESS
# w1
data_fitness <- data_fitness |> 
  group_by(annee_emergence) |> 
  # filter(annee_emergence != 2023) |> 
  mutate(w1 = W1 / mean(W1, na.rm = T)) |> 
  as.data.frame()

data_fitness$p1 <- data_fitness$w1  # relative proportion of phenotype $z$ after episode 1

# w2
data_fitness <- data_fitness |> 
  group_by(as.numeric(as.character(annee_emergence))) |>  
  filter(!is.na(p1)) |>
  mutate(w2 = W2 / (sum(W2 * proportions(p1)))) |> 
  as.data.frame()

data_fitness$p2 <- data_fitness$w1 * data_fitness$w2   # relative proportion of phenotype $z$ after episode 1 and episode 2

# w3
data_fitness <- data_fitness |>
  group_by(as.numeric(as.character(annee_emergence))) |>  #
  filter(!is.na(p2)) |>
  mutate(w3 = W3 / sum(W3 * proportions(p2))) |> 
  as.data.frame()

data_fitness$p3 <- data_fitness$w1 * data_fitness$w2 *  data_fitness$w3

# w4
data_fitness <- data_fitness |> 
  group_by(annee_emergence) |> 
  filter(!is.na(p3)) |>
  mutate(w4 = W4 / sum(W4 * proportions(p3))) |> 
  as.data.frame()

data_fitness$p4 <- data_fitness$w1 * data_fitness$w2 *  data_fitness$w3 * data_fitness$w4

# w5
data_fitness <- data_fitness |> 
  group_by(annee_emergence) |> 
  filter(!is.na(p4)) |>
  mutate(w5 = W5 / sum(W5 * proportions(p4))) |> 
  as.data.frame()

# w total
# w total as durée de dominance
data_fitness <- data_fitness |> 
  group_by(annee_emergence) |> 
  mutate(w = W5 / mean(W5, na.rm = T)) |> 
  as.data.frame()

# as LRS
data_fitness <- data_fitness |> 
  group_by(annee_emergence) |> 
  mutate(LRS_relative = LRS / mean(LRS, na.rm = T)) |> 
  as.data.frame()

2.4 Corrélation entre les proxy de fitness

🧠 Hypothèse
On s’attend à ce que tous les proxy de fitness relatives soient corrélés positivement entre eux.

Les données n’étant pas distribuées normalement, j’ai utilisé la corrélation de Spearman pour calculer l’associaton entre chaque proxy de fitness. A noter que je n’ai pas intégré $w5$ la durée de la dominance car par construction, elle est corrélée à 1 à $w4$ l’accès à la dominance.

##                w1   w2   w3   w4 LRS_relative
## w1                 ***  ***  ***          ***
## w2           0.79       ***  ***          ***
## w3           0.63 0.81       ***          ***
## w4           0.53 0.66 0.81               ***
## LRS_relative 0.32 0.41 0.54 0.65

Tous les proxy de fitness relative sont significativement corrélés positivement entre eux. Donc un individu ayant une bonne / mauvaise fitness relative dans un proxy donné par rapport à ses conspécifiques, aura une bonne / mauvaise fintess relative dans tous les autres proxy.

3 SELECTION GLOBALE

🧠 Hypothèse
Sélection positive et directionnelle sur la masse et la taille sur tous les épisodes de sélection.

3.1 Fonctions

Ci-dessous les fonctions ayant permis de générer les résultats et les plots de différentiel / gradient de sélection de la partie “SELECTION GLOBALE”.

3.1.1 compute_selection()

Fonction permettant de compute les différentiels et les gradients de sélection par trait. Le code est disponible en cliquant sur le bouton Show.

# ──────────────────────────────────────────────────────────────
# Fonction pour calculer les gradients de sélection linéaire et quadratique
# par bootstrap, pour un ou plusieurs traits
# ──────────────────────────────────────────────────────────────

compute_selection <- function(data,
                              trait,               # Nom(s) de(s) trait(s) à tester
                              fitness_proxy,       # Variable de fitness à utiliser ("LRS_relative", "w1", etc.)
                              iteration = 1000,    # Nombre de tirages bootstrap
                              weights = NULL) {    # Variable de pondération (optionnelle)
  
  # Supprimer les lignes avec NA dans les colonnes des traits
  for(trait.temp in trait){
    data <- data[!is.na(data[[trait.temp]]), ]
  }
  
  # Initialiser les matrices pour stocker les coefficients bootstrapés
  cov_boot_linear <- as.data.frame(matrix(
    NA, ncol = length(trait), nrow = iteration,
    dimnames = list(NULL, trait)
  ))
  
  cov_boot_non_linear <- as.data.frame(matrix(
    NA, ncol = length(trait), nrow = iteration,
    dimnames = list(NULL, trait)
  ))
  
  # Formule du modèle incluant les effets quadratiques
  formula_mod <- as.formula(paste0(
    "fitness ~ ",
    paste0(trait, collapse = " + "),
    paste0(" + I(", trait, "^2)", collapse = " + ")
  ))
  
  # Extraire le numéro associé à la fitness 
  k <- gsub(pattern = "w", replacement = "", x = fitness_proxy)
  
  # Identifier chaque ligne pour le bootstrap
  data$row_id <- seq_len(nrow(data))
  
  # BOUCLE BOOTSTRAP
  for (i in 1:iteration) {
    
    # Tirage avec remise, stratifié par année et territoire
    index_random <- data |> 
      group_by(annee_emergence, territoire) |> 
      reframe(row_id = sample(row_id, size = n(), replace = TRUE)) |> 
      pull(row_id)
    
    data_boot_i <- data[index_random, ]
    
    # Sélection de la variable de fitness absolue
    if(fitness_proxy == "LRS_relative"){
      fitness_abs <- sym("LRS")
    } else {
      fitness_abs <- sym(paste0("W", k))
    }

    # ON RECALCULE LA FITNESS RELATIVE PAR ITERATION (car la mean(w_t) doit être égale à 1)
    if (!is.null(weights)) {
      weight_var <- sym(weights)
      data_boot_i <- data_boot_i |>
        group_by(annee_emergence) |>
        filter(!is.na(!!weight_var) & !is.na(!!fitness_abs)) |>
        mutate(fitness = !!fitness_abs / weighted.mean(!!fitness_abs, w = !!weight_var)) |>
        ungroup()
      weights_mod <- proportions(data_boot_i[[weights]])
    } else {
      data_boot_i <- data_boot_i |>
        group_by(annee_emergence) |>
        filter(!is.na(!!fitness_abs)) |>
        mutate(fitness = !!fitness_abs / mean(!!fitness_abs)) |>
        ungroup()
      weights_mod <- NULL
    }
    
    # Ajustement du modèle linéaire
    mod_temp <- lm(
      formula = formula_mod,
      data = data_boot_i,
      weights = weights_mod
    )
    
    # Extraction des coefficients associés aux effets linéaires et quadratiques
    coef_temp <- coef(mod_temp)
    cov_boot_linear[i, ] <- coef_temp[names(coef_temp) %in% trait]
    cov_boot_non_linear[i, ] <- coef_temp[names(coef_temp) %in% paste0("I(", trait, "^2)")]
  }
  
  # Moyennes et intervalles de confiance des coefficients
  mean_boot_linear <- apply(cov_boot_linear, 2, mean, na.rm = TRUE)
  ci_lower_linear <- apply(cov_boot_linear, 2, quantile, probs = 0.025, na.rm = TRUE)
  ci_upper_linear <- apply(cov_boot_linear, 2, quantile, probs = 0.975, na.rm = TRUE)

  mean_boot_non_linear <- apply(cov_boot_non_linear, 2, mean, na.rm = TRUE)
  ci_lower_non_linear <- apply(cov_boot_non_linear, 2, quantile, probs = 0.025, na.rm = TRUE)
  ci_upper_non_linear <- apply(cov_boot_non_linear, 2, quantile, probs = 0.975, na.rm = TRUE)
  
  
  # Fonction d’affichage des niveaux de significativité
  p_values_function <- function(p) {
    ifelse(
      p < 0.001 | p > 0.999, "***",
      ifelse(p < 0.01 | p > 0.99, "**",
             ifelse(p < 0.05 | p > 0.95, "*",
                    ifelse(p < 0.1 | p > 0.9, ".", "")))
    )
  }

  # Calcul des "p-values" bootstrap 
    # values
    p_linear_values <- apply(cov_boot_linear, 2, function(x) {
    prop_neg <- min(sum(x < 0, na.rm = TRUE), sum(x > 0, na.rm = TRUE)) / iteration
  })
  
  p_non_linear_values <- apply(cov_boot_non_linear, 2, function(x) {
    prop_neg <- min(sum(x < 0, na.rm = TRUE), sum(x > 0, na.rm = TRUE)) / iteration
  })
  
  
  # symbol
  p_linear <- apply(cov_boot_linear, 2, function(x) {
    prop_neg <- sum(x < 0, na.rm = TRUE) / iteration
    p_values_function(prop_neg)
  })
  
  p_non_linear <- apply(cov_boot_non_linear, 2, function(x) {
    prop_neg <- sum(x < 0, na.rm = TRUE) / iteration
    p_values_function(prop_neg)
  })
  
  # Préparation du tableau de sortie pour visualisation
  df_gg <- data.frame(
    y = c(mean_boot_linear, mean_boot_non_linear),
    ci_lower = c(ci_lower_linear, ci_lower_non_linear),
    ci_upper = c(ci_upper_linear, ci_upper_non_linear),
    trait = trait,
    p_value = c(p_linear_values, p_non_linear_values),
    p = c(p_linear, p_non_linear),
    method = rep(c("Linear", "Quadratic"), each = length(trait))
  )
  
  return(df_gg)
}

# ──────────────────────────────────────────────────────────────
# Fonction pour itérer compute_selection() sur plusieurs traits
# ──────────────────────────────────────────────────────────────

bind_selection <- function(trait = c("masse_scaled", "tibia_scaled"),
                           data = data,
                           fitness_proxy = fitness_proxy,
                           iteration = 1000,
                           weights = NULL) {
  
  output <- c()  # initialiser le conteneur de résultats
  
  for (trait_temp in trait) {
    output <- rbind(output, compute_selection(
      data = data,
      fitness_proxy = fitness_proxy,
      iteration = iteration,
      weights = weights,
      trait = trait_temp
    ))
  }
  
  return(output)
}

3.1.2 plot_compute_selection()

Fonction permettant de représenter visuellement les différentiels et les gradients de sélection par trait. Le code est disponible en cliquant sur le bouton Show.

plot_compute_selection <- function(data,
                                   ylab,
                                   size_point = 4,
                                   size_es = 1.3,
                                   size_text = 8,
                                   base_size = 25) {
  
  trait_name <- c(
    `masse_scaled` = "Mass",
    `tibia_scaled` = "Tibia",
    `Linear` = "Linear",
    `Quadratic` = "Quadratic"
  )
  
    
  ggplot(data, aes(y = y, x = k, color = method, fill = method, group = trait, shape = method)) +
    scale_shape_manual(values = c("Linear" = 16, "Quadratic" = 15)) +
    facet_grid(method ~ trait, labeller = as_labeller(trait_name), scales = "free_y") +
    geom_hline(
      yintercept = 0,
      linetype = 2,
      alpha = 0.5,
      lwd = 1
    ) +
    geom_text(
      mapping = aes(x = k, y = ci_upper + 0.02, label = p),
      hjust = 0.5,
      size = size_text
    ) +
    
    #  geom_line(color = Fragman:::transp("black", 0.4)) +
    xlab("Episode of selection") +
    ylab(ylab) +
    theme_classic(base_size = base_size) +
    theme(panel.grid.major.y = element_line(
      color =  Fragman:::transp("#474747", 0.3),
      size = 0.5,
      linetype = 2
    )) +
    
    geom_errorbar(
      aes(ymin = ci_lower, ymax = ci_upper),
      width = 0,
      alpha = 0.5,
      lwd = size_es
    ) +
    geom_point(size = size_point) +
  
        # Légende
    scale_color_manual(name = "Modèle",
                       values = c("Linear" = "#27408B", "Quadratic" = "#E76F51")) +
    scale_fill_manual(name = "Intervalle de confiance",
                      values = c("Linear" = "#0000FF", "Quadratic" = "#E76F51")) +
    theme(legend.position = "none")
  
}

3.1.3 plot_overall_fitness_trait()

Fonction représentant visuellement la relation entre les traits et la fitness (différentiel et gradient de sélection) en décomposant la relation linéaire et quadratique.

# ──────────────────────────────────────────────────────────────
# Fonction : plot_overall_fitness_trait
# Description : Trace les courbes de fitness (linéaire et quadratique)
#               en fonction d’un trait donné, avec intervalles de confiance
#               et ajustement selon le type de distribution.
# ──────────────────────────────────────────────────────────────

plot_overall_fitness_trait <- function(fitness_proxy,
                                       trait = "masse_scaled",
                                       data,
                                       significance = 1,
                                       weights,
                                       trait_to_represent = "masse_scaled",
                                       ggtitle,
                                       ylab,
                                       family = "gaussian",
                                       xlab = "Scaled mass") {
  
  # ─── Préparation des données ──────────────────────────────
  
  for(trait_temp in trait){
    
    # Nettoyage des données (suppression des NA)
    data_filtered <- data[!is.na(data[[trait_temp]]), ]
    weights_filtered <- weights[!is.na(data[[trait_temp]])]
    
    # ─── Modélisation ─────────────────────────────────────────
    if(length(trait) == 1){  # univarié
      # Modèle linéaire directionnel
      model_lin <- lm(as.formula(paste(fitness_proxy, "~", trait_temp)),
                      data = data_filtered, weights = weights_filtered)
      
      # Modèle quadratique
      model_quad <- lm(as.formula(paste0(fitness_proxy, " ~ ", trait_temp, " + I(", trait_temp, "^2)")),
                       data = data_filtered, weights = weights_filtered)
      
      # Point pour prédiction à trait == 0
      point_df <- data.frame(x = 0)
      names(point_df)[1] <- trait_temp
      
    } else{  # multivariée
      # Modèle linéaire directionnel
      formula_mod_lin <- as.formula(paste0(
        paste0(fitness_proxy, "~ "),
        paste0(trait, collapse = " + ")
      ))
      model_lin <- lm(formula_mod_lin,
                      data = data_filtered, weights = weights_filtered)
      
      # Modèle quadratique
      formula_mod_quadra <- as.formula(paste0(
        paste0(fitness_proxy, "~ "),
        paste0(trait, collapse = " + "),
        paste0(" + I(", trait, "^2)", collapse = " + ")
      ))
      
      model_quad <- lm(formula_mod_quadra,
                       data = data_filtered, weights = weights_filtered)
      
      # Point pour prédiction à trait == 0
      point_df <- as.data.frame(matrix(0, nrow = 1, ncol = length(trait)))
colnames(point_df) <- trait
    }
   
    
    # Type de prédiction : "response" pour Gaussian, "link" sinon
    type_prediction <- ifelse(family %in% c("binomial", "poisson"), "link", "response")
    
    # ─── Prédictions au point 0 (pour illustrer le pic de fitness) ─────
    point_df$pred_lin <- predict(model_lin, newdata = point_df, type = type_prediction)
    point_df$pred_quad <- predict(model_quad, newdata = point_df, type = type_prediction)
    
    # ─── Prédictions sur un intervalle de la variable x ──────
    # data for prediction
    x_newdata <- seq(min(data_filtered[[trait_temp]], na.rm = TRUE),
                 max(data_filtered[[trait_temp]], na.rm = TRUE), by = 0.05)
    
    
    newdata <- as.data.frame(matrix(0, nrow = length(x_newdata), ncol = length(trait)))
    colnames(newdata) <- trait
    newdata[[trait_temp]] <- x_newdata
    
     # pred
      pred_lin <- predict(model_lin, newdata = newdata, se.fit = TRUE, type = type_prediction)
      pred_quad <- predict(model_quad, newdata = newdata, se.fit = TRUE, type = type_prediction)

    
    # ─── Construction des data.frames pour ggplot ────────────
    df_pred_lin <- data.frame(
      x = newdata[[trait_temp]],
      pred = pred_lin$fit,
      ci_lower = pred_lin$fit - 1.96 * pred_lin$se.fit,
      ci_upper = pred_lin$fit + 1.96 * pred_lin$se.fit,
      model = "Directionnelle"
    )
    
    df_pred_quad <- data.frame(
      x = newdata[[trait_temp]],
      pred = pred_quad$fit,
      ci_lower = pred_quad$fit - 1.96 * pred_quad$se.fit,
      ci_upper = pred_quad$fit + 1.96 * pred_quad$se.fit,
      model = "Quadratique"
    )
    
    df_preds <- rbind(df_pred_lin, df_pred_quad)
    
    # ─── Transformation selon la famille (GLM) ───────────────
    if(family == "binomial"){
      df_preds[, c("pred", "ci_lower", "ci_upper")] <- apply(df_preds[, c("pred", "ci_lower", "ci_upper")], 2, plogis)
      point_df$pred_lin <- plogis(point_df$pred_lin)
      point_df$pred_quad <- plogis(point_df$pred_quad)
    }
    
    if(family == "poisson"){
      df_preds[, c("pred", "ci_lower", "ci_upper")] <- apply(df_preds[, c("pred", "ci_lower", "ci_upper")], 2, exp)
      point_df$pred_lin <- exp(point_df$pred_lin)
      point_df$pred_quad <- exp(point_df$pred_quad)
    }
    
    
    
    
    # ─── Génération du graphique ─────────────────────────────
    # ADD ALPHA FOR SIGNIFICANCE
# FOR PREDICTION
df_preds$significance <- ifelse(significance == 0, 0.15, 1)

# FOR SE
df_preds$significance_se <- ifelse(significance == 0, 0.03, 0.1)

# FOR POINTS AT x == 0
point_df$significance <- ifelse(significance == 0, 0.2, 1)

p <- ggplot() +
  
  # Courbe quadratique + intervalle
  geom_line(
    data = df_preds,
    aes(
      x = x,
      y = pred,
      color = model),
      alpha = df_preds$significance[1],
    size = 1.3
  ) +
  geom_ribbon(
    data = df_preds,
    aes(
      x = x,
      ymin = ci_lower,
      ymax = ci_upper,
      fill = model
    ), alpha = df_preds$significance_se[1]
  ) +
  # Point de prédiction à trait == 0
  geom_point(
    data = data.frame(x = 0, y = point_df$pred_lin),
    aes(x = x, y = y),
    shape = 21,
    size = 4, 
    alpha = point_df$significance[1],
    stroke = 2,
    fill = "white",
    color = "black"
  ) +
  
  theme_classic(base_size = 20) +
  ggtitle(ggtitle) +
  ylab(ylab) +
  xlab(xlab) +
  
  # Légendes
  scale_color_manual(
    name = "Modèle",
    values = c(
      "Directionnelle" = "#27408B",
      "Quadratique" = "#E76F51"
    )
  ) +
  scale_fill_manual(
    name = "Intervalle de confiance",
    values = c(
      "Directionnelle" = "#0000FF",
      "Quadratique" = "#E76F51"
    )
  ) +
  theme(legend.position = "none")
    p
    
    # Assign plot
    assign(x = paste0(trait_temp, "_plot"), value = p)
  }
  
  
  # TRAIT TO REPRESENT
  output <- get(paste0(trait_to_represent, "_plot"))
  
  return(output)
}

3.2 Selection differential

Définition

Il peut être définit comme la covariance entre le trait et la fitness relative : \[ S_k = Cov(\text{trait}, w_k) \] et mesure l’action directe et indirecte de la sélection sur un trait.

Estimation de $S$

Le différentiel de sélection pour chaque trait a été calculé comme le coefficient de régression de la fitness relative à l’épisode de sélection $k$ et du trait seul :

\[ w_{k, i} = \color{blue}{\beta_{k}} \cdot\text{trait}_i + \color{blue}{\gamma _{k}} \cdot\text{trait}^2_i +\epsilon \]

Comme les traits ont été normalisés, le coefficient de régression $\color{blue}{\beta _{k, i}}$ est complètement équivalent à $Cov(\text{trait}_i, w_k)$ car $\sigma_{trait_i} = 1$. J’ai aussi testé une association quadratique entre le trait et la fitness pour détecter un éventuel seuil (e.g. sélection stabilisante).

Indépendance des épisodes de sélection

A noter que la sélection est estimée en traitant chacun des épisodes de sélection $k$ comme indépendendant. Pour cela, la relation entre le trait et la fitness de l’individu $i$ à l’épisode de sélection $k$ a été estimée conditionnellement à l’ensemble des épisodes de sélection précédent $k$. Concrètement, l’idée est de prendre en compte le changement de distribution du trait après l’ensemble des épisodes de sélection précédent. Pour cela, on va utiliser l’argument weights de la fonction lm qui pondère le calcul des moindres carrés pour fit le modèle avec le vecteur précisé dans l’argument :

\[ \sum^n_{i=1}w_i \cdot (y_i - \hat{y}_i)^2 \]

Dans notre cas, weights représente les proportions marginales de la fitness relative (i.e. $w_{i,k-1}$ / $\bar{w}_{k-1}$) de l’épisode de sélection précédent $w_{k - 1}$. En d’autres termes, chaque individu aura un poids proportionnel à sa fitness relative (Arnold & Wade, 1984; Arnold, 2023).

Boostrap procedure

Note importante : ici comme dans le reste du document, les estimations moyennes des paramètres calculés et leurs intervalles de confiance à 95% ont été obtenus via une procédure de bootstrap. Les tirages aléatoires avec remises étaient stratifiés par année et par territoire (i.e. tirage avec remise au sein des cluster année * territoire), pour éviter de sur-représenter les années / territoires avec de plus grands effectifs. La figure ci-dessous décrit la procédure générale de la méthode. A noter également qu’à chaque itération du bootstrap, la fitness relative recalculée (de façon à ce que $\bar{w}_t$ soit toujours 1), condition nécessaire pour que les équations de réponse à la sélection soient valides.

En prenant en compte la sélection directe et indirecte sur les traits, on observe une sélection directionnelle positive sur la masse pour la survie à la première ($w1$) et deuxième année ($w2$), l’accès à la dominance ($w4$) et la fitness totale ($LRS$). Donc plus la masse marmotton est élevée (faibles), plus les chances de survie au premier et deuxième hiver est importante (faible), et plus la probabilité d’accès à la reproduction et le nombre de descendants total est important (faible). On observe également une association non-linéaire entre la masse et la survie à la première année, avec un plateau pour les masses dépassant 2 écart-types au-dessus de la moyenne (voir la partie suivante qui visualise ces relations).

Pour la taille, les résultats montrent une sélection (directe et indirecte) positive pour la survie à la première année ($w1$) et la fitness totale ($LRS$).

3.2.1 Visualisation des différentiels de sélection

J’ai représenté ici différentiel de sélection $S$ estimées précédemment.

Les points blancs correspondent à la prédiction à trait == 0, soit la moyenne de la population. Les graphiques “floutés” correspondent aux effets non-significatifs.

Modèle de sélection directionnelle Modèle de sélection quadratique

3.3 Gradient de sélection

Les gradients de sélection mesure l’effet de la sélection directe sur un trait. On mesure les gradients $\beta$ à l’épisode de sélection $k$ pour les traits $i$ et $j$ avec le modèle suivant :

\[ w_{k,i} = \color{blue}{\beta _{1, k}} \cdot \text{mass}_i + \color{blue}{\gamma _{1, k}} \cdot \text{mass}^2_i + \color{blue}{\beta_{2, k}} \cdot \text{tibia}_i + \color{blue}{\gamma_{2, k}} \cdot \text{tibia}^2_i +\epsilon \]

avec $w_k$ la fitness relative de l’individu $i$ à l’épisode de sélection $k$.

Les résultats de sélection directe sur chaque trait (i.e. en prenant en compte la sélection agissant sur l’autre trait) montrent que c’est surtout la masse qui est sous sélection positive, avec une action directe de la sélection sur $w_1$, $w_4$ et la fitness totale $LRS$. De façon étonnante, on observe une sélection négative de la taille pour l’accès à la dominance $w_4$. Je me demande si cela ne peut pas être un artefact statistique, et si je n’aurais pas dû corriger les p-valeurs pour multi-testing (e.g. avec une procédure de Benjamini et Hochberg, Bonferroni…), car en soit je fais un modèle par proxy de fitness.

3.3.1 Visualisation des gradients de sélection

Les points blancs correspondent à la prédiction à trait == 0, soit la moyenne de la population. Les graphiques “floutés” correspondent aux effets non-significatifs.

Modèle de sélection directionnelle Modèle de sélection quadratique

4 VARIATION SPATIALE DE LA FITNESS ET PHENOLOGIE

[Partie bonus post-terrain…]

On explore si des différences fixes (ou non) de fitness existent entre les différents territoires suivis. Ici est donc représenté la fitness relative moyenne ($\pm$ IC95%) (i.e. pondéré par la fitness de la population à l’annnée $t$) par territoire. Je n’ai considéré ici que les territoires avec plus de $10$ observations. Je n’ai représenté que $w_1$ la survie au premier hiver.

J’ai également représenté les dates d’émergence des différents territoires. Ces dates d’émergences ont été normalisées par année, c’est à dire exprimées en nombre de jour d’émergence relativement à la moyenne de la population à l’année $t$.

J’ai représenté en couleur l’exposition des territoires

Nord Vallée Sud

On remarque des différences fortes de survie entre les territoires avec une forte structuration liée à l’exposition du territoire (Nord, Vallée, Sud). Les territoires exposés Nord et Vallée montrent une survie des marmottons plus élevée que les territoires orientés Sud. Pour les dates d’émergences, les de fortes variabilités existent entre les territoires également (e.g. jusqu’à 10 jours de ‘retard’ pour le territoire $s$ par rapport à la moyenne de la population). Ces différences sont également fortement structurés par l’exposition.

De façon très intéressante, les territoires qui émergent les plus tard sont ceux avec la meilleur survie juvénile. Cela suppose qu’une émergence relative plus précoce n’est pas avantageuse en terme de survie juvénile.

5 VARIATION TEMPORELLE DE LA SELECTION

🧠 Hypothèse
Sélection positive sur les deux traits sur l’ensemble de la période, avec une sélection plus importante pour les années < 2005 comparées à celles > 2005.

Dans cette partie, on va estimer $S$ et $\beta$ en fonction du temps. Je n’ai représenté que les coefficients de sélection directionnelle pour la survie à la première année ($w_1$) et la fitness totale ($w$).

5.1 Fonctions

5.1.1 compute_selection_across_time()

Fonction permettant de computer les différentiels et les gradients de sélection et l’intervalle à 95% par bootstrap, en fonction du temps.

# FUNCTION TO COMPUTE S ACROSS TIME
compute_selection_across_time <- function(data,
                                          trait,
                                          fitness_proxy,
                                          iteration = 1000,
                                          weights,
                                          minimum_sample_size = 30) {

  # -------------------------
  # 1. INITIALISATION
  # -------------------------

  data_boot <- data

  # années comme variables continues
  # data_boot$annee_emergence <- as.numeric(as.character(data_boot$annee_emergence))
  
  # Supprimer les individus avec NA pour les traits sélectionnés
  for(i in seq_along(trait)){
    data_boot <- data_boot[!is.na(data_boot[[trait[i]]]),]
  }

  # Retirer les années avec trop peu d’individus (selon seuil minimum)
  year_sample_size <- data_boot |> 
    group_by(annee_emergence) |> 
    mutate(n = n()) 
  
  rows_to_keep <- year_sample_size$n > minimum_sample_size
  data_boot <- data_boot[rows_to_keep, ]

  # Attribuer un ID unique à chaque ligne pour le bootstrap
  data_boot$row_id <- seq_len(nrow(data_boot))

  # Créer les noms de colonnes pour la matrice de résultats bootstrap
  year <- as.character(sort(unique(data_boot$annee_emergence)))
  colnames_temp <- unlist(lapply(trait, function(t) paste0(t, "_", year)))

  # Initialisation de la matrice pour stocker les coefficients bootstrap
  cov_boot <- as.data.frame(matrix(
    NA,
    ncol = length(colnames_temp),
    nrow = iteration,
    dimnames = list(NULL, colnames_temp)
  ))

  # -------------------------
  # 2. FORMULE DU MODELE
  # -------------------------

  # Crée la formule du modèle : interaction des traits avec l'année
  formula_mod <- as.formula(paste0(
    fitness_proxy,
    "~",
    paste0(trait, collapse = "+"),
    "+",
    paste0(trait, ":annee_emergence", collapse = "+")
  ))


  # ID unique pour chaque ligne
  data_boot$row_id <- seq_len(nrow(data_boot))  # utile pour le bootstrap

  # -------------------------
  # 3. BOUCLE DE BOOTSTRAP
  # -------------------------

  for (i in 1:iteration) {
    # Bootstrap stratifié par année ET territoire
    index_random <- data_boot |>
      group_by(annee_emergence, territoire) |>  # stratification par année / territoire
      reframe(row_id = sample(row_id, size = n(), replace = TRUE)) |>
      pull(row_id)

    # Sous-échantillon tiré avec remise
    data_boot_i <- data_boot[index_random, ]

    # Extraire le numéro de l’épisode de fitness (ex: "1" pour "w1")
    k <- gsub(pattern = "w", replacement = "", x = fitness_proxy)
  
    # Identifier la colonne de fitness absolue (ex: W1)
        # Sélection des individus valides
    if(fitness_proxy == "LRS_relative"){
      fitness_abs <- sym("LRS")
    } else {
      fitness_abs <- sym(paste0("W", k))
    }
    

    # -------------------------
    # 4. CALCUL DE LA FITNESS RELATIVE
    # -------------------------

    if (!is.null(weights)) {
      weight_var <- sym(weights)
      data_boot_i <- data_boot_i |>
        group_by(annee_emergence) |>
        filter(!is.na(!!weight_var) & !is.na(!!fitness_abs)) |>
        mutate(fitness = !!fitness_abs / weighted.mean(!!fitness_abs, w = !!weight_var)) |>
        ungroup()
      
      # Éviter les 0 vrais dans les poids
      weights_mod <- data_boot_i[[weights]]
      weights_mod <- weights_mod / sum(weights_mod, na.rm = TRUE)
      
    } else {
      # Cas sans poids : fitness relative = fitness / moyenne annuelle
      data_boot_i <- data_boot_i |>
        group_by(annee_emergence) |>
        filter(!is.na(!!fitness_abs)) |>
        mutate(fitness = !!fitness_abs / mean(!!fitness_abs)) |>
        ungroup()
      weights_mod <- NULL
    }
    
    # -------------------------
    # 5. MODELISATION
    # -------------------------

    # Modèle linéaire avec interactions trait * année
    mod_temp <- lm(
      formula = formula_mod,
      data = data_boot_i,
      weights = weights_mod
    )

    # -------------------------
    # 6. EXTRACTION DES GRADIENTS DE SELECTION (S et beta)
    # -------------------------

    # Estimation des gradients directionnels par année (i.e., pentes)
    
    params <- c()
    for (j in seq_along(trait)) {
      # emtrends permet d’obtenir la pente de chaque trait pour chaque année
      params <- c(params, as.data.frame(emmeans::emtrends(mod_temp, ~ annee_emergence, var = trait[j]))[, 2])
    }

    # Enregistrement des coefficients dans la matrice bootstrap
    cov_boot[i,] <- params
  }

  # -------------------------
  # 7. STATISTIQUES DE SYNTHESE
  # -------------------------

  # Moyenne des coefficients bootstrap
  mean_boot <- apply(cov_boot, 2, mean, na.rm = TRUE)

  # Intervalles de confiance bootstrap à 95 %
  ci_lower <- apply(cov_boot, 2, quantile, probs = 0.025, na.rm = TRUE)
  ci_upper <- apply(cov_boot, 2, quantile, probs = 0.975, na.rm = TRUE)

  # -------------------------
  # 8. CALCUL DES P-VALEURS BOOTSTRAP
  # -------------------------

  p_values_function <- function(p){
    ifelse(
      p < 0.001 | p > 0.999, "*",
      ifelse(p < 0.01 | p > 0.99, "*",
             ifelse(p < 0.05 | p > 0.95, "*",
                    ifelse(p <= 0.06 | p >= 0.94, ".",
                           ifelse(p < 0.06 | p > 0.94 | is.na(p), yes = "", no = "")))
      )
    )
  }

  # Proportion de coefficients < 0 pour déterminer le signe
  p <- apply(cov_boot, 2, function(x) p_values_function(sum(x < 0, na.rm = TRUE)/iteration))

  # -------------------------
  # 9. FORMATAGE POUR ANALYSE OU VISUALISATION
  # -------------------------

  df_gg <- data.frame(
    y = mean_boot,
    ci_lower = ci_lower,
    ci_upper = ci_upper,
    trait = rep(trait, each = length(year)),     # correspondance trait / année
    year = rep(year, length(trait)),
    p = p
  )

  # Mise au format facteur
  df_gg$trait <- factor(df_gg$trait, levels = trait)
  df_gg$year <- factor(df_gg$year, levels = year)

  # Résultat final
  return(df_gg)
}




# FUNCTION TO BIND SELECTION DIFFERENTIAL OF DIFFERENT TRAIT IN A SAME TABLE 
#' Easier to summarise results
bind_selection_across_time <- function(trait = c("masse_scaled", "tibia_scaled"),
                           data = data,
                           fitness_proxy = fitness_proxy,
                           iteration = iteration,
                           weights = NULL,
                           minimum_sample_size = 5) {
  
  output <- c()
  for (trait_temp in trait) {
    output <- rbind(output, compute_selection_across_time(
      data = data,
      fitness_proxy = fitness_proxy,
      iteration = iteration,
      weights = weights,
      trait = trait_temp
    )
    )
  }
  
  return(output)
}

5.1.2 plot_compute_selection_across_time()

Fonction de plot de la variation temporelle de la sélection.

plot_compute_selection_across_time <- function(data,
                                               title,
                                               ylab,
                                               base_size,
                                               size_point,
                                               size_es, 
                                               size_text = 8) {
  # Facet names
  trait_name <- c(`masse_scaled` = "Mass", `tibia_scaled` = "Tibia")
  
  # Plot
  ggplot(data = data, aes(y = y, x = as.numeric(as.character(year)))) +
    facet_grid(~ trait, labeller = as_labeller(trait_name)) +
    geom_hline(
      yintercept = 0,
      linetype = 2,
      lwd = 1
    ) +
    geom_text(
      mapping = aes(
        x = as.numeric(as.character(year)),
        y = ci_upper + 0.02,
        label = p
      ),
      size = size_text,
      hjust = 0.5,
    ) +
    geom_errorbar(
      aes(ymin = ci_lower, ymax = ci_upper),
      width = 0,
      alpha = 0.3,
      lwd = size_es
    ) +
    geom_point(size = size_point) +
    ggtitle(title) +
    xlab("") +
    ylab(ylab) +
    theme_classic(base_size  = base_size) +
    theme(panel.grid.major.y = element_line(
      color =  Fragman:::transp("#474747", 0.3),
      size = 0.5,
      linetype = 2
    ))
}

5.2 Différentiel de sélection

Comme dans la partie Selection differential, l’idée ici est de calculer $Cov_t(w_k, trait)$ par épisode de sélection $k$, mais ici on calcule $S$ pour l’année $t$. A noter ici que l’année a été modélisée comme une variable catégorielle / facteure. La sélection par année correspond donc au $\beta_{k, t}$ spécifique à l’année $t$.

\[ w_{i,k} = \beta_{0,k} + \beta_{k} \cdot \text{trait}_i + \color{blue}{\beta_{k,t}} \cdot (\text{year}_t : \text{trait}_i) + \epsilon \]

Sample size par année
(En retirant les $NA$ de masse et de tibia)

##    annee_emergence  n
## 1             1996 19
## 2             1997 44
## 3             1998 39
## 4             1999 33
## 5             2000 24
## 6             2001 33
## 7             2002 29
## 8             2003 43
## 9             2004 39
## 10            2005 52
## 11            2006 25
## 12            2007 47
## 13            2008 53
## 14            2009 61
## 15            2010 51
## 16            2011 43
## 17            2012 54
## 18            2013 55
## 19            2014 91
## 20            2015 79
## 21            2016 76
## 22            2017 45
## 23            2018 60
## 24            2019 73
## 25            2020 59

La sélection $S$ est globalement positive pour les deux traits sur l’ensemble de l’année. Le pattern entre les deux périodes n’est en revanche pas si clair, à part pour le LRS et la masse, où la sélection est clairement plus importante en période 1. Les intervalles de confiances sont très importants, il est probable que l’on manque de puissance pour évaluer la covariance trait/fitness à chaque niveau de la variable année…

5.3 Gradient de sélection

On fit le modèle suivant pour déterminer les gradients de sélection par année :

\[ w_{k, i} = \beta_{0,k} + \beta_{1,k} \cdot \text{mass}_i + \beta_{2,k} \cdot \text{tibia}_i + \color{blue}{\beta_{(1,2),k,t} \cdot \text{year}_t :(\text{mass}_i + \text{tibia}_i)} + \epsilon \]

6 FITNESS ET CONTEXTE SOCIALE

J’ai testé ici comment le contexte sociale influence les différents proxy de fitness. Pour cela, j’ai testé un effet linéaire et quadratique pour chaque variable sociale.

\[ \begin{align} w_k &= \beta_{0,k} \\&+ \color{blue}{\beta_{k,s}^{(lineaire)}} \sum_s\text{social}_s \\&+ \color{blue}{\beta_{k,s}^{(non \ lineaire)}} \sum_s\text{social}_s^2 \\&+ \epsilon \end{align} \]

avec $s$ la variable sociale (taille de groupe, taille de portée ou nombre de helpers), $k$ l’épisode de sélection. J’ai ensuite représenté les effets partiels de chaque variable sociale, les autres variables sociales étant fixées à leur moyenne.

Les variables sociales ont été centrées de façon à compute l’effet d’une variable sociale donnée en fonction des valeurs moyennes des autres. Les estimations des paramètres et les p-valeurs ont été obtneues via boostrap.

6.1 Fonctions

6.1.1 partial_effects_social()

partial_effects_social <- function(data,
                                      proxy_fitness,
                                      weights = NULL,
                                      iteration = 100) {
  
  predictors <- c("nb_male_subordinates_scaled", "taille_groupe_scaled", "litter_size_scaled")
  formula_mod <- as.formula(paste0(proxy_fitness, "~ 
    (nb_male_subordinates_scaled + I(nb_male_subordinates_scaled^2)) +
    (taille_groupe_scaled + I(taille_groupe_scaled^2)) +
    (litter_size_scaled + I(litter_size_scaled^2))"
  ))
  
  data$row_id <- seq_len(nrow(data))
  all_preds <- list()
  
  get_significance <- function(p) {
    if (p < 0.001) return("***")
    else if (p < 0.01) return("**")
    else if (p < 0.05) return("*")
    else if (p < 0.1) return(".")
    else return("")
  }
  
  for (pred in predictors) {
    x_vals <- seq(min(data[[pred]], na.rm = TRUE), max(data[[pred]], na.rm = TRUE), length.out = 1000)
    
    newdata <- data.frame(matrix(ncol = length(predictors), nrow = 1000))
    colnames(newdata) <- predictors
    for (v in predictors) {
      newdata[[v]] <- if (v == pred) x_vals else mean(data_mod[[v]],  # MEAN OF THE POPULATION
                                                      na.rm = TRUE)
    }
    for (v in predictors) {
      newdata[[paste0("I(", v, "^2)")]] <- newdata[[v]]^2
    }
    
    pred_matrix <- matrix(NA, nrow = iteration, ncol = 1000)
    coefs_linear <- numeric(iteration)
    coefs_quad <- numeric(iteration)
    
    for (i in 1:iteration) {
      index_random <- data %>%
        group_by(annee_emergence, territoire) %>%
        reframe(row_id = sample(row_id, size = n(), replace = TRUE)) %>%
        pull(row_id)
      
      data_boot_i <- data[index_random, ]
      
      if (proxy_fitness == "LRS_relative") {
        fitness_abs <- sym("LRS")
      } else {
        k <- gsub("w", "", proxy_fitness)
        fitness_abs <- sym(paste0("W", k))
      }
      
      if (!is.null(weights)) {
        weight_var <- sym(weights)
        data_boot_i <- data_boot_i %>%
          group_by(annee_emergence) %>%
          filter(!is.na(!!weight_var) & !is.na(!!fitness_abs)) %>%
          mutate(fitness = !!fitness_abs / weighted.mean(!!fitness_abs, w = !!weight_var)) %>%
          ungroup()
        
        weights_mod <- data_boot_i[[weights]]
        weights_mod <- weights_mod / sum(weights_mod, na.rm = TRUE) + 1e-8
      } else {
        data_boot_i <- data_boot_i %>%
          group_by(annee_emergence) %>%
          filter(!is.na(!!fitness_abs)) %>%
          mutate(fitness = !!fitness_abs / mean(!!fitness_abs)) %>%
          ungroup()
        weights_mod <- NULL
      }
      
      mod_boot <- lm(formula_mod, data = data_boot_i, weights = weights_mod)
      pred_matrix[i, ] <- predict(mod_boot, newdata = newdata)
      
      coefs <- coef(mod_boot)
      coefs_linear[i] <- coefs[pred]
      coefs_quad[i] <- coefs[paste0("I(", pred, "^2)")]
    }
    
    # Significativité globale du prédicteur (linéaire + quadratique)
    p_linear <- 2 * min(mean(coefs_linear > 0), mean(coefs_linear < 0))
    p_quad <- 2 * min(mean(coefs_quad > 0), mean(coefs_quad < 0))
    
    # Choisir le plus petit des deux p-values comme indicateur de l’effet global
    p_min <- min(p_linear, p_quad)
    signif_label <- get_significance(p_min)
    
    pred_df <- data.frame(
      x = x_vals,
      fit = apply(pred_matrix, 2, mean, na.rm = TRUE),
      ci_lower = apply(pred_matrix, 2, quantile, probs = 0.025, na.rm = TRUE),
      ci_upper = apply(pred_matrix, 2, quantile, probs = 0.975, na.rm = TRUE),
      predictor = pred,
      significance = signif_label
    )
    
    all_preds[[pred]] <- pred_df
  }
  
  plot_data <- bind_rows(all_preds)
  
  return(plot_data)
}

6.1.2 plot_partial_effects_social()

plot_partial_effects_social <- function(data, base_size = 23, ylab, period = F){
  
  data$label <- paste0(
    dplyr::recode(data$predictor,
                  "nb_male_subordinates_scaled" = "Helpers",
                  "taille_groupe_scaled" = "Group size",
                  "litter_size_scaled" = "Litter size")
  )
  
  # TRANSPRENCY DEPENDING ON SIGNIFICANCE
  data <- data |> 
    mutate(
      significance_num = case_when(
        significance %in% c("", ".") ~ 0,
        significance %in% c("*", "**", "***") ~ 1
      )
    ) 
  
  # Rename factor
  data$label <- factor(data$label, levels = c("Group size", "Litter size", "Helpers"), labels = c("Group size", "Litter size", "Helpers"))
  
  data <- data |> 
    arrange(label)
  
  
  # FOR PREDICTION
  data$significance_pred <- ifelse(data$significance_num == 0, 0.2, 1)
  
  # FOR SE
  data$significance_se <- ifelse(data$significance_num == 0, 0.1, 0.2)
  
  # PLOT
  plot <- ggplot(data, aes(x = x, y = fit)) +
    geom_hline(yintercept = 1, linetype = 2, alpha = data$significance_pred[!duplicated(data$label)]
    ) + 
    geom_vline(xintercept = 0, linetype = 2, alpha = data$significance_pred[!duplicated(data$label)]
    ) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper),
                fill = "mediumpurple", alpha = data$significance_se) +
    geom_line(color = "mediumpurple", size = 1.5,
              alpha = data$significance_pred) +
    facet_wrap(~ label, scales = "free_x") +
    labs(
      x = "",
      y = ylab
    ) +
    theme_classic(base_size = base_size) +
    annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf) +
    annotate("segment", x = -Inf, xend = -Inf, y = -Inf, yend = Inf) +
    theme(legend.position = "none", legend.title = element_blank())
  
  return(plot)
}

6.2 Fitness across social

Le contexte social a un effet principale sur la survie des individus. En effet, les résultats montrent un effet de la taille de portée et de la taille de groupe. Une grande taille de portée (supérieur à la moyenne [0]) diminue les chances de survie à la première hibernation $w1$. La survie à la deuxième et troisième année n’est pas influencé directement par la taille de portée. La taille de groupe influence la survie à la troisième année, avec un effet linéaire positif sur la probabilité de survie de 2 à 3 ans.

Concernant l’accès à la reproduction et le succès reproducteur, là également on observe un effet du social. La probabilité d’accès à la dominance est modulée par la taille de groupe, qui a un effet négatif sur celle-ci i.e. plus la taille de groupe est élevée, moins la probabilité d’accès à la dominance est importante.

Enfin, la fitness totale $w$ (LRS) est influencé par la taille de groupe et la taille de portée, avec une fitness plus importante dans les grands groupes à la naissance et dans les petites tailles de portée.

A voir si ces effets changent entre la période 1 et la période 2 ?

7 SELECTION ET CONTEXTE SOCIAL

Les estimations des paramètres et les p-valeurs ont été obtenues via boostrap

Dans cette partie, on va tester si la sélection sur la taille et la masse est modulée par le contexte sociale de l’individu. Autrement dit, on teste si la relation trait/fitness varie en fonction de la valeur des variables sociales (e.g. valeur de la sélection à taille de groupe == 1 et taille de groupe == 2).

Pour cela, on va tester l’association trait/fitness en interaction avec les variables sociales. J’ai testé l’effet linéaire ainsi qu’un effet social^2 pour visualiser la présence d’un éventuel seuil de l’effet du social sur la sélection pour la taille de portée, la taille de groupe et le nombre de helpers car il y avait assez de modalités pour fit la variable sociale comme quantitative. Les 3 variables sociales ont été ajoutées simultanément dans le modèle pour rendre compte de leur corrélation.

Pour chaque proxy de fitness, les modèles étaient les suivants :

Différentiel de sélection \[ \begin{align} w_k &= \beta_{0,k} + \beta _{i,k} \cdot \text{trait}_{i} + \sum_s\text{social}_s + \sum_s\text{social}_s^2 \\&+ \sum_s\color{blue}{\gamma_{i,k,s}^{(lineaire)}} \cdot (\text{trait}_i:\text{social}_s) \\&+ \sum_s\color{blue}{\gamma_{i,k,s}^{(non \ lineaire)}} \cdot (\text{trait}_i:\text{social}_s^2) \\&+ \epsilon \end{align} \] avec $k$ l’épisode de sélection $k$, $s$ la variable sociale (taille de groupe, taille de portée et nombre de helpers) et $i$ le trait (masse ou tibia).

Gradient de sélection
Même logique pour le gradient de sélection : \[ \begin{align} w_k &= \beta_{0,k} + \sum_i \beta_{i,k} \cdot \text{trait}_i + \text{social}_s + \text{social}_s^2 \\&+ \sum_s\sum_i \color{blue}{\gamma_{i,k,s}^{(lineaire)}} \cdot (\text{trait}_i : \text{social}_s) \\&+ \sum_s\sum_i \color{blue}{\gamma_{i,k,s}^{(non \ lineaire)}} \cdot (\text{trait}_i : \text{social}_s^2) \\&+ \epsilon \end{align} \]

La significativité de l’effet du social sur la sélection était ensuite testé pour l’effet linéaire de la variable sociale sur le trait ($\gamma_{i,k}^{(lineaire)}$) et pour l’effet non-linéaire ($\gamma_{i,k}^{(non \ lineaire)}$). Pour calculer $\hat{\beta}_{i,k,s}({\text{social} = x})$ la valeur de sélection à une valeur donnée de la variable sociale, j’ai reconstruit le gradient à la valeur sociale $x$ avec la formule :

\[ \hat{\beta}_{i,k,s}({\text{social} = x}) = \beta_{i,k} + x \gamma^{(linéaire)}_{i,k,s} + x^2 \gamma^{(non \ linéaire)}_{i,k,s} \]

Les variables sociales ont été centrées pour (1) éviter une collinéarité trop forte en l’effet linéaire et quadratique et (2) pour pouvoir comparer les $\hat{\beta}$ en fixant les autres effets à 0, 0 correspondant à la moyenne de la variable sociale dans la population (et non pas à e.g. une taille de groupe de 0)

7.1 Fonctions

7.1.1 compute_boot_model_results_social()

Fonction permettant de générer les $n$ itérations de coefficients estimés par le modèle décrit ci-dessus.

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

#                          SELECTION PAR STRUCTURE SOCIALE
#                                 (Significativité)

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

compute_boot_model_results_social <- function(data = data_fitness,
                                              trait,
                                              proxy_fitness,
                                              weights = NULL,
                                              iteration = 50){
 
   # MODEL FORMULA
  formula_mod <- as.formula(paste0(proxy_fitness, "~ 
  (presence_male_subordinates_scaled +
     I(presence_male_subordinates_scaled^2))*(", paste0(trait, collapse = "+"), ") +
  
  (taille_groupe_scaled +
     I(taille_groupe_scaled^2))*(", paste0(trait, collapse = "+"), ") +
  
  (litter_size_scaled +
       I(litter_size_scaled^2))*(", paste0(trait, collapse = "+"), ")"
  )
  )
  
  
  
  # Attribuer un ID unique à chaque ligne (utile pour le bootstrap)
  data$row_id <- seq_len(nrow(data))
  
  # matrix
  mod_temp <- lm(
    formula = formula_mod,
    data = data,
    weights = NULL
  )
  
  
  coef_mod <- coef(mod_temp)
  cov_boot <- matrix(NA, ncol = length(coef_mod), nrow = iteration, dimnames = list(c(1:iteration), names(coef_mod)))
  
  for (i in 1:iteration) {
    
    # 4.1. Échantillonnage bootstrap stratifié par social_variable x territoire
    index_random <- data |>
      group_by(annee_emergence, territoire) |>
      reframe(row_id = sample(row_id, size = n(), replace = TRUE)) |>
      pull(row_id)
    
    data_boot_i <- data[index_random, ]
    
    # 4.2. Définir la variable de fitness absolue
    if(proxy_fitness == "LRS_relative"){
      fitness_abs <- sym("LRS")
    } else {
      k <- gsub("w", "", proxy_fitness)
      fitness_abs <- sym(paste0("W", k))
    }
    
    # 4.3. Calculer la fitness relative (pondérée ou non)
    if (!is.null(weights)) {
      weight_var <- sym(weights)
      data_boot_i <- data_boot_i |>
        group_by(annee_emergence) |>
        filter(!is.na(!!weight_var) & !is.na(!!fitness_abs)) |>
        mutate(fitness = !!fitness_abs / weighted.mean(!!fitness_abs, w = !!weight_var)) |>
        ungroup()
      
      # Normaliser les poids
      weights_mod <- data_boot_i[[weights]]
      weights_mod <- weights_mod / sum(weights_mod, na.rm = TRUE)
      weights_mod <- weights_mod + 1e-8  # éviter les zéros
      
    } else {
      data_boot_i <- data_boot_i |>
        group_by(annee_emergence) |>
        filter(!is.na(!!fitness_abs)) |>
        mutate(fitness = !!fitness_abs / mean(!!fitness_abs)) |>
        ungroup()
      
      weights_mod <- NULL
    }
    
    # 4.4. Ajuster le modèle linéaire
    mod_temp <- lm(
      formula = formula_mod,
      data = data_boot_i,
      weights = weights_mod
    )
    
    # 4.5. Extraire les gradients de sélection (par modalité sociale)
    cov_boot[i,] <- coef(mod_temp)
  }
  
  # GET PARAMETERS
  cov_boot <- as.data.frame(cov_boot)
  cov_boot_df <- pivot_longer(data = cov_boot,
                              cols = 1:ncol(cov_boot), names_to = "params", values_to = "value")
  cov_boot_df$params <- as.factor(cov_boot_df$params)
  
  
  
  #########
  # RESULTS
  #########
  # Significativité ---------------------------------------------------------
  p_values_function <- function(p) {
    ifelse(p < 0.001 | p > 0.999, "***",
           ifelse(p < 0.01 | p > 0.99, "**",
                  ifelse(p < 0.05 | p > 0.95, "*",
                         ifelse(p < 0.1 | p > 0.90, ".", "")))
    )
  }
  
  res_social <- cov_boot_df |> 
    group_by(params) |> 
    summarise(y = mean(value, na.rm = T),
              ci_lower = quantile(value, probs = 0.025, na.rm = T),
              ci_upper = quantile(value, probs = 0.975, na.rm = T),
              p = min(c(sum(value < 0), sum(value > 0))) / iteration,
              significance = p_values_function(p),
              proxy_fitness = proxy_fitness
    ) |> 
    as.data.frame()
  
  res_social
  
  # SORT ONLY INTERACTIVE EFFECT
  res_social <- res_social[grepl(":", res_social$params), ]
  params_name <- res_social$params
  
  # TYPE OF INTERACTION
  res_social$interaction <- NA
  res_social$interaction[grepl(pattern = "I", x = params_name)] <- "Quadratic"
  res_social$interaction[is.na(res_social$interaction)] <- "Linear"
  
  # SOCIAL VARIABLE
  res_social$social <- NA
  
  res_social$social[grepl(pattern = "subordinates", x = params_name)] <- "Nb Helpers"
  res_social$social[grepl(pattern = "group", x = params_name)] <- "Group size"
  res_social$social[grepl(pattern = "litter", x = params_name)] <- "Litter size"
  
  # TRAIT
  res_social$trait <- NA
  for(i in trait){
    res_social$trait[grepl(pattern = i, x = params_name)] <- i
  }
  
  res_social <- res_social |> 
    arrange(trait, interaction, social)
  
  return(list(res_social = res_social, cov_boot = cov_boot, params_name = params_name))
}

7.1.2 plot_social_model_result()

Fonction représentant visuellement les valeurs et la significativité des paramètres d’intéraction linéaire et quadratique entre le trait et la variable sociale.

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

#                          SELECTION PAR STRUCTURE SOCIALE
#                             (Significativité - PLOT)

# -------------------------------------------------------------------------
plot_social_model_result <- function(res_social,
                                     base_size = 20,
                                     proxy_tested,
                                     ylab = "Selection differential",
                                     dodge = 0.3,
                                     size_point = 4,
                                     size_es = 1,
                                     alpha_es = 0.5,
                                     size_text = 8){
 # FORMATTING
  # sort by proxy tested
  res_social <- res_social[res_social$proxy_fitness %in% proxy_tested, ]
  
  # relevel
  res_social$proxy_fitness <- factor(res_social$proxy_fitness, 
                                     levels = unique(res_social$proxy_fitness), labels = proxy_tested)
  
  # relevel
  res_social$trait <- factor(res_social$trait, 
                                     levels = unique(res_social$trait),
                                     labels = c("Mass", "Tibia")
                             )
  
  # relevel significance to format alpha levels
  res_social$significance[res_social$significance == "."] <- ""
  res_social$significance <- factor(res_social$significance, 
                                    levels = c("", "*", "**", "***"),
                                    labels = c("", "*", "*", "*")
                                    )
  
  # res_social$significance <- factor(res_social$significance, 
  #                                   levels = c("", ".", "*", "**", "***"),
  #                                   labels = c("", ".", "*", "*", "*")
  #                                   )
  
  # PLOT
  plot <-  ggplot(data = res_social, aes(x = social,
                                         y = y,
                                         color = interaction, 
                                         alpha = significance)) +
    facet_grid(proxy_fitness ~ trait, scale = "free") +
    geom_hline(yintercept = 0,
               linetype = 2,
               alpha = 0.4) +
     coord_flip() +
     geom_errorbar(aes(ymin = ci_lower,
                      ymax = ci_upper,
                      color = interaction,
                      alpha = significance),
                  width = 0,
                  size = size_es,
                  position = position_dodge(width = dodge)) +
    geom_point(data = res_social, aes(x = social,
                                         y = y,
                                         color = interaction, 
                                         alpha = significance),
               size = size_point,
               position = position_dodge(width = dodge)) +
     # scale_alpha_manual(breaks = c(0.3, 0.5, 1),
     #                    values = c(0.3, 0.5, 1)) +
     # geom_text(aes(label = significance, y = ci_upper+0.1),
     #           position = position_dodge(width = dodge),
     #           size = size_text) +
    theme_classic(base_size = base_size) +
    ylab(ylab) +
  annotate("segment", x=-Inf, xend=Inf, y=-Inf, yend=-Inf) +
  annotate("segment", x=-Inf, xend=-Inf, y=-Inf, yend=Inf) +
    xlab("Interaction terms") +
  scale_color_manual(values =  c("#27408B",  "#E76F51")) +
    scale_fill_manual(values = c("#0000FF",  "#E76F51")) +
    theme(legend.position = "none", legend.title = element_blank())

      return(plot)
}

7.1.3 summary_table()

Fonction permettant de résumer les effets d’interaction significatifs en un tableau

# TABLE
summary_table <- function(df) {
  
  # FORMATTING
  df$trait <- factor(df$trait, levels = c("masse_scaled", "tibia_scaled"), labels = c("Mass", "Tibia"))
  df$proxy_fitness <- factor(df$proxy_fitness, 
                             levels = c("w1", "w2", "w3", "w4", "w5", "LRS_relative"), 
                             labels = c("w1", "w2", "w3", "w4", "w5", "LRS"))

  # DATAFRAME
  df <- df %>%
    mutate(p = as.numeric(p)) %>%
    filter(!is.na(p)) %>%
    filter(p < 0.05) %>%
    select(proxy_fitness, trait, social, interaction) %>%
    rename(
      Fitness = proxy_fitness,
      Trait = trait,
      Social = social,
      Interaction = interaction
    ) %>%
    mutate(Value = "x") %>%
    pivot_wider(names_from = Interaction, values_from = Value, values_fill = "") %>%
    select(Social, Fitness, Trait, Linear, Quadratic) %>%
    arrange(Social, Fitness, Trait)
  
  # Find the last row index of each group (based on column 1)
group_ends <- which(diff(as.integer(factor(df[[1]]))) != 0)
group_ends <- c(group_ends, nrow(df))  # Include the last row

# Start the kable rendering
table_df <- df %>%
   kableExtra::kable(format = "html", align = "c") %>%
   kableExtra::kable_classic(full_width = FALSE, position = "center", font_size = 25, html_font = "Cambria") %>%
   kableExtra::column_spec(1, bold = TRUE) %>%
   kableExtra::column_spec(c(4, 5), bold = TRUE) %>%
   kableExtra::collapse_rows(columns = 1, valign = "middle")

# Add horizontal lines after each collapsed group
for (i in group_ends) {
  table_df <- table_df %>%
     kableExtra::row_spec(i, extra_css = "border-bottom: 2px solid black;")
}
# 
#   # Affichage en kable avec centrage
# table_df <- df %>%
#   kableExtra::kable(format = "html", align = "c") %>%
#   kableExtra::kable_classic(full_width = FALSE, position = "center", font_size = 25, html_font = "Cambria") %>%
#   kableExtra::column_spec(1, bold = TRUE) %>%  # colonne "Social"
#   kableExtra::column_spec(c(4,5), bold = TRUE) %>%  
#   kableExtra::collapse_rows(columns = 1, valign = "middle")

return(table_df)
}

7.1.4 compute_selection_social_level()

Fonction permettant de calculer la pente fitness / trait (sélection) pour chaque valeur de la variable sociale.

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

#                          SELECTION PAR STRUCTURE SOCIALE
#                             (Effet par niveau social)

# -------------------------------------------------------------------------
compute_selection_social_level <- function(cov_boot = cov_boot,
                                           data, 
                                           social_variable,
                                           trait_tested,
                                           proxy_fitness,
                                           significance,
                                           params_name = params_name){
  
  # FORMAT NAME
  params_name <- colnames(cov_boot)
  
    # DEFINE SOCIAL LEVELS ON WHICH COMPUTE SELECTION
    at = sort(na.omit(unique(data[[paste0(social_variable, "_scaled")]])))
    at_original_scale <- sort(na.omit(unique(data[[social_variable]])))
  
  mat_social_level_selection <- matrix(NA,
                                       ncol = length(at),
                                       nrow = iteration,
                                       dimnames = list(c(1:iteration), at_original_scale))
  
  
  # EXTRAIRE LES COEFFICIENTS APPROPRIES
  names_cov_boot <- colnames(cov_boot)
  coef_equation <- names_cov_boot[names_cov_boot == trait_tested |
  (grepl(pattern = trait_tested, names_cov_boot) & grepl(pattern = social_variable, names_cov_boot))]
  
  df_coef <- cov_boot[, coef_equation]
  
  
  # POUR CHAQUE ITERATION BOOTSTRAP, ON CALCULE LES VALEURS DE SELECTION AU LEVEL X
  compute_selection_social <- function(main_effect,
                                       linear_interaction,
                                       non_linear_interaction, 
                                       at = at){
    # y^_{social == x} = \beta + \Gamma^{Linear} x + \Gamma^{Qudra} x^2
    pred_selection <- main_effect + at*linear_interaction + at^2 * non_linear_interaction
    return(pred_selection)
  }
  
  # on applique la fonction sur chaque itération
  res_apply <- apply(df_coef, 1, function(x) compute_selection_social(x[1], x[2], x[3], at = at))
  
  for(i in 1:iteration){
    mat_social_level_selection[i, ] <- res_apply[, i]
  }
  
  df_mat <- pivot_longer(as.data.frame(mat_social_level_selection),
                         cols = 1:ncol(mat_social_level_selection),
                         names_to = "social_variable_levels",
                         values_to = "y")
  
  # compute mean / CI
  df_mat <- df_mat |> 
    group_by(social_variable_levels) |> 
    summarise(mean = mean(y, na.rm = T),
              ci_lower = quantile(y, probs = 0.025, na.rm = T),
              ci_upper = quantile(y, probs = 0.975, na.rm = T)) 
  
  df_mat$social_variable_levels <- factor(df_mat$social_variable_levels,
                                   levels = sort(as.numeric(as.character(df_mat$social_variable_levels))),
                                   labels = sort(as.numeric(as.character(df_mat$social_variable_levels)))
  )
  
  # formatting
  df_mat$social_variable <- social_variable
  df_mat$fitness_proxy <- proxy_fitness
  df_mat$trait <- trait_tested
  df_mat$significance <- significance
  
  return(df_mat)
}

7.1.6 plot_3d()

7.1.5 plot_compute_selection_social_level()

plot_compute_selection_social_level <- function(data, 
                                                social_tested,
                                                proxy_fitness_tested,
                                                ylab = "Selection differential",
                                                xlab = "",
                                                base_size,
                                                size_point,
                                                size_line,
                                                es_size,
                                                alpha_es){
  
  # FILTER DATA
data_filtered <- 
  data[data$social_variable == social_tested &
         data$fitness_proxy == proxy_fitness_tested, ]

# FORMATTING
data_filtered$social_variable <- 
  factor(data_filtered$social_variable,
         levels = c("taille_groupe", "litter_size", "nb_male_subordinates"),
         labels = c("Group Size", "Litter size", "Nb Helpers"))

data_filtered$trait <- factor(data_filtered$trait, levels = c("masse_scaled", "tibia_scaled"), labels = c("Mass", "Tibia"))

# TRANSPRENCY DEPENDING ON SIGNIFICANCE
data_filtered <- data_filtered |> 
  mutate(
    significance_num = case_when(
      significance %in% c("", ".") ~ 0,
      significance %in% c("*", "**", "***") ~ 1
    )
  ) |> 
  mutate(significance_num = factor(significance_num, levels = c(0, 1)))

# PLOT
plot <- ggplot(
    data_filtered,
    aes(x = as.numeric(as.character(social_variable_levels)),
      y = mean,
      alpha = significance_num)
  ) +
  geom_line(aes(x = as.numeric(as.character(social_variable_levels)),
      y = mean, 
      alpha = significance_num),
      col = Fragman::transp("lightsteelblue1",
                             alpha = 0.1),
      linewidth = size_line) +
  scale_alpha_manual(values = c("0" = 0.2, "1" = 1)) +
   # geom_ribbon(mapping = aes(ymin = ci_lower, ymax = ci_upper),
   #              alpha = 0.1, fill = "blue") +
  ylab(ylab) +
  xlab(xlab) +
    facet_grid(fitness_proxy ~ trait, scale = "free") +
      geom_errorbar(aes(ymin = ci_lower,
                        ymax = ci_upper, 
                        alpha = significance_num),
                    width = 0,
                    size = es_size,
                    col = Fragman::transp("darkgrey",
                             alpha = 0.3)) +
    geom_point(size = size_point) +
  geom_hline(yintercept = 0,
             linetype = 2,
             Fragman::transp("black",
                             alpha = 0.4)) +
    theme_classic(base_size = base_size) +
  theme(legend.position = "none")

return(plot)
}

7.2 Selection differential across social

7.2.1 Significance of social effect on selection

Je teste ici la significativité des intéractions linéaire et quadratique entre les variables sociales et la sélection.

Les éléments “floutés” correspondent aux effets non-significatifs. A noter que je n’ai pas intégrer expliciter les effets avec $0.05 < p < 0.1$ dans le graphique, mais ils sont présents dans l’output du modèle avec le symbole “.”. $Aussi, seules les modalités sociales avec $n \geq 20$ ont été conservées.

Effet linéaire Effet quadratique

Résumé des interactions significatives

Social	Fitness	Trait	Linear	Quadratic
Group size	w1	Mass		x
Group size	LRS	Mass	x
Litter size	w1	Tibia	x
	w2	Mass		x
	w2	Tibia		x
	w3	Mass		x
	w5	Tibia	x
	LRS	Mass		x
	LRS	Tibia	x
Nb Helpers	w3	Mass	x	x
	w3	Tibia		x
	w5	Mass	x	x
	LRS	Mass	x

7.2.2 Selection across social level

On compute ici les différentes valeurs du différentiel de sélection $S$ pour chaque valeur de la variable sociale (les autres variables sociales étant fixées à leur moyenne).

Les éléments “floutés” correspondent aux effets non-significatifs.

7.2.2.1 Group size

Lin : interaction linéaire
Qua : interaction quadratique

##    social trait            w interaction     y   ci-   ci+    p sign
## 1   Group  Mass           w1         Lin -0.01 -0.03  0.02 0.34     
## 4   Group  Mass           w1         Qua  0.02  0.01  0.03 0.00   **
## 7   Group Tibia           w1         Lin -0.02 -0.05  0.01 0.12     
## 10  Group Tibia           w1         Qua  0.00 -0.02  0.01 0.28     
## 13  Group  Mass           w2         Lin -0.01 -0.04  0.02 0.25     
## 16  Group  Mass           w2         Qua  0.00 -0.01  0.01 0.42     
## 19  Group Tibia           w2         Lin  0.00 -0.02  0.03 0.41     
## 22  Group Tibia           w2         Qua  0.00 -0.01  0.01 0.42     
## 25  Group  Mass           w3         Lin -0.03 -0.07  0.01 0.05    .
## 28  Group  Mass           w3         Qua -0.01 -0.03  0.01 0.12     
## 31  Group Tibia           w3         Lin  0.00 -0.05  0.05 0.49     
## 34  Group Tibia           w3         Qua -0.01 -0.03  0.01 0.08    .
## 37  Group  Mass           w4         Lin -0.02 -0.07  0.02 0.18     
## 40  Group  Mass           w4         Qua  0.00 -0.01  0.02 0.39     
## 43  Group Tibia           w4         Lin  0.01 -0.04  0.07 0.35     
## 46  Group Tibia           w4         Qua  0.00 -0.02  0.01 0.35     
## 49  Group  Mass           w5         Lin  0.01 -0.04  0.06 0.39     
## 52  Group  Mass           w5         Qua  0.00 -0.02  0.02 0.43     
## 55  Group Tibia           w5         Lin  0.04 -0.03  0.10 0.12     
## 58  Group Tibia           w5         Qua  0.01 -0.02  0.03 0.30     
## 61  Group  Mass LRS_relative         Lin -0.10 -0.21 -0.01 0.02    *
## 64  Group  Mass LRS_relative         Qua -0.01 -0.07  0.05 0.31     
## 67  Group Tibia LRS_relative         Lin -0.05 -0.18  0.08 0.25     
## 70  Group Tibia LRS_relative         Qua -0.03 -0.11  0.04 0.20

Plot

7.2.2.2 Litter size

Lin : interaction linéaire
Qua : interaction quadratique

##    social trait            w interaction     y   ci-   ci+    p sign
## 2  Litter  Mass           w1         Lin -0.01 -0.07  0.05 0.32     
## 5  Litter  Mass           w1         Qua  0.02 -0.01  0.05 0.07    .
## 8  Litter Tibia           w1         Lin -0.04 -0.09  0.01 0.04    *
## 11 Litter Tibia           w1         Qua -0.02 -0.04  0.01 0.07    .
## 14 Litter  Mass           w2         Lin -0.03 -0.09  0.02 0.12     
## 17 Litter  Mass           w2         Qua  0.02  0.00  0.05 0.02    *
## 20 Litter Tibia           w2         Lin -0.05 -0.10  0.01 0.05    .
## 23 Litter Tibia           w2         Qua  0.02  0.00  0.05 0.02    *
## 26 Litter  Mass           w3         Lin -0.02 -0.14  0.10 0.42     
## 29 Litter  Mass           w3         Qua  0.05  0.00  0.09 0.04    *
## 32 Litter Tibia           w3         Lin  0.04 -0.07  0.15 0.23     
## 35 Litter Tibia           w3         Qua  0.03 -0.01  0.08 0.08    .
## 38 Litter  Mass           w4         Lin  0.08 -0.03  0.18 0.06    .
## 41 Litter  Mass           w4         Qua  0.03 -0.03  0.08 0.16     
## 44 Litter Tibia           w4         Lin  0.02 -0.08  0.12 0.35     
## 47 Litter Tibia           w4         Qua  0.02 -0.04  0.07 0.26     
## 50 Litter  Mass           w5         Lin -0.01 -0.15  0.17 0.41     
## 53 Litter  Mass           w5         Qua  0.02 -0.05  0.11 0.25     
## 56 Litter Tibia           w5         Lin -0.17 -0.34 -0.02 0.01    *
## 59 Litter Tibia           w5         Qua  0.04 -0.04  0.13 0.16     
## 62 Litter  Mass LRS_relative         Lin -0.05 -0.30  0.20 0.36     
## 65 Litter  Mass LRS_relative         Qua  0.22  0.08  0.38 0.00   **
## 68 Litter Tibia LRS_relative         Lin -0.29 -0.56 -0.03 0.02    *
## 71 Litter Tibia LRS_relative         Qua  0.11 -0.02  0.25 0.05    .

Plot

7.2.2.3 Number of helpers

Lin : interaction linéaire
Qua : interaction quadratique

##     social trait            w interaction     y   ci-   ci+    p sign
## 3  Helpers  Mass           w1         Lin -0.02 -0.13  0.08 0.33     
## 6  Helpers  Mass           w1         Qua  0.01 -0.06  0.09 0.36     
## 9  Helpers Tibia           w1         Lin  0.02 -0.09  0.12 0.37     
## 12 Helpers Tibia           w1         Qua -0.01 -0.10  0.07 0.42     
## 15 Helpers  Mass           w2         Lin -0.05 -0.15  0.04 0.14     
## 18 Helpers  Mass           w2         Qua  0.01 -0.07  0.08 0.42     
## 21 Helpers Tibia           w2         Lin  0.01 -0.08  0.09 0.43     
## 24 Helpers Tibia           w2         Qua -0.04 -0.11  0.04 0.15     
## 27 Helpers  Mass           w3         Lin  0.16 -0.02  0.35 0.04    *
## 30 Helpers  Mass           w3         Qua -0.13 -0.26 -0.02 0.01    *
## 33 Helpers Tibia           w3         Lin  0.09 -0.08  0.26 0.14     
## 36 Helpers Tibia           w3         Qua -0.13 -0.26  0.01 0.03    *
## 39 Helpers  Mass           w4         Lin  0.04 -0.13  0.23 0.33     
## 42 Helpers  Mass           w4         Qua -0.02 -0.13  0.07 0.35     
## 45 Helpers Tibia           w4         Lin -0.03 -0.21  0.16 0.37     
## 48 Helpers Tibia           w4         Qua  0.02 -0.11  0.13 0.31     
## 51 Helpers  Mass           w5         Lin  0.28  0.07  0.51 0.00   **
## 54 Helpers  Mass           w5         Qua -0.12 -0.24  0.01 0.03    *
## 57 Helpers Tibia           w5         Lin  0.14 -0.07  0.39 0.10     
## 60 Helpers Tibia           w5         Qua -0.09 -0.22  0.08 0.10    .
## 63 Helpers  Mass LRS_relative         Lin  0.41  0.02  0.80 0.02    *
## 66 Helpers  Mass LRS_relative         Qua -0.24 -0.64  0.13 0.11     
## 69 Helpers Tibia LRS_relative         Lin  0.26 -0.17  0.69 0.12     
## 72 Helpers Tibia LRS_relative         Qua -0.25 -0.87  0.34 0.21

Plot

7.3 Selection gradient across social

On étudie maintenant les gradients de sélection $\beta$ avec la même méthode.

7.3.1 Significance of social effect on selection

Je teste ici la significativité des intéractions linéaire et quadratique entre les variables sociales et la sélection.

Les éléments “floutés” correspondent aux effets non-significatifs. Aussi, seules les modalités sociales avec $n \geq 20$ ont été conservées.

Effet linéaire Effet quadratique

Résumé des interactions significatives

Social	Fitness	Trait	Linear	Quadratic
Group size	w1	Mass		x
	w1	Tibia		x
	w2	Mass	x
	w3	Mass	x
	w4	Mass		x
	LRS	Mass		x
Litter size	w1	Mass		x
	w1	Tibia		x
	w5	Tibia	x
	LRS	Tibia	x
Nb Helpers	w2	Mass	x	x
	w2	Tibia	x	x
	w5	Mass	x
	LRS	Mass	x

7.3.2 Selection across social level

On compute ici les différentes valeurs du différentiel de sélection $\beta$ pour chaque valeur de la variable sociale (les autres variables sociales étant fixées à leur moyenne).

Les éléments “floutés” correspondent aux effets non-significatifs.

7.3.2.1 Group size

Model results
Lin : interaction linéaire
Qua : interaction quadratique

##    social trait            w interaction     y   ci-   ci+    p sign
## 1   Group  Mass           w1         Lin  0.01 -0.03  0.05 0.40     
## 4   Group  Mass           w1         Qua  0.04  0.02  0.05 0.00  ***
## 7   Group Tibia           w1         Lin -0.03 -0.07  0.01 0.07    .
## 10  Group Tibia           w1         Qua -0.03 -0.04 -0.01 0.00  ***
## 13  Group  Mass           w2         Lin -0.03 -0.06  0.00 0.03    *
## 16  Group  Mass           w2         Qua  0.00 -0.01  0.02 0.49     
## 19  Group Tibia           w2         Lin  0.03 -0.01  0.06 0.06    .
## 22  Group Tibia           w2         Qua  0.00 -0.01  0.01 0.50     
## 25  Group  Mass           w3         Lin -0.05 -0.10  0.00 0.04    *
## 28  Group  Mass           w3         Qua -0.01 -0.03  0.02 0.33     
## 31  Group Tibia           w3         Lin  0.03 -0.03  0.09 0.19     
## 34  Group Tibia           w3         Qua -0.01 -0.04  0.02 0.20     
## 37  Group  Mass           w4         Lin -0.02 -0.08  0.04 0.28     
## 40  Group  Mass           w4         Qua  0.03  0.00  0.06 0.02    *
## 43  Group Tibia           w4         Lin  0.01 -0.05  0.08 0.46     
## 46  Group Tibia           w4         Qua -0.02 -0.05  0.01 0.06    .
## 49  Group  Mass           w5         Lin -0.04 -0.10  0.03 0.12     
## 52  Group  Mass           w5         Qua  0.02 -0.02  0.05 0.15     
## 55  Group Tibia           w5         Lin  0.04 -0.03  0.11 0.12     
## 58  Group Tibia           w5         Qua -0.01 -0.05  0.02 0.24     
## 61  Group  Mass LRS_relative         Lin -0.08 -0.19  0.02 0.06    .
## 64  Group  Mass LRS_relative         Qua  0.08  0.02  0.16 0.00   **
## 67  Group Tibia LRS_relative         Lin  0.00 -0.15  0.15 0.50     
## 70  Group Tibia LRS_relative         Qua -0.07 -0.16  0.02 0.06    .

Plot

7.3.2.2 Litter size

Model results
Lin : interaction linéaire
Qua : interaction quadratique

##    social trait            w interaction     y   ci-  ci+    p sign
## 2  Litter  Mass           w1         Lin  0.00 -0.07 0.08 0.48     
## 5  Litter  Mass           w1         Qua  0.05  0.01 0.09 0.01    *
## 8  Litter Tibia           w1         Lin  0.00 -0.06 0.05 0.47     
## 11 Litter Tibia           w1         Qua -0.03 -0.07 0.00 0.01    *
## 14 Litter  Mass           w2         Lin -0.01 -0.09 0.07 0.43     
## 17 Litter  Mass           w2         Qua  0.02 -0.03 0.06 0.18     
## 20 Litter Tibia           w2         Lin -0.05 -0.11 0.02 0.07    .
## 23 Litter Tibia           w2         Qua  0.01 -0.03 0.05 0.43     
## 26 Litter  Mass           w3         Lin -0.11 -0.28 0.07 0.14     
## 29 Litter  Mass           w3         Qua -0.03 -0.12 0.06 0.25     
## 32 Litter Tibia           w3         Lin  0.11 -0.03 0.26 0.07    .
## 35 Litter Tibia           w3         Qua  0.06 -0.02 0.15 0.07    .
## 38 Litter  Mass           w4         Lin  0.02 -0.14 0.17 0.43     
## 41 Litter  Mass           w4         Qua -0.03 -0.13 0.08 0.24     
## 44 Litter Tibia           w4         Lin  0.06 -0.10 0.20 0.20     
## 47 Litter Tibia           w4         Qua  0.06 -0.05 0.15 0.12     
## 50 Litter  Mass           w5         Lin  0.14 -0.03 0.30 0.06    .
## 53 Litter  Mass           w5         Qua  0.05 -0.05 0.13 0.14     
## 56 Litter Tibia           w5         Lin -0.16 -0.33 0.00 0.02    *
## 59 Litter Tibia           w5         Qua  0.04 -0.07 0.15 0.23     
## 62 Litter  Mass LRS_relative         Lin  0.23 -0.09 0.59 0.10    .
## 65 Litter  Mass LRS_relative         Qua  0.12 -0.05 0.33 0.09    .
## 68 Litter Tibia LRS_relative         Lin -0.25 -0.56 0.04 0.05    *
## 71 Litter Tibia LRS_relative         Qua  0.05 -0.09 0.20 0.23

Plot

7.3.2.3 Number of helpers

Model results
Lin : interaction linéaire
Qua : interaction quadratique

##     social trait            w interaction     y   ci-   ci+    p sign
## 3  Helpers  Mass           w1         Lin -0.01 -0.14  0.15 0.47     
## 6  Helpers  Mass           w1         Qua -0.05 -0.16  0.05 0.20     
## 9  Helpers Tibia           w1         Lin  0.01 -0.13  0.15 0.46     
## 12 Helpers Tibia           w1         Qua  0.05 -0.07  0.16 0.20     
## 15 Helpers  Mass           w2         Lin -0.14 -0.27 -0.01 0.02    *
## 18 Helpers  Mass           w2         Qua  0.09 -0.01  0.18 0.04    *
## 21 Helpers Tibia           w2         Lin  0.10 -0.01  0.21 0.04    *
## 24 Helpers Tibia           w2         Qua -0.10 -0.18 -0.03 0.00   **
## 27 Helpers  Mass           w3         Lin  0.17 -0.05  0.41 0.07    .
## 30 Helpers  Mass           w3         Qua -0.06 -0.22  0.08 0.20     
## 33 Helpers Tibia           w3         Lin  0.00 -0.22  0.22 0.48     
## 36 Helpers Tibia           w3         Qua -0.10 -0.28  0.07 0.12     
## 39 Helpers  Mass           w4         Lin  0.15 -0.10  0.40 0.11     
## 42 Helpers  Mass           w4         Qua -0.09 -0.26  0.06 0.12     
## 45 Helpers Tibia           w4         Lin -0.11 -0.33  0.10 0.16     
## 48 Helpers Tibia           w4         Qua  0.07 -0.10  0.23 0.16     
## 51 Helpers  Mass           w5         Lin  0.25  0.00  0.50 0.02    *
## 54 Helpers  Mass           w5         Qua -0.10 -0.23  0.04 0.08    .
## 57 Helpers Tibia           w5         Lin  0.00 -0.24  0.27 0.47     
## 60 Helpers Tibia           w5         Qua -0.03 -0.20  0.17 0.35     
## 63 Helpers  Mass LRS_relative         Lin  0.36 -0.04  0.82 0.05    *
## 66 Helpers  Mass LRS_relative         Qua -0.15 -0.55  0.19 0.24     
## 69 Helpers Tibia LRS_relative         Lin  0.07 -0.47  0.55 0.38     
## 72 Helpers Tibia LRS_relative         Qua -0.23 -0.88  0.38 0.24

Plot

SOCIAL AND SELECTION TEMPORAL VARIATION IN THE ALPINE MARMOT OF LA GRANDE SASSIERE

Pierre-Alexandre QUITTET

11-07-2025

SOMMAIRE

1 VARIATION TEMPORELLE DES VARIABLES SOCIALES

1.1 Taille de groupe et taille de portée

1.2 Helpers

1.3 Corrélation entre les variables sociales

2 PROXY DE FITNESS

2.1 Proxy de fitness et scaling

2.2 Distribution de la fitness absolue

2.3 Compute de la fitness relative

2.4 Corrélation entre les proxy de fitness

3 SELECTION GLOBALE

3.1 Fonctions

3.1.1 compute_selection()

3.1.2 plot_compute_selection()

3.1.3 plot_overall_fitness_trait()

3.2 Selection differential

3.2.1 Visualisation des différentiels de sélection

3.3 Gradient de sélection

3.3.1 Visualisation des gradients de sélection

4 VARIATION SPATIALE DE LA FITNESS ET PHENOLOGIE

5 VARIATION TEMPORELLE DE LA SELECTION

5.1 Fonctions

5.1.1 compute_selection_across_time()

5.1.2 plot_compute_selection_across_time()

5.2 Différentiel de sélection

5.3 Gradient de sélection

6 FITNESS ET CONTEXTE SOCIALE

6.1 Fonctions

6.1.1 partial_effects_social()

6.1.2 plot_partial_effects_social()

6.2 Fitness across social

7 SELECTION ET CONTEXTE SOCIAL

7.1 Fonctions

7.1.1 compute_boot_model_results_social()

7.1.2 plot_social_model_result()

7.1.3 summary_table()

7.1.4 compute_selection_social_level()

7.1.6 plot_3d()

7.1.5 plot_compute_selection_social_level()

7.2 Selection differential across social

7.2.1 Significance of social effect on selection

7.2.2 Selection across social level

7.2.2.1 Group size

7.2.2.2 Litter size

7.2.2.3 Number of helpers

7.3 Selection gradient across social

7.3.1 Significance of social effect on selection

7.3.2 Selection across social level

7.3.2.1 Group size

7.3.2.2 Litter size

7.3.2.3 Number of helpers