Regresión Lineal Múltiple

Esta técnica estadística permite evaluar la relación entre una variable dependiente (numérica) y más de una variable predictora (tambien numérica). de modo que evalúa además el efecto de cada predictor en presencia del resto.

Para más información general de regresión vea https://rpubs.com/paraneda/regresion

La ecuaciónde regresión lineal múltiple se deriva de la regresión simple:

\[y = (\beta_{0} + \beta_{1} x_{1i} + \beta_{2} x_{2i} + \beta_{3} x_{3i} + ...+\beta_{n} x_{ni}) + \epsilon_{i}\]

# librerias a cargar
library(mlbench)  # dataset destinados a pruebas de machine learning
library(tidyverse)# 
library(corrplot) # 
library(plotly)   #
library(caret)    # creación de entrenamiento de clasificación y regresión
library(caTools)  # 
library(reshape2) #

Bondad de ajuste

Para determinar que tan bien se ajusta el modelo se utiliza \(R^2\) y el Error estándar residual (RSE)

Boston house

Registros del valor de una casa en los suburbios de Boston (USA) en base a diversas variables. Contiene 506 filas y 14 variables.

Fuente: https://www.kaggle.com/c/boston-housing

# carga de data
data(BostonHousing)
boston <- BostonHousing
str(boston) # mostrar la estructura de los datos.
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : num  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ b      : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

Parametros:

  • crim - Crimen per cápita por ciudad
  • zn - proporción de terrenos residenciales divididos en zonas para lotes de más de 25,000 pies cuadrados
  • indus - proporción de acres de negocios no minoristas por ciudad
  • chas - variable ficticia de Charles River (= 1 si el tramo limita el río, 0 de lo contrario)
  • nox - concentración de óxidos nítricos (partes por 10 millones)
  • rm - número promedio de habitaciones por vivienda
  • age - proporción de unidades ocupadas por sus propietarios construidas antes de 1940
  • dis - Distancias desproporcionadas a cinco centros de empleo de Boston
  • rad - índice de accesibilidad a las autopistas radiales
  • tax - tasa de impuesto a la propiedad de valor completo por USD 10,000
  • ptratio - colegios por localidad
  • b 1000 (B - 0,63)^ 2, donde B es la proporción de negros por ciudad (¿increíble no?)
  • lstat - porcentaje de estado inferior de la población
  • medv - valor mediano de las viviendas ocupadas por sus propietarios en USD 1000
#nombre de las variables
names(boston)
##  [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"    
##  [8] "dis"     "rad"     "tax"     "ptratio" "b"       "lstat"   "medv"
# temas pendientes...
# exploración de datos
# limpieza de datos

Coeficiente de correlación

La librería corrplot perite una visualización gráfica de la relación entre variables

# descartamos la variable chas por no ser numérica.
corrplot(round(cor(subset(boston, select = -chas)), digits = 3), type = "lower")

### Matriz de correlación

# se descartó variable "chas" por no ser numérica
round(cor(select(boston, -chas)), 3)
##           crim     zn  indus    nox     rm    age    dis    rad    tax ptratio
## crim     1.000 -0.200  0.407  0.421 -0.219  0.353 -0.380  0.626  0.583   0.290
## zn      -0.200  1.000 -0.534 -0.517  0.312 -0.570  0.664 -0.312 -0.315  -0.392
## indus    0.407 -0.534  1.000  0.764 -0.392  0.645 -0.708  0.595  0.721   0.383
## nox      0.421 -0.517  0.764  1.000 -0.302  0.731 -0.769  0.611  0.668   0.189
## rm      -0.219  0.312 -0.392 -0.302  1.000 -0.240  0.205 -0.210 -0.292  -0.356
## age      0.353 -0.570  0.645  0.731 -0.240  1.000 -0.748  0.456  0.506   0.262
## dis     -0.380  0.664 -0.708 -0.769  0.205 -0.748  1.000 -0.495 -0.534  -0.232
## rad      0.626 -0.312  0.595  0.611 -0.210  0.456 -0.495  1.000  0.910   0.465
## tax      0.583 -0.315  0.721  0.668 -0.292  0.506 -0.534  0.910  1.000   0.461
## ptratio  0.290 -0.392  0.383  0.189 -0.356  0.262 -0.232  0.465  0.461   1.000
## b       -0.385  0.176 -0.357 -0.380  0.128 -0.274  0.292 -0.444 -0.442  -0.177
## lstat    0.456 -0.413  0.604  0.591 -0.614  0.602 -0.497  0.489  0.544   0.374
## medv    -0.388  0.360 -0.484 -0.427  0.695 -0.377  0.250 -0.382 -0.469  -0.508
##              b  lstat   medv
## crim    -0.385  0.456 -0.388
## zn       0.176 -0.413  0.360
## indus   -0.357  0.604 -0.484
## nox     -0.380  0.591 -0.427
## rm       0.128 -0.614  0.695
## age     -0.274  0.602 -0.377
## dis      0.292 -0.497  0.250
## rad     -0.444  0.489 -0.382
## tax     -0.442  0.544 -0.469
## ptratio -0.177  0.374 -0.508
## b        1.000 -0.366  0.333
## lstat   -0.366  1.000 -0.738
## medv     0.333 -0.738  1.000

Modelo de regresión

