Actividad 3 Modulo 1, Equipo 4
Actividad 3. Desarrollo de Modelos de Regresión Espacial y Machine Learning
En esta actividad, nos enfocaremos en aplicar un modelo de regresión espacial en R, para así poder identificar cuales son las zonas más afectadas por la enfermedad de COVID-19 y su correlación con la capacidad de encontrar servicios de salud por cada municipio.
(La mayoría del código que se usa para depurar o acomodar datos quedará escondido.)
Diferencias entre modelos de regresión global y modelos regresión local (GWR) • Los modelos de regresión global sí identifican autocorrelación espacial, pero no diferencian entre distintas locaciones dentro de la base de datos, a diferencia del GWR que sí lo hace. • Los modelos de regresión global suelen ser mejor para predecir la variable dependiente que los modelos locales, sin embargo, los locales suelen explicar mejor el comportamiento de las variables independientes. • Los modelos de regresión global no explican las relaciones que existen entre distintas áreas, a diferencia del GWR, que sí explica cómo afecta un área a la otra.
En la primera gráfica de mapa identificamos con puntos las zonas pobladas con focos de infección significativos. El mapa muestra que están distribuidos por todo el país, sin embargo, los focos de contagio se muestran más en los municipios de la zona centro del país.
Indentificando densidad poblacional
En este mapa podemos observar la distribución de las personas sin acceso a servicios de salud en México. No se observan cluster determinantes, máss la zona sur de México parece presentar mayores niveles de falta de acceso.
SWM calculo para conexiones municipales
map2 <- map
map2 <- merge(map2,covid,by.x = "IDUNICO",by.y = "cve_ent")
map2$CODELAG <- NULL
map2$CVE_ENT <- NULL
map2$IDUNICO <- NULL
swm_queen <- poly2nb(map_deprctd, queen = TRUE)## although coordinates are longitude/latitude, st_intersects assumes that they
## are planarswm_rook <- poly2nb(map_deprctd, queen = FALSE)## although coordinates are longitude/latitude, st_intersects assumes that they
## are planarQueen Contiguity graphic
### plotting queen contiguity based neighbors maps
plot(map_deprctd, borders = 'lightgrey')
plot(swm_queen, coordinates(map_deprctd), pch = 19, cex = 0.6, add = TRUE, col = "red", xlim =c(0,5), ylim = c(0,2))
title(main = "Queen Contiguity", cex.main = 0.9)
En esta gráfica podemos ver la relación desde los centroides en cada municipio con sus vecinos con un tipo de relación especificada Queen, la cual permite que los municipios desarrollen vecinos en una matrix de conectividad no solo horizontal y verticalmente, sino que también diagonalmente.
Rook Contiguity Graphic
### plotting rook contiguity based neighbors maps
plot(map_deprctd, borders = 'lightgrey')
plot(swm_rook, coordinates(map_deprctd), pch = 19, cex = 0.6, add = TRUE, col = "red")
title(main = "Rook Contiguity", cex.main = 0.9)
Es el mismo tipo de gráfica que la anterior, pero esta vez con una especificación Rook, por lo que desarrolla vecinos unicamente de manera horizontal y vertical. Vale la pena destacar que como hay demasiados municipios, no se puede denotar mucho la correlación en los mapas. Para un análisis más a profundidad, valdría la pena preparar
Obteniendo las coordenadas para los mapas
### computing distance based neighbours
coords <- coordinates(map_deprctd)
head(coords)## [,1] [,2]
## 1001 -102.2958 21.81144
## 1002 -102.0456 22.12651
## 1003 -102.7049 21.90069
## 1004 -102.2970 22.36063
## 1005 -102.4456 21.93212
## 1006 -102.3017 22.10454knn1 <- knn2nb(knearneigh(coords))
knn1_dist <- unlist(nbdists(knn1, coords, longlat = TRUE))Neighbor Links Analysis
dwm <- dnearneigh(coords, 0 ,86, longlat = TRUE)
dwm## Neighbour list object:
## Number of regions: 2456
## Number of nonzero links: 237042
## Percentage nonzero weights: 3.929783
## Average number of links: 96.51547
## 12 regions with no links:
## 12 13 18 19 20 31 34 55 66 203 1807 1943En promedio, la cantidad de vecinos que hay para cada municipio son 96. Además, se encontraron 12 municipios sin conexiones. Todo esto tomando en cuenta una distancia de medición de 86 KM a la redonda.
Mapa: Vecinos a 86 KM
plot(map_deprctd, border = "lightgrey")
plot(dwm, coords, add = TRUE, pch = 19, cex = 0.6)
title(main = "Neighbours within 86 km", cex.main = 0.9)
En esta gráfica a diferencia de la queen o la rook, tiene una menor cantidad de relaciones entre municipios en el norte, pues al limitar los vecinos posibles a 86 KM a la redonda, en el norte, municipios más aislados pierden conexiones, mientras que en el sur, en casos como Oaxaca, los municipios ganan vecinos.
Weighted Queen Neighbors
rswm_queen <- nb2listw(swm_queen, style = "W", zero.policy = TRUE)
rswm_queen## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 2456
## Number of nonzero links: 14392
## Percentage nonzero weights: 0.2385967
## Average number of links: 5.859935
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 2456 6031936 2456 898.7549 10338.48Siguiendo con la lógica anterior, ahora podemos observar como efectivamente, los vecinos creados con especificación queen, son menos que los delimitados por una distancia, pues en este resultado se puede ver como solo hay en promedio 5 vecinos para cada municipio. Cada aclarar que este resultado también incluye el factor weighted. Los pesos aquí creados nos servirán después para hacer el análisis GWR
La parte para predecir la variable principal
# Data cleanup, después lo tengo que mandar llamar otra vez
### lets create a spatial lag of dataset's variable
merged_data$porcentaje_pob_servicios_salud[is.na(merged_data$porcentaje_pob_servicios_salud)] <- 0
merged_data$popden2020[is.na(merged_data$popden2020)] <- 0
merged_data$porcentaje_pob_pobreza[is.na(merged_data$porcentaje_pob_pobreza)] <- 0
merged_data$porcentaje_pob_pobreza_ext[is.na(merged_data$porcentaje_pob_pobreza_ext)] <- 0
merged_data$porcentaje_pob_acceso_ss[is.na(merged_data$porcentaje_pob_acceso_ss)] <- 0
merged_data$GEO_porcentaje_pob_servicios_salud<-lag.listw(rswm_queen,merged_data$porcentaje_pob_servicios_salud,zero.policy=TRUE)
serviciosSalud <- qtm(merged_data, "GEO_porcentaje_pob_servicios_salud")
spatial_lag_serviciosSalud <- qtm(merged_data, "GEO_porcentaje_pob_servicios_salud")Moran test with Weighted values
# Moran's I Test
moran.test(merged_data$porcentaje_pob_servicios_salud, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)##
## Moran I test under randomisation
##
## data: merged_data$porcentaje_pob_servicios_salud
## weights: rswm_queen
##
## Moran I statistic standard deviate = 13.311, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.1618600745 -0.0004073320 0.0001486153Haciendo un análisis de autocorrelación espacial global, pero incluyendo nuestros valores ahora weighted, podemos ver que el resultado es que hay autocorrelación espacial positiva con significancia estadística.
Graphic: Moran’s I correlogram
swm_queen <- poly2nb(map_deprctd, queen = TRUE)
Moran_Correlogram <- sp.correlogram(swm_queen, merged_data$porcentaje_pob_servicios_salud, order = 6, method = "I", style = "B")
plot(Moran_Correlogram)
Por último, se confirma que hay correlación en la dependencia espacial de nuestras variables, y no hay outliers significativos, por lo que podemos continuar.
merged_data$crimen_2018[merged_data$crimen_2018 ==0] <- 0.0001
merged_data$crimen_2019[merged_data$crimen_2019 ==0] <- 0.0001
merged_data$porcentaje_pob_servicios_salud[merged_data$porcentaje_pob_servicios_salud ==0] <- 0.0001
merged_data$popden2020 <- as.numeric(merged_data$popden2020)## Warning: NAs introducidos por coerciónmerged_data$porcentaje_pob_pobreza <- as.numeric(merged_data$porcentaje_pob_pobreza)## Warning: NAs introducidos por coerciónmerged_data$porcentaje_pob_pobreza_ext <- as.numeric(merged_data$porcentaje_pob_pobreza_ext)## Warning: NAs introducidos por coerciónmerged_data$rezago_social <- as.numeric(merged_data$rezago_social)
merged_data$`Entidad federativa` <- as.factor(merged_data$`Entidad federativa`)
merged_data$porcentaje_pob_servicios_salud[is.na(merged_data$porcentaje_pob_servicios_salud)] <- 0
merged_data$popden2020[is.na(merged_data$popden2020)] <- 0
merged_data$porcentaje_pob_pobreza[is.na(merged_data$porcentaje_pob_pobreza)] <- 0
merged_data$porcentaje_pob_pobreza_ext[is.na(merged_data$porcentaje_pob_pobreza_ext)] <- 0
merged_data$porcentaje_pob_acceso_ss[is.na(merged_data$porcentaje_pob_acceso_ss)] <- 0#summary(merged_data)
non_spatial_model <- lm(infectados ~ hospitalespHabitantes + popden2020 + crimen_2018 + crimen_2019 + porcentaje_pob_pobreza + porcentaje_pob_pobreza_ext + porcentaje_pob_acceso_ss + rezago_social + `Entidad federativa`, data = merged_data)
summary(non_spatial_model)##
## Call:
## lm(formula = infectados ~ hospitalespHabitantes + popden2020 +
## crimen_2018 + crimen_2019 + porcentaje_pob_pobreza + porcentaje_pob_pobreza_ext +
## porcentaje_pob_acceso_ss + rezago_social + `Entidad federativa`,
## data = merged_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32685 -1255 -118 681 99120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.303e+03 1.724e+03 4.237 2.35e-05 ***
## hospitalespHabitantes 4.140e+01 1.230e+01 3.366 0.000774 ***
## popden2020 2.202e+00 8.114e-02 27.140 < 2e-16 ***
## crimen_2018 3.033e+00 3.750e+00 0.809 0.418820
## crimen_2019 2.940e+00 3.877e+00 0.758 0.448338
## porcentaje_pob_pobreza -1.399e+01 1.385e+01 -1.010 0.312438
## porcentaje_pob_pobreza_ext 8.787e+01 1.853e+01 4.742 2.24e-06 ***
## porcentaje_pob_acceso_ss -9.475e+03 1.154e+03 -8.211 3.52e-16 ***
## rezago_social -7.472e+02 2.475e+02 -3.019 0.002561 **
## `Entidad federativa`02 1.329e+04 2.840e+03 4.678 3.05e-06 ***
## `Entidad federativa`03 8.378e+03 2.832e+03 2.959 0.003119 **
## `Entidad federativa`04 1.728e+02 2.243e+03 0.077 0.938584
## `Entidad federativa`05 -7.756e+02 1.790e+03 -0.433 0.664759
## `Entidad federativa`06 -5.743e+02 2.314e+03 -0.248 0.804013
## `Entidad federativa`07 -9.008e+02 1.683e+03 -0.535 0.592580
## `Entidad federativa`08 -6.868e+02 1.731e+03 -0.397 0.691546
## `Entidad federativa`09 2.901e+04 2.170e+03 13.368 < 2e-16 ***
## `Entidad federativa`10 6.270e+02 1.793e+03 0.350 0.726617
## `Entidad federativa`11 2.101e+03 1.784e+03 1.178 0.238959
## `Entidad federativa`12 -2.852e+02 1.716e+03 -0.166 0.867999
## `Entidad federativa`13 -7.430e+02 1.692e+03 -0.439 0.660703
## `Entidad federativa`14 -1.096e+03 1.657e+03 -0.661 0.508380
## `Entidad federativa`15 -1.759e+03 1.666e+03 -1.056 0.291024
## `Entidad federativa`16 -9.411e+02 1.676e+03 -0.562 0.574507
## `Entidad federativa`17 -1.605e+03 1.849e+03 -0.868 0.385382
## `Entidad federativa`18 -6.361e+02 1.976e+03 -0.322 0.747497
## `Entidad federativa`19 -1.836e+02 1.756e+03 -0.105 0.916724
## `Entidad federativa`20 -1.779e+03 1.630e+03 -1.091 0.275392
## `Entidad federativa`21 -6.800e+02 1.653e+03 -0.411 0.680866
## `Entidad federativa`22 3.729e+03 1.993e+03 1.871 0.061452 .
## `Entidad federativa`23 3.531e+03 2.299e+03 1.536 0.124735
## `Entidad federativa`24 1.038e+03 1.737e+03 0.598 0.550225
## `Entidad federativa`25 1.128e+03 2.011e+03 0.561 0.575052
## `Entidad federativa`26 -1.358e+03 1.704e+03 -0.797 0.425355
## `Entidad federativa`27 6.083e+03 2.039e+03 2.984 0.002875 **
## `Entidad federativa`28 1.362e+02 1.784e+03 0.076 0.939123
## `Entidad federativa`29 -2.701e+03 1.738e+03 -1.554 0.120371
## `Entidad federativa`30 -1.256e+03 1.640e+03 -0.766 0.443845
## `Entidad federativa`31 -1.850e+03 1.679e+03 -1.102 0.270576
## `Entidad federativa`32 -1.187e+03 1.741e+03 -0.682 0.495544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5242 on 2416 degrees of freedom
## Multiple R-squared: 0.5935, Adjusted R-squared: 0.5869
## F-statistic: 90.45 on 39 and 2416 DF, p-value: < 2.2e-16VIF(non_spatial_model)## GVIF Df GVIF^(1/(2*Df))
## hospitalespHabitantes 1.281813 1 1.132172
## popden2020 1.640840 1 1.280953
## crimen_2018 1.301022 1 1.140624
## crimen_2019 1.310719 1 1.144866
## porcentaje_pob_pobreza 8.280886 1 2.877653
## porcentaje_pob_pobreza_ext 7.141961 1 2.672445
## porcentaje_pob_acceso_ss 2.670169 1 1.634065
## rezago_social 5.491570 1 2.343410
## `Entidad federativa` 9.247524 31 1.036528# detecting spatially autocorrelated regression residuals
moran.test(non_spatial_model$residuals, rswm_queen) ##
## Moran I test under randomisation
##
## data: non_spatial_model$residuals
## weights: rswm_queen
##
## Moran I statistic standard deviate = -1.66, p-value = 0.9515
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.0202115681 -0.0004073320 0.0001423322# convert to sp
municipios.sp = as(map_deprctd, "Spatial")
mx_geodata <-geo_join(map_deprctd,merged_data,'IDUNICO','IDUNICO',how='inner') ### Combining geospatial and non-spatial data## Warning: We recommend using the dplyr::*_join() family of functions instead.El primer modelo tuvo el uso adicional de identificar cuales variables no son significativas, en este caso, son el índice de crímen, pobreza (normal), y la mayoría de las entidades federativas.
Con la función VIF, tenemos un insight adicional a Multicollinearity, y así sabemos qué variables vale la pena quitar para que nuestro modelo sea más viable, en este caso, también deberemos de quitar pobreza y Rezago Social
bw1 <- bw.gwr(infectados ~ log(hospitalespHabitantes) + log(popden2020) + porcentaje_pob_acceso_ss, approach = "AIC", adaptive = T, data = mx_geodata)## Take a cup of tea and have a break, it will take a few minutes.
## -----A kind suggestion from GWmodel development group
## Adaptive bandwidth (number of nearest neighbours): 1518 AICc value: 50203.16
## Adaptive bandwidth (number of nearest neighbours): 946 AICc value: 50089.14
## Adaptive bandwidth (number of nearest neighbours): 591 AICc value: 49990.68
## Adaptive bandwidth (number of nearest neighbours): 373 AICc value: 49860.7
## Adaptive bandwidth (number of nearest neighbours): 237 AICc value: 49714.65
## Adaptive bandwidth (number of nearest neighbours): 154 AICc value: 49641.87
## Adaptive bandwidth (number of nearest neighbours): 101 AICc value: 49478.85
## Adaptive bandwidth (number of nearest neighbours): 70 AICc value: 49343.97
## Adaptive bandwidth (number of nearest neighbours): 49 AICc value: 49320.68
## Adaptive bandwidth (number of nearest neighbours): 38 AICc value: 49360.93
## Adaptive bandwidth (number of nearest neighbours): 58 AICc value: 49303.55
## Adaptive bandwidth (number of nearest neighbours): 61 AICc value: 49317.76
## Adaptive bandwidth (number of nearest neighbours): 53 AICc value: 49317.84
## Adaptive bandwidth (number of nearest neighbours): 58 AICc value: 49303.55# determine the kernel bandwidth
bw2 <- bw.gwr(infectados ~ log(hospitalespHabitantes) + log(popden2020) + porcentaje_pob_acceso_ss, approach = "AIC", adaptive = F, data = mx_geodata)## Take a cup of tea and have a break, it will take a few minutes.
## -----A kind suggestion from GWmodel development group
## Fixed bandwidth: 19.92188 AICc value: 50330.11
## Fixed bandwidth: 12.31486 AICc value: 50306.35
## Fixed bandwidth: 7.613466 AICc value: 50273.75
## Fixed bandwidth: 4.707843 AICc value: 50180.22
## Fixed bandwidth: 2.912069 AICc value: 50081.11
## Error in eval(expr, envir, enclos) : inv(): matrix is singular
## Fixed bandwidth: 1.80222 AICc value: Inf
## Fixed bandwidth: 3.597994 AICc value: 50117.54
## Fixed bandwidth: 2.488145 AICc value: 50046.87
## Fixed bandwidth: 2.226145 AICc value: 50078.21
## Fixed bandwidth: 2.650069 AICc value: 50062.75
## Fixed bandwidth: 2.38807 AICc value: 50342.75
## Fixed bandwidth: 2.549994 AICc value: 50055.18
## Fixed bandwidth: 2.449919 AICc value: 50400.61
## Fixed bandwidth: 2.511769 AICc value: 50051.65
## Fixed bandwidth: 2.473544 AICc value: 50048.32
## Fixed bandwidth: 2.497168 AICc value: 50062.48
## Fixed bandwidth: 2.482568 AICc value: 50049.65
## Fixed bandwidth: 2.491591 AICc value: 50051.51
## Fixed bandwidth: 2.486014 AICc value: 50048.74
## Fixed bandwidth: 2.489461 AICc value: 50051.29
## Fixed bandwidth: 2.487331 AICc value: 50040.96
## Fixed bandwidth: 2.486828 AICc value: 50049.37
## Fixed bandwidth: 2.487642 AICc value: 50049.85
## Fixed bandwidth: 2.487139 AICc value: 50050.59
## Fixed bandwidth: 2.48745 AICc value: 50049.94
## Fixed bandwidth: 2.487258 AICc value: 50050.84
## Fixed bandwidth: 2.487376 AICc value: 50048.54
## Fixed bandwidth: 2.487303 AICc value: 50049.33# fit the GWR model
m.gwr <- gwr.basic(infectados ~ log(hospitalespHabitantes) + log(popden2020) + porcentaje_pob_acceso_ss, adaptive = T, data = mx_geodata, bw = bw1)
m.gwr## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2023-06-13 22:18:15
## Call:
## gwr.basic(formula = infectados ~ log(hospitalespHabitantes) +
## log(popden2020) + porcentaje_pob_acceso_ss, data = mx_geodata,
## bw = bw1, adaptive = T)
##
## Dependent (y) variable: infectados
## Independent variables: hospitalespHabitantes popden2020 porcentaje_pob_acceso_ss
## Number of data points: 2445
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16081 -2249 -456 1145 144379
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5170.06 1046.08 4.942 8.24e-07 ***
## log(hospitalespHabitantes) 400.64 189.30 2.116 0.0344 *
## log(popden2020) 1414.84 84.16 16.811 < 2e-16 ***
## porcentaje_pob_acceso_ss -13815.62 1022.97 -13.505 < 2e-16 ***
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 7164 on 2441 degrees of freedom
## Multiple R-squared: 0.1933
## Adjusted R-squared: 0.1923
## F-statistic: 195 on 3 and 2441 DF, p-value: < 2.2e-16
## ***Extra Diagnostic information
## Residual sum of squares: 125279733715
## Sigma(hat): 7161.083
## AIC: 50352.28
## AICc: 50352.31
## BIC: 47975.3
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: bisquare
## Adaptive bandwidth: 58 (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.
## Intercept -124474.121 -2703.139 -143.984 1215.337
## log(hospitalespHabitantes) -2781.437 40.173 203.678 761.674
## log(popden2020) -399.728 171.539 554.714 1430.157
## porcentaje_pob_acceso_ss -98024.336 -6367.964 -1970.714 -172.436
## Max.
## Intercept 52747
## log(hospitalespHabitantes) 11933
## log(popden2020) 21104
## porcentaje_pob_acceso_ss 25460
## ************************Diagnostic information*************************
## Number of data points: 2445
## Effective number of parameters (2trace(S) - trace(S'S)): 545.9977
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 1899.002
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 49303.55
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 48719.54
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 49082.74
## Residual sum of squares: 54488498271
## R-square value: 0.6491546
## Adjusted R-square value: 0.5482271
##
## ***********************************************************************
## Program stops at: 2023-06-13 22:18:17# Mapping GWR outputs
gwr_sf = st_as_sf(m.gwr$SDF)
# Note: _TV stands for t-values. t-values less than -1.96 and greater than 1.96 are statistically significant.
# local prediction of dependent variable
gwr_sf$y_predicted <- log(gwr_sf$yhat)
coordsDF <- data.frame(x = coords[,1], y = coords[,2])
#mapview(gwr_sf, zcol="y_predicted")
# local prediction of statistically significant explanatory variables
temp <- as.data.frame(gwr_sf$pop_density_TV)
ggplot() +
geom_sf(data = gwr_sf, aes(fill = y_predicted))
ggplot() +
geom_sf(data = gwr_sf, aes(fill = porcentaje_pob_acceso_ss_SE))
# local prediction of R2
ggplot() +
geom_sf(data = gwr_sf, aes(fill = Local_R2))
# local regression residuals
ggplot() +
geom_sf(data = gwr_sf, aes(fill = residual))
Por último y como esperábamos, y_predicted tiene un patrón de colores más brillantes en zonas con población algo densa, pero también rural, es decir, sin mucho acceso a los servicios de salud, de igual manera, el mapa de los residuales no tiene patrón alguno, por lo que podemos declarar que nuestro modelo no tiene problemas de sesgo.
Observando los resultados obtenidos en el modelo GWR, podemos observar como las variables de Hospitales por Habitantes, Densidad Poblacional y Porcentaje de Personas con Acceso a Servicio Social, cuentan con significancia estadística, sin embargo, el AIC mostrado es muy alto, por lo que a pesar de que este modelo explica mejor el comportamiento de las variables, es peor para predecir la variable dependiente.