Bitácora Final Estadística Espacial

Datos de Área

A continuación, se realiza un modelamiento espacial a partir de datos de área, concerniente al modulo impartido, para esto se construyó una base de datos para Cundinamarca, que representa la cantidad de hurtos a personas para el año 2019 que se pueden encontrar en el link Estadística delictiva, además, se integro la información del Censo Nacional de Población y Vivienda - CNPV, teniendo en cuenta los valores de proporción por miseria y proporción de necesidades básicas insatisfechas.

library(rgdal)
CUNDI<-readOGR("~/Estadistica Espacial/Data/ROBOS/robos_cund.shp", layer="robos_cund")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/20151025100/Estadistica Espacial/Data/ROBOS/robos_cund.shp", layer: "robos_cund"
## with 117 features
## It has 10 fields

head(CUNDI@data[2:10])

##   COD_DANE TOTALPOB         MPIO TOTAL CASOS_MPIO RURAL URBANA  PROP_NBI
## 0    25001    10742 AGUA DE DIOS    53         53    10     43  8.079839
## 1    25599     7533        APULO    32         32    11     21 18.406375
## 2    25035    12241     ANAPOIMA    74         73    21     52 10.338005
## 3    25386    29452      LA MESA   267        266    28    237  8.911677
## 4    25307    92903     GIRARDOT  1697       1688    40   1648  7.536184
## 5    25815    13649      TOCAIMA   111        108    30     78 11.171439
##   PROP_MISER
## 0  1.1227329
## 1  3.0810093
## 2  1.3602638
## 3  1.1946689
## 4  0.8859031
## 5  1.4280402

library(leaflet)
pal <- colorNumeric(
  palette = "YlOrRd",
  domain = (CUNDI@data$TOTAL)
)

state_popup <- paste0("<strong>Municipio: </strong>", 
                      CUNDI@data$MPIO, 
                      "<br><strong>Conteo de casos por Hurto a nivel municipal, año 2019: </strong>", 
                      CUNDI@data$TOTAL)

leaflet(CUNDI)%>% 
  #MAPA BASE
  addTiles()%>%addProviderTiles(providers$Esri.NatGeoWorldMap)%>%
  #CAPAS
  addPolygons(color = "#444444", weight = 1, smoothFactor = 0.5,
              popup = state_popup,
              opacity = 1.0, fillOpacity = 0.5,
              fillColor =~colorQuantile("YlOrRd", TOTAL)(TOTAL),
              highlightOptions = highlightOptions(color = "white", weight = 2,bringToFront = TRUE),group = "Hurtos")%>%
  #LEYENDA
  addLegend("bottomleft", pal = pal ,values = ~TOTAL,title = "Total de robos en Cundinamarca año 2019",opacity = 1)%>%
  #CONTROL CAPAS
  addLayersControl(
    #baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
    overlayGroups = c("Hurtos"),
    options = layersControlOptions(collapsed = FALSE)
  ) 

Observando el mapa anterior, se tienen dos datos con valores nulos para la variable respuesta, por lo tanto procedemos a la limpieza de estos para que no afecten la inferencia del modelo.

CUNDI <-CUNDI[is.na(CUNDI@data$MPIO) == F,] #ELIMINANDO VALORES NA

Analísis Exploratorio de datos espaciales (AEDE)

Se procede con el análisis univariante para la variable, mediante estadísticas descriptivas y métodos gráficos como el histograma y diagrama de cajas y bigotes para observar el comportamiento de la variable objeto de estudio.

Para la variable respuesta de total de hurtos en Cundinamarca, se puede observar a partir de sus estadísticas que no cuenta con una distribución normal, dado quue su media es superior a la mediana (quantil 2 - percentil 50), ademásse tienen valores muy altos que pueden asociarse a Bogotá D.C. dado que este cuenta con mayor cantidad de robos respecto a los municipios de Cundinamarca.

