Regresión Lineal con el conjunto de datos de viviendas de Boston
Alexandra Maigua
26 de septiembre de 2018
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)El paquete corrplot es una visualización gráfica de una matriz de correlación, intervalo de confianza. También contiene algunos algoritmos para hacer la reordenación matricial. Además, corrplot es bueno en cuanto a los detalles, como elegir el color, las etiquetas de texto, las etiquetas de color, el diseño, etc.
El paquete mlbench es una colección de bases de datos del mundo real, que incluye, por ejemplo, varios conjuntos de datos del repositorio UCI.
El paquete Amelia es una herramienta que ayuda al tratamiento de datos faltantes en una sola sección transversal (como una encuesta), de una serie temporal (como variables recopiladas para cada año en un país), o de un conjunto de datos de series de tiempo transversales (como los recogidos por años para cada uno de varios países).
La biblioteca gráfica de Plotly hace gráficos interactivos de calidad de publicación en línea y ejemplos de cómo hacer gráficos de líneas, diagramas de dispersión, gráficos de áreas, gráficos de barras, barras de error, diagramas de cajas, histogramas, mapas de calor, subtramas, gráficos de ejes múltiples y gráficos 3D (basados en WebGL).
- El paquete caret (abreviatura de Classification and Regression Training) es un conjunto de funciones que intentan simplificar el proceso para crear modelos predictivos. El paquete contiene herramientas para:
- división de datos
- preprocesamiento
- selección de características
- ajuste de modelo usando remuestreo
- estimación de importancia variable
- así como otras funcionalidades.
Las librerias son actuales.
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,porcentaje de personas analfabetas y con trabajo informal
- medv - valor mediano de las viviendas ocupadas por sus propietarios en USD 1000
Se busca estimar la plusvalia de un barrio con esta base de datos.
Data
Cargue los datos de BostonHousing y asígnelos a la carcasa variable. Para llamar base de datos se utiliza data
data(BostonHousing)
housing <- BostonHousing
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 ...Se prenstan las primeras filas de la base de datos:
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.7El summary de esta base de datos es la siguiente:
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.00Limpieza 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, es decir la existencia de datos perdidos.
missmap(housing,col=c('orange','skyblue'),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í.La correlación es la medida más común entre variables.Correlación sera alta solamente si la relación es lineal
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.0000000medv disminuye cuando hay aumento en crim (medio), indus (alto), nox (bajo), edad (bajo), rad (bajo), impuesto (bajo), ptratio (alto), lstat (alto) y aumenta cuando aumenta 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())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:
Modelo univariado: modelo <- lm (y~x, datos)
Modelo multivariado: modelo <- lm (y~., Datos)
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. Test Data no se debe incluir dentro del modelo.
# establecer una semilla, para hacer uso de la parte aleatoria
set.seed(123)#cambiar la semilla para generación de aleatoreos
#Seccionar los datos , `split ()` asigna un booleano a una nueva columna basada en el SplitRatio especificado.
split <- sample.split(housing,SplitRatio =0.75)# el procentaje de train data
#mismas columnas con sus respectivas filas
train <- subset(housing,split==TRUE)#incluye en el train data
test <- subset(housing,split==FALSE)#parte del test
# train <- select(train,-b) excepto este dato
# 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.
La varibale Pr(>|t|) indica si la variable sea estadisticamente significativa para el modelo, mientras más pequeña sea la probabilidad más significativa es la variable. Asi tenemos una lista de variables que podemos extraer del modelo pero no debemos queitarlas todas en conjunto, si no una a una desde la menos significativa, debido a que nos puede generar un efecto en el resto de variables.
model <- lm(medv ~ -1+rm + tax + lstat, data = train)#especificamos que variables queremos regresiones multilineales, en la iz var dep y der indep, se especifica la data
#sin crim
#model1<-lm(medv~crim,data=train)
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#summary(model1)Ecuación del modelo encontrado
## [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). Para que un modelo sea aceptabble se debe cumpluir que: las variables son estadisticamente significativas R^2 ajustado es >= .5(variable ambigua) Residuo~N(0,1) Residuo no correlacionados es decir no deben seguir un patrón, deben estar dispersos.
res <- residuals(model)
# Convertir residuos en un DataFrame
# para poder graficar
res <- as.data.frame(res)ggplot(res,aes(res)) + geom_histogram(fill='magenta',alpha=0.3)#grafico de los residuos
plot(model)
Observamos la correlación los numeros mostrados son de datos atípicos en el modelo o el modelo no predice bien Ajuste cuantil: la linea recta representa la distribucion normal que se ajustaria a los datos, si los datos se sobreponen tenemos un buen modelo. Residuos estandarizados: quita la escala en los datos, para medir el ajuste de los datos. Generan curvas de nivel.
Predicciones
Probemos nuestro modelo prediciendo en nuestro conjunto de datos de prueba.
test$predicted.medv <- predict(model,test)#compara el real con el esperado
#usamos el 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()
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)#
rmse## [1] 0.7575575Si rmse es cercano a uno el modelo es bueno.
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