Preparación conjunto de datos
Usaremos el conjunto de datos del dataset de R llamado mtcars, con el objetivo de estimar o predecir las millas por galón(mpg) en función a las variables continuas que tiene el dataset. Las variables a utilizar son:
Mpg: Millas por galón
Disp: Desplazamiento
Hp: Potencia Bruta
Drat: engranaje del eje trasero
Wt: Peso en libras
Qsec: 1/4 de milla de tiempo
##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionAnálisis de la relación entre variables
Lo primero es realizar un analisis exploratorio para analizar cuales pueden ser los mejores predictores para el modelo, identificando relaciones lineales y no lineales con las diferentes variables y colinealidad entre predictores.
Para esto utilizaremos la librería GGally, la cual nos permite ver los diagramas de dispersión, los valores de correlación para cada par de variables y la distribución de cada una de las variables en un solo gráfico.## Loading required package: ggplot2##
## Attaching package: 'ggplot2'## The following object is masked from 'datos':
##
## mpg## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2##
## Attaching package: 'GGally'## The following object is masked from 'package:dplyr':
##
## nasaggpairs(datos, lower = list(continuous = "smooth"),
diag = list(continuous = "bar"), axisLabels = "none")## Warning in check_and_set_ggpairs_defaults("diag", diag, continuous =
## "densityDiag", : Changing diag$continuous from 'bar' to 'barDiag'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
De la gráfica, podemos observar es que las 3 variables que mayor relacion lineal tienen con las mpg, son el peso(wt),volumen del motor(disp) y caballos de fuerza(hp).Sin embargo, se debe tener cuidado con disp ya que este a su vez, posee una alta correlacion con wt y hp, por lo que tener todas las variables juntas podrian generar problemas de colinealidad.
Generar el modelo
Para la generación del modelo, lo primero que realizamos es un modelo inicial con todas las variables de interés, para despues seleccionar los mejores predictores.De este, podemos ver como tiene un coeficiente de determinacion alto (0.83) lo que da a entender que explica en un gran porcentaje la variabilidad observada en las mpg.Aun asi, se tienen variables que no estan contribuyendo mucho al modelo, por lo que se podria deshacer de ellas.#Creamos el modelo con las variables de interés iniciales
modelo_parcial <- lm (mpg ~ disp + hp + drat + wt, data = datos)
#Resumen del modelo
summary(modelo_parcial)##
## Call:
## lm(formula = mpg ~ disp + hp + drat + wt, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5077 -1.9052 -0.5057 0.9821 5.6883
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.148738 6.293588 4.631 8.2e-05 ***
## disp 0.003815 0.010805 0.353 0.72675
## hp -0.034784 0.011597 -2.999 0.00576 **
## drat 1.768049 1.319779 1.340 0.19153
## wt -3.479668 1.078371 -3.227 0.00327 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.602 on 27 degrees of freedom
## Multiple R-squared: 0.8376, Adjusted R-squared: 0.8136
## F-statistic: 34.82 on 4 and 27 DF, p-value: 2.704e-10Selección de los mejores predictores
Partiendo del modelo inicial anterior, utilizamos la funcion step() de R para la implementacion de la estrategia de stepwise mixto, con el fin de idenficar los mejores predictores para el modelo. #
##
## Call:
## lm(formula = mpg ~ hp + wt, data = datos)
##
## Coefficients:
## (Intercept) hp wt
## 37.22727 -0.03177 -3.87783En base al resultado anterior, tenemos que el mejor modelo obtenido esta compuesto por las variables hp y wt.
##
## Call:
## lm(formula = mpg ~ hp + wt, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.941 -1.600 -0.182 1.050 5.854
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.22727 1.59879 23.285 < 2e-16 ***
## hp -0.03177 0.00903 -3.519 0.00145 **
## wt -3.87783 0.63273 -6.129 1.12e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.593 on 29 degrees of freedom
## Multiple R-squared: 0.8268, Adjusted R-squared: 0.8148
## F-statistic: 69.21 on 2 and 29 DF, p-value: 9.109e-12En el modelo final, podemos ver como a pesar de utilizar menos variables seguimos conservando un alto coeficiente de determinacion(0.82) y un valor p significativo (9.109e-12), con respecto al generado inicialmente, ya que se preservaron las variables que mas aportaban al modelo.
A demas, tenemos los siguientes intervalos de confianza para cada uno de los coeficientes parciales de regresion, los cuales representan como varia en promedio la variable (Y) cada que un predictor varia y el resto se mantiene constantes, en este caso podemos ver como en terminos generales los predictores son significativos para el modelo.## 2.5 % 97.5 %
## (Intercept) 33.95738245 40.49715778
## hp -0.05024078 -0.01330512
## wt -5.17191604 -2.58374544Validación de condiciones para la regresión múltiple lineal
Relación lineal entre los predictores y la variable respuesta
En este caso, realizaremos la validacion mediante diagramas de dispersión entre cada uno de los predictores y los residuos del modelo,donde los residuos deben de distribuirse aleatoriamente en torno a 0 con una variabilidad constante a lo largo del eje X.
##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combineplot1 <- ggplot(data = datos, aes(hp, mejor_modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos, aes(wt, mejor_modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
como podemos ver, en terminos generales se cumple con esta condicion de linealidad, exceptuando algunos casos donde puede dar indicacion de datos atipicos.
Distribución normal de los residuos
De manera similar al caso anterior, tenemos que en terminos generales los residuos se distribuyen de manera normal a excepcion de algunos posibles valores atipicos, y ademas se tiene la prueba de Shapiro-Wilk que nos da a entender que hay normalidad.
##
## Shapiro-Wilk normality test
##
## data: mejor_modelo$residuals
## W = 0.92792, p-value = 0.03427Variabilidad constante de los residuos (homocedasticidad)
En esta condicion podemos ver que se presentan algunos problemas de homocedasticidad, ya que los residuos frente a los valores ajustados por el modelo, se deben de distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X,pero en este caso vemos mayor dispersión en los extremos, lo que puede significar que la variabilidad tiene cierta dependencia del valor ajustado. # Aun asi, mediante el el test de Breusch-Pagan, que nos da a indicar que se cumple con la homocedasticidad.ggplot(data = datos, aes(mejor_modelo$fitted.values, mejor_modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## studentized Breusch-Pagan test
##
## data: mejor_modelo
## BP = 0.88072, df = 2, p-value = 0.6438No multicolinialidad
Mediante la siguiente matriz de correlación entre predictores, podemos observar que si bien hay cierto porcentaje de correlacion entre estos, no es demasiado significativo para que se presenten problemas de multicolinialidad. #
## corrplot 0.84 loaded
Ademas, mediante el calculo del factor de inflacion de la varianza para cada uno de los predictores, tenemos una medicion mas exacta que nos da a indicar que no se presentan problemas de multicolinialidad.
## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## hp wt
## 1.766625 1.766625Autocorrelacion
Mediante el test de Durbin-Watson, podemos verificar si los residuos del modelo estan correlacionados, en este caso no hay evidencia de correlacion.
## lag Autocorrelation D-W Statistic p-value
## 1 0.2954091 1.362399 0.038
## Alternative hypothesis: rho != 0