Regresión Lineal Múltiple
Lic. Miguel Angel
2024-03-08
Ejemplo de Regresión Lineal Múltiple
Un estudio quiere generar un modelo que permita predecir la esperanza de vida media de los habitantes de una ciudad en función de diferentes variables. Se dispone de información sobre: habitantes, analfabetismo, ingresos, esperanza de vida, asesinatos, universitarios, heladas, área y densidad poblacional.
#install.packages("psych")
#install.packages("GGally")
#install.packages("corrplot")
#install.packages("car")
library(readxl)
library(psych)
library(GGally)
library(dplyr)
library(ggplot2)
library(corrplot)
# Para facilitar su interpretación se renombra y se modifica
datos <- as.data.frame(state.x77)
datos <- rename(habitantes = Population, analfabetismo = Illiteracy,
ingresos = Income, esp_vida = `Life Exp`, asesinatos = Murder,
universitarios = `HS Grad`, heladas = Frost, area = Area,
.data = datos)
datos <- mutate(.data = datos, densidad_pobl = habitantes * 1000 / area)1. Analizar la relación entre las 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 (por lo que no pueden ser incluidas) y para identificar colinealidad entre predictores. A modo complementario, es recomendable representar la distribución de cada variable mediante histogramas.
Las dos formas principales de hacerlo son mediante representaciones gráficas (gráficos de dispersión) y el cálculo del coeficiente de correlación de cada par de variables.
cor_matriz <- round(cor(x = datos, method = "pearson"), 3)
knitr::kable(cor_matriz, caption = "Matriz de Correlación de Pearson")| habitantes | ingresos | analfabetismo | esp_vida | asesinatos | universitarios | heladas | area | densidad_pobl | |
|---|---|---|---|---|---|---|---|---|---|
| habitantes | 1.000 | 0.208 | 0.108 | -0.068 | 0.344 | -0.098 | -0.332 | 0.023 | 0.246 |
| ingresos | 0.208 | 1.000 | -0.437 | 0.340 | -0.230 | 0.620 | 0.226 | 0.363 | 0.330 |
| analfabetismo | 0.108 | -0.437 | 1.000 | -0.588 | 0.703 | -0.657 | -0.672 | 0.077 | 0.009 |
| esp_vida | -0.068 | 0.340 | -0.588 | 1.000 | -0.781 | 0.582 | 0.262 | -0.107 | 0.091 |
| asesinatos | 0.344 | -0.230 | 0.703 | -0.781 | 1.000 | -0.488 | -0.539 | 0.228 | -0.185 |
| universitarios | -0.098 | 0.620 | -0.657 | 0.582 | -0.488 | 1.000 | 0.367 | 0.334 | -0.088 |
| heladas | -0.332 | 0.226 | -0.672 | 0.262 | -0.539 | 0.367 | 1.000 | 0.059 | 0.002 |
| area | 0.023 | 0.363 | 0.077 | -0.107 | 0.228 | 0.334 | 0.059 | 1.000 | -0.341 |
| densidad_pobl | 0.246 | 0.330 | 0.009 | 0.091 | -0.185 | -0.088 | 0.002 | -0.341 | 1.000 |
Otros paquetes permiten representar a la vez los diagramas de dispersión, los valores de correlación para cada par de variables y la distribución de cada una de las variables.
ggpairs(datos, lower = list(continuous = "smooth"),
diag = list(continuous = "barDiag"), axisLabels = "none")
Las variables que tienen una mayor relación lineal con la esperanza de vida son: asesinatos (r= -0.78), analfabetismo (r= -0.59) y universitarios (r= 0.58).
Asesinatos y analfabetismo están medianamente correlacionados (r = 0.7) por lo que posiblemente no sea útil introducir ambos predictores en el modelo.
Las variables habitantes, área y densidad poblacional muestran una distribución exponencial, una transformación logarítmica posiblemente haría más normal su distribución.
2.Generar el modelo
Se va a emplear 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).
modelo <- lm(esp_vida ~ habitantes + ingresos + analfabetismo + asesinatos +
universitarios + heladas + area + densidad_pobl, data = datos )
summary(modelo)##
## Call:
## lm(formula = esp_vida ~ habitantes + ingresos + analfabetismo +
## asesinatos + universitarios + heladas + area + densidad_pobl,
## data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47514 -0.45887 -0.06352 0.59362 1.21823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.995e+01 1.843e+00 37.956 < 2e-16 ***
## habitantes 6.480e-05 3.001e-05 2.159 0.0367 *
## ingresos 2.701e-04 3.087e-04 0.875 0.3867
## analfabetismo 3.029e-01 4.024e-01 0.753 0.4559
## asesinatos -3.286e-01 4.941e-02 -6.652 5.12e-08 ***
## universitarios 4.291e-02 2.332e-02 1.840 0.0730 .
## heladas -4.580e-03 3.189e-03 -1.436 0.1585
## area -1.558e-06 1.914e-06 -0.814 0.4205
## densidad_pobl -1.105e-03 7.312e-04 -1.511 0.1385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7337 on 41 degrees of freedom
## Multiple R-squared: 0.7501, Adjusted R-squared: 0.7013
## F-statistic: 15.38 on 8 and 41 DF, p-value: 3.787e-10El modelo con todas las variables introducidas como predictores tiene un R^2 alta (0.7501), es capaz de explicar el 75,01% de la variabilidad observada en la esperanza de vida. El p-value del modelo es significativo (3.787e-10) por lo que se puede aceptar que el modelo no es por azar, al menos uno de los coeficientes parciales de regresión es distinto de 0. Muchos de ellos no son significativos, lo que es un indicativo de que podrían no contribuir al modelo.
3.Selección de los mejores predictores
En este caso se van a emplear la estrategia de stepwise mixto. El valor matemático empleado para determinar la calidad del modelo va a ser Akaike(AIC).
step(object = modelo, direction = "both", trace = 1)## Start: AIC=-22.89
## esp_vida ~ habitantes + ingresos + analfabetismo + asesinatos +
## universitarios + heladas + area + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - analfabetismo 1 0.3050 22.373 -24.208
## - area 1 0.3564 22.425 -24.093
## - ingresos 1 0.4120 22.480 -23.969
## <none> 22.068 -22.894
## - heladas 1 1.1102 23.178 -22.440
## - densidad_pobl 1 1.2288 23.297 -22.185
## - universitarios 1 1.8225 23.891 -20.926
## - habitantes 1 2.5095 24.578 -19.509
## - asesinatos 1 23.8173 45.886 11.707
##
## Step: AIC=-24.21
## esp_vida ~ habitantes + ingresos + asesinatos + universitarios +
## heladas + area + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - area 1 0.1427 22.516 -25.890
## - ingresos 1 0.2316 22.605 -25.693
## <none> 22.373 -24.208
## - densidad_pobl 1 0.9286 23.302 -24.174
## - universitarios 1 1.5218 23.895 -22.918
## + analfabetismo 1 0.3050 22.068 -22.894
## - habitantes 1 2.2047 24.578 -21.509
## - heladas 1 3.1324 25.506 -19.656
## - asesinatos 1 26.7071 49.080 13.072
##
## Step: AIC=-25.89
## esp_vida ~ habitantes + ingresos + asesinatos + universitarios +
## heladas + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - ingresos 1 0.132 22.648 -27.598
## - densidad_pobl 1 0.786 23.302 -26.174
## <none> 22.516 -25.890
## - universitarios 1 1.424 23.940 -24.824
## + area 1 0.143 22.373 -24.208
## + analfabetismo 1 0.091 22.425 -24.093
## - habitantes 1 2.332 24.848 -22.962
## - heladas 1 3.304 25.820 -21.043
## - asesinatos 1 32.779 55.295 17.033
##
## Step: AIC=-27.6
## esp_vida ~ habitantes + asesinatos + universitarios + heladas +
## densidad_pobl
##
## Df Sum of Sq RSS AIC
## - densidad_pobl 1 0.660 23.308 -28.161
## <none> 22.648 -27.598
## + ingresos 1 0.132 22.516 -25.890
## + analfabetismo 1 0.061 22.587 -25.732
## + area 1 0.043 22.605 -25.693
## - habitantes 1 2.659 25.307 -24.046
## - heladas 1 3.179 25.827 -23.030
## - universitarios 1 3.966 26.614 -21.529
## - asesinatos 1 33.626 56.274 15.910
##
## Step: AIC=-28.16
## esp_vida ~ habitantes + asesinatos + universitarios + heladas
##
## Df Sum of Sq RSS AIC
## <none> 23.308 -28.161
## + densidad_pobl 1 0.660 22.648 -27.598
## + ingresos 1 0.006 23.302 -26.174
## + analfabetismo 1 0.004 23.304 -26.170
## + area 1 0.001 23.307 -26.163
## - habitantes 1 2.064 25.372 -25.920
## - heladas 1 3.122 26.430 -23.877
## - universitarios 1 5.112 28.420 -20.246
## - asesinatos 1 34.816 58.124 15.528##
## Call:
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
## heladas, data = datos)
##
## Coefficients:
## (Intercept) habitantes asesinatos universitarios heladas
## 7.103e+01 5.014e-05 -3.001e-01 4.658e-02 -5.943e-03El mejor modelo resultante del proceso de selección ha sido:
modelo <- (lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
heladas, data = datos))
summary(modelo)##
## Call:
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
## heladas, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## habitantes 5.014e-05 2.512e-05 1.996 0.05201 .
## asesinatos -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## universitarios 4.658e-02 1.483e-02 3.142 0.00297 **
## heladas -5.943e-03 2.421e-03 -2.455 0.01802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12Es recomendable mostrar el intervalo de confianza para cada uno de los coeficientes parciales de regresión:
confint(lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
heladas, data = datos))## 2.5 % 97.5 %
## (Intercept) 6.910798e+01 72.9462729104
## habitantes -4.543308e-07 0.0001007343
## asesinatos -3.738840e-01 -0.2264135705
## universitarios 1.671901e-02 0.0764454870
## heladas -1.081918e-02 -0.0010673977Cada una de las pendientes de un modelo de regresión lineal múltiple (coeficientes parciales de regresión de los predictores) se define del siguiente modo: Si el resto de variables se mantienen constantes, por cada unidad que aumenta el predictor en cuestión, la variable (Y) varía en promedio tantas unidades como indica la pendiente. Para este ejemplo, por cada unidad que aumenta el predictor universitarios, la esperanza de vida aumenta en promedio 0.04658 unidades, manteniéndose constantes el resto de predictores.
4.Validación de condiciones para la regresión múltiple lineal
Relación lineal entre los predictores numéricos y la variable respuesta:
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.
library(ggplot2)
library(gridExtra)
plot1 <- ggplot(data = datos, aes(habitantes, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos, aes(asesinatos, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = datos, aes(universitarios, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot4 <- ggplot(data = datos, aes(heladas, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2, plot3, plot4)
Se cumple la linealidad para todos los predictores
Distribución normal de los residuos:
qqnorm(modelo$residuals)
qqline(modelo$residuals)
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.97935, p-value = 0.525Tanto 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.
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'
library(lmtest)
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 6.2721, df = 4, p-value = 0.1797Con un valor p de 0.1797, no hay suficiente evidencia para rechazar la hipótesis nula de homocedasticidad. Esto significa que no hay indicios estadísticos significativos de heterocedasticidad en los errores del modelo. En otras palabras, no hay pruebas sólidas para afirmar que la varianza de los errores varía con los valores de las variables independientes, y por lo tanto, el supuesto de homocedasticidad no se viola significativamente.
En la práctica, esto indica que, el modelo de regresión satisface el supuesto de homocedasticidad y los errores tienen varianzas aproximadamente constantes. Esto es favorable para la fiabilidad de las estimaciones estándar y las inferencias estadísticas que se derivan del modelo.
No multicolinialidad:
Matriz de correlación entre predictores.
corrplot(cor(dplyr::select(datos, habitantes, asesinatos,universitarios,heladas)),
method = "number", tl.col = "black")
Los valores positivos están en tonos de rojo y los negativos en tonos de azul. La intensidad del color indica la fuerza de la correlación. Las correlaciones más fuertes en esta matriz son negativas y están entre asesinatos y universitarios, y asesinatos y heladas. Las correlaciones positivas más fuertes son entre habitantes y asesinatos, y universitarios y heladas.
No hay correlaciones extremadamente fuertes (cercanas a 1 o -1), lo que indica que ninguna de las variables predictoras está perfectamente predicha por otra. Recordar que la correlación no implica causalidad, y estos resultados solo muestran relaciones lineales.
Análisis de Inflación de Varianza (VIF):
library(car)## Warning: package 'car' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:psych':
##
## logitvif(modelo)## habitantes asesinatos universitarios heladas
## 1.189835 1.727844 1.356791 1.498077El VIF es una medida que cuantifica la intensidad de la multicolinealidad en un modelo de regresión ordinaria.La multicolinealidad ocurre cuando dos o más variables predictoras en el modelo están correlacionadas, y esto puede llevar a una estimación ineficiente de los coeficientes de regresión.
VIF = 1: No hay correlación entre la variable independiente en cuestión y las otras variables. 1 < VIF < 5: Generalmente, se considera un nivel moderado de multicolinealidad. VIF >= 5: Señala una multicolinealidad que puede ser problemática y que podría afectar la precisión de las estimaciones de los coeficientes de la regresión.
Todos los valores de VIF están por debajo de 5, lo cual indica que no hay evidencia de una multicolinealidad excesiva que afecte de manera significativa las estimaciones de los coeficientes.
Autocorrelación:
library(car)
dwt(modelo, alternative = "two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 0.02867262 1.913997 0.736
## Alternative hypothesis: rho != 0El valor p asociado con la prueba de Durbin-Watson. Un valor p alto (generalmente más grande que 0.05) indica que no hay suficiente evidencia para rechazar la hipótesis nula de que no hay autocorrelación (rho=0). Aquí, el valor p es 0.8, lo que es mucho mayor que 0.05, por lo que no rechazamos la hipótesis nula de ausencia de autocorrelación.
Parece que no hay autocorrelación significativa en los residuos del modelo de regresión al nivel de retraso 1. La estadística de Durbin-Watson cercana a 2 y el valor p alto apoyan la suposición de que los residuos son independientes.
5.Identificación de posibles valores atípicos o influyentes
library(dplyr)
datos$studentized_residual <- rstudent(modelo)
ggplot(data = datos, aes(x = predict(modelo), 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))
which(abs(datos$studentized_residual) > 3)## integer(0)El resultado significa que no hay ningún residuo studentizado en el conjunto de datos que tenga un valor absoluto mayor que 3. Esto es significativo porque los valores absolutos de los residuos studentizados que son mayores que 3 son comúnmente considerados como indicativos de valores atípicos potenciales en el modelo de regresión. En términos estadísticos, esto implica que todos los residuos están dentro de un rango que se consideraría “normal” si los residuos se distribuyen aproximadamente de forma normal.No se identifica ninguna observación atípica.
El Modelo Lineal Múltiple:
\[ Esperanza\ de\ vida = 5.014e^{-05} \cdot habitantes - 3.001e^{-01} \cdot asesinatos + 4.658e^{-02} \cdot universitarios - 5.943e^{-03} \cdot heladas \]
es capaz de explicar el 73.6% de la variabilidad observada en la esperanza de vida (R^2: 0.736, R^2-Adjusted: 0.7126). El test F muestra que es significativo (p-value: 1.696e-12). En este caso Se satisfacen todas las condiciones para este tipo de regresión múltiple.