Étude de la démographie française du XIXe siècle à partir de données collaboratives de généalogie

Arthur Charpentier et Ewen Gallic
freakonometrics.hypotheses.org
egallic.fr
Janvier 2018.

---

1 Aspects méthodologique : préparation des données

Le nombre de données relevant quasiment du milliard d’observations, un ordinateur classique se révèle pour l’heure un peu capricieux lors de l’import. La stratégie adoptée consiste à segmenter les observations par morceaux.

1.1 Découpage en tronçons

1.1.1 Une segmentation par tronçons plus petits

Les données dont nous disposons sont partagées en trois tronçons. Malgré ce découpage, le nombre de lignes reste trop conséquent pour l’ordinateur. Nous décidons d’effectuer un découpage plus important.

Pour importer un tronçons, nous utilisons la fonction fread. Nous lisons les observations par petits morceaux de \(20\times10^6\) lignes. Dans un premier temps, il nous est nécessaire de connaître le nombre total d’observations du tronçons sur lequel nous travaillons. L’instruction suivante nous permet de révéler cette information :

```r library(dplyr) library(stringr) library(ggplot2) library(readr) library(tidyr) library(stringi) library(data.table) library(pbapply)

nrow(fread(fichier, select = 1L)) ```

Le nom des colonnes s’obtient facilement :

col_label <- fread(fichier, nrows = 0)
col_label <- colnames(col_label)

Le premier tronçon contient \(412~406~274\) lignes. Nous pouvons créer une variable indiquant les valeurs de débuts et fins de lignes à importer, de manière à lire uniquement \(20\times10^6\) lignes au maximum à chaque importation.

seq_nbs <- seq(1, 412406274, by = 20*10^6)
seq_nbs <- c(seq_nbs, 412406274)
seq_nbs[1] <- 2 # La première ligne correspond aux noms de colonnes
a_parcourir <- data.frame(debut = seq_nbs[-length(seq_nbs)], fin = seq_nbs[-1]-1)

L’importation de chaque tronçon se fait à l’aide de la fonction fread, comme suit (ici, pour le premier tronçon) :

numero_chunk <- 1 # Tronçon 1
val_deb <- a_parcourir$debut[numero_chunk]
val_fin <- a_parcourir$fin[numero_chunk]
val_skip <- val_deb-1
val_nrows <- val_fin-val_deb+1
chunk <- fread(fichier,
encoding = "Latin-1",
header = FALSE, sep = ";",
fill=TRUE, col.names=col_label, quote='',
stringsAsFactors = FALSE,
colClasses=c(rep("character",19)),
na.strings=c("NA","NaN"," ","","0"),
skip = val_skip,
nrows = val_nrows)

Pour ce qui nous intéresse ici, nous pouvons nous permettre de séparer les observations les unes des autres, à condition de conserver celles provenant d’un même arbre d’utilisateur. Nous regroupons chaque observations en fonction du premier caractère alpha numérique des utilistaurs.

# Creation de chunks en fonction de la premiere lettre du nom du prorietaire
chunk[, prem_lettre:=str_to_lower(str_sub(sourcename, 1, 1))]

# Les valeurs différentes pour le nom d'utilisateur
motifs <- unique(chunk$prem_lettre) %>% sort()

# On conserve uniquement les lettres, et on ajoute une valeur pour "non-lettre" (pour les chiffres par exemple)
motifs <- c(motifs[str_detect(motifs, "[:letter:]")], "[^[:letter:]]")

# Sauvegarde dans des tronçons en fonction de la première lettre du nom d'utilisateur
pb2 <- txtProgressBar(min = 1, max = length(motifs), style = 3, char = "@")
for(ii in 1:length(motifs)){
  prem <- motifs[ii]
  if(ii == length(motifs)) prem <- "0"
  chunk_lettre <- chunk[str_detect(prem_lettre, motifs[ii])]
  chunk_lettre[,prem_lettre:=NULL]
  save(chunk_lettre, file = str_c("../data/Geneanet/raw/chunks/chunk_", sprintf("%02s", no_chunks_geneanet), "_", sprintf("%02s", numero_chunk), "_", prem, ".rda"))
  setTxtProgressBar(pb2, ii)
}

Il suffit alors de créer une petite fonction pour faire tourner sur les trois fichiers fournis par Geneanet, et de sauvegarder le résultat dans un fichier de données.

1.1.2 Regroupement des petits tronçons par noms d’utilisateurs

Comme les données initiales sont réparties dans trois fichiers, les petits tronçons obtenus dans l’étape précédente doivent être regroupés (si jamais des données de l’utilisateur toto se trouvent dans au moins deux fichiers différents).

# Chemin vers les tronçons
path_to_files <- "../data/Geneanet/raw/chunks/"
# Liste des fichiers
fichiers <- list.files(path_to_files, pattern = "\\.rda", full.names = TRUE)

motifs <- list.files(path_to_files, pattern = "\\.rda")
motifs <- str_replace(motifs, "chunk_0(1|2|3)_[[:digit:]]{2}_", "") %>% 
  str_replace("\\.rda", "")
motifs <- unique(motifs) %>% sort()
motifs <- c(motifs[str_detect(motifs, "[:letter:]")], "[^[:letter:]]")

# Créé le dossier d'enregistrements des tronçons, si besoin
if(!dir.exists("../data/Geneanet/raw/chunks_letter/")) dir.create("../data/Geneanet/raw/chunks_letter/", recursive = TRUE)

# motif <- "z"
save_chunks_lettre <- function(motif){
  fichiers_charger <- fichiers[str_detect(fichiers, str_c(motif, ".rda"))]
  # Charger les 
  chunks <- pblapply(fichiers_charger, function(x) {load(x) ; unique(chunk_lettre)})
  chunk_lettre <- rbindlist(chunks)
  chunk_lettre <- unique(chunk_lettre)
  if(motif == "[^[:letter:]]") motif <- "0"
  save(chunk_lettre, file = str_c("../data/Geneanet/raw/chunks_letter/chunk_", motif, ".rda"))
}

# Regroupe les tronçons et sauvegarde le résultat
pblapply(motifs, save_chunks_lettre)

1.1.3 Découpage en 5 fichiers pour chaque lettre d’utilisateurs

Afin d’accélérer les étapes suivantes, les tronçons d’utilisateurs sont découpés en 5 parties contenant plus ou moins le même nombre d’utilisateurs différents.

# Les tronçons
N <- list.files("../data/Geneanet/raw/chunks_letter/", pattern = "*.rda", full.names = TRUE)

# Création du dossier d'enregistrement des petits tronçons
if(!dir.exists("../data/Geneanet/raw/chunks_letter_2/")) dir.create("../data/Geneanet/raw/chunks_letter_2/", recursive = TRUE)

# Fonction pour découper des vecteurs en parties à peu près égales
chunk2 <- function(x,n) split(x, cut(seq_along(x), n, labels = FALSE)) 

#' decouper_chunks
#' @i: indice de position du fichier dans N
decouper_chunks <- function(i){
  # Obtenir la lettre des noms d'utilisateurs du tronçon courant
  lettre <- str_replace(N[i], "../data/Geneanet/raw/chunks_letter//chunk_", "")
  lettre <- str_replace(lettre, "\\.rda", "")
  # Chargement du tronçon
  load(N[i], envir = globalenv())
  # Les noms d'utilisateurs
  sourcenames <- unique(chunk_lettre$sourcename)
  # Partage équitable des noms d'utilisateurs (!= des observations)
  chunk_sourcenames <- chunk2(sourcenames, 5)
  
  #' sauvegarder_chunk
  #' @j: indice de chunk_sourcenames à extraire pour traiter
  #'     les utilisateurs qu'il contient
  sauvegarder_chunk <- function(j){
    sourcenames_chunk <- chunk_sourcenames[[j]]
    chunk_lettre_tmp <- chunk_lettre[sourcename %in% sourcenames_chunk]
    save(chunk_lettre_tmp, file = str_c("../data/Geneanet/raw/chunks_letter_2//chunk_lettre_tmp_",
                                        lettre, "_", str_pad(j, width = 2, pad = "0"), ".rda"))
  }# Fin de sauvegarder_chunk()
  
  # Création de 5 petits tronçons pour le fichier courant
  pblapply(1:5, sauvegarder_chunk)
  
}# Fin de decouper_chunks()

# Découper en 5 petits tronçons les fichiers
pblapply(1:length(N), decouper_chunks)

1.2 Regroupement par département et première simplification

Pour déployer les étapes de traitement des données suivantes sur plusieurs machines en même temps, nous regroupons les individus par département. Dans cette même étape, nous regroupons les informations de chaque individu sur une seule ligne. En effet, jusqu’à présent, une ligne correspond à un événement pour un individu dans l’arbre d’un utilisateur.

Pour commencer, nous consignons dans un tableau les différentes régions françaises présentes dans les données :

# Choix d'un tronçon
load("../data/Geneanet/raw/chunks_letter_2/chunk_lettre_tmp_0_05.rda")
liste_depts <-
  chunk_lettre_tmp %>%
  rename(dept = `sous-region`) %>% 
  select(pays, region, dept) %>%
  filter(pays == "FRA") %>%
  unique() %>%
  arrange(region, dept) %>%
  filter(!is.na(dept)) %>%
  tbl_df()

save(liste_depts, file = "liste_depts.rda")

Nous travaillons par tronçons obtenu dans l’étape précédente. Il est de fait nécessaire de les lister.

N <- list.files("../data/Geneanet/raw/chunks_letter_2/", full.names = TRUE)

Le code qui suit permet de traiter un seul tronçon, dont la position dans N est notée i_fichier. Nous prenons le premier fichier en exemple dans le code qui suit. Un simple boucle sur les indices des éléments de N permet de traiter chaque tronçon.

i_fichier <- 1
fichier <- N[i_fichier]
# La première lettre du nom d'utilisateur
lettre <- str_replace(fichier, "../data/Geneanet/raw/chunks_letter_2//chunk_lettre_tmp_", "") %>% 
  str_replace(., "_..\\.rda", "")
# Le numéro de tronçon
num_chunk <- str_replace(fichier, "../data/Geneanet/raw/chunks_letter_2//chunk_lettre_tmp_", "") %>% 
  str_replace(., "._", "") %>% 
  str_replace(., "\\.rda", "")

1.2.1 Quelques fonctions utiles

1.2.1.1 Lieux des événements renseignés

Trois types d’événements sont présents dans la base : naissance, mariage et décès. Nous créons une fonction pour indiquer les informations relatives à ces événements : géographie et date.

#' endroit_acte
#' 
#' Retourne un tableau de donnees indiquant la geographie de l'enregistrement
#' pour un type d'acte donnee
#' 
#' @x: (data.table) informations sur un individu
#' @acte : (string) type d'acte : "N" pour naissance, "M" pour mariage, ou "D" pour deces
#' x <- nai_per ; acte <- "N"
endroit_acte <- function(x, acte){
  if(nrow(x)>0){
    res <- 
      data.table(lieu = x$lieu,
                 dept = x$dept,
                 region = x$region,
                 pays = x$pays,
                 lat = x$latitude,
                 long = x$longitude,
                 stringsAsFactors = FALSE)
    
    if(acte == "N"){
      date <- x$date_naissance
    }else if(acte == "M"){
      date <- NA
    }else{
      date <- x$date_deces
    }
    
    res$date <- date
  }else{
    res <- 
      data.table(lieu = NA, dept = NA, region = NA, pays = NA, lat = NA, long = NA, date = NA, stringsAsFactors = FALSE)
  }
  
  res <- unique(res)
  names(res) <- str_c(names(res), "_", acte)
  res %>% tbl_df()
}# Fin de endroit_acte()

1.2.1.2 Simplification d’une personne

La fonction simplifier_personne() permet de créer un seul enregistrement par individu, pour un arbre d’utilisateur donné. Ainsi, au lieu d’avoir une ligne par événement, nous n’en conservons plus qu’une. Cela a une incidence sur les différents mariages : nous ne pouvons de fait en conserver qu’un seul.

#' simplifier_personne
#' 
#' Retourne les informations relatives a la naissance, le mariage et le deces
#' d'une personne
#' 
#' @un_id_personne: (string) id unique de la personne
#' @dt_source: (data.table) tableau de donnees dans lequel chercher les
#' informations pour cette personne
#' dt_source <- chunk_user ; un_id_personne <- ids_personnes[1]
simplifier_personne <- function(un_id_personne, dt_source){
  
  enregistrements_personne <- dt_source[id_personne == un_id_personne]
  
  sourcename_per <- enregistrements_personne$sourcename %>% unique()
  nom_per <- enregistrements_personne %>% 
    arrange(desc(str_length(nom))) %>% 
    slice(1) %>% 
    .$nom
  prenoms_per <- enregistrements_personne %>% 
    arrange(desc(str_length(prenoms))) %>% 
    slice(1) %>% 
    .$prenoms
  
  sexe_per <- enregistrements_personne$sexe %>% unique()
  
  id_mere_per <- enregistrements_personne$id_mere %>% unique()
  id_pere_per <- enregistrements_personne$id_pere %>% unique()
  
  # Naissance
  nai_per <- enregistrements_personne[stri_detect_regex(event_type, "N")]
  nai <- endroit_acte(nai_per, "N")
  # S'il y a plusieurs lieux de naissance, choisir celui le plus court
  if(nrow(nai)>1){
    nai <- 
      nai %>% 
      filter(!is.na(lieu_N)) %>% 
      arrange(str_length(lieu_N)) %>% 
      slice(1)
  }
  # Mariage
  mar_per <- enregistrements_personne[stri_detect_regex(event_type, "M")]
  mar <- endroit_acte(mar_per, "M")
  # S'il y a plusieurs dates de mariages : on ne retient que la plus ancienne
  if(nrow(mar)>1){
    mar <- 
      mar %>% 
      arrange(date_M) %>% 
      slice(1)
  }
  # Deces
  dec_per <- enregistrements_personne[stri_detect_regex(event_type, "D")]
  dec <- endroit_acte(dec_per, "D")
  # S'il y a plusieurs lieux de deces, choisir celui le plus court
  if(nrow(dec)>1){
    dec <- 
      dec %>% 
      filter(!is.na(lieu_D)) %>% 
      arrange(str_length(lieu_D)) %>% 
      slice(1)
  }
  
  nai$date_N <- enregistrements_personne$date_naissance %>% unique()
  dec$date_D <- enregistrements_personne$date_deces %>% unique()
  
  data.table(sourcename = sourcename_per, id_personne = un_id_personne,
             nom = nom_per, prenoms = prenoms_per, sexe = sexe_per,
             id_mere = id_mere_per, id_pere = id_pere_per, stringsAsFactors = FALSE) %>% 
    cbind(., nai) %>% 
    cbind(., mar) %>% 
    cbind(., dec) %>% 
    unique()
}# Fin de simplifier_personne()

1.2.1.3 Simplification des individus d’un même arbre

La fonction simplifier_proprietaire() permet de simplifier l’ensemble des individus présents dans l’arbre d’un usager de Geneanet. Elle s’appuie sur la fonction simplifier_personne() précédemment définie.

#' simplifier_proprietaire
#' 
#' Pour un proprietaire, simplifie les informations de tous les individus
#' 
#' @id_user: (string) identifiant d'un proprietaire d'arbre
#' @chunk_lettre: (tbl.df) base de donnees contenant les informations des utilisateurs dans le chunk
#' id_user <- "61p"
simplifier_proprietaire <- function(id_user, chunk_lettre){
  # Se restreinde aux individus de l'arbre du proprietaire
  chunk_user <- chunk_lettre[sourcename == id_user]
  
  # Creation d'un identifiant unique par personne
  chunk_user[, id_personne := str_c(sourcename, ref_locale, sep = "@")]
  
  # Ajout de l'identifiant des parents
  identifiants_mere <- chunk_user[, list(ID_num, id_personne)] %>% 
    rename(ID_num_mere = ID_num, id_mere = id_personne)
  
  identifiants_pere <- chunk_user[, list(ID_num, id_personne)] %>% 
    rename(ID_num_pere = ID_num, id_pere = id_personne)
  
  chunk_user <- identifiants_mere[chunk_user, on = "ID_num_mere"]
  chunk_user <- identifiants_pere[chunk_user, on = "ID_num_pere"]
  
  # Suppression de variables qui vont perturber la simplification
  # chunk_user[, c("ref_locale", "ID_num", "ID_num_pere", "ID_num_mere", "ID_num_conjoint") := NULL]
  chunk_user[, c("ID_num", "ID_num_pere", "ID_num_mere", "ID_num_conjoint") := NULL]
  
  # Pour chaque personne, un enregistrement correspond soit :
  # a une naissance (N), un mariage (M), un deces (D),
  # ou un mixe de ces evenements (e.g., NM pour naissance et mariage)
  # Cette distinction est faite parce que les lieux de N, M ou D peuvent varier
  # On va prendre uniquement la date du premier mariage en compte ici.
  ids_personnes <- chunk_user[!is.na(id_personne)]$id_personne %>% unique()
  
  lapply(ids_personnes, simplifier_personne, dt_source = chunk_user) %>% 
    rbindlist
  
}# simplifier_proprietaire

1.2.2 Simplification pour un tronçon

Nous chargeons en mémoire un tronçon précédemment créé, et nous travaillons sur celui-ci. Nous donnons ici les étapes de traitement pour un tronçon ; il est aisé de créer par la suite une fonction et de l’appliquer à chaque tronçon.

# Chargement du chunk
fichier <- N[1]
load(fichier, envir = .GlobalEnv)
chunk_lettre <- chunk_lettre_tmp
rm(chunk_lettre_tmp)

Nous listons les différents usagers (ou propriétaires d’arbre), et procédons à quelques opérations de nommage de variables, afin de faciliter l’écriture du code.

ids_sourcename <- chunk_lettre$sourcename %>% unique()
nbs_obs_ids_sourcename <- length(ids_sourcename)

chunk_lettre <- chunk_lettre %>% rename(dept = `sous-region`)
chunk_lettre <- chunk_lettre %>% rename(ID_num = increment,
                                        ID_num_mere = numero_mere,
                                        ID_num_pere = numero_pere, 
                                        ID_num_conjoint = numero_conjoint,
                                        prenoms = prenom)

Nous avons pris le parti de suivre les descendants d’individus nés entre 1800 et 1804 en France. La base est filtrée en conséquence.

gen_0 <- chunk_lettre
gen_0[, annee_nai := stri_sub(date_naissance, 1, 4)]
gen_0 <- gen_0[annee_nai %in% seq(1800, 1804)]
gen_0 <- gen_0[stri_detect_regex(event_type, "N")]
# Se restreindre uniquement aux abres dont un des individus est né entre 1800 et 1804
chunk_lettre <- chunk_lettre[sourcename %in% unique(gen_0$sourcename)]

1.2.2.1 Pour un tronçon, focus sur un département

Afin d’alléger les opérations à faire subir aux machines, nous travaillons par départements. Nous proposons ici la méthode pour gérer le cas d’un département. À nouveau, une fonction englobant le code qui suit permet de traiter l’ensemble des départements ultérieurement.

i_departement <- "F01"
departement <- liste_depts$dept[i_departement]

# Création du dossier de sauvegarde des résultats
if(!dir.exists(str_c("../data/individuals/migration/", departement, "/"))) dir.create(str_c("../data/individuals/migration/", departement, "/"), recursive = TRUE)

Nous filtrons la base pour se restreindre aux individus de première génération nés dans le département :

gen_0 <- gen_0[dept %in% departement]
# Noms d'utilisateurs de la sous-partie de la base de données
ids_sourcename <- gen_0$sourcename %>% unique()

Il s’agit à présent de se restreindre aux parents et descendants des individus ainsi sélectionnés.

# Initialisation

# La generation courante
gen_n <- gen_0[, list(sourcename, ID_num, ID_num_mere, ID_num_pere)]
nb_obs_gen_0 <- nrow(gen_n)
conserver <- gen_n[, list(sourcename, ID_num)]

# Les generations restantes (on retire la generation courante)
chunk_lettre <- chunk_lettre[sourcename %in% unique(gen_0$sourcename)]
gen_restants <- chunk_lettre[,list(sourcename, ID_num, ID_num_mere, ID_num_pere)]
gen_restants <- fsetdiff(gen_restants, gen_n, all = FALSE)

# Les parents
parents_gen_0 <-
  gen_n %>% 
  select(sourcename, ID_num_mere) %>% 
  rename(ID_num = ID_num_mere) %>% 
  bind_rows(
    gen_n %>% 
      select(sourcename, ID_num_pere) %>% 
      rename(ID_num = ID_num_pere)
  ) %>% 
  unique()

parents_gen_0_complet <- gen_restants[parents_gen_0, on = c("sourcename", "ID_num")]

# On retire les individus trouvés
# Bémol : cette étape retire les individus nés d'une relation incestueuse
gen_restants <- fsetdiff(gen_restants, parents_gen_0_complet, all = FALSE)


# Boucle
# Parmi ces personnes, lesquelles ont pour parent quelqu'un de la generation precedente
compteur <- 0
while(nrow(gen_n) > 0){
  compteur <- compteur+1
  if(compteur>=15) stop("Trop de tours")
  
  parents_m <- data.table(gen_n %>% select(sourcename, ID_num) %>% rename(ID_num_mere = ID_num))
  parents_p <- data.table(gen_n %>% select(sourcename, ID_num) %>% rename(ID_num_pere = ID_num))
  enfants <- data.table(gen_restants %>% select(sourcename, ID_num_mere, ID_num_pere, ID_num))
  
  
  gen_n <- enfants[parents_m, on = c("sourcename", "ID_num_mere")][!is.na(ID_num)] %>% 
    rbind(
      enfants[parents_p, on = c("sourcename", "ID_num_pere")][!is.na(ID_num)]
    ) %>% 
    unique()
  
  conserver <- rbind(conserver, gen_n %>% select(sourcename, ID_num)) %>% unique()
  # Les generations restantes
  gen_restants <- gen_restants[!gen_n, on = c("sourcename", "ID_num")]
}# Fin du while

# Les individus à conserver
conserver <- rbind(parents_gen_0, conserver)
# La base concernant ces individus
chunk_partiel <- chunk_lettre[conserver, on = c("sourcename", "ID_num")]

Reste alors à procéder à la simplification des lignes pour tous ces individus. Le code s’attache à le faire par utilisateur de Geneanet.

res <- pblapply(ids_sourcename, simplifier_proprietaire, chunk_lettre = chunk_partiel)

Comme cette dernière opération informatique prend beaucoup de temps, il est possible de proposer une version parallélisée :

library(parallel)
# Nombre de clusters
ncl <- detectCores()-1
(cl <- makeCluster(ncl))

