# Se cargan las librerías
library(sf)

## Warning: package 'sf' was built under R version 4.2.3

## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(sp)

## Warning: package 'sp' was built under R version 4.2.3

library(tmap)

## Warning: package 'tmap' was built under R version 4.2.3

## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')

library(gstat)

## Warning: package 'gstat' was built under R version 4.2.3

library(spdep)

## Warning: package 'spdep' was built under R version 4.2.3

## 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')`

library(leaflet)

## Warning: package 'leaflet' was built under R version 4.2.3

# se carga la data 
data(columbus)

La base de datos Columbus en la librería spdep de R es un conjunto de datos ampliamente utilizado para la demostración y análisis de modelos espaciales. La base contiene información socioeconómica y de criminalidad para 49 vecindarios de la ciudad de Columbus, Ohio, EE. UU.

Descripción de la base de datos:

Número de observaciones: 49 (correspondientes a los vecindarios). Número de variables: 7.

Variables principales:

CRIME:

Tasa de crímenes per cápita en cada vecindario.

INC: Ingreso promedio per cápita (en miles de dólares).
HOVAL: Valor promedio de las viviendas en miles de dólares.
POPN: Población en cada vecindario.
PLUMB: Porcentaje de viviendas con plomería completa.
DISCBD: Distancia al centro de negocios principal (en millas).
OWNH: Porcentaje de viviendas ocupadas por propietarios

Análisis Descriptivo y exploratorio de los datos para la variable INC (Ingreso Promedio per cápita (en miles de dólares))

# summary de los datos
summary(columbus$INC)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.477   9.963  13.380  14.375  18.324  31.070

A partir del summary se puede observar que el valor mínimo para la variable INC es 4.47 en la ciudad de Columbus, ohio, el promedio del ingreso promedio percápita es de 14.375 y el valor máximo de 31.070. Se puede apreciar tambien que hasta el 25 % que corresponde el primer cuantil, se acumulan 9.963 de ingreso promedio percápita y el segundo cuantil 50 % es cercano al valor de la media 13.380 como valor percápita.

# Boxplot de la variable INC
boxplot(columbus$INC,col = "blue", 
        border = 4, xlab="Ingreso promedio Percápita (INC)")

A partir del boxplot se puede apreciar para la variable INC que los datos están concentrados, aunque se aprecia un dato atípico, la caja del boxplot está mas hacia abajo, es decir el bigote superior del boxplot es mas largo y podrías decir que la mayoría de las observaciones están por encima de la mediana.

Calculo de la coordenada media

data("columbus")
coordinates(columbus) = ~X+Y
coordenada_media1 <- colMeans(coordinates(columbus))
coordenada_media1

##        X        Y 
## 39.46429 32.37265

Interpretación:

La coordenada X (39.46429) y la coordenada Y (32.37265) juntas representan el centroide de los 49 vecindarios de Columbus en el espacio geográfico. Este es el punto promedio donde, en teoría, se ubicaría el “centro” de todos los vecindarios en tu análisis.

Contexto Geográfico:

Este centroide puede ayuda a visualizar cómo están distribuidos los vecindarios respecto a este punto central. Si la mayoría de los vecindarios están cerca de esta coordenada, la distribución es más centralizada; si están dispersos alrededor de ella, la distribución es más esparcida.

Histograma Regional

# histograma
# Se convierte la ada en Data Frame
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

columbus1 <- as.data.frame(columbus)

ggplot(columbus1, aes(x = INC)) +
  geom_histogram(binwidth = 10, fill = "#F8766D", color = "red") +
  labs(title = "Histograma del Ingreso promedio per cápita (en miles de dólares)",
       x = "INC geográficamente",
       y = "Frecuencia") +
  theme_minimal()

A partir del histograma se puede apreciar que la mayoría de las observaciones de INC en la ciudad de Columbus se encuentran dentro del rango de 10 a 20, tambien se puede decir que para los valores mas bajos y altos hay menor frecuencia lo que se podría decir que entre 8 a 25 INC están la mayoría de los datos.

Patrones y tendencias.

Diagrama de Dispersión Básico (INC vs. HOVAL) Muestra la relación entre el ingreso promedio per cápita (INC) y el valor promedio de las viviendas (HOVAL).

# Gráfico de tendencias 
library(ggplot2)

ggplot(columbus1, aes(x = INC, y = HOVAL, col="red")) +
  geom_point() +
  labs(title = "Diagrama de Dispersión: INC vs. HOVAL",
       x = "Ingreso promedio per cápita (INC)",
       y = "Valor promedio de las viviendas (HOVAL)")

