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
  • 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,personas que no tiene trabajo estable
  • 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 leer una base de datos se usa una libreria o se usa el comando data
housing <- BostonHousing # para 
str(housing)
## '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)
##       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','black'),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 (varias relaciones) 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í.Mide la correlacion lineal,es una correlacion buena

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

correlacion = corrplot(cor(select(housing,-chas)))

correlacion
##               crim         zn      indus        nox         rm        age
## crim     1.0000000 -0.2004692  0.4065834  0.4209717 -0.2192467  0.3527343
## zn      -0.2004692  1.0000000 -0.5338282 -0.5166037  0.3119906 -0.5695373
## indus    0.4065834 -0.5338282  1.0000000  0.7636514 -0.3916759  0.6447785
## nox      0.4209717 -0.5166037  0.7636514  1.0000000 -0.3021882  0.7314701
## rm      -0.2192467  0.3119906 -0.3916759 -0.3021882  1.0000000 -0.2402649
## age      0.3527343 -0.5695373  0.6447785  0.7314701 -0.2402649  1.0000000
## dis     -0.3796701  0.6644082 -0.7080270 -0.7692301  0.2052462 -0.7478805
## rad      0.6255051 -0.3119478  0.5951293  0.6114406 -0.2098467  0.4560225
## tax      0.5827643 -0.3145633  0.7207602  0.6680232 -0.2920478  0.5064556
## ptratio  0.2899456 -0.3916785  0.3832476  0.1889327 -0.3555015  0.2615150
## b       -0.3850639  0.1755203 -0.3569765 -0.3800506  0.1280686 -0.2735340
## lstat    0.4556215 -0.4129946  0.6037997  0.5908789 -0.6138083  0.6023385
## medv    -0.3883046  0.3604453 -0.4837252 -0.4273208  0.6953599 -0.3769546
##                dis        rad        tax    ptratio          b      lstat
## crim    -0.3796701  0.6255051  0.5827643  0.2899456 -0.3850639  0.4556215
## zn       0.6644082 -0.3119478 -0.3145633 -0.3916785  0.1755203 -0.4129946
## indus   -0.7080270  0.5951293  0.7207602  0.3832476 -0.3569765  0.6037997
## nox     -0.7692301  0.6114406  0.6680232  0.1889327 -0.3800506  0.5908789
## rm       0.2052462 -0.2098467 -0.2920478 -0.3555015  0.1280686 -0.6138083
## age     -0.7478805  0.4560225  0.5064556  0.2615150 -0.2735340  0.6023385
## dis      1.0000000 -0.4945879 -0.5344316 -0.2324705  0.2915117 -0.4969958
## rad     -0.4945879  1.0000000  0.9102282  0.4647412 -0.4444128  0.4886763
## tax     -0.5344316  0.9102282  1.0000000  0.4608530 -0.4418080  0.5439934
## ptratio -0.2324705  0.4647412  0.4608530  1.0000000 -0.1773833  0.3740443
## b        0.2915117 -0.4444128 -0.4418080 -0.1773833  1.0000000 -0.3660869
## lstat   -0.4969958  0.4886763  0.5439934  0.3740443 -0.3660869  1.0000000
## medv     0.2499287 -0.3816262 -0.4685359 -0.5077867  0.3334608 -0.7376627
##               medv
## crim    -0.3883046
## zn       0.3604453
## indus   -0.4837252
## nox     -0.4273208
## rm       0.6953599
## age     -0.3769546
## dis      0.2499287
## rad     -0.3816262
## tax     -0.4685359
## ptratio -0.5077867
## b        0.3334608
## lstat   -0.7376627
## medv     1.0000000

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).

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

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)

#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)
# test <- select(test,-b)

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)#lm(variable dependiete,variables independientes se separa con +)
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
#Pr(>|t|)probabilidad de que NO sea significativa 
#2 se puede extraer una lista de variables , no se retira todas en conjunto sino una a una desde la menos significativa

Ecuacion del modelo encontrado

paste("medv = ",round(model$coefficients[1],2),"",
      round(model$coefficients[2],2),names(model$coefficients[2]),"+",
      round(model$coefficients[3],2),names(model$coefficients[3]),
      round(model$coefficients[4],2),names(model$coefficients[4])
      
      )
## [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). Un modelo es estadisticamente aceptable Si las variables son significativas Si el R^2 ajustado sea mayor o igual a 0.5 Si el residua sigue una distribucion normal N(0,1) Si los residuos son no correlacionados, es decir, no siguen un patrón

res <- residuals(model)

# Convertir residuos en un DataFrame
# 
res <- as.data.frame(res)
ggplot(res,aes(res)) +  geom_histogram(fill='blue',alpha=0.5)

plot(model)

# 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()
#stat_smooth 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. Si rmse es sercano a 1 es un modelo bueno

error <- test$medv-test$predicted.medv
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