Regresión Lineal con el conjunto de datos de viviendas de Boston
CRISTIAN IVAN SOLORZANO MAYA
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.
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.
data(BostonHousing)#para llamar bases de datos incorporadas
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 ...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.7summary(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.
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 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í.
Nota:
El tamaño implica si hay o no correlacion
Colores frio: \(r\approx 1\), correlación lineal directa
Colores Calidos: \(r\approx -1\), correlación lineal inversa
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 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())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 Siempre se debe cambiar la semilla para los aleatorios
True INcluye en el trainig False Parte del test sample.split cuanto usaremos del training
# 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
Tenemos una lista de variables que podemos extraer del modelo,extraemos una a una desde la menos significativa.
Vamos a construir nuestro modelo teniendo en cuenta que crim, rm, tax, lstat son los principales influyentes en la variable objetivo.
lm lineal model Pr(>|t|): CUal es la probabilidad de que la variable crime no sea estadisticamente significativa para el modelo?
NOta: Extraeremos las variables significativas
#model <- lm(medv ~ rm + tax + lstat, data = train)
model <- lm(medv ~ -1+ rm + tax + lstat, data = train)
#-1 para quitar el intercepto
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-16Ecuacion 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).
Un modelo es acceptable si: -Las variables son estadísticamente significativas -\(R^2\) ajustado$$0.5 -Residuos ~N(0,1) residuos no correlacionados(No hay patron)
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)
Grafico 1: los ajustados Grafico 2: Ver si se ajustan a una normal GRafico 3: Residuos estandarizados Grafico 4: mide los datos dispersos
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()
ggplotly(pl1)valores predichos deben ser lo más cerca posible Interpolamos una posible curva
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)#cercano a uno es buenoEl 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