Regresión Lineal con el conjunto de datos de viviendas de Boston
Gabriela Borja
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 grafica 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 esta base de datos solo la llamamos
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) #las primeras lineas del codigo## 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','pink'),y.at=1,y.labels='',legend=TRUE) #primero:datos perdidos y segundo:datos que si hay
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.
corrplot(cor(select(housing,-chas)))
correlacion<-corrplot(cor(select(housing,-chas))) #para ver los valores de la graifca, calidos correlacion directa y frios inversa, tamaño si hay correlacion 
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). (variable dependiente)
Gráfico de densidad con ggplot
# visualizar la distribución de la variable indepeniente, solo las que hay correlacion
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) #ponemos una semilla para ver que tan aleatorio es, siempre cambiar la semila, debe ver de acuerod a los datos, 3 decimales=3, y asi.
#Seccionar los datos , `split ()` asigna un booleano a una nueva columna basada en el SplitRatio especificado.
#Como aplicar el modelo, tomamos un 80% para ver le modelo y el 20% medimos la precision del modelo, testdata, no en ajuste de modelo, predice eventos que nos ha ocurrido, el train tambien se usa.
split <- sample.split(housing,SplitRatio =0.75) #0.75 cuanto de traindata
train <- subset(housing,split==TRUE) #aleatoriamente estructurada
test <- subset(housing,split==FALSE) #aleatoriamente
# train <- select(train,-b), toda la data excepto estos porque es -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) #lineal model, variable dependediente y demas variables independientes, no es necesario ingresar a la base, especificamos la data. En pr(>|t|),¿cual es la probabilidad de que la variable crim no sea estadisticamente significativa para el modelo?: pr(>|t|) =0.05725, mas pequeña hay mas significancia de la varible y mas grande mas siginificativa.Intercepto al final se le deja, extraemos vairables no signiifcativas.No retiramos todas en conjunto, si no una a una desde la menos significativa.
# de aqui sale una formula
#pongo menos 1 y le quito 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
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).
res <- residuals(model) #model es como guardamos al lm
#un modelo es aceptable:
#las variables son estadisticamente significativas
#si R ajustado es >=0.5, sube cuandos mas variables hay
#si los residuos se ajustan a una distribucion normal de N(0,1), media 0, varianza 1
#residuos no correlacionales, no sigan un patron, deben estar dispersos
# Convertir residuos en un DataFrame
#
res <- as.data.frame(res)ggplot(res,aes(res)) + geom_histogram(fill='pink',alpha=0.5) #si los residuos tiene una N(0,1) debe parecerse a uan campana de gauss, los normaliza de los residuos, ¿Que datos se presenta?: res<-res-mediares/desviacionresiduos (normalizar datos)
plot(model) #se generan 4 graficas, que nos ayudan a ver que los residuos no esten correlacionados y sigan una distribucion normal, mas cercanos a cero, toca ir a buscar lo que esta mal, si tiene datos atipicos
#puntos distribucion empiricas, y los otros teoricos, si mejoramos le modelo no aparece
#raiz, quita la escala de los datos
#generan curvas de nivel que pasan por lo datosPredicciones
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) + #ancho de los puntos
stat_smooth(aes(colour='black')) + #calcula el intervalo de confianza, que tendencia sigues los datos
xlab('Actual value of medv') +
ylab('Predicted value of medv')+
theme_bw() #como se ve el grafico
ggplotly(pl1) #le aplicamos al archivo ggplot y tenemos todas la funciones de ggplotEvaluemos 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 #test, variable dependietne, test, predichos de la variable dependiente
rmse <- sqrt(mean(error)^2) #rmse modelo 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