Ayudantía 01 - Spatial Analytics

Cargamos librerías y bases

library(dplyr)
library(sf)
library(rmapshaper)
library(spdep)
library(chilemapas)
library(data.table)
library(ggplot2)
library(ggpubr)
library(classInt)

path <- 'C:/Users/masie/Dropbox/01 Universidad/13 Extras/2021 S2 Spatial Analytics/Ayudantia 01/'
censo <- fread(paste0(path,"Censo2017_Personas.csv"))
mapa <- readRDS(paste0(path,'zc_censo2017.rds'))

Tasa de la ocupación por sector

Creamos primero el código único identificador de cada zona censal.

censo <- censo %>% 
  mutate(CODIGO_COMUNA   = COMUNA * 10 ^ 6, 
         DISTRITO_CENSAL = DC * 10 ^ 4,
         CODIGO_AREA     = AREA * 10 ^ 3,
         ZONA = as.numeric(CODIGO_COMUNA + DISTRITO_CENSAL + CODIGO_AREA + ZC_LOC))

Posteriormente, calculamos el número de ocupados de zona urbana de la Región Metropolitanala (Santiago) según sector (primario, secundario y tericario/servicios).

censo <- censo %>% 
  filter(P17 %in% c(1,2) & !P18 %in% c(NA,'99','98','Z') & P09 >= 15) %>% 
  mutate(SECTOR := ifelse(P18 %in% c('A','B'),'SECTOR PRIMARIO',
                          ifelse(P18 %in% c('C','D','E','F'),'SECTOR SECUNDARIO',
                                 'SECTOR TERCIARIO'))) %>%
  group_by(ZONA,SECTOR) %>% 
  summarise(OCUPADOS = n())

Para calcular la tasa de participación por sector y zona censal, es necesario transformar la base de datos de formato long a wide. Posible hacerlo con funciones como reshape, spread o dcast.

censo_w <- dcast.data.table(data.table(censo),ZONA ~ SECTOR,value.var = 'OCUPADOS',fill = 0)

Unimos los datos de polígonos de zonas censales a los datos obtenidos de la tasa de ocupación por sector. Centraemos el estudio en la Región Metropolitana.

mapa <- mapa %>% 
  filter(REGION == '13') %>% 
  mutate(provincia = trunc(COD_INE_15/100000000-130)) %>% 
  filter(provincia == 1 & AREA == 1 & 
         !COD_INE_16 %in% c(13119131001,13124071001,13124071002,13124071003,
                            13124071004,13124071005,13124081001)) %>%
  mutate(COD_INE_16 := as.numeric(COD_INE_16))

censo_w <- merge(censo_w,mapa,by.x='ZONA',by.y='COD_INE_16',all.y=T) %>%
  filter(!is.na(AREA) & !is.na(TOTAL))

Observemos la distribución de los valores por sector productivo:

censo %>% 
  group_by(ZONA) %>% 
  filter(ZONA %in% censo_w$COD_INE_15) %>%
  mutate(PORCENTAJE = OCUPADOS/sum(OCUPADOS,na.rm = T)*100) %>%
ggplot(aes(x=PORCENTAJE,color=SECTOR,group=SECTOR)) +
  geom_density() +
  theme_minimal() +
  theme(legend.position = 'bottom') +
  facet_wrap(~SECTOR,ncol = 3,scales = 'free')

Observamos que los sectores secundiario y terciario son los que poseen mayor dispersión, lo cual pudiese ser un indicio de mayor heterogeneidad espacial.

Mapeo de la información

A continuación, obsrvaremos la distribución espacial de la proporción de ocupados en el sector terciario.

censo2 <- censo %>% 
  group_by(ZONA) %>% 
  mutate(PORCENTAJE = OCUPADOS/sum(OCUPADOS,na.rm = T)*100) %>%
  filter(SECTOR == 'SECTOR TERCIARIO')

