1 Introducción

require(AER)
require(tidyverse)
library(psych) #Para obtener los descriptivos

En este caso vamos a analizar el precio de casas en la ciudad de Windsor, Canada, entre los meses de Agosto y Septiembre de 1987. Los datos se encuentran en la librería AER en el dataset HousePrices.

El dataframe contiene 546 observaciones. La variable de interés es predecir el precio price.

data("HousePrices")
str(HousePrices)
## 'data.frame':    546 obs. of  12 variables:
##  $ price     : num  42000 38500 49500 60500 61000 66000 66000 69000 83800 88500 ...
##  $ lotsize   : num  5850 4000 3060 6650 6360 4160 3880 4160 4800 5500 ...
##  $ bedrooms  : num  3 2 3 3 2 3 3 3 3 3 ...
##  $ bathrooms : num  1 1 1 1 1 1 2 1 1 2 ...
##  $ stories   : num  2 1 1 2 1 1 2 3 1 4 ...
##  $ driveway  : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ recreation: Factor w/ 2 levels "no","yes": 1 1 1 2 1 2 1 1 2 2 ...
##  $ fullbase  : Factor w/ 2 levels "no","yes": 2 1 1 1 1 2 2 1 2 1 ...
##  $ gasheat   : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ aircon    : Factor w/ 2 levels "no","yes": 1 1 1 1 1 2 1 1 1 2 ...
##  $ garage    : num  1 0 0 0 0 0 2 0 0 1 ...
##  $ prefer    : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...

La descripción de las variableses:

  1. price: Sale price of a house.

  2. lotsize: Lot size of a property in square feet.

  3. bedrooms: Number of bedrooms.

  4. bathrooms: Number of full bathrooms.

  5. stories: Number of stories excluding basement.

  6. driveway: Factor. Does the house have a driveway?

  7. recreation: Factor. Does the house have a recreational room?

  8. fullbase: Factor. Does the house have a full finished basement?

  9. gasheat: Factor. Does the house use gas for hot water heating?

  10. aircon: Factor. Is there central air conditioning?

  11. garage: Number of garage places.

  12. prefer: Factor. Is the house located in the preferred neighborhood of the city?

2 Análisis Exploratorio

El primer paso es analizar la variable de respuesta y cómo se relaciona con otras variables:

#Density plot
HousePrices %>%
  ggplot(aes(x = price)) +
  geom_density(fill = 'red3', alpha = 0.5) +
  labs(title = 'Desidad del precio de casas en Windsor')

#Obtener los descriptivos del dataset
describe(HousePrices)

El precio promedio de de \(6,8121\) CA con una distribución que cláramente tiene un sesgo hacia la derecha.

2.1 Análisis Relacional

Vamos a comenzar con las variables continuas y la posible relación con el precio:

HousePrices %>%
  ggplot(aes(x = lotsize, y = price)) +
  geom_point() +
  geom_smooth() +
  labs(title = 'Relación precio vs lotsize')
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Acá se puede aprecar una relación positiva entre el tamaño del lote y el precio.

Aún cuando existen otras variables continuas en realidad bedrooms, bathrooms y stories pueden ser consideradas como factores:

#Transformando los datos
House_trf = HousePrices %>%
  mutate(bedrooms = factor(bedrooms),
         bathrooms = factor(bathrooms),
         stories = factor(stories),
         garage = factor(garage))

str(House_trf)
## 'data.frame':    546 obs. of  12 variables:
##  $ price     : num  42000 38500 49500 60500 61000 66000 66000 69000 83800 88500 ...
##  $ lotsize   : num  5850 4000 3060 6650 6360 4160 3880 4160 4800 5500 ...
##  $ bedrooms  : Factor w/ 6 levels "1","2","3","4",..: 3 2 3 3 2 3 3 3 3 3 ...
##  $ bathrooms : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 2 1 1 2 ...
##  $ stories   : Factor w/ 4 levels "1","2","3","4": 2 1 1 2 1 1 2 3 1 4 ...
##  $ driveway  : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ recreation: Factor w/ 2 levels "no","yes": 1 1 1 2 1 2 1 1 2 2 ...
##  $ fullbase  : Factor w/ 2 levels "no","yes": 2 1 1 1 1 2 2 1 2 1 ...
##  $ gasheat   : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ aircon    : Factor w/ 2 levels "no","yes": 1 1 1 1 1 2 1 1 1 2 ...
##  $ garage    : Factor w/ 4 levels "0","1","2","3": 2 1 1 1 1 1 3 1 1 2 ...
##  $ prefer    : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
#Análisis Exploratorio
House_trf[, -2] %>%
  gather(Var, Valor, 2 : ncol(.)) %>%
  ggplot(aes(x = Valor, y = price, fill = Var)) +
  geom_boxplot() +
  facet_wrap(~Var, scales = 'free') +
  theme(legend.position = 'none') +
  labs(title = 'Relación price vs cat.vars')
