Regresión lineal en R

Regresión simple

[sourcecode language=»r» wraplines=»false» collapse=»false»]
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 simple de regresion para predecir el valor de la vivienda (medv) en función del porcentaje de pobreza de la población (lstat)
modelo_simple = lm(data = Boston,formula = medv ~ lstat) # genera un modelo de regresión lineal por mínimos cuadrados en el que la variable respuesta es medv y el predictor lstat.
names(modelo_simple) # para ver el contenido de modelo generado
summary(modelo_simple) # para visualizar los principales parámetros del modelo generado
[/sourcecode]

Datos a observar:

• p-value del estadístico F. Si es bajo (<0.05) significa que hay una relación entre el predictor y 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 el predictor es 0. Es la ordenada en el origen
• Predictor($\beta_1$): Por cada unidad que se incrementa el predictor, indica cuanto vale la variable respuesta. Es la pendiente.

Con los coeficientes de regresión ya pueden hacerse estimaciones, pero esos coeficiente no dejan de ser una estimación con un error estándar y un intervalo de confianza. Para calcularlo se utiliza la función:

[sourcecode language=»r» wraplines=»false» collapse=»false»]
confint(modelo_simple, level = 0.95)
[/sourcecode]

Nos indicará para la ordenada en el origen y la pendiente, los valores extremos al 2.5% y 97.5%

Estimación de datos con el modelo

Para realizar predicciones hay que tener en cuenta que tienen asociado un error y por lo tanto un intervalo. Es importante diferenciar entre dos tipos de intervalo:

• Intervalo de confianza: Devuelve un intervalo para el valor promedio de la variable a predecir con un determinado valor del predictor, supóngase 10.

  [sourcecode language=»r» wraplines=»false» collapse=»false»]
  predict(object = modelo_simple, newdata = data.frame(lstat = c(10)), interval = "confidence", level = 0.95)
  [/sourcecode]

• Intervalo de predicción: Devuelve un intervalo para el valor esperado de un valor a predecir con un determinado valor del predictor.

  [sourcecode language=»r» wraplines=»false» collapse=»false»]
  predict(object = modelo_simple, newdata = data.frame(lstat = c(10)), interval = "prediction", level = 0.95)
  [/sourcecode]

Ambos intervalos estarán centrados en torno al mismo valor. Y, aunque parezcan similares, la diferencia se encuentra en que los intervalos de confianza se aplican al valor promedio que se espera de y para un determinado valor de x, mientras que los intervalos de predicción no se aplican al promedio. Por esta razón, los segundos siempre son más amplios que los primeros.

Representación gráfica

La creación de un modelo de regresión lineal simple suele acompañarse de una representación gráfica superponiendo las observaciones con el modelo. Además de ayudar a la interpretación, es el primer paso para identificar posibles violaciones de las condiciones de la regresión lineal.

[sourcecode language=»r» wraplines=»false» collapse=»false»]
attach(Boston)
plot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30") # Pinta las observaciones
abline(modelo_simple, lwd = 3, col = "red") # Pinta la linea de regresión
[/sourcecode]

Una de las mejores formas de confirmar que las condiciones necesarias para un modelo de regresión lineal simple por mínimos cuadrados se cumplen es mediante el estudio de los residuos del modelo.
En R, los residuos se almacenan dentro del modelo bajo el nombre de residuals. R genera automáticamente los gráficos más típicos para la evaluación de los residuos de un modelo.

Estudio de los residuos del modelo.

Para confirmar que un modelo de regresión lineal simple 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:

[sourcecode language=»r» wraplines=»false» collapse=»false»]
par(mfrow = c(1,2)) # Indicamos que queremos 2 gráficos por pantalla (1 fila y 2 columnas)
plot(modelo_simple) # Pintamos los graficos de residuos
[/sourcecode]

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.

[sourcecode language=»r» wraplines=»false» collapse=»false»]
plot(x = modelo_simple$fitted.values, y = abs(rstudent(modelo_simple)), main = "Residuos absolutes studentizados vs valores predichos", pch = 20, col = "grey30")
abline(h = 3, col = "red") # se pinta la linea límite de 3
[/sourcecode]

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.

[sourcecode language=»r» wraplines=»false» collapse=»false»]
plot(hatvalues(modelo_simple), 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_simple$model)[2]-1 + 1)/dim(modelo_simple$model)[1]), col = "red")
[/sourcecode]