Vers une chaine de traitement reproductible

Exemple d’application



Avant propos et objectifs de séance

Une chaine de traitement reproductible ?

Depuis la collecte des données jusqu’à la valorisation des résultats, produire une analyse reproductible totalement en R

Une chaine de traitement reproductible ?

Depuis la collecte des données jusqu’à la valorisation des résultats, produire une analyse reproductible totalement en R

Utiliser plusieurs packages et leurs fonctions associées


Type d’opération Packages Application
Importer des données readxl Fichiers Excel
sf Données spatiales
Manipuler des données dplyr Manipulation de dataframes (pipes)
stringr Chaînes de caractères
questionr Recodage / transformation de variables
sf Données spatiales
potential Potentiel de Stewart
Visualiser ggplot2 Représentations graphiques de qualité
ggthemes Thèmes pour graphiques ggplot
cartography Cartographie
Valoriser / Documenter rmarkdown Literate programming (code + texte)
spelling Correction fautes d’orthographe

Transposer son analyse

Documenter ses méthodes

Créer un rapport qui associe code, commentaires et sorties graphiques

Se lancer !

Quelques conseils pour débuter en R et mise à disposition des données de l’exercice pour reproduire chez soi les analyses.

Mise en oeuvre d’une chaine de traitement

Données et thématique

Analyse de l’évolution démographique 1968-2017 des communes la métropole du Grand Paris (MGP).

3 fichiers sont utilisés, distribués par l’INSEE et l’IGN.

Population INSEE

Population communale par âge et sexe 1968-2016 (fichier Excel)

Table d’appartenance communale

EPCI de rattachement des communes françaises + recherches pour identifier des EPCI intermédiaires (csv)

Découpages géographiques IGN

Découpages géographiques des communes françaises (par version, shapefiles)

Charger les librairies et les données

Chargement des librairies

## Librairies utiles
library(readxl) # Lire des fichiers Excel

## Manipulation de données
library(dplyr) # Manipulation de données - "pipes"
library(stringr) # Manipulation de chaines de caractères
library(questionr) # Opérations de recodage / transformations de variables

## Visualisation
library(ggplot2) # Créer des visualisations utilisant une grammaire graphique
library(ggthemes) # Thèmes graphiques associés à ggplot

## Spatial 
library(sf) # manipulation de données spatiales
library(cartography) # cartographie
library(potential) # potentiels de Stewart 

Charger les librairies et les données

Les données tabulaires sont importées au format Excel…

# Importer à partir d'un fichier Excel (INSEE) 
library(readxl)

data2017 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", 
                                  skip = 7, 
                                  sheet = "2017")) 
data1990 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", 
                                  skip = 7, 
                                  sheet = "1990")) 
data1968 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", 
                                  skip = 7, 
                                  sheet = "1968")) 

… et au format .csv

# Importer à partir d'un fichier .csv 
metro <- read.csv(file = "input/metropoles.csv", 
                  sep = "\t", 
                  colClasses = c(rep("character", 4)),
                  encoding = "latin1")

Charger les librairies et les données

Le fond de carte (.shp) est finalement importé

# Importer des géométries communales (IGN - Geofla 2016)
library(sf)
com <- st_read(dsn = "input/COMMUNE.shp", stringsAsFactors = F, quiet = TRUE)

par(mar = c(0,0,0,0), bg = "#282828")
plot(com$geometry, col = 2:7, border = NA)

Préparation des données

Regrouper les données

# Jointure : On regroupe les objets data2015, data1990 et data1968 en un seul 
# objet "data"
data <- merge(data2017, data1990, by.x = "COM", by.y = "COM", all.x = T)
data <- merge(data, data1968, by.x = "COM", by.y = "COM", all.x = T) 

# Jointure 2 : EPCI d'appartenance
data <- merge(data, metro, by.x = "COM", by.y = "INSEE_COM", all.y = T)

# Jointure 3 : avec les découpages territoriaux
com <- merge(com, data, by.x = "INSEE_COM", by.y = "COM", 
                  all.y = T)

colnames(com)
##  [1] "INSEE_COM"    "ID_GEOFLA"    "CODE_COM"     "NOM_COM"      "STATUT"      
##  [6] "X_CHF_LIEU"   "Y_CHF_LIEU"   "X_CENTROID"   "Y_CENTROID"   "Z_MOYEN"     
## [11] "SUPERFICIE"   "POPULATION"   "CODE_ARR"     "CODE_DEPT"    "NOM_DEPT"    
## [16] "CODE_REG"     "NOM_REG"      "NCC.x"        "PMUN17"       "NCC.y"       
## [21] "PSDC90"       "NCC"          "PSDC68"       "METRO_LABEL"  "EPCI_SUB"    
## [26] "LIB_EPCI_SUB" "geometry"

Préparation des données

Nettoyer l’espace de travail puis renommer les champs

# Supprimer les objets devenus obsolèthes
rm(data1968, data1990, data2017, metro, data)

# Conserver uniquement certaines variables et les renommer
fields <- c("INSEE_COM", "NOM_COM", "PSDC68", "PSDC90", "PMUN17",
            "METRO_LABEL", "EPCI_SUB", "LIB_EPCI_SUB")
com <- com[,fields]

colnames(com)[1:8] <- c("COM","Nom","POP_1968","POP_1990","POP_2017",
                             "METRO_LABEL", "EPCI_LIB", "EPCI_LABEL")
head(com)
## Simple feature collection with 6 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 1028158 ymin: 6298724 xmax: 1055399 ymax: 6345834
## projected CRS:  RGF93_Lambert_93
##     COM                Nom POP_1968 POP_1990 POP_2017 METRO_LABEL EPCI_LIB
## 1 06006          ASPREMONT      698     1496     2187        Nice       T1
## 2 06009            BAIROLS       37       73      105        Nice       T2
## 3 06011   BEAULIEU-SUR-MER     4050     4013     3715        Nice       T1
## 4 06013          BELVEDERE      501      519      677        Nice       T3
## 5 06020 LA BOLLENE-VESUBIE      243      308      575        Nice       T3
## 6 06021             BONSON      232      456      737        Nice        R
##                    EPCI_LABEL                       geometry
## 1         CU Nice Côte d'Azur MULTIPOLYGON (((1041591 630...
## 2              CC de la Tinée MULTIPOLYGON (((1031138 633...
## 3         CU Nice Côte d'Azur MULTIPOLYGON (((1049887 630...
## 4       CC Vésubie-Mercantour MULTIPOLYGON (((1053577 633...
## 5       CC Vésubie-Mercantour MULTIPOLYGON (((1053577 633...
## 6 Rattachement à la Métropole MULTIPOLYGON (((1037263 631...

Préparation des données

Créer de nouveaux indicateurs

# Variable quantitative : taux de croissance moyens annuels 1968-1990, 1990-2017 et 1968-2017
com$acc_6890 <- 100 * (exp(log(com$POP_1990 / com$POP_1968) / (1990 - 1968))-1)
com$acc_9017 <- 100 * (exp(log(com$POP_2017 / com$POP_1990) / (2017 - 1990))-1)
com$acc_6817 <- 100 * (exp(log(com$POP_2017 / com$POP_1968) / (2017 - 1968))-1)

# Discrétiser une variable quantitative pour créer des classes quali
com[!is.finite(com$acc_6817),] <- NA
brks <- quantile(com$acc_6817, probs = seq(0,1,0.25), na.rm = TRUE)

# Typologie basée sur les quartiles
com$typo6817 <- cut(com$acc_6817, include.lowest = TRUE, 
                    breaks = c(brks[1], brks[2], brks[3], brks[4], brks[5]), 
                    labels = c(paste(round(brks[1],2), round(brks[2],2), sep = " - "),
                               paste(round(brks[2],2), round(brks[3],2), sep = " - "),
                               paste(round(brks[3],2), round(brks[4],2), sep = " - "),
                               paste(round(brks[4],2), round(brks[5],2), sep = " - ")))

# Extraire une chaine de caractère : département d'appartenance de chaque commune
library(stringr)
com$DEP <- str_sub(com$COM,  1, 2)

head(st_set_geometry(com[,c(1:2,10:14)],NULL))
##     COM                Nom    acc_6890   acc_9017   acc_6817     typo6817 DEP
## 1 06006          ASPREMONT  3.52587667  1.4163659  2.3581230  2.31 - 7.03  06
## 2 06009            BAIROLS  3.13702432  1.3554031  2.1514756  1.38 - 2.31  06
## 3 06011   BEAULIEU-SUR-MER -0.04170853 -0.2853713 -0.1760453 -1.39 - 0.53  06
## 4 06013          BELVEDERE  0.16057324  0.9891841  0.6163101  0.53 - 1.38  06
## 5 06020 LA BOLLENE-VESUBIE  1.08327236  2.3390486  1.7733124  1.38 - 2.31  06
## 6 06021             BONSON  3.11927647  1.7940328  2.3869202  2.31 - 7.03  06

Préparation des données

Sélectionner des données

# Quelles métropoles disponibles ? 
levels(as.factor(com$METRO_LABEL))
##  [1] "Bordeaux"         "Brest"            "Clermont-Ferrand" "Dijon"           
##  [5] "Grand Paris"      "Grenoble"         "Lille"            "Lyon"            
##  [9] "Marseille"        "Metz"             "Montpellier"      "Nancy"           
## [13] "Nantes"           "Nice"             "Orléans"          "Rennes"          
## [17] "Rouen"            "Saint-Etienne"    "Strasbourg"       "Toulon"          
## [21] "Toulouse"         "Tours"
# Définir un centre (pour les calculs de distance)
moncentre <- "75101"

