Objectif :

Appréhender une chaîne de traitement de données géographiques simples. Faire des traitements statistiques univariés, bivariés, et la cartographie choroplèthe associée.

Avec un volume horaire de 3 demi-journées (12), l’objectif est centré sur l’écriture d’un programme R puis un document Rmarkdown et d’effectuer des traitements de statistique univariée (valeurs centrales, paramètres de dispersion, histogramme, …) et bivariée (par exemple, corrélation, régression, analyse de variance, test d’égalité des moyennes, tableau de contingence, chi-2). Nous réaliserons ceci en mettant en oeuvre l’ensemble de la chaîne de traitement des données, des tableaux de données jusqu’à la cartographie univariée choroplèthe simple. Afin de vous assurer un accompagnement vers le stage commun avec Carthageo, nous utiliserons à la fois le langage standard R, et les nouvelles évolutions tant pour l’écriture du code (tidyverse) que pour les graphiques (ggplot2). L’objectif est que vous disposiez d’un kit méthodologique pour pouvoir traiter de l’univarié et du bivarié en autonomie, jusqu’à des opérations simples de cartographie.

Support de cours :

Le programme de cet enseignement sera basé sur les outils pédagogiques suivants :

• d’une part les ressources mises à disposition par Claude Grasland sur l’initiation au langage R : https://claudegrasland.github.io/startR/ et celles de Julien Barnier dans son Introduction à R et au Tidyverse : https://juba.github.io/tidyverse/

• d’autre part les ressources de l’UAR RIATE, sur les bases du traitement géomatique avec R : https://rcarto.github.io/geomatique_avec_r/ et https://rcarto.github.io/cartographie_avec_r/

• Nous recommandons, par ailleurs cet ouvrage: https://pur-editions.fr/product/8090/r-pour-la-statistique-et-la-science-des-donnees (une commande collective est en cours pour un fournir un exemplaire aux étudiants qui suivent le stage, nous vous tiendrons informé.es)

• Le Groupe ElementR, commun aux UMR Géographie-cités, UMR PRODIG et UAR RIATE organise la diffusion de modules d’initation et de formation à R, au tidyverse, à la data visualisation et la cartographie : https://elementr.gitpages.huma-num.fr/website/programme.html

Dates et lieu : Les cours auront lieu en salle 310 du bâtiment Olympe de Gouges à l’Université Paris Cité (MAJ Sept 2O25):

• Lundi 15 septembre 2025 : 9h-12 et 13-17h, ODG 310
• Jeudi 18 septembre 2025 ; 9h-12h, ODG 310

Matériel : Les étudiants sont invités à utiliser leur propre ordinateur, et à y installer les dernières versions de R et Rstudio AVANT la première séance. Nous serons en salle informatique et des postes de travail individuels seront également disponibles - https://larmarange.github.io/analyse-R/installation-de-R-et-RStudio.html

Prérequis : Notion de base en statistique univariée et, si possible bivariée.

Documents et syllabus du TD :

Syllabus et données : http://rlegoix.gitpages.huma-num.fr/cours_m2_traitementsstatscarto

L’ensemble des documents sont mis à dispostion sur un dépôt Git : https://gitlab.huma-num.fr/rlegoix/cours_m2_traitementsstatscarto

R: une plateforme statistique, un (des) langage(s) de programmation, un couteau suisse.

De la chaîne de traitement de l’information géographique et la production de documents

La chaîne de traitements d’après C. Grasland

Un éconsystème, ou un environnement, qui évolue vite, en recourrant à des librairies (packages) pour accéder à des fonctions avancées ou plus simples d’utilisation.

Packages et environnements

Objectif: Découvrir l’environnement. Manipuler les principaux objects R.

Travail sur un tableau de données (dataframe)

Importer les données à partir d’une fichier texte (format CSV)

Données utilisées : Données Revenus, pauvreté et niveau de vie en 2021 (Iris) de la base Dispositif Fichier localisé social et fiscal (Filosofi) : https://www.insee.fr/fr/statistiques/8229323/ La création de ce fichier est documentée dans les annexes.

On peut importer directement le fichier, ou contrôler la structure des variables, comme sur la copie d’écran ci-dessous.

Importer un CSV en utilisant l’interface graphique

Notre objectif est de travailler avec les variables

Le tableau de données finales (dataframe) comportera cinq colonnes :

• identifiant géographique : CODE_IRIS
• code département : DEP
• revenus médians 2021 (Médiane du revenu déclaré par unité de consommation (en euros)) : DEC_MED21
• indice de Gini du revenu déclaré par unité de consommation 2021 : DEC_GI21
• part des autres revenus (%) (essentiellement des revenus du patrimoine) : DEC_PAUT21

# Méthode simple
Revenus_IRIS <- read.csv("../public/data/Revenus_IRIS_PC.csv")
## Cette méthode nécessite de contrôler à postériori la structure des variables

# Méthode en utilisant l'interface RStudio qui permet 1. de contrôler la structure des variables que l'on importe ; 2. de copier / coller le code utilisé.
library(readr)
Revenus_IRIS_PC <- read_csv("../public/data/Revenus_IRIS_PC.csv", 
    col_types = cols(IRIS = col_character(), 
    DEC_MED21 = col_number(), DEC_GI21 = col_number()))

## New names:
## • `` -> `...1`

## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
##   dat <- vroom(...)
##   problems(dat)

# Les dimensions du tableau
dim(Revenus_IRIS_PC)

## [1] 2745   28

#> [1] 2529    5

#La fonction class() fournit le type d’un tableau
class(Revenus_IRIS_PC)

## [1] "spec_tbl_df" "tbl_df"      "tbl"         "data.frame"

#> [1] "data.frame"

# La fonction head() fournit les premières lignes`
# Visualiser les colonnes et variables disponibles
head(Revenus_IRIS_PC)

## # A tibble: 6 × 28
##    ...1 IRIS      DEC_PIMP21 DEC_TP6021 DEC_INCERT21 DEC_Q121 DEC_MED21 DEC_Q321
##   <dbl> <chr>     <chr>      <chr>      <chr>        <chr>        <dbl> <chr>   
## 1     1 751010101 ns         ns         so           ns              NA ns      
## 2     2 751010102 ns         ns         so           ns              NA ns      
## 3     3 751010103 ns         ns         so           ns              NA ns      
## 4     4 751010104 ns         ns         so           ns              NA ns      
## 5     5 751010105 nd         nd         so           nd              NA nd      
## 6     6 751010199 nd         nd         so           nd              NA nd      
## # ℹ 20 more variables: DEC_EQ21 <chr>, DEC_D121 <chr>, DEC_D221 <chr>,
## #   DEC_D321 <chr>, DEC_D421 <chr>, DEC_D621 <chr>, DEC_D721 <chr>,
## #   DEC_D821 <chr>, DEC_D921 <chr>, DEC_RD21 <chr>, DEC_S80S2021 <chr>,
## #   DEC_GI21 <dbl>, DEC_PACT21 <chr>, DEC_PTSA21 <chr>, DEC_PCHO21 <chr>,
## #   DEC_PBEN21 <chr>, DEC_PPEN21 <chr>, DEC_PAUT21 <chr>, DEC_NOTE21 <chr>,
## #   DEP <dbl>