library(ggplot2)
z<-as.numeric(CUNDI@data$TOTAL)
summary(z)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##      1.0     10.0     25.0   2262.2    104.5 233277.0

#HISTOGRAMA PARA VER EL COMPORTAMIENTO DE LA VARIABLE DEPENDIENTE
df <- data.frame(PF = z)
ggplot(df, aes(x = PF)) + 
  geom_histogram(aes(y =..density..), 
                 colour = "black", 
                 fill = "white") +
  stat_function(fun = dnorm, args = list(mean = mean(df$PF), sd = sd(df$PF)))

De acuerdo con el histograma de la variable respuesta, se observa que el histograma, es asimétrico positivo, corroborando lo que las estadísticas descriptivas inicialmente descritas.

Además, se observa el peso que tiene la observación atípica dado que aplana la curva de densidad de probabilidad.

boxplot(z, horizontal = T)

A partir del gráfico de Cajas y bigotes, se observan varios valores atípicos a una distancia de 1.5 del rango intercuartilico, estos valores pueden llegar a influenciar la muestra, sin embargo, podemos esclarecer el peso de estas observaciones al realizar el análisis de normalidad.

qqnorm(z)
qqline(z,col="blue")

Si observamos el gráfico de normalidad, se tienen dos observaciones fuera de la línea de ajuste de normalidad, pese a esto, los demás valores se encuentran dentro del ajuste, por lo que es probable que la muestra siga una distribución normal. Para corroborar está información se procede a realizar un test estadístico con una significancia del 5%.

shapiro.test(z)

## 
##  Shapiro-Wilk normality test
## 
## data:  z
## W = 0.077148, p-value < 2.2e-16

Como el (P_value) es menor al 5%, se rechaza la hipótesis nula, por lo tanto, corroboramos que la variable de estudio NO se comporta bajo una distribución normal. Por tal motivo se procede a transformar la variable objeto de estudo Hurtos para que cumpla con el supuesto de normalidad.

library(MASS)
l<-boxcox(lm(z ~ 1))

Se realiza una transformación box-cox de la familia exponencial, para tratar de transformar la variable objeto de estudio como se muestra en la ecuación a continuación. \[z\_trans=\frac{z^{\lambda-1}}{\lambda}\]

Para esto se optimiza el valor de lambda. Y se agrega la variable transformadad al dataset inicial (Shapefile)

lambda <- l$x[which.max(l$y)]

new_z_exact <- (z ^ lambda - 1) / lambda

df <- data.frame(PF = new_z_exact)
ggplot(df, aes(x = PF)) + 
  geom_histogram(aes(y =..density..), 
                 colour = "black", 
                 fill = "white") +
  stat_function(fun = dnorm, args = list(mean = mean(df$PF), sd = sd(df$PF)))

shapiro.test((new_z_exact))

## 
##  Shapiro-Wilk normality test
## 
## data:  (new_z_exact)
## W = 0.99109, p-value = 0.6647

Posterior a la transformación de la variable, se realiza el gráfico asociado al histograma de la variable transformada, y se muestra su función de distribución de probabilidad, en donde se aprecia que hay un mejor ajuste de esta, en comparativa con el histograma inicial. Se corre posteriormente un test de shapiro para esclarecer si la variable sigue una distribución normal.

Como el (P_value) es mayor al 5%, se acepta la hipótesis nula, por lo tanto, corroboramos que la variable de estudio transformada se comporta bajo una distribución normal.

Análisis de la forma funcional

Se realizan diagramas de dispersión en donde se tiene la variable respuesta hurtos (la variable transformada) y se confronta con las demás variables a considerar dentro del modelo.

CUNDI@data$TOTAL_TRANS<-new_z_exact
#RESPECTO A LA POBLACIÓN, SIN TENER EN CUENTA LOS MAYORES A 600.000 HA
ggplot(CUNDI[CUNDI@data$TOTALPOB<600000,]@data,
       aes(log(TOTALPOB),TOTAL_TRANS)) + geom_point()+ geom_smooth(method = "lm", se = FALSE)

