Regresión Lineal con el conjunto de datos de viviendas de Boston
Silvana Romo
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)
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.00Summary nos indica si los datos que se están usando es correcto. ## 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)
#legend=TRUE: la infografia a la derecha del grafico Missmap permite graficar los datos perdidos
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.
En el diagrama de correlación los colores fríos representan una correlacion lineal directa (r=1),los colores cálidos representan una correlacion lineal inversa(r= -1), y los colores tenues indican que no hay una correlacion lineal (r=0).
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#La correlacion es alta o buena solo cuando es lineal cuando se acerca a una recta siendo directa.
#mdv es dependiente y las demas independiente
#Crimen color rojo la corel es buenamedv 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.
# establecer una semilla
set.seed(123) # (123) Puede ir cualquier numero
#Seccionar los datos , `split ()` asigna un booleano a una nueva columna basada en el SplitRatio especificado.
split <- sample.split(housing,SplitRatio =0.75) #Se da la base de datos, y se especifica que porcentaje de train data se asigna.
train <- subset(housing,split==TRUE) # Si incluye en la muestra del Train data
#muestra las columnas de hpusing y solo filas para entrenar el modelo
test <- subset(housing,split==FALSE) # FALSE no incluye en el test data.
# 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)
#-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-16#Funcion lm: lineal model
#Varibles(depend ~ indep + indep)Ecuacion de modelos y predicción
## [1] "medv= 5.15 -0.01 tax + -0.53 lstat NA NA"Como analizar el summary de model: 1. La columna Pr(>|t|)=0.05725 indica cual es la probabilidad de que la variable crime no sea estadísticamente significativo para el modelo. Mientras mas pequeña es la probabilidad mas significativa es la variable.
2. Se tiene una lista de variables que podemos excluir del modelo, pero no se extraen en conjunto si no una a una de la menos significativa, ya que si quitamos una puede tener efecto en las demas.
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 aceptable si las variables son estadísticamente significativas
- \[{R}^2\] (mide la calidad del modelo) ajustado es >=0.5. 3.Los residuos N(0,1) ajusten a una distribucion normal, siendo media=0 y varianza=1, residuos no correlacionados.
- No siguen un patrón deben estar dispersos
- Cálculo de los residuos: \[res=\frac{res-\hat{res}}{\sigma}\]
res <- residuals(model)
#Se calculan los residuales.
# Convertir residuos en un DataFrame
#
res <- as.data.frame(res)ggplot(res,aes(res)) + geom_histogram(fill='blue',alpha=0.5)
Si los residuos siguen una distribucion N(0,1) tienen la forma de la campana de Gauss.
plot(model)
#Grafico 2: Linea reacta se ajusta a los residuos(distribucion teorica), circulo distribucion empírica.
#Gráfico 3: Residuos estandarizados: quitan la distancia de los datosy solo mide el ajuste de los datos, identificando la naturaleza de los residuos.
#
#Gráfico 4: Curvas de nivel: Pasan por los datos, miden la tendencia. Predicciones
Probemos nuestro modelo prediciendo en nuestro conjunto de datos de prueba.
test$predicted.medv <- predict(model,test) #Se usa el test data
pl1 <-test %>%
ggplot(aes(medv,predicted.medv)) +
geom_point(alpha=0.5) + #puntos, alpha mide el ancho de los puntos
stat_smooth(aes(colour='black')) + #Calcula el intervalo de confianza, parte ploma
xlab('Actual value of medv') +
ylab('Predicted value of medv')+
theme_bw() #Malla y color del margen
ggplotly(pl1) #Se aplica plotly a ggplot, para obtener la herramientas dinámicas. 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)Si RMSE 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