ggplot(censo_w) +
  geom_sf(aes(geometry=geometry,fill=`SECTOR TERCIARIO`)) +
  theme_void() +
  theme_minimal() +
  theme(legend.position = 'bottom')

En general, no podemos apreciar mucho de un mapa con una variable continua como “fill”. En consecuencia, es posible modificar dichos mapas. Supongamos que en primera instancia deseamos ver cuanto se alejan los valores de la media o mediana. Observemos primero dichos valores.

q <- quantile(censo_w$`SECTOR TERCIARIO`,na.rm = T,probs = 0.5)
q

FALSE      50% 
FALSE 82.35662

m1 <- mean(censo_w$`SECTOR TERCIARIO`,na.rm = T)
m1

FALSE [1] 81.86938

ggplot(censo_w) +
  geom_sf(aes(geometry=geometry,fill=`SECTOR TERCIARIO`)) +
  theme_void() +
  theme_minimal() +
  scale_fill_gradient2(low='red',mid = 'yellow',high = 'green',
                       midpoint = q) +
  labs(fill=NULL, title = 'Proporción de ocupados en el sector terciario',
       subtitle = paste0('Punto medio en la mediana (',round(q,2),'%)'))

ggplot(censo_w) +
  geom_sf(aes(geometry=geometry,fill=`SECTOR TERCIARIO`)) +
  theme_void() +
  theme_minimal() +
  scale_fill_gradient2(low='red',mid = 'yellow',high = 'green',
                       midpoint = m1) +
  labs(fill=NULL, title = 'Proporción de ocupados en el sector terciario',
       subtitle = paste0('Punto medio en la media (',round(m1,2),'%)'))

Si bien es posible disernir de mejor forma que sectores poseen en mayor o menor intensidad ocupados en el sector terciario, las escalas continuas pueden seguir complicando la interpretación. Pudiese ser mejor agrupar los valores de acuerso a ciertos rangos. Veamos primero según cuartiles, de forma de ver que entre que rangos de tasas se encuentra cada 25% de las zonas censales:

breaks_quantile01 <- classIntervals(censo_w$`SECTOR TERCIARIO`, style = 'quantile',n=4)
censo_w$breaks01 <- cut(censo_w$`SECTOR TERCIARIO`,breaks = breaks_quantile01$brks,
                        include.lowest = T)

ggplot(censo_w) +
  geom_sf(aes(geometry=geometry,fill=breaks01)) +
  theme_void() +
  theme_minimal() +
  scale_fill_viridis_d() +
  labs(fill=NULL, title = 'Proporción de ocupados en el sector terciario',
       subtitle = paste0('Punto medio en la mediana (',round(q,2),'%)'))

Una observación notoria es que un 25% de las zonas censales poseen una tasa de participación en el sector servicios entre el 87.2% al 100%, los cuales se encuentran casi en su totalidad en la zona oriente de la capital. A continuació, repetimos el ejercicos con desviaciones estandar, de forma de observar que tanto se alejan los valores atípicos de la media.

min <- min(censo_w$`SECTOR TERCIARIO`,na.rm = T)
max <- max(censo_w$`SECTOR TERCIARIO`,na.rm = T)
std <- sd(censo_w$`SECTOR TERCIARIO`,na.rm = T)
std.interval01 <- m1+std*seq(-8,4,1)

std.interval01

FALSE  [1]  31.07217  37.42182  43.77147  50.12112  56.47077  62.82043  69.17008
FALSE  [8]  75.51973  81.86938  88.21904  94.56869 100.91834 107.26799

Y en torno a la mediana.

std <- mad(censo_w$`SECTOR TERCIARIO`,na.rm = T)
std.interval02 <- q+std*seq(-8,4,1)

std.interval02

FALSE  [1]  23.38553  30.75692  38.12830  45.49969  52.87108  60.24246  67.61385
FALSE  [8]  74.98523  82.35662  89.72801  97.09939 104.47078 111.84217

