MODELISATION STATISTIQUE SPATIALE DES PRIX MEDIAN DE LOGEMENT A BOSTON
Espoir Zyad KOSSOU
2023-11-13
Introduction
L’analyse spatiale et la modélisation statistique spatiale jouent un rôle crucial dans la compréhension des processus géographiques et dans la prise de décision éclairée dans divers domaines tels que la géographie, l’économie, l’environnement, la santé, etc. Ces méthodes permettent de prendre en compte la dépendance spatiale entre les observations, reconnaissant que les valeurs dans une région peuvent influencer ou être influencées par les valeurs dans les régions voisines.
Parmi les approches couramment utilisées en modélisation statistique spatiale, on trouve les modèles SAR (Spatial Autoregressive), les modèles SEM (Spatial Error), ainsi que d’autres techniques qui intègrent la dimension spatiale dans l’analyse statistique.
Explication des modèles spatiaux
ELHORST 2010 a établi une classification des principaux modèles d’économétrie spatiale, en s’appuyant sur les trois types d’interaction spatiale issus du modèle fondateur de MANSKI 1993 :
— une interaction endogène, lorsque la décision économique d’un agent ou d’une zone géographique va dépendre de la décision de ses voisins;
— une interaction exogène, lorsque la décision économique d’un agent va dépendre des caractéristiques observables de ses voisins;
— une corrélation spatiale des effets liée à de mêmes caractéristiques inobservées.
Ce modèle s’écrit sous forme matricielle :
\[ Y = \rho \cdot WY + X \cdot \beta + W X \cdot \theta + u \] \[ u = \lambda \cdot Wu + \varepsilon \] Avec les paramètres \(\beta\) pour les variables explicatives exogènes, \(\rho\) pour l’effet d’interaction endogène (de dimension 1) dit autorégressif spatial, \(\theta\) pour les effets d’interaction exogène (de dimension égale au nombre de variables exogènes K) et \(\lambda\) pour l’effet de corrélation spatiale des erreurs dit autocorrélation spatiale.
Matériel et Méthode de travail
Matériel
Jeu de données
Le jeu de données “Boston” (aussi appelé “Boston.c”) en R fait référence au “Boston Housing Dataset”. C’est un ensemble de données classique souvent utilisé pour des analyses statistiques, notamment dans le domaine de l’apprentissage automatique et de la régression. Il contient des informations sur le prix médian des maisons dans différentes zones de Boston, ainsi que diverses caractéristiques socio-économiques associées à ces zones.
Ce jeu de données est composé de 506 observations (zones géographiques) et 14 variables. Voici une description des principales variables :
medv : C’est la variable dépendante, représentant le prix médian des logements en milliers de dollars (exprimé en $1000).
crim : Taux de criminalité par habitant.
zn : Proportion de terrains résidentiels pour des lots de plus de 25 000 pieds carrés.
indus : Proportion de superficies commerciales/industrielles par ville.
chas : Variable binaire, égale à 1 si la zone borde le fleuve Charles et 0 sinon.
nox : Concentration d’oxyde nitrique (parties par 10 millions).
rm : Nombre moyen de pièces par logement.
age : Proportion de logements occupés par leur propriétaire construits avant 1940.
dis : Distance pondérée aux cinq centres d’emploi de Boston.
rad : Indice d’accessibilité aux autoroutes radiales.
tax : Taux d’imposition foncière pour 10 000 dollars de valeur.
ptratio : Ratio élèves/enseignants dans les écoles publiques.
black : Proportion de la population afro-américaine par ville.
lstat : Pourcentage de la population à faible statut socio-économique.
Méthode de travail
L’objectif est d’expliquer le prix médian des logements en fonction du Nombre moyen de pièces par logement, du ratio élèves/enseignants dans les écoles publiques, de la Ddistance pondérée aux cinq centres d’emploi de Boston, du taux d’imposition foncière. Pour modéliser nos données spatiales, nous allons faire une regression linéaire. Une matrice de poids spatiale basée sur la contiguïté sera ensuite construite afin de représenter la relation de voisinage entre les observations spatiales. L’étape suivante sera d’éffectuer un test de Moran afin de voir s’il y a autocorrélation spatiale. Si nous la détectons, les modèles spatiales autoregressifs (SAR) et d’autocorrelation spatiale des erreurs (SEM) seront implémentés et le meilleur sera choisi à l’aide du test multiplicatif de Lagrange
La matrice de poids spatiale basée sur la contiguïté \(W\) permet de représenter la relation de voisinage entre les observations spatiales.
Résultats et discussion
Chargement des packages
library(MASS)
library(spdep)
library(spatialreg)
library(spData)
library(stargazer)Chargement du jeu de données
data("Boston")
data = boston.c
attach(data)
Y = cbind(MEDV)
X = cbind(RM, DIS, TAX, PTRATIO)
xy <- cbind(data$X, data$Y)Visualisation du jeu de données
neighbors <- boston.soi
coords <- boston.utm
# Descriptive statistics
summary(Y)## MEDV
## Min. : 5.00
## 1st Qu.:17.02
## Median :21.20
## Mean :22.53
## 3rd Qu.:25.00
## Max. :50.00summary(X)## RM DIS TAX PTRATIO
## Min. :3.561 Min. : 1.130 Min. :187.0 Min. :12.60
## 1st Qu.:5.886 1st Qu.: 2.100 1st Qu.:279.0 1st Qu.:17.40
## Median :6.208 Median : 3.207 Median :330.0 Median :19.05
## Mean :6.285 Mean : 3.795 Mean :408.2 Mean :18.46
## 3rd Qu.:6.623 3rd Qu.: 5.188 3rd Qu.:666.0 3rd Qu.:20.20
## Max. :8.780 Max. :12.127 Max. :711.0 Max. :22.00# Neighbors summary
plot(neighbors, coords)
Modèle de regression linéaire
On considère le modèle de régression linéaire multiple :
\[ y = X\beta + \varepsilon \]
où \(y\) est le vecteur de la variable dépendante, \(X\) la matrice des variables explicatives (\(n\) lignes et \(p\) colonnes), \(\beta\) le vecteur des coefficients de régression et \(\varepsilon\) le vecteur des erreurs.
olsreg = lm(MEDV ~ RM + DIS + TAX + PTRATIO, data = data)
summary(olsreg)##
## Call:
## lm(formula = MEDV ~ RM + DIS + TAX + PTRATIO, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.506 -3.152 -0.455 2.223 41.493
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.360075 4.092513 -0.332 0.740
## RM 7.345490 0.403958 18.184 < 2e-16 ***
## DIS -0.153606 0.147290 -1.043 0.298
## TAX -0.012305 0.002023 -6.084 2.34e-09 ***
## PTRATIO -0.902955 0.141046 -6.402 3.53e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.876 on 501 degrees of freedom
## Multiple R-squared: 0.595, Adjusted R-squared: 0.5918
## F-statistic: 184 on 4 and 501 DF, p-value: < 2.2e-16Interpretation
Le R2 ajusté est de 0.5918, cela signifie que le modèle explique environ 59% de la variance de la variable dépendante MEDV.
Le test F est significatif (p-value < 2.2e-16), le modèle est globalement significatif.
Parmi les variables explicatives :
RM (nombre de pièces) : effet positif significatif sur MEDV, une pièce supplémentaire augmente MEDV de 7.345 toute chose étant égale par ailleurs.
DIS (distance aux centres d’emploi) : effet négatif non significatif sur MEDV
TAX (taux de taxes) : effet négatif significatif, une augmentation du taux de taxes réduit MEDV de 0.012 toute chose étant égale par ailleurs.
PTRATIO (ratio élèves/enseignants) : effet négatif significatif, une augmentation du ratio réduit MEDV de 0.903 toute chose étant égale par ailleurs. L’erreur résiduelle standard est de 5.876, elle représente la variance de MEDV non expliquée par le modèle.
Matrice de Poids basée sur la contiguité
La matrice de poids spatiale basée sur la contiguïté \(W\) permet de représenter la relation de voisinage entre les observations spatiales. Ses éléments \(w_{ij}\) sont définis comme:
\[ w_{ij} = \begin{cases} 1, & \text{si $i$ et $j$ sont contigus} \\ 0, & \text{sinon} \end{cases}\]
Deux observations \(i\) et \(j\) sont contigües si elles partagent une frontière ou un sommet.
Cette matrice décrit donc la topologie du voisinage spatial sous la forme d’une matrice carrée de dimension \(N\), avec \(N\) le nombre d’observations.
Elle est généralement normalisée par ligne pour que les poids de chaque ligne somment à 1.
weights = nb2listw(neighbors)
summary(weights)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 506
## Number of nonzero links: 2152
## Percentage nonzero weights: 0.8405068
## Average number of links: 4.252964
## Link number distribution:
##
## 1 2 3 4 5 6 7 8
## 10 39 93 143 143 52 22 4
## 10 least connected regions:
## 2033 2081 2181 3592 0001 0101 0108 0201 1301 1805 with 1 link
## 4 most connected regions:
## 3564 3744 4196 4223 with 8 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 506 256036 506 261.3891 2071.255Interpretation
Elle contient 506 régions (observations)
Il y a 2152 liens non nuls, c’est à dire 2152 paires de régions considérées comme voisines spatialement
84.05% des poids sont non nuls, cela représente la proportion de paires de régions voisines
En moyenne, chaque région a 4.25 voisins
La distribution du nombre de liens montre que :
10 régions ont seulement 1 voisin
39 régions ont 2 voisins
93 régions ont 3 voisins
etc.
Les 4 régions les plus connectées ont 8 voisins chacune
En résumé, cette sortie donne des informations sur la topologie du graphe spatial en analysant la matrice de poids. Elle permettra de mieux appréhender la structure de dépendance spatiale pour les analyses.
Test de Moran
Le test de Moran permet de tester la présence d’une autocorrélation spatiale.
Soit \(z_i\) la variable d’intérêt en site \(i\) et \(w_{ij}\) les poids spatiaux entre sites \(i\) et \(j\).
La statistique du I de Moran est définie comme:
\[I = \frac{n}{\sum_i \sum_j w_{ij}} \frac{\sum_i \sum_j w_{ij} (z_i - \bar{z})(z_j - \bar{z})}{\sum_i (z_i - \bar{z})^2}\]
Sous l’hypothèse nulle H0 d’absence d’autocorrélation spatiale, on a:
\[E[I] = - \frac{1}{n-1}\] \[Var[I] = \frac{n}{S_0^2} \left( \frac{n-1}{n} \right)\]
où \(S_0 = \sum_i \sum_j w_{ij}\).
La statistique de test est \(I/\sqrt{Var[I]}\) qui suit une loi normale sous H0.
moran.test(Y, weights)##
## Moran I test under randomisation
##
## data: Y
## weights: weights
##
## Moran I statistic standard deviate = 21.599, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.684367739 -0.001980198 0.001009728Interpretation Le test est significatif avec une p-value < 2.2e-16. L’hypothèse nulle d’absence d’autocorrélation spatiale est donc rejetée La statistique I de Moran observée vaut 0.684367739. Il détecte une autocorrélation spatiale significative pour cette variable et cette matrice de poids spatiale. Un modèle prenant en compte cette dépendance spatiale sera nécessaire.
Test multiplicatif de Lagrange
Le test du multiplicateur de Lagrange (LM) permet de tester la présence d’une autocorrélation spatiale des résidus.
On considère un modèle de régression linéaire :
\[y = X\beta + u, \quad u = \lambda W u + \varepsilon\]
Sous l’hypothèse nulle H0 de non autocorrélation spatiale (= 0), la statistique du test LMerr vaut :
\[LM_{err} = \frac{1}{\sigma^2} (e'W'We) \sim \chi^2(1)\]
où \(e\) sont les résidus du modèle OLS.
De même, sous H0, la statistique du test LMlag vaut:
\[LM_{lag} = \frac{1}{\sigma^2} (W y)'(In - X(X'X)^{-1}X')(W y) \sim \chi^2(1)\]
Un rejet de H0 indique la présence d’autocorrélation spatiale des résidus.
lm.LMtests(olsreg, weights, test=c("LMlag", "LMerr"))##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ RM + DIS + TAX + PTRATIO, data = data)
## weights: weights
##
## LMlag = 260.03, df = 1, p-value < 2.2e-16
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ RM + DIS + TAX + PTRATIO, data = data)
## weights: weights
##
## LMerr = 326.17, df = 1, p-value < 2.2e-16Interpretation
Il s’agit des tests du multiplicateur de Lagrange (LM) pour détecter l’autocorrélation spatiale des résidus d’un modèle de régression linéaire (OLS) :
Le test LMlag teste l’autocorrélation de type “spatial lag”, c’est-à-dire sur la variable dépendante Le test LMerr teste l’autocorrélation de type “spatial error”, c’est-à-dire sur le terme d’erreur. Les deux tests sont très significatifs avec une p-value < 2.2e-16. On rejette donc l’hypothèse nulle d’absence d’autocorrélation spatiale des résidus Le test LMerr étant légèrement plus significatif, l’autocorrélation semble plus importante sur le terme d’erreur que sur la variable dépendante En résumé, ces tests confirment la présence d’autocorrélation spatiale des résidus, en particulier de type “spatial error”.
Un modèle de régression spatiale de type “spatial error” serait donc adapté pour prendre en compte cette dépendance spatiale résiduelle.
Modèle SAR
Le modèle autorégressif spatial (SAR) permet de prendre en compte une autocorrélation spatiale de la variable dépendante.
Le modèle s’écrit :
\[y = \rho W y + X \beta + \varepsilon\]
où \(W\) est la matrice de poids spatiaux, \(\rho\) est le coefficient d’autocorrélation spatiale et \(\beta\) les coefficients de régression.
On peut réécrire ce modèle sous la forme réduite:
\[y = (I_n - \rho W)^{-1} X \beta + (I_n - \rho W)^{-1} \varepsilon\]
L’estimateur du maximum de vraisemblance des paramètres \(\beta\) et \(\rho\) est :
\[\hat{\beta} = (X'X)^{-1}X'(I_n - \rho W)y\]
\[\hat{\rho} = \text{solution de l'équation du maximum de vraisemblance}\]
lag <- lagsarlm(MEDV ~ RM + DIS + TAX + PTRATIO, data = data, listw=weights)
summary(lag)##
## Call:lagsarlm(formula = MEDV ~ RM + DIS + TAX + PTRATIO, data = data,
## listw = weights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.55109 -2.28838 -0.38003 1.71112 31.90751
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.5525332 3.0185808 -4.1584 3.205e-05
## RM 5.0942407 0.3270521 15.5762 < 2.2e-16
## DIS -0.2868204 0.1074350 -2.6697 0.0075917
## TAX -0.0061751 0.0015229 -4.0549 5.015e-05
## PTRATIO -0.3620937 0.1065590 -3.3981 0.0006787
##
## Rho: 0.59089, LR test value: 260.53, p-value: < 2.22e-16
## Asymptotic standard error: 0.029891
## z-value: 19.768, p-value: < 2.22e-16
## Wald statistic: 390.79, p-value: < 2.22e-16
##
## Log likelihood: -1481.294 for lag model
## ML residual variance (sigma squared): 18.371, (sigma: 4.2861)
## Number of observations: 506
## Number of parameters estimated: 7
## AIC: 2976.6, (AIC for lm: 3235.1)
## LM test for residual autocorrelation
## test value: 10.603, p-value: 0.0011289Il s’agit de la sortie du modèle de régression spatiale de type “spatial lag” :
Le coefficient Rho est estimé à 0.59089, significatif. Il mesure l’autocorrélation spatiale sur la variable dépendante. Les coefficients des variables explicatives RM, DIS, TAX et PTRATIO sont toujours significatifs, mais leurs valeurs changent par rapport au modèle OLS, notamment l’intercept. Le R2 ajusté n’est pas disponible pour ce type de modèle. Le AIC est plus faible que pour le modèle OLS, indiquant une meilleure qualité d’ajustement. Le test de ratio de vraisemblance (LR) confirme que ce modèle est préférable au modèle OLS. Le test LM des résidus indique qu’il reste une légère autocorrélation spatiale résiduelle. En résumé, ce modèle spatial lag permet de prendre en compte l’autocorrélation détectée précédemment et améliore l’ajustement, mais des améliorations sont encore possibles sur les résidus.
Modèle SEM
Le modèle à erreur spatiale (SEM) permet de prendre en compte une autocorrélation spatiale des erreurs.
Le modèle s’écrit :
\[y = X \beta + u, \quad avec \quad u = \lambda W u + \varepsilon\]
où \(W\) est la matrice de poids spatiaux, \(\lambda\) est le coefficient d’autocorrélation spatiale des erreurs et \(\beta\) les coefficients de régression.
On peut réécrire ce modèle sous la forme réduite :
\[y = X \beta + (I_n - \lambda W)^{-1} \varepsilon\]
L’estimateur du maximum de vraisemblance des paramètres \(\beta\) et \(\lambda\) est :
\[\hat{\beta} = (X'X)^{-1}X'y\]
\[\hat{\lambda} = \text{solution de l'équation du maximum de vraisemblance}\]
err <- errorsarlm(MEDV ~ RM + DIS + TAX + PTRATIO, data = data, listw=weights)
summary(err)##
## Call:errorsarlm(formula = MEDV ~ RM + DIS + TAX + PTRATIO, data = data,
## listw = weights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.08032 -2.33619 -0.46701 1.43650 28.57243
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.3285277 4.1051670 1.5416 0.1232
## RM 5.5794586 0.3416672 16.3301 < 2.2e-16
## DIS -0.1231771 0.2661852 -0.4627 0.6435
## TAX -0.0126354 0.0026203 -4.8222 1.420e-06
## PTRATIO -0.7242637 0.1625443 -4.4558 8.358e-06
##
## Lambda: 0.69027, LR test value: 267.39, p-value: < 2.22e-16
## Asymptotic standard error: 0.033514
## z-value: 20.596, p-value: < 2.22e-16
## Wald statistic: 424.21, p-value: < 2.22e-16
##
## Log likelihood: -1477.865 for error model
## ML residual variance (sigma squared): 17.247, (sigma: 4.153)
## Number of observations: 506
## Number of parameters estimated: 7
## AIC: 2969.7, (AIC for lm: 3235.1)Le coefficient Lambda est estimé à 0.69027, significatif. Il mesure l’autocorrélation spatiale sur le terme d’erreur. Les coefficients sont différents de ceux du modèle OLS. Seuls RM, TAX et PTRATIO sont significatifs. Le R2 ajusté est absent, remplacé par le log-likelihood. Le AIC est inférieur au modèle OLS, indiquant un meilleur ajustement. Le test de ratio de vraisemblance (LR) confirme que ce modèle spatial error est préférable au modèle OLS. Il n’y a pas de test sur les résidus car ce modèle prend en compte l’autocorrélation spatiale des erreurs. En résumé, ce modèle semble plus adapté que le modèle lag pour capturer l’autocorrélation spatiale des erreurs mise en évidence précédemment. Il améliore sensiblement l’ajustement.
Tableau d’illustration de la comparaison entre les modèles
stargazer(olsreg, lag, err, type = "text")##
## ======================================================================
## Dependent variable:
## --------------------------------------------------
## MEDV
## OLS spatial spatial
## autoregressive error
## (1) (2) (3)
## ----------------------------------------------------------------------
## RM 7.345*** 5.094*** 5.579***
## (0.404) (0.327) (0.342)
##
## DIS -0.154 -0.287*** -0.123
## (0.147) (0.107) (0.266)
##
## TAX -0.012*** -0.006*** -0.013***
## (0.002) (0.002) (0.003)
##
## PTRATIO -0.903*** -0.362*** -0.724***
## (0.141) (0.107) (0.163)
##
## Constant -1.360 -12.553*** 6.329
## (4.093) (3.019) (4.105)
##
## ----------------------------------------------------------------------
## Observations 506 506 506
## R2 0.595
## Adjusted R2 0.592
## Log Likelihood -1,481.294 -1,477.865
## sigma2 18.371 17.247
## Akaike Inf. Crit. 2,976.587 2,969.730
## Residual Std. Error 5.876 (df = 501)
## F Statistic 184.010*** (df = 4; 501)
## Wald Test (df = 1) 390.788*** 424.207***
## LR Test (df = 1) 260.535*** 267.392***
## ======================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01