# Vérifier le type des variables avec str()
str(Revenus_IRIS_PC)

## spc_tbl_ [2,745 × 28] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ ...1        : num [1:2745] 1 2 3 4 5 6 7 8 9 10 ...
##  $ IRIS        : chr [1:2745] "751010101" "751010102" "751010103" "751010104" ...
##  $ DEC_PIMP21  : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_TP6021  : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_INCERT21: chr [1:2745] "so" "so" "so" "so" ...
##  $ DEC_Q121    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_MED21   : num [1:2745] NA NA NA NA NA ...
##  $ DEC_Q321    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_EQ21    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_D121    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_D221    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_D321    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_D421    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_D621    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_D721    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_D821    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_D921    : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_RD21    : chr [1:2745] "so" "so" "so" "so" ...
##  $ DEC_S80S2021: chr [1:2745] "so" "so" "so" "so" ...
##  $ DEC_GI21    : num [1:2745] NA NA NA NA NA NA 471 531 503 487 ...
##  $ DEC_PACT21  : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_PTSA21  : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_PCHO21  : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_PBEN21  : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_PPEN21  : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_PAUT21  : chr [1:2745] "ns" "ns" "ns" "ns" ...
##  $ DEC_NOTE21  : chr [1:2745] "so" "so" "so" "so" ...
##  $ DEP         : num [1:2745] 75 75 75 75 75 75 75 75 75 75 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   ...1 = col_double(),
##   ..   IRIS = col_character(),
##   ..   DEC_PIMP21 = col_character(),
##   ..   DEC_TP6021 = col_character(),
##   ..   DEC_INCERT21 = col_character(),
##   ..   DEC_Q121 = col_character(),
##   ..   DEC_MED21 = col_number(),
##   ..   DEC_Q321 = col_character(),
##   ..   DEC_EQ21 = col_character(),
##   ..   DEC_D121 = col_character(),
##   ..   DEC_D221 = col_character(),
##   ..   DEC_D321 = col_character(),
##   ..   DEC_D421 = col_character(),
##   ..   DEC_D621 = col_character(),
##   ..   DEC_D721 = col_character(),
##   ..   DEC_D821 = col_character(),
##   ..   DEC_D921 = col_character(),
##   ..   DEC_RD21 = col_character(),
##   ..   DEC_S80S2021 = col_character(),
##   ..   DEC_GI21 = col_number(),
##   ..   DEC_PACT21 = col_character(),
##   ..   DEC_PTSA21 = col_character(),
##   ..   DEC_PCHO21 = col_character(),
##   ..   DEC_PBEN21 = col_character(),
##   ..   DEC_PPEN21 = col_character(),
##   ..   DEC_PAUT21 = col_character(),
##   ..   DEC_NOTE21 = col_character(),
##   ..   DEP = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

# Sélectionner les colonnes qui nous intéressent (nettoyage des fichiers)
Revenus_IRIS_PC <- Revenus_IRIS_PC[,c("IRIS","DEP","DEC_MED21","DEC_GI21","DEC_PAUT21")]

# Nettoyer les fichiers. On a vu que des lignes comportaient des données numériques manquantes, codées NA. On va les supprimer. Il y a deux solutions

# On supprime les enregistrements qui contiennent des données manquantes
# Supprimer les lignes avec des valeurs NA
Revenus_IRIS_PC <- na.omit(Revenus_IRIS_PC)
# Il doit rester 2529 lignes

# On conserve les enregistrements qui sont complets
Revenus_IRIS_PC <- Revenus_IRIS_PC[complete.cases(Revenus_IRIS_PC), ]
# Il doit rester 2529 lignes (le monde est bien fait)
# 

# Ooops, il semblerait que la variable `DEC_PAUT21` soit toujours une chaîne de caractère. On a besoin de la transformer en variable numérique 
# Ce n'est pas très grave
# Essayer puis voir ce que ça donne.

Revenus_IRIS_PC$RevPatrimoine <- as.numeric(Revenus_IRIS_PC$DEC_PAUT21)

## Warning: NAs introduced by coercion

# On a encore un souci car la décimale est codée avec une virgule.
# On va manipuler la chaîne de caractère pour mettre un point à la place

# On utilise la fonction `gsub` qui permet de remplacer des expressions régulières dans une chaîne de caractère
Revenus_IRIS_PC$RevPatrimoine <- gsub(",", ".", Revenus_IRIS_PC$DEC_PAUT21)

# On peut maintenant converir la variable dans un format numérique
Revenus_IRIS_PC$RevPatrimoine <- as.numeric(Revenus_IRIS_PC$RevPatrimoine)

# Et supprimer une variable devenue inutile
Revenus_IRIS_PC$DEC_PAUT21 <- NULL

summary(Revenus_IRIS_PC)

##      IRIS                DEP          DEC_MED21        DEC_GI21    
##  Length:2529        Min.   :75.00   Min.   : 7330   Min.   :202.0  
##  Class :character   1st Qu.:75.00   1st Qu.:18910   1st Qu.:360.0  
##  Mode  :character   Median :92.00   Median :27770   Median :390.0  
##                     Mean   :86.82   Mean   :29052   Mean   :402.4  
##                     3rd Qu.:93.00   3rd Qu.:37260   3rd Qu.:425.0  
##                     Max.   :94.00   Max.   :77180   Max.   :789.0  
##  RevPatrimoine   
##  Min.   : 0.100  
##  1st Qu.: 2.600  
##  Median : 4.600  
##  Mean   : 6.662  
##  3rd Qu.: 7.800  
##  Max.   :62.700

Variables numériques (variables quantitatives)

Une variable numérique peut faire l’objet d’un ensemble de résumés statistiques à l’aide de fonctions élémentaires

• min() : minimum
• max() : maximum
• mean() : moyenne
• sd() : écart-type (en anglais : standard deviation, soit sd en abrégé)
• sum() : somme

Pour l’exercice on peut créer un vecteur X à partir d’une variable du dataframe Revenus_IRIS_PC.

X <- Revenus_IRIS_PC$DEC_MED21
min(X)

## [1] 7330

max(X)

## [1] 77180

mean(X)

## [1] 29052

sd(X)

## [1] 11998.37

Pour calculer les quantiles on peut utiliser la fonction quantile() en paramétrant la valeur de fréquence cumulée ascendante

• quantile(X,0) : minimum
• quantile(X,0.10) : D1 (premier décile)
• quantile(X,0.25) : Q1 (premier quartile)
• quantile(X,0.5) : Q2 (médiane)
• quantile(X,0.75) : Q3 (troisième quartile)
• quantile(X,0.90) : D9 (dernier décile)
• quantile(X,1) : maximum

X <- Revenus_IRIS_PC$DEC_MED21
quantile(X,0.5)

##   50% 
## 27770

sel<-c(0,0.25,0.5,0.75,1)
quantile(X,sel)

##    0%   25%   50%   75%  100% 
##  7330 18910 27770 37260 77180

sel<-c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)

quantile(X,sel)

##    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
##  7330 14468 17426 20540 24222 27770 31996 35262 39180 44394 77180