Utilizamos classIntervals para crear los intervalor de desviaciones estandar en torno a la media (comparamos lo obtenido en vector std.interval01).

breaks_quantile02 <- classIntervals(censo_w$`SECTOR TERCIARIO`, style = 'sd')
censo_w$breaks02 <- cut(censo_w$`SECTOR TERCIARIO`,breaks = breaks_quantile02$brks, 
                        include.lowest = T)

ggplot(censo_w) +
  geom_sf(aes(geometry=geometry,fill=breaks02)) +
  theme_void() +
  theme_minimal() +
  scale_fill_viridis_d() +
  labs(fill=NULL, title = 'Proporción de ocupados en el sector terciario',
       subtitle = paste0('Desviaciones estandar en torno a la media (',round(m1,2),'%)'))

Repetimos, pero ahora para las desviaciones en torno a la mediana.

breaks_quantile03 <- classIntervals(censo_w$`SECTOR TERCIARIO`, style = 'fixed',
                                    fixedBreaks = std.interval02)
censo_w$breaks03 <- cut(censo_w$`SECTOR TERCIARIO`,breaks = breaks_quantile03$brks, 
                        include.lowest = T)

ggplot(censo_w) +
  geom_sf(aes(geometry=geometry,fill=breaks03)) +
  theme_void() +
  theme_minimal() +
  scale_fill_viridis_d() +
  labs(fill=NULL, title = 'Proporción de ocupados en el sector terciario',
       subtitle = paste0('Desviaciones estandar en torno a la mediana (',round(q,2),'%)')) 

Observamos que existen indicios de ciertos sectores con alta y baja tasa de participación en sector servicios.

Índice global de Morán

Para profundizar en el análisis estudiemos la posibilidad de que existan clusters espaciales con altas y bajas concentración de ocupados en el sector servicios (tasa de participación). Para esto, se comparará el valor de la tasa de cada zona censal con la de sus rezagos espaciales (o vecinos más cercanos). Primero veamos quienes son los vecinos más cercanos de cada zona censal.

crs_planar <- "+proj=longlat +datum=WGS84 +no_defs"
crs_utm <- "+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

censo_sf <- censo_w %>%
  st_as_sf() %>% 
  st_transform(crs_utm) %>% 
  as_Spatial()

censo_queen <- poly2nb(censo_sf)

plot(censo_sf, border = 'lightgrey')
plot(censo_queen, coordinates(censo_sf), add=TRUE, col='red')

Observamos visualmente la relación entre la tasa de la zona censal y el promedio de sus vecinos.

mor <- moran.plot(censo_w$`SECTOR TERCIARIO`,
                  listw = nb2listw(censo_queen, style = "W",zero.policy = T),
                  zero.policy = T,plot = T)

Se aprecia una relación positiva entre estas. Es decir, en la medida que aumenta en las zonas censales el valor de su tasa de ocpación del sector servicios, aumenta (en menor medida) la tasa de sus vecinos. ¿En que medida?

list_q <- nb2listw(censo_queen,zero.policy = T)
globalMoran_q <- moran.test(censo_w$`SECTOR TERCIARIO`,list_q,zero.policy = T)
globalMoran_q[["estimate"]][["Moran I statistic"]]

FALSE [1] 0.6811449

Observemos el ïndice Global de Morán para la tasa de participación de los ocupados en el sector servicios. Testearemos si la tasa se encuetra distribuida aleatoriamente en el espacio (H0) o si existe autocorrelación espacial (H1).

globalMoran_q

FALSE 
FALSE   Moran I test under randomisation
FALSE 
FALSE data:  censo_w$`SECTOR TERCIARIO`  
FALSE weights: list_q    
FALSE 
FALSE Moran I statistic standard deviate = 44.599, p-value <
FALSE 0.00000000000000022
FALSE alternative hypothesis: greater
FALSE sample estimates:
FALSE Moran I statistic       Expectation          Variance 
FALSE      0.6811449322     -0.0007017544      0.0002337364