## Warning: attributes are not identical across measure variables;
## they will be dropped

Acá podemos observar lo siguiente:

  1. aircon: el precio tiene a ser más alto para casas con aire acondicionado

  2. bathrooms: a medida que la casa tiene más baños, el precio se incrementa

  3. bedrooms: a medida que el número de cuartos se incrementa, el precio se incrementa

  4. driveway: cuando la casa tiene entrada para los carros, el precio se incrementa

  5. fullbase: al parecer no hay relación

  6. garage: cuando se tienen más de 1 garage, el precio se incrementa

  7. gasheat: si la casa cuenta con calentador de gas, existe un ligero incremento en el precio

  8. recreation: si existe salón de juegos, el precio tiende a ser mayor

  9. stories: el precio se incrementa de forma significativa a partir de 3 o 4

3 Análisis de Regresión

Ahora vamos a plantear el modelo

3.1 Regresión Lineal

En este caso, el modelo quedaría como:

\(y_i = \beta_0 + \sum_{j=1}^m \beta_jx_{ij} +\epsilon_i\)

#Paso 1: Formula para el modelo ------
f = price ~ lotsize + bedrooms + bathrooms + stories + driveway +
  recreation + fullbase + gasheat + aircon + garage + prefer

#Paso2: Correr la regresion ------
m_lm = lm(f, data = House_trf)
summary(m_lm)
## 
## Call:
## lm(formula = f, data = House_trf)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -39281  -9195   -909   7213  75412 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22507.343  11000.467   2.046 0.041252 *  
## lotsize           3.456      0.360   9.600  < 2e-16 ***
## bedrooms2     -2190.663  11016.988  -0.199 0.842462    
## bedrooms3      1248.530  11053.813   0.113 0.910113    
## bedrooms4      2891.382  11204.557   0.258 0.796466    
## bedrooms5      3913.397  12184.006   0.321 0.748193    
## bedrooms6     11004.913  15566.674   0.707 0.479910    
## bathrooms2    13114.247   1756.304   7.467 3.45e-13 ***
## bathrooms3    29550.820   5157.685   5.729 1.70e-08 ***
## bathrooms4    76993.725  16352.076   4.708 3.20e-06 ***
## stories2       4856.200   1742.991   2.786 0.005527 ** 
## stories3      12505.681   2980.753   4.195 3.20e-05 ***
## stories4      19953.854   3124.777   6.386 3.77e-10 ***
## drivewayyes    6839.223   2068.830   3.306 0.001012 ** 
## recreationyes  4298.723   1909.197   2.252 0.024762 *  
## fullbaseyes    5779.706   1617.102   3.574 0.000384 ***
## gasheatyes    12574.597   3267.652   3.848 0.000134 ***
## airconyes     12335.351   1580.236   7.806 3.23e-14 ***
## garage1        6065.258   1715.040   3.537 0.000441 ***
## garage2        9477.444   1886.412   5.024 6.95e-07 ***
## garage3        2510.257   4851.950   0.517 0.605116    
## preferyes      9270.978   1708.003   5.428 8.73e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15410 on 524 degrees of freedom
## Multiple R-squared:  0.6797, Adjusted R-squared:  0.6668 
## F-statistic: 52.95 on 21 and 524 DF,  p-value: < 2.2e-16
#Pase 3: Siplificar el modelo ------
m_lm2 = step(m_lm)
## Start:  AIC=10551.65
## price ~ lotsize + bedrooms + bathrooms + stories + driveway + 
##     recreation + fullbase + gasheat + aircon + garage + prefer
## 
##              Df  Sum of Sq        RSS   AIC
## - bedrooms    5 1.2545e+09 1.2573e+11 10547
## <none>                     1.2448e+11 10552
## - recreation  1 1.2043e+09 1.2568e+11 10555
## - driveway    1 2.5961e+09 1.2707e+11 10561
## - fullbase    1 3.0346e+09 1.2751e+11 10563
## - gasheat     1 3.5179e+09 1.2800e+11 10565
## - garage      3 7.0010e+09 1.3148e+11 10576
## - prefer      1 6.9990e+09 1.3148e+11 10580
## - stories     3 1.1079e+10 1.3556e+11 10592
## - aircon      1 1.4475e+10 1.3895e+11 10610
## - bathrooms   3 2.2433e+10 1.4691e+11 10636
## - lotsize     1 2.1892e+10 1.4637e+11 10638
## 
## Step:  AIC=10547.12
## price ~ lotsize + bathrooms + stories + driveway + recreation + 
##     fullbase + gasheat + aircon + garage + prefer
## 
##              Df  Sum of Sq        RSS   AIC
## <none>                     1.2573e+11 10547
## - recreation  1 1.2016e+09 1.2693e+11 10550
## - driveway    1 2.2556e+09 1.2799e+11 10555
## - fullbase    1 3.4676e+09 1.2920e+11 10560
## - gasheat     1 3.5174e+09 1.2925e+11 10560
## - garage      3 7.4503e+09 1.3318e+11 10573
## - prefer      1 7.4505e+09 1.3318e+11 10577
## - aircon      1 1.4892e+10 1.4062e+11 10606
## - stories     3 1.6135e+10 1.4187e+11 10607
## - lotsize     1 2.4039e+10 1.4977e+11 10641
## - bathrooms   3 2.6056e+10 1.5179e+11 10644
summary(m_lm2)
## 
## Call:
## lm(formula = price ~ lotsize + bathrooms + stories + driveway + 
##     recreation + fullbase + gasheat + aircon + garage + prefer, 
##     data = House_trf)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -38352  -8988   -762   7511  75479 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.158e+04  2.346e+03   9.198  < 2e-16 ***
## lotsize       3.565e+00  3.545e-01  10.057  < 2e-16 ***
## bathrooms2    1.387e+04  1.702e+03   8.153 2.58e-15 ***
## bathrooms3    3.048e+04  5.071e+03   6.010 3.46e-09 ***
## bathrooms4    7.774e+04  1.630e+04   4.769 2.40e-06 ***
## stories2      6.830e+03  1.484e+03   4.603 5.22e-06 ***
## stories3      1.449e+04  2.812e+03   5.153 3.63e-07 ***
## stories4      2.165e+04  3.005e+03   7.206 2.00e-12 ***
## drivewayyes   6.271e+03  2.036e+03   3.081 0.002173 ** 
## recreationyes 4.292e+03  1.909e+03   2.248 0.024957 *  
## fullbaseyes   6.105e+03  1.598e+03   3.820 0.000150 ***
## gasheatyes    1.252e+04  3.254e+03   3.847 0.000134 ***
## airconyes     1.247e+04  1.576e+03   7.916 1.45e-14 ***
## garage1       5.948e+03  1.711e+03   3.477 0.000549 ***
## garage2       9.896e+03  1.871e+03   5.288 1.81e-07 ***
## garage3       3.348e+03  4.832e+03   0.693 0.488655    
## preferyes     9.489e+03  1.695e+03   5.599 3.47e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15420 on 529 degrees of freedom
## Multiple R-squared:  0.6764, Adjusted R-squared:  0.6667 
## F-statistic: 69.12 on 16 and 529 DF,  p-value: < 2.2e-16