Il peut arriver qu’une fonction soit manquante dans R, comme par exemple le coefficient de variation. Dans ce cas, on peut faire le calcul par des lignes de code ou créer sa propre fonction avec l’instruction function(). La fonction qui est stockée en mémoire apparaît dans la fenêtre Environnement. Lorsqu’on a créé plusieurs fonctions, on peut en faire un programme R qu’on charge en mémoire au début de chaque session. A plus long terme, on peut en faire un package qu’on partagera avec les autres utilisateurs de R (ceux que nous allons bientôt utiliser)

A titre d’exemple, nous créons une fonction cv() qui calcule le rapport entre l’écart-type et la moyenne d’une distribution :

# Calcul du CV avec des lignes de code
X <- Revenus_IRIS_PC$DEC_MED21
sd(X)/mean(X)

## [1] 0.4129963

cv<-function(var) {sd(var)/mean(var)}
cv(X)

## [1] 0.4129963

library(skimr)
skim(Revenus_IRIS_PC)

X<-cut(Revenus_IRIS_PC$DEC_MED21, c(0,10000,20000,30000,40000,50000,60000,70000, 80000,90000, 100000))
str(X)

##  Factor w/ 10 levels "(0,1e+04]","(1e+04,2e+04]",..: 4 5 4 4 5 6 5 5 5 4 ...

table(X)

## X
##     (0,1e+04] (1e+04,2e+04] (2e+04,3e+04] (3e+04,4e+04] (4e+04,5e+04] 
##             8           710           679           674           325 
## (5e+04,6e+04] (6e+04,7e+04] (7e+04,8e+04] (8e+04,9e+04] (9e+04,1e+05] 
##            98            29             6             0             0

plot(X,  main = "Fréquence simple des IRIS par catégorie de revenus", xlab = "classes")

X<-cut(Revenus_IRIS_PC$RevPatrimoine, c(0,20,40,60,80,100))
str(X)

##  Factor w/ 5 levels "(0,20]","(20,40]",..: 1 2 1 1 2 2 3 1 1 1 ...

table(X)

## X
##   (0,20]  (20,40]  (40,60]  (60,80] (80,100] 
##     2389      122       17        1        0

plot(X,  main = "Fréquence simple des IRIS par % de patrimoine dans le revenu des ménages", xlab = "classes")

X<-Revenus_IRIS_PC$DEC_GI21
hist(X,  main = "Histogramme de l'indice de gini (2021)", xlab = "classes")

Variables qualitatives

Les variables qualitatives nominales ou factor sont des objets composés d’une liste de numéros et d’une liste d’étiquettes.

On peut transformer une variable quantitative en facteur avec la fonction cut()

A partir des revenus, nous allons créer 4 catégories de communes en fonction de la médiane, auxquels nous pouvons simplement donner les noms des quartiles :

• Q1 : Iris les plus pauvres “Quartiers plus pauvres”
• Q2: Inférieure à la médiane “Quartiers modestes”
• Q3 : Supérieurs à la médiane “Quartiers aisés”
• Q4 ; 4e quartile, les plus riches médiane “Quartiers riches”

On peut écrire cette transformation des données de la manière suivante, en insérant manuellement les bornes, ce qui créé 4 intervalles :

Revenus_IRIS_PC$RevQuartiles <-cut(Revenus_IRIS_PC$DEC_MED21, breaks=c(7330,    18910.0,    27770.0,    37260.0,    77180.0))

On peut ensuite recoder les classes avec levels()

levels(Revenus_IRIS_PC$RevQuartiles)<-c("Quartiers plus pauvres","Quartiers modestes","Quartiers aisés","Quartiers riches")

Une autre possibilité consiste à utiliser une fonction ce qui créé 4 intervalles statistiques (fonction quantile) :

Revenus_IRIS_PC$RevQuartiles <- cut(Revenus_IRIS_PC$DEC_MED21, breaks=c(quantile(Revenus_IRIS_PC$DEC_MED21, probs = seq(0, 1, by=0.25))))

levels(Revenus_IRIS_PC$RevQuartiles)<-c("Quartiers plus pauvres","Quartiers modestes","Quartiers aisés","Quartiers riches")

Option que l’on peut tester par exemple avec le patrimoine

Revenus_IRIS_PC$CatRevPatrimoine <- cut(Revenus_IRIS_PC$RevPatrimoine, breaks=c(quantile(Revenus_IRIS_PC$RevPatrimoine, probs = seq(0, 1, by=0.25))))

levels(Revenus_IRIS_PC$CatRevPatrimoine)<-c("Faible","Inf Mediane","Sup Mediane","Elevés")

tab <- table(Revenus_IRIS_PC$RevQuartiles)
tab

## 
## Quartiers plus pauvres     Quartiers modestes        Quartiers aisés 
##                    633                    632                    631 
##       Quartiers riches 
##                    632

tab <- table(Revenus_IRIS_PC$RevQuartiles, Revenus_IRIS_PC$CatRevPatrimoine)
tab

##                         
##                          Faible Inf Mediane Sup Mediane Elevés
##   Quartiers plus pauvres    494         111          23      2
##   Quartiers modestes        144         328         137     23
##   Quartiers aisés            13         172         310    136
##   Quartiers riches            0          18         157    457

tab <- addmargins(tab)
tab

##                         
##                          Faible Inf Mediane Sup Mediane Elevés  Sum
##   Quartiers plus pauvres    494         111          23      2  630
##   Quartiers modestes        144         328         137     23  632
##   Quartiers aisés            13         172         310    136  631
##   Quartiers riches            0          18         157    457  632
##   Sum                       651         629         627    618 2525

kable(tab, format = "latex", digits = 2)

x <- Revenus_IRIS_PC$RevQuartiles
y <- Revenus_IRIS_PC$CatRevPatrimoine

tab <- table(x, y)

# Version brute
plot(tab) 

# Version améliorée
plot(tab, main = "Distribution des revenus et du patrimoine par IRIS en 2022", xlab ="Revenus médians (quartiles)", ylab ="% de revenus du patrimoine", las = 1, cex = .5)

test<-chisq.test(x, y)
test

## 
##  Pearson's Chi-squared test
## 
## data:  x and y
## X-squared = 2424.6, df = 9, p-value < 2.2e-16

vcd::mosaic(~ RevQuartiles + CatRevPatrimoine, data = Revenus_IRIS_PC, shade = TRUE, direction = "v", rot_labels=c(45,90,0,0))

library(vcd)

## Warning: package 'vcd' was built under R version 4.3.3

## Loading required package: grid

mosaic(~ RevQuartiles + CatRevPatrimoine, data = Revenus_IRIS_PC, shade = TRUE, 
            labeling = labeling_border(rot_labels = c(90,0, 0, 0), 
                                  offset_label =c(2,5,0, 0),
                                  varnames = c(FALSE, TRUE),
                                  just_labels=c("center","right"),
                                  tl_varnames = FALSE)  ) 

Comparaison d’une variable quantitative et d’une variable qualitative

La comparaison d’une variable quantitative et d’une variable qualitative peut prendre la forme d’une comparaison de la distribution statistique (boxplot) et d’un test de Fischer. On applique cette analyse à l’indice de Gini en fonction des revenus du patrimoine classés par quartiles

