Supuestos del modelo
Wilson Sandoval
10/9/2020
library(MASS)
library(ISLR)
library(psych)
data("Boston")
summary(Boston)
## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00
multi.hist(x = Boston[,1:3], dcol = c("blue","red"), dlty = c("dotted", "solid"), main = "")
modelo_multiple <- lm(formula = medv ~ ., data = Boston)
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
## chas 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-16
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.
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 chas 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.522553
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)
# 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 ***
## chas 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-16
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.
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
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.
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.
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.
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.
library(MASS)
library(carData)
library(car)
qqPlot(rstudent(modelo_multiple))
## [1] 369 372
shapiro.test(rstudent(modelo_multiple))
##
## Shapiro-Wilk normality test
##
## data: rstudent(modelo_multiple)
## W = 0.89411, p-value < 2.2e-16
Test de Kolmogorov-Smirnov
El test de Kolmogorov-Smirnov permite estudiar si una muestra procede de una población con una determinada distribución (media y desviación típica), no está limitado únicamente a la distribución normal. Se ejecuta con la función ks.test().
ks.test(x =rstudent(modelo_multiple) ,"pnorm", mean(rstudent(modelo_multiple)), sd(rstudent(modelo_multiple)))
##
## One-sample Kolmogorov-Smirnov test
##
## data: rstudent(modelo_multiple)
## D = 0.12843, p-value = 1.125e-07
## alternative hypothesis: two-sided
library(lmtest)
bptest(modelo_multiple)
##
## studentized Breusch-Pagan test
##
## data: modelo_multiple
## BP = 59.907, df = 11, p-value = 9.647e-09
#Independencia
## Prueba de Durbin-Watson
dwtest(modelo_multiple)
##
## Durbin-Watson test
##
## data: modelo_multiple
## DW = 1.0779, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
Ejemplo cumpliendo los supuestos
gen_dat <- function(n) {
varianza <- 16
x <- runif(n=n, min=-5, max=6)
media <- 4 - 6 * x
y <- rnorm(n=n, mean=media, sd=sqrt(varianza))
marco_datos <- data.frame(y=y, x=x)
return(marco_datos)
}
datos <- gen_dat(n=500)
mod <- lm(y ~ x, data=datos)
Los gráficos de residuales explicados anteriormente se pueden obtener usando la función plot sobre el modelo ajustado mod.
par(mfrow=c(2, 2))
plot(mod, las=1, col='deepskyblue4', which=1:3)
Ejemplo violando los supuestos
gen_dat <- function(n) {
varianza <- 16
x <- runif(n=n, min=-5, max=6)
media <- 4 - 6 * x + 2 * x^2
y <- rnorm(n=n, mean=media, sd=sqrt(varianza))
marco_datos <- data.frame(y=y, x=x)
return(marco_datos)
}
datos <- gen_dat(n=500)
mod <- lm(y ~ x, data=datos)
Los gráficos de residuales explicados anteriormente se pueden obtener usando la función plot sobre el modelo ajustado mod.
par(mfrow=c(2, 2))
plot(mod, las=1, col='darkseagreen3', which=1:3)
Punto atípico (outlier) y punto influyente
Los conceptos de atípico e influyente son diferentes y se definen así:
- Punto atípico (outlier): es una observación que es numéricamente distante del resto de los datos.
- Punto influyente: punto que tiene impacto en las estimativas del modelo.
Prueba de Bonferroni para detectar outliers
\(H_0:\) la observación i-ésima NO es un outlier frente a
\(H_1:\) la observación i-ésima SI es un outlier.
library(car)
outlierTest(modelo_multiple, cutoff=0.05, n.max=10, order=TRUE, labels=names(rstudent(modelo_multiple)))
## rstudent unadjusted p-value Bonferroni p
## 369 5.893600 7.0113e-09 3.5477e-06
## 372 5.500418 6.0950e-08 3.0841e-05
## 373 5.325354 1.5341e-07 7.7626e-05
Distancia de Cook
Es una medida de cómo influye la observación i-ésima sobre la estimación de
\(\beta\) al ser retirada del conjunto de datos. Una distancia de Cook grande significa que una observación tiene un peso grande en la estimación de \(\beta\)
Una forma sencilla de calcular los valores de la distancia de Cook \(D_i\) es por medio de la expresión. .
\[D_i = \frac{d_i^2}{k+1} \times \frac{h_{ii}}{1-h_{ii}},\]
Son puntos influyentes las observaciones que presenten
\[D_i=\frac{4}{n-p-2}\]
influenceIndexPlot(modelo_multiple, vars="Cook")
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)), ylim=c(0,4), main = "Residuos absolutes studentizados vs valores predichos", pch = 20, col = "grey30")
abline(h = 3, col = "red") # se pinta la linea límite de 3
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")
Colinealidad
Para la colinealidad se recomienda calcular el coeficiente de correlación entre cada par de predictores incluidos en el modelo:
library(corrplot) # importa la libreria para generar graficos de correlación
## corrplot 0.84 loaded
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.
library(car) # importa el paquete car
vif(modelo_multiple) # calcula el factor de inflación de varianza
## crim zn chas nox rm dis rad tax
## 1.789704 2.239229 1.059819 3.778011 1.834806 3.443420 6.861126 7.272386
## ptratio black lstat
## 1.757681 1.341559 2.581984
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.
Evaluación del modelo de regresión
- El objetivo al realizar un modelo de regresión es predecir con precisión un caso desconocido.
- Realizar una evaluación del modelo de regresión después de crear el modelo. Para realizarlo podemos usar el entrenamiento y la prueba en el mismo conjunto de datos.
También se puede realizar una división del conjunto de datos para entrenamiento y otro para prueba. La primera solución consiste en seleccionar una parte de nuestro conjunto de datos para la realización de pruebas. Usamos todo el conjunto de datos para el entrenamiento, y construimos un modelo usando este conjunto de entrenamiento. Luego seleccionamos una pequeña porción del conjunto de datos, que se utiliza como valores reales del conjunto de pruebas. Pasamos este conjunto por nuestro modelo construido, pronosticando valores. Por último, comparamos los valores pronosticados por nuestro modelo con los valores reales.
Modelo Train/Test Split
Una forma de mejorar la precisión fuera de la muestra consiste en utilizar otro enfoque de evaluación denominado “Train/Test Split” (División entre entrenamiento y prueba). En este enfoque, seleccionamos una parte de nuestro conjunto de datos para entrenamiento y el resto se utiliza para probar. El modelo se basa en el conjunto de formación y luego se pasa para la predicción. Finalmente, los valores pronosticados para el conjunto de pruebas se comparan con los valores reales del conjunto de datos.
Evaluación del modelo
Hay diferentes métricas de la evaluación del modelo de regresión, pero la mayoría de ellos se basan en la similitud de los valores pronosticados y reales. Una de las métricas más simples para calcular la precisión de nuestro modelo
Metricas de evaluación
- Error medio absoluto (MAE) es la media del valor absoluto de los errores. Es la medida más fácil de entender, ya que es sólo el error promedio.
- Error cuadrático medio (MSE) es la media de los errores al cuadrado. Es más popular que el error medio absoluto porque el enfoque se orienta más hacia grandes errores. Esto se debe a que el término al cuadrado aumenta exponencialmente los errores más grandes en comparación con los más pequeños.
- Raiz cuadrada del error cuadrático medio (RMSE) es la raíz cuadrada de la anterior medida. Esta es una de las métricas más populares de las métricas de evaluación porque es interpretable en las mismas unidades que el vector de respuesta (o unidades y) haciendo fácil de correlacionar con la información.
- Error absoluto relativo (RAE), también conocido como la suma residual de cuadrado, donde la barra de y es un valor medio de y, toma el error absoluto total y la normaliza dividiendo por el error absoluto total del predictor simple.
- Error cuadrático relativo (RSE) es muy similar a “Error absoluto relativo”, pero se usa para calcular R-cuadrado. R-cuadrado no es un error, si no, una métrica popular para la precisión de su modelo. Representa como de cerca están los valores de los datos de la línea de regresión ajustada. Cuanto más alto sea el R-cuadrado, mejor encaja el modelo a los datos.