Regresión múltiple

require(MASS) # importa la libreria MASS con funciones y métodos estadísticos
require(ISLR) # importa la librería ISLR con datos para un aprender estadística con R
data("Boston") # utiliza el conjunto de datos Boston del paquete MASS. Mediana del valor de la vivienda

# Estudio simple de los datos
require(psych) # importa la librería psych con procedimientos para la investigación psicológica.
summary(Boston) # muestra estadísticos básicos del conjunto de datos
multi.hist(x = Boston[,1:3], dcol = c("blue","red"), dlty = c("dotted", "solid"), main = "") # crea un histograma básico

# Creacion de un modelo múltiple de regresion para predecir el valor de la vivienda (medv) en función del resto de datos
modelo_multiple <- lm(formula = medv ~ ., data = Boston)
# También se pueden especificar una a una 
summary(modelo_multiple) # para visualizar los principales parámetros del modelo generado

Datos a observar:

  • p-value del estadístico F. Si es bajo (<0.05) significa que hay al menos un predictor que tiene relación con la respuesta.
  • RSE (Residual standard error): Son las unidades que se aleja, en promedio, cualquier predicción del modelo. Hay que tener en cuenta el promedio de la variable respuesta, por lo que el RSE se suele medir en porcentaje: \frac{RSE}{promedio variable respuesta}=X\%
  • R^2: Indica el tanto por uno (multiplicar por 100 para saber el %) de la variabilidad explicada por el modelo. La ventaja de R^2 es que es independiente de la escala en la que se mida la variable respuesta, por lo que su interpretación es más sencilla.
  • Intercept(\beta_0): El valor promedio de la variable respuesta cuando todos los predictores son 0. Es la ordenada en el origen en un entorno n-dimensional
  • Predictores(\beta_i): Por cada unidad que se incrementa cada predictor, indica cuanto vale la variable respuesta. Es la pendiente en un entorno n-dimensional. En esta tabla se muestra también los p-values, por lo que puede verse para cada variable si influye o no en el model: p-values muy cercanos a 1 significa que la variable no influye en el modelo y deberian excluirse. Esta exclusión basada en los p-value no es recomendable, en su lugar se recomienda emplear métodos de best subset selection, stepwise selection (forward, backward e hybrid) o Shrinkage/regularization.

Para usar el método de Stepwise existe la funcion “step”: step(object, scope, scale = 0, direction = c(“both”, “backward”, “forward”), trace = 1, keep = NULL, steps = 1000, k = 2, …)

step(modelo_multiple, direction = "both", trace = 0)

La función nos muestra los valores Intercept y de los multiplicadores de cada variable para aquellas que tienen influencia en el odelo, descartando otras que tienen menor influencia. Una vez descartadas algunas de las variables, se ejecuta de nuevo el modelo para ver sus estadísticos:

modelo_multiple <- lm(formula = medv ~ crim + zn + chas +  nox + rm +  dis + rad + tax + ptratio + black + lstat, data = Boston)
# También se pueden indicar todas las variables de un data.frame y exluir algunas
# modelo_multiple <- lm(formula = medv~. -age -indus, data = Boston)
summary(modelo_multiple)

Estudio de los residuos del modelo.

Para confirmar que un modelo de regresión lineal múltiple por mínimos cuadrados cumple su labor, lo mejor es estudiar los residuos del modelo. En R, los residuos se almacenan dentro del modelo, bajo el nombre de “residuals”, además, R puede generar los gráficos más usados para la evaluar los residuos del modelo:

par(mfrow = c(1,2)) # Indicamos que queremos 2 gráficos por pantalla (1 fila y 2 columnas)
plot(modelo_multiple) # Pintamos los graficos de residuos

Aparecen 4 gráficos:

  1. El primer gráfico muestra los residuos para cada predicción. Si hay una distribución lineal de los residuos, entonces podemos pensar que el modelo se ajusta bien. Si se observa cualquier tipo de distribución no lineal ni constante, podemos sospechar que el modelo no se ajusta del todo.
  2. El segundo gráfico es de tipo Normal Q-Q. Compara los residuos estandarizados con la distribución teorica Normal. Si los puntos se alejan de la linea normal, indica que los residuos no se ajustan a una distribución Normal.
  3. El tercer gráfico, Scale-location o también llamado Spread-Location, muestra si los residuos se distribuyen por igual a lo largo de los rangos de predictores. Permite verificar si la varianza es igual (homocedasticidad). Si se ve una línea horizontal con puntos de dispersión iguales (al azar), es una buena señal.
  4. EL último de los gráficos es el Residuals Vs Leverage. Nos ayuda a encontrar casos influyentes. No todos los valores atípicos influyen en el análisis de regresión lineal. Aunque parezcan valores extremos por salirse de lo habitual, igual no influyen en el modelo. De igual forma, valores que aparentemente confirman el modelo, podrían ser muy influyentes y alterar los resultados si los excluimos del análisis. En este gráfico los patrones no son relevantes. Hay que observar valores periféricos en la esquina superior derecha o en la esquina inferior derecha. Esos lugares son los lugares donde los casos pueden influir en una línea de regresión. Hay que buscar casos fuera de la línea, la distancia de Cook. Los resultados de la regresión se alterarán si excluimos esos casos.

Otra forma de identificar las observaciones que puedan ser outliers o puntos con alta influencia (leverage) es emplear las funciones rstudent() y hatvalues().
La función “rstudent” calcula los residuos Studentizados, que vuelve a normalizar los residuos para tener una varianza unitaria. Típicamente, las desviaciones estándar de los residuos en una muestra varían mucho de un punto de datos a otro, incluso cuando todos los errores tienen la misma desviación estándar, particularmente en el análisis de regresión; por lo tanto, no tiene sentido comparar los residuos en diferentes puntos de datos sin primero Studentizarlos. Se considera que los casos que superan un nivel de 3 son outlayers.

