TALLER DE ANÁLISIS DE REGRESIÓN LINEAL

Resumen

En el siguiente taller se desea observar el comportamiento de la variable respuesta precio de venta de la casa, a partir de las variables predictoras: Impuestos, Número de baños, Tamaño del lote, Espacio vital, Número de puestos de garaje, Número de habitaciones, Número de dormitorios, Edad del hogar y Cantidad de chimeneas.

A continuación, se describe cada variable denotando a la variable respuesta, con la letra \(y\), mientras que a cada variable predictora con la letra \(x_{i}\).

\(x_{1}\)= Impuestos (local, escuela, condado)/1000

\(x_{2}\)= Número de baños

\(x_{3}\)= Tamaño del lote (pies cuadrados x 1000 )

\(x_{4}\)= Espacio vital (pies cuadrados x 1000)

\(x_{5}\)= Número de puestos de garaje

\(x_{6}\)= Número de habitaciones (rooms)

\(x_{7}\)= Número de dormitorios (bedrooms)

\(x_{8}\)= Edad del hogar (años)

\(x_{9}\)=Cantidad de chimeneas

\(y\)= precio de venta de la casa/1000

\(\epsilon\)= Error

La ecuación que se requiere estimar por medio del modelo de regresión lineal múltiple tiene la siguiente forma:

\[\hat{y}=\beta_{0}+\beta_{1} x_{1}+\beta_{2} x_{2}+\cdots+\beta_{k} x_{k}+\epsilon\]

Donde:

\(\beta_{0}\)= Representa el intercepto

\(\beta_{1}\)= El coeficiente del predictor \(x_{1}\)

\(\beta_{2}\)= El coeficiente del predictor \(x_{2}\)

\(\beta_{3}\)= El coeficiente del predictor \(x_{3}\)

\(\beta_{k}\)= El coeficiente del predictor \(x_{k}\)

# Resumen estadístico de la base de datos

summary(datos)

##        y               x1              x2              x3       
##  Min.   :25.90   Min.   :3.891   Min.   :1.000   Min.   :2.275  
##  1st Qu.:29.90   1st Qu.:5.057   1st Qu.:1.000   1st Qu.:4.855  
##  Median :33.70   Median :5.974   Median :1.000   Median :5.685  
##  Mean   :34.61   Mean   :6.405   Mean   :1.333   Mean   :6.033  
##  3rd Qu.:38.15   3rd Qu.:7.873   3rd Qu.:2.000   3rd Qu.:7.158  
##  Max.   :45.80   Max.   :9.142   Max.   :2.000   Max.   :9.890  
##        x4              x5               x6            x7              x8       
##  Min.   :0.975   Min.   : 0.000   Min.   :5.0   Min.   :2.000   Min.   : 3.00  
##  1st Qu.:1.161   1st Qu.: 1.000   1st Qu.:6.0   1st Qu.:3.000   1st Qu.:30.00  
##  Median :1.432   Median : 1.000   Median :6.0   Median :3.000   Median :40.00  
##  Mean   :1.384   Mean   : 1.833   Mean   :6.5   Mean   :3.167   Mean   :37.46  
##  3rd Qu.:1.577   3rd Qu.: 2.000   3rd Qu.:7.0   3rd Qu.:3.250   3rd Qu.:48.50  
##  Max.   :1.831   Max.   :15.000   Max.   :8.0   Max.   :4.000   Max.   :62.00  
##        x9      
##  Min.   :0.00  
##  1st Qu.:0.00  
##  Median :0.00  
##  Mean   :0.25  
##  3rd Qu.:0.25  
##  Max.   :1.00

# Nombres de las variables de la base de datos

names(datos)

##  [1] "y"  "x1" "x2" "x3" "x4" "x5" "x6" "x7" "x8" "x9"

# Estructura de la base de datos

str(datos)

