A continuación se realizará una extensa exploración y varios análisis de los datos de paquete sobre los precios de las casas (viviendas) en Columbus, OH, USA, a través de un análisis espacial para identificar y explorar factores asociados con el alto/bajo precio de éstas, mediante el uso de la herramienta de R y los siguientes paquetes y fórmulas.

## Loading required package: sp
## Checking rgeos availability: FALSE
## Please note that 'maptools' will be retired during 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-6, (SVN revision 1201)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.1.8
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
## 
## Loading required package: robustbase
## 
## Loading required package: Rcpp
## 
## Loading required package: spatialreg
## 
## Loading required package: Matrix
## 
## 
## Attaching package: 'Matrix'
## 
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## 
## 
## Attaching package: 'spatialreg'
## 
## 
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
## 
## 
## Welcome to GWmodel version 2.2-9.
## 
## Loading required package: bestglm
## 
## Loading required package: leaps
## 
## Loading required package: VGAM
## 
## Loading required package: stats4
## 
## Loading required package: splines
## 
## 
## Attaching package: 'VGAM'
## 
## 
## The following object is masked from 'package:GWmodel':
## 
##     AICc
## 
## 
## Loading required package: rpart
## 
## Loading required package: randomForest
## 
## randomForest 4.7-1.1
## 
## Type rfNews() to see new features/changes/bug fixes.
## 
## 
## Attaching package: 'randomForest'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
## 
## 
## Important regclass change from 1.3:
## All functions that had a . in the name now have an _
## all.correlations -> all_correlations, cor.demo -> cor_demo, etc.
## 
## 
## Loading required package: viridisLite
## 
## Loading required package: digest
## 
## 
## Attaching package: 'rgeoda'
## 
## 
## The following object is masked from 'package:spdep':
## 
##     skater
## 
## 
## #refugeeswelcome
## 
## 
## Attaching package: 'dlookr'
## 
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## 
## The following object is masked from 'package:base':
## 
##     transform

– Parte 1 –

Ésta puede ser útil al momento de entender distribuciones de variables (socio económicas, financieras, bienestar, etc.) dentro de un espacio (mapa).

La positiva afirma este hallazgo ya que se encuentra cuando los valores similares se encuentran cerca de unos a los otros, mientras que la negativa se muestra cuando no sucede esto y los valores opuestos se encuentran cerca en un espacio.

Distintamente, la no estacionariedad ayuda a entender eventos en donde las propiedades de las variables, cambian dependiendo de su locación/ubicación en el espacio analizado. Ésto significa que las propiedades cambian en diferentes lugares, lo cual puede afectar la precisión de los modelos. Finalmente, se puede observar como una variación de datos y tendencias en el tiempo.

Principalmente, identificar autocorrelación espacial en los residuales de un modelo de regresión estimado enseña sobreestimaciones o subest. de los coeficientes. Además, puede llevar a invalidar las pruebas de hipótesis ya que afecta la precisión y la aplicación de autocorr, espacial ante los residuales. Ésto eventualmente lleva hacia una predicción erronea, ya que impacta a los p-values, el intervalo-de-confz, y entre otros etc. Finalmente, no otorga resultados confiables o creibles.

– Parte 2 –

## Warning: shapelib support is provided by GDAL through the sf and terra packages
## among others

Análisis Espacial de Datos (EDA) - 1

Para iniciar, tenemos que cargar y explorar nuestra base de datos ‘columbus’, entender las variables y describir aquellas principales , crear estadísticos descriptivos y de dispersión, crear histogramas, encontrar correlaciones, transformar éstas y entre otras acciones

view(columbus)
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
columbus <- columbus %>% select (-1)

Correlación de variables

# Identificamos AREA, PERIMETER, HOVAL, INC, CRIME, OPEN, PLUMB, DISCBD como variables
corr_columbus <- columbus %>% select(AREA, PERIMETER, HOVAL, INC, CRIME, OPEN, PLUMB, DISCBD)
correlate(corr_columbus, AREA, PERIMETER, HOVAL, INC, CRIME, OPEN, PLUMB, DISCBD) %>% plot()

