Regresión lineal múltiple
Caso de estudio de análisis del efecto de la movilidad de personas en la contaminación atmosférica
library(pacman)
p_load(rmdformats,readr,readxl,ggplot2,plotly,DT,xfun,gridExtra,leaflet,TSstudio, dplyr, psych, GGally, corrplot)Importar datos
Datos de concentración y movilidad sin fechas
aire <- read_excel("Concentracion.xlsx")
CM <- read_excel("Concentracion_Mov.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 \((O_3)\)
- 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
Los datos de movilidad de las personas, salen de los celulares android de las personas, a través de sus ubicaciones GPS: https://www.google.com/covid19/mobility/
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)
mLa 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
Introducción a 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 (X1, X2, X3…).
Es una extensión de la regresión lineal simple, por lo que es fundamental comprender esta última. 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).
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} \]
$ _{0} $ : es la ordenada en el origen, el valor de la variable dependiente Y cuando todos los predictores son cero.
$ _{i} $ : es el efecto promedio que tiene el incremento en una unidad de la variable predictora Xi sobre la variable dependiente Y manteniéndose constantes el resto de variables. Se conocen como coeficientes parciales de regresión.
$ e_{i} $ : es el residuo o error, la diferencia entre el valor observado y el estimado por el modelo.
Condiciones para la regresión lineal múltiple
Los modelos de correlación lineal múltiple requieren de las mismas condiciones que los modelos lineales simples más otras adicionales.
Parsimonia
Este término hace referencia a que el mejor modelo es aquel capaz de explicar con mayor precisión la variabilidad observada en la variable respuesta empleando el menor número de predictores, por lo tanto, con menos asunciones.
Tamaño de la muestra
No se trata de una condición de por sí pero, si no se dispone de suficientes observaciones, predictores que no son realmente influyentes podrían parecerlo. En el libro Hanbook of biological statistics recomiendan que el número de observaciones sea como mínimo entre 10 y 20 veces el número de predictores del modelo.
La gran mayoría de condiciones se verifican utilizando los residuos, por lo tanto, se suele generar primero el modelo y posteriormente validar las condiciones. De hecho, el ajuste de un modelo debe verse como un proceso iterativo en el que se ajusta el modelo, se evalúan sus residuos y se mejora. Así hasta llegar a un modelo óptimo.
Análisis de correlación
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.000Distribución de los datos
multi.hist(x = aire, dcol = c("blue", "red"), dlty = c("dotted", "solid")
)
## Analizando la dispersión de los datos con ggplot y 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 correlación
modelo <- lm(O3 ~ PM10,data = datos)
summary(modelo)##
## Call:
## lm(formula = O3 ~ PM10, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.0817 -5.0899 -0.6656 5.7951 25.8668
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.04971 1.01381 28.654 < 2e-16 ***
## PM10 -0.12090 0.02436 -4.964 9.87e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.148 on 445 degrees of freedom
## Multiple R-squared: 0.05246, Adjusted R-squared: 0.05033
## F-statistic: 24.64 on 1 and 445 DF, p-value: 9.865e-07La ecuación de la recta de los mínimos cuadrados tomaría la siguiente forma:
\[ y = 29.04971-0.12090x \]
Esta ecuación representa esta recta ajustada a la dispersión de los puntos con los cuales se forma el modelo
plot(datos$PM10, datos$O3)
abline(modelo)
Ahora substituyendo X en nuestro modelo
O3m = 29.04971 - (0.12090*40.76000)
O3m## [1] 24.12183Residuo
$ Datos real - Dato simulado por modelo $
26.91333 - 24.12183## [1] 2.7915En este caso al modelo le “faltan” 2.7915 para llegar a el valor de 26.91333, por lo tanto esta subestimando.
Nota: Pese a que las PM10 NO causan directamente (químicamente) la reacción que da como resultado el O3
Primer modelo de regresión lineal
De que forma las variables predictoras (X) influyen positivamente o negativamente en el valor de SO2
modelo2 <- lm(SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Estaciones_Transito + Lugares_Trabajo + Residencia, data = datos )
summary(modelo2)##
## 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 ( SO2)
Evaluación general 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 = modelo2, 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.035005Con respecto al método AIC, el mejor modelo resultante del proceso de selección ha sido:
modelo3 <- (lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Lugares_Trabajo + Residencia, data = datos))
summary(modelo3)##
## 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:
Intervalo de confianza: Un intervalo de confianza es una técnica de estimación utilizada en inferencia estadística que permite acotar un par o varios pares de valores, dentro de los cuales se encontrará la estimación puntual buscada (con una determinada probabilidad).
Un intervalo de confianza nos va a permitir calcular dos valores alrededor de una media muestral (uno superior y otro inferior). Estos valores van a acotar un rango dentro del cual, con una determinada probabilidad, se va a localizar el parámetro poblacional.
ntervalo de confianza = media +- margen de error
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.026801239Análisis de residuales
Un residual ri es la diferencia entre el valor observado Yi y el valor estimado por la linea de ˆˆ regresión Yi , es decir, ri = Yi − Yi . El residual puede ser considerado como el error aleatorio ei observado. También se acostumbra usar el Residual estandarizado,el cual se obtiene al dividir el residual entre la desviación estándar.
Validació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'
Distribución normal de los residuos
qqnorm(modelo3$residuals)
qqline(modelo3$residuals)
Prueba de normalidad de Shapiro-wilk
shapiro.test(modelo3$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo3$residuals
## W = 0.99227, p-value = 0.02042Tanto el análisis gráfico como es test de hipótesis confirman la normalidad.
Variabilidad constante de los residuos (homocedasticidad):
Al representar los residuos frente a los valores ajustados por el modelo, los primeros se tienen que distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X. Si se observa algún patrón específico, por ejemplo forma cónica o mayor dispersión en los extremos, significa que la variabilidad es dependiente del valor ajustado y por lo tanto no hay homocedasticidad.
corrplot(cor(dplyr::select(datos, O3, PM10, Reactivacion_Comercial, Supermercado_Farmacia, Parques_Centros, Lugares_Trabajo, Residencia)),
method = "number", tl.col = "black")
¿Por qué las variables se correlacionan de esta forma?
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)Análisis de serie de tiempo
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)'))Descarga 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")