Observamos que el valor p es significativamente menor a 0.05, lo cual nos permite rechazar con un 95% de confianza la hipótesis nula. Lo anterior supone que existen indicios de autocorrelación espacial de la tasa de ocupación para el sector servicios entre zonas censales. Luego, extraemos la pendiente que indica la relación entre el valor de la tasa de la zona censal y el lag espacial (valores de vecinos más cercanos).

Índice local de Morán

Dado que existen indicios de existencia de clusters, procedemos a hacer un análisis LISA. El índice local de moral es un estadísticos de autocorrelación espacial a nivel local. En este caso, consdieramos el valor p de cada zona censal y la ubicación en el gráfico de Morán visto peviamente (si los valoes de la tasa local y de vecinos se encuentran por sobre o debajo de la media).

local <- localmoran(x = mor$x,
                    listw = nb2listw(censo_queen,style = "W",zero.policy = T),
                    zero.policy = T)

Posterior a la estimación, indicamos que zonas censales cumplen con los criterios de tolerancia deseados para designar posibles cluster espaciales. Primero veamos el caso de un nivel de significancia del 90%.

quadrant <- vector(mode="numeric",length=nrow(local))

tasa_loc <- mor$x - mean(mor$x,na.rm = T) #Promedio zonas censales
tasa_vec  <- mor$wx - mean(mor$wx,na.rm = T) #Promedio de los vecinos

signif <- 0.1

quadrant[tasa_vec >0 & tasa_loc>0] <- 4 #High - High
quadrant[tasa_vec <0 & tasa_loc<0] <- 1 #Low - Low
quadrant[tasa_vec <0 & tasa_loc>0] <- 2 #High - Low
quadrant[tasa_vec >0 & tasa_loc<0] <- 3 #Low - High
quadrant[local[,5]>signif] <- 0 #P-value > 0.1

censo_w <- cbind(censo_w, quadrant=as.character(quadrant))

Procedemos a repetir el gráfico de Morán, incluyendo la información de aquellas zonas censales que son signifcativas.

moran <- data.frame(cbind(loc = mor$x,
                          vec = mor$wx,
                          Significancia=as.character(quadrant)))

ggplot(moran,aes(x=as.numeric(loc),y=as.numeric(vec),color=Significancia)) +
  geom_vline(aes(xintercept = mean(as.numeric(loc),na.rm = T))) +
  geom_hline(aes(yintercept = mean(as.numeric(vec),na.rm = T))) +
  geom_point(size=3,aes(shape=Significancia)) +
  geom_smooth(method="lm", formula =y ~ x, se=F, col="red",)+
  theme_minimal() +
  theme(legend.position = 'bottom') +
  labs(x='Valor Local',y='Valor promedio vecinos más cercanos',title = 'Grafico de Morán') +
  scale_color_manual(name = "", labels = c('0'="NS",'1'="LL",'2'="HL",'3'="LH",'4'="HH"),
                     values = c("0" = 'black', "1"="blue", "2"= rgb(1,0,0,alpha=0.4), 
                                "3"=rgb(0,0,1,alpha=0.4), "4"= "red")) +
  scale_shape_manual(name = "",labels = c('0'="NS",'1'="LL",'2'="HL",'3'="LH",'4'="HH"),
                     values = c("0" = 1, "1"=16, "2"= 16,"3"=16, "4"= 16))

Notamos que las zonas censales que se ubican cerca de las medias (de valores propios y de los vecinos) no son significativos. A continuación, observamos la distribución espacial de las zonas censales significativas. Lo anterior permite encontrar clusters con altos o bajos niveles de tasa de participación del sector terciario en la ciudad de Santiago.

