1.1 Chargement des données sur le prix de l’immobilier par EPCI

Télécharger les données

Pour rejouer l’analyse, décompressez les données contenues dans l’archive data.zip puis ouvrez le jeu de données donnees_standr.csv et la couche géographique EPCI.shp. La couche REGION.shp est utilisée à des fins d’habillage cartographique.

# On situe le dossier dans lequel se trouve nos données
csv_path <- here("data", "donnees_standr.csv")
# Lecture du CSV dans un dataframe
immo_df <- read.csv2(csv_path)
# Pour visualiser les 10 premières lignes
datatable(format(immo_df, digits = 1, scientific = FALSE))

Ce fichier est composé des 10 variables suivantes (les données datent de 2019 sauf prix médian 2018) :

• SIREN : code SIREN de l’EPCI
• prix_med : prix médian par EPCI à la vente (au m²)
• perc_log_vac : % logements vacants
• perc_maison : % maisons
• perc_tiny_log : % petits logements
• dens_pop : densité de population
• med_niveau_vis : médiane du niveau de vie
• part_log_suroccup : % logements suroccupés, cf. définition de l’INSEE
• part_agri_nb_emploi : % agriculteurs
• part_cadre_profintellec_nbemploi : % cadres et professions intellectuelles

La variable SIREN nous servira de “clé” pour joindre ces données statistiques aux données spatiales, la variable prix_med sera la variable que nous chercherons à expliquer (VD), et toutes les autres seront nos variables explicatives (VI).

1.2 Chargement des données géographiques : les EPCI de France métropolitaine

Ces données proviennent de la base ADMIN-EXPRESS-COG de l’IGN, édition 2021. Le format d’entrée est le shapefile que nous allons lire dans un objet sf :

# Lecture du shapefile en entrée dans un objet sf
shp_path <- here("data", "EPCI.shp")
epci_sf <- st_read(shp_path, quiet = TRUE)
# On vérifie dans les données qu'il existe un champ correspondant au code SIREN
# qui permettra de lier les données de immo_df à cet objet sf
datatable(epci_sf)
## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

Nous allons également charger une couche de régions qui nous servira d’habillage pour les cartes :

shp_path <- here("data", "REGION.shp")
region_sf <- st_read(shp_path, quiet = TRUE)

1.3 Jointure des données géographiques et tabulaires

On constate que les deux jeux de données n’ont pas exactement le même nombre de lignes. En effet, le jeu de données immo_df possède moins de lignes que l’objet sf epci_sf. Cela indique simplement que nous n’avons pas l’indication du prix médian de l’immobilier pour tous les EPCI de France métropolitaine.

nrow(immo_df)
## [1] 1223
nrow(epci_sf)
## [1] 1242

Il peut être intéressant d’identifier et visualiser les EPCI qui n’ont pas de correspondance dans le tableau de données immo_df. Pour ce faire on réalise la jointure on conservant toutes les lignes de epci_sf :

# L'option all.x = TRUE permet de garder toutes les lignes de epci_sf,
# même celles qui n'ont pas de correspondance dans immo_df
data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN", all.x = TRUE)
nrow(data_immo)
## [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),])

plot(data_immo[is.na(data_immo$prix_med),]$geometry)

Cependant, la VD étant prix_med les lignes vides ne nous intéressent pas, nous ne les conserverons pas car elles pourraient poser problème lors de la réalisation de nos analyses. On refait la jointure en ne gardant que les EPCI ayant une correspondance dans le tableau de données :

data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN")
nrow(data_immo)

## [1] 1223

On peut maintenant représenter les données du prix médian de l’immobilier par EPCI sous forme de carte :

mf_map(x = data_immo, 
       var = "prix_med", 
       type = "choro",
       breaks = "jenks",
       nbreaks = 9,
       pal = "Mint",
       lwd = 0.01,
       leg_title = "Discrétisation Jenks",
       leg_box_border = NA,
       leg_val_rnd = 0,
       leg_adj = c(0, 1.5))
mf_map(x = region_sf, 
       add = TRUE,
       type = "base",
       col = NA,
       border = "#34282C",
       lwd = 1)
mf_title("Prix médian de l'immobilier au m² par EPCI")
mf_credits("Sources: Notaires de France 2018, IGN Admin Express")