cor(columbus1$INC,columbus1$HOVAL)

## [1] 0.4998786

Interpretación:

Se puede observar una relación positiva entre éstas dos variables y que a medida que aumenta el ingreso percápita tambien aumenta el valor promedio de las viviendas en la ciudad de Columbus, aunque el valor de la correlación fué de 0.50 lo que quiere decir que no es muy buena pero si existe una correlación entre las variables. Se puede afirmar que los datos atípicos que se observan pueden ser lo que estén ocasionando que el valor de la correlación disminuya.

Realice el mapa.

library(sf)
library(st)

## Warning: package 'st' was built under R version 4.2.3

## Loading required package: sda

## Warning: package 'sda' was built under R version 4.2.3

## Loading required package: entropy

## Loading required package: corpcor

## Loading required package: fdrtool

library(corpcor)
library(fdrtool)
library(sda)
library(entropy)
library(spdep)
columbus <- st_read(system.file("etc/shapes/columbus.shp", package="spdep"))

## Reading layer `columbus' from data source 
##   `C:\Users\User\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

col.gal.nb  <- read.gal(system.file("etc/weights/columbus.gal", package="spdep"))
columbus_sf <- read_sf(system.file("etc/shapes/columbus.shp", package="spdep"))
tm_shape(columbus) + tm_polygons(col='wheat') + 
  tm_style("classic") + 
  tm_text(text='POLYID',size=0.7)

## Warning: Currect projection of shape columbus unknown. Long-lat (WGS84) is
## assumed.

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(columbus) + 
  tm_fill("INC", style="quantile", title = "House Prices (Quantile)") + 
  tm_layout(main.title = "Columbus, Ohio", legend.position = c("left", "top"), 
            legend.title.size = 0.8, legend.text.size = 0.7)

## Warning: Currect projection of shape columbus unknown. Long-lat (WGS84) is
## assumed.

## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.

Mapa de cuartiles

Un mapa de cuantiles o un “box map” es una representación espacial de la variable que destaca las observaciones por cuartiles.

library(tmap)
# Crear una nueva variable con los cuartiles de INC
columbus$INC_quantile <- cut(columbus$INC, 
                             breaks = quantile(columbus$INC, probs = seq(0, 1, 0.25)), 
                             include.lowest = TRUE, 
                             labels = c("1er Cuartil", "2do Cuartil", "3er Cuartil", "4to Cuartil"))

# Mapa de cuartiles de INC
tm_shape(columbus) +
  tm_dots(col = "INC_quantile", palette = "-RdBu", title = "Cuantiles de Ingreso") +
  tm_layout(title = "Box Map: Ingreso en Columbus")

## Warning: Currect projection of shape columbus unknown. Long-lat (WGS84) is
## assumed.

Gráfico de Varianza

library(sf)
library(gstat)
library(sp)
library(spdep)
# Asumimos que columbus_sf es un objeto sf previamente cargado y contiene la variable INC# Convertir el objeto sf a Spatial (si no lo has hecho ya)
columbus_sp <- as(columbus_sf,"Spatial")# Crear una cuadrícula sobre la cual se hace la interpolación
columbus_grid <- st_make_grid(columbus_sf, n =c(100,100), what ="centers")
columbus_grid_sf <- st_sf(geometry = columbus_grid, crs = st_crs(columbus_sf))

# Ajuste de un modelo teórico al variograma nuevo
# psill: Es el valor de la meseta
# range: Es el valor que se le asigna al rango
# nugget: Es el valor de la pepita

# tienes un modelo de variograma
variogram_model <- variogram(log(INC)~1, columbus_sp)

model_fit <- fit.variogram(variogram_model, model = vgm(psill =0.11, model ="Sph",range=0.72, nugget =0.2));model_fit

# Gráfico del modelo ajustado al modelo esférico
library(ggplot2)
ggplot(variogramLine(model_fit, 1.7), aes(x = dist, y = gamma)) +
  geom_path() + 
  geom_point(data = variogram_model, color = "blue") + 
  geom_vline(xintercept = model_fit$range[2], lty = 2) + 
  geom_text(x = model_fit$range[2], y = model_fit$psill[2]/2, label = "range") + 
  theme_bw() +
  geom_hline(yintercept = model_fit$psill[2] + model_fit$psill[1] , lty = 2) + 
  geom_text(x = model_fit$range[2]/2, y = model_fit$psill[2] + model_fit$psill[1],
            label = "psill") +
  geom_text(x = model_fit$range[1], y = model_fit$psill[1], label = "Nugget") +
  ylim(c(0,max(variogram_model$gamma)))

Interpretación:

