Actividad
Evalúa el porcentaje de viviendas de cada AGEB de Cancún con acceso a internet.
df <- read_excel("file_show.xls") %>%
filter(str_detect(MUN, "005"), str_detect(LOC, "0001"), str_detect(NOM_LOC, "Total AGEB urbana")) %>%
dplyr::select(ENTIDAD:AGEB, VIVPAR_HAB, VPH_INTER) %>%
mutate_all(~na_if(., '*'))
df$VIVPAR_HAB[is.na(df$VIVPAR_HAB)] <- 0
df$VIVPAR_HAB <- as.double(df$VIVPAR_HAB)
df$VPH_INTER[is.na(df$VPH_INTER)] <- 0
df$VPH_INTER <- as.double(df$VPH_INTER)
shapefile <- st_read("AGEB_QROO/AGEBS_QROO.shp") %>%
dplyr::select(CVEGEO) %>%
separate(col = CVEGEO, into = c("ENTIDAD", "MUN", "LOC", "AGEB"), sep = c(2, 5, 9)) %>%
filter( MUN == "005" & LOC == "0001")
## Reading layer `AGEBS_QROO' from data source `C:\Users\dovo_\Desktop\Sistemas de informacion geografica\Autocorrelacion espacial\AGEB_QROO\AGEBS_QROO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 812 features and 7 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 3873704 ymin: 762271.2 xmax: 4076865 ymax: 1111848
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs
head(shapefile)
## Simple feature collection with 6 features and 4 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 4057062 ymin: 1091234 xmax: 4065441 ymax: 1095606
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs
## ENTIDAD MUN LOC AGEB geometry
## 1 23 005 0001 3820 POLYGON ((4057319 1092664, ...
## 2 23 005 0001 3816 POLYGON ((4057319 1092664, ...
## 3 23 005 0001 3801 POLYGON ((4057217 1093987, ...
## 4 23 005 0001 3746 POLYGON ((4064699 1094222, ...
## 5 23 005 0001 3661 POLYGON ((4063912 1094215, ...
## 6 23 005 0001 3731 POLYGON ((4064418 1094963, ...
percent_agebs = df$VPH_INTER / df$VIVPAR_HAB * 100
percent_agebs[is.na(percent_agebs)] <- 0
Agebs_Shape <- inner_join(shapefile, df) %>%
st_geometry()
shades <- auto.shading(percent_agebs, n = 7, cols = brewer.pal(7, "Blues"))
choropleth(Agebs_Shape, v = percent_agebs, shades, main="Viviendas de cada AGEB de Cancún con acceso a internet")
Evalúa la media rezagada y construye un mapa para visualizarla.
agebs_nb <- poly2nb(Agebs_Shape, queen=F)
summary(agebs_nb)
## Neighbour list object:
## Number of regions: 309
## Number of nonzero links: 1552
## Percentage nonzero weights: 1.625454
## Average number of links: 5.022654
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12
## 14 25 28 42 70 72 32 14 6 3 2 1
## 14 least connected regions:
## 1 56 72 73 74 77 81 82 84 230 231 290 300 303 with 1 link
## 1 most connected region:
## 136 with 12 links
agebs_nb <- poly2nb(Agebs_Shape)
summary(agebs_nb)
## Neighbour list object:
## Number of regions: 309
## Number of nonzero links: 1618
## Percentage nonzero weights: 1.694578
## Average number of links: 5.236246
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12
## 14 25 25 38 56 73 42 19 11 3 2 1
## 14 least connected regions:
## 1 56 72 73 74 77 81 82 84 230 231 290 300 303 with 1 link
## 1 most connected region:
## 136 with 12 links
plot(Agebs_Shape, border='lightgrey')
plot(agebs_nb, coordinates(as(Agebs_Shape, "Spatial")) , add=T, col='green')
agebs_lw <- nb2listw(agebs_nb)
summary(agebs_lw)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 309
## Number of nonzero links: 1618
## Percentage nonzero weights: 1.694578
## Average number of links: 5.236246
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12
## 14 25 25 38 56 73 42 19 11 3 2 1
## 14 least connected regions:
## 1 56 72 73 74 77 81 82 84 230 231 290 300 303 with 1 link
## 1 most connected region:
## 136 with 12 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 309 95481 309 143.1795 1284.992
weights_l <- lag.listw(agebs_lw, percent_agebs)
choropleth(Agebs_Shape, percent_agebs, shades, main="viviendas de cada AGEB de Cancún con acceso a internet")
choropleth(Agebs_Shape, weights_l, shades, main="Media rezagada")
plot(percent_agebs, weights_l, xlab="Porcentaje de AGEB con Internet", ylab="Media rezagada")
abline(a=0, b=1)
Basándote en la matriz de pesos que se obtiene por medio de spdep (listw2mat), implementa tu evaluación del índice de Moran y prueba con el porcentaje de viviendas con acceso a internet en las AGEBS de Cancún.
P_Moran <- function(Z, Shape){
agebs_nb <- poly2nb(Shape)
agebs_lw <- nb2listw(agebs_nb)
N <- length(agebs_nb)
W_ij <- listw2mat(agebs_lw)
S_Wij <- sum(W_ij)
Zm <- mean(Z)
S_z <- Z - Zm
g <- expand.grid(S_z, S_z)
zizj <- g[,1] * g[,2]
ml <- (N / S_Wij) * sum(W_ij*(zizj)) / (sum((S_z)^2))
return(ml)
}
ml <- P_Moran(percent_agebs, Agebs_Shape)
print(paste("Índice de Moran: ", round(ml, digits=9)))
## [1] "Índice de Moran: 0.052508345"
Compara tus resultados con los obtenidos con la librería spdep.
moran.test(percent_agebs, agebs_lw)
##
## Moran I test under randomisation
##
## data: percent_agebs
## weights: agebs_lw
##
## Moran I statistic standard deviate = 1.4488, p-value = 0.0737
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.052508345 -0.003246753 0.001480988