# Mon espace d'étude
study <- "Grand Paris"

# Ne conserver que les unités géographiques désirées
com <- com[com$METRO_LABEL == study,]

# Supprimer les valeurs manquantes
com <- com[!is.na(com$METRO_LABEL),]

Préparation des données

Aggréger des données

library(dplyr)
# Avec dplyr : communes > Départements 
dep <- com %>% group_by(DEP) %>% summarize(names = first(DEP), .groups = "drop")

# En Rbase : communes > EPCI (+ somme de la population)
epci <- aggregate(x = com[,c("POP_1968", "POP_1990", "POP_2017")], 
                  by = list(EPCI_LABEL = com$EPCI_LABEL),
                  FUN = sum)

head(epci, 5)
## Simple feature collection with 5 features and 4 fields
## Attribute-geometry relationship: 0 constant, 3 aggregate, 1 identity
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 637610.4 ymin: 6853777 xmax: 670929.8 ymax: 6875005
## projected CRS:  RGF93_Lambert_93
##                                   EPCI_LABEL POP_1968 POP_1990 POP_2017
## 1 Association des Communes de l'Est Parisien   445010   462863   508966
## 2                             Boucle Nord 92   401150   384546   440565
## 3                               Est Ensemble   357223   361355   422744
## 4                            Grand-Paris Est   274477   339559   397765
## 5                      Grand Paris Sud Ouest   284832   272160   318893
##                         geometry
## 1 POLYGON ((661417.8 6854301,...
## 2 POLYGON ((650140.2 6866932,...
## 3 POLYGON ((660743.7 6864361,...
## 4 POLYGON ((670317.8 6857208,...
## 5 POLYGON ((646995.3 6857373,...

Préparation des données

Calcul de distances entre objets géographiques

# Extraction du centre prédéfini
center <- com[com$COM == moncentre,] 
# Distance euclidienne du centre au reste des centroides de l'espace d'étude
com$dist <- as.vector(c(st_distance(st_centroid(x = center), 
                                    st_centroid(x = com))))

head(st_set_geometry(com, NULL), 5)
##       COM                      Nom POP_1968 POP_1990 POP_2017 METRO_LABEL
## 777 75101 PARIS-1ER-ARRONDISSEMENT    32332    18360    16266 Grand Paris
## 778 75102  PARIS-2E-ARRONDISSEMENT    35357    20738    20900 Grand Paris
## 779 75103  PARIS-3E-ARRONDISSEMENT    56252    35102    34115 Grand Paris
## 780 75104  PARIS-4E-ARRONDISSEMENT    54029    32226    28088 Grand Paris
## 781 75105  PARIS-5E-ARRONDISSEMENT    83721    61222    58850 Grand Paris
##     EPCI_LIB EPCI_LABEL  acc_6890    acc_9017   acc_6817     typo6817 DEP
## 777       T1      Paris -2.539397 -0.44750439 -1.3922186 -1.39 - 0.53  75
## 778       T1      Paris -2.395960  0.02882413 -1.0672179 -1.39 - 0.53  75
## 779       T1      Paris -2.120750 -0.10557733 -1.0154306 -1.39 - 0.53  75
## 780       T1      Paris -2.321481 -0.50771100 -1.3261858 -1.39 - 0.53  75
## 781       T1      Paris -1.412579 -0.14624395 -0.7168023 -1.39 - 0.53  75
##          dist
## 777    0.0000
## 778  787.5588
## 779 1734.1135
## 780 1811.2153
## 781 2267.9806

Analyses et représentations graphiques

Univarié / résumé statistique

# Pour les graphiques, retirer les géométries des communes et epci
dfcom <- st_set_geometry(com, NULL)
dfepci <- st_set_geometry(epci, NULL)

# Résumé statistique sur certains champs du tableau de données
summary(dfcom[,c("POP_1968", "POP_1990", "POP_2017", "acc_6890", "acc_9017",
                 "acc_6817", "typo6817","dist")])
##     POP_1968         POP_1990         POP_2017         acc_6890       
##  Min.   :   437   Min.   :  1594   Min.   :  1810   Min.   :-2.53940  
##  1st Qu.: 16215   1st Qu.: 17882   1st Qu.: 20228   1st Qu.:-0.45582  
##  Median : 28459   Median : 29472   Median : 33416   Median : 0.01718  
##  Mean   : 44211   Mean   : 42373   Mean   : 47053   Mean   : 0.40161  
##  3rd Qu.: 53809   3rd Qu.: 48460   3rd Qu.: 56501   3rd Qu.: 0.67041  
##  Max.   :244080   Max.   :223940   Max.   :233392   Max.   :10.69107  
##     acc_9017          acc_6817                typo6817        dist      
##  Min.   :-0.7497   Min.   :-1.39222   -1.39 - 0.53:104   Min.   :    0  
##  1st Qu.: 0.1739   1st Qu.: 0.06745   0.53 - 1.38 : 29   1st Qu.: 6904  
##  Median : 0.4716   Median : 0.31633   1.38 - 2.31 : 10   Median :10397  
##  Mean   : 0.5123   Mean   : 0.45908   2.31 - 7.03 :  7   Mean   :10824  
##  3rd Qu.: 0.8253   3rd Qu.: 0.63184                      3rd Qu.:14621  
##  Max.   : 2.4239   Max.   : 4.77500                      Max.   :24753

Analyses et représentations graphiques

Univarié / graphique en bâton n°1

Population des communes de l’espace d’étude (20 classes)

# barplot de la population communale aux 3 recensements avec une boucle
par(mar = c(2,2,2,2), mfrow = c(1,3)) # marge et disposition

for (i in 3:5){ # On prend les colonnes 3 à 5 du dataframe
  popYear <- colnames(dfcom) # récupérer le nom des colonnes
  hist(dfcom[ ,i], # un histogramme
       breaks = 20, # la série statistique est coupée en 20 classes
       main = popYear[i], # titre
       col = "#4286f4", # couleur des barres
       border = "white", # couleur des bordures
       xlab = "Population", # titre de l'axe des X
       ylab = "Fréquence") # titre de l'axe des Y
}

Analyses et représentations graphiques

Univarié / graphique en bâton n°2 (code)

Population cumulée par EPCI aux trois recensements

# Ordonner dfepci de façon croissante sur la pop
dfepci <- dfepci[order(dfepci$POP_2017),] 

par(mar = c(4,4,4,4), mfrow = c(1,1)) # Graphique
barplot(t(as.matrix(dfepci[, c("POP_1968","POP_1990", "POP_2017")])), 
        beside = TRUE, # Barres côte à côte par recensment
        names.arg = str_wrap(dfepci$EPCI_LABEL, width = 15), # gestion labels
        legend.text = c("1968", "1990", "2017"), # Légende
        col = c("#6baed6", "#3182bd","#0d4880"), # Couleur des barres
        las = 2, # Légende des absisses à la verticale
        cex.names = 0.7, # Taille des fontes (X)
        cex.axis = 0.7, # Taille des fontes (Y)
        border = NA, # Pas de bordure autour des barres
        ylim = c(10000,max(dfepci$POP_2017)), # Bornes (Y)
        ylab = "Population résidente totale", # titre (Y)
        log = "y") # Echelle logarithmique

Analyses et représentations graphiques

Univarié / graphique en bâton n°2 (sortie graphique)

Population cumulée par EPCI aux trois recensements

Analyses et représentations graphiques

Quantitatif ~ Quantitatif (statistique)

Effet de la distance au centre sur la croissance démographique ?

# Test corrélation
cor.test(dfcom$dist, dfcom$acc_6890)
## 
##  Pearson's product-moment correlation
## 
## data:  dfcom$dist and dfcom$acc_6890
## t = 12.157, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6164249 0.7788846
## sample estimates:
##       cor 
## 0.7068566
cor.test(dfcom$dist, dfcom$acc_9017)
## 
##  Pearson's product-moment correlation
## 
## data:  dfcom$dist and dfcom$acc_9017
## t = 4.6968, df = 148, p-value = 0.000005982
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2121449 0.4920233
## sample estimates:
##       cor 
## 0.3601614

Analyses et représentations graphiques

Quantitatif ~ Quantitatif (graphique R:base)

Effet de la distance au centre sur la croissance démographique ?

# En quelques instructions on peut obtenir la visualisation des régressions linéaires ! 
par(mar = c(4,4,2,4), mfrow = c(1,2), bg = "white")

plot(dfcom$dist, dfcom$acc_6890)
abline(lm(acc_6890 ~ dist, data = dfcom), col = "red")
plot(dfcom$dist, dfcom$acc_9017)
abline(lm(acc_9017 ~ dist, data = dfcom), col = "red") 

Analyses et représentations graphiques

Quantitatif ~ Quantitatif (ggplot2 - code)

Effet de la distance au centre sur la croissance démographique avec taille et effet d’appartenance?

library(ggplot2)
library(ggthemes)

# Affecter des couleurs selon les départements d'appartenance des EPCI
cols <- c("#EA3531", "#A6CC99", "#92C8E0",  "#7BB6D3" ,"#71A966", "#8CBB80", 
          "#c49edb", "#EF6860", "#458DB3", "#4E9345", "#7BB6D3",  "#F38F84" )

par(mfrow = c(1,1)) # Un graphique par plot 

