# PACKAGES
library(tidyverse)
# LOAD DATA FROM SQL DATA BASE
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/statut_dominance.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$individual_id, individu_bdd$individu_id)
statut_dominance$sexe <- individu_bdd$sexe[index_temp]## individual_id territoire_id statut_social_id type_attribution annee territoire sexe
## 1 15 749 2 deduite 1989 c f
## 2 16 749 2 deduite 1989 c m
## 3 1 2 2 capture 1990 b f
## 4 11 1 2 capture 1990 a m
## 5 12 1 2 capture 1990 a m
## 6 13 1 2 capture 1990 a m
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 :
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\)
Individu est noté parent dans le pédigrée à une année \(t\), mais n’est pas noté comme dominant à
l’année \(t\) dans la table
statut_dominance
# 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) |>
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) |>
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)) |>
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,
type_attribution = "deduite_pedigree",
annee = annee_territoire
) |>
select(individual_id, territoire_id, statut_social_id, type_attribution, annee, territoire, sexe)
# 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)Vérification et statistiques
Nombre de lignes (individual_id × territoire_id) ajoutées : 117
Nombre de territoire_id uniques concernés: 84
# 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) |>
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,
type_attribution = "deduite_sequence",
annee = statut_dominance_filtered_completed$annee,
territoire = statut_dominance_filtered_completed$territoire,
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]
Vérification et statistiques
Nombre de lignes (individual_id × territoire_id) ajoutées : 60
Nombre de territoire_id uniques concernés: 47
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)) |>
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 × 8
## # Groups: territoire, annee [3]
## individual_id territoire_id statut_social_id type_attribution annee territoire sexe couple
## <list> <int> <dbl> <chr> <chr> <chr> <list> <chr>
## 1 <chr [1]> 1 1 capture 1990 a <chr [1]> 7
## 2 <chr [3]> 6 1 capture 1991 a <chr [3]> 32-4-7
## 3 <chr [3]> 6 1 deduite 1991 a <chr [3]> 32-4-7
## 4 <chr [3]> 6 1 capture 1991 a <chr [3]> 32-4-7
## 5 <chr [2]> 15 1 deduite 1992 a <chr [2]> 32-4
## 6 <chr [2]> 15 1 capture 1992 a <chr [2]> 32-4
Généralités sur l’infanticide chez la marmotte
Occurrence de l’infanticide
Dans un premier temps, nous devons caractériser les changements de dominance dont dépendent totalement les occurrences d’infanticide. D’un point de vue opérationnel, je considérerai un changement de dominance lorsque :
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.
# 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))
) |>
select(-cur, -before, -inclus, -is_na) |>
as.data.frame()
La somme totale de changement de dominance :
## [1] 283
Par années
# 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 × 4
## type_changement n proportion pourcentage
## <chr> <int> <dbl> <dbl>
## 1 couple_entier 59 0.211 21.1
## 2 femelle 95 0.339 33.9
## 3 male 126 0.45 45
Détecter l’infanticide
Potentiels effets confondants
On considérera donc comme potentiel infanticide les changements de dominance associés à la présence de marmotton ou de yearlings.
Je vais d’abord construire la variable presence_pup par
année*territoire.
# CONSTRUCT VARIABLE <PRESENCE OF LITTER>
# ... Si le territoire ID est dans le pédigrée, alors une portée cette année...
individu_dominant$presence_pup <-
ifelse(individu_dominant$territoire_id %in% pedigree_data$territoire_id,
T,
F)
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 <- statut_dominance |>
mutate(
annee_emergence = coalesce( # find the first non-missing value
annee_emergence_pedigree,
annee_emergence_pheno
)
) |>
select(-c(annee_emergence_pheno, annee_emergence_pedigree))
# COMPUTE AGE
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$presence_yearling <-
individu_dominant$territoire_id %in% yearling_df$territoire_id
# ON NE GARDE QU'UNE MESURE PAR TERRITOIRE ID
dominance_per_territory <-
individu_dominant |>
filter(!duplicated(territoire_id))
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
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)
Nombre total de territoire-année avec au moins un marmotton
## n_presence_pup prop_presence_pup
## 1 631 0.7211429
Nombre total de territoire-année avec au moins un yearling
## n_presence_pup prop_presence_pup
## 1 338 0.3862857
Nombre total de territoire-année avec au moins un juvénile
## n_presence_pup prop_presence_pup
## 1 687 0.7851429
Nombre total de territoire-année avec potentiel infanticide
## n_presence_pup prop_presence_pup
## 1 198 0.2262857