Inicialmente se realiza el diagrama de dispersión de la variable transformada con el total de la población, sin embargo, se observa que este tiene un comportamiento logaritmico, por lo tanto se le aplica esta función para linealizar y esclarecer la dependencia funcional de la variable independiente sobre la dependiente. En este sentido, se observa que hay una relación positiva entre la variable transformada y el logaritmo de la poblaión.

#RESPECTO A LA PROPORCIÓN DE NBI
ggplot(CUNDI@data,aes(log(PROP_NBI),TOTAL_TRANS))+ geom_point()+ geom_smooth(method = "lm", se = FALSE)

Para la proporción de NBI se observa que se debe linealizar tal cual como el total de la población, sin embargo, la relación de estas es inversa, esto quiere decir que un aumento en la proporción de NBI disminuirá la variable transformada, esto no quiere decir que los casos disminuyan, simplemente su transformación para cumplir con el supuesto de normalidad.

#RESPECTO A LA PROPORCIÓN DE MISERIA
ggplot(CUNDI@data,aes(log(PROP_MISER),TOTAL_TRANS))+ geom_point()+ geom_smooth(method = "lm", se = FALSE)

Para la variable independiente de proporción de miseria se tiene el mismo comportamiento que la variable anteriormente descrita. ## Modelamiento

De acuerdo con el capitulo anterior se tiene planteada la siguiente ecuación que corresponde con la forma funcional para el modelo. \[TOTAL\_TRANS=\beta_0+\beta_1log(TOTALPOB)+\beta_2log(PROP\_NBI)+\beta_3log(PROP\_MISER)+\epsilon\]

modelo_clasico <- lm(CUNDI@data$TOTAL_TRANS~log(CUNDI@data$TOTALPOB)
                     +log(CUNDI@data$PROP_NBI)+log(CUNDI@data$PROP_MISER))
summary(modelo_clasico)

## 
## Call:
## lm(formula = CUNDI@data$TOTAL_TRANS ~ log(CUNDI@data$TOTALPOB) + 
##     log(CUNDI@data$PROP_NBI) + log(CUNDI@data$PROP_MISER))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.65196 -0.17297  0.08965  0.32422  0.81289 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                -0.72260    0.92466  -0.781  0.43619    
## log(CUNDI@data$TOTALPOB)    0.55205    0.04765  11.586  < 2e-16 ***
## log(CUNDI@data$PROP_NBI)   -0.88675    0.27174  -3.263  0.00146 ** 
## log(CUNDI@data$PROP_MISER)  0.30216    0.13390   2.257  0.02599 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4569 on 111 degrees of freedom
## Multiple R-squared:  0.779,  Adjusted R-squared:  0.7731 
## F-statistic: 130.4 on 3 and 111 DF,  p-value: < 2.2e-16

Para la primera estimación se tiene un modelo clásico, dentro de este observamos que las covariables o variables independientes son significativas al 5%, además, la proporción de NBI es significativa al 1% y el total de la población al 0.1%, además, el modelo es significativo dado que el P_VALUE es menor al 1%. El ajuste del modelo es del 77.31%, quiere decir que las variables independientes explican en un 77.31% a la variable transformada del total de hurtos.

Vecindades

Antes de que se pueda modelar formalmente la dependencia que se muestra, primero se debe cubrir cómo los vecindarios están conectados espacialmente entre sí. Para esto se debe definir la conectividad y los pesos de los vecinos.

Se establece un lag de 4 vecinos dada la región, en donde, se establece que una región puede afectar hasta su cuarto orden o vecino más cercano.

#set.ZeroPolicyOption(T)
library(spdep)
towwer <- poly2nb(CUNDI, queen=F)            # Torre
cspc <- sp.correlogram(towwer, CUNDI@data$TOTAL_TRANS, order=4, method="I",
                       zero.policy=TRUE)
plot(cspc)