# Graphique X/Y, avec couleurs et cercles proportionnels
ggplot(dfcom, aes(x = dist, y = acc_6817, color = EPCI_LABEL)) + # Variables
  geom_point(aes(size = POP_2017)) + # Points proportionnels
  theme_wsj() + # Theme "Wall Street Journal" 
  scale_color_manual (values = cols) + # Affectation couleurs
  xlab("Distance au 1er arrondissement de Paris (m)") + 
  ylab("Taux d'accroissement moy. annuel (%)") +
  ggtitle("Croissance démographique 1968 - 2017,\ndistance au centre et appartenance territoriale") +
  labs(subtitle = "Sources: INSEE-IGN, 2020 / Réalisation : R.Ysebaert, 2020") +
  theme(axis.text = element_text(size = 6), 
        axis.title = element_text(size = 10, face = "bold"),
        legend.title = element_text(size = 8, face = "bold"),
        legend.text = element_text(size = 6),
        plot.subtitle = element_text(size = 8),
        plot.title = element_text(size = 15, face = "bold")) 

Analyses et représentations graphiques

Quantitatif ~ Quantitatif (ggplot2 - sortie graphique)

Effet de la distance au centre sur la croissance démographique avec taille et effet d’appartenance?

Analyses et représentations graphiques

Qualitatif ~ Qualitatif (statistique)

Lien entre la catégorie d’évolution démographique et l’appartenance des communes ?

# Test du Chi-2
tab <- table(dfcom$typo6817, dfcom$EPCI_LABEL) # Croisement des modalités
tab.chi2 <- chisq.test(tab) # Chi-2
tab.chi2
## 
##  Pearson's Chi-squared test
## 
## data:  tab
## X-squared = 117.14, df = 33, p-value = 0.0000000000238
library(questionr)
ctab <- cprop(tab, total = FALSE) # Calcul proportion par EPCI des catégories démo

Analyses et représentations graphiques

Qualitatif ~ Qualitatif (graphique - code)

Lien entre la catégorie d’évolution démographique et l’appartenance des communes ?

# Couleurs désirées en fonction des niveaux de la typo
coldf <- data.frame(cols = c("#fc8d59", "#a6d96a","#1a9641", "#085721"),
                    typo6817 = levels(dfcom$typo6817))

# Préparation des données et gestion des couleurs
ctab <- ctab[, order(-ctab[1,], -ctab[2,])] # Ordonner selon catégories croissantes
cols <- coldf[coldf$typo6817 %in% dfcom[,"typo6817"],] # Ne sélectionner que les couleurs correspondantes aux catégories présentes
cols <- cols$cols 

# Diagramme en bâton
par(mar = c(8,4,2,0), mfrow = c(1,1), xpd = TRUE) # Légende en dehors du graphique

barplot(ctab, main = "EPCI de la Métropole du Grand Paris", las = 2, 
        col = cols, 
        border = NA,
        names.arg = str_wrap(colnames(ctab), width = 15),
        ylab = "% des communes", # titre (Y)
        cex.lab = 0.6,
        cex.axis = 0.6, 
        cex.names = 0.6, 
        cex.main = 0.6)

legend(0, -40, title = "Croissance démographique 1968-2017 (% moyen annuel)\nQuartiles de la distribution des communes françaises",
       legend = levels(dfcom$typo6817),
       horiz = TRUE, # Positionnement horizontal
       fill = cols, # Couleur des caissons de légende 
       cex = 0.7, 
       bty = "n")

Analyses et représentations graphiques

Qualitatif ~ Qualitatif (sortie graphique)

Lien entre la catégorie d’évolution démographique et l’appartenance des communes ?

Analyses et représentations graphiques

Quantitatif ~ Qualitatif (statistique et code sortie graphique)

L’EPCI d’appartenance a-t-il un effet sur la croissance démographique ?

# Anova et test de Fisher
anova(lm(acc_6817 ~ EPCI_LABEL, data = dfcom))

# Boxplot taux de croissance 1968-2017 ~ EPCI 
# Ordonner sur la population
dfcom$EPCI_LABEL <- with(dfcom, reorder(EPCI_LABEL, acc_6817, mean,
                                        na.rm = TRUE))

# Couleurs
col <- carto.pal(pal1 = "red.pal", n1 = (nlevels(dfcom$EPCI_LABEL)/2), 
                 pal2 = "green.pal",n2 = (nlevels(dfcom$EPCI_LABEL)/2))

# Boxplot
par(mar = c(4,7,1,1)) # Gestion des marges
boxplot(dfcom$acc_6817 ~ dfcom$EPCI_LABEL, # Les deux variables à mettre en relation
        col = col, # Couleurs
        xlab = "Taux de croissance moyen annuel 1968-2017 (%)", ylab = NA, # Légendes
        names = str_wrap(levels(dfcom$EPCI_LABEL), width = 30), # Labels des EPCI
        horizontal = TRUE, #  Boxplot à l'horizontale
        varwidth = TRUE, # Box proportionnelles au nombre d'observation
        cex.axis = .5, cex.lab = .8, # Paramètres de mise en page
        outline = TRUE,
        las = 1) 

# Lignes de repère
abline (v = seq(-5,5,.5), col = "darkgrey", lwd = .5, lty = 3)

# Représenter les valeurs moyennes
xi <- tapply(dfcom$acc_6817, dfcom$EPCI_LABEL, mean, na.rm = T)
points(x = xi, y = seq(1,nlevels(dfcom$EPCI_LABEL),1), col="white",pch = 19)

Analyses et représentations graphiques

Quantitatif ~ Qualitatif (sortie graphique)

L’EPCI d’appartenance a-t-il un effet sur la croissance démographique ?

## Analysis of Variance Table
## 
## Response: acc_6817
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## EPCI_LABEL  11 72.141  6.5583  19.593 < 2.2e-16 ***
## Residuals  138 46.193  0.3347                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cartographie thématique

Carte choroplèthe (code)

#  2 cartes en vis à vis
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,2), bg = "lightgrey") # Ajuster les marges

# Définition des palettes
cols <- carto.pal(pal1 = "orange.pal", pal2 = "purple.pal",  n1 = 3, n2 = 3)

# 1968 - 1990
choroLayer(x = com, var = "acc_6890", # Objet et variable
           border = NA, # Couleur de bordure des territoires
           col = cols,  # Couleur de fond
           method = "quantile", # Méthode de discrétisation
           nclass = 6,  # Nombre de classes
           legend.pos = "bottom",  # Postion de légende
           legend.horiz = TRUE, #  Légende horizontale
           legend.values.rnd = 1, # Arrondi des valeurs en légende
           legend.title.txt = "Taux de croissance moyen annuel (%),\nQuantiles")

plot(epci$geometry, border = "white", lwd = .5, add = TRUE) # Ajouter les limites d'EPCI

# Mise en page
layoutLayer(title = "Évolution démographique 1968-1990",
            scale = 5,
            col = "lightgrey",
            coltitle = "black")

# 1990 - 2017
choroLayer(x = com, var = "acc_9017", 
           border = NA,
           col = cols, 
           method = "quantile",
           nclass = 6, 
           legend.pos = "bottom", 
           legend.horiz = TRUE,
           legend.values.rnd = 1,
           legend.title.txt = "Taux de croissance moyen annuel (%),\nQuantiles")

plot(epci$geometry, border = "white", lwd = .5, add = TRUE)

layoutLayer(title = "Évolution démographique 1990-2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE,
            horiz = FALSE,
            col = "lightgrey",
            coltitle = "black")

Cartographie thématique

Carte choroplèthe (sortie graphique)

Cartographie thématique

Combiner 2 variables : masse et évolution démographique (code)

# Cercles proportionnels et carte choroplèthe
# 1968 - 1990
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,2), bg = "lightgrey") # Ajuster les marges
plot(st_geometry(com), col = "lightblue", border = NA) # Fond des communes
plot(st_geometry(epci), col = NA, border = "white", lwd = 1, add = T) # Ajouter EPCI

propSymbolsChoroLayer(x = com, var = "POP_1990", var2 = "acc_6890",  # 2 variables
                      col = cols, # couleurs
                      inches = 0.15,  # Taille du cercle maximal...
                      fixmax = max(com$POP_2017), #... Basé sur le max de pop en 2017
                      method = "quantile", # Discrétisation
                      nclass = 6, # Nombre de classes
                      border = "grey50", # Couleur bordures de cercles
                      legend.var.title.txt = "Population totale", # Title légende var1
                      legend.var.values.rnd = -2, # Arrondi valeurs var1
                      legend.var.pos = "bottomleft",  # Positionnement légende var1
                      legend.var2.title.txt =  "Taux de croissance moyen annuel (%)", # Titre légende var2
                      legend.var2.values.rnd = 1, # Aronndis valeurs var2
                      legend.var2.pos = "topleft", # Positionnement légende var2
                      legend.var2.horiz = TRUE, # Positionnement horizontal légende var2
                      add = T) # Ajouter au plot

layoutLayer(title = "Évolution démographique 1968 - 1990",
            author = "", scale = 5,  col = "lightgrey",  coltitle = "black")


# 1990 - 2017
plot(st_geometry(com), col = "lightblue", border = NA) # Fond des communes
plot(st_geometry(epci), col = NA, border = "white", lwd = 1, add = T) # Ajouter EPCI

propSymbolsChoroLayer(x = com, var = "POP_2017", var2 = "acc_9017",  # 2 variables
                      col = cols, # couleurs
                      inches = 0.15,  # Taille du cercle maximal...
                      fixmax = max(com$POP_2017), #... Basé sur le max de pop en 2017
                      method = "quantile", # Discrétisation
                      nclass = 6, # Nombre de classes
                      border = "grey50", # Couleur bordures de cercles
                      legend.var.pos = "n", # Ne pas positionner la légende de var1
                      legend.var2.title.txt =  "Taux de croissance moyen annuel (%)", # Titre légende var2
                      legend.var2.values.rnd = 1, # Aronndis valeurs var2
                      legend.var2.pos = "topleft", # Positionnement légende var2
                      legend.var2.horiz = TRUE, # Positionnement horizontal légende var2
                      add = T) # Ajouter au plot

