Clase 1: Ejemplo de una Regresión lineal
Copertino Quispe1
SUMMARY:
Analizaremos la base de datos Boston Hausing Dataset, la cual recoge el precio de las viviendas ubicadas en Boston, EE.UU. y factores asociados a estas. Cabe señalar que esta base de datos es una muestra del Censo de viviendas de 1970, y su uso es de fin académico.
1. Paquetes en R
A continuación realizaremos una aplicación práctica de una regresión lineal ; para ello utilizaremos los paquete “MASS”, el cual nos permitirá utilizar la base de datos y otras funciones, y “corrplot”, el cual nos ayudará con el gráfico de correlación, y el paquete “lattice”, el cual nos permitirá analizar un gráfico matricial entre el precio de vivienda y sus determinantes.
# Llamamos a los paquetes
library("MASS")
library("corrplot")
library("lattice")2. Cargando base de datos
Una vez instalado los paquetes, cargamos la base de datos, y luego analizaremos la misma.
# Cargamos la base de datos
boston <- data.frame(Boston)
# Visualización de las variables
names(boston) ## [1] "crim" "zn" "indus" "chas" "nox" "rm" "age"
## [8] "dis" "rad" "tax" "ptratio" "black" "lstat" "medv"La base de datos cuenta con 506 observaciones y 14 variables, siendo dos categóricas.
dim(boston)## [1] 506 143. Análisis exploratorio de datos
3.1 Resumen estadístico
Se puede observar como el valor de la media y mediana de algunas variables tienen una gran diferencia; por lo que, se presumiría que algunas variables cuentan con valores atípicos.
summary(boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## 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 black
## 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.003.2 Box plot
A continuación, se realizará un gráfico de cajas con la finalidad de obervar los datos atípicos; para ello se graficarán las variables continuas exceptuando las variables “chas” y “rad” ya que son variables categóricas.
Según los gráficos de cajas, las variables “crim”, “zn”, “rm” y “black” cuentan con valores atípicos.
dropList <- c('chas','rad')
boston_box <- boston[,!colnames(boston) %in% dropList]
# Colores
a <- rainbow(6)
b <- rainbow(6, alpha=0.2)
c <- rainbow(6, v=0.5)
# ---------------
# Grupo 1
par(mfrow = c(1, 4))
boxplot(boston_box$crim, main='crim',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$zn, main='zn',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$indus, main='indus',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$nox, main='nox',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
# ---------------
# Grupo 2
boxplot(boston_box$rm, main='rm',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$age, main='age',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$dis, main='dis',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$tax, main='tax',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
# ---------------
# Grupo 3
boxplot(boston_box$ptratio, main='ptratio',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$black, main='black',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$lstat, main='lstat',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
boxplot(boston_box$medv, main='medv',col=b, boxcol=c, whiskcol=a,
staplecol=c, outcol=c, outbg=c, pch=20, cex=1)
3.3 Análisis de correlación
Al analizar el gráfico de correlación se observa que las variables “lstat” y “medv” se encuentran altamente correlacionadas; mientras la asociación lineal con otras variables es baja.
# Correlación
corrplot(cor(boston), type="upper", tl.col="black", order="hclust", sig.level=0.05)
Al analizar las variables “lstat” y “medv”, se puede observar que existe una relación inversa con los precios; ademas que la función que mejor se ajusta sería una cuadrática.
splom(boston_box[,11:12],col = 'Sky Blue')
4. Regresión lineal
A continuación, se evaluará la relación lineal entre la variable precios (mdev) y su vinculación con los otros determinantes.
En el primer modelo se visualiza que las variables “indus” y “age”, no son significativas, por lo que pasaremos a retirarlas y “correr” otro modelo. De otro lado, el R2 del modelo es igual a 73.38%
mod1 <- lm(medv~.,data=boston)
summary(mod1)##
## Call:
## lm(formula = medv ~ ., data = boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.595 -2.730 -0.518 1.777 26.199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 ***
## crim -1.080e-01 3.286e-02 -3.287 0.001087 **
## zn 4.642e-02 1.373e-02 3.382 0.000778 ***
## indus 2.056e-02 6.150e-02 0.334 0.738288
## chas 2.687e+00 8.616e-01 3.118 0.001925 **
## nox -1.777e+01 3.820e+00 -4.651 4.25e-06 ***
## rm 3.810e+00 4.179e-01 9.116 < 2e-16 ***
## age 6.922e-04 1.321e-02 0.052 0.958229
## dis -1.476e+00 1.995e-01 -7.398 6.01e-13 ***
## rad 3.060e-01 6.635e-02 4.613 5.07e-06 ***
## tax -1.233e-02 3.760e-03 -3.280 0.001112 **
## ptratio -9.527e-01 1.308e-01 -7.283 1.31e-12 ***
## black 9.312e-03 2.686e-03 3.467 0.000573 ***
## lstat -5.248e-01 5.072e-02 -10.347 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
## F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16En el segundo modelo se incluye una función cuadrática de la variable “lstat”; siendo, el R2 del modelo igual a 78.16%. De otro lado, la variable “zn” no es significativa, por lo que esta situación se incluirá en el tercer modelo.
mod2 <- lm(medv~.-age-indus+I(lstat^2),data=boston)
summary(mod2)##
## Call:
## lm(formula = medv ~ . - age - indus + I(lstat^2), data = boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.5603 -2.6709 -0.3071 1.9469 25.0085
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.124187 4.632559 9.093 < 2e-16 ***
## crim -0.149072 0.030006 -4.968 9.34e-07 ***
## zn 0.021281 0.012500 1.703 0.089291 .
## chas 2.589821 0.775307 3.340 0.000900 ***
## nox -13.534522 3.229619 -4.191 3.30e-05 ***
## rm 3.233174 0.372802 8.673 < 2e-16 ***
## dis -1.357892 0.169051 -8.032 7.09e-15 ***
## rad 0.271744 0.057599 4.718 3.11e-06 ***
## tax -0.009546 0.003068 -3.111 0.001970 **
## ptratio -0.790820 0.118089 -6.697 5.82e-11 ***
## black 0.008174 0.002429 3.365 0.000824 ***
## lstat -1.669567 0.119012 -14.029 < 2e-16 ***
## I(lstat^2) 0.033055 0.003198 10.337 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.298 on 493 degrees of freedom
## Multiple R-squared: 0.7868, Adjusted R-squared: 0.7816
## F-statistic: 151.6 on 12 and 493 DF, p-value: < 2.2e-16En el tercer modelo se excluye la variable “zn” y se mantiene la forma funcional anterior, siendo el R2 del modelo igual a 78.08%. Así el modelo es significativa en todas sus variables.
mod3 <- lm(medv~.-age-indus-zn+I(lstat^2),data=boston)
summary(mod3)##
## Call:
## lm(formula = medv ~ . - age - indus - zn + I(lstat^2), data = boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.7575 -2.6846 -0.3668 1.9384 25.1250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.614879 4.632461 9.199 < 2e-16 ***
## crim -0.145772 0.030000 -4.859 1.59e-06 ***
## chas 2.582698 0.776784 3.325 0.000950 ***
## nox -13.976244 3.225361 -4.333 1.78e-05 ***
## rm 3.305519 0.371083 8.908 < 2e-16 ***
## dis -1.212833 0.146289 -8.291 1.07e-15 ***
## rad 0.264051 0.057532 4.590 5.64e-06 ***
## tax -0.008478 0.003009 -2.818 0.005032 **
## ptratio -0.853061 0.112503 -7.583 1.69e-13 ***
## black 0.008169 0.002434 3.357 0.000849 ***
## lstat -1.705182 0.117384 -14.527 < 2e-16 ***
## I(lstat^2) 0.034089 0.003145 10.838 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.306 on 494 degrees of freedom
## Multiple R-squared: 0.7855, Adjusted R-squared: 0.7808
## F-statistic: 164.5 on 11 and 494 DF, p-value: < 2.2e-16UNIVERSIDAD NACIONAL MAYOR DE SAN MARCOS, coperquispe@gmail.com↩︎