[Descargar - RMarkdown]
Regresión espacial
sf
).ggplot2
,
mapview
, tmap
).spatialreg
, antes
spdep
).spatialreg
).dplyr
).Estos ejercicios están basados en ejercicios desarrollados por Luc
Anselin (ver aquí
y aquí
para estudiar los detalles). La especificación de los datos también la
pueden encontrar en la página de GeoDa. Este documento actualiza los
ejericios utilizando spatialreg
en vez de
spdep
.
Leer los paquetes estadísticos:
knitr::opts_chunk$set(message = FALSE, warning=FALSE)
library(sp)
library(sf)
library(spdep)
library(spData)
library(spatialreg)
library(lmtest) # para la prueba Breusch-Pagan de heteroscedasticidad.
library(car) # para la prueba VIG de multicolinearidad
library(tmap)
library(spgwr) # GWR
Datos:
Cargar datos de Columbus:
#Importar datos espaciales y pesos
data(columbus)
Explorar columnas de datos:
columbus$INC
## [1] 19.531 21.232 15.956 4.477 11.252 16.029 8.438 11.337 17.586 13.598
## [11] 7.467 10.048 9.549 9.963 9.873 7.625 9.798 13.185 11.618 31.070
## [21] 10.655 11.709 21.155 14.236 8.461 8.085 10.822 7.856 8.681 13.906
## [31] 16.940 18.942 9.918 14.948 12.814 18.739 17.017 11.107 18.477 29.833
## [41] 22.207 25.873 13.380 16.961 14.135 18.324 18.950 11.813 18.796
Explorar objetos de Columbus:
class(col.gal.nb)
## [1] "nb"
class(columbus)
## [1] "data.frame"
class(polys)
## [1] "polylist"
summary(columbus)
## AREA PERIMETER COLUMBUS. COLUMBUS.I POLYID
## Min. :0.03438 Min. :0.9021 Min. : 2 Min. : 1 Min. : 1
## 1st Qu.:0.09315 1st Qu.:1.4023 1st Qu.:14 1st Qu.:13 1st Qu.:13
## Median :0.17477 Median :1.8410 Median :26 Median :25 Median :25
## Mean :0.18649 Mean :1.8887 Mean :26 Mean :25 Mean :25
## 3rd Qu.:0.24669 3rd Qu.:2.1992 3rd Qu.:38 3rd Qu.:37 3rd Qu.:37
## Max. :0.69926 Max. :5.0775 Max. :50 Max. :49 Max. :49
## NEIG HOVAL INC CRIME
## Min. : 1 Min. :17.90 Min. : 4.477 Min. : 0.1783
## 1st Qu.:13 1st Qu.:25.70 1st Qu.: 9.963 1st Qu.:20.0485
## Median :25 Median :33.50 Median :13.380 Median :34.0008
## Mean :25 Mean :38.44 Mean :14.375 Mean :35.1288
## 3rd Qu.:37 3rd Qu.:43.30 3rd Qu.:18.324 3rd Qu.:48.5855
## Max. :49 Max. :96.40 Max. :31.070 Max. :68.8920
## OPEN PLUMB DISCBD X
## Min. : 0.0000 Min. : 0.1327 Min. :0.370 Min. :24.25
## 1st Qu.: 0.2598 1st Qu.: 0.3323 1st Qu.:1.700 1st Qu.:36.15
## Median : 1.0061 Median : 1.0239 Median :2.670 Median :39.61
## Mean : 2.7709 Mean : 2.3639 Mean :2.852 Mean :39.46
## 3rd Qu.: 3.9364 3rd Qu.: 2.5343 3rd Qu.:3.890 3rd Qu.:43.44
## Max. :24.9981 Max. :18.8111 Max. :5.570 Max. :51.24
## Y AREA NSA NSB
## Min. :24.96 Min. : 1.093 Min. :0.0000 Min. :0.0000
## 1st Qu.:28.26 1st Qu.: 3.193 1st Qu.:0.0000 1st Qu.:0.0000
## Median :31.91 Median : 6.029 Median :0.0000 Median :1.0000
## Mean :32.37 Mean : 6.372 Mean :0.4898 Mean :0.5102
## 3rd Qu.:35.92 3rd Qu.: 7.989 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :44.07 Max. :21.282 Max. :1.0000 Max. :1.0000
## EW CP THOUS NEIGNO
## Min. :0.0000 Min. :0.0000 Min. :1000 Min. :1001
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:1000 1st Qu.:1013
## Median :1.0000 Median :0.0000 Median :1000 Median :1025
## Mean :0.5918 Mean :0.4898 Mean :1000 Mean :1025
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1000 3rd Qu.:1037
## Max. :1.0000 Max. :1.0000 Max. :1000 Max. :1049
## PERIM
## Min. :0.9021
## 1st Qu.:1.4023
## Median :1.8410
## Mean :1.8887
## 3rd Qu.:2.1992
## Max. :5.0775
Cargar datos de Boston:
#Importar datos espaciales y pesos
data(boston)
Explorar datos Boston:
head(boston.utm)
## x y
## 1 338.73 4679.73
## 2 339.23 4683.33
## 3 340.37 4682.80
## 4 341.05 4683.89
## 5 341.56 4684.44
## 6 342.03 4685.09
Regresión Columbus:
lm(CRIME ~ INC + HOVAL, data=columbus)
##
## Call:
## lm(formula = CRIME ~ INC + HOVAL, data = columbus)
##
## Coefficients:
## (Intercept) INC HOVAL
## 68.6190 -1.5973 -0.2739
Observar cómo tenemos que designar la regresión como un objeto para obtener mayor información de la regresión.
columbus.lm <- lm( CRIME ~ INC + HOVAL, data=columbus) # diferencia entre ambos
summary(columbus.lm)
##
## Call:
## lm(formula = CRIME ~ INC + HOVAL, data = columbus)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.418 -6.388 -1.580 9.052 28.649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.6190 4.7355 14.490 < 2e-16 ***
## INC -1.5973 0.3341 -4.780 1.83e-05 ***
## HOVAL -0.2739 0.1032 -2.654 0.0109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.43 on 46 degrees of freedom
## Multiple R-squared: 0.5524, Adjusted R-squared: 0.5329
## F-statistic: 28.39 on 2 and 46 DF, p-value: 9.341e-09
Diagnósticos Columbus:
#Breush-Pagan para heteroscedasticidad (Ho: homoscedasticidad)
bptest(columbus.lm)
##
## studentized Breusch-Pagan test
##
## data: columbus.lm
## BP = 7.2166, df = 2, p-value = 0.0271
#variance inflation factor para beta's
vif(columbus.lm)
## INC HOVAL
## 1.333117 1.333117
Práctica
Ejecuten una regresión básica con en el conjunto de datos de Boston. Especifiquen el valor de la mediana de la vieneda: log(MEDV) o log(CMEDV) (valor medio de la casa) como variable dependiente y usen variables explanatorias utilizando características:
bo1 <- lm(log(MEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) + AGE + log(DIS) + log(RAD) + PTRATIO+B+ log(LSTAT), data=boston.c)
(La regresión hedónica es el uso de un modelo de regresión para estimar la influencia que varios factores tienen en el precio de un bien o, a veces, en la demanda de un bien.)
Tengan en cuenta que para usar las potencias en la fórmula de regresión, deben encerrarlas en un comando I() de identidad. Por ejemplo, como I(NOX^2).
Asegúrense de crear un objeto lm, por ejemplo “bo1”.
Comparen los resultados de la variable dependiente del valor mediano original (MEDV) y el “corregido” (CMEDV).
log(MEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO+B+ log(LSTAT)
Regresión Boston:
Diagnósticos Boston:
Incorporando pesos Columbus:
col.listw <- nb2listw(col.gal.nb)
print(col.gal.nb)
## Neighbour list object:
## Number of regions: 49
## Number of nonzero links: 230
## Percentage nonzero weights: 9.579342
## Average number of links: 4.693878
print(col.listw)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 49
## Number of nonzero links: 230
## Percentage nonzero weights: 9.579342
## Average number of links: 4.693878
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 49 2401 49 23.48489 204.6687
Práctica
Verifiquen la información de contigüidad en el archivo boston.soi. Conviertan el objeto nb en una lista objeto. Comparen la impresión y el resumen de estos objetos.
Pesos Boston:
Para realizar una prueba de Moran sobre los residuos de la regresión
columbus.lm
, pasen tanto el objeto de regresión
(columbus.lm) como la información de contigüidad listw a la función
lm.morantest.
Columbus naive:
col.moran <- lm.morantest(columbus.lm, col.listw)
col.moran
##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: col.listw
##
## Moran I statistic standard deviate = 2.681, p-value = 0.00367
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.212374153 -0.033268284 0.008394853
Un aspecto a tener en cuenta acerca de estos resultados es que la significancia predeterminada se informa para una sola alternativa (“greater”). Para una hipótesis alternativa de dos colas, el significado sería el doble del valor informado. Puede especificar esto con la opción alternativa, como en:
Moran Columbus (2 colas):
col.moran2 <- lm.morantest(columbus.lm,col.listw,alternative="two.sided")
print(col.moran2)
##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: col.listw
##
## Moran I statistic standard deviate = 2.681, p-value = 0.00734
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I Expectation Variance
## 0.212374153 -0.033268284 0.008394853
Hay dos funciones separadas, moran.test, donde la inferencia se basa en un supuesto de normalidad o aleatoriedad, y moran.mc, para una prueba basada en permutaciones.
Cuando randomisation = FALSE, implica que la inferencia está basada en una aproximación de normalidad. Cuando es TRUE, se examina cuando la varianza de I es calculada bajo el supuesto de aleatoriedad.
Moran Columbus forma correcta:
col.e <- resid(columbus.lm)
col.morane <- moran.test(col.e, col.listw, randomisation = FALSE)
print(col.morane)
##
## Moran I test under normality
##
## data: col.e
## weights: col.listw
##
## Moran I statistic standard deviate = 2.4774, p-value = 0.006617
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.212374153 -0.020833333 0.008860962
Autocorrelación residuos MCO columbus:
plot(col.e, lag.listw(col.listw, col.e))
moran.plot(col.e, col.listw)
Estimnando la prueba de Moran con permutaciones, se calcula el rango de la estadística observada relativo a la distribución de referencia de estadísticas para los conjuntos de datos permutados.
El valor semilla de números aleatorios se determina a partir de la hora actual, por lo que ninguna permutación aleatoria será idéntica. Para controlar la semilla, usen la función de R set.seed(seed, kind = NULL) justo antes de invocar el comando moran.mc para obtner el mismo valor cada vez.
Moran Columbus permutaciones:
set.seed(12345678) # Para siempre obtener el mismo resultado de las permutaciones
morpermRES <- moran.mc(col.e, col.listw, 99)
morpermRES
##
## Monte-Carlo simulation of Moran I
##
## data: col.e
## weights: col.listw
## number of simulations + 1: 100
##
## statistic = 0.21237, observed rank = 96, p-value = 0.04
## alternative hypothesis: greater
head(morpermRES$res)
## [1] -0.08887806 -0.08059128 0.07427797 0.03426065 -0.21063126 -0.06661672
morpermRES <- moran.mc(col.e, col.listw, 999)
morpermRES
##
## Monte-Carlo simulation of Moran I
##
## data: col.e
## weights: col.listw
## number of simulations + 1: 1000
##
## statistic = 0.21237, observed rank = 992, p-value = 0.008
## alternative hypothesis: greater
head(morpermRES$res)
## [1] 0.02497607 0.07061868 -0.04435387 0.01270686 -0.15458315 0.03836703
Práctica
Prueben la autocorrelación residual espacial en la regresión de Boston usando boston.soi como los pesos espaciales (nb).
Autocorrelación residual espacial Boston:
La estadística de prueba I de Moran tiene un alto poder estadístico frente a una gama de alternativas espaciales.
Sin embargo, no proporciona mucha ayuda en términos de qué modelo alternativo sería más apropiado.
Repaso:
data(columbus)
data(boston)
columbus.lm <- lm( CRIME ~ INC + HOVAL, data=columbus)
col.listw <- nb2listw(col.gal.nb)
bo1 <- lm(log(MEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) + AGE + log(DIS) + log(RAD) + PTRATIO+B+ log(LSTAT), data=boston.c)
summary(bo1)
##
## Call:
## lm(formula = log(MEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) +
## I(RM^2) + AGE + log(DIS) + log(RAD) + PTRATIO + B + log(LSTAT),
## data = boston.c)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71972 -0.09234 -0.00722 0.09929 0.78515
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5881621 0.1558554 29.439 < 2e-16 ***
## CRIM -0.0123607 0.0012496 -9.892 < 2e-16 ***
## ZN -0.0004927 0.0004823 -1.022 0.307485
## INDUS -0.0031181 0.0021742 -1.434 0.152156
## CHAS1 0.1041558 0.0333501 3.123 0.001895 **
## I(NOX^2) -0.6791523 0.1137241 -5.972 4.49e-09 ***
## I(RM^2) 0.0063361 0.0013265 4.777 2.35e-06 ***
## AGE 0.0001208 0.0005319 0.227 0.820459
## log(DIS) -0.1845076 0.0336954 -5.476 6.95e-08 ***
## log(RAD) 0.0502902 0.0139468 3.606 0.000343 ***
## PTRATIO -0.0350866 0.0049307 -7.116 3.94e-12 ***
## B 0.0003880 0.0001040 3.732 0.000212 ***
## log(LSTAT) -0.3732814 0.0252727 -14.770 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1844 on 493 degrees of freedom
## Multiple R-squared: 0.8013, Adjusted R-squared: 0.7964
## F-statistic: 165.6 on 12 and 493 DF, p-value: < 2.2e-16
bptest(bo1)
##
## studentized Breusch-Pagan test
##
## data: bo1
## BP = 85.506, df = 12, p-value = 3.638e-13
vif(bo1)
## CRIM ZN INDUS CHAS I(NOX^2) I(RM^2) AGE
## 1.715314 1.878820 3.303022 1.065301 3.721287 2.153562 3.328192
## log(DIS) log(RAD) PTRATIO B log(LSTAT)
## 4.907160 2.210203 1.691778 1.338040 3.423939
bos.listw <- nb2listw(boston.soi)
Estadísticas de la prueba del multiplicador de Lagrange para la autocorrelación espacial.
Las estadísticas del multiplicador de Lagrange si permiten una distinción entre modelos de error espacial (spatial error) y modelos de rezago espacial (spatial lag).
Multiplicador de Lagrange de Columbus:
columbus.lagrange <- lm.LMtests(columbus.lm,col.listw,
test=c("LMerr","RLMerr","LMlag","RLMlag"))
columbus.lagrange
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: col.listw
##
## LMerr = 4.6111, df = 1, p-value = 0.03177
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: col.listw
##
## RLMerr = 0.033514, df = 1, p-value = 0.8547
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: col.listw
##
## LMlag = 7.8557, df = 1, p-value = 0.005066
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: col.listw
##
## RLMlag = 3.2781, df = 1, p-value = 0.07021
Para especificar el modelo, primero evaluaríamos cuáles son los LMerr y LMlag significativos. En este caso, ambos lo son, lo que requiere la comparación de sus formas robustas.
Estas pruebas sugieren que el modelo de rezago espacial es la alternativa probable (p = 0.07 en comparación con p =0.85).
Práctica
Intenten correr la autocorrelación espacial en la regresión de Boston y sugieran la alternativa más probable de modelo.
Pueden experimentar con diferentes especificaciones de modelos para tratar de eliminar correlación espacial mediante la inclusión de diferentes variables. Podrían probar esto con el caso Columbus también. De hecho hay algunas especificaciones en las que no queda correlación residual.
Multiplicador de Lagrange de Boston:
La estimación del modelo de rezago espacial se realiza con la función
lagsarlm( )
mediante la estimación de máxima verosimilitud
(Maximum Likelihood). Los argumentos necesarios son una fórmula de
regresión, un conjunto de datos y un objeto de pesos espaciales
listw.
Estimación de spatial lag de Columbus:
columbus.lag <- lagsarlm(formula = CRIME ~ INC + HOVAL,
data = columbus,
listw = col.listw)
summary(columbus.lm)
##
## Call:
## lm(formula = CRIME ~ INC + HOVAL, data = columbus)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.418 -6.388 -1.580 9.052 28.649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.6190 4.7355 14.490 < 2e-16 ***
## INC -1.5973 0.3341 -4.780 1.83e-05 ***
## HOVAL -0.2739 0.1032 -2.654 0.0109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.43 on 46 degrees of freedom
## Multiple R-squared: 0.5524, Adjusted R-squared: 0.5329
## F-statistic: 28.39 on 2 and 46 DF, p-value: 9.341e-09
summary(columbus.lag)
##
## Call:
## lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, listw = col.listw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.4497093 -5.4565567 0.0016387 6.7159553 24.7107978
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 46.851431 7.314754 6.4051 1.503e-10
## INC -1.073533 0.310872 -3.4533 0.0005538
## HOVAL -0.269997 0.090128 -2.9957 0.0027381
##
## Rho: 0.40389, LR test value: 8.4179, p-value: 0.0037154
## Asymptotic standard error: 0.12071
## z-value: 3.3459, p-value: 0.00082027
## Wald statistic: 11.195, p-value: 0.00082027
##
## Log likelihood: -183.1683 for lag model
## ML residual variance (sigma squared): 99.164, (sigma: 9.9581)
## Number of observations: 49
## Number of parameters estimated: 5
## AIC: 376.34, (AIC for lm: 382.75)
## LM test for residual autocorrelation
## test value: 0.19184, p-value: 0.66139
En comparación con los resultados de MCO, todas las variables siguen siendo significativas.
El parámetro autorregresivo espacial (Rho) es altamente significativo, como se indica por el valor p de 0.0008 en una prueba t asintótica (basada en la matriz de varianza asintótica).La prueba de Likelihood Ratio (LR) sobre este parámetro también es altamente significativa (pvalue 0.0037).
También se observa una estadística de prueba LM para la autocorrelación restante del error espacial. Sin embargo, el valor de 0.19 no es significativo (p = 0,66) lo que sugiere que ya no tenemos problemas restantes con la especificación de la dependencia espacial.
El método de estimación de máxima verosimilitud predeterminado utiliza la descomposición de eigenvalues de la matriz de pesos espaciales con el cual se estima el modelo.
También nos interesan los residuales, e.g.,
columbus.lag$residuals
, y los valores de predicción
(valores ajustados) , e.g., columbus.lag$fitted.values
.
Con propósitos ilustrativos, comparamos este modelo con el estimador
de MCO (inconsistente). Primero, necesitamos construir una variable
dependiente rezagada espacialmente. Lo más fácil para manipular los
datos es utilizar el commando attach
. A continuación,
construimos una variable dependiente rezagada espacialmente usando la
función lag.listw
(especificamos un objeto de pesos
espaciales listw y una variable/vector). Finalmente, incluimos el rezago
espacial como variable explicativa en la regresión MCO. Noten cómo la
estimación y el error estándar del coeficiente para lagCRIME es bastante
diferente del resultado de máxima verosimilitud.
MCO minimizan la suma de cudrados de los residuales. ML maximiza la verosimilitud (likelihood) de observar nuestra muestra.
Estimación INCORRECTA de spatial lag de Columbus:
attach(columbus)
lagCRIME <- lag.listw(col.listw,CRIME)
wrong.lag <- lm(CRIME ~ lagCRIME + INC + HOVAL)
summary(wrong.lag)
##
## Call:
## lm(formula = CRIME ~ lagCRIME + INC + HOVAL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.393 -6.541 -0.190 6.799 23.485
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.07773 9.43653 4.247 0.000107 ***
## lagCRIME 0.52957 0.15612 3.392 0.001454 **
## INC -0.91054 0.36314 -2.507 0.015840 *
## HOVAL -0.26877 0.09312 -2.886 0.005970 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.32 on 45 degrees of freedom
## Multiple R-squared: 0.6436, Adjusted R-squared: 0.6198
## F-statistic: 27.08 on 3 and 45 DF, p-value: 3.672e-10
detach(columbus)
Después de ejecutar estos análisis, no olviden hacer
detach(columbus)
porque podríamos encontrarnos con
problemas de memoria. Una vez que hagamos detach un conjunto de datos,
las variables ya no están directamente disponibles por nombre, pero
debemos especificarlas como columnas en una tabla, por ejemplo, como
columbus$CRIME
.
Prueba heteroscedasticidad y normalidad (deberian ser valores similares KB y BP):
bptest.Sarlm(columbus.lag, studentize = TRUE)
##
## studentized Breusch-Pagan test
##
## data:
## BP = 10.855, df = 2, p-value = 0.004393
bptest.Sarlm(columbus.lag, studentize = FALSE)
##
## Breusch-Pagan test
##
## data:
## BP = 25.311, df = 2, p-value = 3.189e-06
Para la interpretación del modelo, tenemos que estudiar los efectos marginales (efecto directo vs efecto indirecto) en SAR ya que es un modelo global.
Adicionalmente, la comparación de los estimadores se puede analizar a través de los efects marginales de las variables explicativas. El efecto de las variables explicativas sobre la variable dependiente consta de dos partes en el modelo de rezago espacial. Uno es el efecto directo que se mide en cada ubicación. El efecto indirecto mide el spillover espacial. Uno de los beneficios de comparar los modelos de rezago espacial es examinar ambas partes explícitamente.
Interpretación Columbus:
#summary(impacts(columbus.lag, listw=col.listw, R=500),zstats=TRUE)
summary(columbus.lag)
##
## Call:
## lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, listw = col.listw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.4497093 -5.4565567 0.0016387 6.7159553 24.7107978
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 46.851431 7.314754 6.4051 1.503e-10
## INC -1.073533 0.310872 -3.4533 0.0005538
## HOVAL -0.269997 0.090128 -2.9957 0.0027381
##
## Rho: 0.40389, LR test value: 8.4179, p-value: 0.0037154
## Asymptotic standard error: 0.12071
## z-value: 3.3459, p-value: 0.00082027
## Wald statistic: 11.195, p-value: 0.00082027
##
## Log likelihood: -183.1683 for lag model
## ML residual variance (sigma squared): 99.164, (sigma: 9.9581)
## Number of observations: 49
## Number of parameters estimated: 5
## AIC: 376.34, (AIC for lm: 382.75)
## LM test for residual autocorrelation
## test value: 0.19184, p-value: 0.66139
Si no establecemos R, el método solo muestra los impactos directos, indirectos y totales de las variables en el modelo.
Efectos Marginales de Columbus:
impacts(columbus.lag, listw=col.listw)
## Impact measures (lag, exact):
## Direct Indirect Total
## INC -1.1225156 -0.6783818 -1.8008973
## HOVAL -0.2823163 -0.1706152 -0.4529315
summary(impacts(columbus.lag, listw=col.listw, R=500),zstats=TRUE)
## Impact measures (lag, exact):
## Direct Indirect Total
## INC -1.1225156 -0.6783818 -1.8008973
## HOVAL -0.2823163 -0.1706152 -0.4529315
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -1.1169 0.31677 0.01417 0.014166
## HOVAL -0.2856 0.09526 0.00426 0.003942
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -1.7080 -1.3298 -1.1136 -0.9063 -0.5131
## HOVAL -0.4668 -0.3503 -0.2843 -0.2195 -0.1008
##
## ========================================================
## Indirect:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -0.7244 0.3893 0.017408 0.018303
## HOVAL -0.1906 0.1151 0.005146 0.005146
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -1.7197 -0.8794 -0.6484 -0.4562 -0.20979
## HOVAL -0.4899 -0.2496 -0.1696 -0.1090 -0.03851
##
## ========================================================
## Total:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -1.8413 0.5911 0.026435 0.026642
## HOVAL -0.4762 0.1880 0.008408 0.007777
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -3.0901 -2.1626 -1.7645 -1.4643 -0.7859
## HOVAL -0.8875 -0.5806 -0.4591 -0.3558 -0.1699
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## INC 0.31677186 0.3892568 0.5910960
## HOVAL 0.09525991 0.1150587 0.1880052
##
## Simulated z-values:
## Direct Indirect Total
## INC -3.526024 -1.860959 -3.115122
## HOVAL -2.998395 -1.656616 -2.533094
##
## Simulated p-values:
## Direct Indirect Total
## INC 0.00042185 0.062750 0.0018387
## HOVAL 0.00271406 0.097597 0.0113061
Si se proporciona R, se utilizarán simulaciones para crear distribuciones para las medidas de impacto, ajustado por la matriz de covarianza de los coeficientes. Las simulaciones se realizan utilizando mvrnorm con los coeficientes y la matriz de covarianza del modelo ajustado.
INTERPRETACION: El efecto de un cambio en una variable de política x en i se extiende más allá de i a sus vecinos, vecinos de vecinos.
Las simulaciones se almacenan como objetos mcmc como se define en el paquete coda; los objetos se utilizan por conveniencia, pero no son generados por un proceso MCMC. Los valores simulados de los coeficientes se verifican para ver que el coeficiente espacial permanezca dentro de su intervalo válido.
Práctica
Especifiquen y estimen el modelo, interpretando sus resultados.
Estimación datos Boston:
Prueba heteroscedasticidad y normalidad de Boston (deberian ser valores similares KB y BP):
Para la interpretación del modelo, tenemos que estudiar los efectos marginales (efecto directo vs efecto indirecto) en SAR ya que es un modelo global.
Efectos Marginales de Boston:
La estimación del modelo de error espacial con máxima verosimilitud
es similar al procedimiento del modelo de rezago espacial y se
implementa con el comando errorsarlm
. Nuevamente, se debe
especificar la fórmula, el conjunto de datos y un objeto de pesos
espaciales listw, como se ilustra a continuación.
Estimación de spatial error de Columbus:
columbus.err <- errorsarlm(CRIME ~ INC + HOVAL,data=columbus, col.listw)
summary(columbus.err)
##
## Call:
## errorsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, listw = col.listw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.45950 -6.21730 -0.69775 7.65256 24.23631
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 61.053618 5.314875 11.4873 < 2.2e-16
## INC -0.995473 0.337025 -2.9537 0.0031398
## HOVAL -0.307979 0.092584 -3.3265 0.0008794
##
## Lambda: 0.52089, LR test value: 6.4441, p-value: 0.011132
## Asymptotic standard error: 0.14129
## z-value: 3.6868, p-value: 0.00022713
## Wald statistic: 13.592, p-value: 0.00022713
##
## Log likelihood: -184.1552 for error model
## ML residual variance (sigma squared): 99.98, (sigma: 9.999)
## Number of observations: 49
## Number of parameters estimated: 5
## AIC: 378.31, (AIC for lm: 382.75)
Observen cómo la estimación del coeficiente autorregresivo espacial es significativa, pero en menor grade que para el modelo de rezago espacial. Una comparación directa entre los dos modelos se puede basar en la maximized log-likelihood (no en medidas de ajuste como la R cuadrada). En este ejemplo, la probabilidad logarítmica del modelo de error de -184.16 es peor (menor que) la del modelo de rezago de -183.17. Sin embargo, los valores de esta estadística no se interpretan y la diferencia entre las probabilidades no se puede comparar de manera más formal en una prueba de Likelihood Ratio porque los dos modelos no están anidados (es decir, uno no se puede obtener del otro estableciendo coeficientes o funciones de coeficientes en cero). Más específicamente, uno no debe usar el método LR.sarlm en spdep para tal comparación, esto solo es válido para modelos anidadas.
Prueba heteroscedasticidad y normalidad (deberian ser valores similares KB y BP):
bptest.Sarlm(columbus.err, studentize = TRUE)
##
## studentized Breusch-Pagan test
##
## data:
## BP = 9.3694, df = 2, p-value = 0.009235
bptest.Sarlm(columbus.err, studentize = FALSE)
##
## Breusch-Pagan test
##
## data:
## BP = 16.285, df = 2, p-value = 0.0002908
Práctica
Estimen un modelo de error espacial para los datos de Boston. Compare las estimaciones de los coeficientes y sus errores estándar con los resultados de MCO y con los resultados del modelo de rezago espacial. ¿Cómo se compara el ajuste del modelo (maximized log-likelihood)?
Estimación de spatial error de Boston:
Prueba heteroscedasticidad y normalidad de Boston (SEM)
Analisis de factor común (common factor analysis)
Una prueba de especificación importante en el modelo de error espacial es la llamada hipótesis del factor común espacial. Esto parte de un modelo de error espacial que también se puede especificar en forma de rezago espacial, incluyendo las variables explicativas espacialmente rezagadas pero con restricciones en los parámetros (las restricciones del factor común).
La forma de rezago espacial del modelo de error también se denomina
modelo de Durbin (SDM). En spdep
, esto se implementa
especificando la opción mixted
para la estimación de rezago
espacial (no es necesario especificar las variables explicativas
rezagados espacialmente).
Este modelo nos ayuda a discriminar entre dependencia sustantiva y residual en una ecuación aparentemente mal especificada.
Estimación de spatial Durbin de Columbus:
columbus.durbin <- lagsarlm(CRIME ~ INC+HOVAL,data=columbus,col.listw,type="mixed")
summary(columbus.durbin)
##
## Call:
## lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, listw = col.listw,
## type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.15904 -6.62594 -0.39823 6.57561 23.62757
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 45.592893 13.128679 3.4728 0.0005151
## INC -0.939088 0.338229 -2.7765 0.0054950
## HOVAL -0.299605 0.090843 -3.2980 0.0009736
## lag.INC -0.618375 0.577052 -1.0716 0.2838954
## lag.HOVAL 0.266615 0.183971 1.4492 0.1472760
##
## Rho: 0.38251, LR test value: 4.1648, p-value: 0.041272
## Asymptotic standard error: 0.16237
## z-value: 2.3557, p-value: 0.018488
## Wald statistic: 5.5493, p-value: 0.018488
##
## Log likelihood: -182.0161 for mixed model
## ML residual variance (sigma squared): 95.051, (sigma: 9.7494)
## Number of observations: 49
## Number of parameters estimated: 7
## AIC: 378.03, (AIC for lm: 380.2)
## LM test for residual autocorrelation
## test value: 0.101, p-value: 0.75063
Nótese cómo el coeficiente autorregresivo espacial de 0.38 es diferente del valor obtenido con el modelo de error puro (0.52).
Asimismo, no hay evidencia de autocorrelación de error espacial restante, como lo indica la prueba significativo de LM (0.75)
Interpretación de spatial Durbin de Columbus:
#summary(impacts(columbus.lag, listw=col.listw, R=500),zstats=TRUE)
impacts(columbus.durbin, listw=col.listw)
## Impact measures (mixed, exact):
## Direct Indirect Total
## INC -1.0418080 -1.4804246 -2.52223256
## HOVAL -0.2836325 0.2302055 -0.05342697
Una prueba sobre la hipótesis del factor común es una prueba sobre las restricciones de que cada coeficiente de una variable explicativa rezagada espacialmente (por ejemplo, lag.INC) es igual al valor negativo del producto del coeficiente autorregresivo espacial (Rho) y el coeficiente de regresión coincidente de la variable no retardada (INC) = \(-\rho\beta\).
En otras palabras, la hipótesis nula es que \(-\rho\beta = \theta\), en la especificación del modelo espacial de Durbin:
\[y=\rho Wy +X\beta+WX\theta+\epsilon\]
En spdep
, esto se puede ilustrar mediante la función
LR.Sarlm
. Esta es una prueba de likelihood ratio, que
compara dos objetos para los que existe una función logLIK
.
En nuestro ejemplo, tenemos columbus.durbin para el modelo “sin
restricciones” (el modelo sin restricciones en los parámetros) y
columbus.err para el modelo “restringido” (las restricciones se han
aplicado). El ajuste (log-likelihood) del modelo sin restricciones (SDM)
siempre será mayor que el del modelo restringido.
NOTA: En estadística, la prueba de log-likelihood (LR test) evalúa el ajuste entre dos modelos estadísticos que compiten en función de la razón de sus probabilidades, específicamente uno maximizando todo el espacio de parámetros y el otro estimado después de imponer alguna restricción.
Diagnóstico de spatial Durbin de Columbus - Likelihood Ratio:
durbin.test1 <- LR.Sarlm(columbus.durbin,columbus.err)
print(durbin.test1)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 4.2782, df = 2, p-value = 0.1178
## sample estimates:
## Log likelihood of columbus.durbin Log likelihood of columbus.err
## -182.0161 -184.1552
Esto sugeriría que la hipótesis del factor común no se rechaza con un valor de p de 0.12 entonces nos quedaríamos con la especificación del modelo de error espacial (SEM).
Si hubiéramos rechazado la prueba de hipótesis (\(\theta=0\)), eso NO implica que tenemos que especificar el modelo de rezago espacial (spatial lag). Puede ser que \(\beta\) o \(\rho\) sean cero. SDM no es una versión anidada del spatial lag.
En otras palabras, la especificación del modelo de error espacial es inconsistente. Esto puede deberse a una mala elección de los pesos espaciales o a una mala especificación del proceso de error (por ejemplo, cuando el verdadero proceso de error no es SAR). Esta evidencia nos muestra que la especificación del modelo debe ser reconsiderada.
Práctica
Comprueben la hipótesis del factor común espacial para el modelo de Boston.
Ejemplo del modelo SLX (Spatial Lag of X), con una variable explanatoria rezagada
COL.SLX <- lmSLX(CRIME ~ INC + HOVAL + I(HOVAL^2), data=columbus, listw=col.listw, Durbin=~INC)
summary(COL.SLX)
##
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
## data = as.data.frame(x), weights = weights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35.720 -6.914 0.665 7.542 25.529
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 83.626625 9.324678 8.968 1.72e-11 ***
## INC -1.085603 0.387235 -2.803 0.00749 **
## HOVAL -0.649836 0.448598 -1.449 0.15454
## I.HOVAL.2. 0.003584 0.004120 0.870 0.38905
## lag.INC -0.977342 0.456277 -2.142 0.03777 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.99 on 44 degrees of freedom
## Multiple R-squared: 0.6042, Adjusted R-squared: 0.5682
## F-statistic: 16.79 on 4 and 44 DF, p-value: 1.994e-08
Aunque spdep no contiene una funcionalidad específica para abordar la heterogeneidad espacial, es posible obtener algunos resultados útiles manipulando el formato de regresión y extrayendo información del objeto de regresión espacial ajustado. Este ejercicio ilustra el análisis de heterogeneidad espacial discreta en forma de análisis espacial de varianza (ANOVA espacial) y pruebas para la presencia de regímenes espaciales.
El análisis espacial de la varianza se puede implementar como una regresión de una variable dummy. En nuestro ejemplo, HR90 es la variable dependiente, con la variable dummy SOUTH como variable explicativa. La estimación de la dummy indica la diferencia entre el régimen SUR = 1 y la media global.
Las regresiones de régimen espacial permiten que los coeficientes del modelo varíen entre subconjuntos de datos espaciales discretos. Implementaremos esto aprovechando la función de interacción de la fórmula en R (indicada por el símbolo : entre dos conjuntos de variables). Específicamente, crearemos una variable dummy para cada régimen (es decir, tomando un valor de uno para las observaciones en el régimen y cero para todas las demás) y luego interactuaremos cada variable explicativa con cada dummy Una prueba de Chow y su extensión a las regresiones espaciales (prueba de Chow espacial) es la base para evaluar la importancia de los regímenes.
Cargar datos de homicidios en Estados Unidos:
nat60 <- st_read("./data/natregimes.geojson")
## Reading layer `natregimes' from data source
## `/Users/irenefarah/Downloads/sp_econ/data/natregimes.geojson'
## using driver `GeoJSON'
## Simple feature collection with 3085 features and 73 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.7314 ymin: 24.95597 xmax: -66.96985 ymax: 49.37173
## CRS: 4326
Crear pesos reina de homicidios en Estados Unidos:
natq <- nb2listw(poly2nb(nat60, row.names = nat60$FIPSNO))
#summary(nat60)
Creamos dos nuevas variables y las agregamos a la tabla nat60. Primero, construimos el complemento de la variable ficticia SUR, NOSOUTH. A continuación, construimos un vector de unos de la misma longitud que las otras variables. Esto se hace fácilmente agregando SOUTH y NOSOUTH sin tener que preocuparnos por especificar la longitud. Llamamos a la nueva tabla de datos nat60a.
Crear variable dicotómica:
nat60$NOSOUTH = 1 - nat60$SOUTH
nat60$ONE = nat60$SOUTH + nat60$NOSOUTH
nat60a <- cbind(nat60,nat60$ONE,nat60$NOSOUTH)
summary(nat60a$NOSOUTH)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.5423 1.0000 1.0000
summary(nat60a$SOUTH)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4577 1.0000 1.0000
Primero ejecutaremos una regresión base de MCO para toda la nación, en la que los coeficientes se mantienen fijos y sin variable dummy SUR. Esto se denominará el modelo restringido (ya que los coeficientes son iguales en todos los regímenes).
MCO
nat60ols <- lm(HR60 ~ RD60 + PS60 + MA60 + DV60 + UE60, data=nat60a)
summary(nat60ols)
##
## Call:
## lm(formula = HR60 ~ RD60 + PS60 + MA60 + DV60 + UE60, data = nat60a)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.999 -2.427 -0.821 1.434 90.954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.51667 0.61064 17.222 < 2e-16 ***
## RD60 2.69490 0.09653 27.919 < 2e-16 ***
## PS60 0.60360 0.09141 6.603 4.72e-11 ***
## MA60 -0.27212 0.01935 -14.063 < 2e-16 ***
## DV60 1.29591 0.09599 13.501 < 2e-16 ***
## UE60 -0.10610 0.03565 -2.976 0.00294 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.84 on 3079 degrees of freedom
## Multiple R-squared: 0.2673, Adjusted R-squared: 0.2661
## F-statistic: 224.7 on 5 and 3079 DF, p-value: < 2.2e-16
Para obtener la regresión sin restricciones, debemos eliminar la constante de la fórmula para asegurarnos de que se estime una intersección diferente para cada régimen. Esto se hace especificando 0 en la fórmula para la especificación de regresión. También necesitamos interactuar cada una de las variables explicativas (incluida la constante, por eso construimos la variable de ONE) con las dummies SOUTH y NOSOUTH.
MCO permitiendo variable dicotómica, incorporando régimen espacial
nat60rega <- lm(HR60 ~ 0 + (ONE + RD60 + PS60 +
MA60 + DV60 + UE60):(SOUTH + NOSOUTH), data=nat60a)
summary(nat60rega)
##
## Call:
## lm(formula = HR60 ~ 0 + (ONE + RD60 + PS60 + MA60 + DV60 + UE60):(SOUTH +
## NOSOUTH), data = nat60a)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.527 -2.182 -0.647 1.339 88.482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ONE:SOUTH 13.21547 0.87183 15.158 < 2e-16 ***
## ONE:NOSOUTH 6.39963 0.99955 6.403 1.76e-10 ***
## RD60:SOUTH 1.76448 0.15369 11.481 < 2e-16 ***
## RD60:NOSOUTH 1.85730 0.23653 7.852 5.60e-15 ***
## PS60:SOUTH 0.29930 0.16611 1.802 0.0717 .
## PS60:NOSOUTH 0.37748 0.11730 3.218 0.0013 **
## MA60:SOUTH -0.27521 0.02951 -9.326 < 2e-16 ***
## MA60:NOSOUTH -0.19537 0.02978 -6.560 6.27e-11 ***
## DV60:SOUTH 1.17945 0.18879 6.247 4.75e-10 ***
## DV60:NOSOUTH 1.11883 0.11463 9.761 < 2e-16 ***
## UE60:SOUTH -0.29186 0.05518 -5.289 1.31e-07 ***
## UE60:NOSOUTH 0.09277 0.04554 2.037 0.0417 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.721 on 3073 degrees of freedom
## Multiple R-squared: 0.5746, Adjusted R-squared: 0.5729
## F-statistic: 345.9 on 12 and 3073 DF, p-value: < 2.2e-16
Tener en cuenta que la R-cuadrada no está estimada correctamente. Los resultados muestran algunas diferencias significativas entre regímenes para algunos de los coeficientes, como PS60 (no significativo en SOUTH pero significativo en NOSOUTH) y UE60 (fuertemente significativo y negativo en SOUTH y positivo pero marginalmente significativo en NOSOUTH).
Una prueba directa para ver si los coeficientes son iguales en todos los regímenes es la prueba de Chow contenida en cualquier prueba econométrica estándar donde er son el vector de residuales de la regresión restringida y eu son los residuales de la regresión no restringida.
Obtener residuales de la regresión
er <- residuals(nat60ols)
#Extraemos grados de libertad
k <- nat60ols$rank
n2k <- nat60ols$df.residual - k
Generar prueba de Chow
chow.test <- function(rest,unrest)
{
er <- residuals(rest)
eu <- residuals(unrest)
er2 <- sum(er^2)
eu2 <- sum(eu^2)
k <- rest$rank
n2k <- rest$df.residual - k
c <- ((er2 - eu2)/k) / (eu2 / n2k)
pc <- pf(c,k,n2k,lower.tail=FALSE)
list(c,pc,k,n2k)
}
chow.test(nat60ols,nat60rega)
## [[1]]
## [1] 27.11882
##
## [[2]]
## [1] 1.141437e-31
##
## [[3]]
## [1] 6
##
## [[4]]
## [1] 3073
El valor del estadístico de prueba es 27,12, que es muy significativo. Por lo tanto, es prueba que los coeficientes del modelo no son constantes entre regímenes, lo que indica heterogeneidad espacial.
Prueba de Multiplicador de Lagrange:
lm.LMtests(nat60rega,natq,test=c("LMerr","RLMerr","LMlag","RLMlag"))
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = HR60 ~ 0 + (ONE + RD60 + PS60 + MA60 + DV60 +
## UE60):(SOUTH + NOSOUTH), data = nat60a)
## weights: natq
##
## LMerr = 97.676, df = 1, p-value < 2.2e-16
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = HR60 ~ 0 + (ONE + RD60 + PS60 + MA60 + DV60 +
## UE60):(SOUTH + NOSOUTH), data = nat60a)
## weights: natq
##
## RLMerr = 9.8075, df = 1, p-value = 0.001738
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = HR60 ~ 0 + (ONE + RD60 + PS60 + MA60 + DV60 +
## UE60):(SOUTH + NOSOUTH), data = nat60a)
## weights: natq
##
## LMlag = 129.19, df = 1, p-value < 2.2e-16
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = HR60 ~ 0 + (ONE + RD60 + PS60 + MA60 + DV60 +
## UE60):(SOUTH + NOSOUTH), data = nat60a)
## weights: natq
##
## RLMlag = 41.325, df = 1, p-value = 1.289e-10
Los residuos de la regresión continúan mostrando una fuerte evidencia de autocorrelación espacial. Tanto las estadísticas de LMerr como las de LMlag son muy significativas, al igual que sus formas robustas, con una ligera ventaja a favor de la alternativa del modelo de rezago espacial.
Los regímenes espaciales se pueden incorporar en los modelos de error espacial y rezago espacial de la misma manera que para MCO, por medio de los términos de interacción en la fórmula del modelo.
Estimación del Spatial Lag con regímenes espaciales
#Se tarda bastante
hr60rlag <- lagsarlm(HR60 ~ 0 + (ONE + RD60 + PS60 + MA60 + DV60 + UE60):(SOUTH + NOSOUTH),
data=nat60a,natq)
summary(hr60rlag)
##
## Call:lagsarlm(formula = HR60 ~ 0 + (ONE + RD60 + PS60 + MA60 + DV60 +
## UE60):(SOUTH + NOSOUTH), data = nat60a, listw = natq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.41352 -2.11312 -0.62028 1.33659 88.38962
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## ONE:SOUTH 9.739325 0.895945 10.8704 < 2.2e-16
## ONE:NOSOUTH 5.273615 0.982158 5.3694 7.899e-08
## RD60:SOUTH 1.420458 0.153772 9.2374 < 2.2e-16
## RD60:NOSOUTH 1.625286 0.231632 7.0167 2.272e-12
## PS60:SOUTH 0.208484 0.162155 1.2857 0.198544
## PS60:NOSOUTH 0.385196 0.114367 3.3681 0.000757
## MA60:SOUTH -0.213140 0.029087 -7.3277 2.340e-13
## MA60:NOSOUTH -0.165157 0.029204 -5.6553 1.556e-08
## DV60:SOUTH 0.979970 0.184224 5.3195 1.041e-07
## DV60:NOSOUTH 0.898503 0.112754 7.9687 1.554e-15
## UE60:SOUTH -0.210014 0.053885 -3.8974 9.722e-05
## UE60:NOSOUTH 0.083777 0.044455 1.8845 0.059494
##
## Rho: 0.25961, LR test value: 106.05, p-value: < 2.22e-16
## Asymptotic standard error: 0.02568
## z-value: 10.11, p-value: < 2.22e-16
## Wald statistic: 102.21, p-value: < 2.22e-16
##
## Log likelihood: -9106.562 for lag model
## ML residual variance (sigma squared): 21.189, (sigma: 4.6032)
## Number of observations: 3085
## Number of parameters estimated: 14
## AIC: 18241, (AIC for lm: 18345)
## LM test for residual autocorrelation
## test value: 38.109, p-value: 6.6916e-10
Aquí nuevamente hay evidencia de una diferencia entre regímenes para los coeficientes de PS60 y UE60. Para evaluar esto de manera más rigurosa, también necesitamos los resultados del modelo de rezago espacial de la especificación restringida.
Estimación del Spatial Lag con regímenes espaciales (especificación restringida)
#Se tarda bastante
hr60lag <- lagsarlm(HR60 ~ RD60 + PS60 + MA60 +DV60 + UE60,data=nat60a,natq)
summary(hr60lag)
##
## Call:
## lagsarlm(formula = HR60 ~ RD60 + PS60 + MA60 + DV60 + UE60, data = nat60a,
## listw = natq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.89718 -2.17000 -0.74723 1.37653 90.09303
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.048613 0.629963 11.1889 < 2.2e-16
## RD60 1.912544 0.106597 17.9419 < 2.2e-16
## PS60 0.456590 0.088178 5.1781 2.242e-07
## MA60 -0.190804 0.019269 -9.9023 < 2.2e-16
## DV60 0.942877 0.093493 10.0850 < 2.2e-16
## UE60 -0.052096 0.034187 -1.5239 0.1275
##
## Rho: 0.32721, LR test value: 192.26, p-value: < 2.22e-16
## Asymptotic standard error: 0.024022
## z-value: 13.621, p-value: < 2.22e-16
## Wald statistic: 185.53, p-value: < 2.22e-16
##
## Log likelihood: -9143.04 for lag model
## ML residual variance (sigma squared): 21.528, (sigma: 4.6399)
## Number of observations: 3085
## Number of parameters estimated: 8
## AIC: 18302, (AIC for lm: 18492)
## LM test for residual autocorrelation
## test value: 50.95, p-value: 9.4769e-13
Cuando los términos de error de regresión se autocorrelacionan espacialmente, la forma simple de la prueba de Chow que se muestra arriba ya no es válida entonces hacemos una Chow espacial.
Prueba de Chow espacial:
spatialchow.test <- function(rest,unrest)
{
lrest <- rest$LL
lunrest <- unrest$LL
k <- rest$parameters - 2
spchow <- - 2.0 * (lrest - lunrest)
pchow <- pchisq(spchow,k,lower.tail=FALSE)
list(spchow,pchow,k)
}
spatialchow.test(hr60lag,hr60rlag)
## [[1]]
## [,1]
## [1,] 72.95666
##
## [[2]]
## [,1]
## [1,] 1.01042e-13
##
## [[3]]
## [1] 6
El resultado muestra un rechazo a la hipótesis nula, sugiriendo coeficientes significativamente diferentes en cada uno de los regímenes, incluso después de corregir por el rezago espacial.
Como primer ejemplo de modelado de heterogeneidad espacial continua, consideramos el método de expansión de Casetti. Este es un caso especial de un modelo con coeficientes variables, donde la variación es generada por una ecuación de expansión. En el ejemplo más simple, la ecuación de expansión es una superficie de tendencia lineal. Esto produce una ecuación final que contiene las variables explicativas originales, así como términos de interacción con las coordenadas X e Y. El cálculo produce valores de coeficientes individuales para cada ubicación, que pueden luego ser mapeados y analizados más a fondo.
Un problema común en la implementación del método de expansión espacial es el alto grado de multicolinealidad que resulta de los términos de interacción. En la práctica, esto requiere el uso de algunos ajustes, como reemplazar los términos originales de interacción por sus componentes principales. En este ejemplo ignoraremos la posible multicolinealidad.
Preparar GWR
data(columbus)
attach(columbus)
Expansión espacial
colex <- lm(CRIME ~ (INC + HOVAL)*(X + Y))
summary(colex)
##
## Call:
## lm(formula = CRIME ~ (INC + HOVAL) * (X + Y))
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.5556 -7.6351 -0.6181 7.8363 30.1948
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.97559 69.36676 1.571 0.1241
## INC -5.82949 3.84408 -1.516 0.1373
## HOVAL 0.27337 0.82049 0.333 0.7407
## X -0.76287 1.13692 -0.671 0.5061
## Y -0.26332 1.21420 -0.217 0.8294
## INC:X -0.01854 0.05396 -0.344 0.7329
## INC:Y 0.13949 0.08004 1.743 0.0891 .
## HOVAL:X 0.03159 0.01549 2.040 0.0480 *
## HOVAL:Y -0.05034 0.02196 -2.293 0.0272 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.04 on 40 degrees of freedom
## Multiple R-squared: 0.6375, Adjusted R-squared: 0.5649
## F-statistic: 8.791 on 8 and 40 DF, p-value: 7.663e-07
Se muestra que muy pocos coeficientes son significativos, aunque la
R-cuadrada es más que aceptable. Este es un síntoma típico de la alta
multicolinealidad. Para nuestros propósitos, ignoraremos este problema
por ahora y calcularemos coeficientes individuales de HOVAL para cada
localidad. Esto se logra extrayendo los coeficientes relevantes del
objeto de regresión usando llamando colex$coeficientes
.
Extraer coeficientes
b <- colex$coefficients
b[[3]]
## [1] 0.2733714
bihoval <- b[3] + b[8] * X + b[9] * Y
bihoval
## [1] -0.71945869 -0.73484522 -0.54173838 -0.61340248 -0.37564117 -0.32192096
## [7] -0.60638062 -0.51564482 -0.16255859 -0.05598467 -0.35829853 -0.37198398
## [13] -0.42336113 -0.39035215 -0.23271508 -0.25443306 0.07333582 -0.29218860
## [19] -0.34176508 0.14326490 -0.18138317 -0.04092354 0.19700015 -0.18584181
## [25] -0.13841490 -0.09558335 0.07053751 0.02287442 -0.01427011 -0.04484941
## [31] -0.34667265 0.35074015 0.13619297 -0.19143232 0.18497546 -0.28527600
## [37] 0.03312495 0.12107370 -0.30416609 0.39765999 0.49266738 -0.14792599
## [43] 0.18759738 0.26427001 0.21424089 -0.21628436 0.61049050 0.27143584
## [49] 0.36488625
summary(bihoval)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.7348 -0.3418 -0.1479 -0.1096 0.1362 0.6105
Vemos que los coeficientes variables oscilan de -0.7348 a 0.6105 (comparando con una estimación de -0.2739 en el modelo de coeficiente fijo).
bi <- data.frame(POLYID,bihoval)
Primero exploraremos los conceptos básicos de la función gwr. Esta función genera una tabla con las estimaciones de coeficientes individuales para cada variable y cada ubicación. Los argumentos de esta función son una fórmula (para la especificación del modelo), un conjunto de datos, una matriz con las coordenadas X, Y de las ubicaciones (o centroides de polígonos) y un “bandwidth”. La función de la kernel por defecto es la gaussiana. El bandwidth es difícil de especificar a priori y el enfoque preferido es realizar una validación cruzada, minimizando el error de predicción del promedio de la raíz cuadrada. Por ahora, probaremos tres valores: 20, 3 y 2. A medida que se amplía el bandwidth, la superficie de los coeficientes estimados se suaviza. Con la regresión CRIME, las coordenadas como las variables X e Y del conjunto de datos incorporado y el bandwidth de 20, esto produce:
attach(columbus)
colg1 <- gwr(CRIME ~ INC + HOVAL,data=columbus, coords=cbind(X,Y),bandwidth=20)
colg1
## Call:
## gwr(formula = CRIME ~ INC + HOVAL, data = columbus, coords = cbind(X,
## Y), bandwidth = 20)
## Kernel function: gwr.Gauss
## Fixed bandwidth: 20
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max. Global
## X.Intercept. 68.41172 68.84934 69.06064 69.24627 69.95574 68.6190
## INC -1.65427 -1.62266 -1.60883 -1.58070 -1.52803 -1.5973
## HOVAL -0.30704 -0.29166 -0.27768 -0.26118 -0.23198 -0.2739
El comando print
proporciona un resumen de los
coeficientes estimados (tenga en cuenta que el summary
no
es útil en este contexto), enumera los parámetros básicos (kernel y
bandwidth) y proporciona estadísticas descriptivas para cada uno de los
coeficientes estimados. Por ejemplo, las estimaciones para HOVAL oscilan
entre -0.3407 y -0.1903 y no están tan lejos de la estimación global de
-0,2739. La razón principal de esto es el gran bandwidth, que no permite
mucha variabilidad. Comparen esto con el rango obtenido por el método de
expansión espacial de -0.7348 a 0.6105.
colg2 <- gwr(CRIME ~ INC + HOVAL,data=columbus, coords=cbind(X,Y),bandwidth=3)
colg2
## Call:
## gwr(formula = CRIME ~ INC + HOVAL, data = columbus, coords = cbind(X,
## Y), bandwidth = 3)
## Kernel function: gwr.Gauss
## Fixed bandwidth: 3
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max. Global
## X.Intercept. 24.823915 58.337575 66.888807 70.734368 77.517040 68.6190
## INC -2.910848 -2.055262 -1.370336 -0.489550 0.605579 -1.5973
## HOVAL -0.938862 -0.370668 -0.072330 -0.016077 0.417796 -0.2739
colg3 <- gwr(CRIME ~ INC + HOVAL, data=columbus, coords=cbind(X,Y),bandwidth=2)
colg3
## Call:
## gwr(formula = CRIME ~ INC + HOVAL, data = columbus, coords = cbind(X,
## Y), bandwidth = 2)
## Kernel function: gwr.Gauss
## Fixed bandwidth: 2
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max. Global
## X.Intercept. 22.991392 51.964036 62.876904 69.090317 81.631781 68.6190
## INC -3.386156 -1.956824 -0.702189 -0.325575 1.339862 -1.5973
## HOVAL -1.209784 -0.372961 -0.115808 0.038842 0.879194 -0.2739
El ancho de banda de 3 produce un rango similar al método de expansión, mientras que el ancho de banda de 2 da como resultado un rango muy amplio. Las tres distribuciones se pueden comparar con la del método de expansión por medio de un diagrama de caja. Para extraer el vector de coeficientes necesitamos acceder al atributo SDF de la clase gwr. Este es un marco de datos que contiene las coordenadas, coeficientes y una medida de error. Por ejemplo, para el segundo ancho de banda (2), esto parece (se han eliminado las columnas R2 y gwr.e):
colg3$SDF
## class : SpatialPointsDataFrame
## features : 49
## extent : 24.25, 51.24, 24.96, 44.07 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 7
## names : sum.w, (Intercept), INC, HOVAL, gwr.e, pred, localR2
## min values : 1.57640620130106, 22.9913921526746, -3.38615619344475, -1.2097842174163, -11.4861625510241, 0.718999007155328, 0.479815888669314
## max values : 6.80180393374637, 81.6317808690283, 1.33986223270747, 0.879193696261171, 11.0890845785619, 57.8029594214381, 0.985577303210009
hovg3 <- colg2$SDF$HOVAL
summary(hovg3)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.93886 -0.37067 -0.07233 -0.16594 -0.01608 0.41780
hovg20 <- colg1$SDF$HOVAL
summary(hovg20)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.3070 -0.2917 -0.2777 -0.2763 -0.2612 -0.2320
hovg2 <- colg3$SDF$HOVAL
summary(hovg2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.20978 -0.37296 -0.11581 -0.14918 0.03884 0.87919
boxplot(bihoval,hovg20,hovg3,hovg2, names=c("Expansion","bw=20","bw=3","bw=2"))
colgwr <- data.frame(POLYID,hovg20,hovg3,hovg2)
#head(colgwr)
Leer datos de Columbus para juntar los polígonos con nuestro análisis:
columbus_sf=st_read("./data/columbus.geojson")
## Reading layer `columbus' from data source
## `/Users/irenefarah/Downloads/sp_econ/data/columbus.geojson'
## using driver `GeoJSON'
## Simple feature collection with 49 features and 20 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 5.874907 ymin: 10.78863 xmax: 11.28742 ymax: 14.74245
## CRS: 4326
col.GWR=merge(columbus_sf, colgwr, by="POLYID")
Generar mapas:
tmap_mode("plot")
tm_shape(col.GWR) +
tm_borders(alpha=0.7)+
tm_fill(col = "hovg20",
alpha = 0.5,
style = "quantile",
title = "GWR Bandwidth 20")
tmap_mode("plot")
tm_shape(col.GWR) +
tm_borders(alpha=0.7)+
tm_fill(col = "hovg3",
alpha = 0.5,
style = "quantile",
title = "GWR Bandwidth 3")
tmap_mode("plot")
tm_shape(col.GWR) +
tm_borders(alpha=0.7)+
tm_fill(col = "hovg2",
alpha = 0.5,
style = "quantile",
title = "GWR Bandwidth 2")