REGRESION LINEAL MULTIPLE
JOSÉ BAEZA D.
30-01-2020
regresión simple y múltiple se han obtenido del libro Introduction to Statistical Learning. El objetivo es mostrar los principales comandos en R para generar modelos lineales. Para obtener un modelo final robusto se tiene que analizar con más detalle cada una de las condiciones que se requieren para estos métodos. Para ver los conceptos teóricos de la regresión lineal múltiple consultar capítulo Introducción a la Regresión Lineal Múltiple. También se incluyen soluciones personales (validadas con otras en la web) de ejercicios propuestos en el libro.
Se van a emplear funciones contenidas en el paquete MASS y datos del paquete ISL
require(MASS)## Loading required package: MASSrequire(ISLR)## Loading required package: ISLRdata("Boston")El dataset Boston del paquete MASS recoge la mediana del valor de la vivienda en 506 áreas residenciales de Boston. Junto con el precio, se han registrado 13 variables adicionales.
crim: ratio de criminalidad per cápita de cada ciudad. zn: Proporción de zonas residenciales con edificaciones de más de 25.000 pies cuadrados. indus: proporción de zona industrializada. chas: Si hay río en la ciudad (= 1 si hay río; 0 no hay). nox: Concentración de óxidos de nitrógeno (partes per 10 millón). rm: promedio de habitaciones por vivienda. age: Proporción de viviendas ocupadas por el propietario construidas antes de 1940. dis: Media ponderada de la distancias a cinco centros de empleo de Boston. rad: Índice de accesibilidad a las autopistas radiales. tax: Tasa de impuesto a la propiedad en unidades de $10,000. ptratio: ratio de alumnos/profesor por ciudad. black: 1000(Bk - 0.63)^2 donde Bk es la proporción de gente de color por ciudad. lstat: porcentaje de población en condición de pobreza. medv: Valor mediano de las casas ocupadas por el dueño en unidades de $1000s.
En primer lugar se realiza un análisis básico de los datos de forma numérica y gráfica.
library(psych)
# La variable chas es una variable categórica por lo que se transforma a factor
Boston$chas <- as.factor(Boston$chas)
summary(Boston)## crim zn indus chas nox
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 0:471 Min. :0.3850
## 1st Qu.: 0.08204 1st Qu.: 0.00 1st Qu.: 5.19 1: 35 1st Qu.:0.4490
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.5380
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.5547
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.6240
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :0.8710
## rm age dis rad
## Min. :3.561 Min. : 2.90 Min. : 1.130 Min. : 1.000
## 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100 1st Qu.: 4.000
## Median :6.208 Median : 77.50 Median : 3.207 Median : 5.000
## Mean :6.285 Mean : 68.57 Mean : 3.795 Mean : 9.549
## 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188 3rd Qu.:24.000
## Max. :8.780 Max. :100.00 Max. :12.127 Max. :24.000
## tax ptratio black lstat
## Min. :187.0 Min. :12.60 Min. : 0.32 Min. : 1.73
## 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38 1st Qu.: 6.95
## Median :330.0 Median :19.05 Median :391.44 Median :11.36
## Mean :408.2 Mean :18.46 Mean :356.67 Mean :12.65
## 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23 3rd Qu.:16.95
## Max. :711.0 Max. :22.00 Max. :396.90 Max. :37.97
## medv
## Min. : 5.00
## 1st Qu.:17.02
## Median :21.20
## Mean :22.53
## 3rd Qu.:25.00
## Max. :50.00# Dado que hay muchas variables, se grafican por grupos de 4, excluyendo las
# categóricas
multi.hist(x = Boston[,1:3], dcol = c("blue","red"), dlty = c("dotted", "solid"),
main = "")
multi.hist(x = Boston[,5:9], dcol = c("blue","red"), dlty = c("dotted", "solid"),
main = "")
multi.hist(x = Boston[,10:14], dcol = c("blue","red"),
dlty = c("dotted", "solid"), main = "")
Regresión Lineal Simple
modelo_simple <- lm(data = Boston,formula = medv ~ lstat)
names(modelo_simple)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"summary(modelo_simple)##
## Call:
## lm(formula = medv ~ lstat, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.168 -3.990 -1.318 2.034 24.500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.55384 0.56263 61.41 <2e-16 ***
## lstat -0.95005 0.03873 -24.53 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.216 on 504 degrees of freedom
## Multiple R-squared: 0.5441, Adjusted R-squared: 0.5432
## F-statistic: 601.6 on 1 and 504 DF, p-value: < 2.2e-16confint(modelo_simple, level = 0.95)## 2.5 % 97.5 %
## (Intercept) 33.448457 35.6592247
## lstat -1.026148 -0.8739505Como era de esperar dado que el p-value del predictor lstat ha resultado significativo para un α=0.05, su intervalo de confianza del 95% no contiene el valor 0.
Una vez generado el modelo, es posible predecir el valor de la vivienda sabiendo el estatus de la población en la que se encuentra. Toda predicción tiene 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 todas las viviendas que se encuentren en una población con un determinado porcentaje de pobreza, supóngase lstat=10.
predict(object = modelo_simple, newdata = data.frame(lstat = c(10)),
interval = "confidence", level = 0.95)## fit lwr upr
## 1 25.05335 24.47413 25.63256Intervalo de predicción: Devuelve un intervalo para el valor esperado de una vivienda en particular que se encuentre en una población con un determinado porcentaje de pobreza.
predict(object = modelo_simple, newdata = data.frame(lstat = c(10)),
interval = "prediction", level = 0.95)## fit lwr upr
## 1 25.05335 12.82763 37.27907Como es de esperar ambos intervalos están centrados en torno al mismo valor. Si bien ambos parecen 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.
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.
attach(Boston)
plot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30")
abline(modelo_simple, lwd = 3, col = "red")
La representación gráfica de las observaciones muestra que la relación entre ambas variables estudiadas no es del todo lineal, lo que apunta a que otro tipo de modelo podría explicar mejor la relación. Aun así la aproximación no es mala.
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.
par(mfrow = c(1,2))
plot(modelo_simple)
par(mfrow = c(1,1))plot(x = modelo_simple$fitted.values, y = abs(rstudent(modelo_simple)),
main = "Absolute studentized residuals vs predicted values", pch = 20,
col = "grey30")
abline(h = 3, col = "red")
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")
En este caso muchos de los valores parecen posibles outliers o puntos con alta influencia porque los datos realmente no se distribuyen de forma lineal en los extremos.
Modelo precio medio vivienda=34.55−0.95lstat
Regresión múltiple
Se desea generar un modelo que permita explicar el precio de la vivienda de una población empleando para ello cualquiera de las variables disponibles en el dataset Boston y que resulten útiles en el modelo.
R permite crear un modelo con todas las variables incluidas en un data.frame de la siguiente forma:
modelo_multiple <- lm(formula = medv ~ ., data = Boston)
# También se pueden especificar una a una
summary(modelo_multiple)##
## Call:
## lm(formula = medv ~ ., data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.595 -2.730 -0.518 1.777 26.199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 ***
## crim -1.080e-01 3.286e-02 -3.287 0.001087 **
## zn 4.642e-02 1.373e-02 3.382 0.000778 ***
## indus 2.056e-02 6.150e-02 0.334 0.738288
## chas1 2.687e+00 8.616e-01 3.118 0.001925 **
## nox -1.777e+01 3.820e+00 -4.651 4.25e-06 ***
## rm 3.810e+00 4.179e-01 9.116 < 2e-16 ***
## age 6.922e-04 1.321e-02 0.052 0.958229
## dis -1.476e+00 1.995e-01 -7.398 6.01e-13 ***
## rad 3.060e-01 6.635e-02 4.613 5.07e-06 ***
## tax -1.233e-02 3.760e-03 -3.280 0.001112 **
## ptratio -9.527e-01 1.308e-01 -7.283 1.31e-12 ***
## black 9.312e-03 2.686e-03 3.467 0.000573 ***
## lstat -5.248e-01 5.072e-02 -10.347 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
## F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16El p-value obtenido para el estadístico F es muy pequeño (< 2.2e-16) lo que indica que al menos uno de los predictores introducidos en el modelo está relacionado con la variable respuesta medv. El modelo es capaz de explicar el 74% de la variabilidad observada en el precio de la vivienda (R2=0.74)
En el summary se puede observar que algunos predictores tienen p-values muy altos, sugiriendo que no contribuyen al modelo por lo que deben ser excluidos, por ejemplo age e indus. La exclusión de predictores basándose en p-values no es aconsejable, en su lugar se recomienda emplear métodos de best subset selection, stepwise selection (forward, backward e hybrid) o Shrinkage/regularization. Para una descripción detallada de cada uno ver capítulo Selección de predictores y mejor modelo: Subset selection, Ridge, Lasso y dimension reduction.
step(modelo_multiple, direction = "both", trace = 0)##
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad +
## tax + ptratio + black + lstat, data = Boston)
##
## Coefficients:
## (Intercept) crim zn chas1 nox rm
## 36.341145 -0.108413 0.045845 2.718716 -17.376023 3.801579
## dis rad tax ptratio black lstat
## -1.492711 0.299608 -0.011778 -0.946525 0.009291 -0.522553La selección de predictores empleando stepwise selection (hybrid/doble) ha identificado como mejor modelo el formado por los predictores crim, zn, chas, nox, rm, dis, rad, tax, ptratio, black, lstat
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)##
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad +
## tax + ptratio + black + lstat, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5984 -2.7386 -0.5046 1.7273 26.2373
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.341145 5.067492 7.171 2.73e-12 ***
## crim -0.108413 0.032779 -3.307 0.001010 **
## zn 0.045845 0.013523 3.390 0.000754 ***
## chas1 2.718716 0.854240 3.183 0.001551 **
## nox -17.376023 3.535243 -4.915 1.21e-06 ***
## rm 3.801579 0.406316 9.356 < 2e-16 ***
## dis -1.492711 0.185731 -8.037 6.84e-15 ***
## rad 0.299608 0.063402 4.726 3.00e-06 ***
## tax -0.011778 0.003372 -3.493 0.000521 ***
## ptratio -0.946525 0.129066 -7.334 9.24e-13 ***
## black 0.009291 0.002674 3.475 0.000557 ***
## lstat -0.522553 0.047424 -11.019 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7348
## F-statistic: 128.2 on 11 and 494 DF, p-value: < 2.2e-16Correlación
require(corrplot)## Loading required package: corrplot## corrplot 0.84 loadedcorrplot.mixed(corr = cor(Boston[,c("crim", "zn", "nox", "rm", "dis", "rad",
"tax", "ptratio", "black", "lstat", "medv")],
method = "pearson"))
Interacción entre predictores
modelo <- lm(medv ~ lstat + age, data = Boston)
summary(modelo)##
## Call:
## lm(formula = medv ~ lstat + age, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.981 -3.978 -1.283 1.968 23.158
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.22276 0.73085 45.458 < 2e-16 ***
## lstat -1.03207 0.04819 -21.416 < 2e-16 ***
## age 0.03454 0.01223 2.826 0.00491 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.173 on 503 degrees of freedom
## Multiple R-squared: 0.5513, Adjusted R-squared: 0.5495
## F-statistic: 309 on 2 and 503 DF, p-value: < 2.2e-16Dado que 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)
Regresión Polinomial: incorporar no-linealidad a los modelos lineales
attach(Boston)## The following objects are masked from Boston (pos = 4):
##
## age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
## rm, tax, znplot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30")
abline(modelo_simple, lwd = 3, col = "red")