Existen algunas correlaciones importantes dentro de la matriz como area y perimetro evidentemente, Income y HOVAL, Crimen e Income, Crimen y HOVAL, PLUMB y DISCBD, entre otras.

Dist. de variables

plot_normality(columbus, AREA, PERIMETER, HOVAL, INC, CRIME, OPEN, PLUMB, DISCBD)

Podemos observar que las variables tienen una distribución normal, menos OPEN y PLUMB, por ello se tienen que realizar las transformaciones de estas dos. Inclusive, se podría realizas una transformación logarítmica a la var. dep. de HOVAL.

#Podríamos realizar la transformación de ambas aqui con: columbus$OPEN <- sqrt(columbus$OPEN) y columbus$PLUMB <- log(columbus$PLUMB), y log(HOVAL), sin embargo lo haremos en cada uno de modelos necesarios, más abajo.
ggplot(columbus, aes(x = CRIME, y = INC, size = HOVAL)) +  geom_point() + scale_x_log10()

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

Entre más alejado del CBD (centro) y mayor INCOME, mayor precio (valor) del housing (por el tamaño de los puntos negros)

Análisis Exploratorio Espacial de los Datos (ESDA) - 2

Conectivity Matrix;

# col.gal.nb <- read.gal(system.file("etc/weights/columbus.gal", package="spdep"))

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

Estandarización de queen (movimiento hacia lados) y rook (vertical, arriba abajo, ambos lados)

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) 

Hay que visualizar la presencia de clusters mediante el spatial lag de los precios de hogares

qtm(columbus_shp, "HOVAL")

qtm(columbus_shp, "sp_HOVAL")

SP - nos explica la dependencia de las zonas con respecto a los vecinos. Podemos observar como en la parte obscura, en el Este, hay más importancia/dependencia de sus municipios vecinos, mientras en la mayoría de la ciudad, los muncp. mismos no afectan el precio de las casas tanto (HOVAL). Vemos ambos la variable dep. como la correlación de vecinos queen, y su dependencia a precios limítrofes.

Identificamos autocorrelación espacial global y local de 3 variables de interés

Global Spatial Analysis

Con respecto el precio de la vivienda (HOVAL), realizamos el análisis de autocorrelación espacial global; nuestra variable dependiente.

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

Obtuvimos autocorr. esp. glob. positiva y significativa Esta variable muestra un valor bajo en autocorrelación espacial

Con respecto ala frecuencia del crimen, realizamos el análisis de autocorrelación espacial global

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

Obtuvimos autocorr. esp. glob. positiva y significativa

Con respecto al Income-ingreso, realizamos el análisis de autocorrelación espacial global

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

Obtuvimos autocorr. esp. glob. positiva y significativa

Con respecto ala distancia del centro-distrito financiero , realizamos el análisis de autocorrelación espacial global

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

Obtuvimos autocorr. esp. glob. positiva y significativa. La distancia al centro tiene alta autocorrelación, muy fuerte con el .8 de valor.

Las tres variables y aquella dependiente cuentan con autocorrelación espacial positiva.

Tres variables con Local spatial a.

Local moran - con mismas variables pero evaluando características individuales y comparadas a sus vecinos, para buscar agrupamiento/dispersión más especificas.

localMoranHOVAL <- as.data.frame(localmoran(columbus_shp$HOVAL, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit))
localMoranCRIME <- as.data.frame(localmoran(columbus_shp$CRIME, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit))
localMoranINC <- as.data.frame(localmoran(columbus_shp$INC, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit))
localMoranDist <- as.data.frame(localmoran(columbus_shp$DISCBD, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit))
map<-read_sf(system.file("etc/shapes/columbus.shp",package="spdep"))
queen_w<-queen_weights(map)

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

map$HOVAL_MEDV<-as.factor(lisa_MEDV$GetClusterIndicators())
levels(map$HOVAL_MEDV)<-lisa_MEDV$GetLabels()
columbus.tr<-sf::st_read(system.file("etc/shapes/columbus.shp",package="spdep")[1])
## Reading layer `columbus' from data source 
##   `/Library/Frameworks/R.framework/Versions/4.2/Resources/library/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.tr<-as(columbus.tr, "Spatial")
columbus_nb<-poly2nb(columbus.tr)
mapview(columbus.tr, zcol="HOVAL") 
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)  

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

