Regresión Lineal

La regresión implica el uso de una o más variables, consideradas como variables independientes, para predecir los valores de otra variable, la variable dependiente. Las variables que están fuertemente correlacionadas con la variable dependiente se usarán para predecir esa variable.

Primero, instalar y cargar los paquetes requeridos.

library(readr)
library(ggplot2)
library(corrplot) 
library(mlbench) 
library(Amelia) 
library(plotly) 
library(reshape2) 
library(caret)  
library(caTools) 
library(dplyr)

Base de datos Boston Housing

Los datos de la vivienda contienen 506 secciones censales de Boston del censo de 1970. La base de datos Boston Housing contiene los datos originales de Harrison y Rubinfeld (1979), el marco de datos BostonHousing 2 es la versión corregida con información espacial adicional.

Esta información esta incluida en la biblioteca mlbench o descargar el conjunto de datos. Los datos tienen las siguientes características, siendo medv la variable de objetivo o independiente:

  • 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 (Que tan lleno esta un lote)
  • 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
  • lstat - porcentaje de estado inferior de la población
  • medv - valor mediano de las viviendas ocupadas por sus propietarios en USD 1000

Data

Cargue los datos de BostonHousing y asígnelos a la carcasa variable.

data(BostonHousing)   #Para llamar datos de una biblioteca (como mlblench), usamos "data()"
housing <- BostonHousing  #Para usar la base de datos le asignamos a la variable "housing"
str(housing) # Obtenemos un data frame con la informacion de la base Bostonhousing
## '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 ...
head(housing)
##      crim zn indus chas   nox    rm  age    dis rad tax ptratio      b
## 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90
## 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90
## 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83
## 4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63
## 5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90
## 6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12
##   lstat medv
## 1  4.98 24.0
## 2  9.14 21.6
## 3  4.03 34.7
## 4  2.94 33.4
## 5  5.33 36.2
## 6  5.21 28.7
summary(housing) #Nos sirve para ver si hay cierta logica en los datos a tratar
##       crim                zn             indus       chas   
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   0:471  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1: 35  
##  Median : 0.25651   Median :  0.00   Median : 9.69          
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14          
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10          
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74          
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio            b         
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

Limpieza de datos

Luego tenemos que limpiar esta información. Hay muchas formas de hacerlo. Utilizaré missmap () del paquete de Amelia.

Compruebe si hay NA en el marco de datos.

missmap(housing,col=c('yellow','blue'),y.at=1,y.labels='',legend=TRUE) #

La gráfica anterior muestra claramente que los datos están libres de NA’s.

Análisis exploratorio de datos

Visualizaciones

Usemos ggplot2, corrplot y plotly para explorar los datos.

Correlación

Correlación y CorrPlots

La correlación se define como:

En las estadística, la dependencia o asociación es una relación estadística, causal o no, entre dos variables aleatorias o dos conjuntos de datos. La correlación es cualquiera de una amplia clase de relaciones estadísticas que involucran dependencia, aunque en el uso común a menudo se refiere a la medida en que dos variables tienen una relación lineal entre sí.

Los diagramas de correlación son una gran forma de explorar datos y ver si hay algún término de interacción.

corrplot(cor(select(housing,-chas))) #USAMOS LA libreria corrplot con la funcion "cor"

INTERPRETACIÓN: medv disminuye con hay aumento en crim (medio), indus (alto), nox (bajo), edad (bajo), rad (bajo), impuesto (bajo), ptratio (alto), lstat (alto) y aumenta con aumento en zn (bajo), rm (alto). “medv es la variale dependiente”.

Los colores fríos indican que si:

  • r=1 correlacion lineal directa

Los colores fríos indican que si:

  • r=-1 correlacion lineal inversa

Si no tiene color:

  • r=0 NO hay correlacion

Gráfico de densidad con ggplot

# visualizar la distribución de la variable indepeniente 

housing %>% 
  ggplot(aes(medv)) +
  stat_density() + 
  theme_bw()

Las visualizaciones anteriores revelan que las densidades máximas de medv están entre 15 y 30.

Gráfico de densidad con plotly

ggplotly(housing %>% 
  ggplot(aes(medv)) +
  stat_density() + 
  theme_bw())

Las visualizaciones anteriores revelan que las densidades máximas de medv están entre 15 y 30.

Efecto de las variables

Veamos el efecto de las variables en la base de datos en medv.

housing %>%
  select(c(crim, rm, age, rad, tax, lstat, medv,indus,nox,ptratio,zn)) %>%
  melt(id.vars = "medv") %>%
  ggplot(aes(x = value, y = medv, colour = variable)) +
  geom_point(alpha = 0.7) +
  stat_smooth(aes(colour = "black")) +
  facet_wrap(~variable, scales = "free", ncol = 2) +
  labs(x = "Variable Value", y = "Median House Price ($1000s)") +
  theme_minimal()