plot(x = modelo_multiple$fitted.values, y = abs(rstudent(modelo_multiple)), main = "Residuos absolutes studentizados vs valores predichos", pch = 20, col = "grey30")
abline(h = 3, col = "red") # se pinta la linea límite de 3

Las observaciones con alto leverage, son observaciones altamente influyentes, puesto que podría estar condicionando el modelo. La eliminación de este tipo de observaciones debe de analizarse con detalle y dependiendo de la finalidad del modelo. Si el fin es predictivo, un modelo sin estas observaciones puede lograr mayor precisión la mayoría de casos. Sin embargo, es muy importante prestar atención a estos valores ya que, de no ser errores de medida, pueden ser los casos más interesantes. El modo adecuado a proceder cuando se sospecha de algún posible valor atípico o influyente es calcular el modelo de regresión incluyendo y excluyendo dicho valor. Observaciones cuyos valores superen 2.5x((p+1)/n), siendo p el número de predictores y n el número de observaciones, deben tenerse en consideración por su posible influencia.

plot(hatvalues(modelo_multiple), main = "Medición de leverage", pch = 20)
# Se añade una línea en el threshold de influencia acorde a la regla 2.5x((p+1)/n)
abline(h = 2.5*((dim(modelo_multiple$model)[2]-1 + 1)/dim(modelo_multiple$model)[1]), col = "red")

En los modelos de regresión lineal con múltiples predictores, además del estudio de los residuos vistos, es necesario descartar colinealidad o multicolinealidad entre variables.

Colinealidad

Para la colinealidad se recomienda calcular el coeficiente de correlación entre cada par de predictores incluidos en el modelo:

require(corrplot) # importa la libreria para generar graficos de correlación
corrplot.mixed(corr = cor(Boston[,c("crim", "zn", "nox", "rm", "dis", "rad", "tax", "ptratio", "black", "lstat", "medv")], method = "pearson")) # grafico con los coeficientes de correlación

El gráfico muestra arriba a la derecha una forma visual de correlación, marcando en azul la correlación positiva y en rojo la negativa. En la parte izquierda/abajo muestra los valores numéricos. Una vez vistas las correlaciones más fuertes, se pueden ver particularmente:

attach(Boston)
par(mfrow = c(2,2))
plot(x = tax, y = rad, pch = 20)
plot(x = tax, y = nox, pch = 20)
plot(x = dis, y = nox, pch = 20)
plot(x = medv, y = rm, pch = 20)

Si la correlación es alta y por lo tanto las variables aportan información redundante, es recomendable analizar si el modelo mejora o no empeora excluyendo alguno de estos predictores.

Multicolinealidad

Para el estudio de la multicolinealidad una de las medidas más utilizadas es el factor de inflación de varianza VIF. Puede calcularse mediante la función vif() del paquete car.

require(car) # importa el paquete car
vif(modelo_multiple) # calcula el factor de inflación de varianza 

Si lo indices VIF son bajos o moderados, valores entre 5 y 10 indican posibles problemas y valores mayores o iguales a 10 se consideran muy problemáticos.

Interacción entre predictores

Una de las asunciones del modelo de regresión lineal múltiple es la de aditividad, según la cual los efectos que causan sobre la variable respuesta Y variaciones en el predictor X_i son independientes del valor que tomen los otros predictores. Se conoce como efecto de interacción cuando el efecto de un predictor varía dependiendo del valor que adquiera otro predictor. Si esto ocurre, el modelo mejorará si se incluye dicha interacción. Si en un modelo se incorpora una interacción, se deben incluir también los predictores individuales que forman la interacción aun cuando por sí solos no sean significativos.

Para simplificar, podemos generar un modelo lineal con dos variables predictoras. El modelo lineal múltiple empleando ambos predictores resulta en:

modelo <- lm(medv  ~ lstat + age, data = Boston) # genera el nuevo modelo
summary(modelo) # muestra los valores del modelo

Como es un modelo con dos predictores continuos se puede representar el plano de regresión:

rango_lstat <- range(Boston$lstat)
nuevos_valores_lstat <- seq(from = rango_lstat[1], to = rango_lstat[2], length.out = 20)
rango_age <- range(Boston$age)
nuevos_valores_age <- seq(from = rango_age[1], to = rango_age[2], length.out = 20)

predicciones <- outer(X = nuevos_valores_lstat, Y = nuevos_valores_age, FUN = function(lstat, age) {
                          predict(object = modelo, newdata = data.frame(lstat, age))
                          })

superficie <- persp(x = nuevos_valores_lstat, y = nuevos_valores_age, z = predicciones, theta = 20, phi = 5,
                    col = "lightblue", shade = 0.1, zlim = range(-10,100), xlab = "lstat", ylab = "age", zlab = "medv",
                    ticktype = "detailed", main = "Predición precio medio ~ lstat y age")

observaciones <- trans3d(Boston$lstat, Boston$age, Boston$medv, superficie)
error <- trans3d(Boston$lstat, Boston$age, fitted(modelo), superficie)
points(observaciones, col = "red", pch = 16)
segments(observaciones$x, observaciones$y, error$x, error$y)

En R se puede generar un modelo con interacción de dos formas: indicando de forma explícita los predictores individuales y entre que predictores se quiere evaluar la interacción

lm(medv ~ lstat + age + lstat:age, data = Boston)

O bien de forma directa

modelo_interaccion  <- lm(medv ~ lstat * age, data = Boston)
summary(modelo_interaccion)

Si el R-squared mejora introduciendo la ineracción, entonces el modelo más adecuado es con la interacción.