layoutLayer(title = "Évolution démographique 1990 - 2017\net population 2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE, horiz = FALSE, col = "lightgrey", coltitle = "black")

Cartographie thématique

Combiner 2 variables : masse et évolution démographique (sortie graphique)

Cartographie thématique

Lissage par calcul de potentiel (code)

library(potential)
com_pt <- st_centroid(com) # Calcul des centres des communes
pts <- create_grid(x = com, res = 500) # Création d'une grille régulière

# Calcul du potentiel de Stewart
pot <- mcpotential(x = com_pt, # Points d'entrée
                   y = pts, # Points de sortie
                   var = c("POP_1968", "POP_1990", "POP_2017"), # Variables
                   fun = "e", # Ajustement exponentiel
                   span = 2500, # Portée de lissage
                   beta = 2,  # Impédence de la fonction
                   limit = 10000) # Limite de prise en compte des valeurs

pot <- as.data.frame(pot) # Tranformer en dataframe
pts <- cbind(pts, pot) # Affecter le calcul de potentiel à la grille régulière

# Calcul des rapports t+1 / t0 * 100
pts$POP_6890 <- pts$POP_1990 / pts$POP_1968 * 100
pts$POP_9017 <- pts$POP_2017 / pts$POP_1990 * 100

# Discrétisation
brks <- c(74, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 363)

# Création de zones d'égales valeurs
equipot90 <- equipotential(pts, var = "POP_6890", 
                           mask = com, # Masquer les résultats en dehors de la zone d'étude
                           breaks = brks) 

equipot17 <- equipotential(pts, var = "POP_9017", 
                           mask = com,
                           breaks = brks) 

# Cartographie
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,2), bg = "lightgrey") # Ajuster les marges
choroLayer(equipot90, var = "center",
           breaks = brks,
           col = carto.pal(pal1 = "red.pal", pal2 = "green.pal",  n1 = 5, n2 = 6),
           border = NA, legend.pos = "topleft", 
           lwd = .2, legend.title.txt = "Potential de population (1990 / 1968), indice 100",
           legend.horiz = T)

plot(epci$geometry, add = TRUE, col = NA, border = "white")

text(x = 635000, y = 6835000, pos = 4, cex = 0.6,
     labels = "Représentation isoplèthe issue d'un calcul de potentiel\n(fonction exponentielle, distance euclidienne, portée = 5km, beta = 2)")

layoutLayer(title = "Évolution démographique 1968 - 1990",
            author = "", scale = 5,  col = "lightgrey",  coltitle = "black")

choroLayer(equipot17, var = "center",
           breaks = brks,
           col = carto.pal(pal1 = "red.pal", pal2 = "green.pal",  n1 = 5, n2 = 6),
           border = NA, legend.pos = "topleft", 
           lwd = .2, legend.title.txt = "Potential de population (2017 / 1990), indice 100",
           legend.horiz = T)

plot(epci$geometry, add = TRUE, col = NA, border = "white")

layoutLayer(title = "Évolution démographique 1990 - 2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE, horiz = FALSE, col = "lightgrey", coltitle = "black")

Cartographie thématique

Lissage par calcul de potentiel (sortie graphique)

En savoir plus sur la méthode des potentiels ? Consultez la vignette du package !

Transposer son analyse

Cibles reproductibles adaptées des fonctions développées par JD Long

L’intérêt de soigner son code

Un script bien construit et un code bien structuré (préparation des données, noms d’objets, séries d’instructions, commentaires etc.) sur l’ensemble de la chaine de traitement permet :

    • Un travail collaboratif plus aisé.
    • D’impacter aisément la mise à jour d’un jeu de données d’entrée.
    • De transposer son analyse et ses visualisations à d’autres dimensions du jeu de données (espace d’étude, dates de recensement, catégories de population…).

Comme toute analyse, cette phase de consolidation succède généralement à celle d’exploration.

Le script pour la MGP…

Ce script comprend l’intégralité des traitements réalisés

##############################################################
#                                                            #
#   MISE EN PLACE D'UNE CHAINE DE TRAITEMENT REPRODUCTIBLE   #
#                 METROPOLE DU GRAND PARIS                   #
#                                                            #
#         RONAN YSEBAERT (RIATE - UNIV. PARIS, CNRS)         # 
#              SEMINAIRE RUSS - DECEMBRE 2020                #
#                                                            #
##############################################################

##############################################################
# 0 - Sélection de l'espace d'étude
##############################################################
# Ceci afin de transposer l'analyse à d'autres métropoles

# Nom de la métropole (une des métropoles contenues dans le fichier métro
study <- "Grand Paris"

# Code INSEE de la commune "centre". 
moncentre <- "75101"

#################################################
# 1 - Imports des librairies et des données
#################################################
## 1.1 Librairies
# Installer les librairies utiles (si manquantes)
## Séparément 
#install.packages("readxl", repos = "http://cran.us.r-project.org")

## Ou tous ensemble
my_packages <- c( "sf", "cartography", "readxl", "stringr", "dplyr", "ggplot2",
                  "ggthemes", "questionr")
missing_packages <- my_packages[!(my_packages %in% installed.packages()[,"Package"])]
if(length(missing_packages)) install.packages(missing_packages, repos = "http://cran.us.r-project.org")

# Lancement des packages
## Séparément... 
library(readxl)
library(sf)
library(cartography)

## Ou tous ensemble 
lapply(my_packages, library, character.only = TRUE)


## 1.2 Données
# Importer à partir d'un fichier Excel (INSEE) 
data2017 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", skip = 7, sheet = "2017")) 
data1990 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", skip = 7, sheet = "1990")) 
data1968 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", skip = 7, sheet = "1968")) 

# Importer des géométries communales (IGN - Geofla 2016)
com <- st_read(dsn = "input/COMMUNE.shp", stringsAsFactors = F)

# Importer l'appartenance intercommunale des communes de la Métropole du Grand Paris
metro <- read.csv(file = "input/metropoles.csv", 
                  sep = "\t", 
                  colClasses = c(rep("character", 4)),
                  encoding = "latin1")


####################################
# 2 - Manipulation de données
####################################

## 2.1 Jointures
# Jointure : On regroupe les objets data2015, data1990 et data1968 en un seul 
# objet "data"
data <- merge(data2017, data1990, by.x = "COM", by.y = "COM", all.x = T)
data <- merge(data, data1968, by.x = "COM", by.y = "COM", all.x = T) 

# Jointure 2 : EPCI d'appartenance
data <- merge(data, metro, by.x = "COM", by.y = "INSEE_COM", all.x = T)

# Jointure 3 : avec les découpages territoriaux
com <- merge(com, data, by.x = "INSEE_COM", by.y = "COM", all.x = T)


## 2.2 Nettoyage
# Supprimer les objets devenus obsolèthes
rm(data1968, data1990, data2017, metro, data)

# Conserver uniquement certaines variables et les renommer
fields <- c("INSEE_COM", "NOM_COM", "PSDC68", "PSDC90", "PMUN17",
            "METRO_LABEL", "EPCI_SUB", "LIB_EPCI_SUB")
com <- com[,fields]

colnames(com)[1:8] <- c("COM","Nom","POP_1968","POP_1990","POP_2017",
                        "METRO_LABEL", "EPCI_LIB", "EPCI_LABEL")


# Retirer les valeurs égales à 0 ou NA (modifications territoriales)
com <- com[com$POP_1968 != 0 & com$POP_1990 != 0 & com$POP_2017 != 0,]

com <- com[!is.na(com$POP_1968) & !is.na(com$POP_1990) & !is.na(com$POP_2017),]


## 2.3 Créer de nouveaux indicateurs à partir de variables pré-existantes
# Variable quantitative : taux de croissance moyens annuels 1975-1990, 1990-2015 et 1990-2015
com$acc_6890 <- 100 * (exp(log(com$POP_1990 / com$POP_1968) / (1990 - 1968))-1)
com$acc_9017 <- 100 * (exp(log(com$POP_2017 / com$POP_1990) / (2017 - 1990))-1)
com$acc_6817 <- 100 * (exp(log(com$POP_2017 / com$POP_1968) / (2017 - 1968))-1)


# Discrétiser une variable quantitative pour créer des classes quali
com[!is.finite(com$acc_6817),] <- NA
brks <- quantile(com$acc_6817, probs = seq(0,1,0.25), na.rm = TRUE)
brks <- round(brks, 2)

# Typologie basée sur les quartiles
com$typo6817 <- cut(com$acc_6817,
                    breaks = c(brks[1], brks[2], brks[3], brks[4], brks[5]), 
                    include.lowest = TRUE, 
                    labels = c(paste(brks[1], brks[2], sep = " - "),
                               paste(brks[2], brks[3], sep = " - "),
                               paste(brks[3], brks[4], sep = " - "),
                               paste(brks[4], brks[5], sep = " - ")))

# Ne conserver que les unités géographiques désirées
com <- com[com$METRO_LABEL == study,]
com <- com[!is.na(com$METRO_LABEL),]

# Extraire une chaine de caractère : département d'appartenance de chaque commune ? 
library(stringr)
com$DEP <- str_sub(com$COM,  1, 2)


## 2.4 Agréger en fonction d'une modalité
# Avec dplyr : communes > Départements 
library(dplyr)
dep <- com %>% group_by(DEP) %>% summarize(names = first(DEP))

# En Rbase : communes > EPCI (+ somme de la population)
epci <- aggregate(x = com[,c("POP_1968", "POP_1990", "POP_2017")], 
                  by = list(EPCI_LABEL = com$EPCI_LABEL),
                  FUN = sum, na.rm = TRUE)

head(epci)


# 2.5 Analyse spatiale - distance au centre (défini plus haut)
center <- com[com$COM == moncentre,] # Extraction du centre prédéfini
# Distance euclidienne du centre au reste des centroides de l'espace d'étude
com$dist <- as.vector(c(st_distance(st_centroid(x = center), 
                                    st_centroid(x = com))))


## 2.6 Visualisation de l'espace d'étude
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,1)) # Ajuster les marges et nb de graphiques en sortie
plot(st_geometry(com), col = "orange", border = "white", lwd = .2) # Le fond communal
plot(st_geometry(epci), col = NA, border = "black", add = T, lwd = 0.5) # Les EPT de la métro du Grand Paris
plot(st_geometry(dep), col = NA, border = "darkred", add = T, lwd = 1) # Les EPT de la métro du Grand Paris


