Laboratorio 3 - Modelos Espaciales 2
Yasna Cortés
11/22/2020
Objetivos del Laboratorio
En este laboratorio, las/los estudiantes deberán utilizar sus conocimientos adquiridos en análisis exploratorio de datos espaciales, para implementar modelos de regresión lineales espaciales. De esta forma, los elementos que consideraremos para la realización de este laboratorio son:
- Exploración de datos espaciales.
- Modelos espaciales de regresión lineal.
Paquetes que se utilizarán
Antes de comenzar, debemos incluir los paquetes estadísticos necesarios para llevar a cabo el análisis exploratorio espacial, la construcción de matrices de peso espacial y medidas de autocorrelación. Preliminarmente, estos paquetes son: (Estos paquetes se incluyen de acuerdo al supuesto de que ya están instalados en R. En el caso de que no estén instalados, usar el comando install.package(“nombre del paquete”). También, es posible instalar los paquetes en el menú Herramientas (tools) –> instalar paquetes)
library(spatialreg)
library(spdep)
library(maptools)
library(RColorBrewer)
library(leaflet)
library(dplyr)
library(ggplot2)
library(tmap)
library(tmaptools)
library(lmtest)
library(normtest)
library(splm)
library(data.table)
library(sphet)Por favor, fijar el directorio de trabajo en un lugar de su computador que sea de fácil acceso.
CASO 1: Determinantes de los Precios de Vivienda en Boston
Información Espacial
En esta sección, usaremos la siguiente base de datos.
data(boston)
ls()## [1] "boston.c" "boston.soi" "boston.utm"data2<-boston.c
attach(data2)Esta base de datos está representada por centroides. Son 506 observaciones de los precios de vivienda de Boston. Con esta información, exploraremos los posibles determinantes de los precios de vivienda en Boston, usando como proxy de los precios de vivienda el precio promedio de 506 census tracts (blocks) y otras variables explicativas medidas al mismo nivel de datos.
Análisis Exploratorio Preliminar de los Datos.
map<-SpatialPointsDataFrame(data=data2, coords = cbind(LON,LAT))
colours=c("dark blue", "blue", "red", "dark red")
spplot(map, "MEDV", cuts=quantile(MEDV), col.regions=colours, cex=0.7,
contour=TRUE)
Housing Prices by Census Tracts in Boston
Otra forma de explorar datos espaciales a través de:
leaflet(data=data2) %>%
addTiles() %>%
addMarkers(clusterOptions = markerClusterOptions()) ## Assuming "LON" and "LAT" are longitude and latitude, respectivelyRegresión Lineal Estándar (SLM - OLS)
Para comenzar con el análisis, primero calcularemos un modelo de regresión lineal a través de OLS. De esta forma, el modelo que construiremos estará determinado por:
\[ y=X\beta +\epsilon \] Donde la variable dependiente \(y\) corresponde al precio de las viviendas, mientras que la matriz \(X\) contiene todos los posibles determinantes de los precios de vivienda.
model_SAR_1<-lm(MEDV ~ CRIM+RM+INDUS+NOX+AGE+DIS+RAD+PTRATIO+B+LSTAT+TAX, data=data2)
summary(model_SAR_1)##
## Call:
## lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS + RAD +
## PTRATIO + B + LSTAT + TAX, data = data2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.7429 -2.8887 -0.7514 1.8144 26.8277
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.308337 5.199690 7.175 2.66e-12 ***
## CRIM -0.103402 0.033339 -3.102 0.002035 **
## RM 4.074379 0.420639 9.686 < 2e-16 ***
## INDUS 0.018212 0.062015 0.294 0.769138
## NOX -17.829176 3.889690 -4.584 5.79e-06 ***
## AGE -0.002647 0.013353 -0.198 0.842957
## DIS -1.210182 0.186123 -6.502 1.94e-10 ***
## RAD 0.304603 0.066878 4.555 6.62e-06 ***
## PTRATIO -1.131146 0.126079 -8.972 < 2e-16 ***
## B 0.009853 0.002735 3.603 0.000346 ***
## LSTAT -0.525072 0.051543 -10.187 < 2e-16 ***
## TAX -0.010901 0.003710 -2.939 0.003452 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.838 on 494 degrees of freedom
## Multiple R-squared: 0.7293, Adjusted R-squared: 0.7233
## F-statistic: 121 on 11 and 494 DF, p-value: < 2.2e-16También, analizamos si el modelo estimado cumple con los supuestos básicos de regresión: normalidad y homogeneidad.
residuals_model_SAR_1<-residuals(model_SAR_1)
jb.norm.test(residuals_model_SAR_1)##
## Jarque-Bera test for normality
##
## data: residuals_model_SAR_1
## JB = 936.74, p-value < 2.2e-16bptest(model_SAR_1)##
## studentized Breusch-Pagan test
##
## data: model_SAR_1
## BP = 59.214, df = 11, p-value = 1.297e-08Análisis de los Residuos
Definimos la Matriz de Pesos Espaciales
De acuerdo a la tipología de los datos, usaremos una matriz de peso espacial basada en distancias. Por tanto, primero obtenemos las coordenadas de esta base de datos, para calcular la matriz W a través del siguiente procedimiento:
coords = cbind(LON,LAT)
dist.mat<-as.matrix(dist(coords, method="euclidean"))
dist.mat.inv<-1/dist.mat
diag(dist.mat.inv)<-0Para obtener la matriz W, basada en distancias usamos la siguiente función:
dist.W1<-mat2listw(dist.mat.inv, row.names = data2$TRACT, style="W")
summary(dist.W1)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 506
## Number of nonzero links: 255530
## Percentage nonzero weights: 99.80237
## Average number of links: 505
## Link number distribution:
##
## 505
## 506
## 506 least connected regions:
## 2011 2021 2022 2031 2032 2033 2041 2042 2043 2044 2045 2046 2047 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068 2069 2070 2071 2072 2081 2082 2083 2084 2091 2092 2101 2102 2103 2104 2105 2106 2107 2108 2109 2111 2112 2113 2114 2121 2141 2151 2161 2171 2172 2173 2174 2175 2176 2181 3301 3302 3311 3312 3313 3321 3322 3323 3324 3331 3332 3333 3334 3335 3336 3341 3342 3343 3344 3351 3352 3353 3354 3361 3362 3363 3364 3371 3372 3373 3381 3382 3383 3384 3385 3391 3392 3393 3394 3395 3396 3397 3398 3399 3400 3401 3411 3412 3413 3414 3415 3416 3417 3418 3419 3421 3422 3423 3424 3425 3426 3427 3501 3502 3503 3504 3505 3506 3507 3508 3509 3510 3511 3512 3513 3514 3515 3521 3522 3523 3524 3525 3526 3527 3528 3529 3530 3531 3532 3533 3534 3535 3536 3537 3538 3539 3540 3541 3542 3543 3544 3545 3546 3547 3548 3549 3550 3561 3562 3563 3564 3565 3566 3567 3571 3572 3573 3574 3575 3576 3577 3578 3581 3583 3584 3585 3586 3587 3591 3593 3602 3611 3612 3613 3651 3652 3661 3662 3671 3672 3681 3682 3683 3684 3685 3686 3687 3688 3689 3690 3691 3701 3702 3703 3704 3731 3732 3733 3734 3735 3736 3737 3738 3739 3740 3741 3742 3743 3744 3745 3746 3747 3748 3821 3822 3823 3824 3825 3826 3831 3832 3833 3834 3835 3836 3837 3838 3839 3840 3851 3852 3861 4001 4002 4003 4004 4005 4006 4007 4008 4009 4010 4011 4012 4021 4022 4023 4024 4025 4031 4032 4033 4034 4035 4041 4042 4043 4044 4051 4061 4071 4091 4111 4112 4113 4121 4122 4123 4131 4132 4133 4134 4135 4141 4142 4143 4151 4152 4153 4161 4162 4163 4164 4171 4172 4173 4174 4175 4176 4177 4178 4179 4180 4181 4182 4191 4192 4193 4194 4195 4196 4197 4198 4201 4202 4203 4211 4212 4221 4222 4223 4224 4225 4226 4227 4228 4231 5001 5011 5012 5021 5022 5031 5041 5051 5052 5061 5062 5071 5081 5082 1 2 3 4 5 6 7 8 101 102 104 105 107 108 201 202 203 301 302 401 402 403 404 405 406 407 501 502 503 504 506 507 508 509 510 511 512 601 602 603 604 605 606 607 608 609 610 612 613 614 702 703 705 706 707 708 709 710 801 802 803 805 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 901 902 903 904 905 906 907 908 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1101 1102 1103 1104 1105 1106 1201 1202 1203 1204 1205 1206 1207 1301 1302 1303 1304 1401 1402 1403 1404 1601 1602 1604 1605 1606 1701 1702 1703 1704 1705 1706 1707 1708 1801 1802 1803 1804 1805 with 505 links
## 506 most connected regions:
## 2011 2021 2022 2031 2032 2033 2041 2042 2043 2044 2045 2046 2047 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068 2069 2070 2071 2072 2081 2082 2083 2084 2091 2092 2101 2102 2103 2104 2105 2106 2107 2108 2109 2111 2112 2113 2114 2121 2141 2151 2161 2171 2172 2173 2174 2175 2176 2181 3301 3302 3311 3312 3313 3321 3322 3323 3324 3331 3332 3333 3334 3335 3336 3341 3342 3343 3344 3351 3352 3353 3354 3361 3362 3363 3364 3371 3372 3373 3381 3382 3383 3384 3385 3391 3392 3393 3394 3395 3396 3397 3398 3399 3400 3401 3411 3412 3413 3414 3415 3416 3417 3418 3419 3421 3422 3423 3424 3425 3426 3427 3501 3502 3503 3504 3505 3506 3507 3508 3509 3510 3511 3512 3513 3514 3515 3521 3522 3523 3524 3525 3526 3527 3528 3529 3530 3531 3532 3533 3534 3535 3536 3537 3538 3539 3540 3541 3542 3543 3544 3545 3546 3547 3548 3549 3550 3561 3562 3563 3564 3565 3566 3567 3571 3572 3573 3574 3575 3576 3577 3578 3581 3583 3584 3585 3586 3587 3591 3593 3602 3611 3612 3613 3651 3652 3661 3662 3671 3672 3681 3682 3683 3684 3685 3686 3687 3688 3689 3690 3691 3701 3702 3703 3704 3731 3732 3733 3734 3735 3736 3737 3738 3739 3740 3741 3742 3743 3744 3745 3746 3747 3748 3821 3822 3823 3824 3825 3826 3831 3832 3833 3834 3835 3836 3837 3838 3839 3840 3851 3852 3861 4001 4002 4003 4004 4005 4006 4007 4008 4009 4010 4011 4012 4021 4022 4023 4024 4025 4031 4032 4033 4034 4035 4041 4042 4043 4044 4051 4061 4071 4091 4111 4112 4113 4121 4122 4123 4131 4132 4133 4134 4135 4141 4142 4143 4151 4152 4153 4161 4162 4163 4164 4171 4172 4173 4174 4175 4176 4177 4178 4179 4180 4181 4182 4191 4192 4193 4194 4195 4196 4197 4198 4201 4202 4203 4211 4212 4221 4222 4223 4224 4225 4226 4227 4228 4231 5001 5011 5012 5021 5022 5031 5041 5051 5052 5061 5062 5071 5081 5082 1 2 3 4 5 6 7 8 101 102 104 105 107 108 201 202 203 301 302 401 402 403 404 405 406 407 501 502 503 504 506 507 508 509 510 511 512 601 602 603 604 605 606 607 608 609 610 612 613 614 702 703 705 706 707 708 709 710 801 802 803 805 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 901 902 903 904 905 906 907 908 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1101 1102 1103 1104 1105 1106 1201 1202 1203 1204 1205 1206 1207 1301 1302 1303 1304 1401 1402 1403 1404 1601 1602 1604 1605 1606 1701 1702 1703 1704 1705 1706 1707 1708 1801 1802 1803 1804 1805 with 505 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 506 256036 506 4.736809 2065.127También, podemos usar la matriz de peso espacial que viene por defecto con la base de datos:
dist.W2<-nb2listw(boston.soi)
summary(dist.W2)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 506
## Number of nonzero links: 2152
## Percentage nonzero weights: 0.8405068
## Average number of links: 4.252964
## Link number distribution:
##
## 1 2 3 4 5 6 7 8
## 10 39 93 143 143 52 22 4
## 10 least connected regions:
## 2033 2081 2181 3592 0001 0101 0108 0201 1301 1805 with 1 link
## 4 most connected regions:
## 3564 3744 4196 4223 with 8 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 506 256036 506 261.3891 2071.255Índice de Moran
lm.morantest(model_SAR_1, listw = dist.W1)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W1
##
## Moran I statistic standard deviate = 20.425, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 5.107215e-02 -4.346662e-03 7.361799e-06lm.morantest(model_SAR_1, listw = dist.W2)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W2
##
## Moran I statistic standard deviate = 15.841, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.4807656728 -0.0152903082 0.0009806596Aquí, debemos observar de que ambas matrices espaciales indican la existencia de autocorrelación espacial en los errores. Esto podría indicarnos que OLS no es el estimador indicado, ya que debemos considerar la existencia de autocorrelación espacial presente en los errores.
Sin embargo, dado que el índice de moran es un contraste de hipótesis abierto, solo nos indica la presencia de autocorrelación espacial pero no nos indica que modelo espacial se ajusta de mejor manera a la estructura de los datos.
Por tanto, el siguiente paso será estimar los contrastes de hipótesis basados en los multiplicadores de lagrange (lagrange multiplier), para determinar que modelo espacial se ajusta de mejor forma a los datos:
lm.LMtests(model_SAR_1, listw = dist.W1, test = c("LMerr", "LMlag","RLMerr", "RLMlag", "SARMA"))##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W1
##
## LMerr = 140.99, df = 1, p-value < 2.2e-16
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W1
##
## LMlag = 37.785, df = 1, p-value = 7.899e-10
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W1
##
## RLMerr = 105.32, df = 1, p-value < 2.2e-16
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W1
##
## RLMlag = 2.1177, df = 1, p-value = 0.1456
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W1
##
## SARMA = 143.11, df = 2, p-value < 2.2e-16lm.LMtests(model_SAR_1, listw = dist.W2, test = c("LMerr", "LMlag","RLMerr", "RLMlag", "SARMA"))##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W2
##
## LMerr = 226.4, df = 1, p-value < 2.2e-16
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W2
##
## LMlag = 146.77, df = 1, p-value < 2.2e-16
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W2
##
## RLMerr = 83.444, df = 1, p-value < 2.2e-16
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W2
##
## RLMlag = 3.8093, df = 1, p-value = 0.05097
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = MEDV ~ CRIM + RM + INDUS + NOX + AGE + DIS +
## RAD + PTRATIO + B + LSTAT + TAX, data = data2)
## weights: dist.W2
##
## SARMA = 230.21, df = 2, p-value < 2.2e-16Recordar que el contraste de hipótesis se plantea de la siguiente forma:
\[ H_{0}=\text{Modelo Restringido} \\ H_{1}=\text{Modelo Sin Restringir} \] Donde el modelo restringido corresponde a un modelo OLS, donde \(\rho\) o \(\lambda\) están restringidos a cero. Mientras que la hipótesis alternativa puede ser cualquier modelo espacial. En este caso, la hipótesis alternativa puede ser el modelo SAR o SEM, según sea el caso.
Según el contraste de hipótesis, en las versiones regulares del LM-test, el modelo SAR y SEM al parecer serían modelos que se ajustan de buena manera a la estructura de datos. Sin embargo, esta conclusión no nos indica que modelo es el mejor. Siguiendo el procedimiento establecido por Anselin (1994), para responder la pregunta, debemos analizar las versiones robustas de ambos contrastes de hipótesis.
Siguiendo con ambos LM-test robustos, podemos encontrar algunas diferencias de acuerdo a la matriz de peso espacial usada. Por ejemplo, usando la matriz basada en distancias, el LM test indica que el modelo del error espacial (SEM) es el que mejor se ajusta a la estructura de los datos. Por otro lado, usando la matriz de peso espacial que viene por defecto, podemos ver que ambos LM test son significativos. Sin embargo, el LM test para el modelo del error espacial (SEM) es más significativo comparado con el modelo del rezago espacial (SAR). Por tanto, podemos llegamos a la misma conclusión obtenida usando la primera matriz de peso espacial. De todas formas, debemos tener cuidado a la hora de elegir este modelo, ya que también podría indicar problemas de especificación del modelo estimado.
Modelo del Error Espacial (SEM)
Estimamos el modelo anterior, considerando los resultados obtenidos de acuerdo al contraste de hipótesis basado en los multiplicadores de lagrange.
Estimaciones Vía ML
boston.sem1<-errorsarlm(MEDV ~ CRIM+RM+INDUS+NOX+AGE+DIS+RAD+PTRATIO+B+LSTAT+TAX, listw=dist.W1, data=data2)
summary(boston.sem1)##
## Call:
## spatialreg::errorsarlm(formula = formula, data = data, listw = listw,
## na.action = na.action, Durbin = Durbin, etype = etype, method = method,
## quiet = quiet, zero.policy = zero.policy, interval = interval,
## tol.solve = tol.solve, trs = trs, control = control)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.36829 -2.65446 -0.78679 1.65060 25.69520
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 31.2946675 7.5098607 4.1671 3.084e-05
## CRIM -0.1183595 0.0320502 -3.6929 0.0002217
## RM 4.2836840 0.4034602 10.6174 < 2.2e-16
## INDUS -0.0042883 0.0620593 -0.0691 0.9449103
## NOX -17.6122695 4.0668263 -4.3307 1.486e-05
## AGE -0.0104953 0.0132409 -0.7926 0.4279859
## DIS -1.1712087 0.1986458 -5.8960 3.725e-09
## RAD 0.3242508 0.0697805 4.6467 3.372e-06
## PTRATIO -0.9966088 0.1305569 -7.6335 2.287e-14
## B 0.0094697 0.0028402 3.3341 0.0008556
## LSTAT -0.4951812 0.0508188 -9.7441 < 2.2e-16
## TAX -0.0111954 0.0036249 -3.0885 0.0020117
##
## Lambda: 0.96026, LR test value: 43.784, p-value: 3.6673e-11
## Asymptotic standard error: 0.026946
## z-value: 35.637, p-value: < 2.22e-16
## Wald statistic: 1270, p-value: < 2.22e-16
##
## Log likelihood: -1487.721 for error model
## ML residual variance (sigma squared): 20.712, (sigma: 4.5511)
## Number of observations: 506
## Number of parameters estimated: 14
## AIC: 3003.4, (AIC for lm: 3045.2)boston.sem2<-errorsarlm(MEDV ~ CRIM+RM+INDUS+NOX+AGE+DIS+RAD+PTRATIO+B+LSTAT+TAX, listw=dist.W2, data=data2)
summary(boston.sem2)##
## Call:
## spatialreg::errorsarlm(formula = formula, data = data, listw = listw,
## na.action = na.action, Durbin = Durbin, etype = etype, method = method,
## quiet = quiet, zero.policy = zero.policy, interval = interval,
## tol.solve = tol.solve, trs = trs, control = control)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.47676 -1.94487 -0.57867 1.29241 22.00163
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 28.7711397 5.3972049 5.3307 9.781e-08
## CRIM -0.0849750 0.0272788 -3.1151 0.0018391
## RM 4.2669725 0.3521893 12.1156 < 2.2e-16
## INDUS -0.0562610 0.0736717 -0.7637 0.4450628
## NOX -14.3128068 5.4079952 -2.6466 0.0081305
## AGE -0.0377138 0.0131799 -2.8615 0.0042170
## DIS -1.3421754 0.2720612 -4.9334 8.083e-07
## RAD 0.2613614 0.0768493 3.4010 0.0006715
## PTRATIO -0.7444100 0.1476810 -5.0407 4.639e-07
## B 0.0119054 0.0029704 4.0081 6.122e-05
## LSTAT -0.3629340 0.0509620 -7.1217 1.066e-12
## TAX -0.0121817 0.0036054 -3.3787 0.0007283
##
## Lambda: 0.62698, LR test value: 191.59, p-value: < 2.22e-16
## Asymptotic standard error: 0.037754
## z-value: 16.607, p-value: < 2.22e-16
## Wald statistic: 275.79, p-value: < 2.22e-16
##
## Log likelihood: -1413.816 for error model
## ML residual variance (sigma squared): 13.843, (sigma: 3.7206)
## Number of observations: 506
## Number of parameters estimated: 14
## AIC: 2855.6, (AIC for lm: 3045.2)Estimación Vía GLS
boston.sem3<-GMerrorsar(MEDV ~ CRIM+RM+INDUS+NOX+AGE+DIS+RAD+PTRATIO+B+LSTAT+TAX, listw=dist.W1, data=data2)
summary(boston.sem3)##
## Call:
## spatialreg::GMerrorsar(formula = formula, data = data, listw = listw,
## na.action = na.action, zero.policy = zero.policy, method = method,
## arnoldWied = arnoldWied, control = control, pars = pars,
## scaleU = scaleU, verbose = verbose, legacy = legacy, se.lambda = se.lambda,
## returnHcov = returnHcov, pWOrder = pWOrder, tol.Hcov = tol.Hcov)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.8806 -3.2507 -1.1079 1.4808 28.3380
##
## Type: GM SAR estimator
## Coefficients: (GM standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 29.3344804 5.4926990 5.3406 9.262e-08
## CRIM -0.1341545 0.0300137 -4.4698 7.830e-06
## RM 4.7032774 0.3790880 12.4068 < 2.2e-16
## INDUS -0.0205716 0.0668562 -0.3077 0.7583114
## NOX -19.4707661 5.1913811 -3.7506 0.0001764
## AGE -0.0328326 0.0136317 -2.4086 0.0160160
## DIS -1.3001899 0.2775214 -4.6850 2.799e-06
## RAD 0.2856166 0.0777420 3.6739 0.0002389
## PTRATIO -0.7819007 0.1432451 -5.4585 4.802e-08
## B 0.0105342 0.0030478 3.4563 0.0005476
## LSTAT -0.4406324 0.0513366 -8.5832 < 2.2e-16
## TAX -0.0094539 0.0035604 -2.6553 0.0079239
##
## Lambda: 3.9323 (standard error): NaN (z-value): NaN
## Residual variance (sigma squared): 17.322, (sigma: 4.162)
## GM argmin sigma squared: 17.322
## Number of observations: 506
## Number of parameters estimated: 14boston.sem4<-GMerrorsar(MEDV ~ CRIM+RM+INDUS+NOX+AGE+DIS+RAD+PTRATIO+B+LSTAT+TAX, listw=dist.W2, data=data2)
summary(boston.sem4)##
## Call:
## spatialreg::GMerrorsar(formula = formula, data = data, listw = listw,
## na.action = na.action, zero.policy = zero.policy, method = method,
## arnoldWied = arnoldWied, control = control, pars = pars,
## scaleU = scaleU, verbose = verbose, legacy = legacy, se.lambda = se.lambda,
## returnHcov = returnHcov, pWOrder = pWOrder, tol.Hcov = tol.Hcov)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7672 -2.7788 -0.8284 1.6586 29.4175
##
## Type: GM SAR estimator
## Coefficients: (GM standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 29.4046169 5.4211169 5.4241 5.825e-08
## CRIM -0.0877069 0.0281650 -3.1140 0.0018454
## RM 4.3056030 0.3640108 11.8282 < 2.2e-16
## INDUS -0.0467449 0.0731892 -0.6387 0.5230275
## NOX -14.5790050 5.2451106 -2.7795 0.0054436
## AGE -0.0349967 0.0134071 -2.6103 0.0090460
## DIS -1.3381876 0.2610185 -5.1268 2.947e-07
## RAD 0.2652296 0.0760499 3.4876 0.0004874
## PTRATIO -0.7925132 0.1470923 -5.3879 7.130e-08
## B 0.0117077 0.0029941 3.9103 9.220e-05
## LSTAT -0.3760408 0.0517244 -7.2701 3.593e-13
## TAX -0.0119977 0.0036738 -3.2658 0.0010916
##
## Lambda: 0.57168 (standard error): 0.071257 (z-value): 8.0228
## Residual variance (sigma squared): 14.786, (sigma: 3.8452)
## GM argmin sigma squared: 15.18
## Number of observations: 506
## Number of parameters estimated: 14Nuevamente, análisis de autocorrelación de los residuos
residuals_sem_1<-residuals(boston.sem1)
moran.test(residuals_sem_1, listw = dist.W1, randomisation = FALSE)##
## Moran I test under normality
##
## data: residuals_sem_1
## weights: dist.W1
##
## Moran I statistic standard deviate = 13.349, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 4.097824e-02 -1.980198e-03 1.035628e-05residuals_sem_2<-residuals(boston.sem2)
moran.test(residuals_sem_2, listw = dist.W2, randomisation = FALSE)##
## Moran I test under normality
##
## data: residuals_sem_2
## weights: dist.W2
##
## Moran I statistic standard deviate = -2.1639, p-value = 0.9848
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.070843610 -0.001980198 0.001012720Caso 2: Los determinantes del crimen en Columbus (OH)
data("columbus")
ls()## [1] "bbs" "boston.c"
## [3] "boston.sem1" "boston.sem2"
## [5] "boston.sem3" "boston.sem4"
## [7] "boston.soi" "boston.utm"
## [9] "col.gal.nb" "colours"
## [11] "columbus" "coords"
## [13] "data2" "dist.mat"
## [15] "dist.mat.inv" "dist.W1"
## [17] "dist.W2" "map"
## [19] "model_SAR_1" "polys"
## [21] "residuals_model_SAR_1" "residuals_sem_1"
## [23] "residuals_sem_2"str("columbus")## chr "columbus"COL.shape<-readShapePoly('columbus.shp')
class(COL.shape)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"plot(COL.shape)
attach(COL.shape@data)COL_OLS<-lm(CRIME ~ HOVAL+INC, data=columbus)
summary(COL_OLS)##
## Call:
## lm(formula = CRIME ~ HOVAL + INC, 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 ***
## HOVAL -0.2739 0.1032 -2.654 0.0109 *
## INC -1.5973 0.3341 -4.780 1.83e-05 ***
## ---
## 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-09Modelo SARAR
W.COL<-nb2listw(col.gal.nb)
summary(W.COL)## 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
## Link number distribution:
##
## 2 3 4 5 6 7 8 9 10
## 7 7 13 4 9 6 1 1 1
## 7 least connected regions:
## 1005 1008 1045 1047 1049 1048 1015 with 2 links
## 1 most connected region:
## 1017 with 10 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 49 2401 49 23.48489 204.6687COL_SARAR1<-sacsarlm(CRIME ~ HOVAL+INC, data=columbus, listw=W.COL, Durbin = FALSE)
summary(COL_SARAR1)##
## Call:spatialreg::sacsarlm(formula = formula, data = data, listw = listw,
## listw2 = listw2, na.action = na.action, Durbin = Durbin,
## type = type, method = method, quiet = quiet, zero.policy = zero.policy,
## tol.solve = tol.solve, llprof = llprof, interval1 = interval1,
## interval2 = interval2, trs1 = trs1, trs2 = trs2, control = control)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.1121 -4.6324 -0.3040 7.0306 24.6929
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 49.051431 10.054986 4.8783 1.07e-06
## HOVAL -0.283114 0.091526 -3.0933 0.001980
## INC -1.068781 0.332839 -3.2111 0.001322
##
## Rho: 0.35326
## Asymptotic standard error: 0.19669
## z-value: 1.796, p-value: 0.072494
## Lambda: 0.13199
## Asymptotic standard error: 0.29905
## z-value: 0.44138, p-value: 0.65894
##
## LR test value: 8.6082, p-value: 0.013513
##
## Log likelihood: -183.0731 for sac model
## ML residual variance (sigma squared): 99.423, (sigma: 9.9711)
## Number of observations: 49
## Number of parameters estimated: 6
## AIC: 378.15, (AIC for lm: 382.75)Sin embargo, un problema de la estimación es que considera que los errores se distribuyen de forma homocedástica. De acuerdo a la estructura de los datos, este supuesto puede ser no realista, ya que las unidades espaciales difieren en términos de tamaño. Por tanto, es muy difícil asumir la hipótesis de que la varianza de los residuos es constante. Esto también implica que nuestras estimaciones ya no sean eficientes.
Para corregir esto, usamos el Generalized Spatial Two-Stage Least Square:
COL_SARAR2<-gstslshet(CRIME ~ HOVAL+INC, data=columbus, listw=W.COL)
summary(COL_SARAR2)##
## Generalized stsls
##
## Call:
## gstslshet(formula = CRIME ~ HOVAL + INC, data = columbus, listw = W.COL)
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -37.468 -5.520 -0.398 -0.005 6.260 24.142
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 44.124087 7.506698 5.8780 4.153e-09 ***
## HOVAL -0.275573 0.176910 -1.5577 0.119305
## INC -0.987477 0.460184 -2.1458 0.031886 *
## lambda 0.452910 0.143270 3.1612 0.001571 **
## rho 0.064822 0.305091 0.2125 0.831742
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Wald test that rho and lambda are both zero:
## Statistics: 3.5907 p-val: 0.058105También, podemos usar una estimación no paramétrica de la matriz de varianza-covarianza, usando HAC (heteroscedastic autocorrelation consistent).
data(coldis)
COL_SARAR3<-stslshac(CRIME ~ HOVAL+INC, data=columbus, listw=W.COL, distance = coldis, type="Epanechnikov")
summary(COL_SARAR3)##
## Stsls with Spatial HAC standard errors
##
## Call:
## stslshac(formula = CRIME ~ HOVAL + INC, data = columbus, listw = W.COL,
## distance = coldis, type = "Epanechnikov")
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -37.83067 -5.43232 -0.44058 6.26723 24.21600
##
## Coefficients:
## Estimate SHAC St.Er. t-value Pr(>|t|)
## Wy 0.45464 0.17250 2.6356 0.008398 **
## (Intercept) 44.11639 8.18951 5.3869 7.167e-08 ***
## HOVAL -0.26950 0.18030 -1.4948 0.134979
## INC -1.00772 0.53790 -1.8734 0.061009 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1