Proyecciones de población Econometría Espacial
Francisco Parra
27 de marzo de 2019
Planteamiento del problema
La población municipal obtenida de los padrones municipales (población observada), se va a estimar en cada municipio con la población municipal que hubiera resultado de su dinámica natural (Poblacion esperada):
\[Población_{t-1}+Nacimientos_t-Defunciones_t\]
Esta población resultado de la dinámica natural de la población de cada municipio es la que se va a proyectar mediante las tasas demográficas (tasa de natalidad por edades y tasas de moratalidad por edades), se trataría por tanto de una población esperada que como resultado de los movimientos residenciales va a diferir en cada municipio de la población real o empadronada. Las diferencias entre ambas poblaciones en los años 2005-2016, tendrán logicamente un correlación espacial, derivada de los movimientos de población entre municipios limitrofes, asociados a los procesos de urbanización y descongestión de los centros urbanos. Estas correlaciones espaciales, estudiadas con metodología de econometría espacial, servirán para pronosticar los saldos migratorios asociados a las proyecciones que se hagan de la población natural en cada municipio.
Leemos los datos:
setwd("~/investigaciones ICANE/INLA proyecciones")
esperados=read.csv("ESPERADOS.csv",header=TRUE,sep=";")
observados=read.csv("OBSERVADOS.csv",header=TRUE,sep=";")Seleccionamos el estadistico que vamos a analizar: promedio del periodo 2005-20015. Para simplificar la notación, asignamos los valores observados y esperados a los objetos O y E respectivamente.
O <- rowMeans(observados[3:13])
E <- rowMeans(esperados[3:13])Activamos la librerias necesarias, leemos y representamos el mapa municipal de Cantabria
library(maptools)## Loading required package: sp## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()library(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: sf## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3library(RColorBrewer)
library(Hmisc)## Loading required package: lattice## Loading required package: survival## Loading required package: Formula## Loading required package: ggplot2##
## Attaching package: 'Hmisc'## The following object is masked from 'package:maptools':
##
## label## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(sp)
library(leaflet)
#library(gplots)
Carto<-readShapePoly("cant_mun_poly_proyecciones.shp")
# Representamos mapa de Cantabria
plot(Carto)
#data.frame(Carto$COD_INE,Carto$NOMBRE)Finalmente, creamos un objeto “data.frame” con todos los datos necesarios, para ejecutar los diferentes modelos:
Datos <- data.frame(Carto$COD_INE,Carto$NOMBRE, O=O, E=E)
head(Datos)## Carto.COD_INE Carto.NOMBRE O E
## 1 39033 Herrerías 670.0000 671.4545
## 2 39095 Val de San Vicente 2757.9091 2738.1818
## 3 39091 Valdáliga 2349.4545 2347.6364
## 4 39090 Udías 854.3636 850.0909
## 5 39024 Comillas 2427.9091 2405.2727
## 6 39012 Cabezón de la Sal 8199.4545 8151.3636El desarrollo del ejercicio esta basado en Sarmiento I (2017).
Econometría Espacial.
Moran (1948) y Geary (1954) son los pioneros en análisis espacial que se populariza con el trabajo de A. D. Cliff and Ord (1973).
•La autocorrelación espacial extiende la idea de autocorrelación de series de tiempo a una estructura geográfica topológica más general, de donde deviene el adjetivo espacial
•Los valores de un mapa pueden ser vistos como una muestra de tamaño 1 con medidas repetidas y cercanamente correlacionadas.
•Estas medidas son tratadas como observaciones, tal que \(N\) es el número de posiciones en un mapa.
•La estructura geográfica topológica general se puede operativizar con una tabla \(N × N\) de variables indicadoras binarias 0-1, para cada posición en un mapa.
•Las celdas son igual a 1 si la fila y la columna de la tabla se definen como vecinos, y son 0 en otro caso.
Dos áreas serán vecinas si dos áreas de una partición comparten un borde de longitud no nula.
Modelo de mínimos cuadrados ordinarios (MCO)
Realizamos una estimación por MCO de una relación lineal simple entre la población observada y esperada de cada municipio:
saldo.ols=lm(O~E,datos=Datos)## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
## extra argument 'datos' will be disregardedsummary(saldo.ols)##
## Call:
## lm(formula = O ~ E, datos = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -244.50 -43.98 -36.43 -5.00 705.93
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.405e+01 1.208e+01 3.645 0.000426 ***
## E 9.989e-01 6.117e-04 1633.055 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 116.9 on 100 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.667e+06 on 1 and 100 DF, p-value: < 2.2e-16plot(saldo.ols)
Análisis gráfico de los residuos del modelo MCO
Carto$saldo.ols.res<-resid(saldo.ols) #residuals ols
spplot(Carto,"saldo.ols.res", at=seq(min(Carto$saldo.ols.res,na.rm=TRUE),max(Carto$saldo.ols.res,na.rm=TRUE),length=12),col.regions=rev(brewer.pal(11,"RdBu")))
Modelos de Dependencia Espacial
La autocorrelación espacial mide el grado en que un fenómeno de interés se correlaciona consigo mismo en el espacio (Cliff y Ord (1973)). En otras palabras, cuando valores similares aparecen próximos entre sí en el espacio, tendremos, autocorrelación espacial positiva, y cuando los valores vecinos sean diferentes, tendremos autocorrelación espacial negativa. La autocorrelación espacial nula indica que el patrón espacial es aleatorio. Siguiendo a Anselin y Bera (1998) podemos expresar la existencia de autocorrelación espacial con la siguiente condición de momento:
\[\begin{equation} Cov(y_{i},y_{j})\neq 0\,\,\,\, for\,\,\,\,i\neq j \end{equation}\]
donde \(y_i\) y \(y_j\) son observaciones en una variable aleatoria en las ubicaciones i y j.
El problema es que necesitamos estimar los términos de covarianza de \(N x N\) directamente para las observaciones de \(N\). Para superar este problema imponemos restricciones sobre la naturaleza de las interacciones. Un tipo de restricción es definir para cada punto de datos un “conjunto de vecindad” relevante. En la econometría espacial, esto se hace operativo a través de la matriz de ponderaciones espaciales. La matriz generalmente denotada por W es una matriz \(N x N\) positiva y simétrica que denota para cada observación (fila) aquellas ubicaciones (columnas) que pertenecen a su conjunto de vecindad como elementos distintos de cero (Anselin y Bera;1998, y Arbia ; 2014).
La especificación del conjunto vecino es bastante arbitraria y hay una amplia gama de sugerencias en la literatura. Los criterios más usados son los siguientes:
• Criterio de torre (rook criterium): dos unidades están cerca una de la otra si comparten un lado.
knitr::include_graphics("rook-1.png")
• Criterio de reina (queen criterium) : dos unidades están cerca si comparten un lado o un borde.
knitr::include_graphics("queen-1.png")
Otro enfoque utilizado es denotar dos observaciones como vecinas aquellas areas que están dentro de una cierta distancia.
Aquí nos centraremos en utilizar el criterio de reina. La función poly2nb de la librería maptools, construye una lista de vecinos, si se especifica la opción queen = TRUE, la construirá utilizando el criterio de queen. Si se especifica FALSE, utilizará el criterio de torre. La opción W fila estandariza la matriz.
vecinos<-poly2nb(Carto, queen=TRUE)
W<-nb2listw(vecinos, style="W", zero.policy=TRUE)
plot(W,coordinates(Carto))
Coeficiente de Moran (MC)
El coeficiente de Moran o I de Moran es una medida de autocorrelación espacial desarrollada por Moran (1948), se formula;
\[I=\frac{N}{\sum_i \sum_j w_{ij}} \frac{\sum_i \sum_j w_{ij} (X_i-\bar X) (X_j-\bar X) }{\sum_i (X_i-\bar X)^2}\] donde \(N\) es el número de unidades espaciales indexadas por i y j, \(X\) es la variable de interés y \(w_{ij}\) es la matriz de pesos espaciales.
Los valores negativos (positivos) indican negativo (positivo) de autocorrelación espacial. Los valores oscilan entre -1 (indicando dispersión perfecta) a 1 (correlación perfecta). Un valor de cero indica un patrón espacial aleatorio. Para las pruebas de hipótesis estadísticas, los valores de la I pueden ser transformados a la Z-score (distribución normal estandarizada) en el que los valores superiores a 1,96 o menor que -1.96 indican autocorrelación espacial a un nivel de significación del 5% \(\alpha=0,05\).
La función R que realiza la prueba I de Moran en lm.morantes de la libreria spdep, la función requiere un objeto de regresión (modelo lm) y un objeto de lista espacial (matriz de vecindades) como argumentos:
moran.lm<-lm.morantest(saldo.ols, W,zero.policy=TRUE)
print(moran.lm)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = O ~ E, datos = Datos)
## weights: W
##
## Moran I statistic standard deviate = 3.203, p-value = 0.00068
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.19589384 -0.01142751 0.00418959Lagrange Multiplier Test
Aunque la prueba I de Moran tiene un gran poder frente a una amplia gama de alternativas (Anselin y Bera ,1998), no guía en la selección de modelos alternativos, ya que dicho contraste no presenta una hipótesis alternativa. Basado en el principio de máxima verosimilitud, la prueba del multiplicador de Lagrange especifica la hipótesis alternativa que nos ayudará con la tarea. Las pruebas LM para la dependencia espacial se incluyen en la función lm.LMtests e incluyen como alternativa la presencia de un retraso espacial y/o la presencia de un retraso espacial en el término de error. Ambas pruebas, así como sus formas robustas se incluyen en la función lm.LMtestsde la libreria spdep . Para ejecutar ambos utilizamos la opción test = “all”. De nuevo, un objeto de regresión y un objeto de lista espacial deben indicarse como argumentos:
LM<-lm.LMtests(saldo.ols, W, test="all",zero.policy=TRUE)
print(LM)##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = O ~ E, datos = Datos)
## weights: W
##
## LMerr = 8.7761, df = 1, p-value = 0.003052
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = O ~ E, datos = Datos)
## weights: W
##
## LMlag = 16.263, df = 1, p-value = 5.514e-05
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = O ~ E, datos = Datos)
## weights: W
##
## RLMerr = 8.5706, df = 1, p-value = 0.003416
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = O ~ E, datos = Datos)
## weights: W
##
## RLMlag = 16.057, df = 1, p-value = 6.145e-05
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = O ~ E, datos = Datos)
## weights: W
##
## SARMA = 24.833, df = 2, p-value = 4.05e-06Modelo de Retraso Autoregresivo Espacial (SAR)
La estimación del modelo SAR se puede abordar asumiemdo la normalidad del término de error y estumando por máxima verosimilitud el siguiente modelo matricial:
\[y=\rho Wy + X\beta + u\]
donde \(y\) es el vector de la dependiente, \(W\) la matriz de vecindades, \(X\) la matriz de de explicativas, e \(u\) un vector de errorese aleatorios normalmente distribuidos.
La función lagsarlm de la librería spdep proporciona la estimación de máxima verosimilitud del modelo SAR.
sar.saldo<-lagsarlm(O~E,data=Datos, W,zero.policy=TRUE)
summary(sar.saldo)##
## Call:
## lagsarlm(formula = O ~ E, data = Datos, listw = W, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -196.0456 -32.6484 -25.0818 -4.3526 666.3916
##
## Type: lag
## Regions with no neighbours included:
## 85
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.2388e+01 1.2009e+01 1.8642 0.06229
## E 9.9813e-01 5.8078e-04 1718.5876 < 2e-16
##
## Rho: 0.0051395, LR test value: 17.532, p-value: 2.8258e-05
## Approximate (numerical Hessian) standard error: 0.0011744
## z-value: 4.3762, p-value: 1.2074e-05
## Wald statistic: 19.152, p-value: 1.2074e-05
##
## Log likelihood: -620.6397 for lag model
## ML residual variance (sigma squared): 11288, (sigma: 106.25)
## Number of observations: 102
## Number of parameters estimated: 4
## AIC: 1249.3, (AIC for lm: 1264.8)Análisis gráfico de los residuos del modelo SAR
Carto$saldo.sar.res<-resid(sar.saldo) #residual sar## Warning: Method residuals.sarlm moved to the spatialreg package## Warning in residuals.sarlm(sar.saldo): install the spatialreg packagespplot(Carto,"saldo.sar.res",at=seq(min(Carto$saldo.sar.res,na.rm=TRUE),max(Carto@data$saldo.sar,na.rm=TRUE), length=12), col.regions=rev(brewer.pal(11,"RdBu")))
La presencia de la matriz de ponderaciones espaciales hace que los efectos marginales sean más ricos y ligeramente más complicados que en el modelo OLS “tradicional”.
Modelo SEM
Estimación de un modelo de errores autorregresivos espaciales simultáneos (SEM), se formula:
\[y = X\beta + u, u = \lambda W u + e\]
Si queremos estimar el modelo de error espacial, tenemos dos enfoques. Primero, podemos usar la Máxima versosimilitud como antes, con la función errorsarlm de la libreria spdep.
errorsalm.saldo<-errorsarlm(O~E, data=Datos, W,zero.policy=TRUE)
summary(errorsalm.saldo)##
## Call:
## errorsarlm(formula = O ~ E, data = Datos, listw = W, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -171.4564 -34.9492 -27.0335 -6.7702 679.2519
##
## Type: error
## Regions with no neighbours included:
## 85
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.6846e+01 1.8003e+01 2.602 0.009267
## E 9.9778e-01 5.9126e-04 1687.568 < 2.2e-16
##
## Lambda: 0.39486, LR test value: 9.0547, p-value: 0.0026203
## Approximate (numerical Hessian) standard error: 0.11514
## z-value: 3.4293, p-value: 0.00060505
## Wald statistic: 11.76, p-value: 0.00060505
##
## Log likelihood: -624.8781 for error model
## ML residual variance (sigma squared): 11834, (sigma: 108.79)
## Number of observations: 102
## Number of parameters estimated: 4
## AIC: 1257.8, (AIC for lm: 1264.8)Gráfica para los residuos de SEM con errorsarlm:
Carto$saldo.sem.res<-resid(errorsalm.saldo) #residual sem## Warning: Method residuals.sarlm moved to the spatialreg package## Warning in residuals.sarlm(errorsalm.saldo): install the spatialreg packagespplot(Carto,"saldo.sem.res",at=seq(min(Carto$saldo.sem.res,na.rm=TRUE),max(Carto@data$saldo.sem.res,na.rm=TRUE), length=12), col.regions=rev(brewer.pal(11,"RdBu")))
Un segundo enfoque es utilizar mínimos cuadrados generalizados factibles (GLS) con la función GMerrorsarde la libreria spdep .
fgls.saldo<-GMerrorsar(O~E, data=Datos, W,zero.policy=TRUE)## Warning: Function GMerrorsar moved to the spatialreg package## Warning in GMerrorsar(O ~ E, data = Datos, W, zero.policy = TRUE): install
## the spatialreg package## Warning: Function as_dgRMatrix_listw moved to the spatialreg package## Warning in as_dgRMatrix_listw(from): install the spatialreg package## Warning: Function as_dsCMatrix_I moved to the spatialreg package## Warning in as_dsCMatrix_I(n): install the spatialreg packagesummary(fgls.saldo)## Warning: Method summary.gmsar moved to the spatialreg package## Warning in summary.gmsar(fgls.saldo): install the spatialreg package## Warning: Method print.summary.gmsar moved to the spatialreg package## Warning in print.summary.gmsar(x): install the spatialreg package##
## Call:
## GMerrorsar(formula = O ~ E, data = Datos, listw = W, zero.policy = TRUE)
##
## Residuals:## Warning: Method residuals.gmsar moved to the spatialreg package## Warning in residuals.gmsar(x): install the spatialreg package## Min 1Q Median 3Q Max
## -71.912 -46.018 -36.635 -1.196 721.726
##
## Type: GM SAR estimator
## Regions with no neighbours included:
## 85
## Coefficients: (GM standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.6819e+01 1.7101e+01 2.7377 0.006186
## E 9.9789e-01 6.0055e-04 1661.6282 < 2.2e-16
##
## Lambda: 0.34849 (standard error): 0.25721 (z-value): 1.3549
## Residual variance (sigma squared): 12272, (sigma: 110.78)
## GM argmin sigma squared: 12404
## Number of observations: 102
## Number of parameters estimated: 4Gráfica para los residuos de SEM con GMerrorsar:
Carto$saldo.fgls.res<-resid(fgls.saldo) #residual fgls## Warning: Method residuals.gmsar moved to the spatialreg package## Warning in residuals.gmsar(fgls.saldo): install the spatialreg packagespplot(Carto,"saldo.fgls.res",at=seq(min(Carto$saldo.fgls.res,na.rm=TRUE),max(Carto@data$saldo.fgls,na.rm=TRUE), length=12), col.regions=rev(brewer.pal(11,"RdBu")))
Pronósticos
Se pueden realizar pronosticos con el modelo SEM (errorsarlm) y la función predict
newdata=data.frame(E=esperados[,14])
E.2016 <- predict(errorsalm.saldo, newdata, listw=lw)
head(E.2016)## [1] 689.4188 2857.6040 2329.7761 925.8936 2332.7695 8385.3289newdata=data.frame(E=esperados[,15])
E.2017 <- predict(errorsalm.saldo, newdata, listw=lw)
head(E.2017)## [1] 696.4033 2841.6395 2326.7828 938.8648 2280.8847 8373.3555newdata=data.frame(E=esperados[,16])
E.2018 <- predict(errorsalm.saldo, newdata, listw=lw)
head(E.2018)## [1] 681.4365 2792.7480 2302.8360 947.8449 2249.9534 8348.4109Datos=cbind(Datos,Est_2016=as.data.frame(E.2016)$fit,Est_2017=as.data.frame(E.2017)$fit,Est_2018=as.data.frame(E.2018)$fit)
head(Datos)## Carto.COD_INE Carto.NOMBRE O E Est_2016 Est_2017
## 1 39033 Herrerías 670.0000 671.4545 689.4188 696.4033
## 2 39095 Val de San Vicente 2757.9091 2738.1818 2857.6040 2841.6395
## 3 39091 Valdáliga 2349.4545 2347.6364 2329.7761 2326.7828
## 4 39090 Udías 854.3636 850.0909 925.8936 938.8648
## 5 39024 Comillas 2427.9091 2405.2727 2332.7695 2280.8847
## 6 39012 Cabezón de la Sal 8199.4545 8151.3636 8385.3289 8373.3555
## Est_2018
## 1 681.4365
## 2 2792.7480
## 3 2302.8360
## 4 947.8449
## 5 2249.9534
## 6 8348.4109write.csv(Datos,file="Datos.csv")Evaluamos los errores medios cuadráticos de cada pronostico
error.2016=sqrt(sum((observados[,14]-E.2016)^2)/105)
error.2016## [1] 100.212error.2017=sqrt(sum((observados[,15]-E.2017)^2)/105)
error.2017## [1] 65.72346error.2018=sqrt(sum((observados[,16]-E.2018)^2)/105)
error.2018## [1] 78.48909Bibliografía
Anselin, Luc. 2003. “An Introduction to Spatial Regression Analysis in R.” Available at: Https://Geodacenter.asu.edu/Drupal_files/Spdepintro.pdf.
Anselin, Luc, and Anil K Bera. 1998. “Spatial Dependence in Linear Regression Models with an Introduction to Spatial Econometrics.” Statistics Textbooks and Monographs 155. MARCEL DEKKER AG: 237–90.
Arbia, Giuseppe. 2014. A Primer for Spatial Econometrics: With Applications in R. Palgrave Macmillan.
Bivand, Roger, and Nicholas Lewin-Koh. 2016. Maptools: Tools for Reading and Handling Spatial Objects. https://CRAN.R-project.org/package=maptools.
Cheng, Joe, and Yihui Xie. 2015. Leaflet: Create Interactive Web Maps with the Javascript ’Leaflet’ Library. http://rstudio.github.io/leaflet/.
Cliff, Andrew David, and J Keith Ord. 1973. Spatial Autocorrelation. Vol. 5. Pion London.
ESRI, Environmental Systems Research Institute. 1998. “ESRI Shapefile Technical Description.” Available at: Https://Www.esri.com/Library/Whitepapers/Pdfs/Shapefile.pdf.
Geary, R.C. (1954) The Contiguity Ratio and Statistical Mapping. The Icorporporated Statistician, 5, 115-145. https://doi.org/10.2307/2986645.
Moran, P. 1948. The Interpretation of Statistical Maps. Journal of the Royal Statistical Society, 10, 243-251.
Neuwirth, Erich. 2014. RColorBrewer: ColorBrewer Palettes. https://CRAN.R-project.org/package=RColorBrewer.
Sarmiento, I: An Introduction to Spatial Econometrics in R.https://www.r-bloggers.com/an-introduction-to-spatial-econometrics-in-r/