El valor aproximado de la peipa es 0.016, de la meseta = 0.13 y del rango respectivamente es 1.3

Ahora se realizará el Krigging

#cuadrícula a SpatialPointsDataFrame
columbus_grid_sp <- as(columbus_grid_sf,"Spatial")

kriging_result <- krige(log(INC)~1, columbus_sp, columbus_grid_sp, model = model_fit)

## [using ordinary kriging]

# Visualización de los resultados de Kriging
spplot(kriging_result["var1.pred"], main ="Kriging ordinario (INC)")

Interpretación:

Se puede ver que los datos espaciales para INC se encuentran bien concentrados con respecto a sus niveles aunque es de aclarar que el gráfico para el krigging generado no tiene la mejor visualización.

Realizar el gráfico de Moran

El gráfico de Moran visualiza la autocorrelación espacial para la variable INC.

# Crear la lista de vecinos y la matriz de pesos espaciales
columbus_sp <- as(columbus,"Spatial")
# Crear una matriz de vecinos basados en contigüidad
nb <- poly2nb(columbus_sp)# Convertir los vecinos a una lista de pesos
listw <- nb2listw(nb, style ="W")
#vecinos <- knn2nb(knearneigh(coordinates(columbus), k=4))
#pesos <- nb2listw(vecinos, style = "W")

# Calcular el I de Moran
morann <- moran.test(columbus_sp$INC, listw)
print(morann)

## 
##  Moran I test under randomisation
## 
## data:  columbus_sp$INC  
## weights: listw    
## 
## 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

Interpretación:

El resultado del test de Moran I indica que hay una autocorrelación espacial positiva significativa en la variable INC en los datos de Columbus. Vamos a desglosar lo que significa cada parte del resultado:

Desglose del Resultado

Moran’s I statistic: 0.4156

Este valor de 0.4156 sugiere una autocorrelación espacial positiva moderada. Un valor cercano a 1 indicaría una fuerte autocorrelación positiva (las áreas cercanas tienen valores similares), mientras que un valor cercano a -1 indicaría autocorrelación negativa (las áreas cercanas tienen valores muy diferentes). Un valor cercano a 0 indicaría que no hay autocorrelación espacial.

Expectation: -0.0208

Este es el valor esperado de Moran’s I bajo la hipótesis nula de que no hay autocorrelación espacial. En este caso, el valor esperado es cercano a 0, lo cual es típico en estos tests.

Variance: 0.0084

La varianza bajo la hipótesis nula, utilizada para calcular el valor p.

Moran I statistic standard deviate: 4.7645

Este valor es una medida de cuántas desviaciones estándar está la estadística de Moran I observada del valor esperado bajo la hipótesis nula. Un valor alto en magnitud indica una evidencia fuerte contra la hipótesis nula.

p-value: 9.467e-07

El valor p es extremadamente bajo, indicando que la autocorrelación espacial positiva observada es altamente significativa. Esto significa que es muy improbable que la autocorrelación observada ocurra por azar.

Alternative hypothesis: greater

La prueba está evaluando si el valor observado de Moran I es mayor que el esperado bajo la hipótesis nula, lo que indica autocorrelación espacial positiva.

Interpretación

En resumen, el test de Moran I sugiere que hay una autocorrelación espacial positiva significativa en la variable INC. Esto significa que las áreas cercanas tienden a tener valores de INC similares, lo cual puede ser relevante para la comprensión de la variable estudiada.

# gráfico
moran.plot(columbus_sp$INC, listw, pch = 16, col = "#2b83ba", 
           main = "Gráfico de Moran: Ingreso promedio percápita", xlab = "INC")

Explicación del Gráfico de Moran

Interpretación del Gráfico

Cuadrantes:

Primer Cuadrante (Arriba a la derecha):

Alto-alto. Zonas con altos valores de INC rodeadas por zonas con altos valores.

Segundo Cuadrante (Arriba a la izquierda): Bajo-alto. Zonas con bajos valores de INC rodeadas por zonas con altos valores.

Tercer Cuadrante (Abajo a la izquierda): Bajo-bajo. Zonas con bajos valores de INC rodeadas por zonas con bajos valores.

Cuarto Cuadrante (Abajo a la derecha): Alto-bajo. Zonas con altos valores de INC rodeadas por zonas con bajos valores.

Interpretación Global:

Si la mayoría de los puntos están en los cuadrantes I y III, esto indica una fuerte autocorrelación espacial positiva (los lugares con valores similares están cerca unos de otros). Si los puntos se distribuyen por todos los cuadrantes, la autocorrelación espacial es más débil o no está presente.