invisible(clusterEvalQ(cl, library(dplyr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(tidyr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringi, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(data.table, warn.conflicts=FALSE, quietly=TRUE)))

# Evnoi des fonctions/données aux clusters
clusterExport(cl, c("simplifier_personne", "endroit_acte"))

res <- pblapply(ids_sourcename, simplifier_proprietaire, cl = cl, chunk_lettre = chunk_partiel)
stopCluster(cl)

Reste à rassembler chaque résultat dans un seul tableau :

res <- 
  res %>% rbindlist

Et enfin de sauvegarder le résultat :

save(res, file = str_c("../data/individuals/migration/", departement,"/chunk_", lettre, "_", num_chunk, ".rda"))

1.3 Rassemblement des arbres et élimination de doublons

Comme chaque usager du site Geneanet peut construire son propre arbre, il existe de nombreux doublons. Nous adoptons une stratégie en plusieurs étapes pour tenter de joindre les arbres entre-eux, tout en évitant les doublons. Du fait de l’important des données manquantes, cette tâche s’avère délicate. Il persiste sûrement au final des doublons, qui ne seront, pour la plupart, pas utilisés dans l’analyse. En effet, s’ils n’ont pas été repérés en tant que doublons, c’est qu’ils sont porteur de très peu d’informations, et seront donc écartés par la suite.

Cette étape dans le processus de nettoyage des données permet non seulement de relier les arbres d’utilisateurs entre eux, mais également de compléter certaines informations qui pourraient manquer dans l’arbre d’un usager mais être présentes dans celui d’un autre.

1.3.1 Quelques fonctions utiles

Nous présentons ici quelques fonctions permettant d’effectuer les regroupements.

1.3.1.1 Valeur la plus probable

La fonction most_probable_value() permet de trouver la valeur la plus probable parmi les différentes proposées. Elle prend en entrée le nom de la variable à laquelle on s’intéresse, le tableau de données contenant les valeurs, et une variable indiquant si des poids dans les observations sont fournis. La valeur finale sera celle dont la fréquence (pondérée le cas échéant) est la plus élevée parmi les propositions.

Cette fonction s’utilise sur des données concernant des individus étant identifiés comme désignant la même personne.

#' most_probable_value
#' Find the most probable value for a variable
#' using the values found within all candidates
#' (excluding NAs)
#' @variable_name: (string) name of the variable
#' @df: (data.table)
#' @weights: (logic) should the simplification consider the weights? (default: FALSE)
#' variable_name <- "date_N"
most_probable_value <- function(df, variable_name, weights = FALSE){
  valeurs <- df[[variable_name]]
  if(!all(is.na(valeurs))){
    if(weights){
      poids <- df[["weight"]]
      res <- data.frame(valeurs, poids, stringsAsFactors = FALSE) %>% 
        group_by(valeurs) %>% 
        summarise(poids = sum(poids)) %>% 
        ungroup()
      
      # Si la variable est une date
      # on va privilegier les informatins relatives aux dates completes
      if(variable_name %in% c("date_N", "date_M", "date_D")){
        tmp <- res %>% 
          mutate(mois = str_sub(valeurs, 5, 6),
                 jour = str_sub(valeurs, 7, 8)) %>% 
          filter(!(mois == "00" | jour == "00"))
        # S'il reste des informations, on se base sur elles
        if(nrow(tmp)>0) res <- tmp
      }
      
      res <- 
        res %>%
        arrange(desc(poids)) %>% 
        slice(1) %>% 
        magrittr::extract2("valeurs")
    }else{# Si pas de poids
      table_freq <- sort(table(valeurs))
      res <- names(table_freq)[1]
    }
  }else{
    res <- NA
  }
  res
}# End of most_probable_value()

1.3.1.2 Valeur la plus probable

La fonction simplifier_personne(), à partir d’individus identifiés comme désignant la même personne, simplifie les valeurs de chaque variable pour n’obtenir qu’une seule observation en sortie. Les valeurs de chaque variable sont obtenues en faisant appel à la fonction most_probable_value() précédemment définie.

# df_corresp <- corresp ; un_groupe_num <- 2 ; weights <- TRUE
# rm(df_corresp, un_groupe_num, groupe_cour, individus_cour, nouvelles_valeurs, weight)
#' @weights: (logic) should the simplification consider the weights? (default: TRUE)
simplifier_personne <- function(df_corresp, un_groupe_num, weights = TRUE){
  groupe_cour <- df_corresp[id_personne_num_groupe %in% un_groupe_num]
  
  individus_cour <- individus_simple[id_personne_num %in% groupe_cour$id_personne_num]
  
  weight <- sum(individus_cour$weight)
  
  if(nrow(individus_cour) == 1){
    nouvelles_valeurs <- individus_cour[, mget(c(variables_names))]
  }else{
    # Les valeurs probables pour les variables d'interet
    nouvelles_valeurs <- 
      variables_names %>% 
      sapply(., most_probable_value, df = individus_cour, weights = weights) %>% 
      t() %>% 
      data.table(stringsAsFactors = FALSE)
  }
  
  
  cbind(id_personne_num = un_groupe_num, nouvelles_valeurs, weight = weight)
}# End of simplifier_personne()

1.3.1.3 Une fonction de minimum

Une fonction permettant de gérer les valeurs manquantes lors de l’utilisation de la fonction min().

my_min <- function(x) ifelse( !all(is.na(x)), min(x, na.rm=T), NA)

1.3.2 Chargement des données

Nous allons retenir uniquement certaines variables pour les individus. Les variables textuelles de lieu sont évincées ici.

# Les vairables qui nous interessent
variables_names <- c("nom", "prenoms",
"sexe",
"lieu_N", "dept_N", "region_N", "pays_N", "lat_N", "long_N", "date_N",
"lieu_M", "dept_M", "region_M", "pays_M", "lat_M", "long_M", "date_M",
"lieu_D", "dept_D", "region_D", "pays_D", "lat_D", "long_D", "date_D",
"prenom")

Il s’agit ensuite de charger les données en mémoire dans la session. Nous repartons des données par département, en se concentrant uniquement sur les arbres dans lesquels se trouvent des individus nés entre 1800 et 1804. Nous ne montrons le traitement des donnée que pour un département. Il est aisé d’enrober le code dans une fonction pour ensuite le déployer sur l’ensemble des départements.

departement <- liste_depts$dept[i_departement]
# L'ensemble des données pour ce département
N <- list.files(str_c("../data/individuals/migration/", departement, "/"), pattern = "*.rda", full.names = TRUE)
# Chargement
individus <- 
  pblapply(N, function(x){load(x) ; res}) %>% 
  rbindlist() %>% 
  unique()

Les observations pour lesquelles les coordonnées de lieux sont 0.00 sont manquantes, donc à interpréter comme non disponibles dans R.

individus <- 
  individus %>% 
  mutate(long_N = ifelse(long_N == "0.00000", yes = NA, no = long_N),
         long_M = ifelse(long_M == "0.00000", yes = NA, no = long_M),
         long_D = ifelse(long_D == "0.00000", yes = NA, no = long_D)) %>% 
  mutate(lat_N = ifelse(lat_N == "0.00000", yes = NA, no = lat_N),
         lat_M = ifelse(lat_M == "0.00000", yes = NA, no = lat_M),
         lat_D = ifelse(lat_D == "0.00000", yes = NA, no = lat_D))

1.3.3 Préparation avant la première simplification

1.3.3.1 Création d’une base d’individus

Nous allons nous appuyer sur des identifiants numériques pour les individus. Nous créons une base individus recensant tous les individus et leur attribuant un identifiant.

# Les individus sont identifies de maniere unique par la colonne
# id_personne. Il s'agit d'une chaine de caracteres qui peut etre assez longue
# On va utiliser un identifiant numerique, ce qui sera plus pratique.
individus <- 
  individus %>% 
  tbl_df() %>% 
  mutate(id_personne_num = row_number()) %>% 
  data.table()

On va creer une variable contenant les prénoms de manière simplifiée, idem pour les noms. Nous retirons les espaces et caractères spéciaux.

individus <- 
  individus %>% 
  tbl_df() %>% 
  mutate(nom = stri_replace_all_regex(nom, "\\((.*?)\\)", "") %>% 
           stri_trans_general(., "Latin-ASCII") %>% 
           stri_replace_all_regex(., "'|\\(|\\)|[[:blank:]]|\\?|_|\\*", "") %>% 
           stri_trim(.) %>% 
           str_to_lower()) %>% 
  mutate(prenom = 
           str_replace_all(prenoms, "\\((.*?)\\)", "") %>% 
           stri_replace_all_regex(., "'|\\(|\\)|\\?|_|\\*|&gt|\\^|\\+|[[:digit:]]|\\.|<|>", "") %>% 
           str_trim() %>% 
           gsub(" .*$", "", .) %>%
           stri_trans_general(., "Latin-ASCII") %>% 
           stri_replace_all_regex(., "-", "") %>% 
           stri_trim(.) %>% 
           str_to_lower()) %>% 
  data.table()

# Les individus uniquement
individus_simple <- copy(individus)

Ensuite, nous ajoutons les informations relatives aux parents de chaque individu, afin de pouvoir identifier les individus en fonction de leurs parents. Si Jeanne Michu dans l’arbre d’un utilisateur a pour parents Jean Michu et Irène Dupond, et que dans un autre arbre, on trouve une personne née au même moment avec des parents portant également les mêmes noms, cela nous permet de rapprocher non seulement Jeanne, mais aussi Jean et Irène dans les deux arbres pour les fusionner.

On prépare de fait des tableaux pour les informations des parents.

# On va creer une base large avec pour chaque ligne,
# les informations sur les parents egalement.
# De fait, chaque ligne representera un triplet (enfant ; mere ; pere)
ind_meres <- match(individus$id_mere, individus$id_personne)
ind_peres <- match(individus$id_pere, individus$id_personne)

# Les meres
individus_meres <- 
  individus[ind_meres,] %>% 
  magrittr::set_colnames(str_c(colnames(individus), "_mere")) %>% 
  select(-sourcename_mere, -id_mere_mere, -id_pere_mere, -id_personne_mere)

# Les peres
individus_peres <- 
  individus[ind_peres,] %>% 
  magrittr::set_colnames(str_c(colnames(individus), "_pere")) %>% 
  select(-sourcename_pere, -id_mere_pere, -id_pere_pere, -id_personne_pere)

Puis on rajoute ces informations aux lignes des individus.

# On ajoute les informations relatives aux parents
individus <- 
  individus %>% 
  cbind(individus_meres) %>% 
  cbind(individus_peres) %>% 
  select(-id_mere, -id_pere) %>% 
  tbl_df()
rm(individus_meres, individus_peres)

Nous supprimons les observations pour lesquelles il y a une anomalie concernant la date de naissance les concernant eux ou leurs parents. Si un des deux parents a son enfant avant 13 ans, nous considérons qu’il s’agit d’une anomalie.

age_min_nenf <- 13

individus <-
  individus %>% 
  mutate(date_N_annee = str_sub(date_N, 1, 4) %>% as.numeric(),
         date_N_annee_mere = str_sub(date_N_mere, 1, 4) %>% as.numeric(),
         date_N_annee_pere = str_sub(date_N_pere, 1, 4) %>% as.numeric()) %>% 
  mutate(anomalie_mere = date_N_annee <= (date_N_annee_mere+age_min_nenf),
         anomalie_mere = ifelse(is.na(anomalie_mere), yes = FALSE, no = anomalie_mere)) %>% 
  mutate(anomalie_pere = date_N_annee <= (date_N_annee_pere+age_min_nenf),
         anomalie_pere = ifelse(is.na(anomalie_pere), yes = FALSE, no = anomalie_pere)) %>% 
  mutate(anomalie = anomalie_mere + anomalie_pere) %>% 
  filter(anomalie == 0) %>% 
  select(-anomalie_mere, -anomalie_pere, -anomalie) %>% 
  data.table()

# On repercute sur la table des individus
individus_simple <- individus_simple[id_personne_num %in% individus$id_personne_num]

individus_simple <- individus_simple[, mget(c("id_personne_num", variables_names))]

Nous rajoutons une colonne indiquant le poids de chaque observation, en prévision des simplifications à venir.

individus$weight <- 1
individus_simple$weight <- 1

1.3.3.2 Création d’un registre

Nous créons un registre indiquant les liens de parentés entre les individus. Ce registre indique pour chaque personne identifiée par son identifiant numérique, les identifiants numériques respectifs de ses parents. Notons que ce registre sera mis à jour par la suite, et se révèlera utile afin d’effectuer les rapprochements entre individus.

registre <- 
  individus %>% 
  select(id_personne_num, id_personne_num_mere, id_personne_num_pere) %>% 
  tbl_df()

1.3.4 Première simplification

Pour commencer, nous allons proceder a un ecremage tres grossier de la base, en supprimant les doublons evidents, c’est-à-dire ayant :

• le même nom de famille ;
• le même prénom ;
• le même sexe ;
• la même date de naissance ;
• le même nom de la mère ;
• le même nom de père.

1.3.4.1 Repérage d’individus identiques

Avant de procéder à la fusion des individus dédoublonés en entités uniques, nous devons les repérer. Le registre devient utile ici. Dès que deux individus sont repérées comme désignant une seule et même personne, l’identifiant numérique le plus petit d’entre les deux devient le nouvel identifiant des deux. Il en va de même pour les parents.

corresp_meres <- 
  individus[, list(id_personne,id_personne_num, nom, prenom, sexe, date_N, nom_mere, prenom_mere, sexe_mere, date_N_mere)] %>% 
  group_by(nom, prenom, sexe, date_N, nom_mere, prenom_mere, sexe_mere, date_N_mere) %>% 
  mutate(id_personne_num_groupe_mere = min(id_personne_num)) %>% 
  ungroup() %>% 
  select(id_personne, id_personne_num, id_personne_num_groupe_mere) %>% 
  data.table()

corresp_peres <- 
  individus[, list(id_personne, id_personne_num, nom, prenom, sexe, date_N, nom_pere, prenom_pere, sexe_pere, date_N_pere)] %>% 
  group_by(nom, prenom, sexe, date_N, nom_pere, prenom_pere, sexe_pere, date_N_pere) %>% 
  mutate(id_personne_num_groupe_pere = min(id_personne_num)) %>% 
  ungroup() %>% 
  select(id_personne, id_personne_num, id_personne_num_groupe_pere) %>% 
  data.table()

corresp <- 
  corresp_meres %>% 
  left_join(corresp_peres, by = c("id_personne_num", "id_personne")) %>% 
  data.table()

corresp <- 
  corresp %>%
  # rowwise() %>% 
  mutate(id_personne_num_groupe = pmin(id_personne_num_groupe_mere, id_personne_num_groupe_pere, na.rm=T)) %>%
  select(-id_personne_num_groupe_mere, -id_personne_num_groupe_pere) %>% 
  data.table()

corresp <- unique(corresp)

Cette première simplification permet peut-être (et c’est le cas) de regrouper davantage d’individus entre eux, en fonction du registre. Nous mettons alors une procédure itérative visant à repérer les fusions qu’il est possible de faire. Dès qu’il n’y a plus de regroupement possibles, la boucle s’achève.

continuer <- TRUE
nb_iterations <- 1

while(continuer){
  registre_new <- 
    registre %>% 
    tbl_df() %>%
    left_join(corresp %>% select(-id_personne) %>% rename(id_personne_num_new = id_personne_num_groupe), by = "id_personne_num") %>% 
    left_join(
      corresp %>% select(-id_personne) %>% rename(id_personne_num_mere = id_personne_num,
                                                  id_personne_num_mere_new = id_personne_num_groupe),
      by = "id_personne_num_mere"
    ) %>%
    left_join(
      corresp %>% select(-id_personne) %>% rename(id_personne_num_pere = id_personne_num,
                                                  id_personne_num_pere_new = id_personne_num_groupe),
      by = "id_personne_num_pere"
    )
  
  
  registre_new <- 
    registre_new %>% 
    mutate(same_id = id_personne_num == id_personne_num_new,
           same_id_mother = id_personne_num_mere == id_personne_num_mere_new,
           same_id_father = id_personne_num_pere == id_personne_num_pere_new) %>% 
    mutate(same_id = ifelse(is.na(same_id), yes = TRUE, no = same_id),
           same_id_mother = ifelse(is.na(same_id_mother), yes = TRUE, no = same_id_mother),
           same_id_father = ifelse(is.na(same_id_father), yes = TRUE, no = same_id_father))
  
  continuer <- any(!registre_new$same_id | !registre_new$same_id_mother | !registre_new$same_id_father)
  
  if(!continuer) next
  
  # On a pu retrouver d'autres associations !
  # Celles des parents d'un enfant.
  
  # Pour les mères
  corresp_meres <-
    registre_new %>% 
    select(id_personne_num_new, id_personne_num_mere) %>% 
    unique() %>% 
    arrange(id_personne_num_new, id_personne_num_mere) %>% 
    filter(!is.na(id_personne_num_mere))
  
  corresp_meres <- data.table(corresp_meres)
  corresp_meres[, id_personne_num_mere_2 := my_min(id_personne_num_mere), by = id_personne_num_new]
  corresp_meres <- 
    corresp_meres %>% select(-id_personne_num_new) %>% 
    rename(id_personne_num = id_personne_num_mere, id_personne_num_groupe = id_personne_num_mere_2) %>% 
    unique()
  
  corresp_meres[, id_personne_num_groupe := my_min(id_personne_num_groupe), by = id_personne_num]
  
  # Idem pour les pères
  corresp_peres <-
    registre_new %>% 
    select(id_personne_num_new, id_personne_num_pere) %>% 
    unique() %>% 
    arrange(id_personne_num_new, id_personne_num_pere) %>% 
    filter(!is.na(id_personne_num_pere)) 
  corresp_peres <- data.table(corresp_peres)
  corresp_peres[, id_personne_num_pere_2 := my_min(id_personne_num_pere), by = id_personne_num_new]
  corresp_peres <- 
    corresp_peres %>% select(-id_personne_num_new) %>% 
    rename(id_personne_num = id_personne_num_pere, id_personne_num_groupe = id_personne_num_pere_2) %>% 
    unique()
  
  corresp_peres[, id_personne_num_groupe := my_min(id_personne_num_groupe), by = id_personne_num]
  
  corresp <-
    corresp %>% 
    tbl_df() %>% 
    left_join(
      corresp_meres %>% 
        rename(id_personne_num_groupe_new = id_personne_num_groupe),
      by = "id_personne_num"
    ) %>% 
    tbl_df() %>% 
    mutate(id_personne_num_groupe = pmin(id_personne_num_groupe, id_personne_num_groupe_new, na.rm=T)) %>% 
    select(-id_personne_num_groupe_new) %>% 
    left_join(
      corresp_peres %>% 
        rename(id_personne_num_groupe_new = id_personne_num_groupe),
      by = "id_personne_num"
    ) %>% 
    mutate(id_personne_num_groupe = pmin(id_personne_num_groupe, id_personne_num_groupe_new, na.rm=T)) %>% 
    select(-id_personne_num_groupe_new) %>% 
    unique()
  
  registre_new <- 
    registre_new %>% 
    mutate(id_personne_num = id_personne_num_new,
           id_personne_num_mere = id_personne_num_mere_new,
           id_personne_num_pere = id_personne_num_pere_new) %>% 
    select(-id_personne_num_new, -id_personne_num_mere_new, -id_personne_num_pere_new,
           -same_id, -same_id_mother, -same_id_father)
  
  registre_new <- data.table(registre_new)
  registre_new <- registre_new[, list(
    id_personne_num_mere = min(id_personne_num_mere, na.rm=T),
    id_personne_num_pere = min(id_personne_num_pere, na.rm=T)), by = id_personne_num] %>% 
    unique()
  
  registre_new$id_personne_num_mere[which(registre_new$id_personne_num_mere == Inf)] <- NA
  registre_new$id_personne_num_pere[which(registre_new$id_personne_num_pere == Inf)] <- NA
  
  # toto <- data.frame(id = c(1,1,2,3,3,4,3), val = c(1,0,2,NA, NA, 2,NA)) %>% data.table()
  # toto[, list(test = pmin(val)), by = id]
  # 
  # 
  # registre_new <- 
  #   registre_new %>% 
  #   tbl_df() %>% 
  #   group_by(id_personne_num) %>% 
  #   summarise(id_personne_num_mere = my_min(id_personne_num_mere),
  #             id_personne_num_pere = my_min(id_personne_num_pere)) %>% 
  #   ungroup()
  
  # registre_new <- unique(registre_new)
  
  registre <- copy(registre_new)
  nb_iterations <- nb_iterations+1
}# Fin de la boucle while

corresp <- data.table(corresp)
registre <- data.table(registre)

1.3.4.2 Fusion des informations des individus

Comme le repérage des dédoublonsages est effectué, il reste à fusionner les informations des doublons. Nous utilisons la fonction simplifier_personne() présentée précédemment pour le faire. Si trois Madame Michu sont repérées comme désignant la même personne, et que deux d’entre-elles ont la date de naissance suivante 1800-01-01 tandis que la troisième renseigne un jour différent, alors, du fait d’une fréquence plus élevée pour la valeur 1800-01-01, nous retiendrons cette dernière comme la date de naissance de Madame Michu.

Nous procédons à la simplifications en fonction des groupes de personnes identiques que l’on a identifiés.

corresp_nb <- 
  corresp %>% 
  tbl_df() %>% 
  group_by(id_personne_num_groupe) %>% 
  tally()

Il n’est pas nécessaire de s’attarder sur les singletons, qui ne nécessitent aucun traitement et qui ralentiraient l’exécution du code. Nous les mettons de côté ; ils seront réintégrés une fois la simplification effectuée.

groupes <- corresp_nb %>% filter(n > 1) %>% magrittr::extract2("id_personne_num_groupe")

Il reste à parcourir chacun de ces groupes afin d’effectuer la simplification. Comme ils sont nombreux, nous parallélisons l’exécution du code.

library(parallel)
ncl <- detectCores()-1
(cl <- makeCluster(ncl))

invisible(clusterEvalQ(cl, library(dplyr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(tidyr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringi, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(data.table, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringdist, warn.conflicts=FALSE, quietly=TRUE)))

# Envoi des fonctions/tableaux 
clusterExport(cl, c("most_probable_value", "variables_names", "individus_simple"), envir=environment())

individus_simple_enfants <- pblapply(groupes, simplifier_personne, df_corresp = corresp, cl = cl)
stopCluster(cl)
individus_simple_enfants <- rbindlist(individus_simple_enfants)

Nous rajoutons les singletons laissés de côté.

groupes_un <- corresp_nb %>% filter(n == 1) %>% magrittr::extract2("id_personne_num_groupe")
individus_simple_un <- individus_simple[id_personne_num %in% groupes_un]
individus_simple <- rbind(individus_simple_un, individus_simple_enfants)
rm(corresp_meres, corresp_peres, individus_simple_enfants, registre_new, groupes, groupes_un, individus_simple_un)

Nous disposons donc de trois tables :

1. La table de registres, qui indique pour chaque personne : son identifiant, ainsi que ceux de ses parents (si disponibles) ;
2. La table des individus, avec leur id et leurs renseignements ;
3. La table de correspondance (n’est plus utile, nous allons en construire une autre).

1.3.5 Seconde simplification

Il reste dans la base de données, des doublons moins évidents. Dans l’esprit, nous procédons de la même manière que pour la première simplification, en rajoutant les informations relatives au parent et en créan un registre, puis en simplifiant les doublons. Dans cette seconde simplification, nous prenons en compte les erreurs d’orthographe dans les noms et prénoms.

1.3.5.1 Repérage d’individus identiques

# On va creer une base large avec pour chaque ligne,
# les informations sur les parents egalement.
# De fait, chaque ligne representera un triplet (enfant ; mere ; pere)
individus_2 <- 
  copy(individus_simple)

# On rajoute les identifiants des parents
individus_2$id_personne_num_mere <- registre$id_personne_num_mere[match(individus_simple$id_personne_num, registre$id_personne_num)]
individus_2$id_personne_num_pere <- registre$id_personne_num_pere[match(individus_simple$id_personne_num, registre$id_personne_num)]


ind_meres <- match(individus_2$id_personne_num_mere, individus_2$id_personne_num)
ind_peres <- match(individus_2$id_personne_num_pere, individus_2$id_personne_num)

# Les meres
individus_meres <- 
  individus_2[ind_meres,] %>% 
  magrittr::set_colnames(str_c(colnames(individus_2), "_mere")) %>% 
  select(-id_personne_num_mere, -id_personne_num_mere_mere, -id_personne_num_pere_mere)

# Les peres
individus_peres <- 
  individus_2[ind_peres,] %>% 
  magrittr::set_colnames(str_c(colnames(individus_2), "_pere")) %>% 
  select(-id_personne_num_pere, -id_personne_num_mere_pere, -id_personne_num_pere_pere)


# On ajoute les informations relatives aux parents
individus_2 <- 
  individus_2 %>% 
  cbind(individus_meres) %>% 
  cbind(individus_peres) %>% 
  # select(-id_mere, -id_pere) %>% 
  tbl_df()

individus_2 <- data.table(individus_2)

individus_2$date_N_annee <- individus_2$date_N %>% str_sub(1,4) %>% as.numeric()

Nous repérons les doublons en tentant de prendre en compte les erreurs d’orthographe dans les noms et prénoms. Pour ce faire, nous nous appuyons sur la fonction stringdist() du package du même nom.

library(stringdisy)
#' calculer_distance
#' help function for ecremer()
# variable <- "nom"
calculer_distance <- function(personne_cour, candidats, variable, threshold){
  stringdist(personne_cour[,get(variable)], candidats[,get(variable)], method = "cosine") < threshold
}# End of calculer_distance()

Nous utilisons la fonction ecremer() définie ci-après pour repérer les doublons potentiels. Nous considérons à présent que deux observations avec à peu près les mêmes nom et prénom, à peu près le même nom de mère ou de père sont potentiellement fusionnables.

#' ecremer
#' Retourne un tableau de candidats possibles
#' @candidats: (data.frame) table dans laquelle chercher les candidats
ecremer <- function(candidats, threshold = 0.1, personne_cour, prenoms = FALSE){
  # Ecremer par le nom
  ind_nom <- which(calculer_distance(personne_cour = personne_cour, candidats = candidats, variable = "nom", threshold = threshold))
  candidats <- candidats[ind_nom,]
  
  if(nrow(candidats)>0){
    # Par le prenom
    if(prenoms == TRUE){
      ind_prenom <- which(calculer_distance(personne_cour = personne_cour, candidats = candidats, variable = "prenoms", threshold = threshold)) 
    }else{
      ind_prenom <- which(calculer_distance(personne_cour = personne_cour, candidats = candidats, variable = "prenom", threshold = threshold))
    }
    candidats <- candidats[ind_prenom,]
    
    if(nrow(candidats)>0 & !is.na(personne_cour$nom_mere)){
      ind_mere_nom <- calculer_distance(personne_cour = personne_cour, candidats = candidats, variable = "nom_mere", threshold = threshold)
      ind_mere_nom <- c(which(is.na(ind_mere_nom)), which(ind_mere_nom))
      candidats <- candidats[ind_mere_nom,]
      
    }# Fin if mere
    
    if(nrow(candidats)>0 & !is.na(personne_cour$nom_pere)){
      ind_pere_nom <- calculer_distance(personne_cour = personne_cour, candidats = candidats, variable = "nom_pere", threshold = threshold)
      ind_pere_nom <- c(which(is.na(ind_pere_nom)), which(ind_pere_nom)) %>% sort()
      candidats <- candidats[ind_pere_nom,]
    }# Fin if pere

  }# Fin du if
  candidats
}# End of ecremer()

Enfin, pour pouvoir effectuer le repérage, nous utilisons la fonction trouver_corresp_indice_naissance() qui s’appuie sur les résultats obtenus par la fonction ecremer() pour définir si un rapprochement peut être effectué entre plusieurs individus. Il faut toutefois que les personnes, outre le fait d’avoir un nom similaire, soient nées dans le même département la même année.

#' trouver_corresp_indice_naissance
#' Pour un numero de ligne du tableau individus
#' retourne sous forme de table les correspondances entre individus.
#' Chaque ligne du resultat correspond a l'identifiant d'une personne
#' et au plus petit identifiant d'une correspondance trouvee dans la base.
trouver_corresp_indice_naissance <- function(i){
  personne_cour <- individus_3[i, ]
  
  resul <- data.table(id_personne_num = personne_cour$id_personne_num,
                      id_personne_num_groupe = personne_cour$id_personne_num)
  
  # Les identifiants des parents
  personne_cour_id_num_mere <- personne_cour$id_personne_num_mere
  personne_cour_id_num_pere <- personne_cour$id_personne_num_pere
  
  # Si les identifiants des deux parents sont manquants, on ne va pas plus loin
  if(!(is.na(personne_cour_id_num_mere) & is.na(personne_cour_id_num_pere))){
    # On a un identifiant pour: mere U pere
    
    id_personne_num <- personne_cour$id_personne_num
    
    # Il faut que les candidats soient de la même année
    # Et du meme sexe
    
    # Il faut que les candidats soient nes la meme annee et soient du meme sexe
    ind <- which(individus_3$date_N_annee %in% personne_cour$date_N_annee & 
                   individus_3$sexe %in% personne_cour$sexe)
    
    candidats <- individus_3[ind]
    
    # individus %>% filter(id_personne == "yvono@le garrec|jean|1") %>% View()
    
    candidats <- candidats[!id_personne_num %in% personne_cour$id_personne_num]
    
    candidats <- ecremer(candidats, personne_cour = personne_cour)
    
    if(nrow(candidats)>0){
      id_retenir <- c(id_personne_num, candidats$id_personne_num)
      resul <- data.table(id_personne_num = id_retenir, id_personne_num_groupe = my_min(id_retenir))
    }# Fin du if()
  }# Fin du if sur les parents
  
  # }# Fin du if sur la date de naissance
  
  resul
}# End of trouver_corresp_indice_naissance()

Nous pouvons utiliser les fonctions précédemment introduite pour repérer les individus à fusionner. Afin de ne pas perdre de temps dans l’exécution du code pour des individus étant rapidement identifiés comme ne pouvant pas être fusionnés, nous repérons a priori ces individus et les mettons de côté.

ind_ne_pas_faire <- which(is.na(individus_2$date_N_annee) | 
                            is.na(individus_2$sexe) | 
                            (is.na(individus_2$id_personne_num_mere) & is.na(individus_2$id_personne_num_pere)))

corresp_naissance_ne_pas_faire <- 
  individus_2[ind_ne_pas_faire,] %>% 
  select(id_personne_num) %>% 
  mutate(id_personne_num_groupe = id_personne_num) %>% 
  data.table()

Reste à exécuter la fonction trouver_corresp_indice_naissance() pour chaque individu restant. À nouveau, cette exécution est parallélisée.

ind_faire <- seq(1, nrow(individus_2))
ind_faire <- ind_faire[-ind_ne_pas_faire]

library(parallel)
ncl <- detectCores()-1
(cl <- makeCluster(ncl))

invisible(clusterEvalQ(cl, library(dplyr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(tidyr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringi, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(data.table, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringdist, warn.conflicts=FALSE, quietly=TRUE)))

individus_3 <- individus_2[ind_faire]

clusterExport(cl, c("individus_3"), envir=environment())
clusterExport(cl, c("calculer_distance", "ecremer", "my_min"), envir=environment())
corresp_naissance <- pblapply(1:nrow(individus_3), trouver_corresp_indice_naissance, cl = cl)
stopCluster(cl)

corresp_naissance <- 
  corresp_naissance %>% 
  rbindlist()

corresp_naissance <- unique(corresp_naissance)

corresp_naissance <- 
  corresp_naissance %>% 
  bind_rows(corresp_naissance_ne_pas_faire) %>% 
  unique()

rm(ind_faire, individus_3)

Nous mettons ensuite la table de correspondance à jour.

corresp_naissance <- corresp_naissance[, list(id_personne_num_groupe = min(id_personne_num_groupe, na.rm=T)), by = id_personne_num] %>% unique()

corresp <- 
  corresp %>% 
  left_join(corresp_naissance %>% rename(id_personne_num_groupe_new = id_personne_num_groupe)) %>% 
  mutate(id_personne_num_groupe = pmin(id_personne_num_groupe, id_personne_num_groupe_new, na.rm=T)) %>% 
  select(-id_personne_num_groupe_new) %>% 
  data.table()

Puis, nous regardons dans cette table, de manière itérative comme lors de la première opération de simplification, s’il n’est pas possible d’effectuer d’autres rapprochements en fonction des liens de parentés.

continuer <- TRUE
nb_iterations <- 1

while(continuer){
  registre_new <- 
    registre %>% 
    tbl_df() %>%
    left_join(corresp %>% select(-id_personne) %>% rename(id_personne_num_new = id_personne_num_groupe), by = "id_personne_num") %>% 
    left_join(
      corresp %>% select(-id_personne) %>% rename(id_personne_num_mere = id_personne_num,
                                                  id_personne_num_mere_new = id_personne_num_groupe),
      by = "id_personne_num_mere"
    ) %>%
    left_join(
      corresp %>% select(-id_personne) %>% rename(id_personne_num_pere = id_personne_num,
                                                  id_personne_num_pere_new = id_personne_num_groupe),
      by = "id_personne_num_pere"
    )
  
  
  registre_new <- 
    registre_new %>% 
    mutate(same_id = id_personne_num == id_personne_num_new,
           same_id_mother = id_personne_num_mere == id_personne_num_mere_new,
           same_id_father = id_personne_num_pere == id_personne_num_pere_new) %>% 
    mutate(same_id = ifelse(is.na(same_id), yes = TRUE, no = same_id),
           same_id_mother = ifelse(is.na(same_id_mother), yes = TRUE, no = same_id_mother),
           same_id_father = ifelse(is.na(same_id_father), yes = TRUE, no = same_id_father))
  
  continuer <- any(!registre_new$same_id | !registre_new$same_id_mother | !registre_new$same_id_father)
  
  if(!continuer) next
  
  corresp_meres <-
    registre_new %>% 
    select(id_personne_num_new, id_personne_num_mere) %>% 
    unique() %>% 
    arrange(id_personne_num_new, id_personne_num_mere) %>% 
    filter(!is.na(id_personne_num_mere))
  
  corresp_meres <- data.table(corresp_meres)
  corresp_meres[, id_personne_num_mere_2 := my_min(id_personne_num_mere), by = id_personne_num_new]
  corresp_meres <- 
    corresp_meres %>% select(-id_personne_num_new) %>% 
    rename(id_personne_num = id_personne_num_mere, id_personne_num_groupe = id_personne_num_mere_2) %>% 
    unique()
  
  corresp_meres[, id_personne_num_groupe := my_min(id_personne_num_groupe), by = id_personne_num]
  
  corresp_peres <-
    registre_new %>% 
    select(id_personne_num_new, id_personne_num_pere) %>% 
    unique() %>% 
    arrange(id_personne_num_new, id_personne_num_pere) %>% 
    filter(!is.na(id_personne_num_pere)) 
  
  corresp_peres <- data.table(corresp_peres)
  corresp_peres[, id_personne_num_pere_2 := my_min(id_personne_num_pere), by = id_personne_num_new]
  corresp_peres <- 
    corresp_peres %>% select(-id_personne_num_new) %>% 
    rename(id_personne_num = id_personne_num_pere, id_personne_num_groupe = id_personne_num_pere_2) %>% 
    unique()
  
  corresp_peres[, id_personne_num_groupe := my_min(id_personne_num_groupe), by = id_personne_num]
  
  corresp <-
    corresp %>% 
    tbl_df() %>% 
    left_join(
      corresp_meres %>% 
        rename(id_personne_num_groupe_new = id_personne_num_groupe),
      by = "id_personne_num"
    ) %>% 
    tbl_df() %>% 
    mutate(id_personne_num_groupe = pmin(id_personne_num_groupe, id_personne_num_groupe_new, na.rm=T)) %>% 
    select(-id_personne_num_groupe_new) %>% 
    left_join(
      corresp_peres %>% 
        rename(id_personne_num_groupe_new = id_personne_num_groupe),
      by = "id_personne_num"
    ) %>% 
    mutate(id_personne_num_groupe = pmin(id_personne_num_groupe, id_personne_num_groupe_new, na.rm=T)) %>% 
    select(-id_personne_num_groupe_new) %>% 
    unique()
  
  registre_new <- 
    registre_new %>% 
    mutate(id_personne_num = id_personne_num_new,
           id_personne_num_mere = id_personne_num_mere_new,
           id_personne_num_pere = id_personne_num_pere_new) %>% 
    select(-id_personne_num_new, -id_personne_num_mere_new, -id_personne_num_pere_new,
           -same_id, -same_id_mother, -same_id_father)
  
  registre_new <- data.table(registre_new)
  registre_new <- registre_new[, list(
    id_personne_num_mere = min(id_personne_num_mere, na.rm=T),
    id_personne_num_pere = min(id_personne_num_pere, na.rm=T)), by = id_personne_num] %>% 
    unique()
  
  registre_new$id_personne_num_mere[which(registre_new$id_personne_num_mere == Inf)] <- NA
  registre_new$id_personne_num_pere[which(registre_new$id_personne_num_pere == Inf)] <- NA
  
  registre <- copy(registre_new)
  nb_iterations <- nb_iterations+1
}# Fin de la boucle while

corresp <- data.table(corresp)
registre <- data.table(registre)

1.3.5.2 Fusion des informations des individus

Le répérage des nouveaux doublons est effectué, il reste alors à fusionner les doublons, comme lors de la première fusion.

individus_simple <- individus_2[, mget(colnames(individus_simple))]

corresp_naissance_nb <- 
  corresp %>% 
  tbl_df() %>% 
  group_by(id_personne_num_groupe) %>% 
  tally()


groupes_naissance <- corresp_naissance_nb %>% filter(n > 1) %>% magrittr::extract2("id_personne_num_groupe")

library(parallel)
ncl <- detectCores()-1
(cl <- makeCluster(ncl))

invisible(clusterEvalQ(cl, library(dplyr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(tidyr, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(stringi, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(data.table, warn.conflicts=FALSE, quietly=TRUE)))

clusterExport(cl, c("most_probable_value", "variables_names", "individus_simple"), envir=environment())

individus_simple_naissance <- pblapply(groupes_naissance, simplifier_personne, df_corresp = corresp, cl = cl, weights=TRUE)
individus_simple_naissance <- rbindlist(individus_simple_naissance)
stopCluster(cl)

Nous rajoutons ensuite les individus mis de côté précédemment.

groupes_un <- corresp_naissance_nb %>% filter(n == 1) %>% magrittr::extract2("id_personne_num_groupe")
individus_simple_naissance_un <- individus_simple[id_personne_num %in% groupes_un]
individus_simple_naissance <- rbind(individus_simple_naissance_un, individus_simple_naissance)
individus_simple <- individus_simple_naissance

1.3.6 Troisième simplification

La troisième simplification ressemble en tous points à la seconde, à l’exception de l’étape de repérage, dans lequel des conditions différentes sont utilisées pour indiquer si deux individus désignent une même personne. Nous ne présentons donc que les deux fonctions qui changent à cette étape, afin de faciliter la lecture.

Cette troisième simplification vise à corriger des erreurs de renseignements dans le sexe des individus.

Pour que deux individus soient désignés identiques, ils doivent avoir au moins un parent en commun, être nés, mariés ou décédés dans le même département. À nouveau, il faut que leur nom soit à peu près identique (nous faisons de nouveau appel à la fonction ecremer()).

#' trouver_corresp_indice_sexe
#' Pour un numero de ligne du tableau individus
#' retourne sous forme de table les correspondances entre individus.
#' Chaque ligne du resultat correspond a l'identifiant d'une personne
#' et au plus petit identifiant d'une correspondance trouvee dans la base.
# i <- 199 ; i <- 32
trouver_corresp_indice_sexe <- function(i){
  personne_cour <- individus_3[i, ]
  
  resul <- data.table(id_personne_num = personne_cour$id_personne_num,
                      id_personne_num_groupe = personne_cour$id_personne_num)
  
  # Les identifiants des parents
  personne_cour_id_num_mere <- personne_cour$id_personne_num_mere
  personne_cour_id_num_pere <- personne_cour$id_personne_num_pere
  
  # Si les identifiants des deux parents sont manquants, on ne va pas plus loin
  if(!(is.na(personne_cour_id_num_mere) & is.na(personne_cour_id_num_pere))){
    # On a un identifiant pour: mere U pere
    id_personne_num <- personne_cour$id_personne_num
    
    # Il faut que les candidats soient du meme departement
    if(!is.na(personne_cour$date_N)){
      ind_N <- which(individus_3$date_N %in% personne_cour$date_N)
    }else{
      ind_N <- NULL
    }
    
    if(!is.na(personne_cour$date_M)){
      ind_M <- which(individus_3$date_M %in% personne_cour$date_M)
    }else{
      ind_M <- NULL
    }
    
    if(!is.na(personne_cour$date_D)){
      ind_D <- which(individus_3$date_D %in% personne_cour$date_D)
    }else{
      ind_D <- NULL
    }
    
    ind <- c(ind_N, ind_M, ind_D) %>% unique()
    
    candidats <- individus_3[ind]
    
    if(!is.na(personne_cour$dept_N)){
      candidats <- 
        candidats[dept_N %in% personne_cour$dept_N]
    }
    
    candidats <- candidats[!id_personne_num %in% personne_cour$id_personne_num]
    
    candidats <- ecremer(candidats, personne_cour = personne_cour, threshold = 0.02)
    
    if(nrow(candidats)>0){
      id_retenir <- c(id_personne_num, candidats$id_personne_num)
      resul <- data.table(id_personne_num = id_retenir, id_personne_num_groupe = my_min(id_retenir))
    }# Fin du if()
  }# Fin du if sur les parents
  
  # }# Fin du if sur la date de naissance
  
  resul
}# End of trouver_corresp_indice_sexe()

Le reste du code est similaire à celui de la seconde simplification.

1.3.7 Quatrième simplification

La quatrième simplification vise à corriger les erreurs concernant les dates. Certaines dates indiquent une année, un mois et un jour. D’autres ne sont pas précises sur le mois ou le jour, et en lieu du nombre correspondant au mois ou au jour, la valeur 00 est indiquée. La simplification présente tente de prendre en compte cet aspect.

À nouveau, la fonction de repérage de doublons est l’élément changeant ici. Nous rajoutons comme contrainte pour le rapprochement, le fait que les individus doivent être nés relativement proches les uns des autres (<5km) tout en restant dans le même département (si ce dernier est indiqué).

#' trouver_corresp_indice_00
#' Pour un numero de ligne du tableau individus
#' retourne sous forme de table les correspondances entre individus.
#' Chaque ligne du resultat correspond a l'identifiant d'une personne
#' et au plus petit identifiant d'une correspondance trouvee dans la base.
trouver_corresp_indice_00 <- function(i){
  personne_cour <- individus_3[i, ]
  
  resul <- data.table(id_personne_num = personne_cour$id_personne_num,
                      id_personne_num_groupe = personne_cour$id_personne_num)
  
  # Si la date de naissance est entiere, on ne va pas chercher plus loin
  if(personne_cour$date_N_mois == 0 | personne_cour$date_N_jour == 0){
    # Les identifiants des parents
    personne_cour_id_num_mere <- personne_cour$id_personne_num_mere
    personne_cour_id_num_pere <- personne_cour$id_personne_num_pere
    
    # Si les identifiants des deux parents sont manquants, on ne va pas plus loin
    if(!(is.na(personne_cour_id_num_mere) & is.na(personne_cour_id_num_pere))){
      # On a un identifiant pour: mere U pere
      
      id_personne_num <- personne_cour$id_personne_num
      
      # Il faut que les candidats soient de la même année et si possible du même département
      candidats <- individus_3[date_N_annee %in% personne_cour$date_N_annee]
      
      if(!is.na(personne_cour$dept_N)){
        candidats <- 
          candidats[dept_N %in% personne_cour$dept_N]
      }
      
      candidats <- candidats[!id_personne_num %in% personne_cour$id_personne_num]
      
      if(nrow(candidats)>0){
        les_distances <- 
          geosphere::distGeo(c(as.numeric(personne_cour$long_N), as.numeric(personne_cour$lat_N)),
                             cbind(as.numeric(candidats$long_N), as.numeric(candidats$lat_N)))
        candidats <- candidats[which(les_distances < 5000),]
        
        candidats <- ecremer(candidats, personne_cour = personne_cour, threshold = 0.01, prenoms = TRUE)
      }
      
      
      if(nrow(candidats)>0){
        id_retenir <- c(id_personne_num, candidats$id_personne_num)
        resul <- data.table(id_personne_num = id_retenir, id_personne_num_groupe = my_min(id_retenir))
      }# Fin du if()
    }# Fin du if sur les parents
  }
  # }
  resul
}# End of trouver_corresp_indice_00()

1.3.8 Cinquième simplification

Dans cette cinquième simplification, nous effectuons le rapprochement par les prénoms, en utilisant le prénom uniquement, et plus le nom de famille. La fonction ecremer() est remplacée par la suivante :

#' ecremer
#' Retourne un tableau de candidats possibles
#' @candidats: (data.frame) table dans laquelle chercher les candidats
ecremer_2 <- function(candidats, threshold = 0.1, personne_cour){
  # Ecremer par le nom
  ind_nom <- which(calculer_distance(personne_cour = personne_cour, candidats = candidats, variable = "nom", threshold = threshold))
  candidats <- candidats[ind_nom,]
  
  if(nrow(candidats)>0){
    # Par le prenom
    pren_cour <- personne_cour$prenom
    ind_prenom <- which(str_detect(string = candidats$prenoms,
                                   pattern = fixed(pren_cour, ignore_case = TRUE)))
    candidats <- candidats[ind_prenom,]
    
    if(nrow(candidats)>0 & !is.na(personne_cour$nom_mere)){
      ind_mere_nom <- calculer_distance(personne_cour = personne_cour, candidats = candidats, variable = "nom_mere", threshold = threshold)
      ind_mere_nom <- c(which(is.na(ind_mere_nom)), which(ind_mere_nom))
      candidats <- candidats[ind_mere_nom,]
    }# Fin if mere
    
    if(nrow(candidats)>0 & !is.na(personne_cour$nom_pere)){
      ind_pere_nom <- calculer_distance(personne_cour = personne_cour, candidats = candidats, variable = "nom_pere", threshold = threshold)
      ind_pere_nom <- c(which(is.na(ind_pere_nom)), which(ind_pere_nom)) %>% sort()
      candidats <- candidats[ind_pere_nom,]
    }# Fin if pere
    
    
  }# Fin du if
  candidats
}# End of ecremer_2()

La fonction de repérage est quant à elle modifiée comme suit. Les personnes doivent être nées le même jour, pas trop loin l’une de l’autre (<5km), leurs parents doivent avoir des noms proches.

#' trouver_corresp_indice_fifth
#' Pour un numero de ligne du tableau individus
#' retourne sous forme de table les correspondances entre individus.
#' Chaque ligne du resultat correspond a l'identifiant d'une personne
#' et au plus petit identifiant d'une correspondance trouvee dans la base.
trouver_corresp_indice_fifth <- function(i){
  personne_cour <- individus_3[i, ]
  
  resul <- data.table(id_personne_num = personne_cour$id_personne_num,
                      id_personne_num_groupe = personne_cour$id_personne_num)
  
  id_personne_num <- personne_cour$id_personne_num
  
  # Il faut que les candidats soient nes le meme jour
  ind <- which(individus_3$date_N %in% personne_cour$date_N)
  candidats <- individus_3[ind]
  
  if(!is.na(personne_cour$dept_N)){
    candidats <- 
      candidats[dept_N %in% personne_cour$dept_N]
  }
  
  candidats <- candidats[!id_personne_num %in% personne_cour$id_personne_num]
  
  if(nrow(candidats)>0){
    les_distances <- 
      geosphere::distGeo(c(as.numeric(personne_cour$long_N), as.numeric(personne_cour$lat_N)),
                         cbind(as.numeric(candidats$long_N), as.numeric(candidats$lat_N)))
    candidats <- candidats[which(les_distances < 5000),]
    
    candidats <- ecremer_2(candidats, personne_cour = personne_cour, threshold = 0.01)
  }
  
  
  if(nrow(candidats)>0){
    id_retenir <- c(id_personne_num, candidats$id_personne_num)
    resul <- data.table(id_personne_num = id_retenir, id_personne_num_groupe = my_min(id_retenir))
  }# Fin du if()
  
  
  # }
  resul
}# End of trouver_corresp_indice_fifth()

1.3.9 Sixième simplification

La dernière simplification s’appuie sur les frères et soeurs.

Nous rajoutons d’abord les informations relatives aux parents des individus.

individus_2 <- individus_simple
individus_2$id_personne_num_mere <- registre$id_personne_num_mere[match(individus_simple$id_personne_num, registre$id_personne_num)]
  individus_2$id_personne_num_pere <- registre$id_personne_num_pere[match(individus_simple$id_personne_num, registre$id_personne_num)]

Nous définissons la fonction trouver_enfants() qui comme son nom l’indique, permet de trouver les identifiants des enfants pour un individu situé à la position i du tableau individus_2.

trouver_enfants <- function(i){
  id_cour <- individus_2$id_personne_num[i]
  ind_enfants <- c(which(individus_2$id_personne_num_mere == id_cour),
                   which(individus_2$id_personne_num_pere == id_cour)) %>%
    unique()
  # ids_enfants <- individus_2$id_personne_num[ind_enfants] %>% str_c(., collapse = "|") %>%
  #   str_c("|", ., "|")
  ids_enfants <- individus_2$id_personne_num[ind_enfants]
  ids_enfants
}# Fin de trouver_enfants()

Nous appliquons cette fonction à chaque observation de la base, en parallélisant l’exécution du code.

library(parallel)
  ncl <- detectCores()-1
  (cl <- makeCluster(ncl))

  invisible(clusterEvalQ(cl, library(dplyr, warn.conflicts=FALSE, quietly=TRUE)))
  invisible(clusterEvalQ(cl, library(stringr, warn.conflicts=FALSE, quietly=TRUE)))
  invisible(clusterEvalQ(cl, library(tidyr, warn.conflicts=FALSE, quietly=TRUE)))
  invisible(clusterEvalQ(cl, library(stringi, warn.conflicts=FALSE, quietly=TRUE)))
  invisible(clusterEvalQ(cl, library(data.table, warn.conflicts=FALSE, quietly=TRUE)))
  invisible(clusterEvalQ(cl, library(stringdist, warn.conflicts=FALSE, quietly=TRUE)))
  
  clusterExport(cl, c("individus_2"), envir=environment())
  
  ids_enfants <- pbsapply(seq(1, nrow(individus_2)), trouver_enfants, cl = cl)
  # stopCluster(cl)

Notons que nous ne fermons pas la connexion aux clusters, nous allons les réutiliser un peu plus loin.

Nous rajoutons à chaque individu de la base les identifiants trouvés pour leurs enfants. Puis, nous recréons une variable pour le prénom des individus, en ne gardant que le premier prénom. “Marie Chantal Jacqueline Michu”" aura comme prénom court “Marie”.

individus_2$ids_enfants <- ids_enfants
  
individus_2 <- 
  individus_2 %>% 
  mutate(date_N_annee = str_sub(date_N, 1, 4) %>% as.numeric(),
         date_D_annee = str_sub(date_D, 1, 4) %>% as.numeric(),
         prenom_court = str_replace_all(prenoms, "\\((.*?)\\)", ""),
         prenom_court = stri_replace_all_regex(prenom_court, "'|\\(|\\)|\\?|_|\\*|&gt|\\^|\\+|[[:digit:]]|\\.|<|>", "") %>% 
           str_trim(),
         prenom_court = gsub(" .*$", "", prenom_court),
         prenom_court = stri_trans_general(prenom_court, "Latin-ASCII"),
         prenom_court = str_extract(prenom_court, '\\w*') %>% 
           stri_trim() %>% 
           str_to_lower()) %>% 
  data.table()

Nous créons ensuite la fonction same_siblings() qui permet de repérer les personnes dédoublonnées, en fonction de leurs frères et soeurs.

#' same_siblings
#' Retourne les freres/soeurs en doublons pour un individu de la base
#' i <- ind_faire[55]
#' i <- 2407
#' rm(individu,ids_parents, parents,ids_siblings,res,siblings,prenom_individu,prenom_individu,candidats, i)
same_siblings <- function(i){
  individu <- individus_2[i, ]
  # Rechercher si cette personne a des freres et soeurs
  ids_parents <- c(individu$id_personne_num_mere, individu$id_personne_num_pere)
  ids_parents <- ids_parents[!is.na(ids_parents)]
  
  parents <- individus_2[id_personne_num %in% ids_parents]
  ids_siblings <- unlist(parents$ids_enfants) %>% unique()
  ids_siblings <- ids_siblings[-which(ids_siblings == individu$id_personne_num)]
  
  res <- NULL
  
  
  # Si la personne a des freres et soeurs
  if(length(ids_siblings)>0){
    siblings <- individus_2[id_personne_num %in% ids_siblings]
    nom_individu <- individu$nom
    prenom_court_individu <- individu$prenom_court
    # Reperer parmi eux ceux qui ont le meme prenom
    candidats <- siblings[prenom_court == prenom_court_individu & nom == nom_individu]
    
    # S'il y a une dates de deces indiquee, elle doit correspondre (meme annee)
    # Si la date est manquante, on emet l'hypothese qu'il peut toujours s'agir de la meme personne
    if(nrow(candidats)>0){
      annee_D_individu <- individu$date_D_annee
      if(is.na(annee_D_individu)){
        annee_D_potentielle <- NA
      }else{
        annee_D_potentielle <- c(NA, annee_D_individu)
      }
      
      candidats <- candidats[date_D_annee %in% annee_D_potentielle]
      
      if(nrow(candidats)>0){
        annee_N_individu <- individu$date_N_annee
        if(is.na(annee_N_individu)){
          annee_N_potentielle <- NA
        }else{
          annee_N_potentielle <- c(NA,annee_N_individu)
        }
        
        candidats <- candidats[date_N_annee %in% annee_N_potentielle]
        
        # A present, il reste des freres ou soeurs avec le meme prenom et la meme annee de naissance et de deces
        if(nrow(candidats)>0){
          res <- c(individu$id_personne_num, candidats$id_personne_num)
        }
        
      }# Fin du if pour la date de naissance
      
    }# Fin du if pour la date de deces
    
    
  }# Fin de s'il y a des siblings
  
  res
}# Fin de same_siblings()

Reste à appliquer cette fonction, en utilisant à nouveau les clusters que nous n’avons pas révoqués.

ind_ne_pas_faire <- c(which(is.na(individus_2$id_personne_num_mere) & is.na(individus_2$id_personne_num_pere)))

ind_faire <- seq(1, nrow(individus_2))
ind_faire <- ind_faire[-ind_ne_pas_faire]

# Calcul en parallele
clusterExport(cl, c("individus_2"), envir=environment())
fusion_siblings <- pblapply(ind_faire, same_siblings, cl = cl)
stopCluster(cl)

# Il faut fusionner ces personnes
fusion_siblings <- fusion_siblings[which(!sapply(fusion_siblings, is.null))]

corresp_siblings <- lapply(fusion_siblings, function(x) data.table(id_personne_num = x, id_personne_num_groupe = min(x))) %>% 
  rbindlist() %>% 
  unique()

corresp_siblings <- corresp_siblings[, list(id_personne_num_groupe = min(id_personne_num_groupe, na.rm=T)), by = id_personne_num] %>% unique()

# On repercute ces nouvelles informations sur corresp
corresp <- 
    corresp %>% 
    left_join(corresp_siblings %>% rename(id_personne_num_groupe_new = id_personne_num_groupe)) %>% 
    mutate(id_personne_num_groupe = pmin(id_personne_num_groupe, id_personne_num_groupe_new, na.rm=T)) %>%
    select(-id_personne_num_groupe_new) %>% 
    data.table()

Reste ensuite à boucler sur le registre et à appliquer la fonction de simplification, comme dans les étapes précédentes.

1.4 Mise en commun des données départementales

Dans la section précédente, nous avons rassemblé les arbres par département. Il est question dans cette nouvelle section de mettre en commun les bases départementales, afin de créer une base française.

1.4.1 Chargement des données

Chargeons la liste des départements que nous avions créée auparavant et rajoutons une colonne avec les numéros du département.

load("liste_depts.rda")
liste_depts <- 
  liste_depts %>% 
  mutate(dep_num = str_sub(dept, 2))

Nous pouvons ajouter le nom de la région, en s’appuyant sur des données contenues dans le package raster.

map_fr_sp <- raster::getData(name = "GADM", country = "FRA", level = 2)
liste_depts <- 
  liste_depts %>% 
  left_join(
    data.frame(dep_num = map_fr_sp$CCA_2, dep_name = map_fr_sp$NAME_2)
  ) %>% 
  tbl_df()

Nous chargeons ensuite les données de généalogie départementales obtenues à l’étape précédente.

#' charger_corresp
#' Charger la table de correspondance, pour un département donné
#' @i_departement: (int) ligne de liste_depts à traiter
charger_corresp <- function(i_departement){
  departement <- liste_depts$dept[i_departement]
  load(str_c("../data/dept/migration/", departement, "/resul.rda"))
  resul$corresp
}# Fin de charger_corresp()

# Charger les données
corresp <- pblapply(seq(1, nrow(liste_depts)), charger_corresp)
corresp <- bind_rows(corresp)

Les individus sont identifiés à l’aide d’une valeur numérique, propre à chaque département, ainsi qu’une chaîne de caractères composée du nom d’identifiant du généalogiste amateur, du nom de famille de l’individu et de ses prénoms. Tandis qu’une personne sera désignée de manière unique par son identifiant textuel, son identidiant numérique pourra être différent en fonction du département. Nous avons en effet travaillé par département, en définissant à ce moment un identifiant pour chaque individu présent dans l’extrait de la base. Maintenant que nous réunissons les départements, nous devons créer un nouvel identifiant numérique. Prenons un exemple : l’individu Pierre Dupond présent dans l’arbre de l’utilisateur batman36 a pour mère une femme née en 1801 dans le Finistère et pour père un homme né en 1800 dans le Morbihan. Aussi, cet individu a reçu un identifiant numérique lorsque nous avons traité le Finistère, et un autre identifiant numérique lorsque nous avons traité le Morbihan. Son identifiant textuel est en revanche identique dans les deux sous-bases. Nous allons ainsi attribuer un nouvel identifiant numérique à ce Pierre Dupond, en nous appuyant sur son identifiant textuel. Nous faisons cela puisque nous allons à nouveau tenter de repérer des associations qui auraient pu nous échaper sur les sous bases, et tenter de compléter certaines informations.

corresp_id <- data.frame(id_personne = unique(corresp$id_personne)) %>% tbl_df()
corresp_id <- corresp_id %>% mutate(id_personne_num_new = row_number())

Maintenant que nous disposons de nouveaux identifiants pour les individus, nous pouvons charger les bases départementales d’individus, mettre à jour leur identifiant numérique, et joindre les bases. Nous le faisons en faisant appel à une fonction que nous appliquons à chaque département.

#' charger_individus_simple
#' Retourne une liste pour un département contenant :
#' - registre : le registre mis à jour
#' - individus_simple : les informations sur les individus
#' - corresp : le tableau de correspondances mis à jour (à chaque personne présente initialement dans la base, l'identifiant numérique correspondant)
#' - departement : le numéro de département
#' @i_departement: (int) ligne de liste_depts à consulter pour récupérer le numéro de département
charger_individus_simple <- function(i_departement){
  departement <- liste_depts$dept[i_departement]
  load(str_c("../data/dept/migration/", departement, "/resul.rda"))
  
  # TABLE DE CORRESPONDANCE #
  
  corresp <- resul$corresp %>% tbl_df()
  # Rajout des identifiants sous forme de texte dans la table de correspondance
  corresp_2 <- 
    corresp %>% 
    left_join(
      corresp %>% 
        select(id_personne, id_personne_num) %>% 
        rename(id_personne_groupe = id_personne, id_personne_num_groupe = id_personne_num) %>% 
        unique(),
      by = "id_personne_num_groupe"
    )
  corresp_2 <- corresp_2 %>% filter(!is.na(id_personne_groupe))
  
  # Rajout des nouveaux identifiants numeriques
  corresp_2 <- 
    corresp_2 %>% left_join(corresp_id %>% mutate(id_personne = as.character(id_personne)), by = "id_personne")
  
  
  corresp_2 <- 
    corresp_2 %>% 
    left_join(
      corresp_id %>% rename(id_personne_groupe = id_personne, id_personne_num_groupe_new = id_personne_num_new) %>% 
        mutate(id_personne_groupe = as.character(id_personne_groupe)),
      by = "id_personne_groupe"
    )
  
  # INDIVIDUS #
  
  individus_simple <- resul$individus_simple
  individus_simple <- individus_simple[weight != 0,]
  individus_simple <- 
    individus_simple %>% tbl_df() %>% 
    left_join(
      corresp_2 %>% 
        select(id_personne_num, id_personne_num_new),
      by = "id_personne_num"
    ) %>% 
    tbl_df()
  
  # REGISTRE DES LIENS DE PARENTE #
  
  registre <- resul$registre %>% data.table()
  registre <- registre[id_personne_num %in% individus_simple$id_personne_num] %>% 
    tbl_df()
  
  registre <- 
    registre %>% 
    left_join(corresp_2 %>% select(id_personne_num, id_personne_num_new), by = "id_personne_num") %>% 
    # Ajout mere
    left_join(corresp_2 %>% select(id_personne_num, id_personne_num_new) %>% 
                rename(id_personne_num_mere = id_personne_num, id_personne_num_mere_new = id_personne_num_new),
              by = "id_personne_num_mere") %>% 
    # Ajout pere
    left_join(corresp_2 %>% select(id_personne_num, id_personne_num_new) %>% 
                rename(id_personne_num_pere = id_personne_num, id_personne_num_pere_new = id_personne_num_new),
              by = "id_personne_num_pere")
  
  registre <- 
    registre %>% 
    select(id_personne_num_new, id_personne_num_mere_new, id_personne_num_pere_new)
  
  corresp_2 <- corresp_2 %>% 
    select(id_personne, id_personne_num_new, id_personne_num_groupe_new) %>% 
    filter(!is.na(id_personne_num_new))
  
  registre <- 
    registre %>% 
    rename(id_personne_num = id_personne_num_new,
           id_personne_num_mere = id_personne_num_mere_new,
           id_personne_num_pere = id_personne_num_pere_new)
  
  individus_simple <- 
    individus_simple %>% 
    select(-id_personne_num) %>% 
    rename(id_personne_num = id_personne_num_new)
  
  corresp_2 <- 
    corresp_2 %>% 
    rename(id_personne_num = id_personne_num_new,
           id_personne_num_groupe = id_personne_num_groupe_new)
  
  list(registre = registre, individus_simple = individus_simple, corresp = corresp_2, departement = departement)
  
}# Fin de charger_individus_simple()

res <- pblapply(seq(1, nrow(liste_depts)), charger_individus_simple)

Nous pouvons extraire et joindre entre-eux les tableaux de correspondance.

# Table des correspondances
corresp <- pblapply(res, function(x) x$corresp)
corresp <- corresp %>% rbindlist() %>% unique()

corresp <- corresp[, list(id_personne_num_groupe = min(id_personne_num_groupe, na.rm=T)), by = id_personne_num] %>% unique()
corresp_id <- 
  corresp_id %>% 
  rename(id_personne_num = id_personne_num_new) %>% 
  data.table()

corresp <- corresp_id[corresp, on = "id_personne_num"]

Idem pour les tableaux contenant les informations relatives à chaque individu.

individus_simple <- pblapply(res, function(x) x$individus_simple)
individus_simple <- individus_simple %>% rbindlist() %>% unique()
individus_simple <- individus_simple[!is.na(id_personne_num)]

1.4.2 Nouvelle simplification

Nous recherchons les associations évidentes qu’il est possible de faire avec ces nouvelles données. Pour cela, nous regardons l’identifiant textuel des individus et voyons si nous pouvons dénombrer des valeurs dupliquées.

corresp_2 <- corresp[, list(id_personne_num, id_personne_num_groupe)] %>% 
  rename(id_personne_num_groupe_new = id_personne_num_groupe,
    id_personne_num_groupe = id_personne_num)

setkey(corresp, id_personne_num_groupe)
setkey(corresp_2, id_personne_num_groupe)
corresp <- corresp_2[corresp]
corresp[, id_personne_num_groupe := pmin(id_personne_num_groupe, id_personne_num_groupe_new, na.rm=T)]
corresp$id_personne_num_groupe_new <- NULL
rm(corresp_2)

Nous faisons de nouveau appel à une fonction pour simplifier les valeurs, basée sur la fréquence observée des valeurs.

#' most_probable_value
#' Find the most probable value for a variable
#' using the values found within all candidates
#' (excluding NAs)
#' @variable_name: (string) name of the variable
#' @df: (data.table)
#' @weights: (logic) should the simplification consider the weights? (default: FALSE)
#' variable_name <- "date_N"
most_probable_value <- function(df, variable_name, weights = FALSE){
  valeurs <- df[[variable_name]]
  # valeurs <- individu_cour_a_fusionner %>% magrittr::extract2(variable_name)
  if(!all(is.na(valeurs))){
    if(weights){
      poids <- df[["weight"]]
      res <- data.frame(valeurs, poids, stringsAsFactors = FALSE) %>% 
        group_by(valeurs) %>% 
        summarise(poids = sum(poids)) %>% 
        ungroup()
      
      # Si la variable est une date
      # on va privilegier les informatins relatives aux dates completes
      if(variable_name %in% c("date_N", "date_M", "date_D")){
        tmp <- res %>% 
          mutate(mois = str_sub(valeurs, 5, 6),
                 jour = str_sub(valeurs, 7, 8)) %>% 
          filter(!(mois == "00" | jour == "00"))
        # S'il reste des informations, on se base sur elles
        if(nrow(tmp)>0) res <- tmp
      }
      
      res <- 
        res %>%
        arrange(desc(poids)) %>% 
        slice(1) %>% 
        magrittr::extract2("valeurs")
    }else{
      table_freq <- sort(table(valeurs))
      res <- names(table_freq)[1]
    }
  }else{
    res <- NA
  }
  res
}# End of most_probable_value()

#' simplifier_personne
#' À partir d'un identifiant numérique, retourne un tableau à une seule ligne contenant les valeurs de chaque variable relative à cette personne.
#' @weights: (logic) should the simplification consider the weights? (default: FALSE)
simplifier_personne <- function(un_groupe_num, weights = TRUE){
  groupe_cour <- df_corresp[id_personne_num_groupe %in% un_groupe_num]
  
  individus_cour <- df_individus_simple[id_personne_num %in% groupe_cour$id_personne_num]
  
  weight <- sum(individus_cour$weight)
  
  if(nrow(individus_cour) == 1){
    nouvelles_valeurs <- individus_cour[, mget(c(variables_names))]
  }else{
    # Les valeurs probables pour les variables d'interet
    nouvelles_valeurs <- 
      variables_names %>% 
      sapply(., most_probable_value, df = individus_cour, weights = weights) %>% 
      t() %>% 
      data.table(stringsAsFactors = FALSE)
  }
  
  
  cbind(id_personne_num = un_groupe_num, nouvelles_valeurs, weight = weight)
}# End of simplifier_personne()

Nous utilisons à nouveau la fonction my_min(), qui permet de retourner la valeur NA plutôt que Inf quand toutes les valeurs comparées sont manquantes.

my_min <- function(x) ifelse( !all(is.na(x)), min(x, na.rm=T), NA)

Les variables que nous souhaitons extraire sont les suivantes :

# Les variables qui nous interessent
variables_names <- c("nom", "prenoms",
                     "sexe",
                     "lieu_N", "dept_N", "region_N", "pays_N", "lat_N", "long_N", "date_N",
                     "lieu_M", "dept_M", "region_M", "pays_M", "lat_M", "long_M", "date_M",
                     "lieu_D", "dept_D", "region_D", "pays_D", "lat_D", "long_D", "date_D",
                     "prenom")

Puis, nous effectuons la simplification des individus, de la même manière que ce que nous avons effectué dans la section précédente.

Reste alors à sauvegarder le résultat :

res <- list(nb_obs = nb_obs,
            corresp = corresp,
            registre = registre,
            individus_simple = individus_simple)

save(res, file = "migration_1800_04_dep.rda")

1.5 Création du tableau de générations

Les étapes précédents ont permis d’obtenir une base de données de généalogie contenant des informations relatives à des individus nés en France entre 1800 et 1804, leur parents et leur descendants. Nous pouvons, à partir de ces données, créer une nouvelle base permettant d’étudier les générations. Cette section donne les étapes utilisées à cet effet et présente les codes utilisés pour générer les graphiques.

Nous commençons par charger certains packages nécessaires.

library(tidyverse)
library(lubridate)
library(stringr)
library(stringi)
library(data.table)
library(pbapply)
library(dplyrExtras)
library(stringdist)
library(sp)

Nous chargeons également la liste des départements en y ajoutant leur nom.

load("liste_depts.rda")
liste_depts <- 
  liste_depts %>% 
  mutate(dep_num = str_sub(dept, 2))

map_fr_sp <- raster::getData(name = "GADM", country = "FRA", level = 2)

liste_depts <- 
  liste_depts %>% 
  left_join(
    data.frame(dep_num = map_fr_sp$CCA_2, dep_name = map_fr_sp$NAME_2)
  ) %>% 
  tbl_df()

Nous récupérons les données de généalogie obtenues à l’issue de l’étape précédente.

load("migration_1800_04_dep.rda")
registre <- res$registre
individus_simple <- res$individus_simple
individus_2 <- individus_simple

À partir des identifiants des parents présents dans chaque ligne du tableau des données individuelles, nous pouvons rajouter à chaque individu les informations relatives à leurs parents, si toutefois elles sont existantes.

# On rajoute les identifiants des parents
individus_2$id_personne_num_mere <- registre$id_personne_num_mere[match(individus_simple$id_personne_num, registre$id_personne_num)]
individus_2$id_personne_num_pere <- registre$id_personne_num_pere[match(individus_simple$id_personne_num, registre$id_personne_num)]

ind_meres <- match(individus_2$id_personne_num_mere, individus_2$id_personne_num)
ind_peres <- match(individus_2$id_personne_num_pere, individus_2$id_personne_num)

# Les meres
individus_meres <- 
  individus_2[ind_meres,] %>% 
  magrittr::set_colnames(str_c(colnames(individus_2), "_mere")) %>% 
  select(-id_personne_num_mere, -id_personne_num_mere_mere, -id_personne_num_pere_mere)

# Les peres
individus_peres <- 
  individus_2[ind_peres,] %>% 
  magrittr::set_colnames(str_c(colnames(individus_2), "_pere")) %>% 
  select(-id_personne_num_pere, -id_personne_num_mere_pere, -id_personne_num_pere_pere)

# On ajoute les informations relatives aux parents
individus_2 <- 
  individus_2 %>% 
  cbind(individus_meres) %>% 
  cbind(individus_peres) %>% 
  # select(-id_mere, -id_pere) %>% 
  data.table()

Concentrons-nous sur des individus nés entre 1680 et 2017.

# On retire les individus nes apres 2017 et avant 1680
individus_2 <- individus_2[date_N_annee > 1680 & date_N_annee <= 2017]

1.5.1 Retrait des anomalies

Il existe des observations avec des valeurs aberrantes. Nous alons les identifier et les retirer.

Une source d’anomalie provient de l’âge auquel une personne peut avoir un enfant. Si nous constations qu’un enfant est né avant qu’un de ses parents a atteint l’âge de 13 ans, nous considérons qu’il s’agit d’une erreur de saisie (ou de fusion, ce qui est plus problématique). Nous considérons qu’un individu est en capacité de se reproduire entre ses 13 et 70 ans.

# On retire les anomalies concernant les dates de naissances
age_min_nenf <- 13
age_max_enf <- 70

individus_2 <-
  individus_2 %>% 
  mutate(date_N_annee = str_sub(date_N, 1, 4) %>% as.numeric(),
         date_D_annee = str_sub(date_D, 1, 4) %>% as.numeric(),
         date_N_annee_mere = str_sub(date_N_mere, 1, 4) %>% as.numeric(),
         date_N_annee_pere = str_sub(date_N_pere, 1, 4) %>% as.numeric()) %>% 
  # Anomalie si la mere a eu son enfant a moins de 13 ans ou a plus de 70 ans
  mutate(anomalie_mere = (date_N_annee <= (date_N_annee_mere+age_min_nenf)) | (date_N_annee >= (date_N_annee_mere + age_max_enf)),
         anomalie_mere = ifelse(is.na(anomalie_mere), yes = FALSE, no = anomalie_mere)) %>% 
  # Idem pour le pere
  mutate(anomalie_pere = (date_N_annee <= (date_N_annee_pere+age_min_nenf)) | (date_N_annee >= (date_N_annee_pere + age_max_enf)),
         anomalie_pere = ifelse(is.na(anomalie_pere), yes = FALSE, no = anomalie_pere)) %>% 
  mutate(anomalie = anomalie_mere + anomalie_pere) %>% 
  filter(anomalie == 0) %>% 
  select(-anomalie_mere, -anomalie_pere, -anomalie) %>% 
  data.table()

Nous retirons les observations pour lesquelles des âges négatifs ou supérieurs à 122 ans sont renseignés.

individus_2 <-
  individus_2 %>% 
  # mutate(age = interval(ymd(date_N), ymd(date_D)) / years(1)) %>% 
  mutate(age = date_D_annee - date_N_annee) %>% 
  mutate(keep = age >= 0 & age <= 122,
         keep = ifelse(is.na(keep), yes = TRUE, no = keep)) %>% 
  filter(keep == TRUE) %>% 
  select(-keep)

Nous retirons les observations pour lesquelles l’année de décès survient avant cette de la naissance

individus_2 <- 
  individus_2 %>% 
  mutate(keep = (str_sub(date_D, 1, 4) %>% as.numeric) >= (str_sub(date_N, 1, 4) %>% as.numeric),
         keep = ifelse(is.na(keep), yes = TRUE, no = keep)) %>% 
  filter(keep == TRUE) %>% 
  select(-keep)

Filtrons les individus de la base individus en ne retenant que ceux qui ont été obtenus dans individus_2.

individus_simple <- 
  individus_simple[id_personne_num %in% individus_2$id_personne_num]

1.5.2 Les générations

1.5.2.1 Aïeux

Nous pouvons à présent construire le tableau de générations, en commençant par les individus à l’origine de l’étude, c’est-à-dire ceux qui sont nés entre 1800 et 1804, en France.

gen_0 <- individus_2[date_N_annee %in% seq(1800, 1804)]
gen_0 <- registre[id_personne_num %in% gen_0$id_personne_num]

Nous allons construire cet tableau génération par génération. Nous définissons à cet effet un tableau contenant la génération courante, qui sera modifé quand on s’intéressera aux enfants, petits-enfants, etc. Il s’agit de l’initialiser.

gen_n <- copy(gen_0)
gen_n$generation <- 0

Nous retirons de notre tableau d’individus les individus de la génération courante.

gen_restants <- registre[id_personne_num %in% individus_2$id_personne_num]
gen_restants <- gen_restants[!id_personne_num %in% gen_n$id_personne_num]

Nous initialisons notre tableau de génération, composé pour l’instant des individus de la génération courante, i.e., la génération 0.

generations <- 
  data.frame(id_personne_num = gen_0$id_personne_num) %>% 
  rename(id_aieul = id_personne_num) %>% 
  data.table()

1.5.2.2 Enfants

Identifionts à présent les enfants de la génération courante.

peres_potentiels <- 
  gen_restants %>% 
  select(id_personne_num, id_personne_num_pere) %>% 
  rename(id_aieul = id_personne_num_pere, id_enfant = id_personne_num) %>% 
  filter(!is.na(id_aieul))

meres_potentiels <- 
  gen_restants %>% 
  select(id_personne_num, id_personne_num_mere) %>% 
  rename(id_aieul = id_personne_num_mere, id_enfant = id_personne_num) %>% 
  filter(!is.na(id_aieul))

parents_potentiels <- 
  peres_potentiels %>% 
  rbind(meres_potentiels)

setkey(parents_potentiels, id_aieul)
setkey(generations, id_aieul)

generations <- parents_potentiels[generations]

1.5.2.3 Petits-enfants

Puis les petits-enfants.

peres_potentiels <- 
  gen_restants %>% 
  select(id_personne_num, id_personne_num_pere) %>% 
  rename(id_enfant = id_personne_num_pere, id_p_enfant = id_personne_num) %>% 
  filter(!is.na(id_enfant))

meres_potentiels <- 
  gen_restants %>% 
  select(id_personne_num, id_personne_num_mere) %>% 
  rename(id_enfant = id_personne_num_mere, id_p_enfant = id_personne_num) %>% 
  filter(!is.na(id_enfant))

parents_potentiels <- 
  peres_potentiels %>% 
  rbind(meres_potentiels)


setkey(parents_potentiels, id_enfant)
setkey(generations, id_enfant)

generations <- parents_potentiels[generations]

1.5.2.4 Arrière-petits-enfants

Et enfin, les arrière-petits-enfants.

peres_potentiels <- 
  gen_restants %>% 
  select(id_personne_num, id_personne_num_pere) %>% 
  rename(id_p_enfant = id_personne_num_pere, id_ap_enfant = id_personne_num) %>% 
  filter(!is.na(id_p_enfant))

meres_potentiels <- 
  gen_restants %>% 
  select(id_personne_num, id_personne_num_mere) %>% 
  rename(id_p_enfant = id_personne_num_mere, id_ap_enfant = id_personne_num) %>% 
  filter(!is.na(id_p_enfant))

parents_potentiels <- 
  peres_potentiels %>% 
  rbind(meres_potentiels)

setkey(parents_potentiels, id_p_enfant)
setkey(generations, id_p_enfant)

generations <- parents_potentiels[generations]

1.5.3 Conversion des départements en polygones

Nous avons besoin d’identifier dans quel département tombent certains points identifiés par des coordonnées (longitude,latitude). Pour cela, une manière simple consiste à utiliser la fonction inside.gpc.poly() du package surveillance. Il faut toutefois fournir un polygone de la classe gpc.poly(). Pour l’heure, nous disposons de données pour tracer les départements sous forme de tableaux de données. Nous convertissons de fait ces tableaux en objets gpc.poly().

library(rgeos)
#' group_to_gpc_poly
#' Returns the polygon from a group within the data frame containing
#' the coordinates of regions
#' @df: data frame with the coordinates of the polygons
#' @group: (string) group ID
group_to_gpc_poly <- function(df, group){
  group_coords <- df[df$group == group, c("long", "lat")]
  as(group_coords, "gpc.poly")
}

#' departement_polygone
#' Returns a French region as a gpc.poly object
#' un_dep: (string) code of the region
departement_polygone <- function(un_dep){
  dep_tmp <-
    france_1800 %>%
    filter(departement == un_dep)
  df_tmp_polys <- lapply(unique(dep_tmp$group), group_to_gpc_poly, df = dep_tmp)
  df_tmp_polys_tot <- df_tmp_polys[[1]]
  if(length(df_tmp_polys)>1){
    for(i in 2:length(df_tmp_polys)){
      df_tmp_polys_tot <- gpclib::union(df_tmp_polys_tot,  df_tmp_polys[[i]])
    }
  }
  df_tmp_polys_tot
}# End of departement_polygone()

france_1800_gpc <- lapply(unique(france_1800$departement), departement_polygone)
names(france_1800_gpc) <- unique(france_1800$departement)
save(france_1800_gpc, file = "france_1800_gpc.rda")

Maintenant que nous disposons des départements sous forme de gpc.poly, nous pouvons créer une fonction permettant d’indiquer si un point (lon,lat) se situe à l’intérieur d’un département donné.

est_dans_poly_naissance <- function(longitude, latitude, departement){
  surveillance::inside.gpc.poly(x = longitude,
                                y = latitude,
                                departement,
                                mode.checked = FALSE)
}

1.5.4 Repérage des départements de naissance des individus

Pour repérer le département de naissance des individus, nous convertissons les coordonnées en valeurs numériques. Elles sont pour l’heure en chaînes de caractères.

individus_simple_reste <- individus_simple
individus_simple_reste <- individus_simple_reste[!is.na(lat_N)]
individus_simple_reste <- 
  individus_simple_reste %>% 
  mutate(long_N = as.numeric(long_N),
         lat_N = as.numeric(lat_N))

individus_simple_reste <- data.table(individus_simple_reste)

Nous allons effectuer une boucle sur les départements pour déterminer, à chaque itération, quels sont les individus nés dans ce département. Comme nous avons un nombre conséquent d’invididus et que le territoire français est relativement grand à l’échelle des départements, nous procédons par un filtrage grossier de nos données à chaque tour. Pour un département donné, nous récupérons les coordonnées d’un rectangle le contenant, à l’aide de la fonction get.bbox(). Nous filtrons ensuite les individus en fonction de ce rectangle, pour écarter ceux dont les coordonnées sont à l’extérieur. Reste alors à appliquer la fonction est_dans_poly_naissance() sur les points filtrés, qui sont bien moins nombreux qu’à l’origine.

dept_n_indiv <- rep(list(NA), length(france_1800_gpc))

pb <- txtProgressBar(min = 1, max = length(france_1800_gpc), style = 3)
for(i in 1:length(france_1800_gpc)){
  poly <- france_1800_gpc[[i]]
  nom_dep <- names(france_1800_gpc)[i]
  
  individus_simple_reste_2 <- 
    individus_simple_reste %>% 
    select(long_N, lat_N, id_personne_num) %>% 
    filter(long_N >= get.bbox(poly)$x[1],
           long_N <= get.bbox(poly)$x[2],
           lat_N >= get.bbox(poly)$y[1],
           lat_N <= get.bbox(poly)$y[2]) %>% 
    rowwise() %>% 
    mutate(in_poly = est_dans_poly_naissance(long_N, lat_N, poly)) %>% 
    select(id_personne_num, in_poly) %>% 
    mutate(dept_N = nom_dep) %>% 
    filter(in_poly == TRUE) %>% 
    select(-in_poly)
  
  
  individus_simple_reste <- 
    individus_simple_reste[! id_personne_num %in% individus_simple_reste_2$id_personne_num]
  
  dept_n_indiv[[i]] <- individus_simple_reste_2
  setTxtProgressBar(pb, i)
  rm(individus_simple_reste_2)
  # On appelle le garbage collector toutes les 10 itérations
  if(i %% 10 == 0) gc()
}

dept_n_indiv <- 
  dept_n_indiv %>% 
  rbindlist()

Reste alors à mettre à jour individus_simple pour répercuter les informations obtenues sur le département de naissance des individus.

dept_n_indiv <- dept_n_indiv %>% rename(dept_N_2 = dept_N)

setkey(dept_n_indiv, id_personne_num)
setkey(individus_simple, id_personne_num)

individus_simple <- dept_n_indiv[individus_simple]
head(individus_simple)

individus_simple <- 
  individus_simple %>% 
  mutate(dept_N = ifelse(!is.na(dept_N_2), yes = dept_N_2, no = dept_N))

individus_simple <- data.table(individus_simple)

1.5.5 Rajout de l’âge des individus

Nous pouvons recalculer l’âge des individus.

individus_simple[, date_N_annee := as.numeric(str_sub(date_N, 1, 4))]
individus_simple[, date_D_annee := as.numeric(str_sub(date_D, 1, 4))]
individus_simple[, age := as.numeric(date_D_annee) - as.numeric(date_N_annee)]

1.5.6 Toutes les informations sur une même ligne

Pour notre tableau de génération, chaque ligne doit permettre de suivre un descendant. Chaque ligne doit contenir les informations d’une personne et de ses aïeux jusqu’à celui né entre 1800 et 1804.

var_recup <- c("id_personne_num", "long_N", "lat_N", "sexe", "date_N", "date_D", "date_M", "dept_N", "date_N_annee", "date_D_annee", "age")

infos_aieul <- 
  individus_simple %>% 
  s_select(var_recup) %>% 
  magrittr::set_colnames(str_c(var_recup, "_aieul")) %>% 
  rename(id_aieul = id_personne_num_aieul)

infos_enfants <- 
  individus_simple %>% 
  s_select(var_recup) %>% 
  magrittr::set_colnames(str_c(var_recup, "_enfants")) %>% 
  rename(id_enfant = id_personne_num_enfants)

infos_p_enfant <- 
  individus_simple %>% 
  s_select(var_recup) %>% 
  magrittr::set_colnames(str_c(var_recup, "_p_enfants")) %>% 
  rename(id_p_enfant = id_personne_num_p_enfants)

infos_ap_enfant <- 
  individus_simple %>% 
  s_select(var_recup) %>% 
  magrittr::set_colnames(str_c(var_recup, "_ap_enfants")) %>% 
  rename(id_ap_enfant = id_personne_num_ap_enfants)

setkey(infos_aieul, id_aieul)
setkey(infos_enfants, id_enfant)
setkey(infos_p_enfant, id_p_enfant)
setkey(infos_ap_enfant, id_ap_enfant)


setkey(generations, id_aieul)
generations <- infos_aieul[generations]

setkey(generations, id_enfant)
generations <- infos_enfants[generations]

setkey(generations, id_p_enfant)
generations <- infos_p_enfant[generations]

setkey(generations, id_ap_enfant)
generations <- infos_ap_enfant[generations]

rm(infos_aieul, infos_enfants, infos_p_enfant, infos_ap_enfant)

generations

1.5.7 Distances avec l’aieul

En fonction des coordonnées GPS fournies, nous sommes à même de calculer la distance séparant un enfant du lieu de naissance de ses aïeux. Passons d’abord toutes les variables de coordonnées en valeurs numériquues.

library(geosphere)
generations <- 
  generations %>% 
  mutate(long_N_aieul = as.numeric(long_N_aieul),
         lat_N_aieul = as.numeric(lat_N_aieul),
         long_N_enfants = as.numeric(long_N_enfants),
         lat_N_enfants = as.numeric(lat_N_enfants),
         long_N_p_enfants = as.numeric(long_N_p_enfants),
         lat_N_p_enfants = as.numeric(lat_N_p_enfants),
         long_N_ap_enfants = as.numeric(long_N_ap_enfants),
         lat_N_ap_enfants = as.numeric(lat_N_ap_enfants))

generations <- data.table(generations)

Nous allons utiliser la fonction distm() du package geosphere pour calculer les distances (en mètres) entre deux points.

calcul_distance <- function(long_1, lat_1, long_2, lat_2){
  distm(c(long_1, lat_1), c(long_2, lat_2), fun = distHaversine) %>% "["(1)
}

Pour chaque ligne du tableau de génération, nous calculons les trois distances suivantes :

1. dist_enfants : distance entre un enfant et l’aïeul de 1800 ;
2. dist_p_enfants : distance entre un petit-enfant et l’aïeul de 1800 ;
3. dist_ap_enfants : distance entre un arrière-petit-enfant et l’aïeul de 1800.

distances_l <- 
  pblapply(1:nrow(generations), function(y){
      x <- generations[y,]
      dist_enfants = calcul_distance(x$long_N_aieul, x$lat_N_aieul, x$long_N_enfants, x$lat_N_enfants)
      dist_p_enfants = calcul_distance(x$long_N_aieul, x$lat_N_aieul, x$long_N_p_enfants, x$lat_N_p_enfants)
      dist_ap_enfants = calcul_distance(x$long_N_aieul, x$lat_N_aieul, x$long_N_ap_enfants, x$lat_N_ap_enfants)
      data.frame(dist_enfants = dist_enfants,
                 dist_p_enfants = dist_p_enfants,
                 dist_ap_enfants = dist_ap_enfants)
    }) %>% 
    rbindlist() %>% 
    data.table()

Il reste à rajouter ces informations au tableau generation.

generations <- 
  generations %>% 
  cbind(distances_l)

Nous pouvons convertir ces distances en kilomètres.

generations <- 
  generations %>% 
  mutate(dist_enfants_km =  dist_enfants / 1000,
         dist_p_enfants_km = dist_p_enfants / 1000,
         dist_ap_enfants_km = dist_ap_enfants / 1000)

Effectuons quelques regroupements de départements, pour refléter l’état approximatif des départements au XIXe siècle.

liste_depts <- 
  liste_depts %>% 
  mutate(dept2 = dept,
         dept2 = ifelse(dept2 %in% c("F78", "F91", "F92", "F95"), yes = "F78", no = dept2),
         dept2 = ifelse(dept2 %in% c("F93", "F94", "F77"), yes = "F77", no = dept2)
  )

Faisons de même sur la base des générations.

generations <- 
  generations %>% 
  mutate(dept_N_aieul = ifelse(dept_N_aieul %in% c("F78", "F91", "F92", "F95"), yes = "F78", no = dept_N_aieul),
         dept_N_enfants = ifelse(dept_N_enfants %in% c("F78", "F91", "F92", "F95"), yes = "F78", no = dept_N_enfants),
         dept_N_p_enfants = ifelse(dept_N_p_enfants %in% c("F78", "F91", "F92", "F95"), yes = "F78", no = dept_N_p_enfants),
         dept_N_ap_enfants = ifelse(dept_N_ap_enfants %in% c("F78", "F91", "F92", "F95"), yes = "F78", no = dept_N_ap_enfants),
         #
         dept_N_aieul = ifelse(dept_N_aieul %in% c("F93", "F94", "F77"), yes = "F77", no = dept_N_aieul),
         dept_N_enfants = ifelse(dept_N_enfants %in% c("F93", "F94", "F77"), yes = "F77", no = dept_N_enfants),
         dept_N_p_enfants = ifelse(dept_N_p_enfants %in% c("F93", "F94", "F77"), yes = "F77", no = dept_N_p_enfants),
         dept_N_ap_enfants = ifelse(dept_N_ap_enfants %in% c("F93", "F94", "F77"), yes = "F77", no = dept_N_ap_enfants)
  )

generations <- data.table(generations)

1.6 Adelphies

Nous pouvons extraire des informations relatives aux adelphies (fratries et sororités).

Le nombre d’enfants par femme de la génération 1 (enfants).

nb_enf_aieux <- 
  generations %>% 
  filter(sexe_aieul == 2) %>% 
  mutate(keep = TRUE) %>% 
  mutate(keep = ifelse(age_aieul < 13, yes = FALSE, no = keep)) %>% 
  mutate(keep = ifelse(is.na(dept_M_aieul), yes = FALSE, no = keep)) %>% 
  filter(keep) %>% 
  select(id_aieul, id_enfant, dept_N_aieul) %>% 
  unique() %>% 
  mutate(nb_enf = ifelse(is.na(id_enfant), yes = 0, no = 1)) %>% 
  group_by(id_aieul, dept_N_aieul) %>% 
  summarise(nb_enf = sum(nb_enf))

summary(nb_enf_aieux$nb_enf)

Le nombre d’enfants par femme de la generation 2 (petits enfants).

nb_enf_enfants <- 
  generations %>% 
  filter(sexe_enfants == 2) %>% 
  mutate(keep = TRUE) %>% 
  mutate(keep = ifelse(age_enfants < 13, yes = FALSE, no = keep)) %>% 
  mutate(keep = ifelse(is.na(dept_M_enfant), yes = FALSE, no = keep)) %>% 
  filter(keep) %>% 
  select(id_enfant, id_p_enfant, dept_N_p_enfants) %>% 
  unique() %>% 
  mutate(nb_enf = ifelse(is.na(id_p_enfant), yes = 0, no = 1)) %>% 
  group_by(id_enfant, dept_N_p_enfants) %>% 
  summarise(nb_enf = sum(nb_enf))

summary(nb_enf_enfants$nb_enf)

Le nombre d’enfants par femme de la generation 3 (arrieres petits enfants).

nb_enf_p_enfants <- 
  generations %>% 
  filter(sexe_p_enfants == 2) %>% 
  mutate(keep = TRUE) %>% 
  mutate(keep = ifelse(age_p_enfants < 13, yes = FALSE, no = keep)) %>% 
  mutate(keep = ifelse(is.na(dept_M_p_enfant), yes = FALSE, no = keep)) %>% 
  filter(keep) %>% 
  select(id_p_enfant, id_ap_enfant, dept_N_ap_enfants) %>% 
  unique() %>% 
  mutate(nb_enf = ifelse(is.na(id_ap_enfant), yes = 0, no = 1)) %>% 
  group_by(id_p_enfant, dept_N_ap_enfants) %>% 
  summarise(nb_enf = sum(nb_enf))

summary(nb_enf_p_enfants$nb_enf)

1.7 Base des individus mobiles

Nous nous intéressons dans l’étude à la mobilité des individus, entre les générations. Il faut de fait se concentrer sur des individus en capacité de se déplacer. Nous plaçons une limite inférieure arbitraire sur les âges pour nos individus, de 16 ans. Nous émettons l’hypothèse qu’avant 16 ans, les individus ne changeront pas de village.

generations_mobilite <- 
  generations %>% 
  filter(age_aieul >= 16,
         age_enfants >= 16,
         age_p_enfants >= 16,
         age_ap_enfants >= 16) %>% 
  data.table()

Enfin, nous sauvegardons les deux bases de données generations et generations_mobilite.

save(generations, file = "generations_migration.rda")
save(generations_mobilite, file = "generations_mobilite.rda")

2 Représentativité des données

Maintenant, nous disposons d’une base nettoyée, bien qu’elle contienne encore sûrement des doublons que les exécutions de nos algorithmes n’ont pas su déceler. Nous pouvons alors l’explorer. Nous nous intéressons dans cette section à quelques comparaisons de statistiques descriptives des données de généalogie avec certains résultats de la littérature.

2.1 Naissances

Un premier élément permettant de caractériser notre échantillon est sa taille. Nous regardons le nombre de naissances par année en différentiant les générations. Puis, nous comparons le nombre de naissances de l’échantillon initial avec des données officielles, pour avoir une idée de la représentativité de la génération à l’origine de notre étude.

2.1.1 Nombre de naissances par génération

Nous chargeons quelques packages, ainsi que les données de généalogie.

library(tidyverse)
library(stringr)
library(data.table)
library(pbapply)
library(dplyrExtras)
library(sp)

load("generations_migration.rda")

Pour chacune des générations, nous calculons le nombre de naissance par année.

nb_nais_aieux <- 
  generations %>% 
  tbl_df() %>% 
  select(id_aieul, date_N_aieul) %>% 
  unique() %>% 
  mutate(date_N_aieul_annee = str_sub(date_N_aieul, 1, 4) %>% as.numeric()) %>% 
  group_by(date_N_aieul_annee) %>% 
  tally() %>% 
  arrange(date_N_aieul_annee)

nb_nais_enfants <- 
  generations %>% 
  tbl_df() %>% 
  select(id_enfant, date_N_enfants) %>% 
  unique() %>% 
  mutate(date_N_enfants_annee = str_sub(date_N_enfants, 1, 4) %>% as.numeric()) %>% 
  group_by(date_N_enfants_annee) %>% 
  tally() %>% 
  arrange(date_N_enfants_annee)

nb_nais_p_enfants <- 
  generations %>% 
  tbl_df() %>% 
  select(id_p_enfant, date_N_p_enfants) %>% 
  unique() %>% 
  mutate(date_N_p_enfants_annee = str_sub(date_N_p_enfants, 1, 4) %>% as.numeric()) %>% 
  group_by(date_N_p_enfants_annee) %>% 
  tally() %>% 
  arrange(date_N_p_enfants_annee)

nb_nais_ap_enfants <- 
  generations %>% 
  tbl_df() %>% 
  select(id_ap_enfant, date_N_ap_enfants) %>% 
  unique() %>% 
  mutate(date_N_ap_enfants_annee = str_sub(date_N_ap_enfants, 1, 4) %>% as.numeric()) %>% 
  group_by(date_N_ap_enfants_annee) %>% 
  tally() %>% 
  arrange(date_N_ap_enfants_annee)

Nous mettons en commun ces observations.

nb_nais <- 
  nb_nais_aieux %>% 
  rename(annee_N = date_N_aieul_annee) %>% 
  mutate(generation = 0) %>% 
  bind_rows(nb_nais_enfants %>% 
              rename(annee_N = date_N_enfants_annee) %>% 
              mutate(generation = 1)
  ) %>% 
  bind_rows(
    nb_nais_p_enfants %>% 
      rename(annee_N = date_N_p_enfants_annee) %>% 
      mutate(generation = 2)
  ) %>% 
  bind_rows(
    
    nb_nais_ap_enfants %>% 
      rename(annee_N = date_N_ap_enfants_annee) %>% 
      mutate(generation = 3)
  ) %>% 
  mutate(generation = factor(generation, levels = seq(0,3), labels = c('Aïeux', "Enfants",
                                                                       "Petits-enfants", "Arrière-petits-enfants")))

Puis, nous pouvons tracer un histogramme, avec une échelle log pour les ordonnées, le nombre de naissance d’aïeux étant relativement explosif vis-à-vis du nombre de naissances annuelles pour les descendants.

p_distrib_naissances_gen <- 
  ggplot(data = nb_nais, aes(x = annee_N, y = log(n))) +
  geom_bar(stat="identity", position = "identity", aes(fill =generation , group = generation), alpha = .5) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10), limits = c(1800, 1940)) +
  xlab("Nombre d'individus") + ylab("Année de naissance") +
  scale_fill_discrete("")

p_distrib_naissances_gen

2.1.2 Comparaison des nombres de naissances avec les données de l’INSEE

Pour avoir une idée de la représentativité des données, nous regardons, par département, combien de naissances sont renseignées dans nos données et comparons le résultat avec les enregistrements de la Statistique Générale de la France (SGF) proposés sur le site de l’INSEE.

Nous chargeons nos données ainsi que celles de la SGF Nous allons nous intéresser seulement à quelques variables parmi les 228 que propose le fichier de la SGF.

load("generations_migration.rda")
pop_insee <- readxl::read_excel("../data/population_insee/TERR_T86.xls", skip = 8, col_names = FALSE)
pop_insee_col <- readxl::read_excel("../data/population_insee/TERR_T86.xls", skip = 6) %>% 
colnames()
colnames(pop_insee) <- pop_insee_col
# V3 : Numéro de département
# V6 : Nom de l'unité d'analyse
# V221 : Naissances totales; légitimes et naturels; garçons, 1800 à 1801
# V222 : Naissances totales; légitimes et naturels: filles 1800 à 1801

Une petite correction de numéro de département, pour pouvoir apparier nos données à celles de la SGF.

pop_insee <- 
pop_insee %>% 
mutate(V3 = ifelse(V3 == "99", yes = "54", no = V3))

Nous retenons uniquement les variables qui nous intéressent.

nb_nais_region <- 
pop_insee %>% 
select(V3, V221, V222) %>% 
rename(num_dep = V3, nb_nais_f = V222, nb_nais_h = V221) %>% 
filter(num_dep != "00")

head(nb_nais_region)

# A tibble: 6 x 3
  num_dep nb_nais_h nb_nais_f
  <chr>       <dbl>     <dbl>
1 01           4892      4811
2 02           6801      6534
3 03           4108      3883
4 04           2813      2739
5 05           2070      1945
6 07           4496      4288

Nous sommons les naissances pour les femmes et les hommes.

nb_nais_1801 <- 
nb_nais_region %>% 
summarise(nb_nais_h = sum(nb_nais_h, na.rm=T), nb_nais_f = sum(nb_nais_f, na.rm=T))

head(nb_nais_1801)

# A tibble: 1 x 2
  nb_nais_h nb_nais_f
      <dbl>     <dbl>
1    464562    439126

Nous calculons à présent le nombre de naissances de filles et de garçons pour chaque génération, pour chaque année, dans les données de généalogie.

nb_naiss_par_annee <- 
generations %>% 
tbl_df() %>% 
select(id_aieul, sexe_aieul, date_N_annee_aieul) %>% 
unique() %>% 
rename(id_personne = id_aieul, sexe = sexe_aieul, date_N_annee = date_N_annee_aieul) %>% 
bind_rows(
generations %>% 
tbl_df() %>% 
select(id_enfant, sexe_enfants, date_N_annee_enfants) %>% 
unique() %>% 
rename(id_personne = id_enfant, sexe = sexe_enfants, date_N_annee = date_N_annee_enfants)
) %>% 
bind_rows(
generations %>% 
tbl_df() %>% 
select(id_p_enfant, sexe_p_enfants, date_N_annee_p_enfants) %>% 
unique() %>% 
rename(id_personne = id_p_enfant, sexe = sexe_p_enfants, date_N_annee = date_N_annee_p_enfants)
) %>% 
bind_rows(
generations %>% 
tbl_df() %>% 
select(id_ap_enfant, sexe_ap_enfants, date_N_annee_ap_enfants) %>% 
unique() %>% 
rename(id_personne = id_ap_enfant, sexe = sexe_ap_enfants, date_N_annee = date_N_annee_ap_enfants)
) %>% 
unique() %>% 
group_by(sexe, date_N_annee) %>% 
tally()

L’idée est ensuite de construire un tableau (representativite) uniquement pour les aïeux nés en 1801, pour comparer avec les valeurs de la SGF. Le tableau que nous construisons indique, pour chaque département, le nombre de naissances enregistrées dans les données de généalogie, pour les femmes et pour les hommes. Nous rajoutons ensuite les données de la SGF contenues dans le tableau nb_nais_region. Il s’agit ensuite de comparer les valeurs, en faisant un ratio du nombre de naissances dans les données de généalogies sur le nombre de naissances reporté par la SGF.

representativite <- 
generations %>% 
tbl_df() %>% 
filter(date_N_annee_aieul == 1801) %>% 
select(id_aieul, sexe_aieul, date_N_annee_aieul, dept_N_aieul) %>% 
unique() %>% 
rename(dept_N = dept_N_aieul,
date_N_annee = date_N_annee_aieul,
sexe = sexe_aieul) %>% 
group_by(sexe, date_N_annee, dept_N) %>% 
summarise(nb_nais = n()) %>% 
left_join(
nb_nais_region %>% gather(sexe, nb_nais_1801, -num_dep) %>% 
mutate(sexe = ifelse(sexe == "nb_nais_h", yes = "1", no = "2")) %>% 
mutate(dept_N = str_c("F", num_dep))
) %>% 
filter(!is.na(sexe)) %>% 
mutate(pct_nais = nb_nais / nb_nais_1801 * 100,
pct_nais = round(pct_nais, 2)) %>% 
ungroup() %>% 
mutate(sexe = factor(sexe, levels = c(2,1), labels = c("Femmes", "Hommes"))) %>% 
arrange(date_N_annee, sexe)

head(representativite)

# A tibble: 6 x 7
  sexe   date_N_annee dept_N nb_nais num_dep nb_nais_1801 pct_nais
  <fct>         <dbl> <chr>    <int> <chr>          <dbl>    <dbl>
1 Femmes         1801 F01       1332 01              4811     27.7
2 Femmes         1801 F02       2341 02              6534     35.8
3 Femmes         1801 F03       1659 03              3883     42.7
4 Femmes         1801 F04        668 04              2739     24.4
5 Femmes         1801 F05        842 05              1945     43.3
6 Femmes         1801 F06        779 <NA>              NA     NA  

Nous retirons les départements pour lesquels la SGF ne communique pas de données.

representativite <- 
representativite %>% 
filter(!is.na(dept_N))

Nous extrayons les numéros du département.

representativite <- 
representativite %>% 
select(dept_N, sexe, nb_nais, nb_nais_1801, pct_nais) %>% 
mutate(num_dep = str_sub(dept_N, 2))

head(representativite)

# A tibble: 6 x 6
  dept_N sexe   nb_nais nb_nais_1801 pct_nais num_dep
  <chr>  <fct>    <dbl>        <dbl>    <dbl> <chr>  
1 F01    Femmes    1332         4811     27.7 01     
2 F01    Hommes    1474         4892     30.1 01     
3 F02    Femmes    2341         6534     35.8 02     
4 F02    Hommes    2836         6801     41.7 02     
5 F03    Femmes    1659         3883     42.7 03     
6 F03    Hommes    1866         4108     45.4 03     

Nous pouvons ensuite proposer une représentation graphique sous forme de carte choroplèthe la représentativité régionale de notre échantillon d’aïeux. Pour ce faire, nous chargons les données permettant de tracer la France. Plus de détails sur l’obtention de ce fichier sont donnés en annexe.

load("france_1800.rda")

representativite_map <- 
france_1800 %>% 
left_join(
representativite %>% 
filter(!is.na(sexe)) %>% 
select(num_dep, sexe, pct_nais) %>% 
spread(sexe, pct_nais)
)

p_repres <- ggplot() +
geom_polygon(data = representativite_map %>% gather(sexe, pct_nais, Femmes, Hommes),
aes(x = long, y = lat, group = group, fill = pct_nais), size = .4) +
coord_quickmap() +
scale_fill_gradientn("%", colours = rev(heat.colors(10)), na.value = "grey") +
facet_wrap(~sexe)

p_repres

2.2 Ratio Femmes/Hommes

On peut également calculer le pourcentages d’hommes à partir des données de généalogie de la manière suivante : on créé un tableau avec tous les individus uniques de l’échantillon, puis ont somme le nombre d’individus par sexe.

sexes <- generations %>% 
  select(id_aieul, sexe_aieul, dept_N_aieul) %>% 
  unique() %>% 
  rename(id_personne = id_aieul, sexe = sexe_aieul, dept_N = dept_N_aieul) %>% 
  mutate(sexe = as.character(sexe)) %>% 
  bind_rows(
    generations %>% 
      select(id_enfant, sexe_enfants, dept_N_enfants) %>% 
      unique() %>% 
      rename(id_personne = id_enfant, sexe = sexe_enfants, dept_N = dept_N_enfants) %>% 
      mutate(sexe = as.character(sexe))
  ) %>% 
  bind_rows(
    generations %>% 
      select(id_p_enfant, sexe_p_enfants, dept_N_p_enfants) %>% 
      unique() %>% 
      rename(id_personne = id_p_enfant, sexe = sexe_p_enfants, dept_N = dept_N_p_enfants) %>% 
      mutate(sexe = as.character(sexe))
  ) %>% 
  bind_rows(
    generations %>% 
      select(id_ap_enfant, sexe_ap_enfants,dept_N_ap_enfants ) %>% 
      unique() %>% 
      rename(id_personne = id_ap_enfant, sexe = sexe_ap_enfants, dept_N = dept_N_ap_enfants) %>% 
      mutate(sexe = as.character(sexe))
  ) %>% 
  unique()

sexes %>% 
  tbl_df() %>% 
  group_by(sexe) %>% 
  tally() %>% 
  ungroup() %>% 
  filter(!is.na(sexe)) %>% 
  mutate(pct = n/sum(n)) %>% 
  head()

# A tibble: 2 x 3
  sexe        n   pct
  <chr>   <int> <dbl>
1 1     1316749 0.538
2 2     1132795 0.462

On peut également calculer le taux de masculinité :

r sexes %>% tbl_df() %>% group_by(sexe) %>% tally() %>% ungroup() %>% filter(!is.na(sexe)) %>% spread(sexe, n) %>% mutate(sex_ratio = `1` / `2` * 100)

# A tibble: 1 x 3 `1` `2` sex_ratio <int> <int> <dbl> 1 1316749 1132795 116

Il n’est pas compliqué de calculer le pourcentage d’hommes par départements :

```r representativite_fh <- sexes %>% filter(str_detect(dept_N, “^F[[:digit:]]{2}”)) %>% filter(!is.na(sexe)) %>% group_by(dept_N, sexe) %>% summarise(n = n()) %>% group_by(dept_N) %>% mutate(pct = n / sum(n)) %>% rename(departement = dept_N)

representativite_fh_sex_ratio <- sexes %>% filter(str_detect(dept_N, “^F[[:digit:]]{2}”)) %>% filter(!is.na(sexe)) %>% group_by(dept_N, sexe) %>% summarise(n = n()) %>% spread(sexe, n) %>% mutate(sex_ratio = 1 / 2 * 100) ```

On peut alors préparer un tableau de données pour effectuer une représentation graphique :

r representativite_map_fh <- france_1800 %>% left_join( representativite_fh %>% filter(sexe == 1) )

p_repres_fh <- ggplot() +
  geom_polygon(data = representativite_map_fh,
               aes(x = long, y = lat, group = group, fill = pct), col = "black", size = .4) +
  coord_quickmap() +
  scale_fill_gradientn("%", colours = rev(heat.colors(10)), na.value = "grey")
p_repres_fh

2.3 Taux de fécondité

Nous allons nous intéresser au taux de fécondité. Aussi, nous devons calculer le nombre d’enfants par femme. Nous le faisons pour chacune des générations.

# Nombre d'enfants par femme : generation 1 (enfants)
nb_enf_aieux <- 
  generations %>% 
  tbl_df() %>% 
  filter(sexe_aieul == 2) %>% 
  mutate(keep = TRUE) %>% 
  mutate(keep = ifelse(age_aieul < 15, yes = FALSE, no = keep)) %>% 
  mutate(keep = ifelse(is.na(dept_M_aieul), yes = FALSE, no = keep)) %>% 
  filter(keep) %>% 
  select(id_aieul, id_enfant, dept_N_aieul) %>% 
  unique() %>% 
  mutate(nb_enf = ifelse(is.na(id_enfant), yes = 0, no = 1)) %>% 
  group_by(id_aieul, dept_N_aieul) %>% 
  summarise(nb_enf = sum(nb_enf))

# Nombre d'enfants par femme : generation 2 (petits enfants)
nb_enf_enfants <- 
  generations %>% 
  filter(sexe_enfants == 2) %>% 
  mutate(keep = TRUE) %>% 
  mutate(keep = ifelse(age_enfants < 15, yes = FALSE, no = keep)) %>% 
  mutate(keep = ifelse(is.na(dept_M_enfants), yes = FALSE, no = keep)) %>% 
  filter(keep) %>% 
  select(id_enfant, id_p_enfant, dept_N_p_enfants) %>% 
  unique() %>% 
  mutate(nb_enf = ifelse(is.na(id_p_enfant), yes = 0, no = 1)) %>% 
  group_by(id_enfant, dept_N_p_enfants) %>% 
  summarise(nb_enf = sum(nb_enf))

# Nombre d'enfants par femme : generation 3 (arrieres petits enfants)
nb_enf_p_enfants <- 
  generations %>% 
  filter(sexe_p_enfants == 2) %>% 
  mutate(keep = TRUE) %>% 
  mutate(keep = ifelse(age_p_enfants < 15, yes = FALSE, no = keep)) %>% 
  mutate(keep = ifelse(is.na(dept_M_p_enfants), yes = FALSE, no = keep)) %>% 
  filter(keep) %>% 
  select(id_p_enfant, id_ap_enfant, dept_N_ap_enfants) %>% 
  unique() %>% 
  mutate(nb_enf = ifelse(is.na(id_ap_enfant), yes = 0, no = 1)) %>% 
  group_by(id_p_enfant, dept_N_ap_enfants) %>% 
  summarise(nb_enf = sum(nb_enf))

summary(nb_enf_aieux$nb_enf)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.000   1.357   2.000  33.000 

summary(nb_enf_enfants$nb_enf)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.000   1.605   2.000  19.000 

summary(nb_enf_p_enfants$nb_enf)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.000   1.615   2.000  21.000 

Ensuite, nous rejoignons toutes ces valeurs dans un seul tableau de données.

meres <- 
  nb_enf_aieux %>% 
  ungroup() %>% 
  rename(id_personne_num = id_aieul, dept_N = dept_N_aieul) %>% 
  mutate(generation = 0) %>% 
  bind_rows(
    nb_enf_enfants %>% 
      ungroup() %>% 
      rename(id_personne_num = id_enfant, dept_N = dept_N_p_enfants) %>% 
      mutate(generation = 1)
  ) %>% 
  bind_rows(
    nb_enf_p_enfants %>% 
      ungroup() %>% 
      rename(id_personne_num = id_p_enfant, dept_N = dept_N_ap_enfants) %>% 
      mutate(generation = 2)
  ) %>% 
  tbl_df()

Le nombre d’enfants moyen, sans se soucier des générations :

mean(meres$nb_enf)

[1] 1.467611

La répartition du nombre d’enfants, par génération (en lignes) est donnée dans le tableau ci-après. Les colonnes indiquent le nombre d’enfants. La colonne intitulée 11 concerne le pourcentage de femmes de l’échantillon a avoir donné naissance à 11 enfants et plus.

meres %>% 
  mutate(nb_enf = ifelse(nb_enf > 10, yes = 11, no = nb_enf)) %>% 
  select(generation, nb_enf) %>% 
  group_by(generation, nb_enf) %>% 
  summarise(nn = n()) %>% 
  ungroup () %>% 
  group_by(generation) %>% 
  mutate(pct_nn = nn / sum(nn)*100) %>% 
  select(generation, nb_enf, pct_nn) %>% 
  spread(nb_enf, pct_nn) %>% 
  knitr::kable(digits=2)

generation	0	1	2	3	4	5	6	7	8	9	10	11
0	46.98	24.54	10.26	6.14	4.02	2.62	1.89	1.20	0.85	0.53	0.36	0.63
1	27.36	39.18	13.18	7.48	4.80	2.86	1.87	1.22	0.84	0.51	0.27	0.43
2	25.81	40.33	13.72	7.88	4.50	2.68	1.89	1.17	0.78	0.47	0.30	0.46

3 Dimension temporelle des données : mortalité

Dans cette section, nous présentons les codes développés pour étudier les aspects temporels relatifs aux données. Nous nous concentrons principalement sur les aïeux de la base, à savoir les individus nés entre 1800 et 1804.

Nous pouvons explorer les questions de mortalité grâce aux dates de naissance et de décès des individus. Nous allons principalement nous concentrer sur les aïeux, mais les méthodes peuvent être appliquées à l’ensemble de la base.

aieux <-
  generations %>%
  select(id_aieul, date_N_annee_aieul, date_D_annee_aieul, age_aieul, sexe_aieul, dept_N_aieul) %>%
  unique() %>%
  filter(!is.na(date_N_annee_aieul), !is.na(date_D_annee_aieul)) %>%
  mutate(end_age = date_D_annee_aieul - date_N_annee_aieul) %>%
  tbl_df()

Nous avons besoin de charger quelques packages.

library(tidyverse)
library(demography)
library(survival)
library(data.table)

3.1 Estimation de la survie des individus

3.1.1 Estimations à partir des données issues de tables de mortalité

Nous avons également besoin de récupérer des données de mortalité. Nous téléchargeons les tables de mortalités de Vallin et al. (2001).

lifetable_fr <- read_csv("lifetable_FRA.csv")

Nous définissons une fonction pour extraire, pour une cohorte donnée, la probabilité de décès instantanné \(\mu_x\) (qui est équivalente à \(q_x\), c’est à dire la probabilité qu’un individu âgé de \(x\) ne survive pas en \(x+1\)).

#' qx_annee_age
#' Obtenir la proba qu'un individu âgé de @annee - @annee_naiss_depart
#' ne survive pas à l'âge de x
#' @annee (int) : Année pour laquelle extraire la table de mortalité
#' @annee_naiss_depart (int) : Annee de naissance des individus que l'on suit
# annee = 1806 ; annee_naiss_depart <- 1804
qx_annee_age <- function(annee, annee_naiss_depart, sexe){
  # Âge des individus que l'on suit
  age_courant <- annee - annee_naiss_depart
  # La table de mortalite en France pour l'année annee et le sexe sexe
  table_mortalite_annee <- lifetable_fr %>% filter(Year1 == annee, Sex == sexe, TypeLT == 1)
  # Extraire la ligne de cette table pour les individus que l'on suit
  res <- 
    table_mortalite_annee %>% 
    filter(Age == age_courant) %>% 
    slice(1) %>% 
    select(`q(x)`, `m(x)`) %>% 
    rename(qx = `q(x)`, mx = `m(x)`)
    # magrittr::extract2("q(x)")
  if(nrow(res) == 0) res <- data.frame(qx = NA, mx = NA)
  res
}# Fin de qx_annee_age()

Nous définissons une autre fonction pour calculer la probabilité qu’un individu âgé de \(x\) atteigne \(x+1\), soit \(p_x\), et la probabilité de survie conditionnelle à l’âge \(S_x(t)\). Ces valeurs sont estimées par cohorte.

#' survie_diag
#' On parcourt chaque année jusqu'à annee_naiss_depart+105
#' Mais on commence à 1806 (pas de données avant)
#' annee_naiss_depart <- 1804 ; sexe <- 1 ; age_min = 6
survie_diag <- function(annee_naiss_depart, sexe = c(1,2), age_min = 6){
  nqx <- lapply(seq(1806, annee_naiss_depart+105), qx_annee_age, annee_naiss_depart = annee_naiss_depart, sexe = sexe) %>% 
    bind_rows()

  ages <- seq(1806, annee_naiss_depart+105) - annee_naiss_depart
  res <- nqx
  res$age <- ages
  res$sexe <- as.character(sexe)
  
  res <- res %>% filter(age >= age_min)
  # Calcul de p_x (proba que (x) atteigne x+1)
  res$px = 1-res$qx
  # Rajout de la survie
  res$sx <- cumprod(res$px)
  res$cohorte <- annee_naiss_depart

  res
}# Fin de survie_diag

On applique la fonction survie_diag() aux années de naissances des individus des cohortes de 1800 à 1804, pour les femmes, puis pour les hommes.

survies_diagonale_1800_1804_f <- 
  pblapply(seq(1800, 1804), survie_diag, sexe = "2") %>% 
  bind_rows()


survies_diagonale_1800_1804_h <- 
  pblapply(seq(1800, 1804), survie_diag, sexe = "1") %>% 
  bind_rows()


survies_diagonale_1800_1804 <- 
  survies_diagonale_1800_1804_f %>% 
  bind_rows(survies_diagonale_1800_1804_h)

head(survies_diagonale_1800_1804)

# A tibble: 6 x 7
       qx      mx   age sexe     px    sx cohorte
    <dbl>   <dbl> <int> <chr> <dbl> <dbl>   <int>
1 0.0138  0.0139      6 2     0.986 0.986    1800
2 0.0118  0.0119      7 2     0.988 0.975    1800
3 0.0101  0.0102      8 2     0.990 0.965    1800
4 0.00808 0.00811     9 2     0.992 0.957    1800
5 0.00666 0.00668    10 2     0.993 0.950    1800
6 0.00591 0.00593    11 2     0.994 0.945    1800

survies_diagonale_1800_1804 %>% 
  mutate(sexe = factor(sexe, levels = c(2,1), labels = c("Femmes", "Hommes"))) %>% 
  ggplot(data = .,
       aes(x = age, y = sx, colour = factor(cohorte), linetype = sexe)) +
  geom_line() +
  ggtitle("Fonctions de survie des femmes et des hommes par cohortes")

Les estimations sont données par cohorte, nous les aggrégeons pour le groupe 1800-1804, en prenant la moyenne des probabilités de survie et de décès instantanné, par âge et sexe.

surv_hist <- 
  survies_diagonale_1800_1804 %>%
  group_by(age, sexe) %>%
  summarise(sx = mean(sx, na.rm=T),
            mx = mean(mx))
head(surv_hist)

# A tibble: 6 x 4
# Groups: age [3]
    age sexe     sx      mx
  <int> <chr> <dbl>   <dbl>
1     6 1     0.986 0.0137 
2     6 2     0.987 0.0135 
3     7 1     0.975 0.0112 
4     7 2     0.976 0.0112 
5     8 1     0.966 0.0102 
6     8 2     0.966 0.00983

3.1.2 Estimations à partir des données de généalogie

Maintenant que nous disposons d’estimations de la probabilité de survie à chaque âge à partir des données issues de tables de mortalité, nous proposons nos estimations réalisées à partir des données de généalogie de Geneanet. Nous définissons une fonction, calcul_survie_geneanet(), qui à partir d’une année de naissance, estime les probabilités de survie, de décès instantanné pour les hommes et les femmes, pour une la cohorte des individus nés cette année là.

#' calcul_survie_geneanet
#' Estime les probabilités de survie, de décès instantanné (en log) pour les hommes et les femmes
#' pour des individus d'une cohorte donnée
#' @une_annee (int) : année de naissance des individus à suivre
calcul_survie_geneanet <- function(une_annee){
  aieux_annee <- aieux %>% 
    tbl_df() %>% 
    filter(date_N_annee_aieul == une_annee) %>% 
    filter(age_aieul >= 6)
  # On compte le nombre de personne encore en vie à chaque année
  calcul_exposition <- function(une_annee, sexe){
    aieux_annee %>% filter(date_D_annee_aieul >= une_annee, sexe_aieul == sexe) %>% nrow()
  }
  
  annees <- seq(une_annee+6, une_annee+110)
  
  exposition_annee_h <- data.frame(annee = annees, exposition = sapply(annees, calcul_exposition,  sexe = "1")) %>% tbl_df()
  exposition_annee_f <- data.frame(annee = annees, exposition = sapply(annees, calcul_exposition, sexe = "2")) %>% tbl_df()
  exposition_annee <- 
    exposition_annee_h %>% rename(exposition_h = exposition) %>% 
    left_join(exposition_annee_f %>% rename(exposition_f = exposition), by = "annee") %>% 
    mutate(age = annees-une_annee)
  
  
  exposition_annee <- 
    exposition_annee %>% 
    mutate(exposition_h = ifelse(exposition_h <= 0, yes = 1, no = exposition_h),
           exposition_f = ifelse(exposition_h <= 0, yes = 1, no = exposition_f)) %>% 
    # Calcul du nombre de décès
    # mutate(expo_moins_1_h = lag(exposition_h),
    #        expo_moins_1_f = lag(exposition_f)) %>% 
    mutate(expo_plus_1_h = lead(exposition_h),
           expo_plus_1_f = lead(exposition_f)) %>%
    mutate(deces_h = exposition_h - expo_plus_1_h,
           deces_f = exposition_f - expo_plus_1_f) %>% 
    mutate(sx_h = (exposition_h - deces_h) / exposition_h,
           sx_f = (exposition_f - deces_f) / exposition_f) %>% 
    select(-expo_plus_1_h, -expo_plus_1_f) %>% 
    # Proba de déceder dans l'année
    mutate(qxt_h = deces_h / exposition_h,
           qxt_f = deces_f / exposition_f) %>% 
    mutate(qxt_h = ifelse(qxt_h <= 0, yes = NA, no = qxt_h),
           qxt_f = ifelse(qxt_f <= 0, yes = NA, no = qxt_f)) %>% 
    # Calcul de la force de mortalité
    mutate(fm_h = log(qxt_h),
           fm_f = log(qxt_f)) %>% 
    mutate(annee_naissance = une_annee)
  
  # Survie
  # exposition_annee$lx_f[1] <- exposition_annee$lx_h[1] <- 1
  exposition_annee$sx_h <- cumprod(exposition_annee$sx_h)
  exposition_annee$sx_f <- cumprod(exposition_annee$sx_f)
  exposition_annee
}# Fin de calcul_survie_geneanet()

On calcule donc les probabilités de survie et de décès instantanné pour les individus des cohortes 1800 à 1804, à partir des données de généalogie :

library(pbapply)
surv_geneanet_tot <- 
  pblapply(1800:1804, calcul_survie_geneanet) %>% 
  bind_rows()

head(surv_geneanet_tot)

# A tibble: 6 x 13
  annee exposition_h exposition_f   age deces_h deces_f  sx_h  sx_f
  <int>        <int>        <int> <int>   <int>   <int> <dbl> <dbl>
1  1806        72587        61392     6     550     570 0.992 0.991
2  1807        72037        60822     7     492     427 0.986 0.984
3  1808        71545        60395     8     376     397 0.980 0.977
4  1809        71169        59998     9     307     277 0.976 0.973
5  1810        70862        59721    10     245     233 0.973 0.969
6  1811        70617        59488    11     272     263 0.969 0.965
# ... with 5 more variables: qxt_h <dbl>, qxt_f <dbl>, fm_h <dbl>, fm_f
#   <dbl>, annee_naissance <int>

Comme pour les données issues des tables de mortalité, nous calculons la moyenne par age et par sexe :

surv_geneanet <- 
  surv_geneanet_tot %>% 
  group_by(age) %>% 
  summarise(fm_h = mean(fm_h, na.rm=T), fm_f = mean(fm_f, na.rm=T),
            sx_h = mean(sx_h, na.rm=T), sx_f = mean(sx_f, na.rm=T))

head(surv_geneanet)

# A tibble: 6 x 5
    age  fm_h  fm_f  sx_h  sx_f
  <int> <dbl> <dbl> <dbl> <dbl>
1     6 -4.77 -4.63 0.992 0.990
2     7 -4.96 -4.82 0.985 0.982
3     8 -5.18 -5.07 0.979 0.976
4     9 -5.32 -5.22 0.974 0.971
5    10 -5.44 -5.32 0.970 0.966
6    11 -5.51 -5.46 0.966 0.962

Nous modifions la forme du dernier tableau pour obtenir une disposition propice à l’utilisation par les fonctions de graphiques de ggplot2.

surv_geneanet_sx <- 
  surv_geneanet %>% 
  select(age, sx_h, sx_f) %>% 
  gather(sexe, sx, sx_h, sx_f) %>% 
  mutate(sexe = factor(sexe, levels = c("sx_h", "sx_f"), labels = c(1,2))) %>% 
  mutate(sexe = as.character(sexe))

Nous mettons en commun les données historiques issues des tables avec celles de notre échantillon de données généalogiques.

df_plot_survie <- 
  surv_hist %>% mutate(type = "Historique") %>%
  bind_rows(
    surv_geneanet_sx %>% 
      mutate(type = "Geneanet")
  ) %>% 
  mutate(sexe = factor(sexe, levels = c("2","1"), labels = c("Femmes", "Hommes")))

Il est alors possible de représenter les fonctions de survie par sexe pour chacune des deux bases.

p <- 
  ggplot(data = df_plot_survie,
         aes(x = age, y = sx, colour = sexe, linetype = type)) +
  geom_line() +
  scale_colour_manual("", values = c("Femmes" = "red", "Hommes" = "blue")) +
  scale_linetype_manual("", values = c("Historique" = "dashed", "Geneanet" = "solid")) +
  xlab("âge") + ylab("Probabilité") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))
p

3.1.3 Estimations à partir des données de généalogie, pour l’ensemble des individus

Cette fois, on peut effectuer les calculs pour l’ensemble des individus, non plus en se restreignant à ceux ayant au moins 6 ans.

#' calcul_survie_geneanet_tot
#' Estime les probabilités de survie, de décès instantanné (en log) pour les hommes et les femmes
#' pour des individus d'une cohorte donnée, dans un département donné
#' @une_annee (int) : année de naissance des individus à suivre
# une_annee <- 1800
calcul_survie_geneanet_tot <- function(une_annee){
  aieux_annee <- 
    aieux %>% 
    tbl_df() %>% 
    filter(date_N_annee_aieul == une_annee)
  
  # On compte le nombre de personne encore en vie à chaque année
  calcul_exposition <- function(une_annee, sexe){
    aieux_annee %>% filter(date_D_annee_aieul >= une_annee, sexe_aieul == sexe) %>% nrow()
  }
  
  annees <- seq(une_annee, une_annee+110)
  
  exposition_annee_h <- data.frame(annee = annees, exposition = sapply(annees, calcul_exposition,  sexe = "1")) %>% tbl_df()
  exposition_annee_f <- data.frame(annee = annees, exposition = sapply(annees, calcul_exposition, sexe = "2")) %>% tbl_df()
  exposition_annee <- 
    exposition_annee_h %>% rename(exposition_h = exposition) %>% 
    left_join(exposition_annee_f %>% rename(exposition_f = exposition), by = "annee") %>% 
    mutate(age = annees-une_annee)
  
  
  exposition_annee <- 
    exposition_annee %>% 
    mutate(exposition_h = ifelse(exposition_h <= 0, yes = 1, no = exposition_h),
           exposition_f = ifelse(exposition_h <= 0, yes = 1, no = exposition_f)) %>% 
    # Calcul du nombre de décès
    # mutate(expo_moins_1_h = lag(exposition_h),
    #        expo_moins_1_f = lag(exposition_f)) %>% 
    mutate(expo_plus_1_h = lead(exposition_h),
           expo_plus_1_f = lead(exposition_f)) %>%
    mutate(deces_h = exposition_h - expo_plus_1_h,
           deces_f = exposition_f - expo_plus_1_f) %>% 
    mutate(sx_h = (exposition_h - deces_h) / exposition_h,
           sx_f = (exposition_f - deces_f) / exposition_f) %>% 
    select(-expo_plus_1_h, -expo_plus_1_f) %>% 
    # Proba de déceder dans l'année
    mutate(qxt_h = deces_h / exposition_h,
           qxt_f = deces_f / exposition_f) %>% 
    mutate(qxt_h = ifelse(qxt_h <= 0, yes = NA, no = qxt_h),
           qxt_f = ifelse(qxt_f <= 0, yes = NA, no = qxt_f)) %>% 
    # Calcul de la force de mortalité
    mutate(fm_h = log(qxt_h),
           fm_f = log(qxt_f)) %>% 
    mutate(annee_naissance = une_annee)
  
  # Survie
  # exposition_annee$lx_f[1] <- exposition_annee$lx_h[1] <- 1
  exposition_annee$sx_h <- cumprod(exposition_annee$sx_h)
  exposition_annee$sx_f <- cumprod(exposition_annee$sx_f)
  exposition_annee
}# Fin de calcul_survie_geneanet_tot()

surv_geneanet_2_tot_c <- pblapply(1800:1804, calcul_survie_geneanet_tot)
surv_geneanet_2_tot_c <- bind_rows(surv_geneanet_2_tot_c)

surv_geneanet_2_tot <- 
  surv_geneanet_2_tot_c %>% 
  group_by(age) %>% 
  summarise(fm_h = mean(fm_h, na.rm=T), fm_f = mean(fm_f, na.rm=T),
            sx_h = mean(sx_h, na.rm=T), sx_f = mean(sx_f, na.rm=T))

3.1.4 Estimations à partir des données de généalogie, par département

Nous voulons également explorer la mortalité au travers de ses différences régionales. Nous adaptons le code précédent pour que les estimations s’effectuent à l’échelle des régions.

#' calcul_survie_geneanet_2
#' Estime les probabilités de survie, de décès instantanné (en log) pour les hommes et les femmes
#' pour des individus d'une cohorte donnée, dans un département donné
#' @une_annee (int) : année de naissance des individus à suivre
#' @un_dep (string) : identifiant du département (lu dans le tableau `aieux`)
calcul_survie_geneanet_2 <- function(une_annee, un_dep){
  aieux_annee <- 
    aieux %>% 
    tbl_df() %>% 
    filter(date_N_annee_aieul == une_annee) %>% 
    filter(dept_N_aieul %in% un_dep)
  
  # On compte le nombre de personne encore en vie à chaque année
  calcul_exposition <- function(une_annee, sexe){
    aieux_annee %>% filter(date_D_annee_aieul >= une_annee, sexe_aieul == sexe) %>% nrow()
  }
  
  annees <- seq(une_annee, une_annee+110)
  
  exposition_annee_h <- data.frame(annee = annees, exposition = sapply(annees, calcul_exposition,  sexe = "1")) %>% tbl_df()
  exposition_annee_f <- data.frame(annee = annees, exposition = sapply(annees, calcul_exposition, sexe = "2")) %>% tbl_df()
  exposition_annee <- 
    exposition_annee_h %>% rename(exposition_h = exposition) %>% 
    left_join(exposition_annee_f %>% rename(exposition_f = exposition), by = "annee") %>% 
    mutate(age = annees-une_annee)
  
  
  exposition_annee <- 
    exposition_annee %>% 
    mutate(exposition_h = ifelse(exposition_h <= 0, yes = 1, no = exposition_h),
           exposition_f = ifelse(exposition_h <= 0, yes = 1, no = exposition_f)) %>% 
    # Calcul du nombre de décès
    mutate(expo_plus_1_h = lead(exposition_h),
           expo_plus_1_f = lead(exposition_f)) %>%
    mutate(deces_h = exposition_h - expo_plus_1_h,
           deces_f = exposition_f - expo_plus_1_f) %>% 
    mutate(sx_h = (exposition_h - deces_h) / exposition_h,
           sx_f = (exposition_f - deces_f) / exposition_f) %>% 
    select(-expo_plus_1_h, -expo_plus_1_f) %>% 
    # Proba de déceder dans l'année
    mutate(qxt_h = deces_h / exposition_h,
           qxt_f = deces_f / exposition_f) %>% 
    mutate(qxt_h = ifelse(qxt_h <= 0, yes = NA, no = qxt_h),
           qxt_f = ifelse(qxt_f <= 0, yes = NA, no = qxt_f)) %>% 
    # Calcul de la force de mortalité
    mutate(fm_h = log(qxt_h),
           fm_f = log(qxt_f)) %>% 
    mutate(annee_naissance = une_annee)
  
  # Survie
  exposition_annee$sx_h <- cumprod(exposition_annee$sx_h)
  exposition_annee$sx_f <- cumprod(exposition_annee$sx_f)
  exposition_annee <- 
    exposition_annee %>% 
    mutate(dept = un_dep)
  exposition_annee
}# Fin de calcul_survie_geneanet_2()

Nous allons parcourir chacun des départements disponibles. Nous dressons la liste de ces derniers.

les_depts <- unique(aieux$dept_N_aieul)
les_depts <- les_depts[!is.na(les_depts)]

Nous parcourons chacun de ces départements, et pour chaque département, nous parcourons chaque cohorte pour lui appliquer la fonction calcul_survie_geneanet_2().

surv_geneanet_2_depts_c <- 
  pblapply(les_depts, function(un_dep){
    lapply(1800:1804, calcul_survie_geneanet_2, un_dep = un_dep) %>% 
      bind_rows()
  })

surv_geneanet_2_depts_c <- surv_geneanet_2_depts_c %>% bind_rows()

head(surv_geneanet_2_depts_c)

# A tibble: 6 x 14
  annee exposition_h exposition_f   age deces_h deces_f  sx_h  sx_f
  <int>        <dbl>        <int> <int>   <dbl>   <int> <dbl> <dbl>
1  1800         2231         2268     0   216       161 0.903 0.929
2  1801         2015         2107     1    67.0      90 0.873 0.889
3  1802         1948         2017     2    31.0      34 0.859 0.874
4  1803         1917         1983     3    29.0      42 0.846 0.856
5  1804         1888         1941     4    22.0      29 0.836 0.843
6  1805         1866         1912     5    13.0      15 0.831 0.836
# ... with 6 more variables: qxt_h <dbl>, qxt_f <dbl>, fm_h <dbl>, fm_f
#   <dbl>, annee_naissance <int>, dept <chr>

Puis, on agrège par âge et département, en calculant les moyennes pour la force de mortalité et pour la fonction de survie.

surv_geneanet_2_depts <- 
  surv_geneanet_2_depts_c %>% 
  group_by(age, dept) %>% 
  summarise(fm_h = mean(fm_h, na.rm=T), fm_f = mean(fm_f, na.rm=T),
            sx_h = mean(sx_h, na.rm=T), sx_f = mean(sx_f, na.rm=T))

On prépare les données pour que le format corresponde à ce qui est attendu par les fonctions de graphiques du package ggplot2.

surv_geneanet_2_dept_sx <- 
  surv_geneanet_2_depts %>% 
  select(dept, age, sx_h, sx_f) %>% 
  gather(sexe, sx, sx_h, sx_f) %>% 
  mutate(sexe = factor(sexe, levels = c("sx_h", "sx_f"), labels = c(1,2))) %>% 
  mutate(sexe = as.character(sexe))

df_plot_survie_2 <- 
  surv_geneanet_2_dept_sx %>% 
  mutate(sexe = factor(sexe, levels = c("2","1"), labels = c("Femmes", "Hommes")))

p <- 
  ggplot(data = df_plot_survie_2,
         aes(x = age, y = sx, colour = sexe, group = interaction(dept, sexe))) +
  geom_line(alpha = .2) +
  xlab("âge") + ylab("Probabilité") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))
p

3.2 Estimation de la force de mortalité

À partir des estimations réalisées dans la section précédente, on peut se pencher sur la force de mortalité, c’est-à-dire la probabilité, pour un individu ayant atteint un âge donné, de décéder dans le courant de l’année.

3.2.1 Pour les données de Vallin & Meslé comparées à celles de Geneanet

df_plot_mx <- 
  surv_hist %>% 
  select(age, sexe, mx) %>% 
  mutate(mx = log(mx)) %>% 
  mutate(type = "Historique") %>% 
  bind_rows(
    surv_geneanet %>% 
      select(age, fm_h, fm_f) %>% 
      gather(sexe, mx, fm_h, fm_f) %>% 
      mutate(type = "Geneanet") %>% 
      mutate(sexe = ifelse(sexe == "fm_h", yes = "1", no = sexe),
             sexe = ifelse(sexe == "fm_f", yes = "2", no = sexe))
  ) %>% 
  mutate(sexe = factor(sexe, levels = c("2","1"), labels = c("Femmes", "Hommes")))

head(df_plot_mx)

# A tibble: 6 x 4
# Groups: age [3]
    age sexe      mx type      
  <int> <fct>  <dbl> <chr>     
1     6 Hommes -4.29 Historique
2     6 Femmes -4.31 Historique
3     7 Hommes -4.49 Historique
4     7 Femmes -4.49 Historique
5     8 Hommes -4.58 Historique
6     8 Femmes -4.62 Historique

p_mx <- 
  ggplot(data = df_plot_mx, aes(x = age, y = mx, colour = sexe, linetype = type)) +
  geom_line() +
  scale_colour_manual("", values = c("Femmes" = "red", "Hommes" = "blue")) +
  scale_linetype_manual("", values = c("Historique" = "dashed", "Geneanet" = "solid")) +
  xlab("âge") + ylab("log(mu_x)") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))
p_mx

Et enfin, on peut tracer ces deux graphiques sur une seule figure :

df_plot <- 
  df_plot_survie %>% select(-mx) %>% rename(value = sx) %>% mutate(variable = "sx") %>% 
  bind_rows(df_plot_mx %>% rename(value = mx) %>% mutate(variable = "mx"))

p_survie <- 
  ggplot(data = df_plot %>% mutate(variable = factor(variable, levels = c("sx", "mx"),
labels = c("Probabilité de survie", "Force de mortalité (en log)"))),
         aes(x = age, y = value, colour = sexe, linetype = type)) +
  geom_line() +
  scale_colour_manual("", values = c("Femmes" = "red", "Hommes" = "blue")) +
  scale_linetype_manual("", values = c("Historique" = "dashed", "Geneanet" = "solid")) +
  xlab("age") + ylab(NULL) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  facet_wrap(~variable, scales = "free_y", nrow = 1)

p_survie

3.2.2 Pour l’ensemble des individus de Geneanet, par département

df_plot_mx_2 <- 
  surv_geneanet_2_depts %>% 
  select(dept, age, fm_h, fm_f) %>% 
  gather(sexe, mx, fm_h, fm_f) %>% 
  mutate(sexe = ifelse(sexe == "fm_h", yes = "1", no = sexe),
         sexe = ifelse(sexe == "fm_f", yes = "2", no = sexe)) %>% 
  mutate(sexe = factor(sexe, levels = c("2","1"), labels = c("Femmes", "Hommes")))



p_mx_2 <- 
  ggplot(data = df_plot_mx_2, aes(x = age, y = mx, colour = sexe, group = interaction(dept, sexe))) +
  geom_line(alpha = .2) +
  scale_colour_manual("", values = c("Femmes" = "red", "Hommes" = "blue")) +
  xlab("\\^age") + ylab("$log(\\mu\\_x$)") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))