x <- Revenus_IRIS_PC$CatRevPatrimoine # Variable qualitative 
y <- Revenus_IRIS_PC$DEC_GI21 # Variable quantitative (indice de Gini)

plot(x,y, 
     col=c("blue","lightblue","pink", "red"),
     xlab ="Indice de Gini",
     ylab = "% des revenus du patrimoine",
     horizontal=T)

# Text de Fischer
mod<-aov(y~x)
summary(mod)

##               Df  Sum Sq Mean Sq F value Pr(>F)    
## x              3 4701027 1567009     497 <2e-16 ***
## Residuals   2521 7949269    3153                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 4 observations deleted due to missingness

x <- as.factor(Revenus_IRIS_PC$DEP) # Variable qualitative 
y <- Revenus_IRIS_PC$DEC_GI21 # Variable quantitative (indice de Gini)

plot(x,y, 
     col=c("blue","lightblue","pink", "red"),
     xlab ="Indice de Gini",
     ylab = "Départements",
     horizontal=T)

# Text de Fischer
mod<-aov(y~x)
summary(mod)

##               Df  Sum Sq Mean Sq F value Pr(>F)    
## x              3 3566909 1188970   330.3 <2e-16 ***
## Residuals   2525 9089080    3600                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparaison de deux variables quantitatives

La comparaison de deux variables qualitatives se fait généralement en traçant un diagramme bivarié, en calculant le coefficient de corrélation (Test de Pearson), et en analysant éventuellement les résultats d’une régression linéaire. On peut faire ces analyses pour deux couples de variables, et analyser la corrélation de l’une avec l’autre.

La régression linéraire permet d’analyser l’hypothèse d’une dépendance entre une variable à expliquer et une variable explicative. La fonction `lm()ù ou lm est l’abbréviation de linear model permet d”effectuer la plupart des modèles de régression linéaire basés sur la méthode des moindres carrés ordinaire. Sa syntaxe est a priori très simple et renvoie les coefficients b et a du modèle de régression.

x <- Revenus_IRIS_PC$DEC_MED21 # Variable quantitative
y <- Revenus_IRIS_PC$RevPatrimoine # Variable quantitative

plot(x,y, xlab="revenus médians",ylab="Revenus du patrimoine")

cor.test(y,x)

## 
##  Pearson's product-moment correlation
## 
## data:  y and x
## t = 56.425, df = 2527, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7288993 0.7634230
## sample estimates:
##       cor 
## 0.7466635

# Modèle
monmodel1<-lm(y~x)
monmodel1

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##  -6.2732415    0.0004453

# Plot avec la régression
plot(x,y, xlab="revenus médians",ylab="Revenus du patrimoine")
abline(monmodel1, col="red",lwd=2)

summary(monmodel1)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.704 -2.563 -0.297  1.398 43.558 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.273e+00  2.480e-01  -25.29   <2e-16 ***
## x            4.453e-04  7.891e-06   56.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.76 on 2527 degrees of freedom
## Multiple R-squared:  0.5575, Adjusted R-squared:  0.5573 
## F-statistic:  3184 on 1 and 2527 DF,  p-value: < 2.2e-16

residus <- monmodel1$residuals
Revenus_IRIS_PC$residus <- residus

x <- Revenus_IRIS_PC$DEC_MED21 # Variable quantitative (revenus médians)
y <- Revenus_IRIS_PC$DEC_GI21 # Variable quantitative (indice de Gini)

plot(x,y, xlab="revenus médians",ylab="indice de Gini")

cor.test(y,x)

## 
##  Pearson's product-moment correlation
## 
## data:  y and x
## t = 25.499, df = 2527, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4208254 0.4828459
## sample estimates:
##       cor 
## 0.4523825

# Modèle
monmodel2<-lm(y~x)
monmodel2

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##   3.249e+02    2.668e-03

# Plot avec la régression
plot(x,y, xlab="revenus médians",ylab="Revenus du patrimoine")
abline(monmodel2, col="blue",lwd=2)

summary(monmodel2)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -203.895  -41.024   -2.148   31.987  311.796 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.249e+02  3.288e+00   98.81   <2e-16 ***
## x           2.668e-03  1.046e-04   25.50   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 63.11 on 2527 degrees of freedom
## Multiple R-squared:  0.2046, Adjusted R-squared:  0.2043 
## F-statistic: 650.2 on 1 and 2527 DF,  p-value: < 2.2e-16

Importer des fichiers cartographiques (ArcGis .SHP)

# Installer le package sf et mapsf

# Ouvrir les fichiers .SHP
library(sf)
COM_PC = st_read("../public/geom/COM_PC.shp")

## Reading layer `COM_PC' from data source 
##   `/Users/renaud/SIG/cours_m2_traitementsstatscarto/public/geom/COM_PC.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 143 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 637302.2 ymin: 6843215 xmax: 671750.3 ymax: 6879246
## Projected CRS: RGF93 v1 / Lambert-93

IRIS_PC = st_read("../public/geom/IRIS_PC.shp")

## Reading layer `IRIS_PC' from data source 
##   `/Users/renaud/SIG/cours_m2_traitementsstatscarto/public/geom/IRIS_PC.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2752 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 637302.2 ymin: 6843215 xmax: 671750.3 ymax: 6879246
## Projected CRS: RGF93 v1 / Lambert-93

##### Visualiser cartes communes et PC #####
###### Méthode avec plot
plot(IRIS_PC)

plot(COM_PC)

###### Méthode avec Mapview. (Utilise OpenStreetMap, nécessite une connexion internet)
library(mapview)

mapview() + mapview (IRIS_PC) + mapview(COM_PC)

Joindre les données

Nous devons joindre les fichiers pour que les variables des revenus par IRIS puissent être jointes au fichier cartographique. Cette opération est importante, car nous avons supprimé des données manquantes lors du nettoyage des fichiers.

Nous pouvons joindre un data.frame à un objet sf en utilisant la fonction merge() et en utiisant les identifiants communs aux deux objets. Attention à l’ordre des arguments, l’objet retourné sera du même type que x. Il n’est pas possible de faire une jointure attributaire en utilisant deux objets sf.

# jointure attributaire 
IRIS_PC <- merge(
  x = IRIS_PC,             # l'objet sf
  y = Revenus_IRIS_PC,     # le data.frame contenant l'information statistique
  by.x = "CODE_IRIS",  # identifiant dans x
  by.y = "IRIS",  # identifiant dans y
  all.x = TRUE         # conserver toutes les lignes du fichier x (ici l'objet sf, ou fond de carte)
)

# Les deux objets ont bien été joints
head(IRIS_PC, 3)