# 2.7 Visualisation des données associées aux géométries
head(com)


#####################################################################
# 3 - Quelques statistiques et graphiques associés 
#####################################################################

# Pour les graphiques, retirer les géométries des communes et epci
dfcom <- st_set_geometry(com, NULL)
dfepci <- st_set_geometry(epci, NULL)

#  3.1 Statistiques de base - résumé statistique
summary(dfcom)


# 3.2 Visualisations simples associés à des types statistiques  
# 3.2.1 - barplot de la population communale aux 3 recensements avec une boucle
par(mar = c(2,2,2,2), mfrow = c(1,3)) # marge et disposition

for (i in 3:5){ # On prend les colonnes 3 à 5 du dataframe
  popYear <- colnames(dfcom) # récupérer le nom des colonnes
  hist(dfcom[ ,i], # un histogramme
       breaks = 20, # la série statistique est coupée en 20 classes
       main = popYear[i], # titre
       col = "#4286f4", # couleur des barres
       border = "white", # couleur des bordures
       xlab = "Population", # titre de l'axe des X
       ylab = "Fréquence") # titre de l'axe des Y
}

# 3.2.2 - barplot à l'échelle des EPCI, échelle log, ordonné avec mise en forme
dfepci <- dfepci[order(dfepci$POP_2017),] # Ordonner dfepci de façon croissante sur la pop
par(mar = c(6,4,0,0), mfrow = c(1,1))
barplot(t(as.matrix(dfepci[, c("POP_1968","POP_1990", "POP_2017")])), 
        beside = TRUE, # Barres côte à côte par recensment
        names.arg = str_wrap(dfepci$EPCI_LABEL, width = 15), # Couper les longs labels
        legend.text = c("1968", "1990", "2017"), # Légende
        col = c("#6baed6", "#3182bd","#0d4880"), # Couleur des barres
        las = 2, # Légende des absisses à la verticale
        cex.names = 0.7, # Taille des fontes (X)
        cex.axis = 0.7, # Taille des fontes (Y)
        border = NA, # Pas de bordure autour des barres
        ylim = c(10000,max(dfepci$POP_2017)), # Bornes (Y)
        ylab = "Population résidente totale", # titre (Y)
        log = "y") # Echelle logarithmique


# 3.3 - DQuantitatif ~ Quantitatif 
# Tests de corrélation
cor.test(dfcom$dist, dfcom$acc_6890) 
cor.test(dfcom$dist, dfcom$acc_9017)

# Régression linéaire
# En quelques instructions on peut obtenir une régression linéaire ! 
par(mar = c(6,4,0,0), mfrow = c(1,2), bg = "white")
plot(dfcom$dist, dfcom$acc_6890)
abline(lm(acc_6890 ~ dist, data = dfcom), col = "red")
plot(dfcom$dist, dfcom$acc_9017)
abline(lm(acc_9017 ~ dist, data = dfcom), col = "red") 


# On peut aussi utiliser ggplot pour réaliser des graphiques
library(ggplot2)
library(ggthemes)

par(mfrow = c(1,1)) # Un graphique par plot 
ggplot(dfcom, aes(x=dist, y=acc_9017, color = EPCI_LABEL)) + 
  geom_point(aes(size=POP_2017)) +
  xlab("Distance au centre (m)") +
  ylab("Taux d'accroissement moyen annuel (1990-2017, %)") +
  ggtitle("Croissance démographique et distance au centre") +
  theme_wsj() +
  scale_fill_wsj(palette = "colors6", "DEP") + 
  labs(subtitle="Sources: INSEE-IGN, 2020 / Réalisation : R.Ysebaert, 2020") +
  theme(axis.text = element_text(size = 6), 
        axis.title = element_text(size = 10, face = "bold"),
        legend.title=element_text(size=12, face = "bold"),
        legend.text=element_text(size=10),
        plot.subtitle=element_text(size=8),
        plot.title = element_text(size = 15, face = "bold")) 


# 3.4  Qualitatif ~ Qualitatif
# Test du Chi-2
tab <- table(dfcom$typo6817, dfcom$EPCI_LABEL)
tab.chi2 <- chisq.test(tab)
tab.chi2

library(questionr)
ctab <- cprop(tab, total = FALSE)

# Diagramme en bâton
# Couleurs désirées en fonction des niveaux de la typo
coldf <- data.frame(cols = c("#fc8d59", "#a6d96a","#1a9641", "#085721"),
                    typo6817 = levels(dfcom$typo6817))

# Préparation des données
ctab <- ctab[, order(-ctab[1,], -ctab[2,])] # Ordonner selon catégories croissantes

# Gestion des couleurs (associer à des modalités)
cols <- coldf[coldf$typo6817 %in% dfcom[,"typo6817"],]
cols <- cols$cols 

# Diagramme en bâton
par(mar = c(8,4,2,0), mfrow = c(1,1), xpd = TRUE) # Légende en dehors du graphique

barplot(ctab, main = "Situation par rapport à l'évolution démographique d'autres communes françaises (1968 - 2017)",
        las = 2, 
        col = cols, 
        border = NA,
        names.arg = str_wrap(colnames(ctab), width = 15),
        ylab = "% des communes", # titre (Y)
        cex.lab = 0.6,
        cex.axis = 0.6, 
        cex.names = 0.5, 
        cex.main = 0.6)

legend(0, -20, 
       title = "Croissance démographique 1968-2017 (% moyen annuel)\nQuartiles de la distribution des communes françaises",
       legend = levels(dfcom$typo6817),
       horiz = TRUE, # Positionnement horizontal
       fill = cols, # Couleur des caissons de légende 
       cex = 0.5, 
       bty = "n")



# 3.5 - Quantitatif ~ Qualitatif
# Anova et test de Fisher
anova(lm(acc_6817 ~ EPCI_LABEL, data = dfcom))

# Boxplot taux de croissance 1968-2017 ~ EPCI 
# Ordonner sur la population
dfcom$EPCI_LABEL <- with(dfcom, reorder(EPCI_LABEL, acc_6817, mean,
                                        na.rm = TRUE))

# Couleurs
col <- carto.pal(pal1 = "red.pal", n1 = (nlevels(dfcom$EPCI_LABEL)/2), 
                 pal2 = "green.pal",n2 = (nlevels(dfcom$EPCI_LABEL)/2))

# Boxplot
par(mar = c(4,7,1,1))
boxplot(dfcom$acc_6817 ~ dfcom$EPCI_LABEL,
        col = col,
        xlab = "Taux de croissance moyen annuel 1968-2017 (%)",
        names = str_wrap(levels(dfcom$EPCI_LABEL), width = 30),
        ylab = NA,
        horizontal = TRUE,
        varwidth = TRUE,
        range = 1,
        cex.axis = .5, cex.lab = .8,
        outline = TRUE,
        las = 1) 

# Lignes de repère
abline (v = seq(-5,5,.5), col = "darkgrey", lwd = .5, lty = 3)

# Représenter les valeurs moyennes
xi <- tapply(dfcom$acc_6817, dfcom$EPCI_LABEL, mean, na.rm = T)
points(x = xi, y = seq(1,nlevels(dfcom$EPCI_LABEL),1), col="white",pch = 19)


############################################
# 4 - Cartographie
############################################

library(cartography)

#  4.1 Cartographie de discontinuités (choro layer +  discLayer)
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,2), bg = "lightgrey") # Ajuster les marges

cols <- carto.pal(pal1 = "orange.pal", pal2 = "purple.pal",  n1 = 2, n2 = 4)

# 1968 - 1990
choroLayer(x = com, var = "acc_6890", 
           border = NA,
           col = cols, 
           method = "quantile",
           nclass = 6, 
           legend.pos = "bottom", 
           legend.horiz = TRUE,
           legend.values.rnd = 1,
           lwd = 0.5,
           legend.title.txt = "Taux de croissance moyen annuel (%),\nQuantiles")

plot(epci$geometry, border = "white", lwd = .5, add = TRUE)

borders <- getBorders(com)

discLayer(x = borders,
          df = com,
          var = "acc_6890", 
          type = "abs", 
          method="equal", 
          nclass= 3, 
          threshold = 0.25, 
          sizemin = 1, 
          sizemax = 6,
          col="red", 
          legend.values.rnd = 1,
          legend.title.txt = "Discontinuités\nabsolues (max-min)",
          legend.pos = "topleft", 
          add = TRUE)

layoutLayer(title = "Discontinuités démographiques 1968-1990",,
            scale = 5,
            col = "lightgrey",
            coltitle = "black")


# 1990 - 2017
choroLayer(x = com, var = "acc_9017", 
           border = NA,
           col = cols, 
           method = "quantile",
           nclass = 6, 
           legend.pos = "bottom", 
           legend.horiz = TRUE,
           legend.values.rnd = 1,
           lwd = 0.5,
           legend.title.txt = "Taux de croissance moyen annuel (%),\nQuantiles")