p_mx_2+theme_paper()


# Les deux graphiques combinés

df_plot_2 <- 
  df_plot_survie_2 %>% rename(value = sx) %>% mutate(variable = "sx") %>% 
  bind_rows(df_plot_mx_2 %>% rename(value = mx) %>% mutate(variable = "mx"))

p_survie_2 <- 
  ggplot(data = df_plot_2 %>% 
           mutate(variable = factor(variable, levels = c("sx", "mx"),
                                    labels = c("Probabilit\\'e de survie", "Force de mortalit\\'e (en log)"))),
         aes(x = age, y = value, colour = sexe, group = interaction(dept, sexe))) +
  geom_line(alpha=.1) +
  scale_colour_manual("", values = c("Femmes" = "red", "Hommes" = "blue")) +
  xlab("\\^age") + ylab(NULL) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  facet_wrap(~variable, scales = "free_y", nrow = 1)


p_survie_2+theme_paper()

3.3 Mortalité régionale par cohorte

df_plot_mx_2 <- 
  surv_geneanet_2_depts %>% 
  select(dept, age, fm_h, fm_f) %>% 
  gather(sexe, mx, fm_h, fm_f) %>% 
  mutate(sexe = ifelse(sexe == "fm_h", yes = "1", no = sexe),
         sexe = ifelse(sexe == "fm_f", yes = "2", no = sexe)) %>% 
  mutate(sexe = factor(sexe, levels = c("2","1"), labels = c("Femmes", "Hommes")))