## Simple feature collection with 3 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 651565.9 ymin: 6861790 xmax: 652179 ymax: 6862522
## Projected CRS: RGF93 v1 / Lambert-93
##   CODE_IRIS INSEE_COM                  NOM_COM IRIS                    NOM_IRIS
## 1 751010101     75101 Paris 1er Arrondissement 0101 Saint-Germain l'Auxerrois 1
## 2 751010102     75101 Paris 1er Arrondissement 0102 Saint-Germain l'Auxerrois 2
## 3 751010103     75101 Paris 1er Arrondissement 0103 Saint-Germain l'Auxerrois 3
##   TYP_IRIS DEP.x DEP.y DEC_MED21 DEC_GI21 RevPatrimoine RevQuartiles
## 1        H    75    NA        NA       NA            NA         <NA>
## 2        A    75    NA        NA       NA            NA         <NA>
## 3        A    75    NA        NA       NA            NA         <NA>
##   CatRevPatrimoine residus                       geometry
## 1             <NA>      NA MULTIPOLYGON (((652096.5 68...
## 2             <NA>      NA MULTIPOLYGON (((651932.7 68...
## 3             <NA>      NA MULTIPOLYGON (((651869 6862...

Cartes choroplèthes

Le code de base pour faire une carte est très simple, avec le package mapsf

mf_map(IRIS_PC, var = "DEC_MED21", type = "choro")

  mf_init(x = COM_PC) # On utilise un fond de carte (ici la carte communale) pour baliser le domaine de travail
# plot population
  mf_map(
  x = IRIS_PC,
  var = "DEC_MED21",
  type = "choro",
  breaks = "quantile",
  nbreaks = 5,
  pal = "Reds",
  border = "white",
  col_na = "lightgrey",
  lwd = 0.01,
  leg_pos = "left",
  leg_title = "Revenus médians 2021", leg_title_cex = 0.5, leg_val_cex = 0.4, add = TRUE
)
      mf_layout(title = "Les revenus des ménages à Paris et petite couronne en 2021",
          credits = "Sources: INSEE, 2021 ; IGN 2022. Réalisation R. Le Goix 2024",   scale = TRUE,
  arrow = TRUE, frame = TRUE)
# Et on rajoute la trame communale pour la lisibilité de la carte
 mf_map(COM_PC, col=NA, add=TRUE)

  mf_init(x = COM_PC) # On utilise un fond de carte (ici la carte communale) pour baliser le domaine de travail
# plot population
  mf_map(
  x = IRIS_PC,
  var = "RevPatrimoine",
  type = "choro",
  breaks = "quantile",
  nbreaks = 5,
  pal = "Blues",
  border = "white",
  col_na = "lightgrey",
  lwd = 0.01,
  leg_pos = "left",
  leg_title = "Revenus du patrimoine 2021 (%)", leg_title_cex   = 0.5, leg_val_cex = 0.4, add = TRUE
)
      mf_layout(title = "Les revenus du patrimoine en % du revenu des ménages en 2021 ",
          credits = "Sources: INSEE, 2021 ; IGN 2022. Réalisation R. Le Goix 2024",   scale = TRUE,
  arrow = TRUE, frame = TRUE)
# Et on rajoute la trame communale pour la lisibilité de la carte
 mf_map(COM_PC, col=NA, add=TRUE)

  mf_init(x = COM_PC) # On utilise un fond de carte (ici la carte communale) pour baliser le domaine de travail
# plot population
  mf_map(
  x = IRIS_PC,
  var = "residus",
  type = "choro",
  breaks = "quantile",
  nbreaks = 5,
  border = "white",
  col_na = "lightgrey",
  lwd = 0.01,
  leg_pos = "left",
  leg_title = "Residus regress", leg_title_cex  = 0.5, leg_val_cex = 0.4, add = TRUE
)
      mf_layout(title = "",
          credits = "Sources: INSEE, 2021 ; IGN 2022. Réalisation R. Le Goix 2024",   scale = TRUE,
  arrow = TRUE, frame = TRUE)
# Et on rajoute la trame communale pour la lisibilité de la carte
 mf_map(COM_PC, col=NA, add=TRUE)

Sauvegarde du nouveau fichier de données géographiques

Afin de pouvoir travailler sur ce fichier finalisé, on recommande d’enregistrer le fichier d’informations géographique contenant les données sélectionnées et transformées, dans un format pratique à exploiter : un geopackage (.gpkg).²

Les lignes de code suivantes exportent les données à l’IRIS dans la couche IRIS_PC du geopackage Paris_PC

st_write(obj = IRIS_PC, dsn = "data/Paris_PC.gpkg", layer = "IRIS_PC")

On peut ajouter d’autres couches dans le même géopackage, par exemple le fond de carte d’habillage à la commune COM_PC

st_write(obj = COM_PC, dsn = "data/Paris_PC.gpkg", layer = "COM_PC")

st_layers("data/Paris_PC.gpkg")

## Driver: GPKG 
## Available layers:
##   layer_name geometry_type features fields              crs_name
## 1    IRIS_PC Multi Polygon     2752     11 RGF93 v1 / Lambert-93
## 2     COM_PC       Polygon      143      1 RGF93 v1 / Lambert-93

RMarkdown : un “pack office killer”

Ce document est rédigé en RMarkdown, qui formate directement une page web.
Le principe du RMarkdown est de réaliser des documents (articles, rapport, pages web simples) à partir d’un document texte simple à formater. Les principaux formats de sortie sont des pages HTML, un document pdf (formaté comme un article), un document word ou ppt, des présentations dynamiques (shiny)³.

Menu de sélection des documents RMarkdown et principaux formats de sortie

Qu’est-ce que le rmarkdown ?

• “Des documents avec l’extension .Rmd : Développez votre code et vos idées côte à côte dans un seul document. Exécutez le code sous forme de morceaux individuels ou dans un document entier.

• Des Documents dynamiques : Rassemblez des graphiques, des tableaux et des résultats avec un texte narratif. Rendu dans une variété de formats tels que HTML, PDF, MS Word ou MS PowerPoint.

• Produire une recherche reproductible : Téléchargez votre rapport, créez un lien vers celui-ci ou joignez-le pour le partager. N’importe qui peut lire ou exécuter votre code pour reproduire votre travail”. Source : https://rstudio.github.io/cheatsheets/html/rmarkdown.html

On insère le code dans le texte, sous forme de morceaux appelés des “chunks”, et le texte est balises insérées dans le texte, pour faire des titres, de caractères gras, des passages en italiques, ou du texte barré. Le principe est assez similaire à celui de LaTex (qui est d’ailleurs le moteur utilisé pour l’exportation en PDF)

Flux de travail

1. Ouvrir un nouveau fichier .Rmd dans l’IDE RStudio en allant à Fichier > Nouveau fichier > R Markdown.
2. Incorporer le code par morceaux. Exécuter le code par ligne, par morceau ou en une seule fois.

Exercice 1 : intégrer dans un nouveau fichier .Rmd le code dejà préparé lors de la séance. Conseil: séparer dans 2 chunks différents la préparation des fichiers, et la réalisation des cartes. Executer le code pour s’assurer du bon fonctionnement de l’ensemble. Voir Rapport Univarié

3. Rédiger du texte et ajouter des tableaux, des figures, des images et des citations. Formater avec la syntaxe Markdown ou l’éditeur visuel Markdown de RStudio.

Exercice 2 : intégrer dans ce Rmarkdown le code pour la réalisation de traitements univariés (tableau univarié, diagramme en bâtons, etc…) de la section précédente. Nb : ces analyses n’ont pas encore été discutées lors de la séance du 13/9/2024, à partir de

4. Définir le(s) format(s) de sortie et les options dans l’en-tête YAML. Personnaliser les thèmes ou ajouter des paramètres pour exécuter ou ajouter de l’interactivité avec Shiny.

