Article 2 - Infanticide
Pierre-Alexandre QUITTET
2026-01-29
CHARGEMENT DES PACKAGES ET DES DONNEES
# PACKAGES
library(tidyverse)
# LOAD DATA FROM SQL DATA BASE
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/statut_dominance_raw.rds") # table statut social
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/pedigree_data.rds") # pédigrée (pour portée)
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/pheno.rds") # morpho
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/individu_bdd.rds") # informations individuelles
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/territoire_ref.rds") # infos territoire
# AJOUT SEXE (à ajouter avec BDD)
index_temp <- match(statut_dominance_raw$individual_id, individu_bdd$individu_id)
statut_dominance_raw$sexe <- individu_bdd$sexe[index_temp]
# REMOVE BLANK SPACE
statut_dominance_raw$dispersant <-
gsub(pattern = " ", replacement = "", statut_dominance_raw$dispersant)
# RENAME LEVELS OF DISPERSANT
statut_dominance <-
statut_dominance_raw |>
mutate(dispersant = case_when(
grepl(pattern = "sortant", x = dispersant) ~ "sortant",
grepl(pattern = "transient", x = dispersant) ~ "transient",
TRUE ~ dispersant
))## individual_id territoire_id statut_social_id dispersant type_attribution annee territoire taille_portee sexe
## 1 15 749 2 <NA> deduite 1989 c NA f
## 2 16 749 2 <NA> deduite 1989 c NA m
## 3 1 2 2 <NA> capture 1990 b 3 f
## 4 11 1 2 <NA> capture 1990 a 0 m
## 5 12 1 2 <NA> capture 1990 a 0 m
## 6 13 1 2 <NA> capture 1990 a 0 m
COMPLETION DES SEQUENCES DE PRESENCE DES INDIVIDUS DOMINANTS
Les données de présence de la table SQL statut_dominance
(contenant les présences par territoire de chaque indiviu, via capture
ou observation) sont probablement incomplètes et quelques données
peuvent probablement être récupérées. Je vois 2 cas de
figures où je pourrais potentiellement récupérer de l’information :
Reconstruction avec les informations du pédigrée - Individu est noté parent dans le pédigrée social (i.e. dominant) à une année \(t\), mais n’est pas noté comme dominant à l’année \(t\) dans la table
statut_dominanceReconstruction avec les séquences \(t\) à \(t+x\) - Individu a été capturé ou aperçu dominant sur un territoire donné de \(t\) à \(t+x\) et ajouter une ligne pour toutes les fois où il n’a pas été noté présent entre \(t\) et \(t+x\)
Reconstruction avec les informations du pédigrée
Pour voir le code associer, cliquer sur le bouton
showà droite.
# FILTRER LES DOMINANCE
statut_dominance_filtered <-
statut_dominance |>
filter(!duplicated(paste0(individual_id, "_", territoire_id, "_", statut_social_id))) |>
arrange(individual_id, territoire_id) |>
filter(statut_social_id == 1)
# CREER UN ID UNIQUE TERRITOIRE_ID (i.e. annee x territoire) x INDIVIDU ID
statut_key <- paste0(statut_dominance_filtered$individual_id, "_", statut_dominance_filtered$territoire_id)
# MERES MANQUANTE PAR TERRITOIRE ID
meres_manquantes <-
pedigree_data |>
filter(!is.na(mother)) |>
mutate(mere_territoire_key = paste0(mother, "_", territoire_id)) |>
filter(!mere_territoire_key %in% statut_key) |>
dplyr::select(individual_id = mother, territoire_id) |>
distinct()
# PERES MANQUANTE PAR TERRITOIRE ID
peres_manquants <- pedigree_data |>
filter(!is.na(father)) |>
mutate(pere_territoire_key = paste0(father, "_", territoire_id)) |>
filter(!pere_territoire_key %in% statut_key) |>
dplyr::select(individual_id = father, territoire_id) |>
distinct()
# COMBINER LES PARENTS MANQUANTS
parents_manquants <- rbind(meres_manquantes, peres_manquants) |>
distinct()
# GET LES INFORMATIONS DU TERRITOIRE
parents_manquants <- parents_manquants |>
left_join(territoire_ref, by = "territoire_id")
# GET SEXE
sexe_ref <- statut_dominance |>
filter(!duplicated(individual_id)) |>
dplyr::select(individual_id, sexe) |>
mutate(individual_id = as.character(individual_id))
parents_manquants <- parents_manquants |>
left_join(sexe_ref, by = "individual_id")
# NOUVELLES LIGNES A AJOUTER
nouvelles_lignes_pedigree <-
parents_manquants |>
filter(!is.na(territoire_id), !is.na(annee_territoire)) |>
mutate(
statut_social_id = 1,
dispersant = NA,
type_attribution = "deduite_pedigree",
taille_portee = NA,
annee = annee_territoire
) |>
dplyr::select(individual_id, territoire_id, statut_social_id, dispersant, type_attribution, annee, territoire, sexe, taille_portee)
# AJOUTER AUX DATAFRAME
statut_dominance <- rbind(
statut_dominance,
nouvelles_lignes_pedigree
) |>
filter(!duplicated(paste0(individual_id, "_", territoire_id, "_", statut_social_id))) |>
arrange(individual_id, territoire_id)Statistiques
Nombre de lignes (individual_id × territoire_id) ajoutées : 132
Nombre de territoire_id uniques concernés: 97
Reconstruction avec les séquences \(t\) à \(t+x\)
Pour voir le code associer, cliquer sur le bouton
showà droite.
# OBTENIR PAR TERRITOIRE, LES SERIES TEMPORELLES DE CHAQUE INDIVIDU
statut_dominance_filtered <-
statut_dominance_filtered |>
mutate(key = paste0(individual_id, territoire, annee))
# GET GUESSABLE SEQUENCES
statut_dominance_filtered_completed <-
statut_dominance_filtered |>
group_by(individual_id, territoire) |>
reframe(annee = min(annee):max(annee),
key = paste0(individual_id, territoire, annee)) |>
ungroup() |>
filter(!key %in% statut_dominance_filtered$key) |>
dplyr::select(-key)
# MERGE LES HISTOIRES RECONSTRUITES AU DATAFRAME
df_temp <- data.frame(
individual_id = statut_dominance_filtered_completed$individual_id,
territoire_id = NA,
statut_social_id = 1,
dispersant = NA,
type_attribution = "deduite_sequence",
annee = statut_dominance_filtered_completed$annee,
territoire = statut_dominance_filtered_completed$territoire,
taille_portee = NA,
sexe = NA
)
statut_dominance_completed <- rbind(statut_dominance, df_temp)
# COMPLETER LA COLONNE TERRITOIRE ID
index_temp <- match(
x = paste0(statut_dominance_completed$territoire, statut_dominance_completed$annee),
table = paste0(territoire_ref$territoire, territoire_ref$annee)
)
statut_dominance_completed$territoire_id <- territoire_ref$territoire_id[index_temp]
# REFILTRER LES DUPLICATIONS
statut_dominance <-
statut_dominance_completed |>
filter(!duplicated(paste0(individual_id, "_", territoire_id, "_", statut_social_id))) |>
arrange(individual_id, territoire_id)
# AJOUT SEXE
index_temp <- match(statut_dominance$individual_id, individu_bdd$individu_id)
statut_dominance$sexe <- individu_bdd$sexe[index_temp]
# SAVE
save(statut_dominance, file = "~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/statut_dominance.rds")Statistiques
Nombre de lignes (individual_id × territoire_id) ajoutées : 60
Nombre de territoire_id (i.e. territoire x) uniques concernés: 47
STATISTIQUES DESCRIPTIVES
On commence par extraire le couple de dominant par territoire par année depuis 1990.
# FILTER DATAFRAME TO KEEP ONLY DOMINANT
individu_dominant <-
statut_dominance |>
filter(statut_social_id %in% c(1)) |> # 1 = dominant
arrange(territoire, annee) |>
group_by(territoire, annee) |>
mutate(individual_id = list(unique(individual_id)),
sexe = list(sexe))
# GET COUPLE
individu_dominant$couple <- sapply(individu_dominant$individual_id, paste, collapse = "-")
# QUICK VIEW
head(individu_dominant)## # A tibble: 6 × 11
## # Groups: territoire, annee [3]
## individual_id territoire_id statut_social_id dispersant type_attribution annee territoire taille_portee sexe litter_id couple
## <list> <int> <dbl> <chr> <chr> <int> <chr> <int> <list> <chr> <chr>
## 1 <chr [1]> 1 1 <NA> capture 1990 a 0 <chr [1]> <NA> 7
## 2 <chr [3]> 6 1 <NA> capture 1991 a 4 <chr [3]> <NA> 32-4-7
## 3 <chr [3]> 6 1 entrant deduite 1991 a 4 <chr [3]> <NA> 32-4-7
## 4 <chr [3]> 6 1 sortant capture 1991 a 4 <chr [3]> <NA> 32-4-7
## 5 <chr [2]> 15 1 <NA> deduite 1992 a 5 <chr [2]> 32_4_15 32-4
## 6 <chr [2]> 15 1 <NA> capture 1992 a 5 <chr [2]> 32_4_15 32-4
Nombre de couple unique
## [1] 447
Nombre de territoires x années suivi
IDENTIFIER LES CAS POTENTIELS D’INFANTICIDE
L’infanticide chez la marmotte Alpine
Généralités sur l’infanticide chez la marmotte
- Changement de dominance mâle comme femelles entraine l’infanticide
- Tous les marmottons et certains des yearlings meurent en cas de changement de dominance
- Plus la saison avance, moins les changements de dominances sont fréquents (Rebecca)
Dans certains cas, les yearlings du sexe opposant au dominant qui saute restent. A voir si cela dépend du poids ? De la date du changement de dominance ? S’ils sont assez gros pour être considérés comme des reproducteurs potentiels…
Occurrence de l’infanticide
- Avant l’arrivée sur le terrain (de l’émergence au début de terrain) [yearlings]
- Pendant le terrain [marmottons et yearlings]
- Après le terrain (de l’après terrain à l’entrée en hibernation) [marmottons et yearlings]
Quantifier les changements de dominance
Dans un premier temps, nous devons caractériser les changements de dominance dont dépendent presque totalement les occurrences d’infanticide. D’un point de vue opérationnel, je considérerai un changement de dominance lorsque :
- Le couple de dominant à \(t\) est différent de celui à \(t-1\)
- Le couple de dominant passe de 1 à 2 (i.e. un des deux n’étaient pas connus)
NOTE DE CODE
Il faut gérer les cas de figures comme à la ligne 8-9 car si “coexistence” de e.g. >3 dominants dans les données (i.e. dominance a été détectée pendant le terrain), alors automatiquement l’année \(t+1\) sera > considérée comme “changement de dominance” alors que le couple était déjà formée à \(t\).Je rajoute donc la condition “si à l’année \(t\) il y a >3 dominants (
n_dom> 3) ET que le couple à \(t+1\) est un sous-ensemble de celui de \(t\), alors ce n’est pas un changement”. Pour cela, j’utilise le packagestringr.
Changement de dominance
# COMPUTE LE CHANGEMENT DE DOMINANCE PAR ANNEE PAR TERRITOIRE
individu_dominant <-
individu_dominant |>
filter(!duplicated(territoire_id)) |>
complete(territoire, annee, fill = list(NA)) |> # on complète les années manquantes par NAs sinon génère artificellement des changements de dominances
arrange(territoire, annee) |>
group_by(territoire) |>
mutate(
cur = strsplit(couple, "-"),
before = strsplit(lag(couple), "-"),
is_na = is.na(couple),
# Tous les dominants de t+1 étaient déjà présents à t
inclus = !is.na(lag(couple)) &
mapply(function(x, y) all(y %in% x), before, cur),
changement_dom = if_else(is_na, NA, if_else(couple != lag(couple) &
!inclus, T, F, missing = F))
) |>
dplyr::select(-cur, -before, -inclus, -is_na) |>
as.data.frame()
La somme totale de changement de dominance :
## [1] 283
Par années
Changement de dominance par sexe
# IDENTIFIER LE TYPE DE CHANGEMENT DE DOMINANCE
# INFORMATIONS PAR SEXE PAR CHANGEMENT DE DOMINANCE
individu_dominant_sexe <-
individu_dominant |>
arrange(territoire, annee) |>
group_by(territoire) |>
mutate(
# Calculer AVANT de passer en rowwise
sexes_couple_avant = lag(sexe),
individual_id_avant = lag(individual_id)
) |>
ungroup()
# IDENTIFIER TYPE DE CHANGEMENT (couple entier, male, femelle)
individu_dominant_sexe <-
individu_dominant_sexe |>
rowwise() |> # permet d'appliquer le mutate par ligne
mutate(
type_changement = if_else(
changement_dom == TRUE,
{
# Extraire les individus et sexes avant et après
ids_avant = unlist(individual_id_avant)
ids_apres = unlist(individual_id)
sexes_av = unlist(sexes_couple_avant)
sexes_ap = unlist(sexe)
# Gérer le cas où ids_avant est NULL (première année)
if (is.null(ids_avant) || length(ids_avant) == 0) {
NA_character_
} else {
# Séparer par sexe
males_avant = ids_avant[sexes_av == "m"]
males_apres = ids_apres[sexes_ap == "m"]
femelles_avant = ids_avant[sexes_av == "f"]
femelles_apres = ids_apres[sexes_ap == "f"]
# Vérifier si chaque sexe a changé
changement_male = !all(males_apres %in% males_avant)
changement_femelle = !all(femelles_apres %in% femelles_avant)
# Classifier le type
if (changement_male & changement_femelle) {
"couple_entier"
} else if (changement_male) {
"male"
} else if (changement_femelle) {
"femelle"
} else {
NA_character_
}
}
},
NA_character_
)
) |>
ungroup()## # A tibble: 3 × 3
## type_changement n proportion
## <chr> <int> <dbl>
## 1 couple_entier 66 0.236
## 2 femelle 92 0.329
## 3 male 122 0.436
Quantifier l’infanticide
Détecter l’infanticide
- Cas 1 : Présence de juvéniles + observation directe de changement de dominance
- Cas 2 : Présence de marmottons à \(t\) + changement de dominance à \(t+1\)
- Cas 3 : Présence de yearling à \(t\) + changement de dominance à \(t+1\)
Il faut distinguer le cas de figure où la portée manquante n’est pas celle du nouveau dominant car si le changement de dominance a lieu tôt dans la saison et que le nouveau couple a eu le temps d’avoir une portée, alors on ne s’attend pas à l’absence de portée à \(t\).
Potentiels effets confondants
- Absence de juvéniles (marmottons + yearlings) meurent pendant l’hibernation + changement de dominance pré-terrain
- Portée manquée + changement de dominance pré-terrain
- Bruce effect / pregnancy block - infanticide pendant la gestation par la propre mère en cas de changement de mâle dominant (mais moins important car ces individus n’apparaissent pas dans la CMR)
On considérera donc comme potentiel infanticide les changements de dominance associés à la présence de marmotton ou de yearlings. Il est également nécessaire de tester si la portée a bien changée… car si c’est la portée du nouveau couple dominant, alors on ne s’attend pas à un infanticide.
Présence de marmottons
Je vais d’abord construire la variable presence_pup par
année*territoire (variable territoire_id).
# CONSTRUCT VARIABLE <PRESENCE OF LITTER>
# ... Si le territoire ID est dans le pédigrée, alors une portée cette année...
individu_dominant_sexe$presence_pup <-
ifelse(individu_dominant_sexe$territoire_id %in% pedigree_data$territoire_id,
T,
F)
# IDENTITE DE LA PORTEE (pour prendre en compte le cas de figure où la portée change)
# Créer une table portée unique par territoire et année
litter_table <- pedigree_data |>
filter(!is.na(litter_id)) |>
mutate(cohort = as.integer(cohort)) |>
distinct(territoire_id, cohort)
# Merge avec individu_dominant_sexe
individu_dominant_sexe <-
individu_dominant_sexe |>
left_join(litter_table,
by = c("territoire_id" = "territoire_id",
"annee" = "cohort"))
individu_dominant_sexe <-
individu_dominant_sexe |>
rename(id_litter = litter_id)
# Variables <changement_litter> qui teste si la portée de t est différente de celle de t-1
individu_dominant_sexe <-
individu_dominant_sexe |>
arrange(territoire, annee) |>
group_by(territoire) |>
mutate(changement_litter = case_when(
is.na(id_litter) & is.na(lag(id_litter)) ~ FALSE, # aucune portée ni avant ni après
is.na(id_litter) | is.na(lag(id_litter)) ~ TRUE, # portée apparue ou disparue
id_litter != lag(id_litter) ~ TRUE, # portée différente
TRUE ~ FALSE
))
Présence de yearlings
Maintenant, on doit ajouter la variable
presence_yearling, aussi victime de l’infanticide en cas de
changement de dominance.
# IDENTIFIER LES YEARLINGS PAR ANNEE PAR TERRITOIRE
# GET EMERGENCE DATE
statut_dominance <-
statut_dominance |>
full_join(pheno[c("individual_id", "age_obs", "annee_emergence")],
by = "individual_id",
multiple = "first") |>
rename(annee_emergence_pheno = annee_emergence)
statut_dominance <-
statut_dominance |>
full_join(pedigree_data[c("individual_id", "cohort")],
by = "individual_id",
multiple = "first") |>
rename(annee_emergence_pedigree = cohort)
statut_dominance$annee_emergence_pheno <- as.numeric(as.character(statut_dominance$annee_emergence_pheno))
statut_dominance <- statut_dominance |>
mutate(
annee_emergence = coalesce( # find the first non-missing value
annee_emergence_pedigree,
annee_emergence_pheno
)
) |>
dplyr::select(-c(annee_emergence_pheno, annee_emergence_pedigree))
# COMPUTE AGE TO IDENTIFY YEARLINGS
statut_dominance$age_annee <-
as.numeric(as.character(statut_dominance$annee)) - as.numeric(as.character(statut_dominance$annee_emergence)) +
statut_dominance$age_obs
# GET ALL YEARLINGS
yearling_df <- statut_dominance |>
filter(age_annee == 1) |>
distinct(territoire_id)
# ADD THE VARIABLE presence_yearling
yearling_df <- statut_dominance |>
filter(age_annee == 1) |>
distinct(individual_id, .keep_all = TRUE)
individu_dominant_sexe$presence_yearling <-
individu_dominant_sexe$territoire_id %in% yearling_df$territoire_id
On ne conserve qu’une seule ligne par territoire-année
# ON NE GARDE QU'UNE MESURE PAR TERRITOIRE ID
dominance_per_territory <-
individu_dominant_sexe |>
filter(!duplicated(territoire_id))
Présence de juvéniles
Je construis la variable presence_juvenile qui regroupe
la présence de marmottons OU de yearlings.
dominance_per_territory$presence_juvenile <- dominance_per_territory$presence_pup | dominance_per_territory$presence_yearling
Potentiel infanticide
Maintenant, on construit la variable
potential_infanticide où les potentiels infanticides sont
tous les cas de figure où un changement de dominance coincide avec la
présence de juvéniles.
dominance_per_territory$potential_infanticide <-
dominance_per_territory$changement_dom &
(dominance_per_territory$presence_juvenile) &
dominance_per_territory$changement_litter
Statistiques descriptives
Nombre total de territoire-année avec au moins un marmotton
## [1] "Somme : 633"## [1] "Proportion : 0.67"
Nombre total de territoire-année avec au moins un yearling
## [1] "Somme : 345"## [1] "Proportion : 0.37"
Nombre total de territoire-année avec au moins un juvénile (marmotton ou yearling)
## [1] "Somme : 691"## [1] "Proportion : 0.73"
Nombre total de territoire-année avec potentiel infanticide
## [1] "Somme : 154"## [1] "Proportion : 0.16"
Répartition des types de changements en cas de potentiels infanticides :
## TYPE DE CHANGEMENT EN CAS DE POTENTIEL INFANTICIDE
## couple_entier femelle male
## 30 50 71
HISTOIRE DE CAPTURE
On considère ici qu’il n’y a pas d’incertitude lié à la connaissance du changement de dominance. On définit les histoires de capture pour le modèle CMR en intégrant 6 évènements :
- 0 : non-capturé sans infos de changement de dominance sur le territoire
- 1 : non-capturé avec changement du dominant de sexe opposé
- 2 : non-capturé avec changement des deux dominants ou du dominant du même sexe
- 3 : capturé sans changement de dominance sur le territoire
- 4 : capturé avec changement du dominant de sexe opposé
- 5 : capturé avec changement des deux dominants ou du dominant du même sexe
L’idée est de croiser une matrice de capture
“simple” contenant 1 pour “capturé” et 0 pour
non-capturé avec une matrice de changement de
dominance contenant des 0 (pas d’information de
changement de dominance), 1 (changement du dominant de sexe
opposé) et 2 (changement des deux dominants ou du dominant
du même sexe). L’intersection des deux matrices donnera une matrice
contenant les évènements pour le modèle (e.g. capturé 1 x
‘perte des deux dominants’ 3 == 5).
Matrice de capture simple
Chargement de la matrice
# LOAD CAPTURE MATRIX FROM histoire_capture_bdd.R script
load("C:/Users/quittet/Documents/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/mat_capture_bdd_infanticide.rds")
# LOAD INFORMATION SUR LES TERRITOIRES
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/territoire_ref.rds")
# AJOUT KEY POUR LES TERRITOIRES POUR RETROUVER L'EQUIVALENCE TERRITOIRE & TERRITOIRE_ID
territoire_ref$key <- paste0(territoire_ref$territoire, "_", territoire_ref$annee_territoire)
# LOAD PEDIGREE
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/pedigree_data.rds")## A1990 A1991 A1992 A1993 A1994 A1995 A1996 A1997 A1998 A1999 A2000 A2001 A2002 A2003 A2004 A2005 A2006 A2007 A2008 A2009 A2010 A2011 A2012 A2013 A2014 A2015 A2016 A2017 A2018 A2019 A2020 A2021
## 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 10 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 100 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 1001 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 1002 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## A2022 A2023 A2024 A2025 classe sexe individual_id
## 1 0 0 0 0 yearling f 1
## 10 0 0 0 0 twoyearsold f 10
## 100 0 0 0 0 adult m 100
## 1000 0 0 0 0 pup f 1000
## 1001 0 0 0 0 pup f 1001
## 1002 0 0 0 0 pup f 1002
Recodage et filtrage
Pour le moment, par simplicité, je ne prends que les individus capturés à la naissance. J’ajouterai les capturés yearlings a posteriori.
# REMPLACE TOUTES LES VALEURS NON NULLES PAR DES 1 (pour "capturé")
mat_capture <- mat_capture_bdd_infanticide |>
mutate(across(1:36, ~ case_when(
.x == 0 ~ 0L,
.x %in% c(1, 2, 3) ~ 1L
)))
# FILTER INDIVIDUALS CAPTURED PUP OR YEARLING
mat_capture <-
mat_capture |>
filter(classe %in% c("pup")) # seuls les pups pour le moment
# AJOUT DE L'INFORMATION TERRITOIRE
bool <- match(mat_capture$individual_id, pedigree_data$individual_id)
mat_capture$territoire_id_naissance <- pedigree_data$territoire_id[bool]
bool <- match(mat_capture$territoire_id_naissance, territoire_ref$territoire_id)
mat_capture$territoire <- territoire_ref$territoire[bool]
Vérification : la matrice ne contient que des 0 et des 1
##
## 0 1
## 60461 3979
Vérification : la matrice ne contient que des individus capturés marmottons
##
## pup
## 1790
Vérification : présence de NAs dans les territoires ?
## [1] 0
Matrice de changement de dominance
On construit la matrice de changement de dominance.
# FORMAT LONG DE LA MATRICE DE CAPTURE
mat_capture_long <-
mat_capture |>
as.data.frame() |>
pivot_longer(
cols = 1:36,
names_to = "annee_col",
values_to = "capture"
) |>
mutate(
annee = as.integer(sub("A", "", annee_col))
) |>
dplyr::select(-annee_col)
# AJOUT TERRITOIRE ID DANS LA MATRICE DE CAPTURE
mat_capture_long$key <- paste0(mat_capture_long$territoire, "_", mat_capture_long$annee)
mat_capture_long <- left_join(mat_capture_long, territoire_ref[c("key", "territoire_id")], by = "key")
# AJOUT DES INFORMATIONS DE DOMINANCE
mat_capture_long <-
left_join(mat_capture_long,
dominance_per_territory[c("territoire_id",
"type_changement")],
by = "territoire_id",
multiple = "first")
Vérification : Quelles fréquences des différents types de changements de dominance ?
##
## couple_entier femelle male <NA>
## 3274 4896 6451 49819Beaucoup de NAs car toute la matrice a été mise en format long, donc beaucoup de 0 (non-capture) car tous les territoires sont comptabilisés de 1990 à 2025 indépendement de leur existence à \(t\). Je mets 0 (évènement) à ces cas de figures.
Matrice des évènements
Il faut considérer les évènements à partir de la naissance uniquement. Je mets des 0 sur tout ce qui se passe AVANT la première capture. Je mets également des 0 APRES les 2 ans de chaque individu (censure des histoires de capture).
################################################################################
# #
# GENERATION DES EVENEMENTS #
# #
################################################################################
# GENERATION DES EVENEMENTS BRUTES
mat_capture_long <-
mat_capture_long |>
# Identifier la première capture par individu
group_by(individual_id) |>
mutate(
premiere_capture = min(annee[capture == 1], na.rm = TRUE),
) |>
ungroup() |>
# Recodage
mutate(
event = case_when(
#AVANT LA PREMIERE CAPTURE
!is.na(premiere_capture) & annee < premiere_capture ~ 0L,
# NON-CAPTURE
# Non-capturé — changement du dominant de même sexe ou des deux
capture == 0 & type_changement == "couple_entier" ~ 2L,
capture == 0 & type_changement == "male" & sexe == "m" ~ 2L,
capture == 0 & type_changement == "femelle" & sexe == "f" ~ 2L,
# Non-capturé — changement du dominant de sexe opposé
capture == 0 & type_changement == "male" & sexe == "f" ~ 1L,
capture == 0 & type_changement == "femelle" & sexe == "m" ~ 1L,
# Non-capturé — sans changement
capture == 0 & is.na(type_changement) ~ 0L,
# CAPTURE
# Capturé — changement du dominant de même sexe ou des deux
capture == 1 & type_changement == "couple_entier" ~ 5L,
capture == 1 & type_changement == "male" & sexe == "m" ~ 5L,
capture == 1 & type_changement == "femelle" & sexe == "f" ~ 5L,
# Capturé — changement du dominant de sexe opposé
capture == 1 & type_changement == "male" & sexe == "f" ~ 4L,
capture == 1 & type_changement == "femelle" & sexe == "m" ~ 4L,
# Capturé — sans changement
capture == 1 & is.na(type_changement) ~ 3L,
TRUE ~ NA_integer_
)
) |>
dplyr::select(-premiere_capture) |>
as.data.frame()
################################################################################
# #
# ALIGNER LES CAPTURES SUR LA GAUCHE #
# #
################################################################################
#' ... Permet d'éviter de censurer les données, on utilisera un effet groupe "année" en covariable à la place pour modéliser les effets temps.
seuil_censure <- 10 # nombre d'évènements à conserver après la naissance
mat_capture_long <-
mat_capture_long |>
# Identifier la première capture et la borne de censure
group_by(individual_id) |>
mutate(
cohort = min(annee[capture == 1], na.rm = TRUE),
borne_censure = cohort + seuil_censure
) |>
ungroup() |>
# Booléen pour supprimer ligne
mutate(row_to_remove = ifelse(!is.na(borne_censure) &
((annee > borne_censure) | (annee < cohort)), T, F)) |>
filter(!row_to_remove) |>
# On centre les années sur la gauche
mutate(
annee = annee - (cohort-1990))Vérification : on est censé avoir 11 évènements de capture par individu
##
## 2 3 4 5 6 7 8 9 10 11
## 79 61 42 85 65 93 63 45 77 1180Quelques individus n’apparaissent que deux fois… Pas sûr de comprendre pourquoi. Problème à régler plus tard. Ils représentent moins de 5% des cas. Je remplace ces lignes par des 0 dans la suite du code.
Vérification : seulement 11 évènements
##
## 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000
## 1790 1790 1711 1650 1608 1523 1458 1365 1302 1257 1180
FORMATAGE
mat_capture_esurge <-
mat_capture_long |>
# Mise en format wide
dplyr::select(classe, sexe, territoire_id_naissance, individual_id, territoire, annee, event, cohort) |>
pivot_wider(
names_from = annee,
values_from = event,
names_prefix = "H:A",
values_fill = 0
) |>
dplyr::select(
starts_with("H:A"),
classe,
sexe,
famille = territoire_id_naissance,
individual_id,
cohort
) |>
as.data.frame()## H:A1990 H:A1991 H:A1992 H:A1993 H:A1994 H:A1995 H:A1996 H:A1997 H:A1998 H:A1999 H:A2000 classe sexe famille individual_id cohort
## 1 3 0 0 0 0 2 0 0 0 0 0 pup f 307 1000 2006
## 2 5 0 0 0 0 0 0 0 0 0 0 pup f 296 1001 2006
## 3 5 0 0 0 0 0 0 0 0 0 0 pup f 296 1002 2006
## 4 3 2 1 2 2 2 0 0 0 0 0 pup m 301 1003 2006
## 5 3 1 2 1 1 1 0 0 0 0 0 pup f 301 1004 2006
## 6 3 5 3 1 1 0 0 0 2 0 0 pup m 32 101 1993
PREPARATION DATA
On va considérer que les changements de dominance se passe à \(t+1\). Je vais donc décaler toutes les informations de changements de dominance de \(+1\).
Histoires de capture impossibles
Ces cas de figure sont :
- un juvénile (marmotton ou yearling) est recapturé (évènements 3, 4 ou 5) l’année après un changement du couple de dominant ou du dominent de même sexe (évènement 2/5)
- un marmotton est recapturé (évènements 3, 4 ou 5) après un changement de dominant de sexe opposé (évènement 4)
Je construis une fonction permettant d’extraire les évènements théoriquement impossible :
check_impossible_case <- function(data = mat_capture_esurge,
event = 5) {
to_check <-
data |>
rowwise() |>
mutate(
to_check = {
list_event <- unlist(across(starts_with("H:A")))
first_event <- which(list_event %in% event)[1]
last_pos <- length(list_event)
if (is.na(first_event)) {
FALSE
} else if (first_event == last_pos) {
FALSE
} else {
any(list_event[(first_event + 1):last_pos] %in% c(3, 4, 5))
}
}
) |>
filter(to_check) |>
as.data.frame()
}
# IDENTIFY IMPOSSIBLE CASES
changement_dom_5 <- check_impossible_case(data = mat_capture_esurge |>
dplyr::select(-c(4:11)), event = c(2, 5))
# Changement d'un des dom..
changement_dom_4 <- check_impossible_case(
data = mat_capture_esurge |>
dplyr::select(-c(3:11)),
event = 4)
total_cases <-
bind_rows(changement_dom_5, changement_dom_4) |>
distinct(individual_id, .keep_all = TRUE)On a 201 cas de figures où un marmotton est recapturé après un cas de figure où un infanticide était censé avoir lieu (133 avec perte de deux dominants / même sexe et 72 avec perte du dominant de sexe opposé en tant que marmotton).
Tous ces cas de figures sont à vérifier à la main.
# EXTRACT PROBLEMATIC CASES
to_check <- total_cases |>
dplyr::select(-c(to_check, classe)) |>
arrange(famille)
# REPARTITION DANS LES ANNEES
barplot(table(to_check$cohort), main = "Répartition des impossibilités dans les années", ylab = "n")
# SAVE
xlsx::write.xlsx(
to_check,
"~/marmot-quantitative-genetics/1_DATA_ANALYSIS/INFANTICIDE/to_check.xlsx",
row.names = FALSE
)
POUR LE MOMENT JE SUPPRIME CES CAS DE FIGURES PROBLEMATIQUES.
Nombre d’occurrences supprimés par évènement :
## AVANT SUPPRESSION
## 0 1 2 3 4 5
## 2619 448 752 2501 368 472## APRES SUPPRESSION
## 0 1 2 3 4 5
## 2468 412 715 2229 257 275
Décalage des histoires de captures en cas de changement de dominance
On considérera que les changements de dominance ont lieu à \(t+1\). Soit :
- Pour l’observation 1 (marmotton)
- si évènement 4 ou 5, je décale tous les évènements à obs2 et les modifie en 1 ou 2 (i.e. infanticide)
- Pour l’observation 2 (yearling)
- Si évènement 4, je vérifie si l’individu a été recapturé après, si oui, je mets 4 en obs3, sinon je mets 1
- Si évènement 5, alors je mets 2 en obs3
- Pour l’observation 3 (adulte) : on ignore l’info
sociale
- Si évènements 0, 1, 2 alors : 0 (non-capturé)
- Si évènements 3, 4, 5 : alors 3 (capturé)
# Sauvegarder les valeurs originales AVANT tout recodage
orig_1990 <- mat_capture_esurge$`H:A1990`
orig_1991 <- mat_capture_esurge$`H:A1991`
orig_1992 <- mat_capture_esurge$`H:A1992`
mat_capture_esurge <-
mat_capture_esurge |>
mutate(
################################################################################
# #
# OBSERVATION 1 (pup) #
# #
################################################################################
# T = 1990 (naissance) : toujours capturé sans changement
`H:A1990_new` = 3L,
################################################################################
# #
# OBSERVATION 2 (yearling) #
# #
################################################################################
# T = 1991 (yearling) : on décale le changement de dominance observé de t=1990
`H:A1991_new` = case_when(
orig_1990 == 5 ~ 2L, # perte des deux dominants / même sexe = infanticide à t+1
orig_1990 == 4 ~ 1L, # perte du dominant de sexe opposé = infanticide à t+1
orig_1991 %in% c(4L, 5L) ~ 3L, # on met 3 en obs2, et on reporte ces cas en obs3
TRUE ~ orig_1991 # on conserve tous les autres cas (i.e. les 3 et les 0)
),
################################################################################
# #
# OBSERVATION 3 (adult) #
# #
################################################################################
# T = 1992 : on décale le changement de dominance observé de t=1991
`H:A1992_new` = case_when(
orig_1991 == 5 ~ 2L, # perte des deux dominants / même sexe = infanticide à t+1
(orig_1991 == 4 & orig_1992 %in% c(0L, 1L, 2L)) ~ 1L,
(orig_1991 == 4 & orig_1992 %in% c(3L, 4L, 5L)) ~ 4L
)) |>
# ... Pour le reste des évènements, on précise simplement l'information de capture / non-capture
mutate(`H:A1992_new` = case_when(
(is.na(`H:A1992_new`) & orig_1992 %in% c(3, 4, 5)) ~ 3,
(is.na(`H:A1992_new`) & orig_1992 %in% c(0, 1, 2)) ~ 0,
TRUE ~ `H:A1992_new`
)) |>
# Supprimer les colonnes originales et renommer les nouvelles
dplyr::select(-`H:A1990`, -`H:A1991`, -`H:A1992`) |>
rename(
`H:A1990` = `H:A1990_new`,
`H:A1991` = `H:A1991_new`,
`H:A1992` = `H:A1992_new`
) |>
# Remettre dans le bon ordre
transmute(
`H:A1990`, `H:A1991`, `H:A1992`,
cohort, classe, sexe, famille, individual_id
)
Verification : on vérifie que des cas théoriquement impossibles ne sont pas présents après le recodage :
# IDENTIFY IMPOSSIBLE CASES
changement_dom_5 <- check_impossible_case(data = mat_capture_esurge, event = c(2, 5))
# Changement d'un des dom..
changement_dom_4 <- check_impossible_case(
data = mat_capture_esurge |>
dplyr::select(-c(3)),
event = 4)Nombre d’histoire impossible : 0
Etat mort après changement de dominance 5 ou 4 - marmotton
Je mets dans l’état mort l’ensemble des juvéniles après un infanticide (i.e. après un 1 ou 2).
mat_capture_esurge <-
mat_capture_esurge |>
mutate(`H:A1992` = if_else(
condition = `H:A1991` %in% c(1, 2),
true = 0,
false = `H:A1992`
))
Vérification : Quelles fréquences des différents évènements ?
## PUPS
## 0 1 2 3
## 489 184 280 636## YEARLINGS
## 0 1 2 3 4
## 1157 45 44 311 32
Vérification : quelles fréquences des différentes transitions d’observation ?
EXPORT DATA
# COLUMNS NAMES
mat_capture_esurge <-
mat_capture_esurge |>
rename(`$COV:classe` = classe,
`$COV:famille` = famille,
`$COV:year` = cohort
) |>
dplyr::select(-c(sexe, individual_id)) |>
as.data.frame()
# E SURGE
# .TXT
write.table(
mat_capture_esurge,
"~/marmot-quantitative-genetics/1_DATA_ANALYSIS/INFANTICIDE/mat_capture_infanticide.txt",
sep = "\t",
eol = "\n",
row.names = FALSE,
quote = FALSE
)
# XLSX
write.csv(
mat_capture_esurge,
"~/marmot-quantitative-genetics/1_DATA_ANALYSIS/INFANTICIDE/mat_capture_esurge.csv",
row.names = FALSE
)