# determinación de modelo con más variables aplicadas
modelo <- lm(medv ~ rm + zn + dis +b, data = boston)
summary(modelo)
## 
## Call:
## lm(formula = medv ~ rm + zn + dis + b, data = boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.768  -2.877  -0.264   2.407  37.908 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -37.446075   2.707528 -13.830  < 2e-16 ***
## rm            8.194770   0.409907  19.992  < 2e-16 ***
## zn            0.062711   0.016134   3.887 0.000115 ***
## dis          -0.238285   0.178453  -1.335 0.182391    
## b             0.024306   0.003134   7.755 4.98e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.129 on 501 degrees of freedom
## Multiple R-squared:  0.5595, Adjusted R-squared:  0.5559 
## F-statistic: 159.1 on 4 and 501 DF,  p-value: < 2.2e-16
modelo2 <- lm(medv ~ rm + zn + dis +b + crim + indus + nox, data = boston)
summary(modelo2)
## 
## Call:
## lm(formula = medv ~ rm + zn + dis + b + crim + indus + nox, data = boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.188  -3.223  -0.534   2.155  36.389 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -10.165095   4.311853  -2.357 0.018785 *  
## rm            7.015433   0.408120  17.190  < 2e-16 ***
## zn            0.071032   0.015150   4.688 3.56e-06 ***
## dis          -1.497164   0.227595  -6.578 1.21e-10 ***
## b             0.014106   0.003135   4.499 8.52e-06 ***
## crim         -0.177686   0.034430  -5.161 3.56e-07 ***
## indus        -0.205569   0.063313  -3.247 0.001245 ** 
## nox         -15.533971   4.010253  -3.874 0.000122 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.702 on 498 degrees of freedom
## Multiple R-squared:  0.621,  Adjusted R-squared:  0.6156 
## F-statistic: 116.5 on 7 and 498 DF,  p-value: < 2.2e-16
# Cómo se distribuye la variable medv
ggplotly(boston %>% 
  ggplot(aes(medv)) +
  stat_density() + 
  theme_bw())

Evaluación del modelo

Establecer una partición de los datos de manera de utilizar una parte para entrenar el valor de predicción y el resto para testear el modelo.

# establecer una semilla
set.seed(123) # para generar num aleatorios , dos decimales seteo dos decimales .. ->ajustar datos

#Seccionar los datos , `split ()` asigna un booleano a una nueva columna basada en el SplitRatio especificado.

split <- sample.split(boston, SplitRatio =0.8)

train <- subset(boston,split==TRUE)#extrae secciones de db housing
test <- subset(boston,split==FALSE)

#modelo
model <- lm(medv ~ -1 + rm + zn + dis + b, data = train) #especifico data para no usar $
summary(model)
## 
## Call:
## lm(formula = medv ~ -1 + rm + zn + dis + b, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.2262  -4.7026  -0.7294   2.1897  28.8332 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## rm   3.060724   0.228885  13.372  < 2e-16 ***
## zn   0.128693   0.020315   6.335 6.49e-10 ***
## dis -0.440750   0.228170  -1.932  0.05412 .  
## b    0.010587   0.004022   2.632  0.00882 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.047 on 394 degrees of freedom
## Multiple R-squared:  0.9159, Adjusted R-squared:  0.915 
## F-statistic:  1073 on 4 and 394 DF,  p-value: < 2.2e-16

Un valor R-squared: 0.9159 indica una alta relación del modelo de regresión.

Visualizando el modelo

res <- residuals(model) 
# Convertir residuos en un DataFrame
# 
res <- as.data.frame(res)
#un modelo es aceptable si las var son significativas
#r^2 ajustado >=0.5
#residuos se ajusta distr normal media=0,varianza=1...-> N)(0,1)
#residuos no corelacionados (No siguen 1 patró, deben estar dispersos)
ggplot(res,aes(res)) +  geom_histogram(fill='blue',alpha=0.5)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Modelo prediccionando datos de prueba.

test$predicted.medv <- predict(model,test)

pl1 <-test %>% 
  ggplot(aes(medv,predicted.medv)) +
  geom_point(alpha=0.7) + 
  geom_smooth(aes(colour='black')) +
  xlab('Actual value of medv') +
  ylab('Predicted value of medv')+
  theme_bw()

ggplotly(pl1)

Evaluando modelo

Para analizar la confiabilidad de un modelo de regresión, además de la determinación del error medio (RMSE) también se realiza un análisis de residuos.

error <- test$medv - test$predicted.medv
rmse <- sqrt(mean(error)^2)
rmse
## [1] 0.1491127

Análisis de residuos

Supuesto 1: Los errores deben seguir una distribución normal

residuos <- rstandard(model) # residuos estándares del modelo ajustado (completo) 
par(mfrow = c(1,3)) # divide la ventana en una fila y tres columnas 
hist(residuos) # histograma de los residuos estandarizados 
boxplot(residuos) # diagrama de cajas de los residuos estandarizados 
qqnorm(residuos) # gráfico de cuantiles de los residuos estandarizados 
qqline(residuos)

Supuesto 2: La varianza de los errores es constante

plot(fitted.values(model), rstandard(model), xlab="Valores ajustados", ylab="Residuos estandarizados")  # gráfico 2D de los valores ajustados vs. los residuos estandarizados 
abline(h=0) # dibuja la recta en cero