Exercice 3 : tester les sorties en HTML ou en PDF

5. Sauvegarder et rendre l’ensemble du document. Faites régulièrement un knit périodiquement pour prévisualiser votre travail au fur et à mesure que vous écrivez.

Un exemple de base de la structure d’un fichier RMarkdown

Le Tidyverse: extension, installation,

tidyverse . Cette section est reprise du tuto de Julien Barnier : https://juba.github.io/tidyverse/06-tidyverse.html

“Le terme tidyverse est une contraction de tidy (qu’on pourrait traduire par “bien rangé”) et de universe. Il s’agit en fait d’une collection d’extensions conçues pour travailler ensemble et basées sur une philosophie commune. Un des objectifs de ces extensions est de fournir des fonctions avec une syntaxe cohérente, qui fonctionnent bien ensemble, et qui retournent des résultats prévisibles.Elles abordent un très grand nombre d’opérations courantes dans R (la liste n’est pas exhaustive) :

• visualisation
• manipulation des tableaux de données
• import/export de données
• manipulation de variables
• extraction de données du Web
• programmation

Un des objectifs de ces extensions est de fournir des fonctions avec une syntaxe cohérente, qui fonctionnent bien ensemble, et qui retournent des résultats prévisibles. Elles sont en grande partie issues du travail d’Hadley Wickham, qui travaille désormais pour RStudio.

tidyverse est également le nom d’une extension qu’on peut installer de manière classique, soit via le bouton Install de l’onglet Packages de RStudio, soit en utilisant la commande :

install.packages("tidyverse")

Cette commande va en fait installer plusieurs extensions qui constituent le “coeur” du tidyverse, à savoir :

• ggplot2 (visualisation)
• dplyr (manipulation des données)
• tidyr (remise en forme des données)
• purrr (programmation)
• readr (importation de données)
• tibble (tableaux de données)
• forcats (variables qualitatives)
• stringr (chaînes de caractères)
• lubridate (manipulation de dates)

De la même manière, charger l’extension avec :

library(tidyverse)

Cette syntaxe s’applique à des jeux de données dits tidy, dont les principes sont les suivants :

• chaque variable est une colonne
• chaque observation est une ligne
• chaque type d’observation est dans une table différente

Source : Barnier, https://juba.github.io/tidyverse/06-tidyverse.html)

Manipulations de bases

Voir Barnier, https://juba.github.io/tidyverse/10-dplyr.html)

Une étude de cas : l’importation et le nettoyage d’une base de données

On utilise les données DVF (Déclarations de Valeur Foncières) pour importer et nettoyer une base de données. C’est une étape importante dans la chaîne de traitement des données. Les traitements réalisés dans cette partie sont adaptés de Mericskay and Demoraes (n.d.). Plusieurs travaux décrivent avec précision la structure de ces fichiers, à la parcelle, et qui nécessitent un traitement pour analyser les prix. En effet, Ces données décrivent des parcelles, et non des transactions immobilières (Casanova Enault, Boulay, and Coulon 2019). Les références de prix sont néanmoins sous réserves de précaution d’emplois, comparables à d’autres sources (Boulay et al. 2020; Casanova et al. 2017; Le Goix et al. 2021).

Les étapes de la préparation des données DVF. Source :Mericskay and Demoraes (n.d.)

On utilise une version adaptée du programme écrit par ces auteurs et disponible sur Github pour réaliser les étapes suivantes de nettoyage des données.

DVF_Script1.Rmd

Exercice : analyser le code et la manière de mettre en oeuvre la chaîne de traitement. Produire le fichier finale DVF nettoyé aggrégeant les départements de Paris et de la petite couronne.

DVFOK <- read.csv("data/DVFOK.csv", encoding="UTF-8", stringsAsFactors=FALSE)

Géographie des prix immobiliers et de l’abordabilité

ggplot(data = DVFOK)
# Ou, équivalent
ggplot(df)

ggplot(DVFOK) +
    geom_histogram()

ggplot(DVFOK) +
    geom_histogram(aes(x = prix))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(DVFOK) +
    geom_point(aes(x = surface, y = prix))

ggplot(DVFOK, aes(x = type, y = prixm2)) +
  geom_boxplot() +                # Scatter plot points
  labs(x = "Type", y = "prix m2", title = "prix au m2")

DVFOK$date <- as.Date(DVFOK$date)

DVFOK <- DVFOK %>%
  mutate(month = format(date, "%m"))  # Extract month as "mm" (numeric format)

ggplot(DVFOK, aes(x = month, y = prixm2)) +
  geom_boxplot() +                # Scatter plot points
  labs(x = "Mois", y = "prix m2", title = "Variation des prix au m2") +
#  scale_x_date(date_labels = "%m") +  # Format date as yyyy-mm-dd
  theme_minimal() 

              # Use a clean theme

df_temp <- DVFOK %>%
            group_by(month) %>%
              summarize(prixm2 = mean(prixm2))

ggplot(df_temp, aes(x = month, y = prixm2)) +
  geom_point() +                # Scatter plot points
  labs(x = "Mois", y = "prix m2", title = "Variation des prix au m2") +
#  scale_x_date(date_labels = "%m") +  # Format date as yyyy-mm-dd
  theme_minimal()  

ggplot(DVFOK) +
    geom_bin2d(aes(x = surface, y = prix))

DVFgeo <- st_as_sf(DVFOK, coords = c("longitude","latitude"), crs = 4326)
plot(st_geometry(DVFgeo), pch = 3, cex = 0.01)

# st_crs(IRIS_PC)
# st_crs(DVFgeo) #Obtenir les informations sur la projection utilisée

DVFgeo <- st_transform(DVFgeo, crs = 2154 ) # On convertit la projection

# Joindre le code IRIS dans le fichier des prix DVF
joined_data <- st_join(DVFgeo, IRIS_PC)

# On ne va garder que le prixm2, donc on nettoie le fichier. 
joined_data <- joined_data %>%
    select("CODE_IRIS","prixm2","prix") %>%  # On garde les deux variables qui nous intéressent
    st_drop_geometry() # On élimine la variable géométrique

# Calcul du prix moyen au m² par IRIS avec RBase
# IRISDVF1 <- aggregate(prixm2 ~ CODE_IRIS, data = joined_data, FUN = mean, na.rm = TRUE)
# 
# IRISDVF2 <- aggregate(prix ~ CODE_IRIS, data = joined_data, FUN = mean, na.rm = TRUE)


# Calcul du prix moyen et prix moyen au m² par IRIS avec tidyverse / dplyr
IRISDVF <- joined_data %>% 
          group_by(CODE_IRIS) %>% 
          summarise(prixm2 = mean(prixm2), prix=mean(prix))

IRIS_PC<-left_join(x=IRIS_PC,y=IRISDVF, by="CODE_IRIS")

  mf_init(x = COM_PC) # On utilise un fond de carte (ici la carte communale) pour baliser le domaine de travail