mapa01 <-
ggplot(censo_w) +
  geom_sf(aes(fill = quadrant,geometry=geometry)) +
  scale_fill_manual(name = "",  labels = c('0'="NS",'1'="LL",'2'="HL",'3'="LH",'4'="HH"),
                    values = c("0" = "white", "1"="blue", "2"= rgb(1,0,0,alpha=0.4), 
                               "3"=rgb(0,0,1,alpha=0.4), "4"= "red")) +
  theme_void() +
  theme(legend.position = 'bottom')

mapa01 + labs(title = 'Indice de Moral Local (p-value <= 0.1)', 
              subtitle = 'Tasa de participación del sector servicios en la ocuapción')

Repitamos el ejercicios nuevamente, aumentando el nivel se significancia a un 95%.

quadrant2 <- vector(mode="numeric",length=nrow(local))

signif <- 0.05

quadrant2[tasa_vec >0 & tasa_loc>0] <- 4 #High - High
quadrant2[tasa_vec <0 & tasa_loc<0] <- 1 #Low - Low
quadrant2[tasa_vec <0 & tasa_loc>0] <- 2 #High - Low
quadrant2[tasa_vec >0 & tasa_loc<0] <- 3 #Low - High
quadrant2[local[,5]>signif] <- 0 #P-value > 0.1

censo_w <- cbind(censo_w, quadrant2=as.character(quadrant2))

Volvemos a observar la distribución espacial de los clusters bajo este nuevo nivel de exigencia, de forma de observa que sectores ven reducidos el número de zonas censales significantes.

mapa02 <-
ggplot(censo_w) +
  geom_sf(aes(fill = quadrant2,geometry=geometry)) +
  scale_fill_manual(name = "",  labels = c('0'="NS",'1'="LL",'2'="HL",'3'="LH",'4'="HH"),
                    values = c("0" = "white", "1"="blue", "2"= rgb(1,0,0,alpha=0.4), 
                               "3"=rgb(0,0,1,alpha=0.4), "4"= "red")) +
  theme_void() +
  theme(legend.position = 'bottom')

mapa02 + labs(title = 'Indice de Moral Local (p-value <= 0.05)', 
              subtitle = 'Tasa de participación del sector servicios en la ocuapción')

Repitamos por última vez el ejercicio, aumentando ahora el nivel se significancia a un 99%.

quadrant3 <- vector(mode="numeric",length=nrow(local))

signif <- 0.01

quadrant3[tasa_vec >0 & tasa_loc>0] <- 4 #High - High
quadrant3[tasa_vec <0 & tasa_loc<0] <- 1 #Low - Low
quadrant3[tasa_vec <0 & tasa_loc>0] <- 2 #High - Low
quadrant3[tasa_vec >0 & tasa_loc<0] <- 3 #Low - High
quadrant3[local[,5]>signif] <- 0 #P-value > 0.1

censo_w <- cbind(censo_w, quadrant3=as.character(quadrant3))

Volvemos a obsevar la distribución espacial de los clusters.

mapa03 <-
ggplot(censo_w) +
  geom_sf(aes(fill = quadrant3,geometry=geometry)) +
  scale_fill_manual(name = "",  labels = c('0'="NS",'1'="LL",'2'="HL",'3'="LH",'4'="HH"),
                    values = c("0" = "white", "1"="blue", "2"= rgb(1,0,0,alpha=0.4), 
                               "3"=rgb(0,0,1,alpha=0.4), "4"= "red")) +
  theme_void() +
  theme(legend.position = 'bottom') 

mapa03 + labs(title = 'Indice de Moral Local (p-value <= 0.01)', 
              subtitle = 'Tasa de participación del sector servicios en la ocuapción')

Por último, agrupamos los tres mapas previos para observar los cambios de mejor forma.

ggarrange(mapa01,mapa02,mapa03,ncol = 3,common.legend = T,
          labels = c('p-value<=0.1','p-value<=0.05','p-value<=0.01'),
          legend = 'bottom')