En este punto hemos obtenido un modelo que explica una \(R^2=0.667\) y ahora vamos a revisar los supuesto de los residuales \(NID(0, \sigma^2)\)

residuals = m_lm2$residuals

#Análisis de residuales
par(mfrow = c(2,2))

#Normalidad
plot(density(residuals), col = 'red', main = 'Desityplot de residuales')  # density plot for 'speed'
polygon(density(residuals), col="red")

qqnorm(residuals)
qqline(residuals, col = 'red')

#Varianza Constante
plot(residuals, main = 'Constant Variance')

#Residuales independientes
acf(residuals)

3.2 Corrección de los Residuales

En particular, cuando revisamos los supuestos, notamos que los errores tienen una dispersión en la cola del lado derecho lo cual los aleja de la normalidad. Cómo podemos corregir esto?

En el módulo de introducción a la regresión lineal https://rpubs.com/blad0914/749636 vimo que podemos transformar la variable \(y\) aplicando \(ln(y)\) para reducir la varianza:

HousePrices %>%
  ggplot(aes(x = log(price))) +
  geom_density(fill = 'red3', alpha = 0.5) +
  labs(title = 'Desidad del precio de casas en Windsor')

lo cual nos aproxima a una distribución normal. Tambíen podemos aplicar \(ln(x_1)\) en donde \(x_1\) es el tamaño de la casa:

HousePrices %>%
  ggplot(aes(x = log(lotsize))) +
  geom_density(fill = 'red3', alpha = 0.5) +
  labs(title = 'Desidad del tamaño de casas en Windsor')

Ahora vamos a correr de nuevo el modelo:

#Paso 1: Formula para el modelo ------
f = log(price) ~ log(lotsize) + bedrooms + bathrooms + stories + driveway +
  recreation + fullbase + gasheat + aircon + garage + prefer

#Paso2: Correr la regresion ------
m_log = lm(f, data = House_trf)
summary(m_log)
## 
## Call:
## lm(formula = f, data = House_trf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68350 -0.12143  0.00641  0.13360  0.62481 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.031886   0.267275  30.051  < 2e-16 ***
## log(lotsize)   0.299964   0.027470  10.920  < 2e-16 ***
## bedrooms2      0.048683   0.150058   0.324 0.745746    
## bedrooms3      0.121459   0.150543   0.807 0.420143    
## bedrooms4      0.128758   0.152603   0.844 0.399197    
## bedrooms5      0.148377   0.165777   0.895 0.371176    
## bedrooms6      0.276717   0.212028   1.305 0.192434    
## bathrooms2     0.168812   0.023937   7.052 5.59e-12 ***
## bathrooms3     0.325218   0.070131   4.637 4.46e-06 ***
## bathrooms4     0.657382   0.222643   2.953 0.003292 ** 
## stories2       0.069899   0.023831   2.933 0.003503 ** 
## stories3       0.193185   0.040602   4.758 2.53e-06 ***
## stories4       0.258982   0.042555   6.086 2.24e-09 ***
## drivewayyes    0.110436   0.028528   3.871 0.000122 ***
## recreationyes  0.054707   0.026119   2.094 0.036696 *  
## fullbaseyes    0.105751   0.022041   4.798 2.09e-06 ***
## gasheatyes     0.166374   0.044507   3.738 0.000206 ***
## airconyes      0.160194   0.021617   7.410 5.08e-13 ***
## garage1        0.079680   0.023398   3.405 0.000711 ***
## garage2        0.115125   0.025671   4.485 8.98e-06 ***
## garage3       -0.004928   0.066157  -0.074 0.940655    
## preferyes      0.122935   0.023160   5.308 1.64e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2099 on 524 degrees of freedom
## Multiple R-squared:  0.6938, Adjusted R-squared:  0.6815 
## F-statistic: 56.53 on 21 and 524 DF,  p-value: < 2.2e-16
#Pase 3: Siplificar el modelo ------
m_log2 = step(m_log)
## Start:  AIC=-1683.04
## log(price) ~ log(lotsize) + bedrooms + bathrooms + stories + 
##     driveway + recreation + fullbase + gasheat + aircon + garage + 
##     prefer
## 
##                Df Sum of Sq    RSS     AIC
## <none>                      23.093 -1683.0
## - bedrooms      5    0.4322 23.526 -1682.9
## - recreation    1    0.1933 23.287 -1680.5
## - gasheat       1    0.6159 23.709 -1670.7
## - driveway      1    0.6604 23.754 -1669.6
## - garage        3    1.1147 24.208 -1663.3
## - fullbase      1    1.0145 24.108 -1661.6
## - prefer        1    1.2418 24.335 -1656.4
## - stories       3    2.0537 25.147 -1642.5
## - aircon        1    2.4202 25.514 -1630.6
## - bathrooms     3    3.0279 26.121 -1621.8
## - log(lotsize)  1    5.2549 28.348 -1573.1
summary(m_log2)
## 
## Call:
## lm(formula = log(price) ~ log(lotsize) + bedrooms + bathrooms + 
##     stories + driveway + recreation + fullbase + gasheat + aircon + 
##     garage + prefer, data = House_trf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68350 -0.12143  0.00641  0.13360  0.62481 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.031886   0.267275  30.051  < 2e-16 ***
## log(lotsize)   0.299964   0.027470  10.920  < 2e-16 ***
## bedrooms2      0.048683   0.150058   0.324 0.745746    
## bedrooms3      0.121459   0.150543   0.807 0.420143    
## bedrooms4      0.128758   0.152603   0.844 0.399197    
## bedrooms5      0.148377   0.165777   0.895 0.371176    
## bedrooms6      0.276717   0.212028   1.305 0.192434    
## bathrooms2     0.168812   0.023937   7.052 5.59e-12 ***
## bathrooms3     0.325218   0.070131   4.637 4.46e-06 ***
## bathrooms4     0.657382   0.222643   2.953 0.003292 ** 
## stories2       0.069899   0.023831   2.933 0.003503 ** 
## stories3       0.193185   0.040602   4.758 2.53e-06 ***
## stories4       0.258982   0.042555   6.086 2.24e-09 ***
## drivewayyes    0.110436   0.028528   3.871 0.000122 ***
## recreationyes  0.054707   0.026119   2.094 0.036696 *  
## fullbaseyes    0.105751   0.022041   4.798 2.09e-06 ***
## gasheatyes     0.166374   0.044507   3.738 0.000206 ***
## airconyes      0.160194   0.021617   7.410 5.08e-13 ***
## garage1        0.079680   0.023398   3.405 0.000711 ***
## garage2        0.115125   0.025671   4.485 8.98e-06 ***
## garage3       -0.004928   0.066157  -0.074 0.940655    
## preferyes      0.122935   0.023160   5.308 1.64e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2099 on 524 degrees of freedom
## Multiple R-squared:  0.6938, Adjusted R-squared:  0.6815 
## F-statistic: 56.53 on 21 and 524 DF,  p-value: < 2.2e-16
#Paso 4: Análisis de residuales