Aquí los clusters significativos de High-High se ubican en el Este, mientras en el centro se estacionan aquellas Low-Low.

Explorar variables;

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

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

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

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

quadr_data <- localMoranHOVAL$quadr

Estimación de Modelos de Predicción - 3

Glob.spatial reg.

LINEAR REGRESSION ANALYSIS

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
# plot(effect("I(RM^2)",lm_sp))

Según nuestro analisis de reg. linear esp, podemos observar que Perimetro, Crimen, Plumb y Distancia centro son significativos y cuentan solo perimetro y crimen con una relación negativa. Además, el modelo lineal es significativo, sin embargo tiene un bajo r-sq con baja accuracy.

sqrt(mean((columbus$HOVAL - lm_sp$fitted.values)^2))
## [1] 12.67496

Aplicamos el log a la variable dependiente y a Plumb, y el sq.square a Open air.

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
moran.test(exp(lm_sp_auto$residuals), rswm_queen)
## 
##  Moran I test under randomisation
## 
## data:  exp(lm_sp_auto$residuals)  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 0.39697, p-value = 0.3457
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.012160279      -0.020833333       0.006907894

Estadístico MoranI es bajo pero positivo; pero es necesario recurrir a modelos espaciales

Durbin model

  • retraso espacial que toma en cuenta las variables independientes y dependientes tomando en cuenta la relación de estás y su impacto posible en ellas mismas. Modelo Global.
spatial_autoregressive <- lagsarlm(log(HOVAL) ~ AREA + PERIMETER + INC + CRIME + OPEN^2 + log(PLUMB) + DISCBD, data = columbus, rswm_queen, type = 'mixed', Durbin = TRUE)
summary(spatial_autoregressive)
## 
## 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

Obtuvimos un p-value muy alto y no signf., pero podemos explicar un cambio hacia el valor de la vivienda de coeficientes negativos como Perimetro, Area, y Distancia al centro b. y crimen, con movimientos de la variable dep.

sqrt(mean((columbus$HOVAL - exp(spatial_autoregressive$fitted.values))^2))
## [1] 11.05289

RMSE-11.05

moran.test(exp(spatial_autoregressive$residuals), rswm_queen)
## 
##  Moran I test under randomisation
## 
## data:  exp(spatial_autoregressive$residuals)  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = -0.036866, p-value = 0.5147
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.024062346      -0.020833333       0.007671804

El Moran stat. muestra disperción y no acumulación de clusters.

spatial_autoregressive2 <- lagsarlm(HOVAL ~ INC + CRIME + PLUMB + DISCBD, data = columbus, listw = rswm_queen, Durbin = TRUE)
summary(spatial_autoregressive2)
## 
## Call:lagsarlm(formula = HOVAL ~ INC + CRIME + PLUMB + DISCBD, data = columbus, 
##     listw = rswm_queen, Durbin = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.4944  -7.0251  -2.1072   3.2795  46.5592 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept)  10.2288609  31.5085142  0.3246 0.745455
## INC           0.7063279   0.4539371  1.5560 0.119707
## CRIME        -0.4490543   0.1783994 -2.5171 0.011832
## PLUMB         1.5622410   0.5843119  2.6736 0.007503
## DISCBD       17.9078488   5.2253501  3.4271 0.000610
## lag.INC      -0.3807637   0.8541255 -0.4458 0.655746
## lag.CRIME     0.4607695   0.4342265  1.0611 0.288632
## lag.PLUMB     0.0082422   0.9928939  0.0083 0.993377
## lag.DISCBD  -14.0859390   7.0118549 -2.0089 0.044550
## 
## Rho: 0.19345, LR test value: 0.99066, p-value: 0.31958
## Asymptotic standard error: 0.18469
##     z-value: 1.0474, p-value: 0.29491
## Wald statistic: 1.0971, p-value: 0.29491
## 
## Log likelihood: -191.8448 for mixed model
## ML residual variance (sigma squared): 146.05, (sigma: 12.085)
## Number of observations: 49 
## Number of parameters estimated: 11 
## AIC: 405.69, (AIC for lm: 404.68)
## LM test for residual autocorrelation
## test value: 0.13504, p-value: 0.71327