Ajuste de los modelos esférico, exponencial, Gaussiano, monómico, y Pepita Pura

Modelo esférico

# Modelo esférico
# tienes un modelo de variograma
variogram_model <- variogram(log(INC)~1, columbus_sp)

model_fit <- fit.variogram(variogram_model, model = vgm(psill =0.11, model ="Sph",range=0.72, nugget =0.2));model_fit

# Gráfico del modelo ajustado al modelo esférico
library(ggplot2)
ggplot(variogramLine(model_fit, 1.7), aes(x = dist, y = gamma)) +
  geom_path() + 
  geom_point(data = variogram_model, color = "blue") + 
  geom_vline(xintercept = model_fit$range[2], lty = 2) + 
  geom_text(x = model_fit$range[2], y = model_fit$psill[2]/2, label = "range") + 
  theme_bw() +
  geom_hline(yintercept = model_fit$psill[2] + model_fit$psill[1] , lty = 2) + 
  geom_text(x = model_fit$range[2]/2, y = model_fit$psill[2] + model_fit$psill[1],
            label = "psill") +
  geom_text(x = model_fit$range[1], y = model_fit$psill[1], label = "Nugget") +
  ylim(c(0,max(variogram_model$gamma)))

Krigging para el modelo esférico

#cuadrícula a SpatialPointsDataFrame
columbus_grid_sp <- as(columbus_grid_sf,"Spatial")

kriging_result <- krige(log(INC)~1, columbus_sp, columbus_grid_sp, model = model_fit)

## [using ordinary kriging]

# Visualización de los resultados de Kriging
spplot(kriging_result["var1.pred"], main ="Kriging ordinario (INC)")

Modelo Exponencial

# Modelo exponencial
# tienes un modelo de variograma
variogram_model <- variogram(log(INC)~1, columbus_sp)

model_fit1 <- fit.variogram(variogram_model, model = vgm(psill =0.11, model ="Exp",range=0.72, nugget =0.2));model_fit1

# Gráfico del modelo ajustado al modelo exponencial
library(ggplot2)
ggplot(variogramLine(model_fit1, 1.7), aes(x = dist, y = gamma)) +
  geom_path() + 
  geom_point(data = variogram_model, color = "blue") + 
  geom_vline(xintercept = model_fit1$range[2], lty = 2) + 
  geom_text(x = model_fit1$range[2], y = model_fit1$psill[2]/2, label = "range") + 
  theme_bw() +
  geom_hline(yintercept = model_fit1$psill[2] + model_fit1$psill[1] , lty = 2) + 
  geom_text(x = model_fit1$range[2]/2, y = model_fit1$psill[2] + model_fit1$psill[1],
            label = "psill") +
  geom_text(x = model_fit1$range[1], y = model_fit1$psill[1], label = "Nugget") +
  ylim(c(0,max(variogram_model$gamma)))

Krigging para el modelo exponencial

#cuadrícula a SpatialPointsDataFrame
#columbus_grid_sp <- as(columbus_grid_sf,"Spatial")

kriging_result1 <- krige(log(INC)~1, columbus_sp, columbus_grid_sp, model = model_fit1)

## [using ordinary kriging]

# Visualización de los resultados de Kriging
spplot(kriging_result1["var1.pred"], main ="Kriging ordinario (INC)")

Modelo Gaussiano

# Modelo Gaussiano
# tienes un modelo de variograma
variogram_model <- variogram(log(INC)~1, columbus_sp)

model_fit2 <- fit.variogram(variogram_model, model = vgm(psill =0.11, model ="Gau",range= 1, nugget =0.2));model_fit2

# Gráfico del modelo ajustado al modelo Gaussiano
library(ggplot2)
ggplot(variogramLine(model_fit2, 1.7), aes(x = dist, y = gamma)) +
  geom_path() + 
  geom_point(data = variogram_model, color = "blue") + 
  geom_vline(xintercept = model_fit2$range[2], lty = 2) + 
  geom_text(x = model_fit2$range[2], y = model_fit2$psill[2]/2, label = "range") + 
  theme_bw() +
  geom_hline(yintercept = model_fit2$psill[2] + model_fit2$psill[1] , lty = 2) + 
  geom_text(x = model_fit2$range[2]/2, y = model_fit2$psill[2] + model_fit2$psill[1],
            label = "psill") +
  geom_text(x = model_fit2$range[1], y = model_fit2$psill[1], label = "Nugget") +
  ylim(c(0,max(variogram_model$gamma)))

Krigging para el modelo Gaussiano

#cuadrícula a SpatialPointsDataFrame
#columbus_grid_sp <- as(columbus_grid_sf,"Spatial")

