[Descargar - RMarkdown]

Regresión espacial

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)

¿Cómo saber qué modelo especificar?

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:

Spatial Lag (Modelo de rezago espacial)

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:

Spatial Error Model (Modelo de error espacial)

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)

Spatial Durbin Model (Modelo de Durbin espacial SDM)

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

Heterogeneidad

Incorporar heterogeneidad espacial (discreta)

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.

Heterogeneidad discreta - régimen espacial

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.

Heterogeneidad continua

Metodo de expansión 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)

GWR

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