p_mx_2 <- 
  ggplot(data = df_plot_mx_2, aes(x = age, y = mx, colour = sexe, group = interaction(dept, sexe))) +
  geom_line(alpha = .2) +
  scale_colour_manual("", values = c("Femmes" = "red", "Hommes" = "blue")) +
  xlab("âge") + ylab("log(mu_x)") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))

p_mx_2

En combinant les deux graphiques :

df_plot_2 <- 
  df_plot_survie_2 %>% rename(value = sx) %>% mutate(variable = "sx") %>% 
  bind_rows(df_plot_mx_2 %>% rename(value = mx) %>% mutate(variable = "mx"))

p_survie_2 <- 
  ggplot(data = df_plot_2 %>% 
           mutate(variable = factor(variable, levels = c("sx", "mx"),
                                    labels = c("Probabilité de survie", "Force de mortalité (en log)"))),
         aes(x = age, y = value, colour = sexe, group = interaction(dept, sexe))) +
  geom_line(alpha=.1) +
  scale_colour_manual("", values = c("Femmes" = "red", "Hommes" = "blue")) +
  xlab("âge") + ylab(NULL) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  facet_wrap(~variable, scales = "free_y", nrow = 1)


p_survie_2

3.4 Espérance de vie à la naissance

À l’aide des estimations des fonctions de survie pour les femmes et les hommes, il est aisé de calculer les espérances de vie à la naissance des individus.