kriging_result2 <- krige(log(INC)~1, columbus_sp, columbus_grid_sp, model = model_fit2)

## [using ordinary kriging]

# Visualización de los resultados de Kriging
spplot(kriging_result2["var1.pred"], main ="Kriging ordinario (INC)")

Modelo Monómico

# Modelo monómico
# tienes un modelo de variograma
variogram_model <- variogram(log(INC)~1, columbus_sp)

model_fit3 <- fit.variogram(variogram_model, model = vgm(psill =0.11, model ="Pow",range= 1, nugget =0.2));model_fit3

# Gráfico del modelo ajustado al modelo Monómico
library(ggplot2)
ggplot(variogramLine(model_fit3, 1.7), aes(x = dist, y = gamma)) +
  geom_path() + 
  geom_point(data = variogram_model, color = "blue") + 
  geom_vline(xintercept = model_fit3$range[2], lty = 2) + 
  geom_text(x = model_fit3$range[2], y = model_fit3$psill[2]/2, label = "range") + 
  theme_bw() +
  geom_hline(yintercept = model_fit3$psill[2] + model_fit3$psill[1] , lty = 2) + 
  geom_text(x = model_fit3$range[2]/2, y = model_fit3$psill[2] + model_fit3$psill[1],
            label = "psill") +
  geom_text(x = model_fit3$range[1], y = model_fit3$psill[1], label = "Nugget") +
  ylim(c(0,max(variogram_model$gamma)))

Krigging para el modelo Monómico

#cuadrícula a SpatialPointsDataFrame
columbus_grid_sp <- as(columbus_grid_sf,"Spatial")

kriging_result3 <- krige(log(INC)~1, columbus_sp, columbus_grid_sp, model = model_fit3)

## [using ordinary kriging]

# Visualización de los resultados de Kriging
spplot(kriging_result3["var1.pred"], main ="Kriging ordinario (INC)")

Modelo Pepita Pura

# Modelo Pepita Pura
# tienes un modelo de variograma
variogram_model <- variogram(log(INC)~1, columbus_sp)

model_fit4 <- fit.variogram(variogram_model, model = vgm( model ="Nug"));model_fit4

# Gráfico del modelo ajustado al modelo Pepita pura
library(ggplot2)
ggplot(variogramLine(model_fit4, 1.7), aes(x = dist, y = gamma)) +
  geom_path() + 
  geom_point(data = variogram_model, color = "blue") + 
  geom_vline(xintercept = model_fit4$range[2], lty = 2) + 
  geom_text(x = model_fit4$range[2], y = model_fit4$psill[2]/2, label = "range") + 
  theme_bw() +
  geom_hline(yintercept = model_fit4$psill[2] + model_fit4$psill[1] , lty = 2) + 
  geom_text(x = model_fit4$range[2]/2, y = model_fit4$psill[2] + model_fit4$psill[1],
            label = "psill") +
  geom_text(x = model_fit4$range[1], y = model_fit4$psill[1], label = "Nugget") +
  ylim(c(0,max(variogram_model$gamma)))

## Warning: Removed 1 rows containing missing values (`geom_vline()`).

## Warning: Removed 200 rows containing missing values (`geom_text()`).

## Warning: Removed 1 rows containing missing values (`geom_hline()`).

## Warning: Removed 200 rows containing missing values (`geom_text()`).

Conclusión:

El mejor modelo, el que se ajusta mucho mejor a los datos es el modelo Gaussiano, se puede apreciar que genera mejores estimaciones en comparación a los otros modelos.

Estadistica Espacial: Análisis de la data Columbus

Javier Daza

2024-08-16

Descripción de la base de datos:

Variables principales:

Análisis Descriptivo y exploratorio de los datos para la variable INC (Ingreso Promedio per cápita (en miles de dólares))

Calculo de la coordenada media

Interpretación:

Histograma Regional

Patrones y tendencias.

Interpretación:

Realice el mapa.

Mapa de cuartiles

Gráfico de Varianza

Interpretación:

Ahora se realizará el Krigging

Interpretación:

Realizar el gráfico de Moran

Interpretación:

Desglose del Resultado

Interpretación

Explicación del Gráfico de Moran

Interpretación del Gráfico

Interpretación Global:

Ajuste de los modelos esférico, exponencial, Gaussiano, monómico, y Pepita Pura

Modelo esférico

Krigging para el modelo esférico

Modelo Exponencial

Krigging para el modelo exponencial

Modelo Gaussiano

Krigging para el modelo Gaussiano

Modelo Monómico

Krigging para el modelo Monómico

Modelo Pepita Pura

Conclusión: