Regresión lineal múltiple
library(pacman)
p_load("prettydoc", "DT", "xfun", "dplyr", "psych", "GGally", "ggplot2","readr")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.
Ejemplo 1. Predictores numericos para la esperanza de vida en ciudades de estados unidos
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 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
Datos
datos <- as.data.frame(state.x77)
datatable(datos)Renombrar y modificar el conjunto de datos
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 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 colinialidad entre predictores. A modo complementario, es recomendable representar la distribución de cada variable mediante histogramas.
round(cor(x = datos, method = "pearson"), 3)## habitantes ingresos analfabetismo esp_vida asesinatos
## habitantes 1.000 0.208 0.108 -0.068 0.344
## ingresos 0.208 1.000 -0.437 0.340 -0.230
## analfabetismo 0.108 -0.437 1.000 -0.588 0.703
## esp_vida -0.068 0.340 -0.588 1.000 -0.781
## asesinatos 0.344 -0.230 0.703 -0.781 1.000
## universitarios -0.098 0.620 -0.657 0.582 -0.488
## heladas -0.332 0.226 -0.672 0.262 -0.539
## area 0.023 0.363 0.077 -0.107 0.228
## densidad_pobl 0.246 0.330 0.009 0.091 -0.185
## universitarios heladas area densidad_pobl
## habitantes -0.098 -0.332 0.023 0.246
## ingresos 0.620 0.226 0.363 0.330
## analfabetismo -0.657 -0.672 0.077 0.009
## esp_vida 0.582 0.262 -0.107 0.091
## asesinatos -0.488 -0.539 0.228 -0.185
## universitarios 1.000 0.367 0.334 -0.088
## heladas 0.367 1.000 0.059 0.002
## area 0.334 0.059 1.000 -0.341
## densidad_pobl -0.088 0.002 -0.341 1.000- Análisis con histogramas
multi.hist( x = datos, dcol = c("blue","red"), dlty = c("dotted", "solid"),
main = ""
)
Ahora representando esto mismo, más la recta 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 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
Como se ha explicado en la introducción, hay diferentes formas de llegar al modelo final más adecuado. En este caso 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 ~ asesinatos, data=datos)
summary(datos)## habitantes ingresos analfabetismo esp_vida
## Min. : 365 Min. :3098 Min. :0.500 Min. :67.96
## 1st Qu.: 1080 1st Qu.:3993 1st Qu.:0.625 1st Qu.:70.12
## Median : 2838 Median :4519 Median :0.950 Median :70.67
## Mean : 4246 Mean :4436 Mean :1.170 Mean :70.88
## 3rd Qu.: 4968 3rd Qu.:4814 3rd Qu.:1.575 3rd Qu.:71.89
## Max. :21198 Max. :6315 Max. :2.800 Max. :73.60
## asesinatos universitarios heladas area
## Min. : 1.400 Min. :37.80 Min. : 0.00 Min. : 1049
## 1st Qu.: 4.350 1st Qu.:48.05 1st Qu.: 66.25 1st Qu.: 36985
## Median : 6.850 Median :53.25 Median :114.50 Median : 54277
## Mean : 7.378 Mean :53.11 Mean :104.46 Mean : 70736
## 3rd Qu.:10.675 3rd Qu.:59.15 3rd Qu.:139.75 3rd Qu.: 81163
## Max. :15.100 Max. :67.30 Max. :188.00 Max. :566432
## densidad_pobl
## Min. : 0.6444
## 1st Qu.: 25.3352
## Median : 73.0154
## Mean :149.2245
## 3rd Qu.:144.2828
## Max. :975.0033Bajo este planteamiento el modelo tomaría la siguiente forma:
\[ y = 72.97356 - 0.28395x \]
- Evaluando el modelo
y = 72.97356 - (0.28395*11.6)
y## [1] 69.67974- Conocer el residuo según este modelo
y - 67.96## [1] 1.71974- Evaluando gráficamente el modelo
plot(datos$asesinatos, datos$esp_vida)
abline(modelo)
Ahora, con múltiples predictores
modelo2 <- (lm(formula = esp_vida ~ habitantes + asesinatos + universitarios + heladas, data = datos))
summary(modelo2) ##
## 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
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.