Análisis de datos de COVID-19 en México, impacto de la movilidad de las personas sobre los contagios, impacto sobre la movilidad urbana sobre la calidad del aire
Paquetes a utilizar
library(pacman)## Warning: package 'pacman' was built under R version 3.6.3p_load(rmdformats, readr, readxl, ggplot2, plotly, DT, xfun, gridExtra, leaflet, GGally, psych, corrplot)Importar datos
CM <- read_excel("Concentracion_Mov.xlsx")
datatable(CM)#Los numeros que salen en reactivación y farmacias y demas son un porcentaje superior o inferior a la media que hay, por eso si es -x es porque hay -x porcentaje de personas en ese sectorDatos de concentracion y movilidad sin fechas
aire <- read_excel("Concentracion.xlsx")
datatable(aire)Caso de estudio
La gestión de la contaminación del aire se ha vuelto importante en el último periodo debido a las siguientes afirmaciones: “A causa del desarrollo industrial y avance tecnológico se estima aproximadamente 1200 millones de personas están expuestas a niveles de dióxido de azufre (SO2), muy por encima de por directrices de la Organización Mundial de la Salud (OMS) y aproximadamente 1400 millones de personas están expuestas a niveles excesivos de humo y material articulado (PM)” (Rico, 2018).
La importancia de respirar aire limpio sin contaminantes es crucial para todo ser viviente. Por lo cual es importante que este sea aire limpio, pero ¿Cómo detectar si el aire que respiramos no tiene contaminantes? ¿Cómo nos afecta a nosotros?. La calidad del aire está directamente relacionada con la movilidad de carros, autobuses y hasta la propia movilidad de la ciudad, donde las personas realizan sus actividades diarias, por lo tanto, entre más movilidad haya en la ciudad, aumentan los contaminantes en el aire.
Los principales contaminantes en el aire han sido creados principalmente por la actividad económica humana, la concentración de estas sustancias es altamente nociva para la salud del ser humano. Este estudio tuvo en cuenta los siguientes contaminantes:
- Material particulado \((PM10)\)
- Ozono \((O3)\)
- Dióxido de nitrógeno \((NO2)\)
- Dióxido de azufre \((SO2)\)
Fuente de de los datos
Los datos de contaminantes atmosfericos provienen de la estacion de calidad del aire de la ERNO del instituto de geologia de la unam ubicado en Hermosillo. Podemos visitar el origen de los datos aqui: http://www.erno.geologia.unam.mx La estación de calidade del aire se encuentra ubicada en el siguiente mapa interactivo:
Origen de los datos
content <- paste(sep = "<br/>",
"<b><a href='https://www.ruoa.unam.mx/index.php?page=estaciones&id=6#datos'>ERNO</a></b>","Lng: -110.9706, Lat: 29.0814")
m <- leaflet() %>%
addTiles() %>%
addMarkers(lng=-110.9706, lat= 29.0814, popup= content)
m- La contaminación del aire es el principal riesgo ambiental para la salud pública en las Américas.
- En todo el mundo, cerca de 7 millones de muertes prematuras fueron atribuibles a la contaminación del aire ambiental en 2016. Alrededor del 88% de estas muertes ocurren en países de ingresos bajos y medios.
- Más de 150 millones de personas en América Latina viven en ciudades que exceden las Guías de Calidad del Aire de la OMS.
- La exposición a altos niveles de contaminación del aire puede causar una variedad de resultados adversos para la salud: aumenta el riesgo de infecciones respiratorias, enfermedades cardiacas, derrames cerebrales y cáncer de pulmón, las cuales afectan en mayor proporción a población vulnerable, niños, adultos mayores y mujeres.
- La contaminación del aire en el hogar se asocia al uso de combustibles y prácticas de cocina ineficiente
Regresión Lineal Múltiple
La regresión lineal múltiple permite generar un modelo lineal en el que el valor de la variable dependiente o respuesta (Y) se determina a partir de un conjunto de variables independientes llamadas predictores \((X_1, X_2, X_3…)\). Los modelos de regresión múltiple pueden emplearse para predecir el valor de la variable dependiente o para evaluar la influencia que tienen los predictores sobre ella (esto último se debe que analizar con cautela para no malinterpretar causa-efecto). (Rodrigo, 2016)
Los modelos lineales múltiples siguen la siguiente ecuación:
\[ Y_{i}=(\beta_{0}+\beta_{1}X_{1i}+\beta_{2}X_{2i}+\cdots+\beta_{n}X_{ni})+e_{i} \]
datos <- as.data.frame(aire)
round(cor(x = datos, method = "pearson"), 3)## O3 SO2 PM10 Reactivacion_Comercial
## O3 1.000 -0.220 -0.229 -0.193
## SO2 -0.220 1.000 0.612 0.404
## PM10 -0.229 0.612 1.000 0.281
## Reactivacion_Comercial -0.193 0.404 0.281 1.000
## Supermercado_Farmacia -0.047 0.443 0.280 0.875
## Parques_Centros -0.182 0.019 0.065 0.814
## Estaciones_Transito -0.174 0.399 0.270 0.944
## Lugares_Trabajo -0.119 0.138 0.100 0.581
## Residencia 0.130 -0.371 -0.238 -0.842
## Supermercado_Farmacia Parques_Centros
## O3 -0.047 -0.182
## SO2 0.443 0.019
## PM10 0.280 0.065
## Reactivacion_Comercial 0.875 0.814
## Supermercado_Farmacia 1.000 0.655
## Parques_Centros 0.655 1.000
## Estaciones_Transito 0.872 0.770
## Lugares_Trabajo 0.492 0.371
## Residencia -0.729 -0.617
## Estaciones_Transito Lugares_Trabajo Residencia
## O3 -0.174 -0.119 0.130
## SO2 0.399 0.138 -0.371
## PM10 0.270 0.100 -0.238
## Reactivacion_Comercial 0.944 0.581 -0.842
## Supermercado_Farmacia 0.872 0.492 -0.729
## Parques_Centros 0.770 0.371 -0.617
## Estaciones_Transito 1.000 0.474 -0.758
## Lugares_Trabajo 0.474 1.000 -0.853
## Residencia -0.758 -0.853 1.000histogramas de los datos
multi.hist(x = aire, dcol = c("blue", "red"), dlty = c("dotted", "solid"),
main = "")
Analizando la dispersion de los datos
ggpairs(aire, lower = list(continuous = "smooth"),
diag = list(continuous = "barDiag"), axisLabels = "none")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Modelo general de correlacion
# lm(variable de respuesta (y) ~ x´s)
modelo <- lm(SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Estaciones_Transito + Lugares_Trabajo + Residencia, data = datos )
summary(modelo)##
## Call:
## lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Estaciones_Transito + Lugares_Trabajo +
## Residencia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51190 -0.09921 -0.00269 0.09870 0.52906
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0187954 0.0380324 -0.494 0.621416
## O3 -0.0057616 0.0009832 -5.860 9.10e-09 ***
## PM10 0.0056324 0.0005408 10.415 < 2e-16 ***
## Reactivacion_Comercial 0.0096304 0.0020097 4.792 2.27e-06 ***
## Supermercado_Farmacia 0.0060624 0.0015668 3.869 0.000126 ***
## Parques_Centros -0.0169489 0.0010860 -15.607 < 2e-16 ***
## Estaciones_Transito -0.0018481 0.0013556 -1.363 0.173465
## Lugares_Trabajo -0.0106817 0.0009837 -10.858 < 2e-16 ***
## Residencia -0.0352798 0.0041747 -8.451 4.32e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1582 on 438 degrees of freedom
## Multiple R-squared: 0.6938, Adjusted R-squared: 0.6882
## F-statistic: 124.1 on 8 and 438 DF, p-value: < 2.2e-16El modelo con todas las variables introducidas como predictores tiene un R2 ajustada (0.6882), es capaz de explicar el 68.82% de la variabilidad observada en el Dióxido de Azufre \((SO_2)\).
Analisis explratorio de datos
Empezamos este análisis utulizando una matriz de diagramas de dispersión
Cual de las variables de movilidad afectan mas a la incidencia de contaminantes atmosféricos?
Nuestra hipótesis es que a medida que las personas aumemtan su estadía en sus casas (y por consecuente no usan su automovil o el transporte público) los contaminantes atmosféricos bajan su concentración
Evaluacion Generacion del modelo
El criterio de información de Akaike (AIC) es una medida de la calidad relativa de un modelo estadístico, para un conjunto dado de datos. Como tal, el AIC proporciona un medio para la selección del modelo. AIC maneja un trade-off entre la bondad de ajuste del modelo y la complejidad del modelo. Se basa en la entropía de información: se ofrece una estimación relativa de la información perdida cuando se utiliza un modelo determinado para representar el proceso que genera los datos. AIC no proporciona una prueba de un modelo en el sentido de probar una hipótesis nula, es decir AIC no puede decir nada acerca de la calidad del modelo en un sentido absoluto. Si todos los modelos candidatos encajan mal, AIC no dará ningún aviso de ello.
# AIC evalua que tanto se relacionan en forma de tendencia, son negativos porque el contaminante baja entre menos movilidad
step(object = modelo, direction = "both", trace = 1)## Start: AIC=-1639.29
## SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Estaciones_Transito + Lugares_Trabajo +
## Residencia
##
## Df Sum of Sq RSS AIC
## - Estaciones_Transito 1 0.0465 11.014 -1639.4
## <none> 10.968 -1639.3
## - Supermercado_Farmacia 1 0.3749 11.343 -1626.3
## - Reactivacion_Comercial 1 0.5750 11.543 -1618.5
## - O3 1 0.8600 11.828 -1607.5
## - Residencia 1 1.7883 12.756 -1573.8
## - PM10 1 2.7161 13.684 -1542.4
## - Lugares_Trabajo 1 2.9523 13.920 -1534.7
## - Parques_Centros 1 6.0997 17.067 -1443.6
##
## Step: AIC=-1639.4
## SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Lugares_Trabajo + Residencia
##
## Df Sum of Sq RSS AIC
## <none> 11.014 -1639.4
## + Estaciones_Transito 1 0.0465 10.968 -1639.3
## - Supermercado_Farmacia 1 0.3298 11.344 -1628.2
## - Reactivacion_Comercial 1 0.6091 11.623 -1617.3
## - O3 1 0.8366 11.851 -1608.7
## - Residencia 1 1.7647 12.779 -1575.0
## - PM10 1 2.7569 13.771 -1541.5
## - Lugares_Trabajo 1 2.9217 13.936 -1536.2
## - Parques_Centros 1 6.0965 17.111 -1444.5##
## Call:
## lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Lugares_Trabajo + Residencia, data = datos)
##
## Coefficients:
## (Intercept) O3 PM10
## -0.020842 -0.005669 0.005668
## Reactivacion_Comercial Supermercado_Farmacia Parques_Centros
## 0.008020 0.005439 -0.016944
## Lugares_Trabajo Residencia
## -0.010417 -0.035005El mejor modelo resultante del proceso de selección ha sido:
modelo <- (lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Lugares_Trabajo + Residencia, data = datos))
summary(modelo)##
## Call:
## lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Lugares_Trabajo + Residencia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5458 -0.1006 -0.0022 0.0964 0.5374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0208416 0.0380399 -0.548 0.584047
## O3 -0.0056694 0.0009818 -5.775 1.46e-08 ***
## PM10 0.0056679 0.0005407 10.482 < 2e-16 ***
## Reactivacion_Comercial 0.0080203 0.0016277 4.927 1.18e-06 ***
## Supermercado_Farmacia 0.0054386 0.0015000 3.626 0.000322 ***
## Parques_Centros -0.0169444 0.0010870 -15.588 < 2e-16 ***
## Lugares_Trabajo -0.0104169 0.0009653 -10.791 < 2e-16 ***
## Residencia -0.0350045 0.0041739 -8.387 6.89e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1584 on 439 degrees of freedom
## Multiple R-squared: 0.6925, Adjusted R-squared: 0.6876
## F-statistic: 141.2 on 7 and 439 DF, p-value: < 2.2e-16Es recomendable mostrar el intervalo de confianza para cada uno de los coeficientes parciales de regresión:
confint(lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Lugares_Trabajo + Residencia, data = datos))## 2.5 % 97.5 %
## (Intercept) -0.095604583 0.053921343
## O3 -0.007599002 -0.003739801
## PM10 0.004605203 0.006730588
## Reactivacion_Comercial 0.004821182 0.011219376
## Supermercado_Farmacia 0.002490524 0.008386774
## Parques_Centros -0.019080816 -0.014808020
## Lugares_Trabajo -0.012314119 -0.008519681
## Residencia -0.043207808 -0.026801239Validación de condiciones para la regresión múltiple lineal
plot1 <- ggplot(data = datos, aes(O3, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos, aes(PM10, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = datos, aes(Reactivacion_Comercial, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot4 <- ggplot(data = datos, aes(Supermercado_Farmacia, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot5 <- ggplot(data = datos, aes(Parques_Centros, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot6 <- ggplot(data = datos, aes(Lugares_Trabajo, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot7 <- ggplot(data = datos, aes(Residencia, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2, plot3, plot4, plot5, plot6, plot7)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Se puede concluir que si se cumple la linealidad para todos los predictores mostrando una buena relación con las concentraciones respecto al movilidad, comprobando su alta correlacción.
Distribución normal de los residuos:
qqnorm(modelo$residuals)
qqline(modelo$residuals)
Prueba de normalidad de los datos residuales
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.99227, p-value = 0.02042Los datos no son normales ya que p < 0.05, así que se aprueba la hipótesis nula porque los datos no son normales, rechazando la hipótesis alternativa.
Matriz de correlación entre predictores.
corrplot(cor(dplyr::select(datos, O3, PM10, Reactivacion_Comercial, Supermercado_Farmacia, Parques_Centros, Lugares_Trabajo, Residencia)),
method = "number", tl.col = "black")
Serie de tiempo de los datos
Para entender como serie de tiempo la movilidad de google vs los contaminantes de la estación ERNO, haremos un gráfico interactivo con GGPLOT2
movilidad <- ggplot(CM) +
geom_line(aes(x=Fecha,y=Reactivacion_Comercial,colour="Recreación y comercio")) +
geom_line(aes(x=Fecha,y=Supermercado_Farmacia,colour="Supermercados y farmacias")) +
geom_line(aes(x= Fecha,y= Parques_Centros, colour= "Parques")) +
geom_line(aes(x= Fecha,y= Estaciones_Transito, colour= "Estaciones de tránsito")) +
geom_line(aes(x= Fecha,y= Lugares_Trabajo,colour= "Lugares de trabajo")) +
geom_line(aes(x= Fecha,y= Residencia,colour= "Lugares residenciales")) +
labs(title= "Reporte de movilidad",x= "Fecha",y= "Procentaje de cambio de movilidad")
ggplotly(movilidad)#Si sube transito O3 sube, si sube la gente en su casa O3 bajaContaminantes atmosfericos para Hermosillo
plot_ly(CM,colors = rainbow(3)) %>%
add_lines(x = ~Fecha, y = ~O3,mode="lines",name = "O3") %>%
add_lines(x = ~Fecha, y = ~SO2,mode="lines", name = "SO2") %>%
add_lines(x = ~Fecha, y = ~PM10,mode="lines", name ="PM10") %>%
rangeslider() %>%
layout(title = 'Contaminantes atmosféricos (concentraciones)',
xaxis = list(title = 'Fecha'),
yaxis = list(title = 'Concentración (ppb para O3 y SO2, ug/m3 para PM10)'))Entender la relacion entre la dispersion cruzada de los datos
Matriz de diagramas de dispersion
#03 esta por partes por billón y es un contaminante traza, lo que quiere decir que es una concentración muy baja con efectos muy adversicos = muy poco Ozono mucho daño
#PM10 son basura que estan por el aire que luego nosotros respiramos y se miden por partes por millón, de ahí su nombre
pairs(CM)
Matriz de coeficientes de correlacion pearson
cor(aire)## O3 SO2 PM10
## O3 1.00000000 -0.22046192 -0.22904955
## SO2 -0.22046192 1.00000000 0.61169793
## PM10 -0.22904955 0.61169793 1.00000000
## Reactivacion_Comercial -0.19313979 0.40436153 0.28084230
## Supermercado_Farmacia -0.04677236 0.44278531 0.28031731
## Parques_Centros -0.18201172 0.01908292 0.06546006
## Estaciones_Transito -0.17411523 0.39928523 0.26992440
## Lugares_Trabajo -0.11913738 0.13779670 0.09976078
## Residencia 0.13008239 -0.37100732 -0.23779333
## Reactivacion_Comercial Supermercado_Farmacia
## O3 -0.1931398 -0.04677236
## SO2 0.4043615 0.44278531
## PM10 0.2808423 0.28031731
## Reactivacion_Comercial 1.0000000 0.87496681
## Supermercado_Farmacia 0.8749668 1.00000000
## Parques_Centros 0.8137739 0.65541398
## Estaciones_Transito 0.9444005 0.87219787
## Lugares_Trabajo 0.5809903 0.49241600
## Residencia -0.8420311 -0.72929468
## Parques_Centros Estaciones_Transito Lugares_Trabajo
## O3 -0.18201172 -0.1741152 -0.11913738
## SO2 0.01908292 0.3992852 0.13779670
## PM10 0.06546006 0.2699244 0.09976078
## Reactivacion_Comercial 0.81377389 0.9444005 0.58099030
## Supermercado_Farmacia 0.65541398 0.8721979 0.49241600
## Parques_Centros 1.00000000 0.7703343 0.37121486
## Estaciones_Transito 0.77033431 1.0000000 0.47409617
## Lugares_Trabajo 0.37121486 0.4740962 1.00000000
## Residencia -0.61722623 -0.7576539 -0.85264903
## Residencia
## O3 0.1300824
## SO2 -0.3710073
## PM10 -0.2377933
## Reactivacion_Comercial -0.8420311
## Supermercado_Farmacia -0.7292947
## Parques_Centros -0.6172262
## Estaciones_Transito -0.7576539
## Lugares_Trabajo -0.8526490
## Residencia 1.0000000Descarga los datos
- Datos de concentracion de contaminantes y de movilidad
xfun::embed_file("Concentracion_Mov.xlsx")Download Concentracion_Mov.xlsx
- Concentracion y movilidad sin fechas
xfun::embed_file("Concentracion.xlsx")