3.4.1 Pour Vallin et Meslé (2010)

Calculons l’espérance de vie à la naissance avec les données de Vallin et Meslé.

# Espérance de vie à la naissance pour la cohorte de 1806, données de Vallin et Meslé.
vallin_1806 <- 
  survie_diag(annee_naiss_depart = 1806, sexe = "1", age_min = 0) %>% 
  bind_rows(survie_diag(annee_naiss_depart = 1806, sexe = "2", age_min = 0))
vallin_1806 %>% 
  group_by(sexe) %>% 
  summarise(e = sum(cumprod(exp(-mx)), na.rm=T))

3.4.2 Pour l’ensemble des aïeux

Nous construisons une fonction effectuant le calcul.

calcul_esperance_vie_naissance_aieux <- function(df, summ = FALSE, age_min = 0){
  res <- 
    df %>% 
    filter(age_min >= age_min) %>% 
    left_join(
      df %>% 
        # group_by(annee_naissance, dept) %>% 
        filter(age == 0) %>% 
        select(annee_naissance, exposition_h, exposition_f) %>% 
        rename(exposition_f1 = exposition_f, exposition_h1 = exposition_h),
      by = c("annee_naissance")
    ) %>% 
    mutate(e_h = exposition_h / exposition_h1,
           e_f = exposition_f / exposition_f1) %>% 
    group_by(annee_naissance) %>% 
    summarise(e_h = sum(e_h)-0.5,
              e_f = sum(e_f)-0.5) %>% 
    ungroup()
  if(summ){
    res <- res %>% summarise(e_h = mean(e_h), e_f = mean(e_f))
  }
  res
}

Nous appliquons cette fonction à notre table de mortalité concernant la totalité des aïeux.

esperance_vie_aieux <- 
  calcul_esperance_vie_naissance_aieux(surv_geneanet_2_tot_c)

esperance_vie_aieux %>% 
  knitr::kable(digit=2)

annee_naissance	e_h	e_f
1800	45.27	43.31
1801	43.66	41.88
1802	43.06	41.32
1803	42.60	40.93
1804	43.12	41.62

esperance_vie_aieux_mean <- 
  esperance_vie_aieux %>% 
  summarise(e_h = mean(e_h), e_f = mean(e_f))
esperance_vie_aieux_mean

# A tibble: 1 x 2
    e_h   e_f
  <dbl> <dbl>
1  43.5  41.8

3.4.3 Par département

calcul_esperance_vie_naissance <- function(df, summ = FALSE, age_min = 0){
  res <- 
    df %>% 
    filter(age_min >= age_min) %>% 
    left_join(
      df %>% 
        # group_by(annee_naissance, dept) %>% 
        filter(age == 0) %>% 
        select(annee_naissance, dept, exposition_h, exposition_f) %>% 
        rename(exposition_f1 = exposition_f, exposition_h1 = exposition_h),
      by = c("annee_naissance", "dept")
    ) %>% 
    mutate(e_h = exposition_h / exposition_h1,
           e_f = exposition_f / exposition_f1) %>% 
    group_by(annee_naissance, dept) %>% 
    summarise(e_h = sum(e_h)-0.5,
              e_f = sum(e_f)-0.5) %>% 
    ungroup()
  if(summ){
    res <- res %>% summarise(e_h = mean(e_h), e_f = mean(e_f))
  }
  res
}

Les espérances de vie à la naissance pour chaque cohorte s’obtient de la manière suivante :

esperance_vie <- 
  calcul_esperance_vie_naissance(surv_geneanet_2_depts_c)
head(esperance_vie)

# A tibble: 6 x 4
  annee_naissance dept    e_h   e_f
            <int> <chr> <dbl> <dbl>
1            1800 F01    42.9  40.3
2            1800 F02    43.8  42.5
3            1800 F03    44.7  38.2
4            1800 F04    44.6  38.3
5            1800 F05    39.2  32.3
6            1800 F06    44.0  48.0

Il suffit alors d’effectuer la moyenne par département pour obtenir une estimation pour l’ensemble des cohortes de 1800 à 1804 :

esperance_vie_mean <- 
  esperance_vie %>% 
  group_by(dept) %>% 
  summarise(e_h = mean(e_h), e_f = mean(e_f))

head(esperance_vie_mean)

# A tibble: 6 x 3
  dept    e_h   e_f
  <chr> <dbl> <dbl>
1 F01    41.9  40.5
2 F02    41.1  39.9
3 F03    41.3  36.7
4 F04    42.7  37.3
5 F05    36.5  31.6
6 F06    42.7  42.5

Une représentation par le biais d’une carte choroplèthe est réalisable, en utilisant à nouveau la carte de france france_1800 dont la création est détaillée en annexe.

esperance_vie_map <- 
  france_1800 %>% 
  left_join(esperance_vie_mean, by = c("departement" = "dept"))

p_map <- 
  ggplot(esperance_vie_map %>% 
           gather(sexe, esperance, e_h, e_f) %>% 
           mutate(sexe = factor(sexe, levels = c("e_f", "e_h"), labels = c("Femmes", "Hommes"))),
         aes(x = long, y = lat, group = group, fill = esperance)) +
  geom_polygon() +
  facet_grid(~sexe) +
  scale_x_continuous(limits = c(-6, 10)) +
  scale_y_continuous(limits = c(41, 51.5)) +
  scale_colour_manual("Sexe", values = c("Femmes" = "blue", "Hommes" = "red"), na.value = "#377eb8") +
  scale_fill_gradient(low = "yellow", high = "red") +
  coord_quickmap() +
  theme(legend.position = "bottom", text = element_text(size = 6),
        legend.key.height   = unit(1, "line"),
        legend.key.width    = unit(1.5, "line"))


p_map

Il est intéressant de comparer les valeurs avec celles qui existent dans la littérature. Nous reportons les valeurs trouvées dans Walle (1973).

esperances_vandewalle <- 
  c('F01', 'Ain', '32.6',
    'F02', 'Aisne', '35.9',
    'F03', 'Allier', '31.7',
    'F04', 'Alpes-de-Haute-Provence', '34.2',
    'F05', 'Hautes-Alpes', '38.3',
    'F06', 'Alpes-Maritimes', NA,
    'F07', 'Ardèche', '41.9',
    'F08', 'Ardennes', '40.7',
    'F09', 'Ariège', '45.2',
    'F10', 'Aube', '34.2',
    'F11', 'Aude', '39.4',
    'F12', 'Aveyron', '44.6',
    'F13', 'Bouches-du-Rhône', NA,
    'F14', 'Calvados', '47.7',
    'F15', 'Cantal', '43.0',
    'F16', 'Charente', '35.8',
    'F17', 'Charente-Maritime', '32.8',
    'F18', 'Cher', '32.3',
    'F19', 'Corrèze', '34.8',
    'F20', 'Corse-du-Sud', '39',
    'F21', "Côte-d'Or", '35.4',
    'F22', "Côtes-d'Armor", '37.9',
    'F23', 'Creuse', '35.2',
    'F24', 'Dordogne', '34.3',
    'F25', 'Doubs', '39.6',
    'F26', 'Drôme', '37.3',
    'F27', 'Eure', '40.6',
    'F28', 'Eure-et-Loir', '30.1',
    'F29', 'Finistère', '29.2',
    'F30', 'Gard', '36.9',
    'F31', 'Haute-Garonne', '38.0',
    'F32', 'Gers', '40.4',
    'F33', 'Gironde', '37.7',
    'F34', 'Hérault', '37.8',
    'F35', 'Ille-et-Vilaine', '32.4',
    'F36', 'Indre', '29.1',
    'F37', 'Indre-et-Loire', '35.2',
    'F38', 'Isère', '37.9',
    'F39', 'Jura', '37.5',
    'F40', 'Landes', '30.7',
    'F41', 'Loir-et-Cher', '24.9',
    'F42', 'Loire', '37.0',
    'F43', 'Haute-Loire', '38.8',
    'F44', 'Loire-Atlantique', '39.4',
    'F45', 'Loiret', '24.7',
    'F46', 'Lot', '35.8',
    'F47', 'Lot-et-Garonne', '41.6',
    'F48', 'Lozère', '45.3',
    'F49', 'Maine-et-Loire', '38.4',
    'F50', 'Manche', '45.9',
    'F51', 'Marne', '38.1',
    'F52', 'Haute-Marne', '33.2',
    'F53', 'Mayenne', '34.7',
    'F54', 'Meurthe-et-Moselle', '36.0',
    'F55', 'Meuse', '35.0',
    'F56', 'Morbihan', '32.1',
    'F57', 'Moselle', '39.2',
    'F58', 'Nièvre', '26.9',
    'F59', 'Nord', '36.3',
    'F60', 'Oise', '35.8',
    'F61', 'Orne', '39.1',
    'F62', 'Pas-de-Calais', '38.0',
    'F63', 'Puy-de-Dôme', '35.9',
    'F64', 'Pyrénées-Atlantiques', '44.2',
    'F65', 'Hautes-Pyrénées', '46.5',
    'F66', 'Pyrénées-Orientales', '32.1',
    'F67', 'Bas-Rhin', '37.1',
    'F68', 'Haut-Rhin', '39.9',
    'F69', 'Rhône', NA,
    'F70', 'Haute-Saône', '34.8',
    'F71', 'Saône-et-Loire', '32.8',
    'F72', 'Sarthe', '32.4',
    'F73', 'Savoie', NA,
    'F74', 'Haute-Savoie', NA,
    'F75', 'Paris', NA,
    'F76', 'Seine-Maritime', '44.0',
    'F77', 'Seine-et-Marne', '28.6',
    'F78', 'Essonne', NA,
    'F79', 'Deux-Sèvres', '41.0',
    'F80', 'Somme', '41.0',
    'F81', 'Tarn', '39.6',
    'F82', 'Tarn-et-Garonne', '37.3',
    'F83', 'Var', '33.9',
    'F84', 'Vaucluse', '33.4',
    'F85', 'Vendée', '36.6',
    'F86', 'Vienne', '46.0',
    'F87', 'Haute-Vienne', '29.5',
    'F88', 'Vosges', '38.3',
    'F89', 'Yonne', '29.2') %>% 
  matrix(ncol = 3, byrow = TRUE) %>% 
  data.frame(stringsAsFactors = FALSE) %>% 
  tbl_df()

colnames(esperances_vandewalle) <- c("dept", "nom_sp", "e_f_walle")

esperances_vandewalle <- 
  esperances_vandewalle %>% 
  mutate(e_f_walle = as.numeric(e_f_walle))

Nous joignons nos estimations régionalles avec celles d’Etienne van de Walle.

esperance_vie_mean_2 <- 
  esperance_vie_mean %>% 
  left_join(esperances_vandewalle)

head(esperance_vie_mean_2)

# A tibble: 6 x 5
  dept    e_h   e_f nom_sp                  e_f_walle
  <chr> <dbl> <dbl> <chr>                       <dbl>
1 F01    41.9  40.5 Ain                          32.6
2 F02    41.1  39.9 Aisne                        35.9
3 F03    41.3  36.7 Allier                       31.7
4 F04    42.7  37.3 Alpes-de-Haute-Provence      34.2
5 F05    36.5  31.6 Hautes-Alpes                 38.3
6 F06    42.7  42.5 Alpes-Maritimes              NA  

L’espérance de vie moyenne à la naissance des femmes vaut, pour notre échantillon et pour celui de van de Walle :

mean(esperance_vie_mean_2$e_f, na.rm=T)

[1] 42.20441

mean(esperance_vie_mean_2$e_f_walle, na.rm=T)

[1] 36.67805

Notre valeur moyenne estimée est supérieure. Toutefois, nous souhaitons voir si les mêmes effets géographiques apparaissent avec nos données. Nous calculons alors les écarts par départements par rapport à la valeur moyenne des départements, pour nos données, ainsi que pour celles de la littérature.

esperance_vie_mean_2 <- 
  esperance_vie_mean_2 %>% 
  mutate(e_f_3 = e_f - mean(e_f, na.rm=T),
         e_f_walle_3 = e_f_walle - mean(e_f_walle, na.rm=T))

head(esperance_vie_mean_2)

# A tibble: 6 x 7
  dept    e_h   e_f nom_sp                  e_f_walle   e_f_3 e_f_walle_3
  <chr> <dbl> <dbl> <chr>                       <dbl>   <dbl>       <dbl>
1 F01    41.9  40.5 Ain                          32.6 - 1.70      - 4.08 
2 F02    41.1  39.9 Aisne                        35.9 - 2.30      - 0.778
3 F03    41.3  36.7 Allier                       31.7 - 5.49      - 4.98 
4 F04    42.7  37.3 Alpes-de-Haute-Provence      34.2 - 4.92      - 2.48 
5 F05    36.5  31.6 Hautes-Alpes                 38.3 -10.6         1.62 
6 F06    42.7  42.5 Alpes-Maritimes              NA     0.325      NA    

Puis, nous représentons le résultats sur des cartes choroplèthes.

esperance_vie_map <- 
  france_1800 %>% 
  left_join(esperance_vie_mean_2 %>% 
              select(dept, e_f_3, e_f_walle_3), by = c("departement" = "dept"))

