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

#Forma de llamar a la base de datos
data(BostonHousing)
#Le llama a la base de datos
housing <- BostonHousing
#Muestra el tipo de variables que estamos manejando
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 ...
#Para ver que variables estamos trabajando
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
#Estadistica por cada una de las variables en función de porcentajes
#Proporciona una idea si los estados estan bien
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ón missmap () del paquete de Amelia.

Compruebe si hay NA en el marco de datos.

#Muestra en negro y amarillo segun los resultados
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íticas 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í.

La correlación que puede tener entre variables, como es el numero de hijos durante la guerra, y que producira esto

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

Hay una correlación directa y buena cuando es lineal, es decir se expresa en linea

#Forma de presentar la correlacion entre variables, la variable dependiente es medv mintras las demás son independientes
corrplot(cor(select(housing,-chas)))

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

Para ver en forma de matriz, lo asociamos a una matriz

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

CORRELACION Si r=1 Correlación lineal directa linea hacia arriba / y de colores fríos Si r=-1 Correlación lineal inversa linea hacia abajo y de colores cálidos Si r=0 no tiene correlacion lineal

El tamaño de la bolita nos indica el nivel de correlación.

Conocemos que variables tienen correlación y asi decidimos como hacer el modelamiento

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

El Test Data no se incluye y predice eventos que no hay ocurrido, es mas viable utilizar el Train

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, para ser aleatorio, el numero debe ser con la dimension de las variables y puede ser cualquiera
set.seed(123)

#Seccionar los datos , `split ()` asigna un booleano a una nueva columna basada en el SplitRatio especificado, le explica el porcentaje de train data.

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

#Extrae secciones de la base de datos
#Destina a que base van a trabajar en funcion del T y F
train <- subset(housing,split==TRUE)
test <- subset(housing,split==FALSE)

# train <- select(train,-b), el b es el excepto de la base de datos
# 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.

Lista de variables que podemos excluir del modelo de observación, no se retiran del todo, si no, una a una desde la menos significativa.

#Regresiones univariables y se van poniendo con el "+"
#Si tomamos el crim, regresion lineal
model <- lm(medv ~ crim + rm + tax + lstat, data = train)
summary(model)
## 
## Call:
## lm(formula = medv ~ crim + rm + tax + lstat, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.266  -3.185  -1.052   2.116  30.121 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.767079   3.573477  -1.054  0.29251    
## crim        -0.070793   0.037113  -1.908  0.05725 .  
## rm           5.580390   0.492854  11.323  < 2e-16 ***
## tax         -0.006392   0.002114  -3.023  0.00268 ** 
## lstat       -0.483836   0.058230  -8.309 2.04e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.183 on 357 degrees of freedom
## Multiple R-squared:  0.6816, Adjusted R-squared:  0.678 
## F-statistic: 191.1 on 4 and 357 DF,  p-value: < 2.2e-16

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]),
       round(model$coefficients[5],2),names(model$coefficients[5])
)
## [1] "medv=  -3.77 -0.07 crim + 5.58 rm -0.01 tax -0.48 lstat"

Para verlo como codigo simple usamos un apostrofe

La parte de arriba del estimate, nos da el coeficiente de la formula del modelo.

La ultima variable Pr(>t) nos indica lo viable de las variables y si tenemos valores, responde a la pregunta “cual es la probabilidad de que la variable crime NO sea estadisticcamente significativa para el modelo”, mientras mas alta es menos viables, mientras mas baja es mas viable

Se extrae a las variables que nos significativas, el intercepto queda para el final.

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)
#Para evidenciar la campana de Gauss en forma de distribucion normal
#Para realizar el diagrama utiliza la normalizacion de los datos
ggplot(res,aes(res)) +  geom_histogram(fill='blue',alpha=0.5)

Un modelo es aceptable si: Las variables son sifnificativas El r^2 ajustado a >= 0.5, residual standard Residuos que se ajustan a una distribución normal donde la media 0-1, se ve como una campana de Gauss reflejada en el histograma, que no estan relacionados, es decir que no sigan un patron deben estar dispersos. Visualmente se comprueba al asociar a una campana.

#La primera grafica es la diferencia entre la realidad y lo ajustado
#Los datos atipicos que se encuentran fuera del ajuste, generalmente son incorrectos
#Lo ideal es que sea una linea recta, ya que estamos viendo que tran cerca esta lo real de lo hecho.

#La grafica 2, son residuos estandarizados, si se suporne la distribucion teorica con la teorica

#La grafica 3, son residuos estandarizados quitan la distancia de la escala, para ver la naturaleza, para medir el ajuste de los datos, lo ideal es la linea recta
plot(model)

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

#ahora ocupamos el test, ya no el train
#la funcion predict
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 sirve para ver la tendecia
#theme_bw es para que muestre la cuadricula 
#ggplotl(pl1) es lo que potencializa y le da herramientas para manipular al observador
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
rmse <- sqrt(mean(error)^2)

El Root Mean Square Error (RMSE)es cercano a uno, significa que nuestro modelo es bueno y viable, 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