diff --git a/figures/bandeau.png b/figures/bandeau.png index cc689e3..65bb52a 100644 Binary files a/figures/bandeau.png and b/figures/bandeau.png differ diff --git a/gwr_rzine.R b/gwr_rzine.R index 9f1b7f8..20098b9 100644 --- a/gwr_rzine.R +++ b/gwr_rzine.R @@ -359,7 +359,7 @@ moran.plot(data_immo$prix_med_z,neighbours_epci_w, # 1- utiliser la fonction queen_weights du package rgeoda pour calculer une matrice de contiguïté de type queen queen_w <- queen_weights(data_immo) # 2- Sortir la variable à étudier dans un vecteur -prix_med_z = data_immo["prix_med_z"] +prix_med_z <- data_immo["prix_med_z"] lisa <- local_moran(queen_w, prix_med_z) # Pour visualiser les résultats des LISA il faut les stocker dans des objets # ou dans des bases de données pour les représenter diff --git a/gwr_rzine.Rmd b/gwr_rzine.Rmd index 2c3ab9b..4f73858 100644 --- a/gwr_rzine.Rmd +++ b/gwr_rzine.Rmd @@ -15,6 +15,7 @@ output: rzine::readrzine: highlight: kate number_sections: true + fig_caption: yes csl: Rzine_citation.csl bibliography: biblio.bib nocite: | @@ -298,16 +299,16 @@ mf_credits("Sources: Notaires de France 2018, IGN Admin Express") Et avoir un aperçu rapide des autres données issues de l'INSEE : -```{r} +```{r fig.cap = "Variables explicatives, discrétisation par quantiles"} # correspondance entre le nom des variables et un intitulé plus lisible noms_variables <- list('perc_log_vac' = '% logements vacants', - 'perc_maison' = '% maisons', - 'perc_tiny_log' = '% petits logements', - 'dens_pop' = 'densité de population', - 'med_niveau_vis' = 'médiane niveau de vie', - 'part_log_suroccup' = '% logements suroccupés', - 'part_agri_nb_emploi' = '% agriculteurs', - 'part_cadre_profintellec_nbemploi' = '% cadres et professions intellectuelles') + 'perc_maison' = '% maisons', + 'perc_tiny_log' = '% petits logements', + 'dens_pop' = 'densité de population', + 'med_niveau_vis' = 'médiane niveau de vie', + 'part_log_suroccup' = '% logements suroccupés', + 'part_agri_nb_emploi' = '% agriculteurs', + 'part_cadre_profintellec_nbemploi' = '% cadres et professions intellectuelles') # les cartes par(mfrow = c(3,3)) for (var in colnames(data_immo)[6:13]) { @@ -1440,12 +1441,24 @@ data_immo$logvac.coef <- mod.gwr_g$SDF$perc_log_vac data_immo$tinylog.coef <- mod.gwr_g$SDF$perc_tiny_log data_immo$suroccup.coef <- mod.gwr_g$SDF$part_log_suroccup data_immo$cadre.coef <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi + +# correspondance entre le nom des variables et un intitulé plus lisible +noms_coefs <- list('agri.coef' = '% agriculteurs', + 'perc_maison.coef' = '% maisons', + 'dens_pop.coef' = 'densité de population', + 'med_vie.coef' = 'médiane niveau de vie', + 'logvac.coef' = '% logements vacants', + 'tinylog.coef' = '% petits logements', + 'suroccup.coef' = '% logements suroccupés', + 'cadre.coef' = '% cadres et professions intellectuelles') ``` Pour réaliser les cartes, avec une discrétisation standardisée centrée sur 0 : -```{r fig.height = 10, fig.cap=paste("Sources : Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")} -par(mfrow = c(4, 2)) +```{r fig.height = 10, fig.cap="Coefficients de la GWR (Sources : Notaires de France 2018, INSEE 2019, IGN Admin Express 2021)"} +# par(mfrow = c(4, 2)) +par(mfrow = c(2, 4)) +par(mai = c(0, 0, 0, 0)) for (var in colnames(data_immo)[22:29]) { # calcul des limites de classe res <- discr(values = data.frame(data_immo)[, var], @@ -1472,12 +1485,14 @@ for (var in colnames(data_immo)[22:29]) { lwd = 0.1, breaks = breaks, pal = palette, - leg_pos = "left", + # leg_pos = "left", + leg_pos = NA, leg_title = NA, leg_box_border = "gray", leg_val_rnd = 0) - var_name <- substr(var, 1, nchar(var) - 5) - mf_title(noms_variables[var_name]) + # var_name <- substr(var, 1, nchar(var) - 5) + mf_title(noms_coefs[var]) + # mf_title(var) } ``` @@ -1499,7 +1514,7 @@ data_immo$cadre.t <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi_TV df <- as.data.frame(data_immo) # On passe les t-values en valeurs absolues pour voir la plus grande # contribution dans un sens sens ou dans l'autre -data_immo$contribmax<- colnames(df[, c(30:37)])[max.col(abs(df[, c(30:37)]),ties.method="first")] +data_immo$contribmax <- colnames(df[, c(30:37)])[max.col(abs(df[, c(30:37)]),ties.method="first")] ``` ```{r} @@ -1511,15 +1526,24 @@ mf_map(x = data_immo, pal = brewer.pal(6,'Set2'), border = "white", lwd = 0.2, - leg_box_border = "white", - leg_adj = c(0, 1.5)) + leg_pos = NA) mf_map(x = region_sf, add = TRUE, type = "base", col = NA, border = "#736F6E", lwd = 1.5) -mf_title('variables contribuant le plus "localement" à l’explication de la variabilité des prix médian par EPCI') +mf_legend(type = "typo", + val = c("Densité de population", + "% logements vacants", + "% maisons", + "médiane du niveau de vie", + "% logements suroccupés", + "% petits logements"), + title = "", + pal = brewer.pal(6,'Set2'), + box_border = "white") +mf_title("Variables contribuant le plus à l’explication de la variabilité du prix médian au niveau de l'EPCI") mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021") ``` @@ -1570,6 +1594,7 @@ mf_map(x = data_immo, border = "gray", lwd = 0.2, leg_box_border = "gray", + leg_title = "Nombre de variables", leg_adj = c(0, 1.5)) mf_map(x = region_sf, add = TRUE, @@ -1609,6 +1634,7 @@ mf_map(x = data_immo, pal= "Reds", border = "gray", lwd = 0.2, + leg_title = "Nombre de betas", leg_box_border = "gray", leg_adj = c(0, 1.5)) mf_map(x = region_sf, @@ -1617,7 +1643,7 @@ mf_map(x = region_sf, col = NA, border = "gray", lwd = 1.5) -mf_title("Nombre des Betas significatifs par EPCI (p-value)") +mf_title("Nombre des betas significatifs par EPCI (p-value)") mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021") ``` @@ -1652,14 +1678,15 @@ mf_map(x = data_immo, lwd = 0.1, leg_box_border = "grey3", pal = mypal2, - leg_title = "Classe P-value") + leg_title = "p-value", + leg_size = 0.8) mf_map(x = region_sf, add = TRUE, type = "base", col = NA, border = "#736F6E", lwd = 1.5) -mf_title("P-value du coefficient de la part d'emploi agriculteurs") +mf_title("P-value du coefficient du % agriculteurs") mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021") # Pour la densité de population @@ -1681,7 +1708,8 @@ mf_map(x = data_immo, lwd = 0.1, leg_box_border = "grey3", pal=mypal2, - leg_title = "Classe P-value") + leg_title = "p-value", + leg_size = 0.8) mf_map(x = region_sf, add = TRUE, type = "base", diff --git a/gwr_rzine.html b/gwr_rzine.html index 3d073bc..d729c9e 100644 --- a/gwr_rzine.html +++ b/gwr_rzine.html @@ -5128,7 +5128,7 @@
Le schéma ci-dessous représente les différents cheminements possibles -pour expliques les variations de Y. Les étapes suivies ici sont -représentées en bleu.
+Le schéma ci-dessous (Frédéric Audard, 2024) représente les +différents cheminements possibles pour expliques les variations de Y. +Les étapes suivies ici sont représentées en bleu.
Ce fichier est composé des 10 variables suivantes (les données datent de 2019 sauf prix médian 2018) :
Nous allons également charger une couche de régions qui nous servira d’habillage pour les cartes :
shp_path <- here("data", "REGION.shp")
@@ -5570,8 +5570,8 @@ 1.3 Jointure des données
## [1] 1242
# filtre des données de la jointure pour ne voir que les epci sans correspondance dans immo_df
datatable(data_immo[is.na(data_immo$prix_med),])
Cependant, la VD étant prix_med
les lignes vides ne nous
@@ -5616,13 +5616,13 @@
Et avoir un aperçu rapide des autres données issues de l’INSEE :
# correspondance entre le nom des variables et un intitulé plus lisible
noms_variables <- list('perc_log_vac' = '% logements vacants',
- 'perc_maison' = '% maisons',
- 'perc_tiny_log' = '% petits logements',
- 'dens_pop' = 'densité de population',
- 'med_niveau_vis' = 'médiane niveau de vie',
- 'part_log_suroccup' = '% logements suroccupés',
- 'part_agri_nb_emploi' = '% agriculteurs',
- 'part_cadre_profintellec_nbemploi' = '% cadres et professions intellectuelles')
+ 'perc_maison' = '% maisons',
+ 'perc_tiny_log' = '% petits logements',
+ 'dens_pop' = 'densité de population',
+ 'med_niveau_vis' = 'médiane niveau de vie',
+ 'part_log_suroccup' = '% logements suroccupés',
+ 'part_agri_nb_emploi' = '% agriculteurs',
+ 'part_cadre_profintellec_nbemploi' = '% cadres et professions intellectuelles')
# les cartes
par(mfrow = c(3,3))
for (var in colnames(data_immo)[6:13]) {
@@ -5637,7 +5637,12 @@ 1.3 Jointure des données
mf_title(noms_variables[var])
mf_credits('Sources: INSEE 2019, IGN Admin Express 2021')
}
Dans le cas où vous préféreriez manipuler vos données sous un format sp (package sp), ou dans le cas où ce format serait requis pour utiliser certains packages ou certaines formules, vous pouvez convertir votre @@ -5789,8 +5794,8 @@
Il est important de réaliser également pour les VI cet histogramme que nous venons de faire pour la VD :
- - + +gtsummary
:
-Pour obtenir les résidus :
- - + +On peut également les visualiser :
par(mfrow=c(1,3))
# diagramme quantile-quantile qui permet de vérifier l'ajustement
@@ -7672,7 +7677,7 @@ 3.4.3 Analyser les
hist(rstudent(mod.lm), breaks = 50, col="darkblue", border="white", main="Analyse visuelle des résidus")
# un graphique pour visualiser l'homoscédasticité des résidus
plot(rstudent(mod.lm))
Si la voie graphique ne vous inspire pas il existe des tests statistiques qui permettent de vérifier la normalité des résidus ou bien leur homoscédasticité.
@@ -7704,8 +7709,8 @@# Pour relancer un nouveau modèle sans l'individu le plus extrême
# Notez que l'on peut en supprimer plusieurs d'un coup avec subset=-c(36,266)
mod.lmx <- update(mod.lm, subset=-266)
@@ -7751,7 +7756,7 @@ 3.4.3 Analyser les
-
+
## [1] 36 180
# Il est possible de comparer les deux modèles et les coefficients
car::compareCoefs(mod.lm, mod.lmx, pvals = TRUE)
@@ -8185,7 +8190,7 @@ 4.2 Niveau local
d’autocorrélation au niveau local, étape nécessaire pour continuer vers
la GWR. Luc Anselin va développer le concept d’indicateur
d’autocorrélation spatiale locale avec les LISA (Local
-Indicators of Spatial Association).
+Indicators of Spatial Association) (Anselin 1995).
D’après Luc Anselin, le LISA de chaque individu statistique indique
l’intensité du regroupement spatial de valeurs similaires autour de
cette individu. Autrement dit, un individu avec un LISA élevé va avoir
@@ -8197,8 +8202,8 @@
4.2 Niveau local
du niveau global par une approche locale. On cherche à la fois à
détecter des clusters qui correspondent à un regroupement significatif
de valeurs identiques dans une zone particulière, et à repérer des zones
-de non stationnarité, c’est-à-dire qui ne suivent pas le processus
-global.
+de non stationnarité spatiale, c’est-à-dire qui ne suivent pas le
+processus global.
Le logiciel GeoDa, développé par Luc Anselin et son équipe pour
étudier l’autocorrélation spatiale et les LISA, constitue une bonne
@@ -8588,7 +8593,7 @@
5.4 Interprétation des
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
-## Program starts at: 2024-07-29 12:15:14.454915
+## Program starts at: 2024-07-29 16:26:46.146628
## Call:
## gwr.basic(formula = formula, data = data, bw = bw, kernel = kernel,
## adaptive = adaptive, p = p, theta = theta, longlat = longlat,
@@ -8676,7 +8681,7 @@ 5.4 Interprétation des
## Adjusted R-square value: 0.9084372
##
## ***********************************************************************
-## Program stops at: 2024-07-29 12:15:57.955583
+## Program stops at: 2024-07-29 16:27:27.148296
Cette visualisation des résultats nous propose d’abord un rappel
complet du modèle linéaire classique. Puis viennent ensuite les
informations concernant la GWR. Le premier indicateur à analyser est le
@@ -8693,8 +8698,8 @@
5.4 Interprétation des
quartiles. Cela permet de rendre compte de la stationnarité de l’effet
ou non. Dans notre cas on voit qu’il y a bien une variation et même dans
certains cas une inversion des signes. Cela laisse supposer une non
-stationnarité des effets : un effet local peut être présent qui ne
-suivrait pas l’effet global.
+stationnarité spatiale des effets : un effet local peut être présent qui
+ne suivrait pas l’effet global.
Par exemple, pour le pourcentage de logement vacant avec un
coefficient global (modèle linéaire) de \(-287\), quand ce pourcentage augmente le
prix médian baisse. En simplifiant le pourcentage baisse d’une unité le
@@ -8741,8 +8746,8 @@
5.4 Interprétation des
#Pour voir à quoi il ressemble dans ce document
datatable(mod.gwr_g$SDF@data)
-
-
+
+
## [1] "Intercept" "perc_log_vac"
@@ -8821,8 +8826,8 @@ 5.4.1 Étude des
5.4.2 Étude des
coefficients
-Pour visualiser la non stationnarité des effets des VI la solution la
-plus efficace est la carte.
+Pour visualiser la non stationnarité spatiale des effets des VI la
+solution la plus efficace est la carte.
# On ajoute à data_immo les coefficients
data_immo$agri.coef <- mod.gwr_g$SDF$part_agri_nb_emploi
data_immo$perc_maison.coef <- mod.gwr_g$SDF$perc_maison
@@ -8831,47 +8836,62 @@ 5.4.2 Étude des
data_immo$logvac.coef <- mod.gwr_g$SDF$perc_log_vac
data_immo$tinylog.coef <- mod.gwr_g$SDF$perc_tiny_log
data_immo$suroccup.coef <- mod.gwr_g$SDF$part_log_suroccup
-data_immo$cadre.coef <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi
+data_immo$cadre.coef <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi
+
+# correspondance entre le nom des variables et un intitulé plus lisible
+noms_coefs <- list('agri.coef' = '% agriculteurs',
+ 'perc_maison.coef' = '% maisons',
+ 'dens_pop.coef' = 'densité de population',
+ 'med_vie.coef' = 'médiane niveau de vie',
+ 'logvac.coef' = '% logements vacants',
+ 'tinylog.coef' = '% petits logements',
+ 'suroccup.coef' = '% logements suroccupés',
+ 'cadre.coef' = '% cadres et professions intellectuelles')
Pour réaliser les cartes, avec une discrétisation standardisée
centrée sur 0 :
-par(mfrow = c(4, 2))
-for (var in colnames(data_immo)[22:29]) {
- # calcul des limites de classe
- res <- discr(values = data.frame(data_immo)[, var],
- center = 0,
- pos_center = "class_center",
- interval = sd(data.frame(data_immo)[, var])*0.5,
- min_nb = 10)
- breaks <- res[[1]]
- # palette de couleurs
- nb_cl_sup0 <- res[[2]]
- nb_cl_inf0 <- res[[3]]
- if (nb_cl_inf0 > 0) {
- palette = mf_get_pal(n = c(nb_cl_inf0, nb_cl_sup0),
- pal = c("Teal", "Peach"),
- neutral = "#f5f5f5")
- } else { # cas de la médiane du niveau de vie où la valeur min est supérieure à 0
- palette = mf_get_pal(n = c(nb_cl_sup0), pal = c("Peach"), rev = TRUE)
- }
- # la carte
- mf_map(x = data_immo,
- var = var,
- type = "choro",
- border = "gray",
- lwd = 0.1,
- breaks = breaks,
- pal = palette,
- leg_pos = "left",
- leg_title = NA,
- leg_box_border = "gray",
- leg_val_rnd = 0)
- var_name <- substr(var, 1, nchar(var) - 5)
- mf_title(noms_variables[var_name])
-}
+# par(mfrow = c(4, 2))
+par(mfrow = c(2, 4))
+par(mai = c(0, 0, 0, 0))
+for (var in colnames(data_immo)[22:29]) {
+ # calcul des limites de classe
+ res <- discr(values = data.frame(data_immo)[, var],
+ center = 0,
+ pos_center = "class_center",
+ interval = sd(data.frame(data_immo)[, var])*0.5,
+ min_nb = 10)
+ breaks <- res[[1]]
+ # palette de couleurs
+ nb_cl_sup0 <- res[[2]]
+ nb_cl_inf0 <- res[[3]]
+ if (nb_cl_inf0 > 0) {
+ palette = mf_get_pal(n = c(nb_cl_inf0, nb_cl_sup0),
+ pal = c("Teal", "Peach"),
+ neutral = "#f5f5f5")
+ } else { # cas de la médiane du niveau de vie où la valeur min est supérieure à 0
+ palette = mf_get_pal(n = c(nb_cl_sup0), pal = c("Peach"), rev = TRUE)
+ }
+ # la carte
+ mf_map(x = data_immo,
+ var = var,
+ type = "choro",
+ border = "gray",
+ lwd = 0.1,
+ breaks = breaks,
+ pal = palette,
+ # leg_pos = "left",
+ leg_pos = NA,
+ leg_title = NA,
+ leg_box_border = "gray",
+ leg_val_rnd = 0)
+ # var_name <- substr(var, 1, nchar(var) - 5)
+ mf_title(noms_coefs[var])
+ # mf_title(var)
+}
-
+
Les cartes des betas vont illustrer la variation des effets en
@@ -8900,7 +8920,7 @@
5.4.2 Étude des
df <- as.data.frame(data_immo)
# On passe les t-values en valeurs absolues pour voir la plus grande
# contribution dans un sens sens ou dans l'autre
-data_immo$contribmax<- colnames(df[, c(30:37)])[max.col(abs(df[, c(30:37)]),ties.method="first")]
par(mfrow = c(1, 1))
# Carte
mf_map(x = data_immo,
@@ -8909,17 +8929,26 @@ 5.4.2 Étude des
pal = brewer.pal(6,'Set2'),
border = "white",
lwd = 0.2,
- leg_box_border = "white",
- leg_adj = c(0, 1.5))
-mf_map(x = region_sf,
- add = TRUE,
- type = "base",
- col = NA,
- border = "#736F6E",
- lwd = 1.5)
-mf_title('variables contribuant le plus "localement" à l’explication de la variabilité des prix médian par EPCI')
-mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")
Au-delà de la non-stationnarité des VI, cette carte met en évidence un autre phénomène : Dans un modèle linéaire classique, la sélection des VI se fait de manière itérative (ascendante, descendante, ou mixte). @@ -8985,16 +9014,17 @@
Il se peut que cela soit plus intéressant d’utiliser les p-value, notamment si vous avez moins de 200 individus.
# Les p-value ne sont pas fournis dans le modèle de la GWR
@@ -9021,17 +9051,18 @@ 5.4.2 Étude des
pal= "Reds",
border = "gray",
lwd = 0.2,
- leg_box_border = "gray",
- leg_adj = c(0, 1.5))
-mf_map(x = region_sf,
- add = TRUE,
- type = "base",
- col = NA,
- border = "gray",
- lwd = 1.5)
-mf_title("Nombre des Betas significatifs par EPCI (p-value)")
-mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")
Sur ces deux dernières cartes (nombre de betas significatifs par EPCI avec les t-value et p-value) nous avons pris volontairement une liberté @@ -9071,45 +9102,47 @@
Ces cartes des p-value sont particulièrement importantes car elles nous donnent les endroits où l’effet est significatif. En effet, on sait que la VI a effet global qui est significatif, qu’elle a en plus une @@ -9173,10 +9206,15 @@