Después de algunas iteraciones, y eliminar variables con p-values por encima de 0.05, obtuvimos el modelo con el menor valor de AIC entre los modelos, significancia, y mayor acierto en comparación a los otros modelos anteriores. Se cuenta con variables negativas explicativas como el crimen, el log-income y de discb, al igual que plumb (+).

Local Spatial Regression Analysis

Geographically Weighted Regression (GWR)

Determinar bandas Kernel

# 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

Obtenemos lista del número óptimo de vecinos parra mejorar el AIC.

GWR se realiza análisis a nivel local

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-05-12 23:57:02 
##    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-05-12 23:57:02

Existe un p-value menor a 0.05 signif., el AIC es alto comparado a los demás modelos.

sqrt(mean((columbus_shp$HOVAL - m.gwr$fitted.values)^2))
## [1] NaN

Muy bajo y negativo, con poca acumulación de vecinos

Aún se tendrá que repetir el modelo sin aquellas variables no signif;

m.gwr <- gwr.basic(HOVAL ~ CRIME + PLUMB, adaptive = T, data = columbus_shp, bw = bw1) 
m.gwr
##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2023-05-12 23:57:02 
##    Call:
##    gwr.basic(formula = HOVAL ~ CRIME + PLUMB, data = columbus_shp, 
##     bw = bw1, adaptive = T)
## 
##    Dependent (y) variable:  HOVAL
##    Independent variables:  CRIME PLUMB
##    Number of data points: 49
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.028  -8.698  -4.721   6.880  59.504 
## 
##    Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
##    (Intercept)  62.2678     4.9697  12.530  < 2e-16 ***
##    CRIME        -0.7681     0.1405  -5.467 1.81e-06 ***
##    PLUMB         1.3327     0.6043   2.205   0.0325 *  
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 14.68 on 46 degrees of freedom
##    Multiple R-squared: 0.3941
##    Adjusted R-squared: 0.3678 
##    F-statistic: 14.96 on 2 and 46 DF,  p-value: 9.892e-06 
##    ***Extra Diagnostic information
##    Residual sum of squares: 9917.329
##    Sigma(hat): 14.52608
##    AIC:  407.2567
##    AICc:  408.1658
##    BIC:  381.3913
##    ***********************************************************************
##    *          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 54.54673 58.51404 62.20620 64.35174 68.6618
##    CRIME     -0.88438 -0.79997 -0.76579 -0.70057 -0.6447
##    PLUMB      1.13947  1.24538  1.33620  1.42991  1.4890
##    ************************Diagnostic information*************************
##    Number of data points: 49 
##    Effective number of parameters (2trace(S) - trace(S'S)): 5.627275 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 43.37273 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 408.3652 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 399.8821 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 364.4335 
##    Residual sum of squares: 9128.596 
##    R-square value:  0.4422831 
##    Adjusted R-square value:  0.3682159 
## 
##    ***********************************************************************
##    Program stops at: 2023-05-12 23:57:02

Se cuenta con un AIC alto en comparativa y un r2 bajo.

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.69331, p-value = 0.2441
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -1.510417e-02     -2.083333e-02      6.828496e-05

Diagnóstico de Resultados Estimados - 4

Non-Spatial Regression Diagnostics

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

Es necesario remover las primeras 2 variables; podemos observar que ambos área como perímetro tienen valores mayores a 5 y presentan multicolinealidad, como se espera.

Lagrange Multiplier Diagnostic for Spatial Dependence (LMlag) - Dependencia y heterogeneidad espacial

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

Valor de p es mayor a 0.05 lo cual significa que el lag espacial de la variable dependiente no es requerida para el modelo.

