Parte 1

Describir autocorrelación espacial, autocorrelación espacial positiva, y autocorrelación espacial negativa.

  • Autocorrelación espacial: La autocorrelación espacial es identificar la presencia de variación espacial de una variable en una zona. r. Abdishakur Hassan 2022
  • Autocorrelación espacial positiva: Es cuando en una zona de vecinos que están cerca, hay valores muy similares/parecidos. Está se divide en alta-alta y baja-baja.
  • Autocorrelación espacial negativa: Es cuando en una zona de vecinos que están cerca, los valores son contrarios. Está se divide en alta-baja y baja-alta.

Describir al menos 3-5 diferencias entre análisis exploratorio de datos (EDA) y análisis exploratorio espacial de datos (ESDA).

  • EDA no toma en cuenta el factor espacial, ESDA sí lo hace.
  • EDA hace análisis globales de los datos, ESDA puede hacer además, análisis locales de los datos.
  • EDA no permite identificar patrones de comportamiento en el espacio, ESDA sí lo hace.

Describir los conceptos de autocorrelación espacial y no estacionariedad en un contexto de análisis espacial de datos.

  • En un contexto de análisis espacial de datos, la autocorrelación espacial permite identificar las relaciones entre vecinos en un territorio, resultando útil para hacer clusters e idetificar patrones de comportamiento. En el caso de no estacionariedad, se refiere a que la variable que cuenta con esto, tiene variabilidad en los diversos vecinos del territorio, es decir, es cambiante en el espacio.

Describir al menos 3-5 diferencias entre la estimación de modelo de regresión no espacial, espacial global, y espacial local.

  • Un modelo no espacial no cuenta con especificaciones que le permiten tomar la autocorrelación espacial en el cálculo de sus predicciones, mientras que uno global sí. *Un modelo espacial global solamente desarrolla un modelo para todos los datos, mientras que uno local, como el GWR, desarrolla un modelo para cada locación de la base de datos. r. Columbia 2023.
  • Por lo tanto, un modelo no global sigue siendo útil cuando los datos no tienen autocorrelación espacial, más si la tienen un modelo espacial global puede ser útil para predecir y uno local para interpretar. r. Columbia 2023.

Describir cómo el proceso de análisis espacial de datos puede mejorar las herramientas de Descriptive Analytics y Predictive Analytics en un contexto de Inteligencia de Negocios.

  • En un contexto de negocios, el ESDA es sumamente útil para entender el comportamiento de factores de interés en un territorio, como lo son clientes, proveedores, competidores, etc. Añade mucho a la inteligencia de negocios pues te permite hacer análisis más especializados para sectores poblacionales o de otro tipo que se pueden llegar a comportar según la primera ley de geografía. Ejemplos de esto son UBER, RAPPI, aerolíneas, etc.

Parte 2

Paso 1 - Análisis Espacial de Datos

Iniciamos con un análisis espacial de datos normales, buscando que las variables en la base datos tengan una distribución normal.

Gráfica: Correlación entre variables

#1. Correlación entre variables
library(dlookr)
## 
## Attaching package: 'dlookr'
## The following object is masked from 'package:tidyr':
## 
##     extract
## The following object is masked from 'package:base':
## 
##     transform
corr_columbus <- columbus %>% select(AREA, PERIMETER, HOVAL, INC, CRIME, OPEN, PLUMB, DISCBD)
correlate(corr_columbus, AREA, PERIMETER, HOVAL, INC, CRIME, OPEN, PLUMB, DISCBD) %>% plot()

Hoval, la cual es la variable dependiente presenta correlación positiva con Income y Distancia, y negativa con Crimen.

Gráficas de Normalidad para todas las variables

plot_normality(columbus, AREA)

plot_normality(columbus, PERIMETER)

plot_normality(columbus, HOVAL)

plot_normality(columbus, INC)

plot_normality(columbus, CRIME)

plot_normality(columbus, OPEN)

plot_normality(columbus, PLUMB)

plot_normality(columbus, DISCBD)

Después de observar las gráficas, podemos determinar que PLUMB y OPEN no siguen una distribución normal, por lo que al hacer los modelos será conveniente transformar las variables con log y sqrt respectivamente. Además la variable dependiente HOVAL cuenta con una inclinación hacia la izquierda, por lo que se podría llegar a transformar con log para mejorar modelos.

Gráfica: Comportamiento de Distancia e Ingreso en el Precio

library(ggplot2)
ggplot(columbus, aes(x = DISCBD, y = INC, size = HOVAL)) +  geom_point() + scale_x_log10()