p_map <- 
  ggplot(esperance_vie_map %>% 
           gather(type, esperance, e_f_3, e_f_walle_3) %>% 
           mutate(type = factor(type, levels = c("e_f_3", "e_f_walle_3"), labels = c("Geneanet", "van de Walle"))),
         aes(x = long, y = lat, group = group, fill = esperance)) +
  geom_polygon(col = "grey") +
  facet_wrap(~type) +
  scale_x_continuous(limits = c(-6, 10)) +
  scale_y_continuous(limits = c(41, 51.5)) +
  scale_fill_gradient2("Écart à la moyenne départementale (années)", low = "red", mid = "white", high = "#71BC78", midpoint = 0, na.value = "grey",
                       guide = guide_colourbar(direction = "horizontal",
                                               barwidth = 30,
                                               title.position = "bottom",
                                               title.hjust = .5)) +
  coord_quickmap() +
  theme(legend.position = "bottom", text = element_text(size = 6),
        legend.key.height   = unit(1, "line"),
        legend.key.width    = unit(1.5, "line")) +
  theme(text = element_text(size = 16), legend.position = "bottom")


p_map

3.5 Intervalles de confiance pour les espérances de vie à la naissance

Nous allons estimer des intervalles de confiance pour les espérances de vie à la naissance. Deux méthodes sont envisagées. Une première est d’utiliser ce que propose Carlo Giovanni Camarda sur son site (https://sites.google.com/site/carlogiovannicamarda/r-stuff/life-expectancy-confidence-interval). Une seconde est fortement inspirée de cette méthode, mais s’avère plus simple à mettre en oeuvre. Les deux méthodes sont fondées sur du bootstrap.

Nous récupérons la fonction lifetable.qx() servant à constuire une table de mortalité à partir des probabilités de décès.

## function for constructing a lifetable starting from probabilities
lifetable.qx <- function(x, qx, sex="M", ax=NULL, last.ax=5.5){
  m <- length(x)
  n <- c(diff(x), NA)
  qx[is.na(qx)] <- 0
  if(is.null(ax)){
    ax <- rep(0,m)
    if(x[1]!=0 | x[2]!=1){
      ax <- n/2
      ax[m] <- last.ax
    }else{    
      if(sex=="F"){
        if(qx[1]>=0.1){
          ax[1] <- 0.350
        }else{
          ax[1] <- 0.05 + 3*qx[1]
        }
      }
      if(sex=="M"){
        if(qx[1]>=0.1){
          ax[1] <- 0.33
        }else{
          ax[1] <- 0.0425 + 2.875*qx[1]
        }
      }
      ax[-1] <- n[-1]/2
      ax[m] <- last.ax
    }
  }
  px  <- 1-qx
  lx  <- cumprod(c(1,px))*100000
  dx  <- -diff(lx)
  Lx  <- n*lx[-1] + ax*dx
  lx <- lx[-(m+1)]
  Lx[m] <- lx[m]*last.ax
  Lx[is.na(Lx)] <- 0 ## in case of NA values
  Lx[is.infinite(Lx)] <- 0 ## in case of Inf values
  Tx  <- rev(cumsum(rev(Lx)))
  ex  <- Tx/lx
  data.frame(x, n, ax, qx, px, lx, dx, Lx, Tx, ex)
}

Cette fonction est appelée par celle qui suit, cette dernière permettant quant à elle de calculer les intervalles de confiance pour l’espérance de vie à un âge donné. Elle est légèrement adaptée des codes proposés par Camarda.

#' ex_ci
#' @description Compute a confidence interval for life expectancy
#' https://sites.google.com/site/carlogiovannicamarda/r-stuff/life-expectancy-confidence-interval
#' @param LT Life table
#' @param which.x Which age for the life expectancy
#' @param level level of the confidence interval
#' @param ns number of samples
#' @param sex 'F' (female) of 'M' (male)
#' @param year_birth year of birth of the cohort
#' @param dept_birth departement of birth
#' which.x <- 0; level <- 0.95; ns <- 1000; sex="M"; year_birth <- 1800 ; dept_birth <- "F29"
ex_ci <- function(LT, which.x=0, level=0.95, ns=1000, sex="M", year_birth = 1800, dept_birth = "F29"){
  ## ages
  x <- LT$age
  if(sex == "M"){
    ## number of deaths
    Dx <- LT$deces_h
    ## estimated probs
    qx <- LT$qxt_h
  }else{
    ## number of deaths
    Dx <- LT$deces_f
    ## estimated probs
    qx <- LT$qxt_f
  }
  
  ## number of ages
  m <- nrow(LT)
  
  ## trials for binomial, rounded
  Ntil <- round(Dx/qx)
  ## ax for last age
  # last.ax <- LT$ax[m]
  ## simulated death counts
  ## from Binomial distribution
  Y <- suppressWarnings(matrix(rbinom(m*ns,
                                      Ntil,
                                      qx),
                               m,ns))
  ## simulated probabilities
  QX <- Y/Ntil
  ## which age?
  wh <- which(x==which.x)
  
  ## for all replicates, compute life expectancy
  fun.ex <- function(qx){
    return(lifetable.qx(x=x, qx, sex=sex)$ex[wh])
  }
  exsim <- apply(QX, 2, fun.ex)
  
  ## confidence interval
  CI <- quantile(exsim,
                 probs = c((1-level)/2,
                           1 - (1-level)/2))
  
  data.frame(ex_mean = mean(exsim),
             ex_inf = CI[[1]],
             ex_sup = CI[[2]],
             age = which.x,
             sexe = sex,
             annee_naissance = year_birth,
             dept = dept_birth,
             stringsAsFactors = FALSE)
  
} # Fin de ex_ci()

Nous introduisons également une fonction permettant d’utiliser la seconde méthode. Cette méthode consiste à considérer une population \(E_0\) née une année donnée (e.g., 1800) ayant une probabilité de décéder dans l’année (\(\widehat{q}_x\), estimée sur nos données pour cette cohorte). Il s’agit de procéder de manière itérative. À l’âge \(x\), s’il reste \(E_x\) personnes en vie, le nombre de décès \(D_x\) est tiré suivant une loi binomiale \(\mathcal{B}(E_x,\widehat{q}_x)\). Ensuite, nous posons \(E_{x+1}=E_x - D_x\), c’est-à-dire que nous considérons le nombre de personnes en vie à l’âge suivant égal au nombre de personnes en vie à l’âge \(x\) auquel nous retranchons le nombre de décès qui vient d’être tiré. Nous tirons alors un nouveau nombre de décès pour l’âge suivant, et ainsi de suite. Ces étapes itératives sont reproduites de manière indépendante dans \(1~000\) échantillons. Dans chacun de ces échantillons, nous calculons l’espérance de vie à la naissance selon les tirages effectués. Les quantiles empiriques d’ordres \(0,025\) et \(0,975\) deviennent alors les bornes de l’intervalle de confiance à \(95\%\) pour l’espérance de vie à la naissance des individus de la cohorte étudiée. La fonction simu() ci-après permet d’effectuer un tirage. Il s’agira de faire appel à simu() de manière itérative pour obtenir différents échantillons.

#' simu
#' @description Estime des IC par bootstrap pour les esperances de vie a la naissance
#' a partir d'une table de vie
#' @param i valeur seed
#' @param LT Table de mortalite
simu <- function(i, LT){
  set.seed(i)
  d_h=d_f=rep(0,111)
  e_h=LT$exposition_h
  e_f=LT$exposition_f
  q_h=LT$qxt_h
  q_f=LT$qxt_f
  
  # Last not NA value :
  ind_max_not_na_h <- max(which(!is.na(q_h)))
  ind_na_h <- which(is.na(q_h))
  q_h[ind_na_h[ind_na_h> ind_max_not_na_h]]=1
  
  ind_max_not_na_f <- max(which(!is.na(q_f)))
  ind_na_f <- which(is.na(q_f))
  q_f[ind_na_f[ind_na_f> ind_max_not_na_f]]=1
  
  if(any(is.na(q_h))){
    loess_h <- loess(q_h~age, data = data.frame(q_h = q_h, age = 1:length(q_h)))
    ind_missing_h <- which(is.na(q_h))
    predicted_h <- predict(loess_h, newdata = data.frame(age = ind_missing_h))
    predicted_h[predicted_h>1] <- 1
    predicted_h[predicted_h<0] <- 0
    q_h[ind_missing_h] <- predicted_h
  }
  
  if(any(is.na(q_f))){
    loess_f <- loess(q_f~age, data = data.frame(q_f = q_f, age = 1:length(q_f)))
    ind_missing_f <- which(is.na(q_f))
    predicted_f <- predict(loess_f, newdata = data.frame(age = ind_missing_f))
    predicted_f[predicted_f>1] <- 1
    predicted_f[predicted_f<0] <- 0
    q_f[ind_missing_f] <- predicted_f
  }
  
  for(t in 1:110){
    D_h = rbinom(1,size=e_h[t],prob=q_h[t])
    d_h[t] = D_h
    e_h[t+1] = e_h[t]-d_h[t]
    D_f = rbinom(1,size=e_f[t],prob=q_f[t])
    d_f[t] = D_f
    e_f[t+1] = e_f[t]-d_f[t]  
  }
  s_h=e_h/e_h[1]
  s_f=e_f/e_f[1]
  
  list(E_h=sum(s_h)-.5,E_f=sum(s_f)-.5)
}# Fin de simu()

3.5.1 Pour Vallin et Meslé

LT_h <- lifetable_fr %>% filter(Year1 == 1806, Sex == 1, TypeLT == 1)
LT_f <- lifetable_fr %>% filter(Year1 == 1806, Sex == 2, TypeLT == 1)

#'
#'@param LT Life table
#'@param ns number of samples to draw
#'@param which.x age for which to compute life expectancy
#'@param level level for the confidence interval
ex_ci_vallin <- function(LT, ns=1000, which.x=0, level=0.95){
  ## sex
  sex <- LT$Sex %>% unique()
  ## ages
  x <- LT$Age
  ## number of deaths
  Dx <- LT$`d(x)`
  ## estimated probs
  qx <- LT$`q(x)`
  
  ## number of ages
  m <- nrow(LT)
  
  ## trials for binomial, rounded
  Ntil <- round(Dx/qx)
  ## ax for last age
  # last.ax <- LT$ax[m]
  ## simulated death counts
  ## from Binomial distribution
  Y <- suppressWarnings(matrix(rbinom(m*ns,
                                      Ntil,
                                      qx),
                               m,ns))
  ## simulated probabilities
  QX <- Y/Ntil
  ## which age?
  wh <- which(x==which.x)
  
  ## for all replicates, compute life expectancy
  fun.ex <- function(qx){
    return(lifetable.qx(x=x, qx, sex=sex)$ex[wh])
  }
  exsim <- apply(QX, 2, fun.ex)
  
  ## confidence interval
  CI <- quantile(exsim,
                 probs = c((1-level)/2,
                           1 - (1-level)/2))
  
  data.frame(ex_mean = mean(exsim),
             ex_inf = CI[[1]],
             ex_sup = CI[[2]],
             age = which.x,
             sexe = sex,
             annee_naissance = which.x,
             stringsAsFactors = FALSE)
  
}# End of ex_ci_vallin()

ex_vallin <- ex_ci_vallin(LT = LT_h) %>% 
  bind_rows(
    ex_ci_vallin(LT = LT_f)
  )

ex_vallin

3.5.2 Pour les aïeux

Nous proposons une fonction, e_ic_cohorte(), qui à partir d’une année de naissance, estime par bootstrap les intervalles de confiance pour l’espérance de vie à la naissance de la cohorte d’individus nés cette année.

#' e_ic_cohorte
#' @param year_birth annee de naissance de la cohorte
#' @param ns nombre d'échantillons
#' #' @param level niveau de l'intervalle de confiance
e_ic_cohorte <- function(year_birth, ns=1000, level=0.95){
  LT <- surv_geneanet_2_tot_c %>% filter(annee_naissance == year_birth)
  E=lapply(1:ns,simu, LT=LT)
  BE=matrix(unlist(E),nrow=2)
  
  ex_mean <- apply(BE,1,mean)
  ex_inf <- apply(BE,1,function(x) quantile(x,(1-level)/2))
  ex_sup <- apply(BE,1,function(x) quantile(x,1 - (1-level)/2))
  
  data.frame(ex_mean = ex_mean,
             ex_inf = ex_inf,
             ex_sup = ex_sup,
             sexe = c("M", "F"),
             annee_naissance = year_birth,
             stringsAsFactors = FALSE)
}# Fin de e_ic_cohorte()

ex_aieux <- pblapply(1800:1804, e_ic_cohorte)
ex_aieux <- bind_rows(ex_aieux)

ex_aieux_avg <- 
  ex_aieux %>% 
  group_by(sexe) %>% 
  summarise(ex_mean = mean(ex_mean),
            ex_inf = mean(ex_inf),
            ex_sup = mean(ex_sup))

ex_aieux_avg

# A tibble: 2 x 4
  sexe  ex_mean ex_inf ex_sup
  <chr>   <dbl>  <dbl>  <dbl>
1 F        41.8   41.6   42.0
2 M        43.5   43.3   43.7

3.5.3 Pour les aïeux, par département

Ici, nous proposons d’utiliser les deux méthodes d’estimation des intervalles de confiance.

3.5.3.1 Première méthode

Nous appliquons la méthode de Carlo Giovanni Camarda (https://sites.google.com/site/carlogiovannicamarda/r-stuff/life-expectancy-confidence-interval).

Nous bouclons sur les années de naissance des aïeux, pour chaque sexe, et pour chaque département.

annees <- 1800:1804
sexes <- c("F", "M")
depts <- sort(unique(surv_geneanet_2_depts_c$dept))

n_tot <- length(annees)*length(sexes)*length(depts)
pb <- txtProgressBar(min = 1, max = n_tot, style = 3)
i <- 1
ex_depts <- NULL
for(annee in annees){
 for(sexe in sexes) {
   for(dept in depts){
     LT <- surv_geneanet_2_depts_c %>% filter(annee_naissance == annee, dept == dept)
     tmp <- ex_ci(LT=LT, year_birth = annee, sex = sexe, dept_birth = dept)
     ex_depts <- bind_rows(ex_depts, tmp)
     rm(tmp, LT)
     setTxtProgressBar(pb, i)
     i <- i+1
   }
 }
}

head(ex_depts)

   ex_mean   ex_inf   ex_sup age sexe annee_naissance dept
1 40.27501 38.32497 42.34023   0    F            1800  F01
2 42.42850 40.84798 43.98664   0    F            1800  F02
3 38.24195 36.38821 40.12682   0    F            1800  F03
4 38.34529 34.82395 41.78034   0    F            1800  F04
5 32.26114 29.72819 35.05179   0    F            1800  F05
6 47.84448 44.62942 51.02613   0    F            1800  F06

3.5.3.2 Seconde méthode

En utilisant notre méthode à présent.

#' ex_ci
#' @description Compute a confidence interval for life expectancy
#' From an initial population $E_0$ born in `year_birth`, et a probability of dying within the current year
#' $\widehat{q}_x$ previously estimated, proceeds as follows:
#' For age $x$, if $E_x$ people are still alive, the number of deaths $D_x$ is drawn from a Binomial $B(E_x,q_x)$ distribution
#' and we then posit $E_{x+1}=E_x-D_x$.
#' @param level level of the confidence interval
#' @param ns number of samples
#' @param year_birth year of birth of the cohort
#' @param dept_birth departement of birth
#' which.x <- 0; level <- 0.95; ns <- 1000; sex="M"; year_birth <- 1800 ; dept_birth <- "F29"
ex_ci_2 <- function(level=0.95, ns=1000, year_birth = 1800, dept_birth = "F29"){
  
  LT <- surv_geneanet_2_depts_c %>% filter(annee_naissance == year_birth, dept == dept_birth)
  
  E=lapply(1:ns,simu, LT=LT)
  BE=matrix(unlist(E),nrow=2)
  
  ex_mean <- apply(BE,1,mean)
  ex_inf <- apply(BE,1,function(x) quantile(x,(1-level)/2))
  ex_sup <- apply(BE,1,function(x) quantile(x,1 - (1-level)/2))
  
  data.frame(ex_mean = ex_mean,
             ex_inf = ex_inf,
             ex_sup = ex_sup,
             sexe = c("M", "F"),
             annee_naissance = year_birth,
             dept = dept_birth,
             stringsAsFactors = FALSE)
  
} # Fin de ex_ci_2()

Comme cette méthode effectue les estimations pour chaque sexe dans la fonction simu(), nous itérons uniquement sur les années de naissance et sur les départements.

n_tot <- length(annees)*length(depts)
pb <- txtProgressBar(min = 1, max = n_tot, style = 3)
i <- 1
ex_depts_2 <- NULL
for(annee in annees){
  for(dept in depts){
    tmp <- ex_ci_2(year_birth = annee, dept_birth = dept)
    ex_depts_2 <- bind_rows(ex_depts_2, tmp)
    rm(tmp)
    setTxtProgressBar(pb, i)
    i <- i+1
  }
}

head(ex_depts_2)

   ex_mean   ex_inf   ex_sup sexe annee_naissance dept
1 41.02150 39.07937 42.96377    M            1800  F01
2 39.66530 37.69183 41.63957    F            1800  F01
3 43.74457 42.35724 45.04229    M            1800  F02
4 42.40019 40.97947 43.99872    F            1800  F02
5 44.09319 42.46796 45.83167    M            1800  F03
6 37.91267 36.07012 39.75813    F            1800  F03

ex_depts_avg <- 
  ex_depts_2 %>% 
  group_by(dept, sexe) %>% 
  summarise(ex_mean = mean(ex_mean),
            ex_inf = mean(ex_inf),
            ex_sup = mean(ex_sup))

head(ex_depts_avg)

# A tibble: 6 x 5
# Groups: dept [3]
  dept  sexe  ex_mean ex_inf ex_sup
  <chr> <chr>   <dbl>  <dbl>  <dbl>
1 F01   F        39.7   37.6   41.8
2 F01   M        40.9   38.9   42.9
3 F02   F        39.9   38.3   41.5
4 F02   M        41.0   39.5   42.4
5 F03   F        35.0   33.0   37.1
6 F03   M        39.4   37.4   41.4

Nous pouvons alors grapher les estimations des espérances de vie à la naissance accompagnées des intervalles de confiance.

esperance_vie_mean_ic <- 
esperance_vie_mean_2 %>% 
  left_join(
    ex_depts_avg %>% 
      filter(sexe == "F") %>% 
      rename(ex_mean_F = ex_mean, ex_inf_F = ex_inf, ex_sup_F = ex_sup) %>% 
      select(-sexe) %>% 
      left_join(
        ex_depts_avg %>% 
          filter(sexe == "M") %>% 
          rename(ex_mean_M = ex_mean, ex_inf_M = ex_inf, ex_sup_M = ex_sup) %>% 
          select(-sexe)
      )
  )

esperance_vie_mean_ic %>% 
  mutate(outside = e_f_walle <= ex_inf_F,
         inside = e_f_walle > ex_inf_F & e_f_walle < ex_sup_F) %>% 
  summarise(outside = sum(outside, na.rm=T),
            inside = sum(inside, na.rm=T))

p_error_bar_e <- 
  esperance_vie_mean_ic %>% 
  # mutate(departement = str_c(nom_sp, " (", str_sub(dept, 2), ")")) %>% 
  mutate(departement = str_sub(dept, 2)) %>% 
  select(departement, e_h, e_f, ex_inf_M, ex_sup_M, ex_inf_F, ex_sup_F) %>% 
  unite(values_f, e_f, ex_inf_F, ex_sup_F, sep = "@") %>% 
  unite(values_h, e_h, ex_inf_M, ex_sup_M, sep = "@") %>% 
  gather(sexe, value, values_h, values_f) %>% 
  separate(value, c("e", "e_inf", "e_sup"), sep = "@", convert = TRUE) %>% 
  mutate(echantillon = "Geneanet") %>% 
  filter(!is.na(departement)) %>% 
  ggplot(data = ., aes(x = departement, y = e, colour = sexe)) +
  geom_point() +
  geom_errorbar(aes(ymin = e_inf, ymax = e_sup), width=1, size = .8) +
  scale_colour_manual("Sexe", values = c("values_f" = "red", "values_h" = "blue"),
                      labels = c("values_f" = "Femmes", "values_h" = "Hommes")) +
  geom_point(data = 
               esperance_vie_mean_ic %>% 
                   # mutate(departement = str_c(nom_sp, " (", str_sub(dept, 2), ")")) %>% 
                   mutate(departement = str_sub(dept, 2)) %>% 
                   select(departement, e_f_walle) %>% 
                   mutate(sexe = "values_f") %>% 
                   rename(e = e_f_walle) %>% 
               mutate(echantillon = "van de Walle"),
             aes(x = departement, y = e, shape = echantillon), colour = "black") +
  scale_shape_manual("", values = c("van de Walle" = 17)) +
  theme(axis.text.x = element_text(angle = 90, size = .6*20)) +
  xlab("Département") + ylab("Espérance de vie à la naissance")

p_error_bar_e

Nous produisons également une représentation sous forme de carte.

esperance_vie_map_ic <- 
  france_1800 %>% 
  left_join(esperance_vie_mean_2 %>% 
              select(dept, e_f, e_h), by = c("departement" = "dept"))

Pour pouvoir ajouter les intervalles de confiance relatifs aux estimations de chaque département, nous récupérons les coordonnées de barycentres de chaque département.

# Barycentres
library(sp)
obtenir_spa_poly_sp <- function(un_dep){
  dep_tmp <- france_1800 %>%
    filter(departement == un_dep)
  
  obtenir_poly_groupe <- function(un_groupe){
    coords_group <- dep_tmp[which(dep_tmp$group %in% un_groupe), c("long", "lat")]
    return(Polygons(list(Polygon(coords_group)), un_groupe))
  }
  dep_polys <- lapply(unique(dep_tmp$group), obtenir_poly_groupe)
  dep_polys_sp <- SpatialPolygons(dep_polys, seq(1, length(dep_polys)))
  dep_polys_sp
}

france_1800_sp <- lapply(unique(france_1800$departement), obtenir_spa_poly_sp)
names(france_1800_sp) <- unique(france_1800$departement)
france_1800_barycentres <- lapply(france_1800_sp, function(x) rgeos::gCentroid(x))
france_1800_barycentres <- do.call("rbind", lapply(france_1800_barycentres, function(x) x@coords))
france_1800_barycentres <- cbind(france_1800_barycentres, departement = unique(france_1800$departement))
france_1800_barycentres <- data.frame(france_1800_barycentres, stringsAsFactors = FALSE) %>% tbl_df()
rownames(france_1800_barycentres) <- NULL
france_1800_barycentres <-
  france_1800_barycentres %>%
  mutate(long = as.numeric(x), lat = as.numeric(y)) %>%
  select(-x,-y)

La carte de base est la suivante :

p_map <- 
  ggplot() +
  geom_polygon(data = esperance_vie_map_ic %>% 
                 gather(sexe, esperance, e_f, e_h) %>% 
                 mutate(sexe = factor(sexe, levels = c("e_f", "e_h"), labels = c("Femmes", "Hommes"))),
               aes(x = long, y = lat, group = group, fill = esperance)) +
  facet_grid(~sexe) +
  scale_x_continuous(limits = c(-6, 10)) +
  scale_y_continuous(limits = c(41, 51.5)) +
  scale_fill_gradient("Espérance de vie à la naissance (années)", low = "yellow", high = "red") +
  coord_quickmap() +
  theme(legend.position = "bottom", text = element_text(size = 6),
        legend.key.height   = unit(1, "line"),
        legend.key.width    = unit(1.5, "line"))

Nous rajoutons les intervalles de confiance :

val_longueur <- 1
val_corresp <- .1

longueur_ic <- 
  ex_depts_avg %>% 
  left_join(france_1800_barycentres, by = c("dept" = "departement")) %>% 
  mutate(longueur_ic_tx = ex_sup - ex_inf) %>% 
  mutate(
    lat_beg = lat - (val_corresp * longueur_ic_tx/val_longueur)/2,
    lat_2 = lat + (val_corresp * longueur_ic_tx/val_longueur)/2) %>% 
  mutate(sexe = factor(sexe, levels = c("F", "M"), labels = c("Femmes", "Hommes")))

legende_x <- -4.5
legende_y <- 43

legende_min <- 5
legende_max <- 10

p_map <- 
  p_map +
  geom_errorbar(data = longueur_ic,
                aes(x = long, ymin = lat_beg, ymax = lat_2), colour = "black", size = 1.5, width=0.2) +
  geom_errorbar(data = longueur_ic,
                aes(x = long, ymin = lat_beg, ymax = lat_2), colour = "white", size = 1, width=0.1) +
  geom_errorbar(data = data.frame(long = legende_x, lat = legende_y, longueur_ic_tx = c(legende_min, legende_max)) %>%
                  mutate(lat_2 = lat + val_corresp * longueur_ic_tx/val_longueur),
                aes(x = long, ymin = lat, ymax = lat_2), colour = "black", size = 1.5, width=0.2) +
  geom_errorbar(data = data.frame(long = legende_x, lat = legende_y, longueur_ic_tx = c(legende_min, legende_max)) %>%
                  mutate(lat_2 = lat + val_corresp * longueur_ic_tx/val_longueur),
                aes(x = long, ymin = lat, ymax = lat_2), colour = "white", size = 1, width=0.1) +
  geom_text(data = data.frame(long = legende_x + .5, lat = legende_y, longueur_ic_tx = c(legende_min, legende_max)) %>% 
              mutate(lat_2 = lat + val_corresp * longueur_ic_tx/val_longueur),
            aes(x = long, y = lat_2, label = longueur_ic_tx))

p_map

4 Dimension spatiale des données : sédentarité et migration intérieure en France

Nous pouvons explorer les aspects spatiaux qu’offrent les données de généalogie dans plus de détail. On commence par charger quelques packages, ainsi que des cartes de France et les données de généalogie spécifiques au sous-échantillon d’individus susceptibles d’être en âge de migrer.

library(tidyverse)
library(stringr)
library(data.table)
library(pbapply)
library(dplyrExtras)
library(sp)

load("liste_depts.rda")
load("france.rda")
load("france_1800.rda")

liste_depts <- 
  liste_depts %>% 
  mutate(dept2 = dept,
         dept2 = ifelse(dept2 %in% c("F78", "F91", "F92", "F95"), yes = "F78", no = dept2),
         dept2 = ifelse(dept2 %in% c("F93", "F94", "F77"), yes = "F77", no = dept2)
  )

# Les données de généalogie
load("generations_mobilite.rda")
generations <- generations_mobilite
rm(generations_mobilite)

Nous allons avoir besoin de définir dans quel département certains points repérés par leurs coordonnées (longitude,latitude) se trouvent. Nous redéfinissons deux fonctions permettant de transformer des tableaux de données contenant les contours des départements en polygones.

library(rgeos)
#' group_to_gpc_poly
#' Returns the polygon from a group within the data frame containing
#' the coordinates of regions
#' @df: data frame with the coordinates of the polygons
#' @group: (string) group ID
group_to_gpc_poly <- function(df, group){
  group_coords <- df[df$group == group, c("long", "lat")]
  as(group_coords, "gpc.poly")
}

#' departement_polygone
#' Returns a French region as a gpc.poly object
#' un_dep: (string) code of the region
departement_polygone <- function(un_dep){
  dep_tmp <-
    france_1800 %>%
    filter(departement == un_dep)
  df_tmp_polys <- lapply(unique(dep_tmp$group), group_to_gpc_poly, df = dep_tmp)
  df_tmp_polys_tot <- df_tmp_polys[[1]]
  if(length(df_tmp_polys)>1){
    for(i in 2:length(df_tmp_polys)){
      df_tmp_polys_tot <- gpclib::union(df_tmp_polys_tot,  df_tmp_polys[[i]])
    }
  }
  df_tmp_polys_tot
}# End of departement_polygone()

france_1800_gpc <- lapply(unique(france_1800$departement), departement_polygone)
names(france_1800_gpc) <- unique(france_1800$departement)
save(france_1800_gpc, file = "france_1800_gpc.rda")

Nous créons une fonction s’appuyant sur une autre, à savoir inside.gpc.poly() du package surveillance, pour déterminer, à partir de vecteurs de coordonnées (longitude, latitude) et d’un département sous forme de polygone, quels sont les points à l’intérieur du département.

est_dans_poly_naissance <- function(longitude, latitude, departement){
  surveillance::inside.gpc.poly(x = longitude,
                                y = latitude,
                                departement,
                                mode.checked = FALSE)
}

generations %>% 
  filter(dist_enfants_km >0) %>% 
  .$dist_enfants_km %>% 
  summary()

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
    0.004     6.452    23.126   201.593   219.673 16895.501 

generations %>% 
  filter(dist_p_enfants_km >0) %>% 
  .$dist_p_enfants_km %>% 
  summary()

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
    0.020     6.523    27.075   269.284   242.784 19190.069 

generations %>% 
  filter(dist_ap_enfants_km >0) %>% 
  .$dist_ap_enfants_km %>% 
  summary()

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
    0.031     7.881    40.540   316.859   273.826 19190.069 

4.1 Migrations de courtes et longues distances

Nous définissons deux variables relatives à la distance entre les aïeux et leurs descendants : une indiquant si la distance séparant les lieux de naissance entre un aïeul et ses descendants est inférieure ou non à 10 km, et une autre pour savoir si cette distance est inférieure ou non à 20km.

generations <- 
  generations %>% 
  mutate(moins_10_enfants = dist_enfants_km <= 10,
         moins_20_enfants = dist_enfants_km <= 20,
         moins_10_p_enfants = dist_p_enfants_km <= 10,
         moins_20_p_enfants = dist_p_enfants_km <= 20,
         moins_10_ap_enfants = dist_ap_enfants_km <= 10,
         moins_20_ap_enfants = dist_ap_enfants_km <= 20)
generations <- generations %>% data.table()

Nous regardons, pour chaque génération, le pourcentage de descendants nés à moins de 10 ou 20km du lieu de naissance de leur aïeul.

tableau_distances_1 <- 
  generations %>% 
  select(id_aieul, id_enfant, moins_10_enfants, moins_20_enfants, dist_enfants_km) %>% 
  rename(moins_10 = moins_10_enfants,
         moins_20 = moins_20_enfants,
         dist_km = dist_enfants_km) %>% 
  unique() %>% 
  summarise(
    moins_10_not_na = sum(!is.na(moins_10)),
    moins_10 = sum(moins_10, na.rm=T),
    moins_20_not_na = sum(!is.na(moins_20)),
    moins_20 = sum(moins_20, na.rm=T),
    nas = sum(is.na(dist_km), na.rm=T),
    nb = n()
  ) %>% 
  mutate(generation = 1,
         moins_10_pct = (moins_10 / moins_10_not_na * 100) %>% round(2),
         moins_20_pct = (moins_20 / moins_20_not_na * 100) %>% round(2))

tableau_distances_2 <- 
  generations %>% 
  select(id_aieul, id_p_enfant, moins_10_p_enfants, moins_20_p_enfants, dist_p_enfants_km) %>% 
  rename(moins_10 = moins_10_p_enfants,
         moins_20 = moins_20_p_enfants,
         dist_km = dist_p_enfants_km) %>% 
  unique() %>% 
  summarise(
    moins_10_not_na = sum(!is.na(moins_10)),
    moins_10 = sum(moins_10, na.rm=T),
    moins_20_not_na = sum(!is.na(moins_20)),
    moins_20 = sum(moins_20, na.rm=T),
    nas = sum(is.na(dist_km), na.rm=T),
    nb = n()
  ) %>% 
  mutate(generation = 2,
         moins_10_pct = (moins_10 / moins_10_not_na * 100) %>% round(2),
         moins_20_pct = (moins_20 / moins_20_not_na * 100) %>% round(2))

tableau_distances_3 <- 
  generations %>% 
  select(id_aieul, id_ap_enfant, moins_10_ap_enfants, moins_20_ap_enfants, dist_ap_enfants_km) %>% 
  rename(moins_10 = moins_10_ap_enfants,
         moins_20 = moins_20_ap_enfants,
         dist_km = dist_ap_enfants_km) %>% 
  unique() %>% 
  summarise(
    moins_10_not_na = sum(!is.na(moins_10)),
    moins_10 = sum(moins_10, na.rm=T),
    moins_20_not_na = sum(!is.na(moins_20)),
    moins_20 = sum(moins_20, na.rm=T),
    nas = sum(is.na(dist_km), na.rm=T),
    nb = n()
  ) %>% 
  mutate(generation = 3,
         moins_10_pct = (moins_10 / moins_10_not_na * 100) %>% round(2),
         moins_20_pct = (moins_20 / moins_20_not_na * 100) %>% round(2))

On agrègr ces trois tableaux en un seul :

tableau_distances <- 
  tableau_distances_1 %>% 
  bind_rows(tableau_distances_2) %>% 
  bind_rows(tableau_distances_3)
tableau_distances

   moins_10_not_na moins_10 moins_20_not_na moins_20  nas    nb generation
1:           26384    20596           26384    21754 1132 27516          1
2:           32793    19265           32793    21899 1705 34498          2
3:           66878    29065           66878    34990 7038 73916          3
   moins_10_pct moins_20_pct
1:        78.06        82.45
2:        58.75        66.78
3:        43.46        52.32

4.1.1 Distribution de la distance séparant les lieux de naissance des descendants avec ceux des aïeux

Pour avoir une idée de la distribution des distances séparant les lieux de naissances des descendants de celui de leur aïeul, nous préparons des tableaux, pour chaque génération, indiquant pour chaque individu les distances en kilomètres, le sexe, ainsi que la génération à laquelle chaque individu appartient.

enfants <- 
  generations %>% 
  select(id_enfant, long_N_enfants, lat_N_enfants, sexe_enfants,
         dept_N_enfants, moins_10_enfants, moins_20_enfants,
         dist_enfants_km) %>% 
  mutate(generation = 1) %>% 
  unique()

colnames(enfants) <- str_replace(colnames(enfants), "\\_enfants?", "")
enfants <- data.table(enfants)

p_enfants <- 
  generations %>% 
  select(id_p_enfant, long_N_p_enfants, lat_N_p_enfants, sexe_p_enfants,
         dept_N_p_enfants, moins_10_p_enfants, moins_20_p_enfants,
         dist_p_enfants_km) %>% 
  mutate(generation = 2) %>% 
  unique()

colnames(p_enfants) <- str_replace(colnames(p_enfants), "\\_p_enfants?", "")
p_enfants <- data.table(p_enfants)

ap_enfants <- 
  generations %>% 
  select(id_ap_enfant, long_N_ap_enfants, lat_N_ap_enfants, sexe_ap_enfants,
         dept_N_ap_enfants, moins_10_ap_enfants, moins_20_ap_enfants,
         dist_ap_enfants_km) %>% 
  mutate(generation = 3) %>% 
  unique()

colnames(ap_enfants) <- str_replace(colnames(ap_enfants), "\\_ap_enfants?", "")
ap_enfants <- data.table(ap_enfants)

Nous fusionnons ces trois tableaux dans un seul :

generations_2 <- 
  enfants %>% 
  bind_rows(p_enfants) %>% 
  bind_rows(ap_enfants)

Nous retirons les valeurs manquantes. Par ailleurs, nous créons une variable de distance exprimée en logarithme. Cela permet d’’observer d’avoir une idée de la distribution des distances pour ceux qui migrent, les sédentaires étant présents dans une écrasante majorité.

generations_2 <- generations_2[!is.na(dist_km)]
generations_2 <- generations_2[, dist_km_2 := dist_km+1]
generations_2 <- generations_2[, dist_km_log := log(dist_km+1)]

Nous mettons en forme le tableau de données pour que ces dernières soient lisibles par les fonctions du package ggplot2.

stat_des_distances <- 
  generations_2 %>% 
  select(generation, dist_km) %>% 
  tbl_df() %>% 
  mutate(type = ifelse(dist_km == 0, yes = "Sédentaire", no = NA),
         type = ifelse(dist_km<=20 & dist_km >0, yes = "Courte", no = type),
         type = ifelse(dist_km>20 , yes = "Longue", no = type))

Un aperçu de la répartition entre sédentaires, migrants de courte distance et migrants de longue distance s’obtient à l’aide des instructions suivantes :

stat_des_distances %>% 
  group_by(type) %>% 
  summarise(nb_obs = n()) %>% 
  mutate(pct = nb_obs / sum(nb_obs)*100) %>% 
  knitr::kable(digits=2)

type	nb_obs	pct
Courte	30207	25.62
Longue	45496	38.59
Sédentaire	42203	35.79

Et en différentiant en fonction de la génération (en colonnes) :

stat_des_distances %>% 
  filter(!is.na(type)) %>% 
  group_by(generation, type) %>% 
  summarise(nb_obs = n()) %>% 
  group_by(generation) %>% 
  mutate(pct = nb_obs / sum(nb_obs) * 100) %>% 
  select(-nb_obs) %>% 
  spread(generation, pct) %>% 
  knitr::kable(digits=2)

type	1	2	3
Courte	19.43	27.70	27.06
Longue	18.40	34.24	48.77
Sédentaire	62.17	38.06	24.17

Nous voulons visualiser la distribution des distances, par génération. Dans un premier temps, nous nous appuyons sur les paramètres d’échelle utilisés lors de l’appel à la fonction hist() du package graphics.

hist_dist <- hist(generations_2$dist_km_log, plot=FALSE)

Dans un second temps, nous créons une fonction retournant pour chacune des barres de l’histogramme, les valeurs en abscisse et en ordonnées, ces première reflétant les distances, ces secondes représentant le pourcentage d’observations pour chaque distance donnée.

#' hist_gen
#' Valeurs pour l'histogramme des distances, pour une generation
#' @une_gen (int) : numéro de la génération
hist_gen <- function(une_gen){
  hist_tmp <- 
    generations_2 %>% filter(generation == une_gen) %>% 
    magrittr::extract2("dist_km_log") %>% 
    hist(plot=F, breaks = hist_dist$breaks)
  
  data.frame(x = hist_tmp$mids, counts = hist_tmp$counts) %>% 
    mutate(pct = counts / sum(counts),
           generation = une_gen) %>% 
    tbl_df()
}# Fin de hist_gen()

Nous appliquons cette fonction à chacune des trois générations à notre disposition.

distances_gen_hist <- 
  lapply(unique(generations_2$generation), hist_gen) %>% 
  bind_rows()

Enfin, nous pouvons créer et afficher les histogrammes.

p_dist_histo <- 
  ggplot(data = distances_gen_hist %>% 
           mutate(generation = factor(generation,
                                      levels = c(0, 1, 2, 3),
                                      labels = c('Aïeux', "Enfants", "Petits-enfants", "Arrière-petits-enfants"))),
       aes(x = exp(x), y = pct)) +
  geom_bar(stat = "identity") +
  facet_wrap(~generation, nrow = 2, scales = "free_x") +
  scale_x_log10(breaks=c(0, 1, 2, 5, 10, 25, 50, 100, 500, 2500, 10000)) +
  ylab(NULL) + xlab("Distance (km)")

p_dist_histo

Nous pouvons aussi nous intéresser aux déplacements intergénérationnels en France en fonction du genre. Pour ce faire, nous proposons de regarder cela sur un tableau indiquant les distances en lignes, les générations et le sexe en colonne.

# Définition des classes de distances
cut_dist <- c(0:25, 50, 75, 100, 150, 200, +Inf)
cut(seq(0, 5), cut_dist, ordered_result = TRUE, include.lowest = TRUE)

[1] [0,1] [0,1] (1,2] (2,3] (3,4] (4,5]
31 Levels: [0,1] < (1,2] < (2,3] < (3,4] < (4,5] < (5,6] < ... < (200,Inf]

generations_dist_table <- 
  generations_2 %>% 
  filter(!is.na(dist_km), !is.na(sexe), !is.na(generation)) %>% 
  mutate(dist_km_rounded = round(dist_km)) %>% 
  select(generation, sexe, dist_km_rounded) %>% 
  mutate(dist_km_categ = cut(dist_km_rounded, cut_dist, ordered_result = TRUE, include.lowest=TRUE)) %>% 
  tbl_df() %>% 
  group_by(generation, sexe, dist_km_categ) %>% 
  tally() %>% 
  group_by(generation, sexe) %>% 
  mutate(pct = n/sum(n) * 100,
         pct_cumul = cumsum(pct)) %>% 
  mutate(pct = round(pct, 2),
         pct_cumul = round(pct_cumul, 2))

generations_dist_table_tex <- 
  generations_dist_table %>% 
  select(generation, sexe, dist_km_categ, pct, pct_cumul) %>% 
  ungroup() %>% 
  mutate(sexe = factor(sexe, levels = c(2,1), labels = c("Femmes", "Hommes"))) %>% 
  mutate(value = str_c(pct, " (", pct_cumul, ")")) %>% 
  select(-pct, -pct_cumul) %>% 
  dcast(., dist_km_categ~generation + sexe, value.var = "value")

generations_dist_table_tex %>% 
  knitr::kable(digits=2)

dist_km_categ	1_Femmes	1_Hommes	2_Femmes	2_Hommes	3_Femmes	3_Hommes
[0,1]	63.74 (63.74)	62.24 (62.24)	39.07 (39.07)	38.54 (38.54)	25.21 (25.21)	24.32 (24.32)
(1,2]	1.82 (65.55)	1.55 (63.78)	1.88 (40.94)	1.88 (40.42)	1.82 (27.03)	1.64 (25.96)
(2,3]	2.15 (67.7)	2.26 (66.04)	2.71 (43.66)	2.66 (43.08)	2.26 (29.29)	2.2 (28.16)
(3,4]	2.41 (70.11)	2.15 (68.19)	2.89 (46.55)	3.19 (46.27)	2.66 (31.95)	2.72 (30.88)
(4,5]	2.08 (72.19)	2.1 (70.28)	2.67 (49.21)	2.7 (48.97)	2.39 (34.33)	2.59 (33.47)
(5,6]	1.51 (73.7)	1.88 (72.16)	2.57 (51.78)	2.39 (51.35)	2.34 (36.67)	2.38 (35.85)
(6,7]	1.37 (75.07)	1.5 (73.66)	1.91 (53.69)	2.18 (53.54)	2.05 (38.72)	2.04 (37.89)
(7,8]	1.25 (76.31)	1.24 (74.9)	1.61 (55.29)	1.75 (55.28)	1.8 (40.52)	1.62 (39.52)
(8,9]	0.92 (77.23)	1.09 (75.99)	1.46 (56.76)	1.54 (56.82)	1.6 (42.12)	1.55 (41.07)
(9,10]	0.82 (78.05)	1.11 (77.1)	1.3 (58.05)	1.53 (58.36)	1.47 (43.59)	1.29 (42.36)
(10,11]	0.64 (78.69)	0.78 (77.88)	1.16 (59.21)	1.14 (59.5)	1.29 (44.88)	1.34 (43.7)
(11,12]	0.69 (79.39)	0.49 (78.36)	0.91 (60.12)	0.97 (60.47)	1.11 (45.99)	1.01 (44.71)
(12,13]	0.49 (79.87)	0.52 (78.88)	0.98 (61.1)	0.93 (61.4)	1.14 (47.13)	1.03 (45.74)
(13,14]	0.35 (80.22)	0.42 (79.3)	0.8 (61.9)	0.95 (62.35)	0.84 (47.97)	0.89 (46.63)
(14,15]	0.27 (80.49)	0.44 (79.74)	0.67 (62.57)	0.79 (63.14)	0.85 (48.82)	0.81 (47.44)
(15,16]	0.43 (80.92)	0.39 (80.13)	0.8 (63.37)	0.6 (63.74)	0.89 (49.71)	0.78 (48.22)
(16,17]	0.33 (81.25)	0.38 (80.51)	0.77 (64.14)	0.73 (64.47)	0.69 (50.4)	0.68 (48.9)
(17,18]	0.35 (81.6)	0.32 (80.83)	0.57 (64.71)	0.63 (65.11)	0.66 (51.07)	0.73 (49.63)
(18,19]	0.26 (81.86)	0.31 (81.14)	0.5 (65.21)	0.67 (65.78)	0.65 (51.71)	0.69 (50.31)
(19,20]	0.27 (82.13)	0.32 (81.45)	0.49 (65.69)	0.5 (66.28)	0.64 (52.35)	0.57 (50.88)
(20,21]	0.26 (82.39)	0.22 (81.67)	0.45 (66.14)	0.43 (66.71)	0.5 (52.85)	0.55 (51.43)
(21,22]	0.24 (82.64)	0.25 (81.92)	0.46 (66.6)	0.35 (67.06)	0.44 (53.29)	0.52 (51.95)
(22,23]	0.24 (82.88)	0.28 (82.19)	0.34 (66.94)	0.34 (67.4)	0.41 (53.7)	0.54 (52.49)
(23,24]	0.16 (83.04)	0.29 (82.48)	0.46 (67.4)	0.33 (67.73)	0.43 (54.13)	0.4 (52.89)
(24,25]	0.09 (83.13)	0.19 (82.67)	0.28 (67.67)	0.25 (67.98)	0.35 (54.48)	0.32 (53.22)
(25,50]	3 (86.13)	2.51 (85.19)	5.36 (73.03)	5.07 (73.05)	6.97 (61.46)	7.07 (60.29)
(50,75]	1.27 (87.41)	1.57 (86.75)	2.7 (75.74)	2.84 (75.89)	3.95 (65.41)	3.68 (63.97)
(75,100]	1.11 (88.51)	1 (87.75)	1.82 (77.56)	1.84 (77.72)	2.71 (68.11)	2.63 (66.61)
(100,150]	1.59 (90.1)	1.72 (89.47)	2.78 (80.34)	2.83 (80.55)	4.08 (72.19)	4.23 (70.83)
(150,200]	1.34 (91.44)	1.48 (90.95)	2.49 (82.83)	2.82 (83.38)	3.5 (75.7)	3.75 (74.58)
(200,Inf]	8.56 (100)	9.05 (100)	17.17 (100)	16.62 (100)	24.3 (100)	25.42 (100)

4.2 Transitions entre petites et grandes villes

Nous allors faire un point sur la migration entre les villes, selon le type de ces dernières. La Statistique Générale de la France (SGF) fournit des statistiques pour certaines communes : https://www.insee.fr/fr/statistiques/2591293?sommaire=2591397. Nous décidons de séparer les communes en plusieurs catégories, en reprenant le principe de Fleury & Henry (1958) et Blayo (1975), qui consiste à segmenter les communes en fonction de leur nombre d’habitants. Nous utilisons les statistiques de population de 1836 pour définir, parmi les communes présentes dans les données de la SGF :

• les grandes villes (plus de 50 000 habitants) ;
• les villes moyennes (entre 10 000 et 50 000 habitants) ;
• les petites villes (moins de 10 000 habitants).

Nous considérons que les communes n’étant pas présentes dans les données de la SGF sont à classer dans une catégorie à part, que nous appelons “très petites villes”.

Nous chargeons les données de la SGF :

library(stringi)
pop_chef_lieux <- readxl::read_excel("../data/population_insee/TERR_T116.xls", skip=7)

pop_chef_lieux <- 
  pop_chef_lieux %>% 
  filter(`Code de l'unité géographique` == 3) %>% 
  select(`Code du département`, `Nom de l'unité d'analyse`,
         `Population totale des villes chefs-lieux d'arrondissement en 1836`) %>% 
  rename(code_dep = `Code du département`, nom_insee = `Nom de l'unité d'analyse`,
         pop_1836 = `Population totale des villes chefs-lieux d'arrondissement en 1836`) %>% 
  mutate(nom = stri_replace_all_regex(nom_insee, "\\((.*?)\\)", "") %>% 
  stri_trans_general(., "Latin-ASCII") %>% 
  stri_replace_all_regex(., "'|\\(|\\)|[[:blank:]]|\\?|_|\\*", "") %>% 
  stri_trim(.) %>% 
  str_to_lower()) %>% 
  mutate(code_dep = ifelse(code_dep %in% c("78", "91", "92", "95"), yes = "78", no = code_dep),
         code_dep = ifelse(code_dep %in% c("93", "94", "77"), yes = "77", no = code_dep),
         code_dep = ifelse(code_dep %in% c("2A", "2B"), yes = "20", no = code_dep))

head(pop_chef_lieux)

# A tibble: 6 x 4
  code_dep nom_insee       pop_1836 nom            
  <chr>    <chr>              <dbl> <chr>          
1 01       BELLEY              3970 belley         
2 01       BOURG               9528 bourg          
3 01       NANTUA              3696 nantua         
4 01       TREVOUX             2559 trevoux        
5 02       CHATEAU-THIERRY     4761 chateau-thierry
6 02       LAON                8230 laon           

Nous rajoutons Nice à la main…

pop_chef_lieux <- 
  pop_chef_lieux %>% 
  bind_rows(data.frame(code_dep = "06", nom_insee = "NICE", pop_1836 = 33811, nom = "nice", stringsAsFactors = FALSE))

Nous allons ajouter des données complémentaires, qui pourraient par la suite être utilisées, notamment la surface des communes. Ces données sont proposées par Open Street Map (OSM), une annexe permet de comprendre comment nous les obtenons.

load("../shp_OSM/communes_sans_cercle.rda")
communes_osm <- lapply(communes_sans_cercle, function(x) x@data) %>% bind_rows() %>% tbl_df()
communes_osm <- 
  communes_osm %>%
  mutate(code_dep = str_sub(insee, 1, 2)) %>% 
  rename(nom_osm = nom) %>% 
  mutate(nom = stri_replace_all_regex(nom_osm, "\\((.*?)\\)", "") %>% 
         stri_trans_general(., "Latin-ASCII") %>% 
         stri_replace_all_regex(., "'|\\(|\\)|[[:blank:]]|\\?|_|\\*", "") %>% 
         stri_trim(.) %>% 
         str_to_lower()) %>% 
  mutate(code_dep = ifelse(code_dep %in% c("78", "91", "92", "95"), yes = "78", no = code_dep),
         code_dep = ifelse(code_dep %in% c("93", "94", "77"), yes = "77", no = code_dep),
         code_dep = ifelse(code_dep %in% c("2A", "2B"), yes = "20", no = code_dep))

head(communes_osm)

Le point épineux ici consiste à apparier les bases de la SGF avec celles d’OSM. Le plus simple serait d’utiliser les codes INSEE, mais ils ne sont pas fournis par la SGF. Nous décidons alors d’utiliser le code département et le nom de la commune. Cela dit, l’appariement n’est pas total, certaines orthographes diffèrent selon la source des données.

pop_chef_lieux %>%
    anti_join(communes_osm) %>% 
  head()

# A tibble: 6 x 4
  code_dep nom_insee          pop_1836 nom               
  <chr>    <chr>                 <dbl> <chr>             
1 01       BOURG                  9528 bourg             
2 03       MOULINS-SUR-ALLIER    15231 moulins-sur-allier
3 04       DIGNE                  6365 digne             
4 07       TOURNON                4174 tournon           
5 08       MEZIERES               4083 mezieres          
6 08       ROCROY                 3682 rocroy            

Nous créons alors une fonction pratique pour repérer les potentiels appariements possibles.

chercher_match_commune <- function(i){
  (tmp <- pop_chef_lieux %>% 
     anti_join(communes_osm) %>% 
     slice(i)
  )
  
  cat(str_c("Je cherche pour : ", tmp$nom, "\n"))
  print(tmp)
  
  res_tmp <-
    communes_osm %>% 
    filter(code_dep == tmp$code_dep) %>% 
    filter(str_sub(nom, 1, 4) == str_sub(tmp$nom, 1, 4)) %>% 
    arrange(nom)
  
  if(nrow(res_tmp) == 0){
    print("Pas trouvé... Je cherche plus loin\n")
    res_tmp <- communes_osm %>% 
      filter(code_dep == tmp$code_dep) %>% 
      filter(str_detect(nom, tmp$nom)) %>% 
      arrange(nom)
  }
  res_tmp
}

L’utilisation de cette fonction est simple :

chercher_match_commune(1)

Je cherche pour : bourg
# A tibble: 1 x 4
  code_dep nom_insee pop_1836 nom  
  <chr>    <chr>        <dbl> <chr>
1 01       BOURG         9528 bourg

# A tibble: 2 x 6
  insee nom_osm                wikipedia       surf_ha code_dep nom       
  <chr> <chr>                  <chr>             <dbl> <chr>    <chr>     
1 01053 Bourg-en-Bresse        fr:Bourg-en-Br…    2393 01       bourg-en-…
2 01054 Bourg-Saint-Christophe fr:Bourg-Saint…     900 01       bourg-sai…

Nous effectuons les corrections suivantes, en utilisant successivement la fonction chercher_match_commune() :

a_remplacer_noms <- 
  c("bourg", "bourg-en-bresse",
    "moulins-sur-allier", "moulins",
    "digne", "digne-les-bains",
    "tournon", "tournon-sur-rhone",
    "mezieres", "charleville-mezieres",
    "rocroy", "rocroi",
    "milhau", "millau",
    "aix", "aix-en-provence",
    "arles-sur-rhone", "arles",
    "vire", "virenormandie",
    "barbezieux", "barbezieux-saint-hilaire",
    "st-jean-dangely", "saint-jean-dangely",
    "st-amand-montrond", "saint-amand-montrond",
    "brive", "brive-la-gaillarde",
    "semur", "semur-en-auxois",
    "guincamp", "guingamp",
    "sarlat", "sarlat-la-caneda",
    "beaume-les-dames", "baume-les-dames",
    "montelimart", "montelimar",
    "alais", "ales",
    "lesparre", "lesparre-medoc",
    "saint-pons", "saint-pons-de-thomieres",
    "montfort", "montfort-sur-meu",
    "romorantin", "romorantin-lanthenay",
    "lepuy", "lepuy-en-velay",
    "yssengeaux", "yssingeaux",
    "florac", "florac-trois-rivieres",
    "bauge", "bauge-en-anjou",
    "beaupreau", "beaupreau-en-mauges",
    "segre", "segre-en-anjou-bleu",
    "cherbourg", "cherbourg-en-cotentin",
    "mortain", "mortain-bocage",
    "chalons-sur-marne", "chalons-en-champagne",
    "vitry", "vitry-le-francois",
    "vassy", "wassy",
    # "briey", "valdebriey",
    "cosne", "cosne-cours-sur-loire",
    "avesnes", "avesnes-sur-helpe",
    "domfront", "domfront-en-poiraie",
    "mortagne", "mortagne-au-perche",
    "montreuil-sur-mer", "montreuil",
    "saint-pol", "saint-pol-sur-ternoise",
    "mauleon", "mauleon-licharre",
    "oloron", "oloron-sainte-marie",
    "argeles", "argeles-gazost",
    "bagneres", "bagneres-de-bigorre",
    "schelestadt", "selestat",
    "mantes", "mantes-la-jolie",
    "castel-sarrazin", "castelsarrasin",
    "saint-yrieix", "saint-yrieix-la-perche",
    "saint-die", "saint-die-des-vosges",
    "corbeil", "corbeil-essonnes",
    "briey", "valdebriey"
    ) %>% 
  matrix(byrow=TRUE, ncol=2) %>% 
  data.frame(stringsAsFactors = FALSE)

colnames(a_remplacer_noms) <- c("nom", "nouveau_nom")

Et nous répercutons ces modifications sur pop_chef_lieux :

pop_chef_lieux <- 
  pop_chef_lieux %>% 
  left_join(a_remplacer_noms) %>% 
  mutate(nom = ifelse(!is.na(nouveau_nom), yes = nouveau_nom, no = nom)) %>% 
  select(-nouveau_nom)

Nous modifions également les informations relatives à certains départements créés récemment.

a_remplacer_noms_dep <- 
  c("briey", "54",
    "saint-denis", "77",
    "sceaux", "78",
    "corbeil", "78",
    "etampes", "78",
    "pontoise", "78",
    "grasse", "06",
    "chateau-salins", "57",
    "luneville", "54",
    "nancy", "54",
    "sarrebourg", "57",
    "toul", "54") %>% 
  matrix(byrow=TRUE, ncol=2) %>% 
  data.frame(stringsAsFactors = FALSE)

colnames(a_remplacer_noms_dep) <- c("nom", "nouveau_dep")

pop_chef_lieux <- 
  pop_chef_lieux %>% 
  left_join(a_remplacer_noms_dep) %>% 
  mutate(code_dep = ifelse(!is.na(nouveau_dep), yes = nouveau_dep, no = code_dep)) %>% 
  select(-nouveau_dep)

L’appariement entre les données d’OSM et de la SGF peut ainsi s’effectuer :

villes_francaises <- 
  pop_chef_lieux %>% 
  left_join(communes_osm) %>% 
  filter(!is.na(pop_1836))

villes_francaises %>% arrange(-pop_1836) %>% data.frame() %>% slice(1:12)

# A tibble: 12 x 8
   code_dep nom_insee  pop_1836 nom        insee nom_osm wikipedia surf_ha
   <chr>    <chr>         <dbl> <chr>      <chr> <chr>   <chr>       <dbl>
 1 75       PARIS        909126 paris      75056 Paris   fr:Paris    10537
 2 69       LYON         150814 lyon       69123 Lyon    fr:Lyon      4798
 3 13       MARSEILLE    146239 marseille  13055 Marsei… fr:Marse…   24189
 4 33       BORDEAUX      98705 bordeaux   33063 Bordea… fr:Borde…    4984
 5 76       ROUEN         92083 rouen      76540 Rouen   fr:Rouen     2139
 6 31       TOULOUSE      77372 toulouse   31555 Toulou… fr:Toulo…   11802
 7 44       NANTES        75895 nantes     44109 Nantes  fr:Nantes    6579
 8 59       LILLE         72005 lille      59350 Lille   fr:Lille     3483
 9 67       STRASBOURG    57885 strasbourg 67482 Strasb… fr:Stras…    7828
10 80       AMIENS        46129 amiens     80021 Amiens  fr:Amiens    5007
11 30       NIMES         43036 nimes      30189 Nîmes   fr:Nîmes    16170
12 57       METZ          42793 metz       57463 Metz    fr:Metz      4176

Nous calculons la densité de population (nombre d’habitants au kilomètre carré), et définissons les trois types de villes parmi les communes citées par la SGF.

villes_francaises <- 
  villes_francaises %>% 
  # Densite au km2
  mutate(densite_pop_1836 = pop_1836 / (surf_ha * 0.01)) %>% 
  mutate(type_commune = "petite",
         type_commune = ifelse(pop_1836 >= 10000, yes = "moyenne", no = type_commune),
         type_commune = ifelse(pop_1836 >= 50000, yes = "grande", no = type_commune)
         ) %>% 
  filter(!is.na(insee))

Nous utilisons cette classification pour segmenter les communes de naissance des individus composant notre échantillon :

# Aïeux
types_communes_aieul <- 
  villes_francaises %>% 
  select(insee, type_commune) %>% 
  rename(code_insee_N_aieul = insee, type_commune_aieul = type_commune) %>% 
  data.table()

setkey(types_communes_aieul, code_insee_N_aieul)
setkey(generations, code_insee_N_aieul)
generations <- types_communes_aieul[generations]

# Enfants
types_communes_enfants <- 
  villes_francaises %>% 
  select(insee, type_commune) %>% 
  rename(code_insee_N_enfants = insee, type_commune_enfants = type_commune) %>% 
  data.table()

setkey(types_communes_enfants, code_insee_N_enfants)
setkey(generations, code_insee_N_enfants)
generations <- types_communes_enfants[generations]

# Petits-enfants
types_communes_p_enfants <- 
  villes_francaises %>% 
  select(insee, type_commune) %>% 
  rename(code_insee_N_p_enfants = insee, type_commune_p_enfants = type_commune) %>% 
  data.table()

setkey(types_communes_p_enfants, code_insee_N_p_enfants)
setkey(generations, code_insee_N_p_enfants)
generations <- types_communes_p_enfants[generations]

# Arrière-petits-enfants
types_communes_ap_enfants <- 
  villes_francaises %>% 
  select(insee, type_commune) %>% 
  rename(code_insee_N_ap_enfants = insee, type_commune_ap_enfants = type_commune)%>% 
  data.table()


setkey(types_communes_ap_enfants, code_insee_N_ap_enfants)
setkey(generations, code_insee_N_ap_enfants)
generations <- types_communes_ap_enfants[generations]

Ensuite, nous rajoutons la modalité “très petite ville” aux communes n’étant pas citées dans les données de la SGF.

generations$type_commune_aieul[which(is.na(generations$type_commune_aieul))] <- "autres"
generations$type_commune_enfants[which(is.na(generations$type_commune_enfants))] <- "autres"
generations$type_commune_p_enfants[which(is.na(generations$type_commune_p_enfants))] <- "autres"
generations$type_commune_ap_enfants[which(is.na(generations$type_commune_ap_enfants))] <- "autres"

generations <- 
  generations %>% tbl_df() %>% 
  mutate(type_commune_aieul = factor(type_commune_aieul, levels = c("grande", "moyenne", "petite", "autres")),
         type_commune_enfants = factor(type_commune_enfants, levels = c("grande", "moyenne", "petite", "autres")),
         type_commune_p_enfants = factor(type_commune_p_enfants, levels = c("grande", "moyenne", "petite", "autres")),
         type_commune_ap_enfants = factor(type_commune_ap_enfants, levels = c("grande", "moyenne", "petite", "autres"))) %>% 
  data.table()

Nous pouvons dès lors commencer à sortir quelques statistiques. Nous regardons la proportion d’individus dans chaque type de villes, par génération :

prop_villes_gen <- 
  generations %>% 
  tbl_df() %>% 
  select(id_aieul, type_commune_aieul) %>% 
  unique() %>% 
  group_by(type_commune_aieul) %>% 
  tally() %>%
  mutate(generation = "0") %>% 
  mutate(pct = n / sum(n)) %>% 
  rename(type_commune = type_commune_aieul) %>% 
  bind_rows(
    generations %>% 
      tbl_df() %>% 
      select(id_enfant, type_commune_enfants) %>% 
      unique() %>% 
      group_by(type_commune_enfants) %>% 
      tally() %>%
      mutate(generation = "1") %>% 
      mutate(pct = n / sum(n)) %>% 
      rename(type_commune = type_commune_enfants)
  ) %>% 
  bind_rows(
    generations %>% 
      tbl_df() %>% 
      select(id_p_enfant, type_commune_p_enfants) %>% 
      unique() %>% 
      group_by(type_commune_p_enfants) %>% 
      tally() %>%
      mutate(generation = "2") %>% 
      mutate(pct = n / sum(n)) %>% 
      rename(type_commune = type_commune_p_enfants)
  ) %>% 
  bind_rows(
    generations %>% 
      tbl_df() %>% 
      select(id_ap_enfant, type_commune_ap_enfants) %>% 
      unique() %>% 
      group_by(type_commune_ap_enfants) %>% 
      tally() %>%
      mutate(generation = "3") %>% 
      mutate(pct = n / sum(n)) %>% 
      rename(type_commune = type_commune_ap_enfants)
  )

Le résultat est le suivant :

prop_villes_gen %>% select(-n) %>%  
  mutate(pct = pct*100) %>% spread(generation, pct) %>% 
  knitr::kable(digits=2)

type_commune	0	1	2	3
grande	3.40	3.92	6.31	7.63
moyenne	5.18	5.25	6.55	7.86
petite	3.43	3.56	3.73	3.81
autres	87.99	87.27	83.41	80.70

Nous regardons ensuite les déplacements transitoires conditionnels entre petites, moyennes et grandes villes, en pourcentages.

# Aieux --> enfants
villes_aieul_enfant <- 
  generations %>% 
  select(id_aieul, id_enfant, type_commune_aieul, type_commune_enfants) %>% 
  unique()
villes_aieul_enfant <- 
  xtabs(~type_commune_aieul + type_commune_enfants, villes_aieul_enfant)

# Aieux --> petits-enfants
villes_aieul_p_enfant <- 
  generations %>% 
  select(id_aieul, id_p_enfant, type_commune_aieul, type_commune_p_enfants) %>% 
  unique()
villes_aieul_p_enfant <- 
  xtabs(~type_commune_aieul + type_commune_p_enfants, villes_aieul_p_enfant)

# Aieux --> arrière-petits-enfants
villes_aieul_ap_enfant <- 
  generations %>% 
  select(id_aieul, id_ap_enfant, type_commune_aieul, type_commune_ap_enfants) %>% 
  unique()
villes_aieul_ap_enfant <- 
  xtabs(~type_commune_aieul + type_commune_ap_enfants, villes_aieul_ap_enfant)


table_villes_aieul_enfant <- ftable(villes_aieul_enfant)
table_villes_aieul_enfant <- table_villes_aieul_enfant/rowSums(table_villes_aieul_enfant) * 100

table_villes_aieul_p_enfant <- ftable(villes_aieul_p_enfant)
table_villes_aieul_p_enfant <- table_villes_aieul_p_enfant/rowSums(table_villes_aieul_p_enfant) * 100

table_villes_aieul_ap_enfant <- ftable(villes_aieul_ap_enfant)
table_villes_aieul_ap_enfant <- table_villes_aieul_ap_enfant/rowSums(table_villes_aieul_ap_enfant) * 100

Les pourcentages se lisent par ligne.

table_villes_aieul_enfant

                   type_commune_enfants     Paris  Couronne     Autre
type_commune_aieul                                                   
Paris                                   74.794521  3.561644 21.643836
Couronne                                 4.360465 79.651163 15.988372
Autre                                         NaN       NaN       NaN

table_villes_aieul_p_enfant

                   type_commune_p_enfants     Paris  Couronne     Autre
type_commune_aieul                                                     
Paris                                     55.040323  9.072581 35.887097
Couronne                                   6.766917 68.421053 24.812030
Autre                                           NaN       NaN       NaN

table_villes_aieul_ap_enfant

                   type_commune_ap_enfants     Paris  Couronne     Autre
type_commune_aieul                                                      
Paris                                      43.447038  8.976661 47.576302
Couronne                                   11.786600 58.188586 30.024814
Autre                                            NaN       NaN       NaN

4.3 Le cas parisien

Nous pouvons nous pencher davantage sur le cas parisien. Pour ce faire, nous créons une base comprenant les descendants nés dans le département de Paris ainsi que leurs aïeux. Les aïeux ne sont pas forcément nés dans le département de Paris, et c’est ce que nous allons étudier pour commencer.

enfants_75 <- 
  generations %>% 
  filter(dept_N_enfants == "F75") %>% 
  select(id_aieul, id_enfant, dept_N_aieul) %>% 
  mutate(generation = 1) %>% 
  unique() %>% 
  select(dept_N_aieul, generation)

p_enfants_75 <- 
  generations %>% 
  filter(dept_N_p_enfants == "F75") %>% 
  select(id_aieul, id_p_enfant, dept_N_aieul) %>% 
  mutate(generation = 2) %>% 
  unique() %>% 
  select(dept_N_aieul, generation)


ap_enfants_75 <- 
  generations %>% 
  filter(dept_N_ap_enfants == "F75") %>% 
  select(id_aieul, id_ap_enfant, dept_N_aieul) %>% 
  mutate(generation = 3) %>% 
  unique() %>% 
  select(dept_N_aieul, generation)

generations_paris <- 
  enfants_75 %>% 
  bind_rows(p_enfants_75) %>% 
  bind_rows(ap_enfants_75) %>% 
  tbl_df()

Nous calculons, pour chaque génération de descendants (enfants, petits-enfants et arrière-petits-enfants), le nombre d’aïeux nés dans chacun des départements français. Puis, nous calculons la part relative que représente chaque département de naissance des aïeux, en excluant Paris, pour voir si certaines régions sont davantage représentées que d’autres dans l’origine des descendants nés à Paris.

nb_obs_depts <- 
  generations_paris %>% 
  filter(str_detect(dept_N_aieul, "^F[[:digit:]]{2}")) %>% 
  filter(dept_N_aieul != "F75") %>% 
  group_by(generation, dept_N_aieul) %>% 
  summarise(nb_obs = n()) %>% 
  group_by(generation) %>% 
  mutate(pct = nb_obs / sum(nb_obs) * 100)

nb_obs_depts %>% arrange(-pct) %>% head()

# A tibble: 6 x 4
# Groups: generation [3]
  generation dept_N_aieul nb_obs   pct
       <dbl> <chr>         <int> <dbl>
1       1.00 F78              23  8.98
2       3.00 F67             267  8.54
3       2.00 F67              82  8.34
4       1.00 F77              17  6.64
5       2.00 F57              62  6.31
6       2.00 F77              59  6.00

Nous voulons proposer une représentation graphique de ces valeurs, sur une carte choroplèthe.

france_paris <- 
  france_1800 %>% 
  left_join(
    nb_obs_depts %>% 
      select(generation, pct, dept_N_aieul) %>% 
      spread(generation, pct) %>% 
      mutate(num_dep = str_sub(dept_N_aieul, 2)),
    by = "num_dep")

france_paris <- 
  france_paris %>% 
  gather(generation, pct, `1`,`2`, `3`) %>% 
  mutate(pct = ifelse(num_dep == "75", yes = NA, no = pct)) %>% 
  mutate(generation = factor(generation,
                             levels = c(0, 1, 2, 3),
                             labels = c('Aïeux', "Enfants", "Petits-enfants", "Arrière-petits-enfants")))

max_pct <- france_paris$pct %>% max(na.rm=T)
max_pct <- ceiling(max_pct)

p <- ggplot() +
  geom_polygon(data = france_paris, aes(x = long, y = lat, group = group, fill = pct), col = "grey", size = .4) +
  geom_polygon(data = france_paris %>% filter(departement == "F75"),
               aes(x = long, y = lat, group = group), fill = NA, size = .4, col = "grey50", linetype = "dashed") +
  coord_quickmap() +
  scale_fill_gradientn(
    "Pourcentage",
    colours = c("#FFFFFF", rev(heat.colors(10))), na.value = "white",
                       limits = c(0, max_pct),
                       guide = guide_colourbar(direction = "horizontal",
                                               barwidth = 30,
                                               title.position = "bottom",
                                               title.hjust = .5)) +
  facet_wrap(~generation) + 
  theme(text = element_text(size = 16), legend.position = "bottom")

p

On peut également regarder le pourcentage de descendants, dans chaque département, à naître à Paris.

pct_naiss_enfants <- 
  generations %>% 
  select(id_aieul, id_enfant, dept_N_aieul, dept_N_enfants) %>% 
  filter(!is.na(dept_N_aieul), !is.na(dept_N_enfants)) %>% 
  unique() %>% 
  tbl_df() %>% 
  mutate(naiss_paris = dept_N_enfants == "F75") %>% 
  group_by(dept_N_aieul) %>% 
  summarise(n = n(), n_paris = sum(naiss_paris), pct_paris = n_paris / n * 100)

pct_naiss_p_enfants <- 
  generations %>% 
  select(id_aieul, id_p_enfant, dept_N_aieul, dept_N_p_enfants) %>% 
  filter(!is.na(dept_N_aieul), !is.na(dept_N_p_enfants)) %>% 
  unique() %>% 
  tbl_df() %>% 
  mutate(naiss_paris = dept_N_p_enfants == "F75") %>% 
  group_by(dept_N_aieul) %>% 
  summarise(n = n(), n_paris = sum(naiss_paris), pct_paris = n_paris / n * 100)


pct_naiss_ap_enfants <- 
  generations %>% 
  select(id_aieul, id_ap_enfant, dept_N_aieul, dept_N_ap_enfants) %>% 
  filter(!is.na(dept_N_aieul), !is.na(dept_N_ap_enfants)) %>% 
  unique() %>% 
  tbl_df() %>% 
  mutate(naiss_paris = dept_N_ap_enfants == "F75") %>% 
  group_by(dept_N_aieul) %>% 
  summarise(n = n(), n_paris = sum(naiss_paris), pct_paris = n_paris / n * 100)

pct_naiss_paris <-
  pct_naiss_enfants %>% mutate(generation = 1) %>% 
  bind_rows(pct_naiss_p_enfants %>% mutate(generation = 2)) %>% 
  bind_rows(pct_naiss_ap_enfants %>% mutate(generation = 3))

stat_des_pct_naiss_paris <- 
  pct_naiss_paris %>% 
  group_by(generation) %>% 
  summarise(pct_paris_avg = mean(pct_paris),
            pct_paris_min = min(pct_paris),
            pct_paris_max = max(pct_paris)
            )

stat_des_pct_naiss_paris

# A tibble: 3 x 4
  generation pct_paris_avg pct_paris_min pct_paris_max
       <dbl>         <dbl>         <dbl>         <dbl>
1       1.00          1.69         0              73.3
2       2.00          3.30         0              53.6
3       3.00          5.11         0.256          42.2

La carte, qui montre pour chaque département, le pourcentage de descendants à naître à Paris, pour chaque département, s’obtient de la manière suivante :

france_paris <- 
  france_1800 %>% 
  left_join(
    pct_naiss_paris %>% 
      select(generation, pct_paris, dept_N_aieul) %>% 
      spread(generation, pct_paris) %>% 
      mutate(num_dep = str_sub(dept_N_aieul, 2)),
    by = "num_dep")

france_paris <- 
  france_paris %>% 
  gather(generation, pct, `1`,`2`, `3`) %>% 
  mutate(pct = ifelse(num_dep == "75", yes = NA, no = pct)) %>% 
  mutate(generation = factor(generation,
                             levels = c(0, 1, 2, 3),
                             labels = c('A\\"\\\'eux', "Enfants", "Petits-enfants", "Arri\\`ere-petits-enfants")))

max_pct <- france_paris$pct %>% max(na.rm=T)
max_pct <- ceiling(max_pct)

p <- ggplot() +
    geom_polygon(data = france_paris, aes(x = long, y = lat, group = group, fill = pct), col = "grey", size = .4) +
    geom_polygon(data = france_paris %>% filter(departement == "F75"),
                 aes(x = long, y = lat, group = group), fill = NA, size = .4, col = "grey50", linetype = "dashed") +
    coord_quickmap() +
    scale_fill_gradientn(
      "Pourcentage",
      colours = c("#FFFFFF", rev(heat.colors(10))), na.value = "white",
      limits = c(0, max_pct),
      guide = guide_colourbar(direction = "horizontal",
                              barwidth = 30,
                              title.position = "bottom",
                              title.hjust = .5)) +
    facet_wrap(~generation) + 
    theme(text = element_text(size = 16), legend.position = "bottom")

p

4.3.1 Proportion de descendants à Paris : un effet distance avec la capitale ?

La carte montrant les proportions de descendants à naître à Paris semble refléter un effet de distance : plus le département d’origine des aïeux est loin de Paris, moins le pourcentage de descendants à y naître semble élevé. Nous creusons alors dans ce sens.

Il nous faut une mesure de la distance entre chaque département et Paris. Nous retenons celle qui sépare les barycentres. Nous devons de facto obtenir le barycentre de chaque région.

library(rgeos)
obtenir_spa_poly_sp <- function(un_dep){
  dep_tmp <- france_1800 %>%
    filter(departement == un_dep)
  
  obtenir_poly_groupe <- function(un_groupe){
    coords_group <- dep_tmp[which(dep_tmp$group %in% un_groupe), c("long", "lat")]
    return(Polygons(list(Polygon(coords_group)), un_groupe))
  }
  dep_polys <- lapply(unique(dep_tmp$group), obtenir_poly_groupe)
  dep_polys_sp <- SpatialPolygons(dep_polys, seq(1, length(dep_polys)))
  dep_polys_sp
}

france_1800_sp <- lapply(unique(france_1800$departement), obtenir_spa_poly_sp)
names(france_1800_sp) <- unique(france_1800$departement)
france_1800_barycentres <- lapply(france_1800_sp, function(x) gCentroid(x))
france_1800_barycentres <- do.call("rbind", lapply(france_1800_barycentres, function(x) x@coords))
france_1800_barycentres <- cbind(france_1800_barycentres, departement = unique(france_1800$departement))
france_1800_barycentres <- data.frame(france_1800_barycentres, stringsAsFactors = FALSE) %>% tbl_df()
rownames(france_1800_barycentres) <- NULL
france_1800_barycentres <- 
  france_1800_barycentres %>% 
  mutate(long = as.numeric(x), lat = as.numeric(y)) %>% 
  select(-x,-y)

En particulier, les coordonnées du barycentre du département de Paris sont les suivantes :

(coords_paris <- france_1800_barycentres %>% filter(departement == "F75"))

# A tibble: 1 x 3
  departement  long   lat
  <chr>       <dbl> <dbl>
1 F75          2.37  48.9

Nous encapsulons la fonction distm() du package geosphere dans une autre, pour pouvoir calculer la distance séparant deux points du globe.

calcul_distance <- function(long_1, lat_1, long_2, lat_2){
  geosphere::distm(c(long_1, lat_1), c(long_2, lat_2), fun = geosphere::distHaversine) %>% "["(1)
}

Nous appliquons la fonction fraîchement créée pour obtenir les distances (en km) séparant chaque département de Paris :

distances_paris <- 
  france_1800_barycentres %>% 
  mutate(long_paris = coords_paris$long, lat_paris = coords_paris$lat) %>% 
  rowwise() %>% 
  mutate(distance_paris = calcul_distance(long, lat, long_paris, lat_paris) / 1000) %>% 
  ungroup() %>% 
  select(departement, distance_paris)

head(distances_paris)

# A tibble: 6 x 2
  departement distance_paris
  <chr>                <dbl>
1 F01                    380
2 F02                    118
3 F03                    282
4 F04                    607
5 F06                    657
6 F07                    482

Puis, nous ajoutons ces distances à nos données représentant le pourcentage de descendants dans chaque département à être né à Paris.

pct_naiss_paris <- 
  pct_naiss_paris %>% 
  left_join(distances_paris, by = c("dept_N_aieul" = "departement")) %>% 
  filter(dept_N_aieul != "F75") %>% 
  mutate(generation = factor(generation,
                             levels = c(0, 1, 2, 3),
                             labels = c('Aïeux', "Enfants", "Petits-enfants", "Arrière-petits-enfants")))

Pour attribuer les couleurs à chaque génération en accord avec celles utilisées précédemment, nous faisons appel à une fonction permettant de les retrouver (ce sont les couleurs qui sont utilisée automatiquement par les fonctions du package ggplot2).

gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}

couleurs <- gg_color_hue(4)

p_prop_paris_distance <- ggplot(data = pct_naiss_paris,
                                aes(x = distance_paris, y = pct_paris,
                                    colour = generation, group = generation
                                    )) +
  geom_point() +
  geom_smooth(se = FALSE) +
  xlab("Distance à Paris (km)") + ylab("Pourcentage") +
  scale_colour_manual("", values = c('AÏeux' = couleurs[1], "Enfants" = couleurs[2],
                                     "Petits-enfants" = couleurs[3], "Arrière-petits-enfants" = couleurs[4])) +
  scale_linetype_manual("", values = c("Enfants" = "solid",
                                       "Petits-enfants" = "dotted",
                                       "Arrière-petits-enfants" = "dashed")) +
  scale_shape_manual("", values = c("Enfants" = 15,
                                       "Petits-enfants" = 16,
                                       "Arrière-petits-enfants" = 17))

p_prop_paris_distance

4.3.2 Déplacements transitoires entre Paris et sa couronne

Nous explorons les mouvements entre la commune de Paris et sa couronne. Nous récupérons les codes des communes dans un rayon de 20km de Paris, que nous considérons comme la couronne. Pour ce faire, nous réutilisons l’objet communes_proches_20km, dont l’obtention est détaillée en annexe.

codes_insee_proches_20km <- communes_proches_20km[["75056"]]$limitrophes_20
codes_retenir <- c(codes_insee_proches_20km, "75056")

Pour chaque génération, nous calculons le pourcentage de descendants d’aïeux parisiens à naître à Paris, dans sa couronne ou ailleurs.

# Aieux --> enfants
villes_aieul_enfant <- 
  generations %>% 
  filter(code_insee_N_aieul %in% codes_retenir) %>% 
  filter(!is.na(code_insee_N_aieul), !is.na(code_insee_N_enfants)) %>% 
  mutate(type_commune_aieul = ifelse(code_insee_N_aieul == "75056", yes = "Paris", no = "Autre")) %>% 
  mutate(type_commune_aieul = ifelse(code_insee_N_aieul %in% codes_insee_proches_20km, yes = "Couronne", no = type_commune_aieul)) %>% 
  mutate(type_commune_enfants = ifelse(code_insee_N_enfants == "75056", yes = "Paris", no = "Autre")) %>% 
  mutate(type_commune_enfants = ifelse(code_insee_N_enfants %in% codes_insee_proches_20km, yes = "Couronne", no = type_commune_enfants)) %>% 
  select(id_aieul, id_enfant, type_commune_aieul, type_commune_enfants) %>% 
  unique() %>% 
  mutate(type_commune_aieul = factor(type_commune_aieul, levels = c("Paris", "Couronne", "Autre")),
         type_commune_enfants = factor(type_commune_enfants, levels = c("Paris", "Couronne", "Autre")))

villes_aieul_enfant <- 
  xtabs(~type_commune_aieul + type_commune_enfants, villes_aieul_enfant)

# Aieux --> petits-enfants
villes_aieul_p_enfant <- 
  generations %>% 
  filter(code_insee_N_aieul %in% codes_retenir) %>% 
  mutate(type_commune_aieul = ifelse(code_insee_N_aieul == "75056", yes = "Paris", no = "Autre")) %>% 
  mutate(type_commune_aieul = ifelse(code_insee_N_aieul %in% codes_insee_proches_20km, yes = "Couronne", no = type_commune_aieul)) %>% 
  mutate(type_commune_p_enfants = ifelse(code_insee_N_p_enfants == "75056", yes = "Paris", no = "Autre")) %>% 
  mutate(type_commune_p_enfants = ifelse(code_insee_N_p_enfants %in% codes_insee_proches_20km, yes = "Couronne", no = type_commune_p_enfants)) %>% 
  select(id_aieul, id_p_enfant, type_commune_aieul, type_commune_p_enfants, code_insee_N_p_enfants) %>% 
  unique() %>% 
  mutate(type_commune_aieul = factor(type_commune_aieul, levels = c("Paris", "Couronne", "Autre")),
         type_commune_p_enfants = factor(type_commune_p_enfants, levels = c("Paris", "Couronne", "Autre")))

villes_aieul_p_enfant <- 
  xtabs(~type_commune_aieul + type_commune_p_enfants, villes_aieul_p_enfant)

# Aieux --> arrière-petits-enfants
villes_aieul_ap_enfant <- 
  generations %>% 
  filter(code_insee_N_aieul %in% codes_retenir) %>% 
  mutate(type_commune_aieul = ifelse(code_insee_N_aieul == "75056", yes = "Paris", no = "Autre")) %>% 
  mutate(type_commune_aieul = ifelse(code_insee_N_aieul %in% codes_insee_proches_20km, yes = "Couronne", no = type_commune_aieul)) %>% 
  mutate(type_commune_ap_enfants = ifelse(code_insee_N_ap_enfants == "75056", yes = "Paris", no = "Autre")) %>% 
  mutate(type_commune_ap_enfants = ifelse(code_insee_N_ap_enfants %in% codes_insee_proches_20km, yes = "Couronne", no = type_commune_ap_enfants)) %>% 
  select(id_aieul, id_ap_enfant, type_commune_aieul, type_commune_ap_enfants, code_insee_N_ap_enfants) %>% 
  unique() %>% 
  mutate(type_commune_aieul = factor(type_commune_aieul, levels = c("Paris", "Couronne", "Autre")),
         type_commune_ap_enfants = factor(type_commune_ap_enfants, levels = c("Paris", "Couronne", "Autre")))

villes_aieul_ap_enfant <- 
  xtabs(~type_commune_aieul + type_commune_ap_enfants, villes_aieul_ap_enfant)


table_villes_aieul_enfant <- ftable(villes_aieul_enfant)
table_villes_aieul_enfant <- table_villes_aieul_enfant/rowSums(table_villes_aieul_enfant) * 100

table_villes_aieul_p_enfant <- ftable(villes_aieul_p_enfant)
table_villes_aieul_p_enfant <- table_villes_aieul_p_enfant/rowSums(table_villes_aieul_p_enfant) * 100

table_villes_aieul_ap_enfant <- ftable(villes_aieul_ap_enfant)
table_villes_aieul_ap_enfant <- table_villes_aieul_ap_enfant/rowSums(table_villes_aieul_ap_enfant) * 100

5 Annexes

5.1 Carte de France en 1800

Pour obtenir une carte de France (approximative) en 1800, nous nous appuyons d’abord sur les tracés actuels, fournis par le package maps.

france <- map_data(map = "france") %>% rename(nom_map = region) %>% tbl_df()
head(france)

Nous récupérons également les tracés de la base de données GADM. Le package raster propose une fonction pour télécharger les données.

map_fr_sp <- raster::getData(name = "GADM", country = "FRA", level = 2)

Les données de la base GADM permettent d’obtenir des noms mieux orthographiés que celles du package maps. Nous les extrayons, et procédons à quelques manipulations afin d’effectuer l’appariement entre les deux sources.

noms_departements <- map_fr_sp@data %>% select(CCA_2, NAME_2) %>% unique() %>% tbl_df() %>% 
  rename(num_dep = CCA_2, nom_sp = NAME_2) %>% 
  mutate(nom_jointure = iconv(nom_sp, to='ASCII//TRANSLIT')) %>% 
  mutate(nom_jointure = str_replace_all(nom_jointure, "[^[:alnum:]]", "")) %>% 
  mutate(nom_jointure = str_to_lower(nom_jointure)) %>% 
  arrange(nom_jointure)

Nous pouvons dès lors obtenir une jolie base de données pour tracer les frontières des départements français, que nous sauvegardons.

france <- 
  france %>% 
  mutate(nom_jointure = str_to_lower(nom_map)) %>% 
  mutate(nom_jointure = str_replace_all(nom_jointure, "[^[:alnum:]]", "")) %>% 
  tbl_df() %>% 
  arrange(nom_jointure) %>% 
  left_join(noms_departements)

save(france, file="france.rda")

Reste alors à effectuer quelques fusions de départements pour mieux refléter l’état de la France au XIXe siècle. Nous devons regrouper quelques départements, principalement ceux ayant été créé récemment, notamment dans le bassin parisien.

Nous créons une fonction permettant de transformer un département représenté par les coordonnées de son polygone en objet gpc.poly. Cette étape est utile pour pouvoir utiliser les fonctions d’union de polygones par la suite.

#' group_to_gpc_poly
#' Returns the polygon from a group within the data frame containing
#' the coordinates of regions
#' @df: data frame with the coordinates of the polygons
#' @group: (string) group ID
group_to_gpc_poly <- function(df, group){
  group_coords <- df[df$group == group, c("long", "lat")]
  as(group_coords, "gpc.poly")
}

Nous créons une autre fonction permettant d’effectuer la fusion de départements en un seul. Cette fonction procède de la manière suivante : nous mettons sous forme de polygone les départements à fusionner, créons l’union des ces polygones, puis récupérons les coordonnées de l’union des polygones.

#' fusion_depts
#' A partir d'un tableau contenant les coordonnées des frontières des départements et d'une liste d'identifiants
#' de départements, fusionne les départements cités en un seul polygone.
#' @depts_a_grouper (string) : vecteur de caractères indiquant les identifiants des départements à fusionner
#' @nouveau_nom (string) : nouveau nom du département issu de la fusion
#' @nouveau_departement (string) : nouvel identifiant du département issu de la fusion
#' @df : tableau à double entrée contenant les coordonnées des frontières des départements
fusion_depts <- function(depts_a_grouper, nouveau_nom, nouveau_departement, df){
  df_tmp <- df %>% filter(departement %in% depts_a_grouper)
  df_tmp_polys <- lapply(unique(df_tmp$group), group_to_gpc_poly, df = df_tmp)
  df_tmp_polys_tot <- df_tmp_polys[[1]]
  if(length(df_tmp_polys)>1){
    for(i in 2:length(df_tmp_polys)){
      df_tmp_polys_tot <- gpclib::union(df_tmp_polys_tot,  df_tmp_polys[[i]])
    }
  }
  coords <- data.frame(long = df_tmp_polys_tot@pts[[1]]$x, lat = df_tmp_polys_tot@pts[[1]]$y)
  noms_col <- colnames(df_tmp)
  res <- 
    df_tmp %>% 
    select(-long, -lat) %>% 
    slice(1)
  res <- 
    cbind(coords, res) %>% tbl_df() %>% 
    s_select(noms_col) %>% 
    mutate(nom_map = nouveau_nom,
           departement = nouveau_departement)
  res
}# Fin de fusion_depts()

Nous effectuons les trois fusions suivantes :

seine_oise <- fusion_depts(c("F78", "F91", "F92", "F95"), nouveau_nom = "Seine-et-Oise", nouveau_departement = "F78", df = france)
seine_oise$num_dep <- "78"

seine_marne <- fusion_depts(c("F93", "F94", "F77"), nouveau_nom = "Seine-et-Marne", nouveau_departement = "F77", df = france)
seine_marne$num_dep <- "77"

corse <- fusion_depts(c("F2A", "F2B"), nouveau_nom = "Corse", nouveau_departement = "F20", df = france)
corse$num_dep <- "20"

Nous retirons ces départements du tableau de données initial, puis rajoutons les coordonnées de la fusion, et sauvegardons le résultat.

france_1800 <- 
  france %>% 
  filter(! departement %in% c("F78", "F91", "F92", "F95", "F93", "F94", "F77", "F2A", "F2B")) %>% 
  bind_rows(seine_oise) %>% 
  bind_rows(seine_marne) %>% 
  bind_rows(corse) %>% 
  mutate(num_dep)

save(france_1800, file="france_1800.rda")

l_fusion <- c("F78", "F91", "F92", "F95", "F93", "F94", "F77", "F2A", "F2B")
france %>% 
  mutate(dept_a_fusionner = departement %in% l_fusion) %>% 
  ggplot(data = .,
         aes(x=long,y=lat,group=group)) +
  geom_polygon(col = "white", aes(fill = dept_a_fusionner)) +
  coord_quickmap()

l_fusionnes <- c("F78", "F77", "F20")
france_1800 %>% 
  mutate(depts_fusionnes = departement %in% l_fusionnes) %>% 
  ggplot(data = ., aes(x=long,y=lat,group=group)) +
  geom_polygon(col = "white", aes(fill = depts_fusionnes)) +
  coord_quickmap()

5.2 Carte des communes

Cette annexe propose d’indiquer comment récupérer les limites géographiques des communes à partir des données d’Open Street Map (OSM), puis des les réutiliser pour obtenir une extension des frontières des communes dans un rayon donné. L’obtention des extensions de frontières permettent, entre autre, d’identifier quelles communes se trouvent dans un voisinage de \(x\) kilomètres d’une autre commune. En choisissant une faible valeur de \(x\)x, cela permet d’identifier les communes limitrophes.

5.2.1 Récupération des données

Nous récupérons les données d’OSM via la plateforme ouverte des données publiques françaises data.gouv.fr. En particulier, nous téléchargeons le shapefile des communes le plus récent (Janvier 2018 au moment du téléchargement). La suite se passe avec R.

D’abord, nous chargeons quelques packages.

library(tidyverse)
library(lubridate)
library(stringr)
library(stringi)
library(data.table)
library(pbapply)
library(dplyrExtras)
library(stringdist)
library(sp)
library(rgeos)
library(maptools)
library(rgdal)

Nous chargeons les données :

communes <- readOGR(dsn="communes-20180101-shp/", layer="communes-20180101")
communes <- spTransform(communes, CRS( "+init=epsg:2154" ))

Nous allons extraire les informations de chaque commune de l’objet communes. Il faut dans ce cas faire attention à un détail technique qui concerne la taille de chaque objet contenant les informations d’une seule commune. Si l’on se contente d’extraire telle quelle une commune, dans l’état actuel, l’objet créé occupera une lourde place en mémoire, puisqu’il contiendra de nombreuses (très nombreuses) informations inutiles. En effet, le slot data de l’objet communes contient des informations : les codes communes INSEE, les noms et le liens Wikipedia. Or, ces informations sont stockées sous forme de facteurs : R considère donc chaque valeur comme des entiers, et se réfère à un dictionnaire indiquant les niveaux correspondants. Lorsque l’on extrait un facteur d’un vecteur de facteurs en R, on récupère un sous-élément de ce vecteur… et le dictionnaire au complet ! Aussi, lorsque ce dictionnaire est très volumineux, on perd en efficacité. Ici, comme chaque ligne du slot data contient un code INSEE unique, un nom unique et un lien wikipedia unique, il est sous-optimal de stocker ces informations sous forme de facteurs, de simples chaînes de caractères suffisent et se révèlent nettement plus efficace par la suite, lors des extractions de communes.

codes_insee <- unique(communes$insee) %>% as.character()
communes@data <- 
  communes@data %>% 
  mutate(insee = as.character(insee),
         nom = as.character(nom),
         wikipedia = as.character(wikipedia))

5.2.2 Agrandissement des limites des communes

Nous allons agrandir les polygones limitant chaque commune, selon une distance donnée. Pour ce faire, nous utilons la fonction gBuffer() du package rgeos. Nous choisissons d’étendre les frontières des communes de 20 km.

distance <- 20000 # en metres

Nous créons une fonction qui se chargera d’agrandir les frontières d’une commune, pour une distance donnée. Cette fonction retourne une liste contenant 4 objets : les coordonnées d’un rectangle délimitant les limites de la commune, celles d’un rectangle délimitant les limites étendues de la commune, puis les objets spatiaux contenant les coordonnées de la commune et de la commune étendue, respectivement.

#' communes_buffer
#' Obtenir la surface etendue de la commune
#' avec une distance donnee
#' @code_insee: (string) code INSEE de la commune
#' @distance: (num) distance pour etendre
communes_buffer <- function(code_insee, distance){
  tmp <- communes[communes$insee == code_insee,]
  tmp_buffer <- gBuffer(tmp, width = distance, byid = TRUE)
  bbox_commune <- bbox(tmp)
  bbox_commune_buffer <- bbox(tmp_buffer)
  tmp_buffer <- spTransform(tmp_buffer, CRS("+proj=longlat +datum=WGS84")) 
  tmp <- spTransform(tmp, CRS("+proj=longlat +datum=WGS84")) 
  list(bbox_commune = bbox_commune, bbox_commune_buffer = bbox_commune_buffer, tmp = tmp, tmp_buffer = tmp_buffer)
}# Fin de communes_buffer()

Voici un exemple du résultat pour une commune en particulier, Rennes, avec un facteur d’agrandissement de 1km.

res_rennes <- communes_buffer(code_insee = "35238", distance = 1000)

Les coordonnées du cadre délimitant la commune, et celles du cadre délimitant la commune étendue :

res_rennes$bbox_commune

        min       max
x  346353.8  356295.1
y 6785457.4 6793920.0

res_rennes$bbox_commune_buffer

        min     max
x  345356.3  357295
y 6784457.4 6794920

Les limites de la commune et l’agrandissement :

plot(res_rennes$tmp_buffer, border = "red")
plot(res_rennes$tmp, add=TRUE)

Pour permettre à l’ordinateur d’avoir à gérer de moins gros objets, nous séparons en 20 tronçons les 36 000 codes INSEE, appliquons la fonction communes_buffer() sur chaque code INSEE des tronçons, et sauvegardons le résultat de chaque tronçon.

chunk2 <- function(x,n) split(x, cut(seq_along(x), n, labels = FALSE)) 
a_parcourir <- chunk2(1:length(codes_insee), 20)

if(!(dir.exists(str_c("communes/")))) dir.create(str_c("communes/"), recursive = TRUE)
for(i in 1:length(a_parcourir)){
  communes_cercles_tmp <- 
    pblapply(codes_insee[a_parcourir[[i]]], communes_buffer, distance = distance)
  save(communes_cercles_tmp, file = str_c("communes/communes_cercles_tmp_", i, ".rda"))
  rm(communes_cercles_tmp)
}

Reste alors à charger les 20 résultats intermédiaires, pour obtenir les limites étendues de chaque communes :

communes_cercles <- 
  lapply(1:length(a_parcourir), function(i){
    load(str_c("communes/communes_cercles_tmp_", i, ".rda"))
    lapply(communes_cercles_tmp, function(x) x$tmp_buffer)
  })

communes_cercles <- unlist(communes_cercles)
names(communes_cercles) <- codes_insee

Puis celles des communes non étendues :

communes_sans_cercle <- 
  lapply(1:length(a_parcourir), function(i){
    load(str_c("communes/communes_cercles_tmp_", i, ".rda"))
    lapply(communes_cercles_tmp, function(x) x$tmp)
  })

communes_sans_cercle <- unlist(communes_sans_cercle)
names(communes_sans_cercle) <- codes_insee

Et enfin les rectangles délimitant les frontières des communes et des communes étendues :

communes_cercles_bbox <- 
  lapply(1:length(a_parcourir), function(i){
    load(str_c("communes/communes_cercles_tmp_", i, ".rda"))
    lapply(communes_cercles_tmp, function(x) x$bbox_commune_buffer)
  })

communes_cercles_bbox <- unlist(communes_cercles_bbox, recursive=FALSE)
names(communes_cercles_bbox) <- codes_insee

communes_bbox <- 
  lapply(1:length(a_parcourir), function(i){
    load(str_c("communes/communes_cercles_tmp_", i, ".rda"))
    lapply(communes_cercles_tmp, function(x) x$bbox_commune)
  })

communes_bbox <- unlist(communes_bbox, recursive=FALSE)
names(communes_bbox) <- codes_insee

Nous transformons ensuite en tableaux de données les objets spatiaux contenant les limites des communes.

options(warn=-1)
communes_cercles_df <- 
  pblapply(communes_cercles, function(x){
    suppressMessages(broom::tidy(x, region = "insee"))
  }) %>% 
  bind_rows()
options(warn=1)

Nous faisons de même pour les communes :

communes <- spTransform(communes, CRS("+proj=longlat +datum=WGS84"))
communes_df <- broom::tidy(communes, region = "insee")
communes_df <- tbl_df(communes_df)

5.2.3 Communes proches

À présent, nous pouvons utiliser les limites des communes étendues pour identifier, pour chacune des communes, les autres proches dans un rayon de 20km. Nous créons une fonction qui fonctionne en deux temps, pour une commune donnée. Dans un premier, nous utilisons les bounding box des communes pour réaliser un écrémage rapide des communes potentiellement proches de notre commune de référence. Cette étape vise à accélérer la seconde étape qui consiste à utiliser la fonction gIntersects() du package rgeos. Cette fonction, qui n’est pas des plus rapides à s’exécuter, indique si deux polygones s’intersectent. Elle nous permet donc d’identifier les communes en intersection avec la commune dont les limites ont été élargies de 20km.

#' trouver_intersection_commune
#' Pour la commune i de communes_cercles, retourne
#' l'indice Insee de cette commune et les indices Insee des
#' communes dans un rayon de 20km de cette commune
#' @i (int) : indice de la commune
trouver_intersection_commune <- function(i){
  comm_courante <- communes_cercles[[i]]
  comm_restantes <- communes_sans_cercle[-i]
  
  # On fait un premier ecremage à l'aide des box
  bbox_courante <- communes_cercles_bbox[[i]]
  bbox_restantes <- communes_bbox[-i]
  
  box_se_touchent <- function(x){
    # Est-ce que les box se touchent
    touche <- 
      bbox_courante["x", "min"] <= x["x", "max"] & bbox_courante["x", "max"] >= x["x", "min"] &
      bbox_courante["y", "min"] <= x["y", "max"] & bbox_courante["y", "max"] >= x["y", "min"]
    touche
  }# Fin de box_se_touchent()
  
  touchent <- sapply(bbox_restantes, box_se_touchent)
  
  inter <- sapply(comm_restantes[touchent], function(x){
    gIntersects(x, comm_courante)
  })
  
  insee_intersection <- names(comm_restantes)[which(touchent)[which(inter)]]
  
  list(insee = names(communes_cercles[i]), limitrophes_20 = insee_intersection)
}

Nous appliquons cette fonction à toutes les communes. Pour accélérer les choses, nous parallélisons l’exécution.

library(parallel)
ncl <- detectCores()-1
(cl <- makeCluster(ncl))

invisible(clusterEvalQ(cl, library(tidyverse, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(geosphere, warn.conflicts=FALSE, quietly=TRUE)))
invisible(clusterEvalQ(cl, library(rgeos, warn.conflicts=FALSE, quietly=TRUE)))
clusterExport(cl, c("communes_cercles", "communes_sans_cercle"), envir=environment())
clusterExport(cl, c("communes_cercles_bbox", "communes_bbox"), envir=environment())

communes_proches_20km <- pblapply(1:length(communes_cercles), trouver_intersection_commune, cl = cl)
names(communes_proches_20km) <- names(communes_cercles)

stopCluster(cl)

Reprenons l’exemple de la commune de Rennes pour observer le résultat.

ind_rennes <- which(names(communes_cercles) == "35238")
proche_rennes_20 <- trouver_intersection_commune(i = ind_rennes)
proche_rennes_20

$insee
[1] "35238"

$limitrophes_20
 [1] "35289" "35149" "35312" "35340" "35188" "35012" "35221" "35231"
 [9] "35212" "35030" "35322" "35218" "35333" "35066" "35266" "35206"
[17] "35352" "35055" "35334" "35001" "35039" "35203" "35143" "35144"
[25] "35146" "35346" "35265" "35233" "35193" "35079" "35274" "35197"
[33] "35296" "35251" "35007" "35003" "35315" "35067" "35085" "35286"
[41] "35165" "35338" "35264" "35283" "35087" "35252" "35229" "35052"
[49] "35141" "35057" "35175" "35127" "35155" "35099" "35176" "35173"
[57] "35069" "35060" "35169" "35309" "35282" "35208" "35035" "35196"
[65] "35198" "35242" "35164" "35022" "35139" "35135" "35245" "35180"
[73] "35134" "35089" "35258"
 [ reached getOption("max.print") -- omitted 102 entries ]

Regardons le résultat sur une carte :

map_rennes <- 
  ggplot(data = communes_df %>% 
           filter(id %in% unlist(proche_rennes_20)) %>% 
           mutate(limitrophe = ifelse(id %in% proche_rennes_20$limitrophes_20, yes = "limitrophe", no = "non"),
                  limitrophe = ifelse(id == proche_rennes_20$insee, yes = "focus", no = limitrophe))) +
  geom_polygon(data = map_data("france"), aes(x = long, y = lat, group = group), fill = NA, col = "white") +
  geom_polygon(aes(x = long, y= lat, group = group, fill = limitrophe)) +
  geom_polygon(data = communes_cercles[[ind_rennes]] %>% 
                 broom::tidy(region = "insee") %>% 
                 tbl_df(),
               aes(x = long, y=lat, group = group), fill = NA, col = "red", linetype = "dashed") +
  scale_fill_manual("", values = c("limitrophe" = "dodgerblue", "non" = "white", "focus" = "red"), guide = FALSE) +
  coord_quickmap(xlim = c(-5,0),
                 ylim = c(47.5,48.5))

map_rennes

Nous pouvons sauvegarder les résultats.

save(communes, file = "communes.rda")
save(communes_df, file = "communes_df.rda")
save(communes_cercles, file = "communes_cercles.rda")
save(communes_sans_cercle, file = "communes_sans_cercle.rda")
save(communes_cercles_bbox, file = "communes_cercles_bbox.rda")
save(communes_bbox, file = "communes_bbox.rda")
save(communes_proches_20km, file = "communes_proches_20km.rda")

5.3 Rivières

L’ajout des rivières sur certaines cartes peut être utile. Cette annexe propose une méthode pour le faire, en s’appuyant sur les données d’Open Street Map (OSM), dont une extraction est proposée sur la plateforme ouverte des données publiques françaises. Les données relatives aux cours d’eau sont disponibles à l’adresse suivante : https://www.data.gouv.fr/fr/datasets/carte-des-departements-2/. Comme elles sont proposées par département, nous allons scraper les liens vers les fichiers de chaque département, puis les télécharger à l’aide d’une fonction.

library(tidyverse)
library(lubridate)
library(stringr)
library(stringi)
library(data.table)
library(pbapply)
library(dplyrExtras)
library(stringdist)
library(sp)
library(rgeos)
library(maptools)
library(rgdal)
library(rvest)

if(!(dir.exists(paste0("osm_datagouv")))) dir.create(paste0("osm_datagouv"), recursive = TRUE)

page_datagrouv <- read_html("https://www.data.gouv.fr/fr/datasets/carte-des-departements-2/")

liens <- 
  page_datagrouv %>% 
  html_nodes(".resources-list") %>% 
  .[[1]] %>% 
  html_nodes("h4 a") %>% 
  html_attr("href")

telecharger_fichiers <- function(lien){
  nom_fichier <- str_replace(string = lien, pattern = "^.*?departement-", replacement = "")
  nom_dep <- str_replace(string = nom_fichier, "\\.zip", "")
  download.file(lien, destfile = str_c("osm_datagouv/", nom_fichier))
  unzip(str_c("osm_datagouv/", nom_fichier), exdir = str_c("osm_datagouv/", nom_dep))
  file.remove(str_c("osm_datagouv/", nom_fichier))
}

Nous téléchargeons les fichiers avec l’instruction suivante :

pblapply(liens[1:2], telecharger_fichiers)

Nous allons stocker les données, par département, dans un dossier spécifique :

if(!(dir.exists(paste0("rivieres_osm")))) dir.create(paste0("rivieres_osm"), recursive = TRUE)

Nous récupérons la liste des fichiers d’OSM :

dossiers_departements <- list.files("osm_datagouv", full.names = TRUE)

Puis, nous créons une fonction qui importe les fichiers shapefile récupérés, les simplifie légèrement et retourne le résultat sous forme d’un tableau de données.

riviere_rda <- function(un_dep){
  num_dep <- str_replace(un_dep,"osm_datagouv/", "") %>% 
    str_replace("-(.*?)$", "")
  chemin <- str_c(un_dep, "/departement-", num_dep)
  rivieres <- readOGR(dsn=chemin, layer="waterways")
  rivieres <- gSimplify(rivieres, tol = 0.0001)
  
  rivieres@data$osm_id <- as.character(rivieres@data$osm_id)
  rivieres@data$name <- as.character(rivieres@data$name)
  rivieres@data$type <- as.character(rivieres@data$type)
  rivieres@data$width <- as.character(rivieres@data$width)
  rivieres@data$id <- rownames(rivieres@data)
  rivieres@data$dep <- num_dep
  rivieres_df <- broom::tidy(rivieres, region = "id")
  rivieres_df <- rivieres_df %>% left_join(rivieres@data)
  save(rivieres_df, file = str_c("rivieres_osm/", num_dep, ".rda"))
}

Enfin, nous appliquons la fonction riviere_rda() pour traiter tous les départements.

pblapply(dossiers_departements, riviere_rda)

References

Blayo, Y. (1975). Mouvement naturel de la population française de 1740 a 1829. Population, 30, 15. https://doi.org/10.2307/1530644

Fleury, M., & Henry, L. (1958). Pour connaitre la population de la France depuis Louis XIV. Plan de travaux par sondage. Population, 13(4), 663. https://doi.org/10.2307/1525088

Walle, E. van de. (1973). La mortalité des départements français ruraux au XIXe siècle. Annales de démographie historique, 1973(1), 581‑589.