plot(epci$geometry, border = "white", lwd = .5, add = TRUE)

discLayer(x = borders,
          df = com,
          var = "acc_9017", 
          type = "abs", 
          method="equal", 
          nclass= 3, 
          threshold = 0.25, 
          sizemin = 1, 
          sizemax = 6,
          col="red", 
          legend.values.rnd = 1,
          legend.title.txt = "Discontinuités\nabsolues (max-min)",
          legend.pos = "topleft", 
          add = TRUE)

layoutLayer(title = "Discontinuités démographiques 1990-2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE,
            horiz = FALSE,
            col = "lightgrey",
            coltitle = "black")


# 4.2 Cercles proportionnels et carte choroplèthe
# 1968 - 1990
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,2), bg = "lightgrey") # Ajuster les marges
plot(st_geometry(com), col = "lightblue", border = NA)
plot(st_geometry(epci), col = NA, border = "white", lwd = 1, add = T)

propSymbolsChoroLayer(x = com, var = "POP_1990", var2 = "acc_6890", 
                      col = cols,
                      inches = 0.15, 
                      method = "quantile",
                      nclass = 6,
                      border = "grey50",
                      fixmax = max(com$POP_2017),
                      legend.var.values.rnd = -2,
                      legend.var.pos = "n", 
                      legend.var2.values.rnd = 1,
                      legend.var2.pos = "bottom",
                      legend.var2.horiz = TRUE,
                      legend.var2.title.txt = 
                        "Taux de croissance moyen annuel (%)",
                      add = T)

layoutLayer(title = "Évolution démographique 1968 - 1990\n et population 1990",
            author = "", scale = 5,  col = "lightgrey",  coltitle = "black")


# 1990 - 2017
plot(st_geometry(com), col = "lightblue", border = NA)
plot(st_geometry(epci), col = NA, border = "white", lwd = 1, add = T)

propSymbolsChoroLayer(x = com, var = "POP_2017", var2 = "acc_9017", 
                      col = cols,
                      inches = 0.15, 
                      method = "quantile",
                      nclass = 6,
                      border = "grey50",
                      fixmax = max(com$POP_2017),
                      legend.var.values.rnd = -2,
                      legend.var.pos = "topleft", 
                      legend.var2.values.rnd = 1,
                      legend.var2.pos = "bottom",
                      legend.var2.horiz = TRUE,
                      legend.var2.title.txt = 
                        "Taux de croissance moyen annuel (%)",
                      legend.var.title.txt = "Population totale",
                      add = T)

layoutLayer(title = "Évolution démographique 1990 - 2017\net population 2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE, horiz = FALSE, col = "lightgrey", coltitle = "black")


#  Carte de potentiel - évolution démographique 1968-1990 / 1990-2017
library(potential)
com_pt <- st_centroid(com) # Calcul des centres des communes
pts <- create_grid(x = com, res = 500) # Création d'une grille régulière

# Calcul du potentiel dans un voisinage de 2.5 km aux trois dates de recensement
# Ajustement exponentiel, beta = 2
pot <- mcpotential(
  x = com_pt, y = pts,
  var = c("POP_1968", "POP_1990", "POP_2017"),
  fun = "e", span = 2500,
  beta = 2, limit = 10000)

pot <- as.data.frame(pot) 
pts <- cbind(pts, pot) # Affecter le calcul de potentiel à la grille régulière

# Calcul des rapports t+1 / t0 * 100
pts$POP_6890 <- pts$POP_1990 / pts$POP_1968 * 100
pts$POP_9017 <- pts$POP_2017 / pts$POP_1990 * 100

# Discrétisation
brks <- c(74, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 363)

# Création de zones d'égales valeurs
equipot90 <- equipotential(pts, var = "POP_6890", mask = com,
                           breaks = brks) 

equipot17 <- equipotential(pts, var = "POP_9017", mask = com,
                           breaks = brks) 

# Cartographie
choroLayer(equipot90, var = "center",
           breaks = brks,
           col = carto.pal(pal1 = "red.pal", pal2 = "green.pal",  n1 = 5, n2 = 6),
           border = NA, legend.pos = "bottom", 
           lwd = .2, legend.title.txt = "Potential de population (1990 / 1968), indice 100",
           legend.horiz = T)

plot(epci$geometry, add = TRUE, col = NA, border = "white")

layoutLayer(title = "Évolution démographique 1968 - 1990",
            author = "", scale = 5,  col = "lightgrey",  coltitle = "black")

choroLayer(equipot17, var = "center",
           breaks = brks,
           col = carto.pal(pal1 = "red.pal", pal2 = "green.pal",  n1 = 5, n2 = 6),
           border = NA, legend.pos = "bottom", 
           lwd = .2, legend.title.txt = "Potential de population (2017 / 1990), indice 100",
           legend.horiz = T)

plot(epci$geometry, add = TRUE, col = NA, border = "white")

layoutLayer(title = "Évolution démographique 1990 - 2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE, horiz = FALSE, col = "lightgrey", coltitle = "black")
# ------------------------------------------------------------ #

… Et son adaptation marseillaise !

Vous pouvez jouer aux “7 différences”, mais seule la partie 0 est mise à jour !

##############################################################
#                                                            #
#   MISE EN PLACE D'UNE CHAINE DE TRAITEMENT REPRODUCTIBLE   #
#        TRANSPOSITION A LA METROPOLE D'AIX-MARSEILLE        #
#                                                            #
#         RONAN YSEBAERT (RIATE - UNIV. PARIS, CNRS)         # 
#              SEMINAIRE RUSS - DECEMBRE 2020                #
#                                                            #
##############################################################

##############################################################
# 0 - Sélection de l'espace d'étude
##############################################################
# Ceci afin de transposer l'analyse à d'autres métropoles

# Nom de la métropole (une des métropoles contenues dans le fichier métro
study <- "Marseille"

# Code INSEE de la commune "centre". 
moncentre <- "13055"

#################################################
# 1 - Imports des librairies et des données
#################################################
## 1.1 Librairies
# Installer les librairies utiles (si manquantes)
## Séparément 
#install.packages("readxl", repos = "http://cran.us.r-project.org")

## Ou tous ensemble
my_packages <- c( "sf", "cartography", "readxl", "stringr", "dplyr", "ggplot2",
                  "ggthemes", "questionr")
missing_packages <- my_packages[!(my_packages %in% installed.packages()[,"Package"])]
if(length(missing_packages)) install.packages(missing_packages, repos = "http://cran.us.r-project.org")

# Lancement des packages
## Séparément... 
library(readxl)
library(sf)
library(cartography)

## Ou tous ensemble 
lapply(my_packages, library, character.only = TRUE)


## 1.2 Données
# Importer à partir d'un fichier Excel (INSEE) 
data2017 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", skip = 7, sheet = "2017")) 
data1990 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", skip = 7, sheet = "1990")) 
data1968 <- data.frame(read_excel("input/Fichier_poplegale_6817.xls", skip = 7, sheet = "1968")) 

# Importer des géométries communales (IGN - Geofla 2016)
com <- st_read(dsn = "input/COMMUNE.shp", stringsAsFactors = F)

# Importer l'appartenance intercommunale des communes de la Métropole du Grand Paris
metro <- read.csv(file = "input/metropoles.csv", 
                  sep = "\t", 
                  colClasses = c(rep("character", 4)),
                  encoding = "latin1")


####################################
# 2 - Manipulation de données
####################################

## 2.1 Jointures
# Jointure : On regroupe les objets data2015, data1990 et data1968 en un seul 
# objet "data"
data <- merge(data2017, data1990, by.x = "COM", by.y = "COM", all.x = T)
data <- merge(data, data1968, by.x = "COM", by.y = "COM", all.x = T) 

# Jointure 2 : EPCI d'appartenance
data <- merge(data, metro, by.x = "COM", by.y = "INSEE_COM", all.x = T)

# Jointure 3 : avec les découpages territoriaux
com <- merge(com, data, by.x = "INSEE_COM", by.y = "COM", all.x = T)


## 2.2 Nettoyage
# Supprimer les objets devenus obsolèthes
rm(data1968, data1990, data2017, metro, data)

# Conserver uniquement certaines variables et les renommer
fields <- c("INSEE_COM", "NOM_COM", "PSDC68", "PSDC90", "PMUN17",
            "METRO_LABEL", "EPCI_SUB", "LIB_EPCI_SUB")
com <- com[,fields]

colnames(com)[1:8] <- c("COM","Nom","POP_1968","POP_1990","POP_2017",
                        "METRO_LABEL", "EPCI_LIB", "EPCI_LABEL")


# Retirer les valeurs égales à 0 ou NA (modifications territoriales)
com <- com[com$POP_1968 != 0 & com$POP_1990 != 0 & com$POP_2017 != 0,]

com <- com[!is.na(com$POP_1968) & !is.na(com$POP_1990) & !is.na(com$POP_2017),]


## 2.3 Créer de nouveaux indicateurs à partir de variables pré-existantes
# Variable quantitative : taux de croissance moyens annuels 1975-1990, 1990-2015 et 1990-2015
com$acc_6890 <- 100 * (exp(log(com$POP_1990 / com$POP_1968) / (1990 - 1968))-1)
com$acc_9017 <- 100 * (exp(log(com$POP_2017 / com$POP_1990) / (2017 - 1990))-1)
com$acc_6817 <- 100 * (exp(log(com$POP_2017 / com$POP_1968) / (2017 - 1968))-1)


# Discrétiser une variable quantitative pour créer des classes quali
com[!is.finite(com$acc_6817),] <- NA
brks <- quantile(com$acc_6817, probs = seq(0,1,0.25), na.rm = TRUE)
brks <- round(brks, 2)

# Typologie basée sur les quartiles
com$typo6817 <- cut(com$acc_6817,
                    breaks = c(brks[1], brks[2], brks[3], brks[4], brks[5]), 
                    include.lowest = TRUE, 
                    labels = c(paste(brks[1], brks[2], sep = " - "),
                               paste(brks[2], brks[3], sep = " - "),
                               paste(brks[3], brks[4], sep = " - "),
                               paste(brks[4], brks[5], sep = " - ")))

# Ne conserver que les unités géographiques désirées
com <- com[com$METRO_LABEL == study,]
com <- com[!is.na(com$METRO_LABEL),]

# Extraire une chaine de caractère : département d'appartenance de chaque commune ? 
library(stringr)
com$DEP <- str_sub(com$COM,  1, 2)


## 2.4 Agréger en fonction d'une modalité
# Avec dplyr : communes > Départements 
library(dplyr)
dep <- com %>% group_by(DEP) %>% summarize(names = first(DEP))

# En Rbase : communes > EPCI (+ somme de la population)
epci <- aggregate(x = com[,c("POP_1968", "POP_1990", "POP_2017")], 
                  by = list(EPCI_LABEL = com$EPCI_LABEL),
                  FUN = sum, na.rm = TRUE)

head(epci)


# 2.5 Analyse spatiale - distance au centre (défini plus haut)
center <- com[com$COM == moncentre,] # Extraction du centre prédéfini
# Distance euclidienne du centre au reste des centroides de l'espace d'étude
com$dist <- as.vector(c(st_distance(st_centroid(x = center), 
                                    st_centroid(x = com))))


## 2.6 Visualisation de l'espace d'étude
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,1)) # Ajuster les marges et nb de graphiques en sortie
plot(st_geometry(com), col = "orange", border = "white", lwd = .2) # Le fond communal
plot(st_geometry(epci), col = NA, border = "black", add = T, lwd = 0.5) # Les EPT de la métro du Grand Paris
plot(st_geometry(dep), col = NA, border = "darkred", add = T, lwd = 1) # Les EPT de la métro du Grand Paris