Al observar el correlograma a partir de la vecindad de tipo torre, se tiene que a mayor orden disminuye el indice de moran, eliminando la autocorrelación o asociación e incidencia que puede tener la variable transformada.

queen <- poly2nb(CUNDI, queen=T)              #reina
cspc <- sp.correlogram(queen, CUNDI@data$TOTAL_TRANS, order=4, method="I",
                       zero.policy=TRUE)
plot(cspc)

Al observar el correlograma a partir de la vecindad de tipo reina, se tiene que a mayor orden disminuye el indice de moran, eliminando la autocorrelación o asociación e incidencia que puede tener la variable transformada.

coordenadas<-coordinates(CUNDI)    #k-vecinos
k1 <- knn2nb(knearneigh(coordenadas,4))
cspc <- sp.correlogram(k1, CUNDI@data$TOTAL_TRANS, order=4, method="I",
                       zero.policy=TRUE)
plot(cspc)

Otro método común para definir vecinos es k-vecinos más cercanos. Este método encuentra las k observaciones más cercanas para cada observación de interés, donde k es algún entero. En este caso, para hacer el simil con los métodos anteriormente mencionados, se deja un k=4.

A partir de estos correlogramas se puede observar que sin importar el método, se tiene el mismo comportamiento por lo cual a mayor aumento en el orden, se disminuye la correlación, por lo tanto la contigüidad de orden 1 parece ser la más adecuada sin importar el método.

I de Moran

Los estadísticos globales de autocorrelación constituyen la aproximación más tradicional al efecto de dependencia espacial, permitiendo contrastar la presencia o ausencia de un esquema de dependencia espacial a nivel univariante, es decir, contrastar si se cumple la hipótesis de que una variable se encuentra distribuida de forma totalmente aleatoria en el espacio o si por el contrario existe una asociación significativa de valores similares o disímiles entre regiones vecinas.

A partir de la matriz de vecindad de orden 1 con efecto reina, se tiene la matriz de pesos para estimar el I de Moran.

queen <- poly2nb(CUNDI, queen=T)              #reina
queen<-nb2listw(queen, style="W", glist = NULL)

Se plotea el diagrama de dispersión del índice de Moran, en donde se puede apreciar que existe una asociación positiva o autocorrelación positiva, además, existen varios municipios que pueden estar correlacionados de forma negativa.

moran.plot(CUNDI@data$TOTAL_TRANS, queen, zero.policy=NULL)

Para determinar si existe correlación de forma global, se da apartir de un test estadístico, allí se encuentra el valor detallado que observabamos en el correlograma.

moran.test(CUNDI@data$TOTAL_TRANS, queen)

## 
##  Moran I test under randomisation
## 
## data:  CUNDI@data$TOTAL_TRANS  
## weights: queen    
## 
## Moran I statistic standard deviate = 7.6183, p-value = 1.285e-14
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.435140440      -0.008771930       0.003395318

Como se puede apreciar el valor del estadístico de moran es de 0.43, además, el P_VALUE es inferior a una significancia del 5%, por lo tanto se tiene autocorrelación espacial.

