Regresión lineal múltiple
A continuación se presenta un caso de estudio sobre el efecto de la movilidad de las personas y su relación con el material particulado PM10, el cual actúa como contaminante en el aire.
library(pacman)
p_load("rmdformats","readr","readxl","ggplot2","plotly","DT","xfun","gridExtra","leaflet","TSstudio", "dplyr", "psych", "GGally", "corrplot", "tidyverse", "lubridate")Datos de concentración y movilidad
Los datos de contaminantes se obtuvieron de la estacion de calidad del aire de la ERNO del instituto de geologia de la unam ubicado en Hermosillo.
aire <- read_excel("Concentracion.xlsx")
CM <- read_excel("Concentracion_Mov.xlsx")
datatable(aire)Los datos de movilidad de las personas, salen de los celulares android de las personas, a través de sus ubicaciones GPS:
datatable(CM)Análisis de relación entre variables
A
datos <- as.data.frame(aire)
round(cor(x = datos, method="pearson"), 3)## PM10 Reactivacion_Comercial Supermercado_Farmacia
## PM10 1.000 0.281 0.280
## Reactivacion_Comercial 0.281 1.000 0.875
## Supermercado_Farmacia 0.280 0.875 1.000
## Parques_Centros 0.065 0.814 0.655
## Estaciones_Transito 0.270 0.944 0.872
## Lugares_Trabajo 0.100 0.581 0.492
## Residencia -0.238 -0.842 -0.729
## Parques_Centros Estaciones_Transito Lugares_Trabajo
## PM10 0.065 0.270 0.100
## Reactivacion_Comercial 0.814 0.944 0.581
## Supermercado_Farmacia 0.655 0.872 0.492
## Parques_Centros 1.000 0.770 0.371
## Estaciones_Transito 0.770 1.000 0.474
## Lugares_Trabajo 0.371 0.474 1.000
## Residencia -0.617 -0.758 -0.853
## Residencia
## PM10 -0.238
## Reactivacion_Comercial -0.842
## Supermercado_Farmacia -0.729
## Parques_Centros -0.617
## Estaciones_Transito -0.758
## Lugares_Trabajo -0.853
## Residencia 1.000Distribución de los datos
multi.hist(x = aire, dcol = c("green", "purple"), dlty = c("dotted", "solid")
)
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`.
Al analizar los datos y mostrarlos de forma gráfica, se puede observar que la variable que más correlación tiene con los niveles de PM10 es la de Reactivación Comercial, con un valor de 0.281, también se puede observar que Reactivación Comercial está altamente relacionada con variables como Supermercado-Farmacia (0.875), Parques-Centros (0.814), Estaciones-Tránsito (0.944) y Residencia (0.842), por lo cual conviene no agregarlos al análisis si se quiere seguir el principio de Parsimonia.
Primer modelo
Se probará un primer modelo utilizando únicamente la relación directa entre los niveles de PM10 y la Reactivación Comercial.
modelo <- lm(PM10 ~ Reactivacion_Comercial,data = datos)
summary(modelo)##
## Call:
## lm(formula = PM10 ~ Reactivacion_Comercial, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.189 -9.592 -2.002 8.228 59.423
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46.68846 1.50967 30.926 < 2e-16 ***
## Reactivacion_Comercial 0.26157 0.04237 6.173 1.51e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.22 on 445 degrees of freedom
## Multiple R-squared: 0.07887, Adjusted R-squared: 0.0768
## F-statistic: 38.1 on 1 and 445 DF, p-value: 1.512e-09Se obtiene una R cuadrada del 0.07887, lo cual indica que la relación entre las variables no es tan fuerte, sin emmbargo, el valor de significancia es de 1.5e-09, debido a que es menor a 0.05, se puede determinar que el resultado de y no es completamente al azar con los valores de x.
Por medio de esta validación, se obtiene la siguiente fórmula:
\[ y=46.68846+0.2515x \]
Representación grafica
plot(datos$Reactivacion_Comercial, datos$PM10)
abline(modelo)
Sustituyendo X en el modelo
y =46.68846+(0.2515*7)
y## [1] 48.44896Conocer el residuo del modelo
y-39.2058333## [1] 9.243127Primer modelo de regresión lineal
Se hará un analisis de un modelo involucrando todas las variabls para descartar completamente aquellas que no aportan un ajuste necesario al modelo
modelo2 <- lm(PM10 ~ Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Estaciones_Transito + Lugares_Trabajo + Residencia, data = datos )
summary(modelo2)##
## Call:
## lm(formula = PM10 ~ Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Estaciones_Transito + Lugares_Trabajo +
## Residencia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -65.754 -9.009 -1.347 8.087 52.265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.90828 2.09937 18.057 < 2e-16 ***
## Reactivacion_Comercial 0.66279 0.17806 3.722 0.000223 ***
## Supermercado_Farmacia 0.03773 0.13692 0.276 0.782998
## Parques_Centros -0.62999 0.09377 -6.719 5.7e-11 ***
## Estaciones_Transito -0.08419 0.12275 -0.686 0.493135
## Lugares_Trabajo -0.33070 0.08685 -3.808 0.000160 ***
## Residencia -0.85596 0.37221 -2.300 0.021935 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.37 on 440 degrees of freedom
## Multiple R-squared: 0.1881, Adjusted R-squared: 0.177
## F-statistic: 16.99 on 6 and 440 DF, p-value: < 2.2e-16Ha aumentado la R cuadrada con relación al anterior modelo, y el p value de 2.2e-16 indica que si es representativo. Sin embargo, solo puede explcar el 18.81% de las variables.
Evaluación general del modelo y selección de predictores
Se utilizará la estrategia de stepwise mixto, determinando la calidad del modelo con Akaike (AIC).
step(object = modelo2, direction = "both", trace = 1)## Start: AIC=2389.62
## PM10 ~ Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros +
## Estaciones_Transito + Lugares_Trabajo + Residencia
##
## Df Sum of Sq RSS AIC
## - Supermercado_Farmacia 1 15.7 90882 2387.7
## - Estaciones_Transito 1 97.2 90964 2388.1
## <none> 90867 2389.6
## - Residencia 1 1092.1 91959 2393.0
## - Reactivacion_Comercial 1 2861.5 93728 2401.5
## - Lugares_Trabajo 1 2994.5 93861 2402.1
## - Parques_Centros 1 9321.9 100189 2431.3
##
## Step: AIC=2387.7
## PM10 ~ Reactivacion_Comercial + Parques_Centros + Estaciones_Transito +
## Lugares_Trabajo + Residencia
##
## Df Sum of Sq RSS AIC
## - Estaciones_Transito 1 82.9 90965 2386.1
## <none> 90882 2387.7
## + Supermercado_Farmacia 1 15.7 90867 2389.6
## - Residencia 1 1088.1 91971 2391.0
## - Lugares_Trabajo 1 2990.1 93872 2400.2
## - Reactivacion_Comercial 1 3314.2 94197 2401.7
## - Parques_Centros 1 9917.2 100800 2432.0
##
## Step: AIC=2386.11
## PM10 ~ Reactivacion_Comercial + Parques_Centros + Lugares_Trabajo +
## Residencia
##
## Df Sum of Sq RSS AIC
## <none> 90965 2386.1
## + Estaciones_Transito 1 82.9 90882 2387.7
## + Supermercado_Farmacia 1 1.5 90964 2388.1
## - Residencia 1 1070.4 92036 2389.3
## - Lugares_Trabajo 1 2909.1 93874 2398.2
## - Reactivacion_Comercial 1 6151.6 97117 2413.4
## - Parques_Centros 1 9850.2 100816 2430.1##
## Call:
## lm(formula = PM10 ~ Reactivacion_Comercial + Parques_Centros +
## Lugares_Trabajo + Residencia, data = datos)
##
## Coefficients:
## (Intercept) Reactivacion_Comercial Parques_Centros
## 38.0294 0.5961 -0.6324
## Lugares_Trabajo Residencia
## -0.3199 -0.8469El modelo que arrojó AIC toma en cuenta aquel que relaciona los niveles de PM10 con las variables Reactivación Comercial, Parques-Centros, Lugares-Trabajo y Residencia, eliminando en este caso Estaciones-Tránsito y Supermercado-Farmacia.
Evaluación de modelo utilizando variables relevantes
modelo3 <- (lm(PM10 ~ Reactivacion_Comercial + Parques_Centros +
Lugares_Trabajo + Residencia, data = datos))
summary(modelo3)##
## Call:
## lm(formula = PM10 ~ Reactivacion_Comercial + Parques_Centros +
## Lugares_Trabajo + Residencia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -65.711 -9.082 -1.442 7.900 53.059
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.02938 1.88544 20.170 < 2e-16 ***
## Reactivacion_Comercial 0.59610 0.10903 5.467 7.66e-08 ***
## Parques_Centros -0.63242 0.09141 -6.918 1.61e-11 ***
## Lugares_Trabajo -0.31989 0.08509 -3.760 0.000193 ***
## Residencia -0.84686 0.37134 -2.281 0.023050 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.35 on 442 degrees of freedom
## Multiple R-squared: 0.1872, Adjusted R-squared: 0.1798
## F-statistic: 25.45 on 4 and 442 DF, p-value: < 2.2e-16El p-value sigue mostrando significancia del modelo, y en este caso se eliminaron dos variables sin comprometer de manera significativa la R cuadrada, puesto que ésta pasó de ser 0.1881 a 0.1872.
Intervalo de confianza
confint(lm(formula = PM10 ~ Reactivacion_Comercial + Parques_Centros + Lugares_Trabajo + Residencia, data = datos))## 2.5 % 97.5 %
## (Intercept) 34.3238269 41.7349256
## Reactivacion_Comercial 0.3818139 0.8103830
## Parques_Centros -0.8120734 -0.4527584
## Lugares_Trabajo -0.4871172 -0.1526720
## Residencia -1.5766584 -0.1170517Se puede observar cuánto afecta a los niveles de PM10 cada una de las variables si las demás se mantuveran constantes.
##Validación de condiciones para la regresión múltiple lineal
plot1 <- ggplot(data = datos, aes(PM10, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos, aes(Reactivacion_Comercial, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = datos, aes(Parques_Centros, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot4 <- ggplot(data = datos, aes(Lugares_Trabajo, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot5 <- 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)## `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'
Existe una linealidad general de todas las variables.
Distribución normal de los residuos
qqnorm(modelo3$residuals)
qqline(modelo3$residuals)
Se puede observar que los cuantiles de muestra y los cuantiles teóricos coinciden bastante en un rango de -2.5 y 2.
Prueba de Shapiro Wilk para determinar normalidad
shapiro.test(modelo3$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo3$residuals
## W = 0.95956, p-value = 9.746e-10El valor de W se acerca a 1, por lo cual se llega a la conclusión de 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'
Conclusión
Se puede observar que hay un aumento de la contaminación PM10 junto con la reactivación económica, lo cual puede ser debido a que la circulación de vehículos aumenta al volver a abrir negocios, lo cual en un periodo de cuarentena no sucede.