Los resultados del gráfico anterior están en correlación con el corrplot.

Construcción de modelos y predicción

El modelo de regresión lineal general en R:

Train Y Test Data

SIRVE PARA HACER UNA DE MANERA ALETORIA EL ENTRETAMIENTO(Calibracion) y eL TEST (Validacion del modelo) Permite dividir los datos en entrenamento y prueba de datos utilizando la biblioteca caTools. Lets split the data into train and test data using caTools library.

# establecer una semilla
set.seed(123) # Al usar siempre la misma semilla se genera un patron en la generacion de numeros aleatorios por eso es importante asignar una nueva semilla, con la misma de cantidad de digitos de los que se tiene en la base de datos.

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

split <- sample.split(housing,SplitRatio =0.75)

train <- subset(housing,split==TRUE)
test <- subset(housing,split==FALSE)

# train <- select(train,-b) #selccionar todos los datos excepto b.
# test <- select(test,-b)   #se quita el criterio personal del modelo y se deja todo aleactorio

Entrenando nuestro modelo

Vamos a construir nuestro modelo teniendo en cuenta que crim, rm, tax, lstat son los principales influyentes en la variable objetivo.

model <- lm(medv ~ -1 + rm + tax + lstat, data = train) #Con una variable es MRL univariante y con el resto de variables en multivariante
summary(model)             
## 
## Call:
## lm(formula = medv ~ -1 + rm + tax + lstat, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.148  -3.240  -1.262   2.192  29.844 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## rm     5.153675   0.101739  50.656  < 2e-16 ***
## tax   -0.008066   0.001940  -4.157 4.03e-05 ***
## lstat -0.534921   0.045159 -11.845  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.198 on 359 degrees of freedom
## Multiple R-squared:  0.9542, Adjusted R-squared:  0.9538 
## F-statistic:  2494 on 3 and 359 DF,  p-value: < 2.2e-16
#Para evitar el uso continuo de $ podemos incluir la data
#Se debe hacer un analisis de los residuales
#La variable Pr(>|t|) significa que: ¿Cual es la probabilidad que un a variable no sea estadistiacmente significativa para el modelo?
#En resumen mientras mayor sea la probabilidad no es significativa
#las mas significativas son los que tienen el Pr mas pequeño.
#El siguiente paso es extraer las variables que no son significativas, no retiramos todas en conjunto, si no uan a una desde la menos significativa.

Ecuacion del modelo encontrado

Vamos a copiar la ecuacion del modelo encontrado con la funcion model

## [1] "medv= 5.15 + -0.01 tax -0.53 lstat + NA NA"

Visualizando nuestro modelo

Permite visualizar nuestro modelo de regresión lineal trazando los residuos. La diferencia entre el valor observado de la variable dependiente (y) y el valor predicho (y) se denomina residual (e).

res <- residuals(model)

# Convertir residuos en un DataFrame
# 
res <- as.data.frame(res)

Un modelo es aceptable si: * Las variables son estadisticamente significativas * Si el R^2>=0.5 * Residuos que se ajusten a una N(0,1) distribucion normal de 0 a 1, (Varianza=1) * Residuos no correlacionados, deben estar diversos.

ggplot(res,aes(res)) +  geom_histogram(fill='blue',alpha=0.5)

plot(model)

#Los puntos con 373, 370, 371 son llamados valores atipicos,o que no hay una buena prediccion del modelo
#la segunda grafica muestra para que ajuste normal se adaptan mejor los datos
# el tecer grafico quita la escala y solo muestra la dispersion de los datos

CONCLUSION DE LOS RESULTADOS DEL MODELO

Predicciones

Probemos nuestro modelo prediciendo en nuestro conjunto de datos de prueba.

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

pl1 <-test %>% 
  ggplot(aes(medv,predicted.medv)) +
  geom_point(alpha=0.5) + 
  stat_smooth(aes(colour='black')) +
  xlab('Actual value of medv') +
  ylab('Predicted value of medv')+
  theme_bw()
#Calcula el intervalo de confianza 

ggplotly(pl1)

Evaluemos nuestro modelo

usando Root Mean Square Error, una medida estandarizada de cuán lejos estábamos con nuestros valores predichos.

error <- test$medv-test$predicted.medv  #Calculamos el error entre si es cercano a uno es un mododelo bueno
rmse <- sqrt(mean(error)^2)

El Root Mean Square Error (RMSE) para nuestro modelo es 0.7602882 y los resultados pueden mejorarse aún más utilizando la extracción de variables y entrenando el modelo.

Evaluando nuestro modelo