Et avoir un aperçu rapide des variables explicatives issues de l’INSEE. Ces cartes suivent une 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')
# Création des cartes
par(mfrow = c(3,3))
for (var in colnames(data_immo)[6:13]) {
  mf_map(x = data_immo,
         var = var,
         type = "choro",
         breaks = "quantile",
         nbreaks = 7,
         pal = "Purples",
         lwd = 0.01,
         leg_pos = NA)
  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 objet sf en sp à l’aide de la ligne de commande suivante (nous en aurons besoin pour la suite) :

data_immo_sp <- as(data_immo, "Spatial")

2.1 Voisinage

Nous ne développerons pas ici ce qu’implique la définition d’un voisinage. Pour cela, nous vous renvoyons vers les travaux de Sébastien Oliveau (Oliveau 2011).

Ici, il est important de savoir qu’un voisinage peut être de 3 types :

• Basé sur la contiguïté
• Basé sur la distance
• Basé sur la proximité

Lorsque nous travaillons avec des polygones (comme c’est le cas ici), le plus souvent l’usage consiste à se baser sur une matrice de contiguïté. Il faut par ailleurs savoir qu’il existe plusieurs types de voisinages basé sur la contiguïté. Dans un cas classique nous utiliserons celui de type queen. Queen est une référence à la reine aux échecs, qui peut se déplacer dans toutes les directions. Nous considérons ici les voisins contigus à notre polygone de tous côtés. Il s’oppose au type rook qui fait référence à la tour, les voisins seront donc définis à partir des mouvements de cette pièce sur l’échiquier (dans toutes les directions sauf en diagonale).

Figure 2.8 du manuel INSEE Codifier la structure de voisinage (De Bellefon et al. 2018)

R permet assez simplement de définir le voisinage.

# Création de la liste des voisins : avec l'option queen = TRUE, 
# sont considérés comme voisins 2 polygones possédant au moins 1 sommet commun
neighbours_epci <- poly2nb(data_immo, queen = TRUE)
# Obtention des coordonnées des centroïdes
coord <- st_coordinates(st_centroid(data_immo))

Voici la représentation graphique du voisinage :

# Faire un graphe de voisinage
par(mfrow = c(1,1))
par(mai = c(0, 0, 0, 0)) # Pour réduire les marges de la figure
plot.nb(x = neighbours_epci, coords = coord, points = FALSE)

Pour comprendre ce que contient neighbours_epci :

# Si on prend le 1er élément de neighbours_epci, 
# on voit qu'il a pour voisins les polygones 62, 74 etc.
neighbours_epci[[1]]

## [1]   62   74  338 1135 1136 1137 1140

# Ce qu'on peut vérifier sur la carte :
neighbours1 <- data_immo[c(1,62,74,338,1135,1136,1137,1140),]
neighbours1$index <- rownames(neighbours1)
mf_map(x = neighbours1, border = "white")
# Pour voir les numéros
mf_label(x = neighbours1, var = "index", halo = TRUE)
# Pour ajouter les liens entre voisins, il faut convertir neighbours_epci en objet sf
neighbours_epci_sf <- spdep::nb2lines(neighbours_epci, coords=coord, proj4string="EPSG:2154", as_sf=T)
mf_map(x = neighbours_epci_sf[neighbours_epci_sf$i == 1,], col = "grey30", add = TRUE)

Nous précisons qu’un voisinage peut aussi tout à fait se calculer lorsque l’on n’a pas de polygones mais simplement des coordonnées (des points). Les matrices de distances sont alors souvent plus adaptées. Pour définir le voisinage il faut utiliser les fonctions knearneigh() et knn2nb().

2.2 Création de la matrice de voisinage

Une fois le voisinage défini nous pouvons créer une matrice de voisinage, qui permettra d’attribuer un poids à chaque voisin.

# La fonction nb2listw attribue un poids à chaque voisin
# par exemple si un polygone a 4 voisins, ils auront chacun un poids de 1/4 = 0.25
#help("nb2listw")
neighbours_epci_w <- nb2listw(neighbours_epci)
# Les poids sont stockés dans le 3e élément de neighbours_epci_w
# par exemple si on veut connaître les poids des voisins du 1er élément :
neighbours_epci_w[[3]][1]

## [[1]]
## [1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

# Cet élément a 7 voisins qui ont donc un poids de 1/7 soit ~0,14

3.1 Exploration des variables

Cette première étape très importante permet d’étudier la distribution des données et d’identifier d’éventuels individus extrêmes qui pourrait venir perturber les résultats.

Ici par exemple nous serions en droit de nous poser la question de conserver dans notre jeu de données l’entité spatiale avec un prix médian à plus de 10 000 qui se détache très clairement de tous les autres EPCI.

Pour visualiser la distribution d’une variable quantitative l’histogramme est une bonne solution. Pour le réaliser nous utilisons dans ce cas le package plotly qui permet l’interactivité de la figure.

# Distribution de la variable dépendante :
add_histogram(plot_ly(data_immo, x = ~prix_med))

Il est important de réaliser également pour les VI cet histogramme que nous venons de faire pour la VD :

# Distribution des variables indépendantes :
a <- add_histogram(plot_ly(data_immo, x = ~(perc_log_vac), name = "perc_log_vac"))
b <- add_histogram(plot_ly(data_immo, x = ~(perc_maison), name = "perc_maison"))
c <- add_histogram(plot_ly(data_immo, x = ~(perc_tiny_log), name = "perc_tiny_log"))
d <- add_histogram(plot_ly(data_immo, x = ~(dens_pop), name = "dens_pop"))
e <- add_histogram(plot_ly(data_immo, x = ~(med_niveau_vis), name = "med_niveau_vis"))
f <- add_histogram(plot_ly(data_immo, x = ~(part_log_suroccup), name = "part_log_suroccup"))
g <- add_histogram(plot_ly(data_immo, x = ~(part_agri_nb_emploi), name = "part_agri_nb_emploi"))
h <- add_histogram(plot_ly(data_immo, x = ~(part_cadre_profintellec_nbemploi), name = "part_cadre_profintellec_nbemploi"))
fig = subplot(a, b, c, d, e, f, g, h, nrows = 2)
fig

3.2 Etude des corrélations

Très rapidement, la corrélation permet d’étudier le lien, la relation entre deux variables. La corrélation repose sur la covariance entre les variables.

Quand 2 variables covarient, pour un même individu, un écart à la moyenne pour une variable sera accompagné par un écart à la moyenne pour l’autre variable. Si les 2 écarts se font dans le même sens, on parle de corrélation positive, sinon il s’agit d’une corrélation négative.

La conception d’un modèle statistique doit absolument être le fruit d’une réflexion sur le choix des variables indépendantes (explicatives) et le choix de la méthode de régression. Pour définir un modèle de régression certaine règles doivent être respectées.

L’étude des corrélations peut apporter une aide précieuse dans cette réflexion. Elle pourra nous aider à choisir les variables à intégrer au modèle et dans le même temps à vérifier certaines des conditions de réalisation de notre régression.

Ainsi, une analyse des corrélations pourra vérifier :

• l’existence d’un lien entre les variables indépendantes et la variable à étudier. En effet dans une régression linéaire, il est nécessaire d’avoir une relation linéaire entre la VD et les différentes VI.
• la multicolinéarité des variables indépendantes. Les corrélations ne doivent pas être trop fortes entre les VI. Un coefficient > 0,7 en valeur absolue doit entraîner la suppression des variables concernées. Cela peut aussi être vérifié très efficacement avec le VIF (Variance Inflation Factor) mais peut se faire seulement après avoir lancé le modèle.
• que toutes les variables d’intérêt sont bien présentes dans le modèle (sauf dans le cas où cela induirait une trop grande multicolinéarité).

Pour calculer une matrice de corrélation :

# On commence par créer un dataframe identique à immo_df mais sans la colonne SIREN
data_cor <- immo_df %>% select(-SIREN)
# On peut ensuite créer la matrice de corrélation
immo_cor <- correlation(data_cor, redundant = TRUE)
# Et afficher cette matrice
summary(immo_cor)

## # Correlation Matrix (pearson-method)
## 
## Parameter                        | part_cadre_profintellec_nbemploi | part_agri_nb_emploi | part_log_suroccup | med_niveau_vis | dens_pop | perc_tiny_log | perc_maison | perc_log_vac | prix_med
## -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
## prix_med                         |                          0.55*** |            -0.45*** |           0.60*** |        0.59*** |  0.52*** |       0.49*** |    -0.60*** |     -0.68*** |  1.00***
## perc_log_vac                     |                         -0.31*** |             0.39*** |          -0.28*** |       -0.46*** | -0.21*** |      -0.19*** |     0.35*** |      1.00*** |         
## perc_maison                      |                         -0.58*** |             0.54*** |          -0.79*** |       -0.24*** | -0.46*** |      -0.75*** |     1.00*** |              |         
## perc_tiny_log                    |                          0.75*** |            -0.51*** |           0.87*** |        0.22*** |  0.62*** |       1.00*** |             |              |         
## dens_pop                         |                          0.59*** |            -0.29*** |           0.62*** |        0.18*** |  1.00*** |               |             |              |         
## med_niveau_vis                   |                          0.43*** |            -0.33*** |           0.15*** |        1.00*** |          |               |             |              |         
## part_log_suroccup                |                          0.65*** |            -0.43*** |           1.00*** |                |          |               |             |              |         
## part_agri_nb_emploi              |                         -0.57*** |             1.00*** |                   |                |          |               |             |              |         
## part_cadre_profintellec_nbemploi |                          1.00*** |                     |                   |                |          |               |             |              |         
## 
## p-value adjustment method: Holm (1979)

On peut aussi visualiser les corrélations de manière plus lisible à l’aide d’un corrélogramme, en utilisant le package corrplot :

mat_cor_comp <- summary(immo_cor, redundant = TRUE)

# Nom des lignes = valeurs de la première colonne ("Parameter")
rownames(mat_cor_comp ) <- mat_cor_comp[,1]

# Transformation du dataframe en objet matrice (+ suppression de la première colonne)
mat_cor <- as.matrix(mat_cor_comp[,-1])

# Calcul du nombre total d'individus
nb <- nrow(data_cor)

# Calcul des matrices de p-values et des intervalles de confiance
p.value <- cor_to_p(mat_cor, n = nb, method = "auto")

# Extraction de la matrice des p-values uniquement
p_mat <- p.value$p

# Affichage du corrélogramme
corrplot(mat_cor, 
         p.mat = p_mat, 
         type = "upper", 
         order = "hclust", 
         addCoef.col = "white", 
         tl.col = "gray",
         number.cex = 0.5,
         tl.cex= 0.7,
         tl.srt = 45, 
         col=brewer.pal(n = 8, name = "PRGn"), 
         sig.level = 0.000001, 
         insig = "blank", 
         diag = FALSE, )

L’étude des corrélations nous permet ici de confirmer une relation entre la variable à expliquer et toutes les variables explicatives définies. Par exemple, elle est négativement corrélée avec le pourcentage de logements vacants.

Elle met aussi à jour de très fortes multicolinéarités, ce qui pose problème. Par exemple, la part de logements suroccupés est fortement corrélée au pourcentage de petits logements. Dans le cadre de la définition d’un modèle linéaire classique, il faudrait sortir du modèle les variables explicatives qui sont trop fortement corrélées. Dans le cadre de cet article nous faisons le choix de conserver toutes les VI, cela fera sens au moment de la GWR.

3.3 Régression linéaire ou méthode des moindres carrés ordinaire (MCO)

# Dans le fonctionnement sur R il est important de stocker la régression dans un objet.
# Pour lancer la régression on va utiliser la fonction lm() 
# dont les 2 lettres sont l'acronyme de linear model
mod.lm <- lm(formula = prix_med ~ 
               perc_log_vac + 
               perc_maison + 
               perc_tiny_log + 
               dens_pop + 
               med_niveau_vis + 
               part_log_suroccup + 
               part_agri_nb_emploi + 
               part_cadre_profintellec_nbemploi, 
             data = data_immo)
# On affiche les principaux résultats avec la fonction summary
summary(mod.lm)

## 
## Call:
## lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
##     dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
##     part_cadre_profintellec_nbemploi, data = data_immo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1566.8  -220.2   -27.2   174.0  3277.3 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1592.873     11.204 142.171  < 2e-16 ***
## perc_log_vac                     -287.337     14.134 -20.330  < 2e-16 ***
## perc_maison                      -128.255     20.041  -6.400 2.22e-10 ***
## perc_tiny_log                    -261.556     27.934  -9.363  < 2e-16 ***
## dens_pop                          173.942     15.272  11.389  < 2e-16 ***
## med_niveau_vis                    288.017     13.860  20.781  < 2e-16 ***
## part_log_suroccup                 388.982     27.352  14.221  < 2e-16 ***
## part_agri_nb_emploi               -20.785     15.059  -1.380    0.168    
## part_cadre_profintellec_nbemploi   -7.841     19.904  -0.394    0.694    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 391.8 on 1214 degrees of freedom
## Multiple R-squared:  0.7784, Adjusted R-squared:  0.7769 
## F-statistic: 532.9 on 8 and 1214 DF,  p-value: < 2.2e-16

mod.lm %>%
  tbl_regression(intercept = TRUE)

GGally::ggcoef_model(mod.lm)

3.4 Diagnostic du modèle linéaire

vif(mod.lm)

##                     perc_log_vac                      perc_maison 
##                         1.577984                         3.204058 
##                    perc_tiny_log                         dens_pop 
##                         6.221631                         1.864236 
##                   med_niveau_vis                part_log_suroccup 
##                         1.532403                         5.977951 
##              part_agri_nb_emploi part_cadre_profintellec_nbemploi 
##                         1.812419                         3.162576

mod.lm %>%
  tbl_regression(intercept = TRUE) %>% add_vif()

score_vif <- vif(mod.lm)
score_vif

##                     perc_log_vac                      perc_maison 
##                         1.577984                         3.204058 
##                    perc_tiny_log                         dens_pop 
##                         6.221631                         1.864236 
##                   med_niveau_vis                part_log_suroccup 
##                         1.532403                         5.977951 
##              part_agri_nb_emploi part_cadre_profintellec_nbemploi 
##                         1.812419                         3.162576

# Création du graphique
par(mar=c(3, 12, 3, 1), mgp=c(1, 0.4, 0)) # pour pouvoir lire les noms des variables
x = barplot(height = score_vif, main = "VIF Values", horiz = TRUE, col = "steelblue", las = 1)
# Ajout du seuil de 4
abline(v = 4, lwd = 3, lty = 2)
# Et de la limite de 3
abline(v = 3, lwd = 3, lty = 2)

mod.lm2 <- lm(formula = prix_med ~ 
                perc_log_vac + 
                perc_maison + 
                dens_pop + 
                med_niveau_vis + 
                part_log_suroccup + 
                part_agri_nb_emploi + 
                part_cadre_profintellec_nbemploi, 
              data = data_immo)

summary(mod.lm2)
## 
## Call:
## lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + 
##     med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
##     part_cadre_profintellec_nbemploi, data = data_immo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1502.2  -226.3   -42.3   173.5  3535.9 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1592.549     11.597 137.329  < 2e-16 ***
## perc_log_vac                     -328.301     13.911 -23.601  < 2e-16 ***
## perc_maison                      -100.010     20.507  -4.877 1.22e-06 ***
## dens_pop                          161.783     15.750  10.272  < 2e-16 ***
## med_niveau_vis                    280.586     14.322  19.591  < 2e-16 ***
## part_log_suroccup                 237.069     22.792  10.401  < 2e-16 ***
## part_agri_nb_emploi                 2.696     15.370   0.175    0.861    
## part_cadre_profintellec_nbemploi  -77.744     19.098  -4.071 4.99e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 405.5 on 1215 degrees of freedom
## Multiple R-squared:  0.7623, Adjusted R-squared:  0.761 
## F-statistic: 556.8 on 7 and 1215 DF,  p-value: < 2.2e-16
vif(mod.lm2)
##                     perc_log_vac                      perc_maison 
##                         1.426783                         3.131469 
##                         dens_pop                   med_niveau_vis 
##                         1.850756                         1.527378 
##                part_log_suroccup              part_agri_nb_emploi 
##                         3.874580                         1.762159 
## part_cadre_profintellec_nbemploi 
##                         2.717634

mod.lm2 %>%
  tbl_regression(intercept = TRUE) %>% add_vif()

GGally::ggcoef_model(mod.lm2)

# L'argument "both" permet d'utiliser les deux méthodes : ascendante et ascendante
step(mod.lm2, direction = "both")

## Start:  AIC=14696.68
## prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + 
##     part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi
## 
##                                    Df Sum of Sq       RSS   AIC
## - part_agri_nb_emploi               1      5060 199817567 14695
## <none>                                          199812507 14697
## - part_cadre_profintellec_nbemploi  1   2725360 202537867 14711
## - perc_maison                       1   3911245 203723752 14718
## - dens_pop                          1  17351111 217163618 14796
## - part_log_suroccup                 1  17792190 217604698 14799
## - med_niveau_vis                    1  63121115 262933622 15030
## - perc_log_vac                      1  91601521 291414028 15156
## 
## Step:  AIC=14694.71
## prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + 
##     part_log_suroccup + part_cadre_profintellec_nbemploi
## 
##                                    Df Sum of Sq       RSS   AIC
## <none>                                          199817567 14695
## + part_agri_nb_emploi               1      5060 199812507 14697
## - part_cadre_profintellec_nbemploi  1   3211704 203029272 14712
## - perc_maison                       1   4155413 203972980 14718
## - dens_pop                          1  17567092 217384659 14796
## - part_log_suroccup                 1  18007058 217824626 14798
## - med_niveau_vis                    1  63117884 262935451 15028
## - perc_log_vac                      1  94962190 294779757 15168

## 
## Call:
## lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + 
##     med_niveau_vis + part_log_suroccup + part_cadre_profintellec_nbemploi, 
##     data = data_immo)
## 
## Coefficients:
##                      (Intercept)                      perc_log_vac  
##                          1592.55                           -327.82  
##                      perc_maison                          dens_pop  
##                           -99.01                            162.05  
##                   med_niveau_vis                 part_log_suroccup  
##                           280.55                            237.44  
## part_cadre_profintellec_nbemploi  
##                           -78.93

mod.lm3 <- lm(formula = prix_med ~ 
                perc_log_vac + 
                perc_maison + 
                dens_pop + 
                med_niveau_vis + 
                part_log_suroccup + 
                part_cadre_profintellec_nbemploi, 
              data = data_immo)

res_modlm <- mod.lm$residuals
datatable(cbind("Nom EPCI" = data_immo$NOM, data.frame("Résidu" = res_modlm))
          %>% format(digits = 2, scientific = FALSE))

par(mfrow=c(1,3))
# Diagramme quantile-quantile qui permet de vérifier l'ajustement
# d'une distribution à un modèle théorique, ici la loi normale
qqPlot(mod.lm)

## [1]  36 266

# Histogramme pour donner une autre indication sur la normalité
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))