Lagrange Multiplier Diagnostic for Spatial Error Dependence (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

Valor de p es mayor a 0.05 lo cual significa que el lag del error es requerida para el modelo.

Autocorrelación Espacial de los residuales estimados (εi)

columbus_shp$non_spatial_regression_residuals <-lm_sp$residuals
moran.test(columbus_shp$non_spatial_regression_residuals, rswm_queen) 
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$non_spatial_regression_residuals  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 0.97495, p-value = 0.1648
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.064083081      -0.020833333       0.007586055

El estadístico Moran I es bajo pero positivo, muestra que no hay concentrationa o clusters de los residuales de la regresion OLS.

Selección de Modelo - 5

Predicción local de variable dep. y mapeo

gwr_sf = st_as_sf(m.gwr$SDF)
summary(columbus$HOVAL)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   17.90   25.70   33.50   38.44   43.30   96.40

Vis. de predicción local de variables expl. estd. significantes

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=12, title="Rate range") +
   tm_layout(title= 'Precios de viviendas Columbus',  title.position = c('right', 'top'))
## Warning: Currect projection of shape gwr_sf unknown. Long-lat (WGS84) is
## assumed.

Podemos observar como nuestra variable dependiente de HOVAL (precios vivienda) muestra concentraciones de rangos altos en diversas áreas de la ciudad. Especialmente en los municipios de punta norte, y en el este, hay una concentración de los precios de las casas (altos rangos de la leyenda).

A continuación tenemos unos mapas extra de algunas variables explic. en donde se busca la realción tal cual con predición,.

tm_shape(gwr_sf) +
  tm_polygons(col = "CRIME", palette="PuBu", style="quantile", n=8, title="t-statistic") +
  tm_layout(title= 'Crimen',  title.position = c('right', 'top'))

Curiosamente crimen en el norte hay una mayor concentración de crimen, quizá atribuido al precio de casas (con mayor valor al hurtar quizá).

#tm_shape(gwr_sf) +
  #tm_polygons(col = "INC", palette="Greens", style="quantile", n=8, title="R2") +
 # tm_layout(title= 'Ingreso',  title.position = c('right', 'top'))

Es interesante ver además como habría concentración de ingreso en el sur, más no tanto valor de las casas.

Conclusión modelo y análisis - hallazgos

Finalmente, podemos concluir que el mejor modelo para la predicción del comportamiento de los precios de viviendas en Columbus, OH, USA, es el segundo model durbin - global. El modelo mantiene una buena cantidad de variables significativas y logicamente y estadisticamente explicativas ante cambios de la dependiente. Explica el AIC más bajo de los demás; se compara por más de dos: 403 a diferencia de 408 y 406. Con respecto al análisis: Podemos observar que el mapa predictivo muestra un incremento en los precios de viviendas alejadas del centro. También, aquellas concentradas en el Este y Norte son más costosas. Ambos tests de LaGrange muestran que el modelo lineal muestra dependencia espacial entre vecinos mediante distancias queen (no rook). El acceso al derenaje y el índice de crimen son significativas hacia los precios, en los modelos locales detallados. La zona con más dependencia de sus vecinos son las mismas de alta concentración, Este y Norte. Como se exploro en la gráfica de arriba, el norte tiene el mayor indice de crimen, mientras que el income en el sur. Un último hallazgo es que el distrito de la ciudad (financiero-business) no cuenta con viviendas de alto precio.



Referencias

Abdishakur. (2022, 8 abril). What is Exploratory Spatial Data Analysis (ESDA)? - Towards Data Science. Medium. https://towardsdatascience.com/what-is-exploratory-spatial-data-analysis-esda-335da79026ee How Spatial Autocorrelation (Global Moran’s I) works—ArcGIS Pro | Documentation. (s. f.). https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-statistics/h-how-spatial-autocorrelation-moran-s-i-spatial-st.htm Regression analysis basics—ArcGIS Pro | Documentation. (s. f.). https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-statistics/regression-analysis-basics.htm#:~:text=Regression%20analysis%20allows%20you%20to,factors%20behind%20observed%20spatial%20patterns What is Exploratory Data Analysis? | IBM. (s. f.). https://www.ibm.com/topics/exploratory-data-analysis What is Spatial Analysis? Definition and Examples. (s. f.). Spiceworks. https://www.spiceworks.com/tech/artificial-intelligence/articles/what-is-spatial-analysis/ What is Spatial Regression Analysis | IGI Global. (s. f.). https://www.igi-global.com/dictionary/spatial-regression-analysis/57225#:~:text=1.,variable%20and%20independent%20spatial%20variables