Caso de estudio de análisis del efecto de la movilidad de personas en la contaminación atmosférica y la incidencia de casos de contagio de COVID-19
Paquetes
library(pacman)
p_load(rmdformats,readr,readxl,ggplot2,plotly,DT,xfun,gridExtra,leaflet,TSstudio, dplyr, psych, GGally)Importar datos
Datos de concentracion y movilidad sin fechas
aire <- read_excel("Concentracion.xlsx")
datatable(aire)setwd("~/Documents/ESTADISTICA")
CM <- read_excel("Concentracion_Mov.xlsx")
datatable(CM)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 \((O_3)\)
- Dióxido de nitrógeno \((NO_2)\)
- Dióxido de azufre \((SO_2)\)
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:
Ubicación de donde se obtuvieron 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
Par que esto funcione necesitamos tener instalada y actividad la libreria pysch
multi.hist(x = aire, dcol = c("blue", "red"), dlty = c("dotted", "solid"),
main = "")
Analizando la dispersion de los datos
Esto es realizado con la libreria GGally
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
Intenta explicar la correlacion que existe del SO2 (Y) dependiente de las demas variables X
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\)
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.
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
library(gridExtra)
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.02042Matriz de correlación entre predictores.
library(corrplot)## corrplot 0.88 loadedcorrplot(cor(dplyr::select(datos, O3, PM10, Reactivacion_Comercial, Supermercado_Farmacia, Parques_Centros, Lugares_Trabajo, Residencia)),
method = "number", tl.col = "black")
Movilidad de personas en Sonora
Para esto estamos usando datos de google mobility report, que pueden ser encontrados con su documentacion aqui: https://www.google.com/covid19/mobility/
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)Contaminantes 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
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
- Concentracion y movilidad
xfun::embed_file("Concentracion_Mov.xlsx")Download Concentracion_Mov.xlsx
- Concentracion y movilidad sin fechas
xfun::embed_file("Concentracion.xlsx")