# plot population
  mf_map(
  x = IRIS_PC,
  var = "prixm2",
  type = "choro",
  breaks = "quantile",
#  pal = "Reds",
  nbreaks = 5,
  border = "white",
  col_na = "lightgrey",
  lwd = 0.01,
  leg_pos = "left",
  leg_title = "2023 (prix / m2)", leg_title_cex = 0.5, leg_val_cex = 0.4, add = TRUE
)
      mf_layout(title = "Prix au m2 moyen",
          credits = "Sources: INSEE, 2021 ; IGN 2022. Réalisation R. Le Goix 2024",   scale = TRUE,
  arrow = TRUE, frame = TRUE)
# Et on rajoute la trame communale pour la lisibilité de la carte
 mf_map(COM_PC, col=NA, add=TRUE)

Calculer et cartographier un nouvel indicateur : carte de l’abordabilité

• A partir des valeurs calculées dans DVF de prix au m^2, on détermine un indice d’abordabilité, qui peut être calculé de deux manières différentes :

  • un ratio \(prix/revenus\)
  • un ratio plus intuitif à comprendre, le nombre de mois de revenus par m2 acheté

• En tidyverse, les deux opérations peuvent être combinées en une seule commande.

## Calcul du ratio prix / revenus brut

IRIS_PC <- IRIS_PC %>%
  mutate(PIR_brut = prix/DEC_MED21,           ##  ratio prix / revenus brut
         PIR_par_mois = prixm2/(DEC_MED21/12) )  ##  ratio prix / revenus par mois

• On peut travailler la visualisation en construisant directement l’ensemble des histogrammes qui nous intéressent. Cette action nécessite un changement de format, d’un format “large” à un format long, qui va être réalisé en utilisant la commande gather.

# Créer un vecteur avec les variables que je veux sélectionner
selected_vars <- c("prixm2","PIR_brut","PIR_par_mois") 

# Etude de la commande gather
# From https://stackoverflow.com/questions/1181060
stocks <- tibble(
  time = as.Date("2009-01-01") + 0:9,
  X = rnorm(10, 0, 1),
  Y = rnorm(10, 0, 2),
  Z = rnorm(10, 0, 4)
)

toto <- stocks %>% gather("stock", "price", -time)

toto <- IRIS_PC %>% 
  gather(Attributes, value, selected_vars)

## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(selected_vars)
## 
##   # Now:
##   data %>% select(all_of(selected_vars))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Application à la production de l'histogramme
IRIS_PC %>% 
  gather(Attributes, value, selected_vars) %>% 
  ggplot(aes(x=value)) +
  geom_histogram(fill = "lightblue2", color = "black") + 
  facet_wrap(~Attributes, scales = "free_x") +
  labs(x = "Value", y = "density")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 835 rows containing non-finite outside the scale range
## (`stat_bin()`).

  mf_init(x = COM_PC) # On utilise un fond de carte (ici la carte communale) pour baliser le domaine de travail
# plot population
  mf_map(
  x = IRIS_PC,
  var = "PIR_par_mois",
  type = "choro",
  breaks = "quantile",
  pal = "Blues",
  nbreaks = 5,
  border = "white",
  col_na = "lightgrey",
  lwd = 0.01,
  leg_pos = "left",
  leg_title = "nb de mois pour acheter 1 m2)", leg_title_cex    = 0.5, leg_val_cex = 0.4, add = TRUE
)
      mf_layout(title = "Indice d'abordabilité de l'achat d'un logement en 2023",
          credits = "Sources: INSEE, 2021 ; IGN 2022. Réalisation R. Le Goix 2024",   scale = TRUE,
  arrow = TRUE, frame = TRUE)
# Et on rajoute la trame communale pour la lisibilité de la carte
 mf_map(COM_PC, col=NA, add=TRUE)

Données sur la pratique de la course à pieds (Source : STRAVA)

Pour réaliser cet exercice, on utilie un nouveau jeu de données, décrivant des segments parcourus par des sportifs, extraits en utilisant une API (Application Programming Interface) qui permet de construire des requêtes et d’extraires les données crowsourcées de Strava. Strava est un réseau social d’enregistrement et de suivi d’activité sportive, qui permet de type de requêtes de manière limitée et contrôlées. Voir Annexe sur la définition des objets et l’accès aux sources Données des segments Strava.

Le chunk ci-dessous importe, et optionnellement permet de contrôler visuellement le fichier d’information géographique créé par requêtes à partir de Strava, sur Paris et la Petite Couronne. Par précaution, on vérifie la nature de la projection géographique du fichier, afin de correspondre aux autres fonds de carte utilisés.

Quel est le type de projection utilisé par le fichier Source “Strava”. Dans quelle projection le transformer pour pouvoir être utilisé avec les fonds communaux et IRIS qui sont en Lambert93. Comment obtenir les informations sur les codes à utiliser pour vérifier et modifier ces projections cartographique ?

sdf_result <- st_read("geom/sdf_result.geojson")

## Reading layer `sdf_result' from data source 
##   `/Users/renaud/SIG/cours_m2_traitementsstatscarto/public/geom/sdf_result.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 4718 features and 11 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 2.09952 ymin: 48.65597 xmax: 2.64408 ymax: 49.02369
## Geodetic CRS:  WGS 84

sdf_result <- st_make_valid(sdf_result)# The error is due to one vertex with no edage in Polyline

#Obtenir les informations sur la projection utilisée
# st_crs(IRIS_PC)
# st_crs(sdf_result) 


sdf_result <- st_transform(sdf_result, crs = 2154 ) # On convertit la projection

mapview(COM_PC)+mapview(sdf_result) # Carto rapide pour vérifier le fichier

Calculer une densité

Des travaux ont montré un lien entre l’environnement urbain, les prix immobiliers, et les inégalités sociales, dans les pratiques de mobilité et de sports de plein air. On fait l’hypothèse qu’à Paris et Petite Couronne, la pratique de la course à pieds est corrélée avec la division sociale de l’espace.

“Crowdsourced data now represent a valuable asset to document urban mobility. Despite limitations, smartphones crowdsourced apps became an increasingly valuable tool for research, and many are used in urban studies to analyze travel behaviors and planning, with topic ranging from urban amenities, services, to prices. Using GPS-equipped participants and diaries, a study in Shanghaï found walking and cycling behaviors differ according to land-use, commercial facilities being less walkable than compact neighborhoods, although pollution in compact neighborhoods is also found detrimental to active travel (Zhang et al. 2023). Conrow, Mooney, and Wentz (2020) used Strava data to study bicycle ridership, in relation with the economic value of housing market.” (Le Goix and Lucas 2024)

Notre objectif méthodologique est donc d’analyser la densité des pratiques de course à pieds on veut pour cela calculer une densité de segments de course à pieds par km2.

• Procéder à une jointure spatiale déterminer la longueur des segments, et préparer un fichier qui calcule la somme des longues des segments par IRIS

### Km of strava segments by IRIS ###
# Perform a spatial join between line segments and IRIS polygons

joined_data <- st_join(sdf_result, IRIS_PC)

# Calculate the length of each line segment within each census tract
joined_data$segment_length <- st_length(joined_data$geometry)

# Group by census tract and calculate the total length of segments each IRIS
IRIS_summary <- joined_data %>%
                  group_by(CODE_IRIS) %>%
                  summarize(total_segment_length = sum(segment_length))

