Actividad 2 Modulo 1, Equipo 4
Situación Problema
De acuerdo a la Cámara Nacional de la Industria Farmacéutica, (CANIFARMA) en México las personas en situación de pobreza se caracterizan por tener una probabilidad 5 veces mayor de fallecer por COVID-19 que las personas con relativamente mayor nivel de ingresos (Arceo-Gómez, et al., 2021)2. Además de la falta de acceso a servicios de salud y posibles cormobilidades, otro factor relevante en incrementar dicha probabilidad es el perfil socioeconómico (Arceo-Gómez, et al., 2021). A partir de la pandemia por COVI19, la firma de consultoría XYZ (México) establece que “Las organizaciones que su principal actividad de negocios es brindar servicios de salud requieren soluciones específicas e innovadoras, para aprovechar oportunidades, afrontar retos, así como favorecer su consolidación y crecimiento”. Algunos de los servicios enfocados por parte de la firma es detectar las necesidades y potencial del crecimiento del sector salud a partir de Analítica de Datos.
¿Cómo se relaciona ESDA con el Location Intelligence
Location Intelligence es el uso de datos geo-espaciales para el desarrollo de información que ayude a la toma de decisiones en un ambiente de negocios. ESDA, es una forma de transformar datos geo-espaciales a información, por lo que ESDA, es una forma de hacer Location Intelligence y debe ser visto como una herramienta para eso, más no la totalidad del proceso.
Carga de las bases de datos
Uniendo las bases de datos.
#Preparing the variable 'clee' to match 'cve_ent' to do the joint
hospitales$cve_ent <- substr(hospitales$clee, 1, 5)
covid$cve_ent <- as.character(covid$cve_ent)
covid$cve_ent <- ifelse(str_length(covid$cve_ent) == 4, paste("0",covid$cve_ent, sep = ""),covid$cve_ent)#Deleting all the health centers that aren't considered as hospitals or related to the treatment of Covid or related.
hospitalesFiltered <- filter(hospitales, codigo_act != 624191 & codigo_act !=
623311 & codigo_act !=
621331 & codigo_act !=
624198 & codigo_act !=
621311 & codigo_act !=
623991 & codigo_act !=
623992 & codigo_act !=
624221 & codigo_act !=
624112 & codigo_act !=
624111 & codigo_act !=
624311 & codigo_act !=
622311 & codigo_act !=
624411 & codigo_act !=
624412 & codigo_act !=
624122 & codigo_act !=
624222 & codigo_act !=
624312 & codigo_act !=
621411 & codigo_act !=
623312 & codigo_act !=
624121 & codigo_act !=
621412 & codigo_act !=
621312)#Creating a Dataset that counts the number of hospitals in México
hospitalesMexico <- hospitalesFiltered %>% count(cve_ent,entidad)
hospitalesMexico <- hospitalesMexico[-c(2),]#Adding up the cases of Covid 19 into 2 columns
covid$infectados2020 <- rowSums(covid[,c("feb_2020", "march_2020","april_2020","may_2020","june_2020","july_2020","august_2020","sept_2020","oct_2020","nov_2020","dic_2020")])
covid$infectados2021 <- rowSums(covid[,c("jan_2021", "feb_2021", "mar_2021","april_2021","may_2021","june_2021","july_2021","august_2021","sept_2021","oct_2021","nov_2021","dic_2021")])densidad <- readxl::read_excel("C:\\Users\\art191127\\Downloads\\Población_07.xlsx")
superficie <- readxl::read_excel("C:\\Users\\art191127\\Downloads\\inegi_superficie_territorial.xlsx",sheet = "Distribución territorial")Creando la variable Hospitales x Habitantes
superficie <- na.omit(superficie)
colnames(superficie) <- c("Entidad federativa", "Municipio", "Superficie (Km2)", "Porcentaje de la superficie estatal","Demsidad de población (hab./km2)","Total de municipios o demarcaciones territoriales / localidades")
superficie$`Entidad federativa` <- substr(superficie$`Entidad federativa`,1,2)
superficie$Municipio <- substr(superficie$Municipio,1,3)
superficie$cve_ent <- paste(superficie$`Entidad federativa`,superficie$Municipio,sep = "")
#Adding population density into hospitalesMexico dataset
hospitalesMexico <- merge(hospitalesMexico,superficie,by.x = "cve_ent", by.y = "cve_ent")
hospitalesMexico <- merge(hospitalesMexico, covid %>% select(cve_ent, poblacion_2022), by.x = "cve_ent", by.y = "cve_ent")
#Deleting extra municipalities
hospitalesMexico <- hospitalesMexico[-c(2457, 2458), ]
#Creating the columns for KM2 extension and population density
hospitalesMexico$hospitalesKm2 <- hospitalesMexico$n/hospitalesMexico$`Superficie (Km2)`
hospitalesMexico$hospitalespHabitantes <- (hospitalesMexico$n/hospitalesMexico$poblacion_2022)*10000map <- st_read("C:\\Users\\art191127\\Documents\\Tec de Monterrey\\Semestre 8\\Clase 2\\Saucedo\\shp_mx_mpios\\mx_mpios.shp")## Reading layer `mx_mpios' from data source
## `C:\Users\art191127\Documents\Tec de Monterrey\Semestre 8\Clase 2\Saucedo\shp_mx_mpios\mx_mpios.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2456 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -118.4076 ymin: 14.5321 xmax: -86.71041 ymax: 32.71865
## Geodetic CRS: WGS 84map_deprctd <- readShapePoly("C:\\Users\\art191127\\Documents\\Tec de Monterrey\\Semestre 8\\Clase 2\\Saucedo\\shp_mx_mpios\\mx_mpios.shp", IDvar = "IDUNICO", proj4string = CRS("+proj=longlat"), sf_use_s2(FALSE))## Shapefile type: Polygon, (5), # of Shapes: 2456Ajustando variables
# Transforming selected variables to plot them
covid <- covid %>% mutate_at(c("hogrem2015", "hogremjefmuj2015", "popnoafmed2015", "gini2015", "porcentaje_pob_pobreza", "porcentaje_pob_pobreza_ext", "porcentaje_pob_servicios_salud", "porcentaje_pob_acceso_ss", "popden2020"), as.numeric)## Warning: There were 9 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `hogrem2015 = .Primitive("as.double")(hogrem2015)`.
## Caused by warning:
## ! NAs introducidos por coerción
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 8 remaining warnings.#covid1 <- na.omit(covid1)
# Dividing the percentage variables between 100 for better interpretation
covid$porcentaje_pob_pobreza <- covid$porcentaje_pob_pobreza/100
covid$porcentaje_pob_acceso_ss <- covid$porcentaje_pob_acceso_ss/100
covid$porcentaje_pob_pobreza_ext <- covid$porcentaje_pob_pobreza_ext/100
covid$porcentaje_pob_servicios_salud <- covid$porcentaje_pob_servicios_salud/100
covid$cve_ent <- as.character(covid$cve_ent)Modelo de pre-elección
lm_model <- lm(log1p(infectados2020) ~ porcentaje_pob_pobreza_ext + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + log(popden2020) + hogrem2015 + hogremjefmuj2015 + gini2015 + inclusion_fin_2019 + crimen_2019 + grado_rs + popnoafmed2015 + porcentaje_pob_acceso_ss, data = covid)
summary(lm_model)##
## Call:
## lm(formula = log1p(infectados2020) ~ porcentaje_pob_pobreza_ext +
## porcentaje_pob_pobreza + porcentaje_pob_servicios_salud +
## rezago_social + log(popden2020) + hogrem2015 + hogremjefmuj2015 +
## gini2015 + inclusion_fin_2019 + crimen_2019 + grado_rs +
## popnoafmed2015 + porcentaje_pob_acceso_ss, data = covid)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.8128 -0.8423 -0.0292 0.7749 3.8881
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0148033 0.3685754 -0.040 0.96797
## porcentaje_pob_pobreza_ext -0.1040968 0.4547244 -0.229 0.81895
## porcentaje_pob_pobreza -2.8767105 0.2922149 -9.845 < 2e-16 ***
## porcentaje_pob_servicios_salud 0.8247289 0.2828436 2.916 0.00358 **
## rezago_social -0.2718712 0.0989677 -2.747 0.00606 **
## log(popden2020) 0.4984695 0.0180212 27.660 < 2e-16 ***
## hogrem2015 -0.0141398 0.0033839 -4.179 3.04e-05 ***
## hogremjefmuj2015 -0.0014831 0.0044526 -0.333 0.73909
## gini2015 9.7601555 0.7310990 13.350 < 2e-16 ***
## inclusion_fin_2019 0.6264764 0.0362137 17.299 < 2e-16 ***
## crimen_2019 0.0047722 0.0008101 5.891 4.37e-09 ***
## grado_rsBajo 0.4555337 0.1745176 2.610 0.00910 **
## grado_rsMedio 0.3433603 0.1264246 2.716 0.00666 **
## grado_rsMuy alto 0.1685609 0.1699354 0.992 0.32134
## grado_rsMuy bajo 0.4469024 0.2253911 1.983 0.04750 *
## popnoafmed2015 -0.0205013 0.0044851 -4.571 5.10e-06 ***
## porcentaje_pob_acceso_ss -1.1300702 0.2713014 -4.165 3.22e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.227 on 2438 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.6959, Adjusted R-squared: 0.6939
## F-statistic: 348.7 on 16 and 2438 DF, p-value: < 2.2e-16Antes de pasar a los modelos que se van a comparar, hacemos un modelo rápido con todas las variables y determinamos el comportamiento de algunas, así como su relevancia, para solo quedarnos con las relevantes.
Mergin the data frames
colnames(map)[3] ="cve_ent"
map$cve_ent <- as.character(map$cve_ent)
# Deleting the extra row
map <- map[-2456,]
map$cve_ent <- as.character(map$cve_ent)
# Matching the key of the dataset 'map' with the key of 'covid'
map$cve_ent <- ifelse(str_length(map$cve_ent) == 4, paste("0",map$cve_ent, sep = ""),map$cve_ent)
connectivityM <- merge(map, covid, by = "cve_ent")
# Selecting only the interest variables
vector <- as.vector(c("cve_ent","gini2015", "porcentaje_pob_pobreza", "porcentaje_pob_servicios_salud", "rezago_social", "inclusion_fin_2019", "porcentaje_pob_acceso_ss", "popden2020","infectados2020", "popnoafmed2015", "crimen_2019", "geometry"))
connectivityM <- connectivityM[,vector] hospitalesMexico$cve_ent <- as.character(hospitalesMexico$cve_ent)
hospitalesMexico$cve_ent <- ifelse(str_length(hospitalesMexico$cve_ent) == 4, paste("0",hospitalesMexico$cve_ent, sep = ""),hospitalesMexico$cve_ent)
connectivityM <- merge(connectivityM, hospitalesMexico, by = "cve_ent")
vector <- as.vector(c("gini2015", "porcentaje_pob_pobreza", "porcentaje_pob_servicios_salud", "rezago_social", "inclusion_fin_2019", "porcentaje_pob_acceso_ss", "popden2020","hospitalespHabitantes", "infectados2020", "popnoafmed2015", "crimen_2019", "geometry"))
connectivityM <- connectivityM[,vector]
#map_deprctd <- map_deprctd[-2457,]
connectivityM <- na.omit(connectivityM)Explorando la variable dependiente
Antes de hacer modelos, es necesario identificar el comportamiento de la variable de infectados, pues al igual que con la variable de población, es muy posible que cuente con outliers muy fuertes. Por lo tanto hacemos un boxplot.
boxplot(connectivityM$infectados2020)
Como se esperaba, la variable tiene un comportamiento problemático, por lo que para obtener una distribución más normal, haremos una transformación log1p, que es una versión de log que maneja mucho mejor los valores 0 o cercanos a este.
boxplot(log1p(connectivityM$infectados2020))
Después de la transformación, la variable se ve mucho mejor, y aún cuenta con municipios que resaltan por sus contagios, más ya no afectan como outliers.
connectivityM$popden2020 <- log(connectivityM$popden2020)
connectivityM$infectados2020 <- log1p(connectivityM$infectados2020)Creating train and test set
set.seed(257)
partition <- createDataPartition(y = connectivityM$infectados2020, p=0.7, list=F)
train <- connectivityM[partition, ]
test <- connectivityM[-partition, ]Modelo: Regresión lineal simple
lm_model <- lm(infectados2020 ~ gini2015 + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + popden2020 + inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes + popnoafmed2015 + crimen_2019, data = train)
summary(lm_model)##
## Call:
## lm(formula = infectados2020 ~ gini2015 + porcentaje_pob_pobreza +
## porcentaje_pob_servicios_salud + rezago_social + popden2020 +
## inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes +
## popnoafmed2015 + crimen_2019, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6922 -0.7837 -0.0520 0.7290 4.4283
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.797251 0.382134 2.086 0.0371 *
## gini2015 9.345167 0.841825 11.101 < 2e-16 ***
## porcentaje_pob_pobreza -2.181556 0.266597 -8.183 5.35e-16 ***
## porcentaje_pob_servicios_salud 1.379920 0.317367 4.348 1.45e-05 ***
## rezago_social -0.482311 0.049751 -9.694 < 2e-16 ***
## popden2020 0.476322 0.019395 24.559 < 2e-16 ***
## inclusion_fin_2019 0.663356 0.041236 16.087 < 2e-16 ***
## porcentaje_pob_acceso_ss -1.995072 0.295892 -6.743 2.12e-11 ***
## hospitalespHabitantes -0.029365 0.003140 -9.353 < 2e-16 ***
## popnoafmed2015 -0.020816 0.005110 -4.074 4.84e-05 ***
## crimen_2019 0.004732 0.000894 5.293 1.36e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.193 on 1708 degrees of freedom
## Multiple R-squared: 0.7177, Adjusted R-squared: 0.7161
## F-statistic: 434.3 on 10 and 1708 DF, p-value: < 2.2e-16Dado el cuidado que se tuvo con la selección de variables, podemos ver como todas parecen ser significativas dentro del modelo.
Gráfica: Efecto de la inclusión en infectados
plot(effect("inclusion_fin_2019",lm_model))
Ahora utilizamos la variable inclusión, que esperábamos no fuera relevante para el modelo, en está gráfica para ver su relación con la variable dependiente. Como se puede ver se sigue una tendencia positiva.
Gráfica: Predicted vs Observed Infectados
ggplot(train, aes(x = exp(lm_model$fitted.values), y = train$infectados2020)) +
geom_point() +
stat_smooth() +
labs(x='Predicted Values', y='Actual Values', title='OLS Predicted vs. Actual Values')## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Al visualizar está gráfica, nos podemos dar cuenta que nuestro modelo está lejos de ser perfecto, pues al predecir, obtuvo un valor muy grande en comparación al real, creando un outlier por su propia cuenta.
Multicollinearity test
VIF(lm_model)## gini2015 porcentaje_pob_pobreza
## 1.058609 4.175998
## porcentaje_pob_servicios_salud rezago_social
## 1.869221 3.083361
## popden2020 inclusion_fin_2019
## 1.416338 1.538518
## porcentaje_pob_acceso_ss hospitalespHabitantes
## 2.391673 1.179617
## popnoafmed2015 crimen_2019
## 1.695696 1.025160Ahora probamos el modelo para ver si cuenta con multicollinearity, pero todas las variables se encuentran entre los valores 1 y 5, por lo que podemos decir que no hay.
Homoscedasticity Test
bptest(lm_model)##
## studentized Breusch-Pagan test
##
## data: lm_model
## BP = 131.48, df = 10, p-value < 2.2e-16El p-value es mejor a 0.05, y por lo tanto no hay homoscedasticity, pero si hay heteroscedasticity.
Gráfica: Residuals vs Fit Plot
plot(fitted(lm_model), resid(lm_model), main="Linear Regression Residual vs. Fitted Values", xlab="Fitted Values", ylab="Residuals")
abline(0,0)
Los residuales no se comportan de manera horizontal a la línea de la regresión estima, por lo tanto no hay Homoscedasticity.
Gráfica: Normality of residuals
hist(lm_model$residuals, xlab="Estimated Regression Residuals", main='Distribution of OLS Estimated Regression Residuals', col='lightblue', border="white")
Los residuales tienen una distribución normal como era de esperarse para una regresión lineal, acercando así los errores al 0.
RMSE: Modelo Regresión Lineal.
prediction_lm_model <- lm_model %>% predict(test)
RMSE(prediction_lm_model, test$infectados2020)## [1] 1.24692Modelo: WLS
Para evitar la heteroscedasticity creamos un modelo Weighted Least Squares
wt <- 1 / lm(abs(lm_model$residuals) ~ lm_model$fitted.values)$fitted.values^2
wls_model<-lm(infectados2020 ~ gini2015 + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + popden2020 + inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes + popnoafmed2015 + crimen_2019, data = train, weights=wt)
summary(wls_model)##
## Call:
## lm(formula = infectados2020 ~ gini2015 + porcentaje_pob_pobreza +
## porcentaje_pob_servicios_salud + rezago_social + popden2020 +
## inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes +
## popnoafmed2015 + crimen_2019, data = train, weights = wt)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -3.7726 -0.8339 -0.0486 0.7878 5.6383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6287739 0.3804088 1.653 0.0985 .
## gini2015 9.5126148 0.8330972 11.418 < 2e-16 ***
## porcentaje_pob_pobreza -2.4620156 0.2557329 -9.627 < 2e-16 ***
## porcentaje_pob_servicios_salud 1.3536627 0.3052202 4.435 9.79e-06 ***
## rezago_social -0.4232849 0.0456177 -9.279 < 2e-16 ***
## popden2020 0.4637735 0.0196288 23.627 < 2e-16 ***
## inclusion_fin_2019 0.7404021 0.0453906 16.312 < 2e-16 ***
## porcentaje_pob_acceso_ss -1.6165166 0.2919605 -5.537 3.56e-08 ***
## hospitalespHabitantes -0.0269722 0.0027475 -9.817 < 2e-16 ***
## popnoafmed2015 -0.0206731 0.0048529 -4.260 2.16e-05 ***
## crimen_2019 0.0042365 0.0008773 4.829 1.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.268 on 1708 degrees of freedom
## Multiple R-squared: 0.7046, Adjusted R-squared: 0.7028
## F-statistic: 407.3 on 10 and 1708 DF, p-value: < 2.2e-16Todas las variables son significativas.
bptest(wls_model)##
## studentized Breusch-Pagan test
##
## data: wls_model
## BP = 159.9, df = 10, p-value < 2.2e-16Gráfica: Residual vs Fitted
plot(fitted(wls_model), resid(wls_model), main="Residual vs. Fitted Values", xlab="Fitted Values", ylab="Residuals")
abline(0,0)
Los residuales no se comportan de manera horizontal a la línea de la regresión estima, por lo tanto no hay Homoscedasticity.
RMSE: WLS model
prediction_wls_model <- wls_model %>% predict(test)
RMSE(prediction_wls_model, test$infectados2020) ## [1] 1.244346Modelo: SVM
fit.svm <- svm(infectados2020 ~ gini2015 + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + popden2020 + inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes + popnoafmed2015 + crimen_2019, data = train, type = 'eps-regression', kernel = 'radial')
summary(fit.svm)##
## Call:
## svm(formula = infectados2020 ~ gini2015 + porcentaje_pob_pobreza +
## porcentaje_pob_servicios_salud + rezago_social + popden2020 +
## inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes +
## popnoafmed2015 + crimen_2019, data = train, type = "eps-regression",
## kernel = "radial")
##
##
## Parameters:
## SVM-Type: eps-regression
## SVM-Kernel: radial
## cost: 1
## gamma: 0.1
## epsilon: 0.1
##
##
## Number of Support Vectors: 1418Gráfica: Residual vs Fitted
plot(fit.svm$fitted, fit.svm$residuals, main="SVM Residual vs. Fitted Values", xlab="Fitted Values", ylab="Residuals")
abline(0,0)
Los residuales se comportan de manera horizontal respecto a la línea estimada de regresión, por lo que se intuye que hay homoscedasticity.
RMSE: SVM model
predicted_dv=predict(fit.svm, newdata = test)
rmse(predicted_dv, test$infectados2020)## [1] 1.029832Gráfica: Predicted vs Actual Values
dv_svm<-data.frame(exp(fit.svm$fitted),train$infectados2020)
ggplot(dv_svm, aes(x = exp.fit.svm.fitted., y = train.infectados2020)) +
geom_point() +
stat_smooth() +
labs(x='Predicted Values', y='Actual Values', title='SVM Predicted vs. Actual Values')## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
A pesar de que mejoro el RMSE considerablemente, aún hay problemas con la generación de un outlier.
RANDOM FOREST
random_forest<-randomForest(infectados2020 ~ gini2015 + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + popden2020 + inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes + popnoafmed2015 + crimen_2019, data = train, proximity=TRUE)
print(random_forest)##
## Call:
## randomForest(formula = infectados2020 ~ gini2015 + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + popden2020 + inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes + popnoafmed2015 + crimen_2019, data = train, proximity = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 3
##
## Mean of squared residuals: 0.9460134
## % Var explained: 81.12RMSE: Random Forest Model
rf_prediction_test_data <-predict(random_forest,test)
rmse(rf_prediction_test_data, test$infectados2020)## [1] 0.9953379Gráfica: Variable Importance
varImpPlot(random_forest, n.var = 10, main = "Top 10 - Variable")
En este modelo, la variable inclusión es la más importante, mientras que la Población sin Acceso a Servicios de Salud es la que menos.
importance(random_forest) ## IncNodePurity
## gini2015 268.3932
## porcentaje_pob_pobreza 1191.1899
## porcentaje_pob_servicios_salud 281.6057
## rezago_social 1500.8790
## popden2020 877.4876
## inclusion_fin_2019 2337.9893
## porcentaje_pob_acceso_ss 707.2152
## hospitalespHabitantes 301.7978
## popnoafmed2015 199.2086
## crimen_2019 758.0459En este resultado, se puede ver que tanto el error incrementa si se elimina cada variable.
SPATIAL REGRESSION ANALYSIS
Moran Test Linear Model
Usaremos un modelo lineal simple para hacer el Moran Test.
lm_model_alt <- lm(infectados2020 ~ gini2015 + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + popden2020 + inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes + popnoafmed2015 + crimen_2019, data = connectivityM)
summary(lm_model_alt)##
## Call:
## lm(formula = infectados2020 ~ gini2015 + porcentaje_pob_pobreza +
## porcentaje_pob_servicios_salud + rezago_social + popden2020 +
## inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes +
## popnoafmed2015 + crimen_2019, data = connectivityM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.1716 -0.7845 -0.0351 0.7458 4.3301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5942866 0.3214904 1.849 0.0646 .
## gini2015 9.9165119 0.7068047 14.030 < 2e-16 ***
## porcentaje_pob_pobreza -2.3974276 0.2259255 -10.612 < 2e-16 ***
## porcentaje_pob_servicios_salud 1.0949504 0.2682924 4.081 4.62e-05 ***
## rezago_social -0.4794235 0.0427411 -11.217 < 2e-16 ***
## popden2020 0.4716068 0.0166314 28.356 < 2e-16 ***
## inclusion_fin_2019 0.6812239 0.0352210 19.341 < 2e-16 ***
## porcentaje_pob_acceso_ss -1.7199649 0.2511186 -6.849 9.36e-12 ***
## hospitalespHabitantes -0.0312555 0.0027317 -11.442 < 2e-16 ***
## popnoafmed2015 -0.0180319 0.0043498 -4.145 3.51e-05 ***
## crimen_2019 0.0044926 0.0007897 5.689 1.43e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.208 on 2442 degrees of freedom
## Multiple R-squared: 0.7107, Adjusted R-squared: 0.7095
## F-statistic: 599.8 on 10 and 2442 DF, p-value: < 2.2e-16connectivityM.link_a_rook<-poly2nb(connectivityM,queen=T)
connectivityM.linkW_a_rook<-nb2listw(connectivityM.link_a_rook, style="W")lm.morantest(lm_model_alt, nb2listw(connectivityM.link_a_rook))##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = infectados2020 ~ gini2015 + porcentaje_pob_pobreza
## + porcentaje_pob_servicios_salud + rezago_social + popden2020 +
## inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes +
## popnoafmed2015 + crimen_2019, data = connectivityM)
## weights: nb2listw(connectivityM.link_a_rook)
##
## Moran I statistic standard deviate = 33.142, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.4001735307 -0.0020647314 0.0001473046Al tener un p-value significativo, y un Moran positivo de 0.4001, se puede confirmar la existencia de autocorrelación espacial positiva, por lo tanto si se pueden y deben hacer modelos autoregresivos.
SPATIAL AUTOGRESSIVE MODEL
lmat_c<-coordinates(map_deprctd)
map.centroid_c<-coordinates(map_deprctd)
connectivityM.link_a_rook<-poly2nb(connectivityM,queen=T)## although coordinates are longitude/latitude, st_intersects assumes that they
## are planarconnectivityM.linkW_a_rook<-nb2listw(connectivityM.link_a_rook, style="W")spatial_lag_model <- lagsarlm(infectados2020 ~ gini2015 + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + popden2020 + inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes + popnoafmed2015 + crimen_2019, data = connectivityM, nb2listw(connectivityM.link_a_rook), method="Matrix")
summary(spatial_lag_model)##
## Call:lagsarlm(formula = infectados2020 ~ gini2015 + porcentaje_pob_pobreza +
## porcentaje_pob_servicios_salud + rezago_social + popden2020 +
## inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes +
## popnoafmed2015 + crimen_2019, data = connectivityM, listw = nb2listw(connectivityM.link_a_rook),
## method = "Matrix")
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.075238 -0.687010 -0.012965 0.640372 4.425443
##
## Type: lag
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.18197954 0.29360979 -4.0257 5.681e-05
## gini2015 9.51085665 0.62309173 15.2640 < 2.2e-16
## porcentaje_pob_pobreza -1.01701970 0.20722786 -4.9077 9.213e-07
## porcentaje_pob_servicios_salud 1.03149807 0.23870720 4.3212 1.552e-05
## rezago_social -0.30853176 0.03831501 -8.0525 8.882e-16
## popden2020 0.33411032 0.01558350 21.4400 < 2.2e-16
## inclusion_fin_2019 0.62510035 0.03100494 20.1613 < 2.2e-16
## porcentaje_pob_acceso_ss -1.92462703 0.22185106 -8.6753 < 2.2e-16
## hospitalespHabitantes -0.01326906 0.00250801 -5.2907 1.219e-07
## popnoafmed2015 -0.01759362 0.00385365 -4.5654 4.984e-06
## crimen_2019 0.00347585 0.00069528 4.9992 5.756e-07
##
## Rho: 0.41206, LR test value: 552.02, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.016118
## z-value: 25.565, p-value: < 2.22e-16
## Wald statistic: 653.58, p-value: < 2.22e-16
##
## Log likelihood: -3663.047 for lag model
## ML residual variance (sigma squared): 1.1235, (sigma: 1.0599)
## Number of observations: 2453
## Number of parameters estimated: 13
## AIC: 7352.1, (AIC for lm: 7902.1)Todas las variables salieron significativas.
RMSE: SLM
sqrt(mean((connectivityM$infectados2020 - spatial_lag_model$fitted.values)^2))## [1] 1.059941SPATIAL ERROR MODEL
spatial_error_model <- errorsarlm(infectados2020 ~ gini2015 + porcentaje_pob_pobreza + porcentaje_pob_servicios_salud + rezago_social + popden2020 + inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes + popnoafmed2015 + crimen_2019, data = connectivityM, nb2listw(connectivityM.link_a_rook), method="Matrix")
summary(spatial_error_model)##
## Call:errorsarlm(formula = infectados2020 ~ gini2015 + porcentaje_pob_pobreza +
## porcentaje_pob_servicios_salud + rezago_social + popden2020 +
## inclusion_fin_2019 + porcentaje_pob_acceso_ss + hospitalespHabitantes +
## popnoafmed2015 + crimen_2019, data = connectivityM, listw = nb2listw(connectivityM.link_a_rook),
## method = "Matrix")
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8789254 -0.6122848 -0.0059539 0.5947503 3.6752538
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.41814465 0.32148022 -1.3007 0.193366
## gini2015 8.85576027 0.63796512 13.8813 < 2.2e-16
## porcentaje_pob_pobreza -1.66641960 0.25127804 -6.6318 3.317e-11
## porcentaje_pob_servicios_salud 1.37409168 0.23406132 5.8706 4.341e-09
## rezago_social -0.33485185 0.04590027 -7.2952 2.982e-13
## popden2020 0.59732790 0.02205653 27.0817 < 2.2e-16
## inclusion_fin_2019 0.48287548 0.02980776 16.1997 < 2.2e-16
## porcentaje_pob_acceso_ss -1.90591591 0.25646738 -7.4314 1.075e-13
## hospitalespHabitantes -0.00784171 0.00251769 -3.1146 0.001842
## popnoafmed2015 -0.01283189 0.00401535 -3.1957 0.001395
## crimen_2019 0.00202139 0.00069007 2.9293 0.003398
##
## Lambda: 0.69519, LR test value: 856.54, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.018216
## z-value: 38.163, p-value: < 2.22e-16
## Wald statistic: 1456.4, p-value: < 2.22e-16
##
## Log likelihood: -3510.786 for error model
## ML residual variance (sigma squared): 0.9202, (sigma: 0.95927)
## Number of observations: 2453
## Number of parameters estimated: 13
## AIC: 7047.6, (AIC for lm: 7902.1)Todas las variables son significativas, menos el intercepto.
sqrt(mean((connectivityM$infectados2020 - spatial_error_model$fitted.values)^2))## [1] 0.9592708jtools::export_summs(lm_model_alt, spatial_lag_model, spatial_error_model)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 0.59 | -1.18 *** | -0.42 |
| (0.32) | (0.29) | (0.32) | |
| gini2015 | 9.92 *** | 9.51 *** | 8.86 *** |
| (0.71) | (0.62) | (0.64) | |
| porcentaje_pob_pobreza | -2.40 *** | -1.02 *** | -1.67 *** |
| (0.23) | (0.21) | (0.25) | |
| porcentaje_pob_servicios_salud | 1.09 *** | 1.03 *** | 1.37 *** |
| (0.27) | (0.24) | (0.23) | |
| rezago_social | -0.48 *** | -0.31 *** | -0.33 *** |
| (0.04) | (0.04) | (0.05) | |
| popden2020 | 0.47 *** | 0.33 *** | 0.60 *** |
| (0.02) | (0.02) | (0.02) | |
| inclusion_fin_2019 | 0.68 *** | 0.63 *** | 0.48 *** |
| (0.04) | (0.03) | (0.03) | |
| porcentaje_pob_acceso_ss | -1.72 *** | -1.92 *** | -1.91 *** |
| (0.25) | (0.22) | (0.26) | |
| hospitalespHabitantes | -0.03 *** | -0.01 *** | -0.01 ** |
| (0.00) | (0.00) | (0.00) | |
| popnoafmed2015 | -0.02 *** | -0.02 *** | -0.01 ** |
| (0.00) | (0.00) | (0.00) | |
| crimen_2019 | 0.00 *** | 0.00 *** | 0.00 ** |
| (0.00) | (0.00) | (0.00) | |
| rho | 0.41 *** | ||
| (0.02) | |||
| lambda | 0.70 *** | ||
| (0.02) | |||
| N | 2453 | 2453 | 2453 |
| R2 | 0.71 | 0.78 | 0.82 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | |||
En todos los modelos, todas las variables fueron significativas
Seleccionando el mejor modelo
Estimation_Method<-c('OLS','WLS','SVM', 'RF', 'SAR', 'SEM')
RMSE<-c(1.24692, 1.244346, 1.029832, 0.9953379, 0.9592708, 1.059941)
RMSE_df<-data.frame(Estimation_Method,RMSE)
RMSE_df<-RMSE_df[order(RMSE),]
# knitr::kable(RMSE_df, format = "html")
export_table(RMSE_df, format = "md")| Estimation_Method | RMSE |
|---|---|
| SAR | 0.96 |
| RF | 1.00 |
| SVM | 1.03 |
| SEM | 1.06 |
| WLS | 1.24 |
| OLS | 1.25 |
Utilizando como medida el RMSE, concluimos que el mejor modelo es el SPATIAL ERROR MODEL.
Conclusiones y preguntas detonantes
¿Cuáles son los principales factores socioeconómicos que propician un incremento / Disminución de los casos confirmados de COVID-19? R: Hay muchas variables que explican fuertemente los casos de covid, así que tomando los resultados del modelo seleccionado, concluimos que las variables de mayor importancia son el GINI, Porcentaje de Pobreza, Porcentaje de Población sin Acceso a Servicios de Salud, Porcentaje de la Población sin Acceso a Seguro Social y la Tasa de Inclusión.
¿Cuáles son los principales factores socioeconómicos que propician un incremento / Disminución de los casos confirmados de COVID-19? R: Como era de esperarse, el porcentaje de pobreza afecta fuertemente a la tasa de COVID, sin embargo, al contrario de las expectativas, afecta negativamente, es decir, entre mayor sea el porcentaje de pobreza, menor serán los contagios de COVID. Así mismo pasa con el Porcentaje de Acceso a Seguro Social. Las variables que se comportan según lo esperado son la de GINI y Acceso a Servicios de Salud.