Regresión lineal múltiple
Caso de estudio de análisis del efecto de la movilidad de personas en la contaminación atmosférica
setwd("~/Universidad en general/4to semestre/Estadistica aplicada/ea9am")
library(pacman)
p_load(rmdformats,readr,readxl,ggplot2,plotly,DT,xfun,gridExtra,leaflet,TSstudio, dplyr, psych, GGally, corrplot,xlsx,lmtest,corrplot,car)Caso de estudio
Caso de estudio de análisis del efecto de la movilidad de personas en la contaminación atmosférica y su relacion con el Dioxido de azufre
Actualmente gracias a la propagación de información el cuidado del medio ambiente ha sido un tema que se empezó a mencionar cada vez más, con la intencion de que la gente se de cuenta de los daños que puede ocacionar cada cosita que hacemos en nuestro dia a dia, aparte de las grandes empresas que ocacionan una cantidad increible de daño al medio ambiente con el fin de continuar con los productos que consumimos en nuestra vida diaria al momemento de realizar los productos generan una gran cantidad de gases tales como el dióxido de azufre(SO2) y los excesos de humo y material articulado(PM).
Estos gases nos afectan de una forma muy dañina a nosotras las personas, ya que la contaminación genera problemas tales como cáncer de pulmón, asma entre otras enfermedades respiratorias por el simple hecho de continuar realizando nuestra vida diaria, por lo que si se reduce el uso de multiples artefactos y se mantiene un control de contaminación de empresas la vida del ser humano será mejor en el ámbito de la salud respiratoria.
Los principales contaminantes de los cuales informaremos y nos causan una gran cantidad de daño a nosotros y que han sido creados a causa de las actividades diarias y económicas de los seres humanos son:
- 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:
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. ## Datos
Datos de concentración y movilidad sin fechas
aire <- read_excel("~/Universidad en general/4to semestre/Estadistica aplicada/ea9am/Concentracion.xlsx")
CM <- read_excel("~/Universidad en general/4to semestre/Estadistica aplicada/ea9am/Concentracion_Mov.xlsx")
datatable(aire)Análisis de correlación entre variables
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 en histogramas
multi.hist(x = aire, dcol = c("purple", "pink"), 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`.
Los datos analizados hasta el momento muestran que la variable que mas correlacion tiene con el SO2 es PM10(r=0.612), reactivacion comercial(r=0.404) y supermercado farmacia (r=0.443)
Modelo general de correlación
En esta parte se genera el primer modelo de regresion lineal simple con la ecuacion recta de los minimos cuadrados. Para realizar este primer modelo utilizamos la variable que mas relacion tiene con el SO2 siendo PM10.
modelo <- lm(SO2 ~PM10,data=datos)
summary(modelo)##
## Call:
## lm(formula = SO2 ~ PM10, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.68827 -0.18494 0.03467 0.15767 0.89115
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2339600 0.0279254 -8.378 7.1e-16 ***
## PM10 0.0109437 0.0006709 16.311 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2244 on 445 degrees of freedom
## Multiple R-squared: 0.3742, Adjusted R-squared: 0.3728
## F-statistic: 266.1 on 1 and 445 DF, p-value: < 2.2e-16La ecuación de la recta de los mÃnimos cuadrados tomarÃa la siguiente forma:
\[ y = -0.2339600 - 0.0109437x \] Esta formula resultante representa la recta que se ajusta a la dispersion
Ejemplo grafico del modelo
plot(datos$PM10, datos$SO2)
abline(modelo)
Ahora substituyendo X en nuestro modelo
y = -0.2339600 - (0.0109437*39.2058333)Residuo
$ Datos real - Dato simulado por modelo $
y-0.19791667## [1] -0.8609335Primer modelo de regresión lineal
De que forma las variables predictoras (X) 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.
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 el 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
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.026801239Se puede analizar que tanto afecta al CO2 las demas variables con una certeza del 97.5%
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.02042El valor es cercano a uno por lo cual se concluye que el modelo muestra normalidad
Variabilidad constante de los residuos
ggplot(data = datos, aes(modelo$fitted.values, modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'Linea múltiple
bptest(modelo3)##
## studentized Breusch-Pagan test
##
## data: modelo3
## BP = 19.372, df = 7, p-value = 0.0071Matriz de correlación entre predictores para el modelo 3
corrplot(cor(dplyr::select(datos,O3,PM10,Reactivacion_Comercial,Supermercado_Farmacia,Parques_Centros,Lugares_Trabajo,Residencia)),
method = "number", tl.col = "black")
vif(modelo3)## O3 PM10 Reactivacion_Comercial
## 1.198005 1.304100 13.623986
## Supermercado_Farmacia Parques_Centros Lugares_Trabajo
## 4.815401 3.618397 5.084395
## Residencia
## 11.372743Autocorrelacion
dwt(modelo, alternative = "two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 0.8228273 0.3457909 0
## Alternative hypothesis: rho != 0##Identificacion de los posibles valores atipicos o influyentes
datos$studentized_residual <- rstudent(modelo3)
ggplot(data = datos, aes(x = predict(modelo3), y = abs(studentized_residual))) +
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
# se identifican en rojo observaciones con residuos estandarizados absolutos > 3
geom_point(aes(color = ifelse(abs(studentized_residual) > 3, 'red', 'black'))) +
scale_color_identity() +
labs(title = "Distribución de los residuos studentized",
x = "predicción modelo") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5))
Se presentan unas cuantas observaciones atipicas
which(abs(datos$studentized_residual) > 3)## [1] 27 159 295 296Cuantificacion de las influencias de los predictores
summary(influence.measures(modelo))## Potentially influential observations of
## lm(formula = SO2 ~ PM10, data = datos) :
##
## dfb.1_ dfb.PM10 dffit cov.r cook.d hat
## 3 0.17 -0.23 -0.26_* 0.99 0.03 0.01
## 5 0.16 -0.22 -0.24_* 1.00 0.03 0.01_*
## 27 -0.10 0.09 -0.10 1.01_* 0.00 0.01
## 33 -0.04 0.03 -0.04 1.02_* 0.00 0.01
## 56 0.02 -0.02 0.02 1.02_* 0.00 0.01
## 105 0.14 -0.21 -0.26_* 0.97_* 0.03 0.01
## 107 -0.05 0.01 -0.11 0.98_* 0.01 0.00
## 108 0.01 -0.06 -0.14 0.98_* 0.01 0.00
## 137 0.03 -0.07 -0.13 0.99_* 0.01 0.00
## 219 0.17 -0.22 -0.23_* 1.01 0.03 0.02_*
## 228 0.29 -0.36 -0.38_* 1.01 0.07 0.03_*
## 255 0.21 -0.26 -0.28_* 1.00 0.04 0.02_*
## 256 0.13 -0.16 -0.17 1.02_* 0.01 0.02_*
## 260 0.06 -0.08 -0.08 1.02_* 0.00 0.02_*
## 263 0.03 -0.05 -0.05 1.01_* 0.00 0.01
## 277 0.06 -0.08 -0.08 1.03_* 0.00 0.02_*
## 278 0.04 -0.05 -0.05 1.02_* 0.00 0.02_*
## 279 0.01 -0.01 -0.01 1.01_* 0.00 0.01
## 291 0.01 -0.02 -0.02 1.02_* 0.00 0.02_*
## 292 -0.01 0.01 0.01 1.02_* 0.00 0.01
## 294 0.06 -0.02 0.12 0.98_* 0.01 0.00
## 295 0.12 -0.08 0.16 0.97_* 0.01 0.00
## 296 -0.01 0.07 0.15 0.97_* 0.01 0.00
## 298 -0.05 0.07 0.08 1.02_* 0.00 0.01
## 314 -0.23 0.28 0.29_* 1.02_* 0.04 0.03_*
## 315 -0.10 0.13 0.15 1.01 0.01 0.01_*
## 329 -0.03 0.04 0.05 1.02_* 0.00 0.01
## 335 -0.08 0.10 0.11 1.02_* 0.01 0.02_*
## 340 0.02 -0.03 -0.03 1.01_* 0.00 0.01
## 371 0.03 0.02 0.12 0.98_* 0.01 0.00
## 372 0.01 0.04 0.12 0.98_* 0.01 0.00
## 378 0.15 -0.19 -0.19 1.04_* 0.02 0.04_*
## 379 -0.04 0.05 0.05 1.03_* 0.00 0.03_*
## 397 0.73 -0.71 0.73_* 0.96_* 0.26 0.03_*