# Get rid of geometry
st_geometry(IRIS_summary) <- NULL

• Dans un fichier temporaire, calculer la surface des IRIS.

## Surface of IRIS in km2
IRIS_PC$area <- st_area(IRIS_PC$geometry)
IRIS_PC$area <- units::set_units(IRIS_PC$area, km^2)


IRIS_PC<-left_join(x=IRIS_PC,y=IRIS_summary, by="CODE_IRIS")

#pop_data$total_gpx_length <- units::set_units(pop_data$total_gpx_length, km) #Set units to km2
IRIS_PC$total_segment_length <- units::set_units(IRIS_PC$total_segment_length, km) #Set units to km2

#pop_data <- pop_data%>%
#  mutate(segment_density = total_segment_length/area,
#         gpx_density = total_gpx_length/area)

IRIS_PC <- IRIS_PC%>%
  mutate(segment_density = total_segment_length/area) # without GPX

#hist(pop_data$segment_density)

• Dans le cadre d’une pratique d’entraînement, de fitness et de compétition relayée sur l’application de Strava, la carte produite correspond à une densité de segments parcourus par des coureurs, exprimée en nb de km / km^2.

  mf_init(x = COM_PC) #Permet de définir un cadre de la carte
# plot Densité
  mf_map(
  x = IRIS_PC,
  var = "segment_density",
  type = "choro",
  breaks = "quantile",
  nbreaks = 5,
  pal = "Reds",
  border = "white",
  lwd = 0.01,
  leg_pos = "topright",
  leg_title = "Densité de segments", leg_title_cex  = 0.5, leg_val_cex = 0.4, add = TRUE
)
      mf_layout(title = "Course à pieds: densité de segments / km2",
          credits = "Sources: IGN 2024 ; Strava API 2024 | R. Le Goix, 2024",   scale = TRUE,
  arrow = TRUE, frame = TRUE)
# Et on rajoute la trame communale pour la lisibilité de la carte
 mf_map(COM_PC, col=NA, add=TRUE)

## Réalisation d’une série d’histogrammes

# Créer un vecteur avec les variables que je veux sélectionner
selected_vars <- c("DEC_MED21","DEC_GI21","RevPatrimoine","segment_density") 

IRIS_PC %>% 
  gather(Attributes, value, selected_vars) %>% 
  ggplot(aes(x=value)) +
  geom_histogram(fill = "lightblue2", color = "black") + 
  facet_wrap(~Attributes, scales = "free_x") +
  labs(x = "Value", y = "binwidth")

## Warning: attributes are not identical across measure variables; they will be
## dropped

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 1669 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Matrice de corrélation globale

Le bout de code ci-dessous permet de réaliser une matrice de corrélation d’une sélection de variable, en utilisant les fonctions de `ggplot2``.

####### 4.4 HEAT MAP CORRELATION MATRIX  #######

# Créer un vecteur avec les variables que je veux sélectionner
selected_vars <- c("DEC_MED21","DEC_GI21","RevPatrimoine","prixm2","PIR_par_mois","segment_density") 

toto <- IRIS_PC %>%
  drop_na(all_of(selected_vars)) %>%  # Ne garder que les cas complets des variables
  select(all_of(selected_vars)) %>%   # Sélectionner les colonnes dont on a besoin
  st_drop_geometry()                  # Nettoyer en enlevant la géométrie


toto <- toto[complete.cases(toto),]
### FUN Calc pct ACQ
total_col = apply(toto, 1, sum)
pcts = lapply(toto, function(x) {
  x / total_col
})
### pcts as df
pcts = as.data.frame(pcts)

corr <- round(cor(pcts), 2)

df <- reshape2::melt(corr)

gg <- ggplot(df, aes(x=Var1, y=Var2, fill=value, label=value)) + 
  geom_tile() + 
  theme_bw() + 
  labs(title="Correlation Matrix of selected data") +
  theme(text=element_text(size=10), legend.position="bottom",axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_gradient2()

gg

Préparation des données revenus et sources des données

# Importer fichier CSV INSEE
Revenus_IRIS <- read_delim("../public/data/BASE_TD_FILO_IRIS_2021_DEC_CSV/BASE_TD_FILO_IRIS_2021_DEC.csv", delim = ";", escape_double = FALSE, trim_ws = TRUE)

# Créer une nouvelle variable dans un DF avec les 2 caractères de gauche du code IRIS
Revenus_IRIS$DEP <- substr(Revenus_IRIS$IRIS, 1, 2)

# Selectionner les départements de Paris et petite couronne
Revenus_IRIS_PC <- Revenus_IRIS[Revenus_IRIS$DEP %in% c("75", "92", "93", "94"), ]

Revenus_IRIS_Paris <- Revenus_IRIS[Revenus_IRIS$DEP == "75", ]
write.csv(Revenus_IRIS_Paris, "../public/data/Revenus_IRIS_Paris.csv")
write.csv(Revenus_IRIS_PC, "../public/data/Revenus_IRIS_PC.csv")

#On supprime les fichiers dont on n'a plus besoin
rm(Revenus_IRIS)

IRIS <- st_read("../CONTOURS-IRIS_2-1_SHP_LAMB93_FXX-2022/CONTOURS-IRIS.shp")

##### Sélectionner Paris et PC ####
# Créer une nouvelle variable dans un DF avec les 2 caractères de gauche du code IRIS
IRIS$DEP <- substr(IRIS$INSEE_COM, 1, 2)

# Selectionner les départements petite couronne
IRIS_PC <- IRIS[IRIS$DEP %in% c("75", "92", "93", "94"), ]

##### Créer carte communes ####

# Agréger les géométries IRIS pour reconstruire les géométries communales
COM_PC <- aggregate(IRIS_PC[, "geometry"], 
                           by = list(IRIS_PC$INSEE_COM), 
                           FUN = st_union)



##### Visualiser cartes communes et PC #####
###### Méthode avec plot
plot(IRIS_PC)
plot(COM_PC)

###### Méthode avec Mapview. (Utilise OpenStreetMap, nécessite une connexion internet)
#library(mapview)

#mapview() + mapview (IRIS_PC) + mapview(COM_PC)

### On vérifie la projection cartographique de ces fichiers
st_crs(IRIS_PC)

##### Exporter les fichiers
st_write(IRIS_PC, "../public/geom/IRIS_PC.shp")
st_write(COM_PC, "../public/geom/COM_PC.shp")

rm(IRIS)

Note sur la préparation d’une bibliographie

On profite du passage ci-dessus présentant les data Strava pour noter la manière dont on peut utiliser RMarkdown pour gérer les éléments d’une bibliographie, à partir d’un fichier .BIB formaté selon le format BibTex utilisé dans Latex. Il est très aisé d’utiliser Zotero pour préparer ce fichier bibliographique. Il convient d’ajouter le chemin vers le fichier bibliographique dans l’en-tête du document RMarkdown, comme ci-dessous :

---
title: "Stage d'initiation à la statistique et cartographie avec R"
author: "Renaud Le Goix, Université Paris Cité"
date: "`r Sys.Date()`"
output:
  rmdformats::readthedown:
    highlight: kate
    toc: 4
 bibliography: Bibliographie.bib
---