# Pour étudier la normalité on peut utiliser le test de Shapiro-Wilk
shapiro.test(mod.lm$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mod.lm$residuals
## W = 0.90792, p-value < 2.2e-16

# Pour évaluer l'homoscédasticité on peut utiliser le test de Breusch-Pagan
# Le package `car` propose une fonction pour le réaliser
ncvTest(mod.lm)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 1151.844, Df = 1, p = < 2.22e-16

# Pour visualiser les individus concernés
datatable(data_immo[c(36, 266),])

# 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)

# Etudier le nouveau modèle
summary(mod.lmx)

## 
## Call:
## lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
##     dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
##     part_cadre_profintellec_nbemploi, data = data_immo, subset = -266)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1622.6  -215.0   -37.5   170.0  2767.7 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1587.246     10.526 150.797  < 2e-16 ***
## perc_log_vac                     -302.182     13.317 -22.692  < 2e-16 ***
## perc_maison                      -132.159     18.814  -7.024 3.57e-12 ***
## perc_tiny_log                    -227.353     26.356  -8.626  < 2e-16 ***
## dens_pop                           -4.725     19.978  -0.237    0.813    
## med_niveau_vis                    284.947     13.012  21.899  < 2e-16 ***
## part_log_suroccup                 404.159     25.701  15.725  < 2e-16 ***
## part_agri_nb_emploi               -17.545     14.138  -1.241    0.215    
## part_cadre_profintellec_nbemploi   23.384     18.841   1.241    0.215    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 367.8 on 1213 degrees of freedom
## Multiple R-squared:  0.7823, Adjusted R-squared:  0.7809 
## F-statistic: 544.9 on 8 and 1213 DF,  p-value: < 2.2e-16

vif(mod.lmx)

##                     perc_log_vac                      perc_maison 
##                         1.589305                         3.168447 
##                    perc_tiny_log                         dens_pop 
##                         6.073006                         2.089258 
##                   med_niveau_vis                part_log_suroccup 
##                         1.531885                         5.726830 
##              part_agri_nb_emploi part_cadre_profintellec_nbemploi 
##                         1.811307                         3.111740

par(mfrow=c(1,2))
plot(rstudent(mod.lmx)) 
qqPlot(mod.lmx)

## [1]  36 180

# Il est possible de comparer les deux modèles et les coefficients
car::compareCoefs(mod.lm, mod.lmx, pvals = TRUE)

## Calls:
## 1: lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
##   dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
##   part_cadre_profintellec_nbemploi, data = data_immo)
## 2: lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
##   dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
##   part_cadre_profintellec_nbemploi, data = data_immo, subset = -266)
## 
##                                  Model 1 Model 2
## (Intercept)                       1592.9  1587.2
## SE                                  11.2    10.5
## Pr(>|z|)                         < 2e-16 < 2e-16
##                                                 
## perc_log_vac                      -287.3  -302.2
## SE                                  14.1    13.3
## Pr(>|z|)                         < 2e-16 < 2e-16
##                                                 
## perc_maison                       -128.3  -132.2
## SE                                  20.0    18.8
## Pr(>|z|)                         1.6e-10 2.1e-12
##                                                 
## perc_tiny_log                     -261.6  -227.4
## SE                                  27.9    26.4
## Pr(>|z|)                         < 2e-16 < 2e-16
##                                                 
## dens_pop                          173.94   -4.72
## SE                                 15.27   19.98
## Pr(>|z|)                         < 2e-16    0.81
##                                                 
## med_niveau_vis                     288.0   284.9
## SE                                  13.9    13.0
## Pr(>|z|)                         < 2e-16 < 2e-16
##                                                 
## part_log_suroccup                  389.0   404.2
## SE                                  27.4    25.7
## Pr(>|z|)                         < 2e-16 < 2e-16
##                                                 
## part_agri_nb_emploi                -20.8   -17.5
## SE                                  15.1    14.1
## Pr(>|z|)                            0.17    0.21
##                                                 
## part_cadre_profintellec_nbemploi   -7.84   23.38
## SE                                 19.90   18.84
## Pr(>|z|)                            0.69    0.21
## 

lm.morantest(model = mod.lm, 
             listw = neighbours_epci_w)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prix_med ~ perc_log_vac + perc_maison +
## perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup +
## part_agri_nb_emploi + part_cadre_profintellec_nbemploi, data =
## data_immo)
## weights: neighbours_epci_w
## 
## Moran I statistic standard deviate = 21, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.359707408     -0.003517610      0.000299161

# Attention : pour avoir le  coefficient il faut faire 1-"Résultat test de Geary"
# (soit ici le coefficient est 0,67)
# Le coefficient de Geary s'étend de 0 à 2, 1 signifiant qu'il n'y a aucune corrélation
# Par ailleurs, un score < 1 implique une corrélation positive 
# et un score > 1 une corrélation négative.
geary(x = data_immo$prix_med, 
      listw = neighbours_epci_w,
      n = length(neighbours_epci), 
      n1 = length(neighbours_epci)-1, 
      S0 = Szero(neighbours_epci_w))

## $C
## [1] 0.3324317
## 
## $K
## [1] 18.24222

data_immo$res_reg <- mod.lm$residuals

# Pour appeler la fonction discr contenue dans le fichier discr.R
source("discr.R")

res_residus <- discr(values = data_immo$res_reg, 
                     center = 0, 
                     pos_center = "class_center", 
                     interval = sd(data_immo$res_reg)*0.5, 
                     min_nb = 10)
breaks_residus <- res_residus[[1]]
nb_cl_sup0_res <- res_residus[[2]]
nb_cl_inf0_res <- res_residus[[3]]

# On fonce légèrement la couleur de fond pour mieux voir la classe centrale
mf_theme("default", bg = "#dedede")
# Création d'une palette divergente avec une couleur neutre centrale
palette = mf_get_pal(n = c(nb_cl_inf0_res, nb_cl_sup0_res), 
                     pal = c("Teal", "Peach"), 
                     neutral = "#f5f5f5")
# Création de la carte
par(mfrow = c(1,1))
mf_map(x = data_immo, 
       var = "res_reg", 
       type = "choro", 
       border = NA,
       lwd = 0.1,
       breaks = breaks_residus,
       pal = palette,
       leg_title = "Valeur centrale =  0\nIntervalle = σ / 2",
       leg_box_border = NA,
       leg_val_rnd = 0,
       leg_adj = c(0, 1.5))
mf_map(x = region_sf, 
       add = TRUE,
       type = "base",
       col = NA,
       border = "#B6B6B4",
       lwd = 1.5)
mf_title("Résidus de régression linéaire classique")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

# réinitialisation du thème
mf_theme("default")

4.1 Niveau global

Pour mesurer l’autocorrélation comme nous l’avons vu précédemment les deux outils les plus utilisés sont les test de Moran et de Geary.

L’indice de Moran va considérer les variances et covariances en tenant compte de la différence entre chaque et la moyenne de toutes les observations. L’indice de Geary de son côté prend en compte la différence entre les observations voisines.

Plus concrètement, l’indice de Moran est une retranscription directe du coefficient de corrélation qu’on peut avoir entre une VI et une VP ; on opère simplement un remplacement de la VP par la valeur moyenne des voisins de la VI. Il s’interprète donc comme un coefficient de corrélation.

Le nombre de voisins étant défini par une distance, il peut être variable en fonction des contextes. On appelle ça le biais de l’effet de bord (Anselin 1995). Le calcul présenté dans cet article prend en compte ce biais en standardisant la matrice de voisinage en ligne.

Concernant l’indice de Geary, il met en rapport d’une part l’écart entre la VI et la moyenne des voisins et d’autre part l’écart entre la VI et la moyenne de l’ensemble des observations. Un indice proche de 0 traduit ainsi une autocorrélation positive, un score proche de 2 une autocorrélation négative. Un score de 1 signifie qu’aucune autocorrélation n’est constatée.

Pourquoi utiliser un indice plutôt qu’un autre ? Ils traduisent sensiblement la même chose, mais diffèrent dans leur construction : Moran est une agrégation de logique locale quand Geary se rapporte à des écarts au global. Cela dit le coût de réalisation n’est pas très élevé et rien n’empêche de faire l’un et l’autre pour voir si les résultats sont cohérents.

Commençons par représenter sur une carte la variable du prix médian des logements par EPCI :

mf_map(x = data_immo, 
       var = "prix_med", 
       type = "choro",
       breaks = "quantile",
       nbreaks = 7,
       pal = "Mint",
       lwd = 0.01,
       leg_title = "Discrétisation par quantile",
       leg_box_border = NA,
       leg_val_rnd = 0,
       leg_adj = c(0, 1.5))
mf_map(x = region_sf, 
       add = TRUE,
       type = "base",
       col = NA,
       border = "#34282C",
       lwd = 1.5)
mf_title("Prix médian de l'immobilier au m² par EPCI")
mf_credits("Sources: Notaires de France 2018, IGN Admin Express 2021")

Une structure spatiale apparaît de manière assez prégnante, ce qui est vérifié par les indicateurs statistiques :

# Pour l'occasion on va standardiser le prix médian
# Cela permettra par la suite de le comparer à d'autres variables 
# si d'autres analyses d'autocorrélation spatiale sont réalisées
data_immo$prix_med_z <- (data_immo$prix_med-mean(data_immo$prix_med))/sd(data_immo$prix_med)

Nous pouvons commencer par réaliser le test de Geary :

# Attention à la lecture particulière des résultats de l'indice de Geary
geary.test(data_immo$prix_med_z, 
           neighbours_epci_w, 
           zero.policy = TRUE, 
           randomisation = FALSE) 

## 
##  Geary C test under normality
## 
## data:  data_immo$prix_med_z 
## weights: neighbours_epci_w   
## 
## Geary C statistic standard deviate = 35.972, p-value < 2.2e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.3324317242      1.0000000000      0.0003444045

On réalise ensuite le test de Moran :

On indique dans un premier temps la variable que l’on souhaite analyser puis la matrice de voisinage.

• L’argument zero.policy = TRUE permet de préciser que l’on souhaite intégrer à l’analyse les entités spatiales qui n’auraient pas de voisins
• L’argument randomisation = FALSE est l’instruction utile pour préciser que nous supposons que la distribution est normale ; dans le cas contraire on devrait partir sur une solution de type Monte-Carlo qui repose sur la randomisation

moran.test(data_immo$prix_med_z, 
           neighbours_epci_w, 
           zero.policy = TRUE, 
           randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  data_immo$prix_med_z  
## weights: neighbours_epci_w    
## 
## Moran I statistic standard deviate = 41.143, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.7154340217     -0.0008183306      0.0003030652

Ce test d’autocorrélation se lit exactement comme un test de corrélation classique. On va donc s’intéresser au signe, à la grandeur du coefficient et à la p-value du test.

Ici on a donc une autocorrélation spatiale positive et importante.

Qu’implique l’existence de cette autocorrélation ?

Comme nous l’avons mentionné, l’autocorrélation spatiale, positive ou négative, décrit la relation d’une variable avec elle-même du fait de la localisation spatiale des observations.

Cette autocorrélation spatiale peut donc être :

• Structurelle : les valeurs s’influencent dans l’espace et créent ainsi une structure spatiale sous-jacente. On contrevient ainsi à l’hypothèse d’indépendance statistique des individus. Cette autocorrélation peut être constatée autant sur une VI que sur des VP.
• Contextuelle : le résultat d’un phénomène inobservé (un aléa) ou qu’on ne peut mesurer efficacement. Ce phénomène s’inscrit dans l’espace et de ce fait crée une structure spatiale. Les phénomènes de diffusion par exemple : il existe différents phénomènes sociaux de la sorte comme par exemple des phénomènes d’interaction ou de diffusion (comme les phénomènes de diffusion technologique). Ces phénomènes peuvent produire de l’autocorrélation spatiale.
• Symptomatique : dans la définition du modèle statistique, la mesure de l’autocorrélation spatiale peut être envisagée comme un outil de diagnostic et de détection d’un “mauvais” (du point de vue statistique) modèle (variables non intégrées spatialement corrélées, erreurs sur le choix de l’échelle à laquelle le phénomène spatial est analysé, etc.)

Pour visualiser rapidement la structure spatiale on peut aussi réaliser un diagramme de Moran qui est complémentaire au test statistique.

moran.plot(data_immo$prix_med_z, neighbours_epci_w, 
           labels = TRUE, 
           xlab = "prix medians centrés réduits par epci" , 
           ylab = "moyenne du prix médian centré réduit par epci des voisins")

Comment lire ce diagramme ?

Il présente le nuage de points entre les valeurs d’une variable et les valeurs moyennes des voisins de chaque valeur de la variable initiale. Si les individus sont répartis aléatoirement dans l’espace alors c’est le signe d’une absence d’autocorrélation, la pente de la droite de régression sera donc de 0. Si au contraire la pente est non nulle, comme c’est le cas ici, c’est le signe de la présence d’une autocorrélation spatiale.

Un aspect important de ce diagramme est qu’il permet d’ores et déjà de caractériser les coefficients d’autocorrélation spatiale. Comme ce graphique présente des valeurs centrées-réduites, les valeurs sont centrées sur la moyenne qui est de 0. De ce fait, tous les points à droite de l’axe des ordonnées pour x = 0 auront une moyenne > 0 et ceux à gauche < 0. Cette réflexion s’applique également à l’axe des abscisses. Par convention on désigne les individus avec une moyenne > 0 par le terme high et ceux avec une moyenne < 0 par le terme low, au sens de supérieur ou inférieur à la moyenne.

Ainsi on peut découper ce diagramme en 4 quadrants. Les quadrants en haut à droite et en bas à gauche correspondent aux individus ayant une autocorrélation spatiale positive, c’est-à-dire une valeur proche de celle de leurs voisins.

• Pour le quadrant en haut à droite on parle du quadrant High-High composé d’individus ayant une valeur de la variable plus élevée que la moyenne entourés d’autres individus qui leur ressemblent.
• Pour le quadrant en bas à gauche on parle du quadrant Low-Low composé d’individus avec un score plus faible que la moyenne avec des voisins avec un score similaire.
• Le quadrant en bas à droite est considéré comme le quadrant High-Low. On y retrouve des individus avec un score plus élevé que la moyenne sur la variable du prix médian mais avec un voisinage qui ne lui ressemble pas : autocorrélation spatiale négative mais score élevé à la variable.
• Enfin, en haut à gauche on retrouve à l’inverse les individus avec une valeur du prix médian plus faible que la moyenne et un indice d’autocorrélation spatiale négatif. C’est le quadrant Low-High.

Diagramme de Moran avec découpage en 4 quadrants

4.2 Niveau local

Les indices de Geary et Moran présentent une limite importante : ils partent du principe d’un processus spatial stationnaire, autrement dit global. Ainsi, l’effet de dépendance serait le même dans tout notre espace et la mesure fournie ne permet pas de différencier les contextes spatiaux locaux. De plus, l’effet de la dimension serait le même dans tout notre espace ; ce qui pose question et ce d’autant plus que l’emprise géographique augmente.

La compréhension de ces contextes locaux passe par la réalisation 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) (Anselin 1995).

D’après Luc Anselin, le LISA de chaque individu statistique indique l’intensité du regroupement spatial de valeurs similaires autour de cet individu. Autrement dit, un individu avec un LISA élevé va avoir une concentration autour de lui de voisins avec des valeurs similaires (pour nous par exemple un regroupement d’individus avec un prix particulièrement élevé ou à l’inverse particulièrement bas).

Le LISA le plus utilisé est le I de Moran local.

L’idée est de faire appel ici au LISA pour compléter la compréhension 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é spatiale, c’est-à-dire qui ne suivent pas le processus global.

Nous allons calculer le I de Moran local sur nos données, grâce au package rgeoda développé également par Luc Anselin pour réaliser sur R les traitements de GeoDa.

# Calcul I de Moran local

# Pour utiliser la fonction local_Moran du package rgeoda, 2 pré-requis :

# 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"]

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
lisa_colors <- lisa_colors(lisa)
lisa_labels <- lisa_labels(lisa)
lisa_clusters <- lisa_clusters(lisa)
lisa_value <- lisa_values(lisa)
lisa_pvalue <- lisa_pvalues(lisa)

Pour illustrer les Moran locaux on peut réaliser une carte :

# On récupère les moran locaux dans data_immo
data_immo$moranlocalvalue <- lisa_values(lisa)
# Discrétisation standardisée avec des classes de un demi-écart-type
res_moranloc <- discr(values = data_immo$moranlocalvalue, 
                      center = 0, 
                      pos_center = "class_break", 
                      interval = sd(data_immo$moranlocalvalue)/2, 
                      min_nb = 10)
breaks_moranloc <- res_moranloc[[1]]
nb_cl_sup0_moranloc <- res_moranloc[[2]]
nb_cl_inf0_moranloc <- res_moranloc[[3]]
# Création d'une palette associée
palette = mf_get_pal(n = c(nb_cl_inf0_moranloc, nb_cl_sup0_moranloc), 
                     pal = c("Teal", "Reds"))
# Création de la carte
mf_map(x = data_immo, 
       var = "moranlocalvalue", 
       type = "choro",
       breaks = breaks_moranloc,
       border = "gray",
       lwd = 0.2,
       pal = palette,
       leg_title = "Discrétisation standardisée\nvaleur centrale = 0\nintervalle = σ / 2", 
       leg_box_border = "gray",
       leg_val_rnd = 1,
       leg_adj = c(0, 1.5))
mf_map(x = region_sf, 
       add = TRUE,
       type = "base",
       col = NA,
       border = "gray",
       lwd = 1.5)
mf_title("Carte des LISA du prix médian du logement (Moran local)")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

Si le score du Moran local est > 0 cela implique que l’on a un regroupement de valeurs similaires, plus faibles ou plus élevées que la moyenne. En revanche si le score est < 0 alors on a un regroupement de valeurs dissimilaires, par exemple des valeurs plus faibles entourées de valeurs plus fortes (centre de la Corse).

L’étude des p-values associées est également importante pour l’interprétation : une p-value < 0,05 des LISA apparaît comme significative et renvoie à des clusters de valeurs (similaires ou dissimilaires) qui sont plus marqués que ce l’on pourrait observer si la répartition spatiale du phénomène était considérée comme aléatoire.

# Carte des p-values des Moran locaux
data_immo$moranlocalpvalue<- lisa_pvalues(lisa)

# Pour plus de lisibilité dans la carte on va faire des classes des p-values
data_immo <- data_immo %>% mutate(lisapvalue_fac = case_when(moranlocalpvalue <= 0.002 ~ "[0.001 - 0.002[",
                           moranlocalpvalue <= 0.01 ~ "[0.002 - 0.01[",
                           moranlocalpvalue <= 0.05 ~ "[0.01 - 0.05[",
                           TRUE ~ "[0.05;0.5]")) %>%
  mutate(lisapvalue_fac = factor(lisapvalue_fac,
                                 levels = c("[0.001 - 0.002[", 
                                            "[0.002 - 0.01[", 
                                            "[0.01 - 0.05[", 
                                            "[0.05 - 0.5]")))

mypal <- mf_get_pal(n = 4, palette = "Reds")

mf_map(x = data_immo, 
       var = "lisapvalue_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.08, 
       pal = mypal, 
       leg_title = "p-value Moran local",
       leg_box_border = "grey3",
       leg_adj = c(0, 1.5))
mf_map(x = region_sf, 
       add = TRUE,
       type = "base",
       col = NA,
       border = "gray",
       lwd = 1.5)
mf_title("Carte de significativité des LISA")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

Les regroupements que l’on observe vont pouvoir se rapprocher des 4 types que nous avions sur le diagramme de Moran (High-High, Low-Low, High-Low et Low-High).

Cette carte nous apporte donc une information supplémentaire sur les LISA en différenciant, parmi les clusters, les valeurs plus faibles ou plus fortes que la moyenne. On voit, sur cette carte, se dégager la diagonale du vide qui présente un cluster important des prix médians de l’immobilier faibles, qui se différencie des la Région parisienne, de la Côte d’Azur, de la Corse ou des Alpes du nord qui présentent des clusters de prix médians élevés..

data_immo$testmoran <- sapply(lisa_clusters, function(x){return(lisa_labels[[x+1]])})

colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")

mf_map(x = data_immo, 
       var = "testmoran", 
       type = "typo", 
       border = "black", 
       lwd = 0.1, 
       pal= colors,
       val_order = c("Non significatif","Low-Low","Low-High","High-Low","High-High"),
       leg_title = "Lisa cluster",
       leg_adj = c(0, 1.5))
mf_map(x = region_sf, 
       add = TRUE,
       type = "base",
       col = NA,
       border = "#34282C",
       lwd = 1.5)
mf_title("Regroupement des LISA")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

Si on est en présence d’une autocorrélation spatiale au niveau global, les LISA pourront aussi servir d’indicateur d’instabilité locale. Ils indiquent à la fois des clusters locaux qui vont avoir un impact fort sur le processus spatial global (un score d’autocorrélation locale plus marqué que l’autocorrélation globale) ou à l’inverse des zones qui se démarquent du processus global (plus faible autocorrélation).

En revanche, s’il n’y a pas d’autocorrélation spatiale au niveau global les LISA vont nous permettre de détecter des zones où des valeurs semblables se regroupent de façon significative. Ils font émerger des structures au niveau local au sein desquelles les liens entre voisins seront particulièrement importants. On peut imaginer certaines configurations où l’autocorrélation spatiale globale ne paraisse pas significative, mais où ce score masque localement des clusters marqués qui restreignent la qualité des modèles de régression classique.

5.1 Calcul de la matrice des distances

La première étape est de calculer la distance entre toutes les observations. Pour ce faire nous utiliserons la fonction gw.dist().

# Le package GWmodel n'est pas compatible avec le format sf il a besoin d'un objet sp
# (contrairement à spgwr qui peut travailler avec un format sf)

# Pour construire la matrice de distances entre centroïdes des EPCI :
dm.calib <- gw.dist(dp.locat = coordinates(data_immo_sp))

5.2 Définition de la bande passante

La bande passante est une distance au-delà de laquelle le poids des observations est considéré comme nul. Le calcul de cette distance est très important car la valeur de la bande passante pourra fortement influencer notre modèle. La définition de la bande passante renvoie au type de pondération que nous souhaitons appliquer. Heureusement la fonction bw.gwr va choisir pour nous le résultat optimal…

Pour ce faire la fonction va se baser sur un critère statistique que l’utilisateur devra définir : le CV (Cross Validation ou validation croisée) ou le AIC (Akaike Information Criterion ou critère d’information d’Akaike). Elle reposera aussi sur un type de noyau qu’il faudra également définir : gaussien, exponentiel, bicarré, tricube ou encore Boxcar. Enfin, il sera également nécessaire de savoir si ce noyau pourra être adaptatif ou fixe.

Voici quelques informations pour guider nos choix :

• Le critère CV a pour objectif de maximiser le pouvoir prédictif du modèle, le critère AIC va chercher un compromis entre le pouvoir prédictif du modèle et son degré de complexité. En général, le critère AIC est privilégié.
• Avec un noyau fixe l’étendue du noyau est déterminée par la distance au point d’intérêt et il est identique en tout point de l’espace. Un noyau fixe est adapté si la répartition des données est homogène dans l’espace, l’unité de la bande passante sera donc une distance. Avec un noyau adaptatif l’étendue du noyau est déterminée par le nombre de voisins. Il est donc plus adapté à une répartition non homogène, l’unité sera alors le nombre de voisins.

Concernant la forme des noyaux :

• Les noyaux gaussiens et exponentiels vont pondérer toutes les observations avec un poids qui tend vers zéro avec la distance au point estimé.
• Les noyaux bisquare et tricube (dont les formes sont très proches) accordent également aux observations un poids décroissant avec la distance, mais par contre ce poids est nul au-delà de la distance définie par la bande passante.
• Le noyau Boxcar traite un phénomène continu de façon discontinue.

Manuel de géographie quantitative (Feuillet et al. 2019)

Sachant que, sur la forme du noyau, il est tout à fait possible de comparer deux pondérations et deux modèles de GWR.

Ici, nous allons tester avec un noyau gaussien, ce qui sera justifié un petit peu plus bas.

Pour éviter les répétitions la formule est stockée dans une variable au préalable.

formula <- prix_med ~ 
  perc_log_vac +
  perc_maison + 
  perc_tiny_log + 
  dens_pop + 
  med_niveau_vis + 
  part_log_suroccup + 
  part_agri_nb_emploi + 
  part_cadre_profintellec_nbemploi

# Définition de la bande passante (bandwidth en anglais) :
bw_g <- bw.gwr(data = data_immo_sp, 
              approach = "AICc", 
              kernel = "gaussian", 
              adaptive = TRUE, 
              dMat = dm.calib,
              formula = formula)

## Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17959.04 
## Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17884.14 
## Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17798.07 
## Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17715.24 
## Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17622.76 
## Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17526.89 
## Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 17422.28 
## Adaptive bandwidth (number of nearest neighbours): 45 AICc value: 17353.68 
## Adaptive bandwidth (number of nearest neighbours): 33 AICc value: 17283.09 
## Adaptive bandwidth (number of nearest neighbours): 29 AICc value: 17261.83 
## Adaptive bandwidth (number of nearest neighbours): 23 AICc value: 17250.23 
## Adaptive bandwidth (number of nearest neighbours): 22 AICc value: 17247.23 
## Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79 
## Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79

bw_g

## [1] 19

La bande passante est donc ici de 19 voisins, ce qui implique que les EPCI au-delà de cette “distance” auront un poids ramené à 0 et ne joueront donc plus de rôle dans la description de la relation statistique.

5.3 Estimation du modèle

Une fois la bande passante définie on peut lancer l’estimation du modèle de GWR.

# un des problèmes de la GWR est de gérer des individus "aberrants" au niveau local
# 2 méthodes ont été définies pour gérer cela :
# Méthode 1 (filtered = TRUE) on filtre en fonction des individus standardisés
# L'objectif est de détecter les individus dont les résidus sont très élevés et de les exclure.
# Méthode 2 (filtered = FALSE) on diminue le poids des observations aux résidus élevés.
mod.gwr_g <- gwr.robust(data = data_immo_sp, 
                        dMat = dm.calib,
                        bw = bw_g,
                        kernel = "gaussian",
                        filtered = FALSE,
                        adaptive = TRUE,
                        formula = formula)

Si on souhaite comparer deux modèles car nous avons un doute sur les paramètres c’est tout à fait possible. Par exemple ici nous souhaitons comparer deux formes de noyau :

bw_tri <- bw.gwr(data = data_immo_sp, 
              approach = "AICc", 
              kernel = "tricube", 
              adaptive = TRUE, 
              dMat = dm.calib,
              formula = formula)

## Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17701.17 
## Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17586.94 
## Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17422.51 
## Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17294.75 
## Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17254.1 
## Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17279.33 
## Adaptive bandwidth (number of nearest neighbours): 154 AICc value: 17259.05 
## Adaptive bandwidth (number of nearest neighbours): 112 AICc value: 17257.17 
## Adaptive bandwidth (number of nearest neighbours): 138 AICc value: 17254.8 
## Adaptive bandwidth (number of nearest neighbours): 122 AICc value: 17253.75 
## Adaptive bandwidth (number of nearest neighbours): 117 AICc value: 17255.04 
## Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57 
## Adaptive bandwidth (number of nearest neighbours): 126 AICc value: 17254.25 
## Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57

mod.gwr_tri <- gwr.robust(data = data_immo_sp, 
                   dMat = dm.calib,
                   bw = bw_tri,
                   kernel = "gaussian",
                   filtered = FALSE,
                   adaptive = TRUE,
                   formula = formula)

Best_gwr <- cbind(
  rbind(bw_g, bw_tri),
  rbind(mod.gwr_g$GW.diagnostic$gw.R2,mod.gwr_tri$GW.diagnostic$gw.R2),
  rbind(mod.gwr_g$GW.diagnostic$AIC,mod.gwr_tri$GW.diagnostic$AIC)) %>% 
  `colnames<-`(c("Nb Voisins","R2","AIC")) %>% 
  `rownames<-`(c("GAUSSIAN","TRICUBE"))

Best_gwr

##          Nb Voisins        R2      AIC
## GAUSSIAN         19 0.9324328 16849.26
## TRICUBE         123 0.8577370 17567.05

Le modèle utilisant un noyau gaussien a un meilleur pouvoir explicatif (\(R^2\)) et un meilleur score de qualité (\(AIC\)). Ce modèle est donc plus performant et dans notre situation c’est plutôt ce modèle qu’il faut privilégier.

5.4 Interprétation des premiers résultats

Comme pour le modèle linéaire classique, l’objet qui contient la GWR est composé de plusieurs éléments. Pour obtenir les résultats il suffit d’appeler l’objet.

# Pour voir les différents éléments qui composent notre modèle de GWR
summary(mod.gwr_g)

##               Length Class                    Mode
## GW.arguments    11   -none-                   list
## GW.diagnostic    8   -none-                   list
## lm              14   lm                       list
## SDF           1223   SpatialPolygonsDataFrame S4  
## timings          5   -none-                   list
## this.call       15   -none-                   call
## Ftests           0   -none-                   list

# Pour accéder aux résultats
mod.gwr_g

##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2024-09-25 17:54:21.062216 
##    Call:
##    gwr.basic(formula = formula, data = data, bw = bw, kernel = kernel, 
##     adaptive = adaptive, p = p, theta = theta, longlat = longlat, 
##     dMat = dMat, F123.test = F123.test, cv = T, W.vect = W.vect, 
##     parallel.method = parallel.method, parallel.arg = parallel.arg)
## 
##    Dependent (y) variable:  prix_med
##    Independent variables:  perc_log_vac perc_maison perc_tiny_log dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi part_cadre_profintellec_nbemploi
##    Number of data points: 1223
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##     Min      1Q  Median      3Q     Max 
## -1566.8  -220.2   -27.2   174.0  3277.3 
## 
##    Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
##    (Intercept)                      1592.873     11.204 142.171  < 2e-16 ***
##    perc_log_vac                     -287.337     14.134 -20.330  < 2e-16 ***
##    perc_maison                      -128.255     20.041  -6.400 2.22e-10 ***
##    perc_tiny_log                    -261.556     27.934  -9.363  < 2e-16 ***
##    dens_pop                          173.942     15.272  11.389  < 2e-16 ***
##    med_niveau_vis                    288.017     13.860  20.781  < 2e-16 ***
##    part_log_suroccup                 388.982     27.352  14.221  < 2e-16 ***
##    part_agri_nb_emploi               -20.785     15.059  -1.380    0.168    
##    part_cadre_profintellec_nbemploi   -7.841     19.904  -0.394    0.694    
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 391.8 on 1214 degrees of freedom
##    Multiple R-squared: 0.7784
##    Adjusted R-squared: 0.7769 
##    F-statistic: 532.9 on 8 and 1214 DF,  p-value: < 2.2e-16 
##    ***Extra Diagnostic information
##    Residual sum of squares: 186354789
##    Sigma(hat): 390.6721
##    AIC:  18086.13
##    AICc:  18086.31
##    BIC:  16985.31
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: gaussian 
##    Adaptive bandwidth: 19 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: A distance matrix is specified for this model calibration.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                                            Min.     1st Qu.      Median
##    Intercept                         1109.13394  1309.99018  1516.79508
##    perc_log_vac                      -919.82407  -224.60028  -163.40316
##    perc_maison                       -898.16365  -258.64481  -110.20592
##    perc_tiny_log                    -1210.96882  -206.68578   -92.91874
##    dens_pop                          -411.20172    72.82448   217.03593
##    med_niveau_vis                      81.29015   287.78992   358.32914
##    part_log_suroccup                 -543.68600    42.42839   142.48859
##    part_agri_nb_emploi               -498.74706   -65.81443   -32.88823
##    part_cadre_profintellec_nbemploi  -347.37935   -48.57329     0.58811
##                                         3rd Qu.    Max.
##    Intercept                         1696.07367 2330.78
##    perc_log_vac                      -113.83881  388.64
##    perc_maison                         76.01079  896.95
##    perc_tiny_log                       32.81544  773.30
##    dens_pop                           268.35651  659.84
##    med_niveau_vis                     413.36560  853.98
##    part_log_suroccup                  247.33650  700.50
##    part_agri_nb_emploi                 -1.15602  120.33
##    part_cadre_profintellec_nbemploi    49.93194  388.50
##    ************************Diagnostic information*************************
##    Number of data points: 1223 
##    Effective number of parameters (2trace(S) - trace(S'S)): 320.246 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 902.754 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 17201.79 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 16849.26 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 17068.01 
##    Residual sum of squares: 56808914 
##    R-square value:  0.9324328 
##    Adjusted R-square value:  0.9084372 
## 
##    ***********************************************************************
##    Program stops at: 2024-09-25 17:54:32.206493

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 \(R^2 ajusté\) de la GWR, qui est nettement meilleur que celui de la régression linéaire multiple. On passe de \(77\ %\) de variance expliquée à \(91\ %\) avec la GWR. Ce \(R^2\), pour la GWR, est en fait une moyenne calculée des \(R^2\) de toutes les régressions locales réalisées par la GWR.

La seconde information qui nous intéresse particulièrement est les coefficients associés aux VI. Nous voyons qu’ils ne sont pas présentés de la même manière que ceux de la régression linéaire. En effet, chaque VI va avoir des coefficients en fonction du minimum, du maximum et 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é 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, quand le pourcentage baisse d’une unité le prix médian augmente de \(287€\). Dans le cas de la densité de population on a un coefficient global positif donc une relation positive. La densité augmente donc le prix médian augmente. Ici, au global, quand la densité augmente d’une unité le prix médian augmente de \(173€\).

Les résultats de la GWR peuvent donc être lus à une échelle globale pour mesurer la pertinence du modèle ; mais également à des échelles locales : les résultats illustrent ainsi comment les coefficients varient en fonction des unités spatiales. Gardons l’exemple de la densité de population. Dans les lieux où le coefficient est à son minimum, c’est-à-dire \(-411\), on a une relation négative. Dans ces espaces, quand la densité augmente d’un écart-type le prix médian baisse de \(411€\).

Ensuite nous pouvons constater une inversion du signe du coefficient. Ainsi dans les EPCI du dernier quartile où le prix médian du logement est le plus élevé (par ex. Paris) le coefficient est positif. À son maximum une augmentation d’une unité entraîne une augmentation du prix de \(659€\). On a donc très clairement un effet de la densité qui ne sera pas du tout le même en fonction du lieu.

Nous pouvons également étudier l’intervalle interquartile. Ainsi, toujours pour la densité, ce résultat implique que pour 50 % des unités spatiales (EPCI entre quartile 1 et 3), une augmentation d’une unité de la densité va impliquer une augmentation du prix médian entre \(72€\) et \(268€\).

La cartographie va être la meilleure manière de représenter les betas (coefficients) et les différents indicateurs fournis avec la GWR, cela nous permet de décrire plus finement et plus précisément les phénomènes observés.

L’ensemble des données est stocké dans le sous-objet SDF du modèle. Il contient l’ensemble des informations du modèle associé à chaque donnée spatiale.

On peut le convertir en un dataframe pour le visualiser plus facilement. À l’origine il est au format “SpatialPointsDataFrame”.

# L'ensemble des données est stocké dans le sous objet SDF de notre modèle
# Pour voir à quoi il ressemble avec les noms des EPCI
datatable(cbind("Nom EPCI" = data_immo$NOM, mod.gwr_g$SDF@data)
          %>% format(digits = 2, scientific = FALSE))

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

# Pour voir les variables qui le constituent
names(mod.gwr_g$SDF@data)

##  [1] "Intercept"                           "perc_log_vac"                       
##  [3] "perc_maison"                         "perc_tiny_log"                      
##  [5] "dens_pop"                            "med_niveau_vis"                     
##  [7] "part_log_suroccup"                   "part_agri_nb_emploi"                
##  [9] "part_cadre_profintellec_nbemploi"    "y"                                  
## [11] "yhat"                                "residual"                           
## [13] "CV_Score"                            "Stud_residual"                      
## [15] "Intercept_SE"                        "perc_log_vac_SE"                    
## [17] "perc_maison_SE"                      "perc_tiny_log_SE"                   
## [19] "dens_pop_SE"                         "med_niveau_vis_SE"                  
## [21] "part_log_suroccup_SE"                "part_agri_nb_emploi_SE"             
## [23] "part_cadre_profintellec_nbemploi_SE" "Intercept_TV"                       
## [25] "perc_log_vac_TV"                     "perc_maison_TV"                     
## [27] "perc_tiny_log_TV"                    "dens_pop_TV"                        
## [29] "med_niveau_vis_TV"                   "part_log_suroccup_TV"               
## [31] "part_agri_nb_emploi_TV"              "part_cadre_profintellec_nbemploi_TV"
## [33] "E_weigts"                            "Local_R2"

# Intercept : c'est la constante c'est-à-dire le prix médian de référence
# nom de la variable : estimation du coefficient, du beta associé à la VI en chaque point.
# y : valeurs de la VD
# yhat : valeur de y prédite.
# residual, Stud_residual : résidu et résidu standardisé
# CV_score : score de validation croisée
# _SE : erreur standard de l'estimation du coefficient
# _TV : t-value de l'estimation du coefficient
# E_weight : poids des observations dans la régression robuste
# Local_R2 : R ²au niveau de chaque unité spatiale

# On récupère les résidus dans data_immo
res_gwr <- mod.gwr_g$SDF$Stud_residual
data_immo$res_gwr <- res_gwr
# calcul des limites de classes avec la fonction discr, centrées sur 0
res_resgwr <- discr(values = data_immo$res_gwr, 
                    center = 0, 
                    pos_center = "class_center", 
                    interval = sd(data_immo$res_gwr)*0.5, 
                    min_nb = 10)
breaks_gwr <- res_resgwr[[1]]
nb_cl_sup0_gwr <- res_resgwr[[2]]
nb_cl_inf0_gwr <- res_resgwr[[3]]
# Création de la palette correspondante
palette = mf_get_pal(n = c(nb_cl_inf0_gwr, nb_cl_sup0_gwr), 
                     pal = c("Teal", "Peach"), 
                     neutral = "#f5f5f5")
# Création de la carte des résidus
mf_map(x = data_immo, 
       var = "res_gwr", 
       type = "choro", 
       border = "gray", 
       lwd = 0.2, 
       breaks = breaks_gwr,
       pal = palette, 
       leg_title = "Discrétisation standardisée :\nvaleur centrale = 0\nintervalle = σ / 2",
       leg_box_border = "gray",
       leg_val_rnd = 1,
       leg_adj = c(0, 1.5))
mf_map(x = region_sf, 
       add = TRUE,
       type = "base",
       col = NA,
       border = "gray",
       lwd = 1.5)
mf_title("Résidus de la GWR")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

moran.test(data_immo$res_gwr,
           neighbours_epci_w,
           zero.policy = TRUE,
           randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  data_immo$res_gwr  
## weights: neighbours_epci_w    
## 
## Moran I statistic standard deviate = 1.7079, p-value = 0.04383
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0289140229     -0.0008183306      0.0003030652

5.4.2 Étude des coefficients

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
data_immo$dens_pop.coef <- mod.gwr_g$SDF$dens_pop
data_immo$med_vie.coef <- mod.gwr_g$SDF$med_niveau_vis
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 :

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)
  mf_title(noms_coefs[var])
}