## tibble [24 x 10] (S3: tbl_df/tbl/data.frame)
##  $ y : num [1:24] 25.9 29.5 27.9 25.9 29.9 29.9 30.9 28.9 35.9 31.5 ...
##  $ x1: num [1:24] 4.92 5.02 4.54 4.56 5.06 ...
##  $ x2: num [1:24] 1 1 1 1 1 1 1 1 1 1 ...
##  $ x3: num [1:24] 3.47 3.53 2.27 4.05 4.46 ...
##  $ x4: num [1:24] 0.998 1.5 1.175 1.232 1.121 ...
##  $ x5: num [1:24] 1 2 1 1 1 1 1 0 2 1 ...
##  $ x6: num [1:24] 7 7 6 6 6 6 7 6 6 6 ...
##  $ x7: num [1:24] 4 4 3 3 3 3 3 3 3 3 ...
##  $ x8: num [1:24] 42 62 40 54 42 56 51 32 32 30 ...
##  $ x9: num [1:24] 0 0 0 0 0 0 1 0 0 0 ...

1. Ajuste un modelo de regresión múltiple que relacione el precio de venta con las nueve variables regresoras.

# Ajuste del modelo1

modelo1 <-lm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9,data=datos)

# Resumen del modelo1

summary(modelo1)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9, 
##     data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6230 -1.2920 -0.3015  0.9506  4.5351 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 20.23164    5.39006   3.754  0.00214 **
## x1           1.24538    1.06496   1.169  0.26176   
## x2           5.22728    2.27619   2.297  0.03760 * 
## x3           0.15553    0.45629   0.341  0.73827   
## x4           1.89401    4.06618   0.466  0.64853   
## x5           0.59810    0.28330   2.111  0.05322 . 
## x6          -0.91380    2.22767  -0.410  0.68787   
## x7           0.92566    3.07345   0.301  0.76771   
## x8          -0.07754    0.06406  -1.211  0.24611   
## x9           2.76934    1.81362   1.527  0.14904   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.757 on 14 degrees of freedom
## Multiple R-squared:  0.8716, Adjusted R-squared:  0.7891 
## F-statistic: 10.56 on 9 and 14 DF,  p-value: 7.694e-05

