Regresión lineal múltiple
setwd("~/ea9AM")
library(pacman)
p_load("prettydoc","DT","xfun", "dplyr", "psych", "GGally", "ggplot2", "readr", "gridExtra", "rmdformats", "readxl", "plotly", "leaflet","TSstudio", "lmtest","car", "corrplot")
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 X1X2X3.
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 para el caso de estudio (Se muestran los contaminantes PM10 // O3 // SO2)
En la actualidad, el cuidado del medio ambiente ha sido un tema en tendencia gracias a la propagación masiva de información en redes sociales, con la intención de concientizar a la población acerca de este tema y así provocar un cambio en la huella ecológica de todas las personas.
Sin embargo, es importante considerar que la contaminación del aire es un riesgo medioambiental latente para la salud del ser humano. A través de la reducción de los niveles de contaminación del aire, los países son capaces de disminuir la carga de morbilidad respecto a accidentes cerebrovasculares, cáncer de pulmón y neumopatías crónicas y agudas, por ejemplo el asma. Mientras más bajos sean los niveles de contaminación del aire, será mejor la salud cardiovascular y respiratoria de la población, en corto y largo plazo.
Al momento de hacer mediciones y evaluaciones para determinar el impacto de la contaminación del aire en la población y en los recursos naturales, es muy necesario tener sistemas, redes y programas pertinentes para medición de la calidad del aire.
Los principales contaminantes se han generado debido a la huella ecologica que la humanidad ha dejado en el planeta tierra, la concentración de estas sustancias es altamente nociva para nuestra salud y por eso en este trabajo realizaremos el estudio de los siguientes contaminantes los cuales son los más preocupantes:
Ozono (O3)
Dióxido de azufre (SO2)
Particulado (PM10)
Tabla de datos de los contaminantes
aire <- read_excel("Concentracion.xlsx")
CM <- read_excel("Concentracion_Mov.xlsx")
datos <- as.data.frame(aire)
datatable(datos)
1. Se analiza la relación entre las variables
Lo primero que se tiene que hacer a la hora de establecer un modelo lineal múltiple es estudiar la relación que existe entre variables. Esta de suma importancia a la hora de identificar cuáles pueden ser los mejores predictores para el modelo, qué variables presentan relaciones de tipo no lineal y para identificar colinealidad entre predictores.
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.000
- Análisis con histogramas
multi.hist(x = aire, dcol = c("blue", "red"), dlty = c("dotted", "solid"))
- Ahora representando esto mismo, utilizando ggplot y ggally
ggpairs(datos, lower = list(continuous ="smooth"),
diag = list (continuos = "barDiag"), axisLabels = "none")
Entonces, de los análisis realizados hasta el momento, podemos obtener las siguientes conclusiones preliminares:
SO2 y PM10 están medianamente relacionados (r = 0.612) por lo que posiblemente no sea útil introducir ambos predictores en el modelo.
Las variables que tienen una mayor relación lineal con el O3 son: SO2 (r= -0.22), PM10 (r= -0.229) y Reactivación comercial (r= -0.193).
2. Generar el primer modelo de regresión lineal simple basado en la ecuación de la recta de mínimos cuadrados para una relación de X y Y
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-07
Bajo esta misma premisa el modelo tomaría forma de la siguiente manera:
\[ y = 29.04971 - 0.12090x \]
Evaluamos el modelo
y= 29.04971 - (0.12090*42.335417)
y
## [1] 23.93136
Se conoce el residuo segun este modelo
y- 24.91250
## [1] -0.9811419
Evaluamos gráficamente el modelo
plot(datos$PM10, datos$O3)
abline(modelo)
En este caso el modelo se pasa por 0.9811419 al llegar a el valor de 24.91250, por lo tanto esta sobrestimado.
Se utilizan multiples predictores
Aquí se implementa el método mixto iniciando el modelo con todas las variables como predictores y realizando la selección de los mejores predictores con la medición Akaike (AIC).
modelo2<-lm(O3 ~ SO2 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Estaciones_Transito + Lugares_Trabajo + Residencia, data = datos )
summary(modelo2)
##
## Call:
## lm(formula = O3 ~ SO2 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Estaciones_Transito + Lugares_Trabajo +
## Residencia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.9388 -4.4283 -0.0429 4.8415 23.5113
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.18847 1.46454 14.468 < 2e-16 ***
## SO2 -12.61924 2.15337 -5.860 9.10e-09 ***
## PM10 -0.05151 0.02816 -1.829 0.0681 .
## Reactivacion_Comercial -0.13450 0.09628 -1.397 0.1631
## Supermercado_Farmacia 0.47829 0.07098 6.738 5.07e-11 ***
## Parques_Centros -0.27200 0.06205 -4.383 1.47e-05 ***
## Estaciones_Transito -0.11122 0.06335 -1.756 0.0799 .
## Lugares_Trabajo -0.27806 0.05014 -5.546 5.06e-08 ***
## Residencia -1.06194 0.20450 -5.193 3.18e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.406 on 438 degrees of freedom
## Multiple R-squared: 0.2296, Adjusted R-squared: 0.2156
## F-statistic: 16.32 on 8 and 438 DF, p-value: < 2.2e-16
Seleccionando los mejores predictores
step(object= modelo2, direction="both", trace= 1)
## Start: AIC=1798.93
## O3 ~ SO2 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Estaciones_Transito + Lugares_Trabajo +
## Residencia
##
## Df Sum of Sq RSS AIC
## - Reactivacion_Comercial 1 107.04 24129 1798.9
## <none> 24022 1798.9
## - Estaciones_Transito 1 169.03 24191 1800.1
## - PM10 1 183.50 24206 1800.3
## - Parques_Centros 1 1053.76 25076 1816.1
## - Residencia 1 1478.93 25501 1823.6
## - Lugares_Trabajo 1 1686.94 25709 1827.3
## - SO2 1 1883.52 25906 1830.7
## - Supermercado_Farmacia 1 2490.04 26512 1841.0
##
## Step: AIC=1798.92
## O3 ~ SO2 + PM10 + Supermercado_Farmacia + Parques_Centros + Estaciones_Transito +
## Lugares_Trabajo + Residencia
##
## Df Sum of Sq RSS AIC
## <none> 24129 1798.9
## + Reactivacion_Comercial 1 107.04 24022 1798.9
## - PM10 1 190.65 24320 1800.4
## - Estaciones_Transito 1 568.24 24698 1807.3
## - Residencia 1 1372.22 25502 1821.6
## - Lugares_Trabajo 1 1762.50 25892 1828.4
## - Parques_Centros 1 2051.53 26181 1833.4
## - SO2 1 2257.88 26387 1836.9
## - Supermercado_Farmacia 1 2391.87 26521 1839.2
##
## Call:
## lm(formula = O3 ~ SO2 + PM10 + Supermercado_Farmacia + Parques_Centros +
## Estaciones_Transito + Lugares_Trabajo + Residencia, data = datos)
##
## Coefficients:
## (Intercept) SO2 PM10
## 21.40785 -13.37428 -0.05249
## Supermercado_Farmacia Parques_Centros Estaciones_Transito
## 0.45202 -0.31897 -0.16388
## Lugares_Trabajo Residencia
## -0.28339 -0.98924
Según la llamada final de este proceso (AIC) tenemos que el mejor modelo es:
modelo3<-lm(O3 ~ SO2 + PM10 + Supermercado_Farmacia + Parques_Centros + Estaciones_Transito + Lugares_Trabajo + Residencia, data = datos )
summary(modelo3)
##
## Call:
## lm(formula = O3 ~ SO2 + PM10 + Supermercado_Farmacia + Parques_Centros +
## Estaciones_Transito + Lugares_Trabajo + Residencia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.2494 -4.3746 -0.1753 4.9050 24.1825
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.40785 1.45767 14.686 < 2e-16 ***
## SO2 -13.37428 2.08670 -6.409 3.78e-10 ***
## PM10 -0.05249 0.02818 -1.862 0.0632 .
## Supermercado_Farmacia 0.45202 0.06852 6.597 1.21e-10 ***
## Parques_Centros -0.31897 0.05221 -6.109 2.21e-09 ***
## Estaciones_Transito -0.16388 0.05097 -3.215 0.0014 **
## Lugares_Trabajo -0.28339 0.05005 -5.663 2.70e-08 ***
## Residencia -0.98924 0.19798 -4.997 8.44e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.414 on 439 degrees of freedom
## Multiple R-squared: 0.2262, Adjusted R-squared: 0.2139
## F-statistic: 18.33 on 7 and 439 DF, p-value: < 2.2e-16
Es recomendable mostrar el intervalo de confianza para cada uno de los coeficientes parciales de regresión:
confint(lm(formula= O3 ~ SO2 + PM10 + Supermercado_Farmacia + Parques_Centros + Estaciones_Transito + Lugares_Trabajo + Residencia, data=datos))
## 2.5 % 97.5 %
## (Intercept) 18.5429587 24.272734653
## SO2 -17.4754428 -9.273108691
## PM10 -0.1078856 0.002902609
## Supermercado_Farmacia 0.3173485 0.586692242
## Parques_Centros -0.4215821 -0.216357958
## Estaciones_Transito -0.2640599 -0.063709120
## Lugares_Trabajo -0.3817479 -0.185032682
## Residencia -1.3783591 -0.600127926
Aquí se puede observar la relacion que guardan los sectores que estudiamos con el O3, el intervalo de 97.5 nos indica que el valor se encuentra en un determinado rango de valores con un 97.5% de certeza.
Se validan las condiciones para la regresión múltiple lineal.
La condicion tiene distintas formas de validarse, en este caso se validará mediante los diagramas de dispersión entre la variable dependiente y entre cada uno de los predicctores, o tambien se puede validar con diagramas de dispersión entre cada uno de los predictores asi como los residuos del modelo. Si la relacion es linea, los residuos deben de distribuirse de forma aleatoria en torno a 0 con una variabilidad constante a lo largo del eje x, esta ultima opción de validación suele ser la mas indicada y usada por la mayoría, ya que este permite identificar posibles datos anormales o atipicos.
plot1 <- ggplot(data= datos, aes(SO2, 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(Supermercado_Farmacia, modelo$residuals)) + geom_point() + geom_smooth(color= "firebrick") + geom_hline(yintercept= 0) + theme_bw()
plot4 <- ggplot(data= datos, aes(Parques_Centros, modelo$residuals)) + geom_point() + geom_smooth(color= "firebrick") + geom_hline(yintercept= 0) + theme_bw()
plot5 <- ggplot(data= datos, aes(Estaciones_Transito, 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 cumple la linealidad para todos los predictores.
Distribución normal de los residuos:
qqnorm(modelo3$residuals)
qqline(modelo3$residuals)
- Prueba de shapiro wilk para determinar normalidad
shapiro.test(modelo$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.9947, p-value = 0.1267
Despues de haber analizado toda la información gráfica que pudimos deducir, se determinó que la normalidad de los datos y la linealidad para los predictores, esto es debido a que se tiene una variabilidad constante respecto a estos residuos dominantes, estos en uso debido a los valores que se representan de manera continua a lo largo del eje x, mostrando a su vez una variabilidad no atipica, sino que mayormente constante y sin cambios visibles.
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'
Probando el modelo de regresión múltiple
bptest(modelo)
##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 10.687, df = 1, p-value = 0.001079
Matriz de correlación entre las variables predictoras
corrplot(cor(dplyr::select(datos,Supermercado_Farmacia,Residencia,
Reactivacion_Comercial,Parques_Centros,Lugares_Trabajo)),
method = "number", tl.col = "black")
> Podemos encontrar que se puede encontrar una relacion entre predictores asi mismo residencia tiene correlacion con los predictores. por otro lado podemos notar que hay una relacion en un aspecto positivo entre la reactivasion comercial con las supermecado-farmacias y parques-centros.
vif(modelo3)
## SO2 PM10 Supermercado_Farmacia
## 2.837797 1.617476 4.586790
## Parques_Centros Estaciones_Transito Lugares_Trabajo
## 3.810335 7.096867 6.237850
## Residencia
## 11.680499
Las personas al momento de estar en su casa, mantienen un estado de sedentarismo donde no hay lugar para moverse en lugares como supermercados, plazas y parques. La reactivación comercial provoca el incremento en las actividades de los sitios anteriores.
Autocorrelación
dwt(modelo, alternative = "two.sided")
## lag Autocorrelation D-W Statistic p-value
## 1 0.7529511 0.4940213 0
## Alternative hypothesis: rho != 0
Identificación de los posibles valores atípicos o influyentes
datos$studentized_residual <- rstudent(modelo)
ggplot(data = datos, aes(x = predict(modelo), y = abs(studentized_residual))) +
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
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))
Como lo marcan los puntos rojos podemos observar que si existen valores atípicos
which(abs(datos$studentized_residual) > 3)
## [1] 441
Cuantificación de las influencias de los predictores
summary(influence.measures(modelo))
## Potentially influential observations of
## lm(formula = O3 ~ PM10, data = datos) :
##
## dfb.1_ dfb.PM10 dffit cov.r cook.d hat
## 3 0.05 -0.06 -0.07 1.02_* 0.00 0.01
## 5 0.04 -0.06 -0.06 1.02_* 0.00 0.01_*
## 27 -0.22 0.20 -0.22_* 1.00 0.02 0.01
## 33 -0.03 0.03 -0.03 1.02_* 0.00 0.01
## 56 0.08 -0.07 0.08 1.01_* 0.00 0.01
## 199 -0.25 0.21 -0.25_* 0.98_* 0.03 0.01
## 201 -0.24 0.20 -0.25_* 0.97_* 0.03 0.01
## 206 -0.14 0.10 -0.16 0.98_* 0.01 0.00
## 207 -0.16 0.13 -0.17 0.99_* 0.01 0.01
## 213 -0.11 0.07 -0.14 0.98_* 0.01 0.00
## 218 -0.10 0.05 -0.15 0.97_* 0.01 0.00
## 219 0.13 -0.17 -0.19 1.01 0.02 0.02_*
## 220 -0.10 0.05 -0.15 0.97_* 0.01 0.00
## 224 -0.03 -0.01 -0.11 0.98_* 0.01 0.00
## 225 -0.10 0.07 -0.13 0.98_* 0.01 0.00
## 227 -0.06 0.02 -0.11 0.98_* 0.01 0.00
## 228 0.24 -0.30 -0.31_* 1.02_* 0.05 0.03_*
## 229 -0.01 0.01 0.01 1.01_* 0.00 0.01
## 255 -0.05 0.06 0.07 1.02_* 0.00 0.02_*
## 256 -0.03 0.04 0.04 1.03_* 0.00 0.02_*
## 260 0.04 -0.05 -0.06 1.02_* 0.00 0.02_*
## 263 0.00 0.00 0.00 1.02_* 0.00 0.01
## 277 0.04 -0.04 -0.05 1.03_* 0.00 0.02_*
## 278 0.06 -0.08 -0.08 1.02_* 0.00 0.02_*
## 281 0.05 -0.06 -0.07 1.01_* 0.00 0.01
## 291 0.04 -0.06 -0.06 1.02_* 0.00 0.02_*
## 292 0.04 -0.05 -0.05 1.01_* 0.00 0.01
## 298 0.03 -0.03 -0.04 1.02_* 0.00 0.01
## 314 0.00 0.00 0.00 1.03_* 0.00 0.03_*
## 315 0.05 -0.07 -0.07 1.02_* 0.00 0.01_*
## 329 0.04 -0.06 -0.06 1.02_* 0.00 0.01
## 335 0.04 -0.06 -0.06 1.02_* 0.00 0.02_*
## 378 -0.05 0.06 0.06 1.04_* 0.00 0.04_*
## 379 0.05 -0.06 -0.06 1.03_* 0.00 0.03_*
## 397 -0.39 0.38 -0.39_* 1.01_* 0.08 0.03_*
## 415 -0.02 0.07 0.13 0.98_* 0.01 0.00
## 432 -0.05 0.10 0.15 0.98_* 0.01 0.00
## 433 0.06 -0.03 0.10 0.99_* 0.01 0.00
## 440 0.01 0.04 0.14 0.97_* 0.01 0.00
## 441 0.09 -0.03 0.15 0.96_* 0.01 0.00
## 442 0.00 0.05 0.13 0.98_* 0.01 0.00
## 443 -0.11 0.18 0.22_* 0.98_* 0.02 0.01