# 2.7 Visualisation des données associées aux géométries
head(com)


#####################################################################
# 3 - Quelques statistiques et graphiques associés 
#####################################################################

# Pour les graphiques, retirer les géométries des communes et epci
dfcom <- st_set_geometry(com, NULL)
dfepci <- st_set_geometry(epci, NULL)

#  3.1 Statistiques de base - résumé statistique
summary(dfcom)


# 3.2 Visualisations simples associés à des types statistiques  
# 3.2.1 - barplot de la population communale aux 3 recensements avec une boucle
par(mar = c(2,2,2,2), mfrow = c(1,3)) # marge et disposition

for (i in 3:5){ # On prend les colonnes 3 à 5 du dataframe
  popYear <- colnames(dfcom) # récupérer le nom des colonnes
  hist(dfcom[ ,i], # un histogramme
       breaks = 20, # la série statistique est coupée en 20 classes
       main = popYear[i], # titre
       col = "#4286f4", # couleur des barres
       border = "white", # couleur des bordures
       xlab = "Population", # titre de l'axe des X
       ylab = "Fréquence") # titre de l'axe des Y
}

# 3.2.2 - barplot à l'échelle des EPCI, échelle log, ordonné avec mise en forme
dfepci <- dfepci[order(dfepci$POP_2017),] # Ordonner dfepci de façon croissante sur la pop
par(mar = c(6,4,0,0), mfrow = c(1,1))
barplot(t(as.matrix(dfepci[, c("POP_1968","POP_1990", "POP_2017")])), 
        beside = TRUE, # Barres côte à côte par recensment
        names.arg = str_wrap(dfepci$EPCI_LABEL, width = 15), # Couper les longs labels
        legend.text = c("1968", "1990", "2017"), # Légende
        col = c("#6baed6", "#3182bd","#0d4880"), # Couleur des barres
        las = 2, # Légende des absisses à la verticale
        cex.names = 0.7, # Taille des fontes (X)
        cex.axis = 0.7, # Taille des fontes (Y)
        border = NA, # Pas de bordure autour des barres
        ylim = c(10000,max(dfepci$POP_2017)), # Bornes (Y)
        ylab = "Population résidente totale", # titre (Y)
        log = "y") # Echelle logarithmique


# 3.3 - DQuantitatif ~ Quantitatif 
# Tests de corrélation
cor.test(dfcom$dist, dfcom$acc_6890) 
cor.test(dfcom$dist, dfcom$acc_9017)

# Régression linéaire
# En quelques instructions on peut obtenir une régression linéaire ! 
par(mar = c(6,4,0,0), mfrow = c(1,2), bg = "white")
plot(dfcom$dist, dfcom$acc_6890)
abline(lm(acc_6890 ~ dist, data = dfcom), col = "red")
plot(dfcom$dist, dfcom$acc_9017)
abline(lm(acc_9017 ~ dist, data = dfcom), col = "red") 


# On peut aussi utiliser ggplot pour réaliser des graphiques
library(ggplot2)
library(ggthemes)

par(mfrow = c(1,1)) # Un graphique par plot 
ggplot(dfcom, aes(x=dist, y=acc_9017, color = EPCI_LABEL )) + 
  geom_point(aes(size=POP_2017)) +
  xlab("Distance au centre (m)") +
  ylab("Taux d'accroissement moyen annuel (1990-2017, %)") +
  ggtitle("Croissance démographique et distance au centre") +
  theme_wsj() +
  scale_fill_wsj(palette = "colors6", "DEP") + 
  labs(subtitle="Sources: INSEE-IGN, 2020 / Réalisation : R.Ysebaert, 2020") +
  theme(axis.text = element_text(size = 6), 
        axis.title = element_text(size = 10, face = "bold"),
        legend.title=element_text(size=12, face = "bold"),
        legend.text=element_text(size=10),
        plot.subtitle=element_text(size=8),
        plot.title = element_text(size = 15, face = "bold")) 


# 3.4  Qualitatif ~ Qualitatif
# Test du Chi-2
tab <- table(dfcom$typo6817, dfcom$EPCI_LABEL)
tab.chi2 <- chisq.test(tab)
tab.chi2

library(questionr)
ctab <- cprop(tab, total = FALSE)

# Diagramme en bâton
# Couleurs désirées en fonction des niveaux de la typo
coldf <- data.frame(cols = c("#fc8d59", "#a6d96a","#1a9641", "#085721"),
                    typo6817 = levels(dfcom$typo6817))

# Préparation des données
ctab <- ctab[, order(-ctab[1,], -ctab[2,])] # Ordonner selon catégories croissantes

# Gestion des couleurs (associer à des modalités)
cols <- coldf[coldf$typo6817 %in% dfcom[,"typo6817"],]
cols <- cols$cols 

# Diagramme en bâton
par(mar = c(8,4,2,0), mfrow = c(1,1), xpd = TRUE) # Légende en dehors du graphique

barplot(ctab, main = "Situation par rapport à l'évolution démographique d'autres communes françaises (1968 - 2017)",
        las = 2, 
        col = cols, 
        border = NA,
        names.arg = str_wrap(colnames(ctab), width = 15),
        ylab = "% des communes", # titre (Y)
        cex.lab = 0.6,
        cex.axis = 0.6, 
        cex.names = 0.5, 
        cex.main = 0.6)

legend(0, -20, 
       title = "Croissance démographique 1968-2017 (% moyen annuel)\nQuartiles de la distribution des communes françaises",
       legend = levels(dfcom$typo6817),
       horiz = TRUE, # Positionnement horizontal
       fill = cols, # Couleur des caissons de légende 
       cex = 0.5, 
       bty = "n")



# 3.5 - Quantitatif ~ Qualitatif
# Anova et test de Fisher
anova(lm(acc_6817 ~ EPCI_LABEL, data = dfcom))

# Boxplot taux de croissance 1968-2017 ~ EPCI 
# Ordonner sur la population
dfcom$EPCI_LABEL <- with(dfcom, reorder(EPCI_LABEL, acc_6817, mean,
                                        na.rm = TRUE))

# Couleurs
col <- carto.pal(pal1 = "red.pal", n1 = (nlevels(dfcom$EPCI_LABEL)/2), 
                 pal2 = "green.pal",n2 = (nlevels(dfcom$EPCI_LABEL)/2))

# Boxplot
par(mar = c(4,7,1,1))
boxplot(dfcom$acc_6817 ~ dfcom$EPCI_LABEL,
        col = col,
        xlab = "Taux de croissance moyen annuel 1968-2017 (%)",
        names = str_wrap(levels(dfcom$EPCI_LABEL), width = 30),
        ylab = NA,
        horizontal = TRUE,
        varwidth = TRUE,
        range = 1,
        cex.axis = .5, cex.lab = .8,
        outline = TRUE,
        las = 1) 

# Lignes de repère
abline (v = seq(-5,5,.5), col = "darkgrey", lwd = .5, lty = 3)

# Représenter les valeurs moyennes
xi <- tapply(dfcom$acc_6817, dfcom$EPCI_LABEL, mean, na.rm = T)
points(x = xi, y = seq(1,nlevels(dfcom$EPCI_LABEL),1), col="white",pch = 19)


############################################
# 4 - Cartographie
############################################

library(cartography)

#  4.1 Cartographie de discontinuités (choro layer +  discLayer)
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,2), bg = "lightgrey") # Ajuster les marges