names(modelo1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

# Residuales del modelo1

modelo1$df.residual

## [1] 14

El modelo de regresión multiple tiene la siguiente forma:

\[\hat{y}= 1.24538 x_{1} + 5.22728 x_{2}+ 0.15553 x_{3} + 1.89401x_{4} + 0.59810x_{5}-0.91380x_{6}+ 0.92566 x_{7} -0.07754x_{8} +2.76934x_{9}+20.23164\].

2. Prueba de signícancia de la regresión. ¿Qué conclusiones puedes sacar?.

El Residual standar error del modelo (RSE) indica que cualquier predicción se aleja 2.757 unidades del valor verdadero y el coeficiente de determinación \(R^2\) establece que 87.16% de la variabilidad observada en la variable respuesta \(y\), es explicado por las variables predictoras \(x_{1}\), \(x_{2}\), \(x_{3}\),….,\(x_{9}\) con un \(R_{ajustado}^2 =78.91%\).
El p-value del modelo es significativo \(7.694e-5<0.05\), por lo que rechazamos la hipotesis nula, esto quiere decir que algun coeficiente es distinto de cero. Esto implica que al menos una de las variables predictoras contribuye de forma significativa a la explicación de la variable respuesta, por lo que el modelo es útil.

3. Utilice las pruebas t para evaluar la contribución de cada variable regresora al modelo. Discute tus hallazgos.

La variable predictora \(x_{2}\) (Número de baños) tiene una relación estadisticamente significativa con la variable respuesta \(y\) (precio de venta de la casa/1000), porque tiene un p-value <0.05. Las variables \(x_{1}\), \(x_{3}\), \(x_{4}\), \(x_{5}\), \(x_{6}\), \(x_{7}\), \(x_{8}\) y \(x_{9}\) no tienen una relación estadisticamente significativa con la variable respuesta \(y\), dado que tienen un p-value mayor a 0.05, estas variables aportan poca información al momento de explicar la variabilidad de la variable respuesta \(y\), por lo tanto se pueden eliminar algunos predictores del modelo.

4. Encuentre el subconjunto de variables predictoras en el conjunto de datos que da como resultado el “mejor” modelo según algún criterio de información.

Selección de los mejores predictores por medio del método stepwise y AIC.

Método paso a paso (stepwise): emplea criterios matemáticos para decidir qué predictores contribuyen significativamente al modelo y en qué orden se introducen. Dentro de este método se diferencias tres estrategias:

Dirección forward: El modelo inicial no contiene ningún predictor, solo el parámetro β0. A partir de este se generan todos los posibles modelos introduciendo una sola variable de entre las disponibles. Aquella variable que mejore en mayor medida el modelo se selecciona. A continuación se intenta incrementar el modelo probando a introducir una a una las variables restantes. Si introduciendo alguna de ellas mejora, también se selecciona. En el caso de que varias lo hagan, se selecciona la que incremente en mayor medida la capacidad del modelo. Este proceso se repite hasta llegar al punto en el que ninguna de las variables que quedan por incorporar mejore el modelo.

Dirección backward: El modelo se inicia con todas las variables disponibles incluidas como predictores. Se prueba a eliminar una a una El método paso a paso requiere de algún criterio matemático para determinar si el modelo mejora o empeora con cada incorporación o extracción. Existen varios parámetros empelados, de entre los que destacan el Cp, AIC, BIC y R2ajustado, cada uno de ellos con ventajas e inconvenientes. El método Akaike(AIC) tiende a ser más restrictivo e introducir menos predictores que el R2-ajustado. Para un mismo set de datos, no todos los métodos tienen porque concluir en un mismo modelo.

# Selección de los mejores predictores

step(modelo1, direction = "both", trace = 1)

## Start:  AIC=55.75
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
## 
##        Df Sum of Sq    RSS    AIC
## - x7    1     0.690 107.12 53.901
## - x3    1     0.883 107.31 53.944
## - x6    1     1.279 107.71 54.032
## - x4    1     1.649 108.08 54.115
## <none>              106.43 55.746
## - x1    1    10.396 116.82 55.982
## - x8    1    11.140 117.57 56.135
## - x9    1    17.725 124.15 57.443
## - x5    1    33.884 140.31 60.379
## - x2    1    40.092 146.52 61.418
## 
## Step:  AIC=53.9
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x8 + x9
## 
##        Df Sum of Sq    RSS    AIC
## - x6    1     0.675 107.79 52.051
## - x3    1     0.862 107.98 52.093
## - x4    1     2.018 109.14 52.349
## <none>              107.12 53.901
## - x1    1     9.706 116.82 53.982
## - x8    1    10.457 117.57 54.136
## + x7    1     0.690 106.43 55.746
## - x9    1    19.724 126.84 55.957
## - x5    1    33.194 140.31 58.379
## - x2    1    42.154 149.27 59.865
## 
## Step:  AIC=52.05
## y ~ x1 + x2 + x3 + x4 + x5 + x8 + x9
## 
##        Df Sum of Sq    RSS    AIC
## - x3    1     1.149 108.94 50.306
## - x4    1     1.464 109.25 50.375
## <none>              107.79 52.051
## - x1    1     9.710 117.50 52.121
## - x8    1    17.647 125.44 53.690
## + x6    1     0.675 107.12 53.901
## - x9    1    19.340 127.13 54.012
## + x7    1     0.085 107.71 54.032
## - x5    1    34.429 142.22 56.704
## - x2    1    46.114 153.91 58.599
## 
## Step:  AIC=50.31
## y ~ x1 + x2 + x4 + x5 + x8 + x9
## 
##        Df Sum of Sq    RSS    AIC
## - x4    1     2.439 111.38 48.837
## <none>              108.94 50.306
## - x1    1    13.780 122.72 51.164
## + x3    1     1.149 107.79 52.051
## + x6    1     0.962 107.98 52.093
## + x7    1     0.197 108.74 52.263
## - x8    1    20.456 129.40 52.436
## - x9    1    26.213 135.15 53.481
## - x5    1    33.467 142.41 54.735
## - x2    1    45.150 154.09 56.628
## 
## Step:  AIC=48.84
## y ~ x1 + x2 + x5 + x8 + x9
## 
##        Df Sum of Sq    RSS    AIC
## <none>              111.38 48.837
## + x4    1     2.439 108.94 50.306
## + x3    1     2.124 109.25 50.375
## - x8    1    18.259 129.64 50.481
## + x6    1     0.184 111.19 50.798
## + x7    1     0.017 111.36 50.834
## - x9    1    24.345 135.72 51.582
## - x1    1    24.575 135.95 51.622
## - x5    1    32.666 144.04 53.010
## - x2    1    60.752 172.13 57.285

## 
## Call:
## lm(formula = y ~ x1 + x2 + x5 + x8 + x9, data = datos)
## 
## Coefficients:
## (Intercept)           x1           x2           x5           x8           x9  
##    19.65737      1.34545      5.69737      0.57455     -0.07865      2.53594

Las variables \(x_{3}\), \(x_{4}\), \(x_{6}\) y \(x_{7}\) son las variables eliminadas en el proceso de selección. Reajustamos el modelo excluyendo dicha variables:

El mejor modelo resultante del proceso de selección ha sido:

# Regresión del mejor modelo

modelo2 <-lm(y~x1+x2+x5+x8+x9,data=datos)

# Resumen del modelo2

summary(modelo2)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x5 + x8 + x9, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3423 -1.4806 -0.3722  1.2254  4.3120 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.65737    3.75142   5.240 5.54e-05 ***
## x1           1.34545    0.67512   1.993  0.06166 .  
## x2           5.69737    1.81828   3.133  0.00574 ** 
## x5           0.57455    0.25006   2.298  0.03379 *  
## x8          -0.07865    0.04579  -1.718  0.10298    
## x9           2.53594    1.27849   1.984  0.06277 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.488 on 18 degrees of freedom
## Multiple R-squared:  0.8657, Adjusted R-squared:  0.8283 
## F-statistic:  23.2 on 5 and 18 DF,  p-value: 2.899e-07

# Residuales del modelo2

modelo2$df.residual

## [1] 18

Tras excluir las variables \(x_{3}\), \(x_{4}\), \(x_{6}\) y \(x_{7}\), el \(R^2\) apenas se ha modificado y el valor del \(R_{ajustado}^2\) aumento, lo que indica que las variables excluidas aportan poca información al momento de explicar la variabilidad de la variable respuesta y.

5. Interprete los parámetros estimados del nuevo modelo.

En el nuevo modelo (modelo2) que se obtuvo al eliminar las variables que aportaban poca información en el modelo1, observamos que los coeficientes de regresión \(\beta_0\), \(\beta_2\) y \(\beta_5\) son significativos porque ambos tienen un p-value menor a 0.05.
Los coeficientes \(\beta_1\), \(\beta_8\) y \(\beta_9\) tienen un menor grado de significancia debido a que tienen un p-value mayor a 0.05.
El coeficiente \(\beta_2\) de la variable predictora \(x_{2}\), es el parámetro que aporta mas información al modelo2.
Manteniendo las demas variables constantes, por cada unidad que se incrementa el predictor \(x_{2}\) (número de baños), el precio de venta de la casa aumenta en promedio 5.69737 unidades.
Por cada unidad que aumente el valor del predictor \(x_{8}\) (Edad del hogar), el precio de venta de la casa disminuira en 0.07865 unidades. Es decir el valor de la vivienda se desvalua.

6.Construya e interprete un gráfico de residuales contra la respuesta estimada. Es sensato pensar que el modelo cumple los supuestos? Verifique esto a un nivel de signícancia del 5%

# gráfico de residuales contra la respuesta estimada.

plot1 <- ggplot(data = datos, aes(modelo2$residuals, y)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()

plot1

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

En el gráfico anterior se puede observar que existe linealidad entre la variable respuesta con cada uno de los predictores.
Se puede evidenciar en el gráfico que no se sigue un patrón ajustado y por tanto podemos inferir que hay dependecia entre la variable respuesta con los predictores.

# Análisis gráfico de residuos 

par(mfrow = c(2, 2))
plot(modelo2)

## Warning in sqrt(crit * p * (1 - hh)/hh): Se han producido NaNs

## Warning in sqrt(crit * p * (1 - hh)/hh): Se han producido NaNs

Supuestos del modelo2

1. Linealidad de los datos

plot(modelo2, 1)

En el gráfico de residuos se puede observar que no se sigue un patrón ajustado y que la linea roja se distribuye de manera aproximadamente horizontal en torno a cero, manteniendo la misma variabilidad a lo largo del eje x. Por lo tanto podemos asumir una relación lineal entre los predictores y la variable respuesta.

2. Normalidad de los residuos

plot(modelo2, 2)

Como todos los puntos caen aproximadamente a lo largo de la línea recta de referencia, podemos asumir la normalidad de los errores residuales.

Prueba de normalidad

\(H_{0}\): Los residuales tienen distribución normal
\(H_{1}\): Los residuales no tienen distribución normal

# Prueba de normalidad

shapiro.test(rstandard(modelo2))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(modelo2)
## W = 0.9491, p-value = 0.259

# Test de hipótesis para el análisis de normalidad de los residuos

shapiro.test(modelo2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo2$residuals
## W = 0.94818, p-value = 0.2474

Como el p-value >0.05, rechazamos la hipotesis alternativa. Por lo tanto, los residuos tienen una distribución normal.
Tanto el análisis gráfico como el test de hipótesis confirman la normalidad.

3. Homogeneidad de la varianza de los residuos

plot(modelo2, 3)

bptest(modelo2)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 3.1764, df = 5, p-value = 0.6728

 # gráfico 2D de los valores ajustados vs. los residuos estandarizados

plot(fitted.values(modelo2),rstandard(modelo2), xlab="Valores ajustados", ylab="Residuos estandarizados")  

abline(h=0) # dibuja la recta en cero

En los gráficos se observa que los residuos son homocedásticos. Por lo tanto podemos suponer que los residuos tienen una varianza constante (homocedasticidad).

4. Independencia

Residuales<-resid(modelo2)

plot(Residuales)

Para que se cumpla el supuesto de independencia debe haber un patrón disperso de puntos, es decir que el residual y el tiempo se han independientes. En este caso no se cumple el supuesto y podemos observar en el gráfico que existe un patrón de dependencia y posibles datos atípicos.

Prueba de durbin Watson para la prueba de independencia

\(H_{0}\): Independencia
\(H_{1}\): Autocorrelación (+ o -)

dwtest(modelo2,alternative="two.sided")

## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 2.0973, p-value = 0.8485
## alternative hypothesis: true autocorrelation is not 0

Existe autocorrelación positiva, por lo tanto, para un nivel de significancia del 0.05, rechazamos la hipotesis nula de la independencia.
Podemos observar que se cumplen los supuestos de linealidad, homocedasticidad y normalidad, pero no se cumple el supuesto de independencia.

7. Realice un análisis de diagnóstico e identífique (si es que existen) observaciones con alto leverage, influyentes y extremas en la respuesta.

Identificando valores atípicos, puntos de alto apalancamiento y observaciones extremas.

# gráfico de Residuos vs Apalancamiento

plot(modelo2, 5)

## Warning in sqrt(crit * p * (1 - hh)/hh): Se han producido NaNs

## Warning in sqrt(crit * p * (1 - hh)/hh): Se han producido NaNs

La gráfica anterior resalta los 3 puntos más extremos (17, 22 y 24), con un residuo estandarizado por debajo de 2.5.
24 es una observación atípica con un residuo estandarizado por encima de 2.
En la gáfica se puede observar que no hay valores atípicos que superen las 3 desviaciones estandar de los residuos.
hay observaciones que presentan un alto punto de apalancamiento por encima de 0.5. \(2(p + 1)/n = 2*(5+1)/24=12/24 = 0.5\).
Las observaciones 15 y 17 tienen un alto punto de apalancamiento.
La observación 17 es un valor influyente.

# Estadistica de apalancamiento 

hatvalues(modelo2)

##          1          2          3          4          5          6          7 
## 0.08843826 0.23011848 0.11966997 0.14545038 0.08366199 0.20949159 0.24520385 
##          8          9         10         11         12         13         14 
## 0.11397612 0.07718464 0.10234289 0.12358852 0.07872581 0.23452563 0.30457068 
##         15         16         17         18         19         20         21 
## 0.54270793 0.23639078 0.96570325 0.22184195 0.20738478 0.32434041 0.29233878 
##         22         23         24 
## 0.25769278 0.48547608 0.30917443

plot(hatvalues(modelo2))

Observaciones extremas respecto a x: valor de distancia

datos$vd<-hatvalues(modelo2)

which(datos$vd>=2*(5+1)/24)

## 15 17 
## 15 17

Las observaciones 15 y 17 son extremas respecto a x.

Observaciones extremas respecto a y: Residual estudentizado

datos$RS<-rstudent(modelo2)

which(abs(datos$RS)>=2)

## 24 
## 24

La observación 24 es extrema respecto a la variable de respuesta y.

Observaciones influyentes: Distancia cook

datos$DC<-cooks.distance(modelo2)

which(datos$DC>1)

## 17 
## 17

# Identificación y visualización de valores influyentes

influencePlot(modelo2)

##       StudRes        Hat      CookD
## 9   1.8080492 0.07718464 0.04046920
## 15  0.2648128 0.54270793 0.01462631
## 17 -0.9176812 0.96570325 3.98702338
## 24  2.3274683 0.30917443 0.32444747

summary(influence.measures(modelo2))

## Potentially influential observations of
##   lm(formula = y ~ x1 + x2 + x5 + x8 + x9, data = datos) :
## 
##    dfb.1_  dfb.x1  dfb.x2 dfb.x5  dfb.x8 dfb.x9 dffit   cov.r   cook.d  hat    
## 15  0.17   -0.19    0.19   0.11   -0.20   0.18   0.29    3.01_*  0.01    0.54  
## 17 -0.15    0.40   -0.21  -3.66_*  0.30  -0.06  -4.87_* 30.74_*  3.99_*  0.97_*
## 23 -0.01   -0.09    0.10   0.08    0.07   0.01  -0.17    2.71_*  0.01    0.49  
## 24 -1.02_*  1.08_* -0.32  -0.82    0.48  -0.73   1.56    0.39    0.32    0.31

# Cooks distance
plot(modelo2, 4, id.n =5)

La observación 17 es influyente y afecta las estimaciones del nuevo modelo.

Observaciones atípicas

# Identificación de los residuos estudentizados > 2 considerados como outliers

datos$studentized_residual <- rstudent(modelo2)

which(rstudent(modelo2) > 2)

## 24 
## 24

outlierTest(modelo2)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 24 2.327468           0.032549      0.78117

24 es una observación atípica con un residuo estandarizado por de bajo de 3.

Conclusión

El análisis muestra que la observación 17 es influyente y excede los límites para los valores de Leverages o Distancia Cook, por lo tanto dicho valor altera los resultado del analisis de regresión del nuevo modelo.
El modelo de regresión lineal multiple (modelo2), determina que el 86.57% de la variabilidad observada en la variable respuesta y (precio de venta de la casa/1000).
El test F del modelo2 muestra que es significativo (p-value: 2.899e-07).
Es recomendable realizar un modelo sin la observación 17 y observar si mejoran los resultados.