A6U1
Equipo 4
01/3/2022
Regresión líneal múltiple
library(pacman)
p_load("prettydoc","DT","xfun", "dplyr", "psych", "GGally", "ggplot2", "readr", "gridExtra", "rmdformats", "readxl", "plotly", "leaflet","TSstudio", "corrplot", "lmtest", "car")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.
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 nivel de la superficie, el ozono en altas concentraciones es un contaminante del aire que provoca efectos nocivos en la salud humana, las plantas y los animales. Tiene además una contribución al calentamiento global” (Dirección de Monitoreo Atmosférico, 2016).
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 los datos
Los datos de contaminantes atmosféricos provienen de la estación de calidad del aire de la ERNO del instituto de geología de la unam ubicado en Hermosillo. Podemos visitar el origen de los datos aqui: http://www.erno.geologia.unam.mx
Datos relacionados con el Ozono como un contaminante aquí: http://www.aire.cdmx.gob.mx/descargas/noticias/que-es-ozono/que-es-ozono.pdf
Datos
aire <- read_excel("Concentracion.xlsx")
CM <- read_excel("Concentracion_Mov.xlsx")
datos <- as.data.frame(aire)
datatable(datos)1. Analizar la relación entre variables
El primer paso a la hora de establecer un modelo lineal múltiple es estudiar la relación que existe entre variables. Esta información es crítica 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:
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).
SO2 y PM10 están medianamente relacionados (r = 0.612) por lo que posiblemente no sea útil introducir ambos predictores en el modelo.
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-07Bajo este planteamiento el modelo tomaría la siguiente forma:
\[ y = 29.04971 - 0.12090x \]
Evaluando el modelo
y= 29.04971 - (0.12090*42.335417)
y## [1] 23.93136Conocer el residuo según este modelo
y- 24.91250## [1] -0.9811419Evaluando 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.
Ahora, con múltiples predictores
A continuación, se implementará 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-16Seleccionando 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.98924Segú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-16Es 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.600127926Con la tabla anterior se puede observar la relación que guardan los sectores de estudio 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.
Validación de condiciones para la regresión múltiple lineal.
Esta condición se puede validar bien mediante diagramas de dispersión entre la variable dependiente y cada uno de los predictores (como se ha hecho en el análisis preliminar) o con diagramas de dispersión entre cada uno de los predictores y los residuos del modelo. Si la relación es lineal, los residuos deben de distribuirse aleatoriamente en torno a 0 con una variabilidad constante a lo largo del eje X. Esta última opción suele ser más indicada ya que permite identificar posibles datos atípicos.
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.1267Gracias al análisis gráfico se pudo determinar que se confirma la normalidad de los datos y la linealidad para los predictores, debido a que se tiene una variabilidad constante respecto a los residuos. Esto dandose debido a que los valores se representan de manera continua a lo largo del eje x, mostrando a la vez una varibilidad mayormente constante.
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'
bptest(modelo3)##
## studentized Breusch-Pagan test
##
## data: modelo3
## BP = 35.112, df = 7, p-value = 1.066e-05El test de Breusch-Pagan anterior nos muestra que el valor de p es menor a 0.05 por lo que hay evidencia de falta de homocedasticidad, por lo que decimos que la heterocedasticidad está presente en los residuos, por ende los resultados se vuelven poco fiables.
Matriz de correlación entre predictores para el modelo 3:
corrplot(cor(dplyr::select(datos, SO2, PM10, Supermercado_Farmacia, Parques_Centros, Estaciones_Transito, Lugares_Trabajo, Residencia)),
method = "number", tl.col = "black")
En la matriz se puede observar que residencia guarda una relación fuerte con otros predictores como Supermercado_Farmacia, Parques_Centros, Estaciones_Transito, Lugares_Trabajo.
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.680499Como se puede observar también la matriz el factor de residencia muestra una alta correlación línea e inflación de varianza.
dwt(modelo, alternative = "two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 0.7529511 0.4940213 0
## Alternative hypothesis: rho != 0Identificación de los posibles valores atípicos 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") + 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))
Existen observaciones muy atípicas
which(abs(datos$studentized_residual)>3)## [1] 12Cuantifiació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