Regresión Lineal con el conjunto de datos de viviendas de Boston
Andrea Conde
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)variable dependiente mdv
- 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,variable factor)
- 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 de una libreria, solo la estamos llamando
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.00#el summary me ayuda a ver si la base est{a bien}, es dcir si hay un sueldo de -10 hay que ocrregirLimpieza 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)
#me dice si hay datos perdidos
#legend=TRUE->para lo de la parte derecha con los cuadritosLa 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.
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 buena solo cuando es LINEAL, sin embargo hay otras formas de medir la correlación
#colores cálidos -> relación lineal inversa
#colores frios ->relacion lineal directamedv 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, como consejo se debe ir cmbaindo la semilla para no hacer un patron
set.seed(123)
#Seccionar los datos , `split ()` asigna un booleano a una nueva columna basada en el SplitRatio especificado, ES PARA SELECCIONAR DE MANERA ALEATORIA LAS OBSERVACIONES Y ASI PODER HACER EL ANALISIS DE LA BASE DE DATOS
split <- sample.split(housing,SplitRatio =0.75)#SplitRatio =0.75->que porcentaje vamos a usar de Trainning data
train <- subset(housing,split==TRUE) #subset->para dividir en los que tienen split=TRUE
test <- subset(housing,split==FALSE)
#-b, TODO EXCEPTO ESO
# train <- select(train,-b)
# test <- select(test,-b)Entrenando nuestro modelo
lista de variables que podemos excluir del modelo
observaion!! no retiramo stodas sino una a una , desde la mens significativa
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#Pr(>|t|) ->la probabilidad de que la variable CRIME no sea estadisticamente significativa para el modelo, mas pequeña->mas significativa
#AL INTERCEPTO SE LO DEJA AL FINAL, ES DECIR LUEGO DE QUITAR LAS VARIABLES QUE NO SON SIGNIFICATIVAS
#para la formula del modeloecuacion del modelo encontrado
#objeto espacial es decir una lists
model$coefficients## rm tax lstat
## 5.153674759 -0.008065944 -0.534921396paste("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]))## [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).
res <- residuals(model)
# Convertir residuos en un DataFrame
#
res <- as.data.frame(res)Un modelo es aceptable si:
1.Las variables son estadisticamente significativa. 2.R^2 ajustado es >=0.5 3.residuos ~N(0,1) 4.NO PRESENTREN PATRONES(DEBEN ESTAR DISPERSOS), ES DECIR NO DEBEN ESTAR CORRELACIONADOS ENTRE SI
ggplot(res,aes(res)) + geom_histogram(fill='blue',alpha=0.5)
plot(model)
#1°grafico(AJIUSTADOS vs REISUDOS)->las numeritos son las posiciones de ñas observaciones que no s ajustan==DATOS ATIPICOS
#3°RESIUDOS ESTANDARIZADOS(quitan la escala de la base de datos)
#4°no interesa mucho CURVAS DE NIVELPredicciones
Probemos nuestro modelo prediciendo en nuestro conjunto de datos de prueba.
#ajusta
test$predicted.medv <- predict(model,test)
#aes(medv,predicted.medv)->que va en el eje x y en el y
#stat_smooth(aes(colour='black'))->calcula el intervalo de confianza, que es la linea ploma
# theme_bw()->para que salga la cuarilla, el color del margen,etc
#ggplotly(pl1)->par que aparezca los valores con el mouse, y se le aplica a un objeto
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(valores reales) con nuestros valores predichos.
error <- test$medv-test$predicted.medv
rmse <- sqrt(mean(error)^2)#raiz media cuadratica
#rmse->si es cercano a uno entonces el modelo 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