Coefficients de la GWR (Sources : Notaires de France 2018, INSEE 2019, IGN Admin Express 2021)

Les cartes des betas vont illustrer la variation des effets en fonctions des entités spatiales et de leur voisinage. Dans notre cas on verra quels sont les EPCI où l’effet du coefficient est négatif et ceux où il est positif, c’est-à-dire dans quel EPCI la VI va entraîner une augmentation du prix médian et dans quel autre au contraire une diminution, toutes choses égales par ailleurs. Sachant que dans notre cas toutes les VI sont significatives, elles ont donc toutes un effet qui varie localement.

Au-delà de cette visualisation VI par VI, il peut être intéressant de savoir par EPCI quelle variable sera la plus explicative dans la relation à la VD et laquelle a l’impact le plus important. Nous avons donc réalisé une carte des contributions maximales par EPCI. Pour la réaliser nous nous sommes basés sur le t-value.

data_immo$agri.t <- mod.gwr_g$SDF$part_agri_nb_emploi_TV
data_immo$maison.t <- mod.gwr_g$SDF$perc_maison_TV
data_immo$dens.t <- mod.gwr_g$SDF$dens_pop_TV
data_immo$medvie.t <- mod.gwr_g$SDF$med_niveau_vis_TV
data_immo$logvac.t <- mod.gwr_g$SDF$perc_log_vac_TV
data_immo$tinylog.t <- mod.gwr_g$SDF$perc_tiny_log_TV
data_immo$suroccup.t <- mod.gwr_g$SDF$part_log_suroccup_TV
data_immo$cadre.t <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi_TV     