residuals = m_log2$residuals

#Análisis de residuales
par(mfrow = c(2,2))

#Normalidad
plot(density(residuals), col = 'red', main = 'Desityplot de residuales')  # density plot for 'speed'
polygon(density(residuals), col="red")

qqnorm(residuals)
qqline(residuals, col = 'red')

#Varianza Constante
plot(residuals, main = 'Constant Variance')

#Residuales independientes
acf(residuals)

Ahora, vemos que el supuesto de normalidad se ha corregido. Sin embargo, al parecer los errores están autocorrelacionados tal cual se muestra en la función de autocorrelación o ACF.

En el contexto del problema, la autocorrelación no hace sentido a menos que los datos hayan sido recolectados en forma secuenciales. En este caso, no sabemos el orden en el cual se han tomado las observaciones y es posible que esta autocorrelación se deba a un efecto raro o extraño.

Para verificar la independencia podemos graficar los residuales vs lotsize:

plot(y = residuals, x = log(House_trf$lotsize))

Lo que observamos, es una nube de puntos distribuidos en forma aleatoria lo cual confirma la independencia en los residuales.

Otra forma de verificar esto es reordenar los datos en forma aleatoria y hacer el análisis:

House_trf$residuals = residuals
House_trf$aleatorio = runif(nrow(House_trf))
#Reordenar
House_trf = House_trf[order(House_trf$aleatorio), ]


#Analisis de Residuales
acf(House_trf$residuals)

Ahora lo que observamos es que la autocorrelación desaparece corroborando lo anterior mencionado referente a que dado que los datos no han sido recolectados de manera secuencial en el tiempo, la autocorrelación en este caso no hace sentido.