Siendo el Precio de las casas el tamaño, se determina que entre mayor ingreso y mayor distancia al centro, mayor es el precio de las casas.

swm_queen <- poly2nb(columbus_shp, queen = TRUE)
swm_rook  <- poly2nb(columbus_shp, queen = FALSE)

## Standardized Queen & Rook Spatial Connectivity Matrices
rswm_queen <- nb2listw(swm_queen, style = "W", zero.policy = TRUE)
rswm_rook  <- nb2listw(swm_rook, style = "W", zero.policy = TRUE)

columbus_shp$sp_HOVAL<-lag.listw(rswm_queen,columbus_shp$HOVAL,zero.policy=TRUE)

Gráfica: Precios en Columbus

qtm(columbus_shp, "HOVAL")

Confirmando la gráfica anterior, entre más alejados del centro, los precios son mayores.

Gráfica: Precios en Columbus con spatial lag

qtm(columbus_shp, "sp_HOVAL")

Al igual que el anterior mapa, entre más lejos del mapa, mayores los precios, pero al ser este un mapa con lag, los precios son menores.

Paso 2 - Modelos de predicción

Análisis de Autocorrelación Espacial Global: HOVAL

moran.test(columbus_shp$HOVAL, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$HOVAL  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 2.2071, p-value = 0.01365
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.180093114      -0.020833333       0.008287783

Análisis de Autocorrelación Espacial Global: CRIME

moran.test(columbus_shp$CRIME, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$CRIME  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 5.5894, p-value = 1.139e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.500188557      -0.020833333       0.008689289

Análisis de Autocorrelación Espacial Global: INCOME

moran.test(columbus_shp$INC, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$INC  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 4.7645, p-value = 9.467e-07
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.415628778      -0.020833333       0.008391926

Análisis de Autocorrelación Espacial Global: DISCBD (distancia al centro)

moran.test(columbus_shp$DISCBD, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$DISCBD  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 8.7904, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.800562617      -0.020833333       0.008731396

Todas las variables analizadas tienen autocorrelación espacial positiva estádisticamente significativa.

Local Spatial Analysis

Ahora iniciamos con el análisis local, utilizando la connectivity matrix para luego crear las gráficas para hacer el análisis local.

## Reading layer `columbus' from data source 
##   `C:\Users\art191127\AppData\Local\R\win-library\4.2\spdep\etc\shapes\columbus.shp' 
##   using driver `ESRI Shapefile'
## 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:           NA
columbus_map_sf<-read_sf(system.file("etc/shapes/columbus.shp",package="spdep"))
queen_w<-queen_weights(columbus_map_sf)
columbus_map_centroid<-coordinates(columbus.tr) 
columbus_map.linkW<-nb2listw(columbus_nb, style="W")   
plot(columbus.tr,border="blue",axes=FALSE,las=1, main="Columbus Spatial Connectivity Matrix")
plot(columbus.tr,col="grey",border=grey(0.9),axes=T,add=T) 
plot(columbus_map.linkW,coords=columbus_map_centroid,pch=19,cex=0.1,col="red",add=T)

Mapa: Local Clusters HOVAL

lisa_MEDV<-local_moran(queen_w, columbus_map_sf["HOVAL"]) 

columbus_map_sf$HOVAL_MEDV<-as.factor(lisa_MEDV$GetClusterIndicators())
levels(columbus_map_sf$HOVAL_MEDV)<-lisa_MEDV$GetLabels()

ggplot(data=columbus_map_sf) +
  geom_sf(aes(fill=HOVAL_MEDV)) + 
  ggtitle(label = "Median Values of Housing Units", subtitle = "Columbus Housing Market")

Para la variable HOVAL no notamos más que unos pocos municipios con los 4 tipos de autocorrelación, sin embargo sin un patrón claro.

Mapa: Local Clusters INCOME

lisa_MEDV<-local_moran(queen_w, columbus_map_sf["INC"]) 

columbus_map_sf$INC_MEDV<-as.factor(lisa_MEDV$GetClusterIndicators())
levels(columbus_map_sf$INC_MEDV)<-lisa_MEDV$GetLabels()

ggplot(data=columbus_map_sf) +
  geom_sf(aes(fill=INC_MEDV)) + 
  ggtitle(label = "Median Income", subtitle = "Columbus Housing Market")

En este caso, se pueden notar 2 agrupaciones de autocorrelación espacial positiva, más de diferentes tipos, siendo la agrupación del noroeste de valores altos, mientras que la del sureste valores bajos.

Mapa: Local Clusters CRIME

lisa_MEDV<-local_moran(queen_w, columbus_map_sf["CRIME"]) 

columbus_map_sf$CRIME_MEDV<-as.factor(lisa_MEDV$GetClusterIndicators())
levels(columbus_map_sf$CRIME_MEDV)<-lisa_MEDV$GetLabels()

ggplot(data=columbus_map_sf) +
  geom_sf(aes(fill=CRIME_MEDV)) + 
  ggtitle(label = "Median Crime", subtitle = "Columbus Housing Market")

En el caso de crimen, podemos ver como también cuenta con agrupaciones de municipios con correlación espacial positiva, siendo en el centro valores altos, lo cual era de esperarse en una ciudad urbana, mientras que en el noreste, del lado donde también había alto ingreso en la anterior gráfica, hay valores bajos de crimen.

Mapa: Local Clusters DISTANCIA

lisa_MEDV<-local_moran(queen_w, columbus_map_sf["DISCBD"]) 

columbus_map_sf$DISCBD_MEDV<-as.factor(lisa_MEDV$GetClusterIndicators())
levels(columbus_map_sf$DISCBD_MEDV)<-lisa_MEDV$GetLabels()

ggplot(data=columbus_map_sf) +
  geom_sf(aes(fill=DISCBD_MEDV)) + 
  ggtitle(label = "Median Distance", subtitle = "Columbus Housing Market")

En este caso, las correlaciones son bastantes claras, pues al ser la variable de distancia al centro es normal que la auotocorrelación espacial positiva de valores bajos este en el centro y la de valores altos en las orillas.

Mapa: Local Clusters PLUMB

lisa_MEDV<-local_moran(queen_w, columbus_map_sf["PLUMB"]) 

columbus_map_sf$PLUMB_MEDV<-as.factor(lisa_MEDV$GetClusterIndicators())
levels(columbus_map_sf$PLUMB_MEDV)<-lisa_MEDV$GetLabels()

ggplot(data=columbus_map_sf) +
  geom_sf(aes(fill=PLUMB_MEDV)) + 
  ggtitle(label = "Median Plumb", subtitle = "Columbus Housing Market")

En el caso del porcentaje de casas sin tuberías, notamos un patrón similar al de las anteriores gráficas, donde en el centro hay correlación espacial positiva alta, por lo que hay más casas sin tuberías, mientras que en el sureste hay baja, por lo que hay más tuberías.

Paso 3 - Modelos de Predicción

Regresión lineal simple

Para los modelos de predicción inicio con una regresión simple.

lm_sp <- lm(HOVAL ~ AREA + PERIMETER + INC + CRIME + OPEN + PLUMB + DISCBD, data = columbus)
summary(lm_sp)
## 
## Call:
## lm(formula = HOVAL ~ AREA + PERIMETER + INC + CRIME + OPEN + 
##     PLUMB + DISCBD, data = columbus)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.184  -7.958  -2.442   4.922  54.291 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  34.3767    17.1915   2.000  0.05220 . 
## AREA          2.7221     1.2578   2.164  0.03632 * 
## PERIMETER   -13.8862     7.3959  -1.878  0.06757 . 
## INC           0.2750     0.5112   0.538  0.59346   
## CRIME        -0.3513     0.2168  -1.621  0.11276   
## OPEN          0.5714     0.4494   1.272  0.21069   
## PLUMB         1.8889     0.6674   2.830  0.00717 **
## DISCBD        5.3576     2.4311   2.204  0.03321 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.86 on 41 degrees of freedom
## Multiple R-squared:  0.5191, Adjusted R-squared:  0.4369 
## F-statistic: 6.321 on 7 and 41 DF,  p-value: 4.618e-05
sqrt(mean((columbus$HOVAL - lm_sp$fitted.values)^2))
## [1] 12.67496
AIC(lm_sp)
## [1] 405.9396

Regresión lineal simple mod

Ahora incluimos a la regresión lineal simple los hallazgos realizados en el EDA.

lm_sp_auto <- lm(log(HOVAL) ~ AREA + PERIMETER + INC + CRIME + OPEN^2 + log(PLUMB) + DISCBD, data = columbus)
summary(lm_sp_auto)
## 
## Call:
## lm(formula = log(HOVAL) ~ AREA + PERIMETER + INC + CRIME + OPEN^2 + 
##     log(PLUMB) + DISCBD, data = columbus)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37428 -0.19566 -0.00598  0.13564  0.94605 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.066010   0.378133   8.108 4.71e-10 ***
## AREA         0.053575   0.025402   2.109 0.041090 *  
## PERIMETER   -0.293780   0.149767  -1.962 0.056625 .  
## INC          0.032541   0.011467   2.838 0.007031 ** 
## CRIME       -0.008801   0.004434  -1.985 0.053868 .  
## OPEN         0.012236   0.009255   1.322 0.193497    
## log(PLUMB)   0.239541   0.062413   3.838 0.000421 ***
## DISCBD       0.175039   0.053884   3.248 0.002319 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2835 on 41 degrees of freedom
## Multiple R-squared:  0.6329, Adjusted R-squared:  0.5702 
## F-statistic:  10.1 on 7 and 41 DF,  p-value: 2.87e-07
sqrt(mean((columbus$HOVAL - exp(lm_sp_auto$fitted.values))^2))
## [1] 12.2514
AIC(lm_sp_auto)
## [1] 24.79238

Nos da un mejor RMSE e AIC, por lo que se repetiran las transformaciones a HOVAL, OPEN y PLUMB en los siguientes modelos

Moran Test para el Modelo No Espacial

moran.test(lm_sp_auto$residuals, rswm_queen)
## 
##  Moran I test under randomisation
## 
## data:  lm_sp_auto$residuals  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 0.55256, p-value = 0.2903
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.029113040      -0.020833333       0.008170386

Dado que el p-value no es significativo, no hay concentraciones de Clusters en los resiudales del modelo no espacial, por lo tanto, el modelo está bien especificado.

Spatial Models

DURBIN

spatial_durbin <- lagsarlm(log(HOVAL) ~ AREA + PERIMETER + INC + CRIME + OPEN^2 + log(PLUMB) + DISCBD, data = columbus, rswm_queen, type = 'mixed', Durbin = TRUE)
summary(spatial_durbin)
## 
## Call:lagsarlm(formula = log(HOVAL) ~ AREA + PERIMETER + INC + CRIME + 
##     OPEN^2 + log(PLUMB) + DISCBD, data = columbus, listw = rswm_queen, 
##     Durbin = TRUE, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.382796 -0.175914 -0.026944  0.125061  0.754857 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     2.6746983  0.9251331  2.8911 0.0038384
## AREA            0.0552639  0.0254916  2.1679 0.0301643
## PERIMETER      -0.3697178  0.1547896 -2.3885 0.0169165
## INC             0.0347711  0.0095507  3.6407 0.0002719
## CRIME          -0.0095583  0.0036701 -2.6044 0.0092038
## OPEN            0.0151493  0.0082976  1.8257 0.0678900
## log(PLUMB)      0.1941515  0.0607883  3.1939 0.0014037
## DISCBD          0.4360162  0.1188223  3.6695 0.0002430
## lag.AREA       -0.0114703  0.0504131 -0.2275 0.8200151
## lag.PERIMETER   0.1666277  0.3383121  0.4925 0.6223471
## lag.INC         0.0040684  0.0190855  0.2132 0.8311954
## lag.CRIME       0.0043599  0.0095626  0.4559 0.6484384
## lag.OPEN        0.0304823  0.0241902  1.2601 0.2076297
## lag.log(PLUMB)  0.0718721  0.1297700  0.5538 0.5796870
## lag.DISCBD     -0.3184185  0.1644086 -1.9368 0.0527758
## 
## Rho: 0.029334, LR test value: 0.018026, p-value: 0.8932
## Asymptotic standard error: 0.20191
##     z-value: 0.14528, p-value: 0.88449
## Wald statistic: 0.021106, p-value: 0.88449
## 
## Log likelihood: 3.013575 for mixed model
## ML residual variance (sigma squared): 0.051764, (sigma: 0.22752)
## Number of observations: 49 
## Number of parameters estimated: 17 
## AIC: 27.973, (AIC for lm: 25.991)
## LM test for residual autocorrelation
## test value: 1.5004, p-value: 0.22061
sqrt(mean((columbus$HOVAL - exp(spatial_durbin$fitted.values))^2))
## [1] 11.05289
moran.test(spatial_durbin$residuals, rswm_queen)
## 
##  Moran I test under randomisation
## 
## data:  spatial_durbin$residuals  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = -0.072712, p-value = 0.529
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.027474370      -0.020833333       0.008341805

Nuevamente, el Moran Test para los residuales del modelo dio un p-value no significativo, por lo que no hay autocorrelación espacial en los residuales, por lo que el modelo está bien especificado.

GWR

# Kernel bandwidth
bw1 <- bw.gwr(HOVAL ~ AREA + PERIMETER + INC + CRIME + OPEN + PLUMB + DISCBD, 
             approach = "AIC", adaptive = T, data=columbus_shp)
## Adaptive bandwidth (number of nearest neighbours): 37 AICc value: 435.0207 
## Adaptive bandwidth (number of nearest neighbours): 31 AICc value: 451.2796 
## Adaptive bandwidth (number of nearest neighbours): 42 AICc value: 425.3363 
## Adaptive bandwidth (number of nearest neighbours): 44 AICc value: 423.0324 
## Adaptive bandwidth (number of nearest neighbours): 46 AICc value: 420.7616 
## Adaptive bandwidth (number of nearest neighbours): 47 AICc value: 420.1263 
## Adaptive bandwidth (number of nearest neighbours): 48 AICc value: 418.7132 
## Adaptive bandwidth (number of nearest neighbours): 48 AICc value: 418.7132
m.gwr <- gwr.basic(log(HOVAL) ~ AREA + PERIMETER + INC + CRIME + OPEN^2 + log(PLUMB) + DISCBD, adaptive = T, data = columbus_shp, bw = bw1) 
m.gwr
##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2023-06-13 20:56:01 
##    Call:
##    gwr.basic(formula = log(HOVAL) ~ AREA + PERIMETER + INC + CRIME + 
##     OPEN^2 + log(PLUMB) + DISCBD, data = columbus_shp, bw = bw1, 
##     adaptive = T)
## 
##    Dependent (y) variable:  HOVAL
##    Independent variables:  AREA PERIMETER INC CRIME OPEN PLUMB DISCBD
##    Number of data points: 49
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42664 -0.18131 -0.01919  0.13215  0.90479 
## 
##    Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
##    (Intercept)  3.457573   0.459983   7.517  3.1e-09 ***
##    AREA         2.070143   1.250545   1.655 0.105481    
##    PERIMETER   -0.363342   0.224519  -1.618 0.113263    
##    INC          0.025569   0.011986   2.133 0.038934 *  
##    CRIME       -0.013356   0.004665  -2.863 0.006577 ** 
##    OPEN         0.014023   0.009522   1.473 0.148444    
##    log(PLUMB)   0.246176   0.064681   3.806 0.000463 ***
##    DISCBD       0.157577   0.053799   2.929 0.005532 ** 
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 0.289 on 41 degrees of freedom
##    Multiple R-squared: 0.6185
##    Adjusted R-squared: 0.5534 
##    F-statistic: 9.497 on 7 and 41 DF,  p-value: 5.969e-07 
##    ***Extra Diagnostic information
##    Residual sum of squares: 3.424192
##    Sigma(hat): 0.2699169
##    AIC:  26.6692
##    AICc:  31.28458
##    BIC:  29.72196
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: bisquare 
##    Adaptive bandwidth: 48 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: Euclidean distance metric is used.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                     Min.    1st Qu.     Median    3rd Qu.    Max.
##    Intercept   3.0030257  3.2550108  3.3798423  3.5085636  3.5998
##    AREA        0.2073765  1.0300755  1.6081441  1.8822719  2.3409
##    PERIMETER  -0.4411890 -0.3459993 -0.2943410 -0.2171521 -0.0967
##    INC         0.0259320  0.0284559  0.0299024  0.0322012  0.0419
##    CRIME      -0.0176501 -0.0155647 -0.0131996 -0.0110549 -0.0067
##    OPEN        0.0045670  0.0089048  0.0127941  0.0156084  0.0182
##    log(PLUMB)  0.1569920  0.2083373  0.2353698  0.2652081  0.3241
##    DISCBD      0.0828695  0.1191438  0.1431939  0.1658951  0.2139
##    ************************Diagnostic information*************************
##    Number of data points: 49 
##    Effective number of parameters (2trace(S) - trace(S'S)): 15.23485 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 33.76515 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 37.02456 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 10.33914 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): -1.701527 
##    Residual sum of squares: 2.729495 
##    R-square value:  0.6959227 
##    Adjusted R-square value:  0.5545355 
## 
##    ***********************************************************************
##    Program stops at: 2023-06-13 20:56:01
sqrt(mean((columbus_shp$HOVAL - m.gwr$fitted.values)^2))
## [1] NaN
gwr_data <- m.gwr$SDF
moran.test(exp(gwr_data$residual), rswm_queen)
## 
##  Moran I test under randomisation
## 
## data:  exp(gwr_data$residual)  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 0.17362, p-value = 0.4311
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.006429839      -0.020833333       0.006882518

Nuevamente, el p-value sale no significativo, por lo que no hay autocorrelación espacial y el modelo está bien especificado.

Paso 4

Análisis de Multicolinealidad

VIF(lm_sp_auto)
##       AREA  PERIMETER        INC      CRIME       OPEN log(PLUMB)     DISCBD 
##   7.607459   7.340353   2.554205   3.286746   1.114715   3.676390   3.612775

Dentro de las variables del modelo, solamente Area y Perimetro cuentan con valores por arriba de 5 en esta prueba, por lo tanto, son las unicas que cuentan con multicolinealidad dada su naturaleza de ser ambas del mismo origen.

Lagrange Multiplier: LMlag

lm.LMtests(lm_sp_auto,rswm_queen,test=c("RLMlag"))
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = log(HOVAL) ~ AREA + PERIMETER + INC + CRIME +
## OPEN^2 + log(PLUMB) + DISCBD, data = columbus)
## weights: rswm_queen
## 
## RLMlag = 0.22497, df = 1, p-value = 0.6353

Dado que el p-value es no significativo (arriba de 0.05), no es necesario especificar el spatial lag de la variable dependiente en el modelo.

Lagrange Multiplier: LMerr

lm.LMtests(lm_sp_auto,rswm_queen,test=c("RLMerr"))
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = log(HOVAL) ~ AREA + PERIMETER + INC + CRIME +
## OPEN^2 + log(PLUMB) + DISCBD, data = columbus)
## weights: rswm_queen
## 
## RLMerr = 0.083756, df = 1, p-value = 0.7723

Nuevamente, el p-value de la prueba no es significativo, por lo que no es necesario especificar el spatial lag del error en el modelo.

Paso 5

Criterio de selección del modelo

Para seleccionar el mejor modelo, se utilizará AIC como criterio de selección dado que existe para todos los modelos y además es una métrica útil para encontrar un modelo no muy complejo pero tampoco muy simple.

Selección del modelo

##   model_names aic_values
## 1         OLS   24.79238
## 2         SDM   27.97300
## 3         GWR   10.33914

Como se puede ver, el modelo con menor AIC, es el modelo GWR, y por lo tanto es el modelo seleccionado como el mejor.

Principales hallazgos

  • Las variables de HOVAL, PLUMB y OPEN se comportan mejor en los modelos después de transformarlas con log y sqrt para que su distribución se asemeje más a una normal.
  • Existe autocorrelación espacial global positiva y estadísticamente significativa para todas las variables seleccionadas.
  • Existen patrones de agrupaciones de autocorrelación espacial positiva en el análisis local para las variables INCOME, DISTANCIA, CRIME Y PLUMB.
  • Estas mismas variables son las únicas estadísticamente significativas en el modelo seleccionado GWR.
  • El modelo seleccionado no cuenta con autocorrelación espacial en sus residuales, por lo que está bien especificado.
  • En las pruebas de Lagrange Multiplier, el modelo no espacial obtuvo como resultado que no hay necesidad de especificar el spatial lag, tanto para la variable dependiente como el error.
  • Se puede notar un patrón donde las casas con menor crimen, mejor tubería y mayor income en la zona están lejos del centro, donde las casas están en una situación contraría.
  • Por lo tanto, las mejores casas con los mayores precios se encuentran en las orillas de Columbus, mientras que las peores están en el centro.

Mapa: Valores predichos para precios en Columbus

gwr_sf = st_as_sf(m.gwr$SDF)
gwr_sf$y_predicted <- exp(gwr_sf$yhat)
#mapview(gwr_sf, zcol="y_predicted", col.regions=brewer.pal(5, "Oranges"))
tm_shape(gwr_sf) +
  tm_polygons(col = "yhat", palette="YlOrRd", style="quantile", n=5, title="Prices") +
   tm_layout(title= 'Housing Price',  title.position = c('right', 'top'))

qtm(columbus_shp, "HOVAL")

Al comparar los mapas, tanto de los valores predichos como el de los valores reales, podemos ver como en general, tienen un comportamiento similar, donde los precios altos de las casas se encuentran a las orillas, mientras que el centro tiene precios más bajos. Siguen un patrón muy similar, sin embargo, los valores predichos tienen cambios más bruscos entre zonas y zonas. En general el modelo gwr hace un buen trabajo el predecir los precios dado su especificación espacial.

Referencias: