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.
Describir al menos 3-5 diferencias entre análisis exploratorio de datos (EDA) y análisis exploratorio espacial de datos (ESDA). Una de las principales y más notorias diferencias entre un análisis exploratorio de datos (EDA) y análisis exploratorio espacial de datos (ESDA), es en el enfoque ya que uno se aproxima más ala exploración de datos y tendencias, y el otro de ello con referencia e impacto a la ubicación geográfica y visual, asi como la distribución y casua espacial de variables. Una segunda diferencia es que el análisis exploratorio de datos no toma en cuenta la autocorrelación espacial, pero la ESDA si lo hace dándole importancia ala relación de las variables por ésto mismo. Una tercer diferencia es que el EDA busca entender y explorar histogramas, tendencias y dispersiones de datos a través del tiempo, mientras que el ESDA busca un impacto geospacial de vecinos, clusters, e impactos explicativos de una causa en un país, ciudad, municipio, etc. Similarmente, el EDA utiliza histogramas y gráficos de dispersión unicamente para explorar y sustentar resutlados de los datos, mientras que el ESDA heat maps, diagramas de dispersión espacial, matrices de correlacion en regiones y más.
Describir al menos 3-5 diferencias entre GWR y GRF. El modelo Geographically Weighted Regression realiza análisis de ubicaciones geospaciales, mientras que el Geostatistical Regression Framework utiliza una concentración especificada ya sea de regiones, municipios, etc. Se puede ver como global y local, respectivamente. El primero se ajusta la ubicación geospacial, sin embargo GRF se enfoca en una región especifica. El GWR aplica pesos (wieghted) hacia variables cercanas, mientras que el GRF una matriz de corr. espacial.
Describir cuáles son las 3-5 principales consecuencias de identificar autocorrelación espacial en los residuales de un modelo de regresión estimado.
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
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)
# 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.
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)
# 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
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.
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
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
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.
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
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.
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.
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.
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.
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
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.
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