# Définir contribrution maxixmale
df <- as.data.frame(data_immo)
# On passe les t-values en valeur absolue pour voir la plus grande
# contribution dans un 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))
# Création de la carte
mf_map(x = data_immo, 
       var = "contribmax", 
       type = "typo", 
       pal = brewer.pal(6,'Set2'),
       border = "white",
       lwd = 0.2,
       leg_pos = NA)
mf_map(x = region_sf, 
       add = TRUE,
       type = "base",
       col = NA,
       border = "#736F6E",
       lwd = 1.5)
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")

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). L’inclusion d’une VI a pour conséquence d’en exclure d’autres qui présenteraient trop de multicolinéarité. Le phénomène mis en évidence ici repose sur le fait qu’en fonction du lieu, de la hiérarchie des VI et donc de leur pertinence au sein du modèle, n’est pas constante. Cette visualisation des données ouvre une perspective intéressante pour la suite. Il serait tout à fait pertinent de développer une méthode permettant d’effectuer un stepwise préalable à toute GWR valable pour chaque lieu individuellement. Un tel modèle traiterait dans sa totalité le problème de non-stationnarité.

Nous pouvons également cartographier les \(R^2\) locaux, ce qui fournit une indication sur les zones où la variabilité est la mieux expliquée.

data_immo$r2local <- mod.gwr_g$SDF$Local_R2

