Regresión Práctica
Oscar Ramírez-Alán
19 de octubre de 2016
Regresión Lineal simple
- Utilizaremos la base de datos “trees”.
Observe que la base de datos es interactiva, para visualizar la cantidad de datos que usted requiera.
datatable(trees) # Interactive table Recomiendo hacer inspección básica, por ejemplo, revisar el supuesto de que las variables se ajustan a un modelo de regresión lineal.
scatterplotMatrix(trees)
Note que Volume y Height pareciera no ajustarse de manera lineal
Modelo de regresión
reg1 <- lm(Volume ~ Height + Girth, data=trees)
summary(reg1)##
## Call:
## lm(formula = Volume ~ Height + Girth, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4065 -2.6493 -0.2876 2.2003 8.4847
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57.9877 8.6382 -6.713 2.75e-07 ***
## Height 0.3393 0.1302 2.607 0.0145 *
## Girth 4.7082 0.2643 17.816 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared: 0.948, Adjusted R-squared: 0.9442
## F-statistic: 255 on 2 and 28 DF, p-value: < 2.2e-16Revisamos los supuestos
par(mfrow=c(2,2))
plot(reg1)
Parece haber ciertos problemas en el cumplimiento de los supuestos!
El primer gráfico enfrenta los errores residuales frente a sus valores ajustados. Los residuos deben estar distribuidos al azar alrededor de la línea horizontal que representa un error residual de cero; es decir, no debe haber una tendencia clara en la distribución de puntos. Una tendencia en la variabilidad de los residuos sugiere que la varianza está relacionada con la media, violando el supuesto de varianza constante. Si el gráfico tiene forma de embudo, es decir, si los puntos parecen estar más o menos extendidos a lo largo del gráfico, entonces lo más probable es que exista heterocedastididad en los datos Así mismo, si se presentan agrupados los residuos, es posible que la aleatoriedad no se cumpla.
El segundo gráfico: residuos tipificados se trazan contra los cuantiles de una distribución normal estándar. Si los residuos se distribuyen normalmente los datos se deben situarse a lo largo de la línea.
El tercer gráfico escala-ubicación en el que los residuos están estandarizados por sus desviaciones estándar estimadas. Esta gráfica se utiliza para detectar si la difusión de los residuos es constante en el rango de valores ajustados. Por lo tanto una variación de la tendencia, es que existen valores altos muestran una gran variación, y se puede traducir como el no cumplimiento de la homogeneidad de las varianzas
El cuarto gráfico muestra el valor leverage de cada punto, que se relaciona con la medida de su importancia en la determinación del modelo de regresión. Están representados los datos que ejercen mayor influencia.
A veces es bueno correr algunas pruebas alrternativas para identificar esos valores influyentes de la regresión: Ejm
library(car)
influencePlot(reg1, id.n = 2)
## StudRes Hat CookD
## 18 -1.859901 0.1434615 0.1775359
## 20 -1.105744 0.2112366 0.1082855
## 31 2.765603 0.2270585 0.6052326#Observe que los valores son la posición dentro de los datosDe momento nos queda primeramente solucionar el cumplimeinto de los supuestos.
Por lo que ejecutaremos una transformación de los datos.
reg2 <- lm(log(Volume) ~ log(Height) + log(Girth), data=trees)
summary(reg2)##
## Call:
## lm(formula = log(Volume) ~ log(Height) + log(Girth), data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.168561 -0.048488 0.002431 0.063637 0.129223
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.63162 0.79979 -8.292 5.06e-09 ***
## log(Height) 1.11712 0.20444 5.464 7.81e-06 ***
## log(Girth) 1.98265 0.07501 26.432 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08139 on 28 degrees of freedom
## Multiple R-squared: 0.9777, Adjusted R-squared: 0.9761
## F-statistic: 613.2 on 2 and 28 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(reg2)
Note como mejoran los supuestos
library(visreg)
par(mfrow=c(2,2))
visreg(reg2)
Forma la parte gráfica intereactiva en ggplot (se requieren de ciertos paquetes)
library(ggplot2)#ggplot2 es una extension poderosa para graficar
plot1<-ggplot(trees, aes(x=log(Height), y=log(Volume))) +
geom_point(shape=1) + # genera circulos en el grafico
geom_smooth(method=lm)
ggplotly(plot1) # Make the plot interactive## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`## Warning in plyr::split_indices(scale_id, n): '.Random.seed' is not an
## integer vector but of type 'NULL', so ignoredBuscando valores influyentes
library(car)
avPlots(reg2, id.n=2, id.cex=0.7)
identificamos valores extremos
outlierTest(reg2)##
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 18 -2.32572 0.027789 0.86146Extraemos el valor extremos y probamos un nuevo modelo
trees1<-trees[-c(3,15,17,18,20),]
trees1## Girth Height Volume
## 1 8.3 70 10.3
## 2 8.6 65 10.3
## 4 10.5 72 16.4
## 5 10.7 81 18.8
## 6 10.8 83 19.7
## 7 11.0 66 15.6
## 8 11.0 75 18.2
## 9 11.1 80 22.6
## 10 11.2 75 19.9
## 11 11.3 79 24.2
## 12 11.4 76 21.0
## 13 11.4 76 21.4
## 14 11.7 69 21.3
## 16 12.9 74 22.2
## 19 13.7 71 25.7
## 21 14.0 78 34.5
## 22 14.2 80 31.7
## 23 14.5 74 36.3
## 24 16.0 72 38.3
## 25 16.3 77 42.6
## 26 17.3 81 55.4
## 27 17.5 82 55.7
## 28 17.9 80 58.3
## 29 18.0 80 51.5
## 30 18.0 80 51.0
## 31 20.6 87 77.0reg3 <- lm(log(Volume) ~ log(Height) + log(Girth), data=trees1)
summary(reg3)##
## Call:
## lm(formula = log(Volume) ~ log(Height) + log(Girth), data = trees1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.152211 -0.057899 0.001432 0.059596 0.115497
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.95432 0.95074 -7.315 1.92e-07 ***
## log(Height) 1.20408 0.24134 4.989 4.79e-05 ***
## log(Girth) 1.96470 0.07111 27.630 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07021 on 23 degrees of freedom
## Multiple R-squared: 0.9842, Adjusted R-squared: 0.9829
## F-statistic: 718.1 on 2 and 23 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(reg3)
Nota: no se recomienda hacer mas de una extracción al modelo.
Comparación de los modelos
AIC(reg2,reg3) #El modelo reg1, fue descartado por incumplir supuestos, sin embargo hay que ponerle atención a la advertencia, que es un error muy clásico en la comparación de modelos.## Warning in AIC.default(reg2, reg3): models are not all fitted to the same
## number of observations## df AIC
## reg2 4 -62.71125
## reg3 4 -59.53002Revisión de Supuestos con cálculos estadísticos
Normalidad
shapiro.test(residuals(reg2))##
## Shapiro-Wilk normality test
##
## data: residuals(reg2)
## W = 0.95922, p-value = 0.2782Homogeneidad de las varianzas
library(car)
ncvTest(reg2)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.2266356 Df = 1 p = 0.6340298Explicación
summary(reg2)##
## Call:
## lm(formula = log(Volume) ~ log(Height) + log(Girth), data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.168561 -0.048488 0.002431 0.063637 0.129223
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.63162 0.79979 -8.292 5.06e-09 ***
## log(Height) 1.11712 0.20444 5.464 7.81e-06 ***
## log(Girth) 1.98265 0.07501 26.432 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08139 on 28 degrees of freedom
## Multiple R-squared: 0.9777, Adjusted R-squared: 0.9761
## F-statistic: 613.2 on 2 and 28 DF, p-value: < 2.2e-16Observemos el Adjusted R-squared:0.98, que mide la variabilidad de los datos explicada por el modelo. En nuestro caso decir que el 98% de la variabilidad de los datos fue recogida por el modelo, en otras palabras el modelo es bueno, además de ser significativo.
Así mismo la ecuación de la regresión múltiple quedaría escrita de la siguiente manera:
log(Volumen)= -6.63162 + (log(Height) * (1.11712)) + (log(Girth) * (1.98265))
Por lo que si quisieramos hacer una predicción, sustituimos los valores en la ecuación.
Vamos a predecir por ejemplo, cuando la atura de un árbol es de 65 y el “girth” de 10.
(-6.63162)+(log(65)*1.11712) + (log(10)*1.98265)## [1] 2.596892Por lo que el valor esperado del volumen sería “2.596892”. Note que este valor está basado en la transformación logaritmica, y para transformarlo a su escala original, invertimos el logaritmo aplicando “exp” (función exponencial)
exp(2.596892)## [1] 13.42196Por lo que 13.421 sería el valor en las unidades originales en las que se encuentran los datos.
La forma correcta de aplicar la predicción en R es utilizando la función “predict.lm”
Y como nuestro modelo es bueno para predecir, aplicamos para un un conjunto los valores los siguiente.
predict.lm(reg2,data.frame(Height=c(65,72,75,90), Girth=c(10,15,20,25)))## 1 2 3 4
## 2.596908 3.515062 4.131038 4.777129**2) Utilizaremos la base de datos “mtcars”“**
a) genere un modelo de regresión utilizando únicamente las filas 1:15
pairs(mtcars2,pch=8)
b) Genere su modelo de regresión utilizando como variable dependiente a “mpg” y como varibles independientes a “hp+wt”.
Con esto, vamos a explicar el consumo (mpg) en función de la potencia (hp) y del peso (wt):
##
## Call:
## lm(formula = mpg ~ hp + wt, data = mtcars2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9656 -0.9031 -0.1894 0.9988 1.8708
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.556402 1.846469 18.173 4.25e-10 ***
## hp -0.040534 0.009174 -4.419 0.000838 ***
## wt -2.668977 0.692562 -3.854 0.002294 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.313 on 12 degrees of freedom
## Multiple R-squared: 0.8929, Adjusted R-squared: 0.8751
## F-statistic: 50.03 on 2 and 12 DF, p-value: 1.507e-06Esto se interpreta: Cuanto más potente es el coche, menos millas recorre (de ahí el signo negativo de su coeficiente), y cuanto más pesa, menos millas recorre.
El R-squared es del 87.5% , lo que quiere decir que esas dos variables explican bastante bien el consumo.
library(visreg)
visreg(regmt1)
c) Predecir, las millas recorridas por galón de un coche que tiene 150 caballos y pesa 2.5 (x 1000 lbs)
## 1
## 20.80379Interacción
regmt2<-lm(mpg~hp*wt, data=mtcars2)
summary(regmt2)##
## Call:
## lm(formula = mpg ~ hp * wt, data = mtcars2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1949 -0.6462 -0.2770 0.8963 1.8546
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.676528 7.741540 3.704 0.00348 **
## hp -0.010234 0.047549 -0.215 0.83352
## wt -1.162528 2.423634 -0.480 0.64086
## hp:wt -0.009041 0.013908 -0.650 0.52899
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.346 on 11 degrees of freedom
## Multiple R-squared: 0.8969, Adjusted R-squared: 0.8688
## F-statistic: 31.89 on 3 and 11 DF, p-value: 1.009e-05En este caso no hubo un interación sigmnificativa, lo que se interpreta que ninguna de las variables independientes se comporta de manera diferente con respecto a las variable dependiente.