moran.mc(CUNDI@data$TOTAL_TRANS, queen, nsim=999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  CUNDI@data$TOTAL_TRANS 
## weights: queen  
## number of simulations + 1: 1000 
## 
## statistic = 0.43514, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

C de Geary

geary.test(CUNDI@data$TOTAL_TRANS, queen)

## 
##  Geary C test under randomisation
## 
## data:  CUNDI@data$TOTAL_TRANS 
## weights: queen 
## 
## Geary C statistic standard deviate = 6.2555, p-value = 1.981e-10
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.60135355        1.00000000        0.00406115

I de moran local

El Índice Local de Asociación Espacial (LISA), permite la identificación de patrones locales de asociación espacial, descomponiendo el Índice Moran para evaluar la influencia de ubicaciones individuales en la estadística global que amplía las capacidades de visualización de los valores analizados a través del uso de Sistemas de Información Geográfica.

Este índice se encarga de representar aquellas localizaciones con valores significativos en indicadores estadísticos de asociación espacial local, alertando así de la presencia de puntos calientes hot spots o atípicos espaciales, cuya intensidad depende de la significativa asociada de los datos estadísticos analizados.

locali<-localmoran(CUNDI@data$TOTAL_TRANS, queen)
library(dplyr)
CUNDI@data <- CUNDI@data %>%
              mutate(localmi = locali[,1], localz = locali[,4])

CUNDI@data <- CUNDI@data %>%
                mutate(mcluster = cut(localz, breaks = c(min(localz),-1.96, 1.96, max(localz)), include.lowest = TRUE, labels = c("Negative Correlation", "Not Significant", "Positive Correlation")))

Para esto se estableció la significancia del LISA que se encuentra dentro de los límites de la distribución normal con valor \(z\sim+-1.96\), así pues, como se muestra en el siguiente mapa.

factpal <- colorFactor(heat.colors(3), CUNDI@data$mcluster)


state_popup <- paste0("<strong>Municipio: </strong>", 
                      CUNDI@data$MPIO, 
                      "<br><strong>Significancia: </strong>", 
                      CUNDI@data$mcluster)

leaflet(CUNDI)%>% 
  #MAPA BASE
  addTiles()%>%addProviderTiles(providers$Esri.NatGeoWorldMap)%>%
  addPolygons(color = "#444444", weight = 1, smoothFactor = 0.5,
              popup = state_popup,
              opacity = 1.0, fillOpacity = 0.5,
              fillColor =~colorQuantile("YlOrRd", TOTAL)(TOTAL),
              highlightOptions = highlightOptions(color = "white", weight = 2,bringToFront = TRUE),group = "Hurtos")%>%#CAPAS
  addPolygons(weight = 1, smoothFactor = 0.5,
              popup = state_popup,
              opacity = 1.0, fillOpacity = 0.5,
              color = ~factpal(mcluster))%>%
  #CONTROL CAPAS
  addLayersControl(
    #baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
    overlayGroups = c("I local","Hurtos"),
    options = layersControlOptions(collapsed = FALSE)
  ) 

Si observamos el mapa previo, podemos observar que la Sabana de Bogotá cuenta con un cluster de significancia, donde existe correlación positiva o hot Spots, mientras que al norte de la sabana podemos evidenciar que el municipio de Pacho tiene una correlación negativa. Además existen ciertas conglomeraciones a los extremos que representan asociación espacial positiva, allí se mantiene una relación de bajos rodeados de bajos.

Modelamiento Espacial

El último paso después de haber encoentrado autocorrelación espacial, es modelar esta, en función de la variable objeto de estudio previamente transformada. Para esto se tiene en cuenta dos modelos, uno refrente al retraso espacial, que hace referencia cuando las observaciones vecinas se afectan entre sí. Tales dependencias pueden conducir a estimaciones inconsistentes y sesgadas en un modelo OLS (previamente calculado). E incluso si no le importa el “espacio” en un sentido geográfico.

modelo_lags<-lagsarlm(CUNDI@data$TOTAL_TRANS~log(CUNDI@data$TOTALPOB)+log(CUNDI@data$PROP_NBI)+
                        log(CUNDI@data$PROP_MISER),listw=queen)
summary(modelo_lags)

## 
## Call:lagsarlm(formula = CUNDI@data$TOTAL_TRANS ~ log(CUNDI@data$TOTALPOB) + 
##     log(CUNDI@data$PROP_NBI) + log(CUNDI@data$PROP_MISER), listw = queen)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.452471 -0.175381  0.092225  0.242778  0.622677 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                -1.72262    0.88079 -1.9558  0.05049
## log(CUNDI@data$TOTALPOB)    0.52943    0.04460 11.8705  < 2e-16
## log(CUNDI@data$PROP_NBI)   -0.65271    0.25535 -2.5562  0.01058
## log(CUNDI@data$PROP_MISER)  0.26953    0.12398  2.1740  0.02971
## 
## Rho: 0.27219, LR test value: 11.926, p-value: 0.00055347
## Asymptotic standard error: 0.067849
##     z-value: 4.0116, p-value: 6.0301e-05
## Wald statistic: 16.093, p-value: 6.0301e-05
## 
## Log likelihood: -65.10448 for lag model
## ML residual variance (sigma squared): 0.17891, (sigma: 0.42298)
## Number of observations: 115 
## Number of parameters estimated: 6 
## AIC: 142.21, (AIC for lm: 152.14)
## LM test for residual autocorrelation
## test value: 1.2202, p-value: 0.26932

Dentro de este modelo podemos observar que todas las variables resultan ser significativas al 5% de confiabilidad. Allí se puede mencionar que un aumento del 1% del total de la población aumenta en 0.52% la variable transformada bajo box-cox. Además, esta se encuentra relacionada de forma directamente proporcional con la variable de interés.

En cuanto a la variable de proporción de NBI, un aumento de esta en 1%, representa la disminución de la variable transformada en 0.65%. Así mismo, se puede aseverar que hay una relación inversamente proporcional entre esta y la variable de interés.

Finalmente para la variable independiente de proporción de miseria, se tiene que un aumento de esta en 1%, aumenta en 0.29% la variable transformada.

El modelo surge de la presencia de dependencia espacial en el término de error de una unidad espacial y las correspondientes unidades vecinas. Ocurre en el caso de que algunas variables influyan valor de la variable dependiente pero excluyendo en el modelo la correlación entre las unidades espaciales.

modelo_error <- errorsarlm(CUNDI@data$TOTAL_TRANS~log(CUNDI@data$TOTALPOB)+log(CUNDI@data$PROP_NBI)+
                             log(CUNDI@data$PROP_MISER),listw=queen)
summary(modelo_error)

## 
## Call:errorsarlm(formula = CUNDI@data$TOTAL_TRANS ~ log(CUNDI@data$TOTALPOB) + 
##     log(CUNDI@data$PROP_NBI) + log(CUNDI@data$PROP_MISER), listw = queen)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.716106 -0.179537  0.077197  0.306379  0.670537 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                            Estimate Std. Error z value  Pr(>|z|)
## (Intercept)                -0.87541    0.91446 -0.9573  0.338422
## log(CUNDI@data$TOTALPOB)    0.53044    0.04586 11.5666 < 2.2e-16
## log(CUNDI@data$PROP_NBI)   -0.72796    0.27408 -2.6560  0.007907
## log(CUNDI@data$PROP_MISER)  0.20976    0.13329  1.5737  0.115564
## 
## Lambda: 0.33438, LR test value: 5.9933, p-value: 0.01436
## Asymptotic standard error: 0.12283
##     z-value: 2.7222, p-value: 0.006484
## Wald statistic: 7.4106, p-value: 0.006484
## 
## Log likelihood: -68.07098 for error model
## ML residual variance (sigma squared): 0.18683, (sigma: 0.43224)
## Number of observations: 115 
## Number of parameters estimated: 6 
## AIC: 148.14, (AIC for lm: 152.14)

Para este modelo se tienen solamente las variables de población y proporción de necesidades básicas insatisfechas como, variables significativas al 5%. Se mantienen las relaciones de dependencia al igual que el modelo SPATIAL LAG, la población esta relacionada directamente y la proporción de NBI negativamente.

anova(modelo_lags,modelo_error)            

##              Model df    AIC  logLik
## modelo_lags      1  6 142.21 -65.104
## modelo_error     2  6 148.14 -68.071

Como podemos observar a partir del críterio del AIC se puede inferir que el modelo SPATIAL LAG es mucho mejor que el modelo SPATIAL ERROR dado que cuenta con un menor valor del críterio. Además, este modelo contempla las tres variables de interés y resultan ser significativas al 5% de confiabilidad.