mf_map(x = data_immo, 
       var = "r2local", 
       type = "choro",
       breaks = "quantile",
       nbreaks = 11,
       pal= "Reds",
       border = "gray",
       lwd = 0.2,
       leg_title = "Discrétisation par quantile",
       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("R² locaux")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

Nous avons choisi d’utiliser une palette de valeurs large pour représenter au mieux les différences de \(R^2\) locaux. Toutefois, cette carte peut paraître trompeuse ; tous les EPCI obtiennent une explication satisfaisante avec un R² local minimum de 0.68.

À partir des t-values on peut aussi étudier la significativité des effets sur le territoire. On peut ainsi calculer et cartographier un indicateur qui représenterait le nombre de VI dont l’effet est significatif sur chaque unité spatiale. Cela donne une bonne idée de la complexité du phénomène sur un espace donné (en effet sur un EPCI on peut avoir toutes les variables significatives, elle jouent donc sur cet espace toutes un rôle) et souligne l’importance d’avoir une carte par coefficient. Cette carte montre également qu’un modèle parfaitement adaptatif gagnerait en sobriété et donnerait un AIC plus satisfaisant en ne sélectionnant localement que les variables réellement significatives dans la relation.

# Pour rappel si on a plus de 200 individus et le t-value > |1.96| 
# on pourra considérer le coefficient comme significatif au seuil de 0.05 
# (95 % chances que ce ne soit pas dû au hasard)
data_immo$nbsignif_t <- rowSums(abs(df[, c(30:37)]) > 1.96)

mf_map(x = data_immo, 
       var = "nbsignif_t", 
       type = "typo",
       pal = "Reds",
       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,
       type = "base",
       col = NA,
       border = "gray",
       lwd = 1.5)
mf_title("Nombre de variables expliquant significativement la variabilité des prix médians par EPCI (t-value)")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

Il se peut que cela soit plus intéressant d’utiliser les p-values, notamment si vous avez moins de 200 individus. Dans le cas présent, les 2 cartes sont identiques puisque nous avons 1 223 individus.

# Les p-values ne sont pas fournis dans le modèle de la GWR
# On pourrait les calculer à partir de t-value et de l'erreur standard 
# mais le package GWmodel propose une fonction pour les obtenir
pvalue <- gwr.t.adjust(mod.gwr_g)

# On ajoute les p-value à notre fichier
data_immo$agri.p <- pvalue$SDF$part_agri_nb_emploi_p 
data_immo$maison.p <- pvalue$SDF$perc_maison_p
data_immo$dens.p <- pvalue$SDF$dens_pop_p
data_immo$medvie.p <- pvalue$SDF$med_niveau_vis_p
data_immo$logvac.p <- pvalue$SDF$perc_log_vac_p
data_immo$tinylog.p <- pvalue$SDF$perc_tiny_log_p
data_immo$suroccup.p <- pvalue$SDF$part_log_suroccup_p
data_immo$cadre.p <- pvalue$SDF$part_cadre_profintellec_nbemploi_p

df <- as.data.frame(data_immo)
data_immo$nbsignif_p <- rowSums(df[, c(41:48)] < 0.05)

mf_map(x = data_immo, 
       var = "nbsignif_p", 
       type = "typo",
       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, 
       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")

Dans ce cadre, il est possible de réaliser une collection de cartes des p-values (ou t-values) comme ce qui a été fait pour les coefficients. L’intérêt est de voir où l’effet de la VI est significatif et où il ne l’est pas.

# Ici nous représenterons les p-values avec un découpage par classe de significativité
# et seulement les p-values de 2 VI
par(mfrow = c(1, 2))
# Par exemple les p-values des coefficients de la variable part de l'emploi agriculteur
data_immo <- data_immo %>%  mutate(agri.p_fac = case_when(agri.p<= 0.002 ~ "[0;0.002[",
                           agri.p <= 0.01 ~ "[0.002;0.01[",
                           agri.p <= 0.05 ~ "[0.01;0.05[",
                           agri.p <= 0.1 ~ "[0.05;0.1[",
                           TRUE ~ "[0.1;1]")) %>%
  mutate(agri.p_fac = factor(agri.p_fac,
                        levels = c("[0;0.002[", 
                                   "[0.002;0.01[",
                                   "[0.01;0.05[",
                                   "[0.05;0.1[", 
                                   "[0.1;1]")))

mypal2 <- mf_get_pal(n = 5, palette = "OrRd")
mf_map(x = data_immo, 
       var = "agri.p_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       leg_box_border = "grey3",
       pal = mypal2, 
       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 du % agriculteurs")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

# Pour la densité de population
data_immo <- data_immo %>%  mutate(dens.p_fac = case_when(dens.p <= 0.002 ~ "[0;0.002[",
                           dens.p <= 0.01 ~ "[0.002;0.01[",
                           dens.p <= 0.05 ~ "[0.01;0.05[",
                           dens.p <= 0.1 ~ "[0.05;0.1[",
                           TRUE ~ "[0.1;1]")) %>%
  mutate(dens.p_fac = factor(dens.p_fac,
                        levels = c("[0;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.1[", "[0.1;1]")))

mypal2 <- mf_get_pal(n = 5, palette = "OrRd")
mf_map(x = data_immo, 
       var = "dens.p_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       leg_box_border = "grey3",
       pal=mypal2, 
       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 densité de population")
mf_credits("Sources: Notaires de France 2018, INSEE 2019, IGN Admin Express 2021")

Ces cartes des p-values sont particulièrement importantes car elles nous donnent les endroits où l’effet est significatif. En effet, on sait que la VI a un effet global significatif, qu’elle a en plus une variabilité locale or localement elle n’est pas partout significative. Pour la part d’agriculteurs dans l’emploi, l’effet est significatif quasiment uniquement dans le Sud-Est.