cols <- carto.pal(pal1 = "orange.pal", pal2 = "purple.pal",  n1 = 2, n2 = 4)

# 1968 - 1990
choroLayer(x = com, var = "acc_6890", 
           border = NA,
           col = cols, 
           method = "quantile",
           nclass = 6, 
           legend.pos = "bottom", 
           legend.horiz = TRUE,
           legend.values.rnd = 1,
           lwd = 0.5,
           legend.title.txt = "Taux de croissance moyen annuel (%),\nQuantiles")

plot(epci$geometry, border = "white", lwd = .5, add = TRUE)

borders <- getBorders(com)

discLayer(x = borders,
          df = com,
          var = "acc_6890", 
          type = "abs", 
          method="equal", 
          nclass= 3, 
          threshold = 0.25, 
          sizemin = 1, 
          sizemax = 6,
          col="red", 
          legend.values.rnd = 1,
          legend.title.txt = "Discontinuités\nabsolues (max-min)",
          legend.pos = "topleft", 
          add = TRUE)

layoutLayer(title = "Discontinuités démographiques 1968-1990",,
            scale = 5,
            col = "lightgrey",
            coltitle = "black")


# 1990 - 2017
choroLayer(x = com, var = "acc_9017", 
           border = NA,
           col = cols, 
           method = "quantile",
           nclass = 6, 
           legend.pos = "bottom", 
           legend.horiz = TRUE,
           legend.values.rnd = 1,
           lwd = 0.5,
           legend.title.txt = "Taux de croissance moyen annuel (%),\nQuantiles")

plot(epci$geometry, border = "white", lwd = .5, add = TRUE)

discLayer(x = borders,
          df = com,
          var = "acc_9017", 
          type = "abs", 
          method="equal", 
          nclass= 3, 
          threshold = 0.25, 
          sizemin = 1, 
          sizemax = 6,
          col="red", 
          legend.values.rnd = 1,
          legend.title.txt = "Discontinuités\nabsolues (max-min)",
          legend.pos = "topleft", 
          add = TRUE)

layoutLayer(title = "Discontinuités démographiques 1990-2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE,
            horiz = FALSE,
            col = "lightgrey",
            coltitle = "black")


# 4.2 Cercles proportionnels et carte choroplèthe
# 1968 - 1990
par(mar = c(0.5,0.5,1.5,0.5), mfrow = c(1,2), bg = "lightgrey") # Ajuster les marges
plot(st_geometry(com), col = "lightblue", border = NA)
plot(st_geometry(epci), col = NA, border = "white", lwd = 1, add = T)

propSymbolsChoroLayer(x = com, var = "POP_1990", var2 = "acc_6890", 
                      col = cols,
                      inches = 0.15, 
                      method = "quantile",
                      nclass = 6,
                      border = "grey50",
                      fixmax = max(com$POP_2017),
                      legend.var.values.rnd = -2,
                      legend.var.pos = "n", 
                      legend.var2.values.rnd = 1,
                      legend.var2.pos = "bottom",
                      legend.var2.horiz = TRUE,
                      legend.var2.title.txt = 
                        "Taux de croissance moyen annuel (%)",
                      add = T)

layoutLayer(title = "Évolution démographique 1968 - 1990\n et population 1990",
            author = "", scale = 5,  col = "lightgrey",  coltitle = "black")


# 1990 - 2017
plot(st_geometry(com), col = "lightblue", border = NA)
plot(st_geometry(epci), col = NA, border = "white", lwd = 1, add = T)

propSymbolsChoroLayer(x = com, var = "POP_2017", var2 = "acc_9017", 
                      col = cols,
                      inches = 0.15, 
                      method = "quantile",
                      nclass = 6,
                      border = "grey50",
                      fixmax = max(com$POP_2017),
                      legend.var.values.rnd = -2,
                      legend.var.pos = "topleft", 
                      legend.var2.values.rnd = 1,
                      legend.var2.pos = "bottom",
                      legend.var2.horiz = TRUE,
                      legend.var2.title.txt = 
                        "Taux de croissance moyen annuel (%)",
                      legend.var.title.txt = "Population totale",
                      add = T)

layoutLayer(title = "Évolution démographique 1990 - 2017\net population 2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE, horiz = FALSE, col = "lightgrey", coltitle = "black")


#  Carte de potentiel - évolution démographique 1968-1990 / 1990-2017
library(potential)
com_pt <- st_centroid(com) # Calcul des centres des communes
pts <- create_grid(x = com, res = 500) # Création d'une grille régulière

# Calcul du potentiel dans un voisinage de 2.5 km aux trois dates de recensement
# Ajustement exponentiel, beta = 2
pot <- mcpotential(
  x = com_pt, y = pts,
  var = c("POP_1968", "POP_1990", "POP_2017"),
  fun = "e", span = 2500,
  beta = 2, limit = 10000)

pot <- as.data.frame(pot) 
pts <- cbind(pts, pot) # Affecter le calcul de potentiel à la grille régulière

# Calcul des rapports t+1 / t0 * 100
pts$POP_6890 <- pts$POP_1990 / pts$POP_1968 * 100
pts$POP_9017 <- pts$POP_2017 / pts$POP_1990 * 100

# Discrétisation
brks <- c(74, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 363)

# Création de zones d'égales valeurs
equipot90 <- equipotential(pts, var = "POP_6890", mask = com,
                           breaks = brks) 

equipot17 <- equipotential(pts, var = "POP_9017", mask = com,
                           breaks = brks) 

# Cartographie
choroLayer(equipot90, var = "center",
           breaks = brks,
           col = carto.pal(pal1 = "red.pal", pal2 = "green.pal",  n1 = 5, n2 = 6),
           border = NA, legend.pos = "bottom", 
           lwd = .2, legend.title.txt = "Potential de population (1990 / 1968), indice 100",
           legend.horiz = T)

plot(epci$geometry, add = TRUE, col = NA, border = "white")

layoutLayer(title = "Évolution démographique 1968 - 1990",
            author = "", scale = 5,  col = "lightgrey",  coltitle = "black")

choroLayer(equipot17, var = "center",
           breaks = brks,
           col = carto.pal(pal1 = "red.pal", pal2 = "green.pal",  n1 = 5, n2 = 6),
           border = NA, legend.pos = "bottom", 
           lwd = .2, legend.title.txt = "Potential de population (2017 / 1990), indice 100",
           legend.horiz = T)

plot(epci$geometry, add = TRUE, col = NA, border = "white")

layoutLayer(title = "Évolution démographique 1990 - 2017",
            author = "Ronan Ysebaert, 2020", sources = "INSEE, IGN, 2020",
            scale = FALSE, horiz = FALSE, col = "lightgrey", coltitle = "black")
# ------------------------------------------------------------ #

… Qui fonctionne !

Evolution démographique par EPCI 1968-2017

Effet de la distance au centre sur la croissance démographique

… Qui fonctionne !

Evolution démographique 1968-2017 des communes de chaque EPCI au regard des quartiles de la distribution des communes françaises

Dispersion statistique de la croissance démographique par EPCI

… Qui fonctionne !

Évolution démographique communale comparée (1968 - 1990 et 1990-2017) et discontinuités territoriales

… Qui fonctionne !

Évolution démographique communale comparée (1968 - 1990 et 1990-2017), mise en relation des potentiels de poulation (voisinage 5 km)

Documenter son code

Literate programming

Literate programming : paradigme de programmation où est associé au programme informatique une explication de sa logique dans un mode d’expression naturel.

    • Un processus d’explication de la pensée : Expliquer à des êtres humains ce que nous voulons que l’ordinateur fasse (Knuth, 1970).

    • Des programmes qui gagnent en qualité : une explication intintérompue de la logique, à la manière d’un essai.

    • Rend les mauvaises décisions de conception plus évidentes.

    • Permet de reprendre des analyses ultérieurement sans se perdre, même longtemps après.

Un gain général en reproductiblité de la démarche et en ouverture des méthodes scientifiques.

R Mardown

Dans RStudio, on peut aisément importer un modèle R Markdown pour associer texte et code.

R Mardown

Le document est généré en “tricotant” (knit) le fichier .Rmd

Plusieurs formats de sortie sont disponibles (presentations, rapports, site Web, ouvrages…)

Un R Markdown de synthèse ?

Le résumé commenté de cette séance et de la chaine de traitement mobilisée est disponible ici.

Il a été généré à partir de ce fichier .Rmd.

Reproduisez l’analyse !


Téléchargez les données

https://frama.link/RUSS_IntroR_data

Ouvrez le projet et le script

Il est nécessaire d’installer prélablement R (Linux, Windows ou Mac) puis RStudio

Exécutez le script

Exécuter une ligne de commande : CTRL + ENTRÉE

Comment a été créé le rapport ?

Le fichier report.rmd génère le fichier report.html. A un stade plus avancé, il est possible d’adapter la mise en forme en utilisant des modèles (templates) ou en paramétrant un .css

Quelques conseils

    • L’apprentissage de R reste coûteux, il faut y consacrer du temps…

    • Ne pas être effrayé à la vue d’un warning ou d’une erreur, généralement liés à des erreurs de syntaxe (virgule, casse, “orthographe”).

    • On n’est pas seuls ! Consultez les forums, comme stack overflow. Préférez les recherches en anglais (make barplot with R).

    • Exécuter les exemples reproductibles proposés par des blogs, données d’exemple des packages, voire des vignette. Bien regarder la nature des objets utilisés en entrées (class(x)), le type statistique des variables (str(x)) et à quoi ressemblent les données (head(x)) pour comprendre leur utilisation-type.

A vous de jouer !

Remerciements


Timothée Giraud (CNRS)

Hugues Pécout (CNRS)