House prices: Lasso, XGBoost, and a detailed EDA
Jesús Alba
#Introduccion
Para el trabajo de Open Data he elegido estos datos de Kaggle sobre el precio de la vivienda. Dejo el link abajo al código original.
He elegido estos datos por una situación familiar en la que estamos buscando comprarnos una casa. Puesto que es prioridad para mi este tema en estos momentos me ha parecido curioso.
Contestar a las preguntas de tipo: ¿Cuáles son las variables que se tienen en cuenta para el precio? ¿Cuáles son más importantes o tienen más peso?
1 Cargar y explorar los datos
##Cargar las librerias necesarias en R
library(knitr)
library(ggplot2)
library(plyr)
library(dplyr)
library(corrplot)
library(caret)
library(gridExtra)
library(scales)
library(Rmisc)
library(ggrepel)
library(randomForest)
library(psych)
library(xgboost)Leemos los datos tanto de test como de entrenamiento.
train <- read.csv("train.csv", stringsAsFactors = F)
test <- read.csv("test.csv", stringsAsFactors = F)##Tamaño de los datos y estructura.
El conjunto de datos del tren consiste en variables de caracteres y enteras. La mayoría de las variables de caracteres son en realidad factores (ordinales), pero he optado por leerlas en R como cadenas de caracteres, ya que la mayoría de ellas requieren una limpieza y/o una ingeniería de características primero. En total, hay 81 columnas/variables, de las cuales la última es la variable de respuesta (SalePrice). A continuación, sólo muestro una parte de las variables. Todas ellas se comentan con más detalle a lo largo del documento.
dim(train)## [1] 1460 81str(train[,c(1:10, 81)]) #mostrar las 10 primeras variables y la variable de respuesta## 'data.frame': 1460 obs. of 11 variables:
## $ Id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ MSSubClass : int 60 20 60 70 60 50 20 60 50 190 ...
## $ MSZoning : chr "RL" "RL" "RL" "RL" ...
## $ LotFrontage: int 65 80 68 60 84 85 75 NA 51 50 ...
## $ LotArea : int 8450 9600 11250 9550 14260 14115 10084 10382 6120 7420 ...
## $ Street : chr "Pave" "Pave" "Pave" "Pave" ...
## $ Alley : chr NA NA NA NA ...
## $ LotShape : chr "Reg" "Reg" "IR1" "IR1" ...
## $ LandContour: chr "Lvl" "Lvl" "Lvl" "Lvl" ...
## $ Utilities : chr "AllPub" "AllPub" "AllPub" "AllPub" ...
## $ SalePrice : int 208500 181500 223500 140000 250000 143000 307000 200000 129900 118000 ...#Deshacerse de los IDs pero mantener los IDs de las pruebas en un vector. Estos son necesarios para componer el archivo de presentación
test_labels <- test$Id
test$Id <- NULL
train$Id <- NULLtest$SalePrice <- NA
all <- rbind(train, test)
dim(all)## [1] 2919 80Sin los IDs el dataset tiene 79 predictores y una sola variable de salida, el precio de la vivienda.
#Explorando las más importante
##La variable respuesta; SalePrice
Como se puede ver, los precios de venta están muy sesgados. Esto era de esperar, ya que pocas personas pueden permitirse casas muy caras. Esto, el autor, lo toma en cuenta a la hora de modelar.
ggplot(data=all[!is.na(all$SalePrice),], aes(x=SalePrice)) +
geom_histogram(fill="blue", binwidth = 10000) +
scale_x_continuous(breaks= seq(0, 800000, by=100000), labels = comma)
summary(all$SalePrice)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 34900 129975 163000 180921 214000 755000 1459##Los predictores numéricos más importantes
Las variables de carácter necesitan algo de trabajo antes de que pueda utilizarlas. Para hacerme una idea del conjunto de datos, decidí ver primero qué variables numéricas tienen una alta correlación con el Precio de venta.
###Correlación con SalePrice
En total, hay 10 variables numéricas con una correlación de al menos 0,5 con SalePrice. Todas esas correlaciones son positivas.
numericVars <- which(sapply(all, is.numeric)) #index vector numeric variables
numericVarNames <- names(numericVars) #saving names vector for use later on
cat('There are', length(numericVars), 'numeric variables')## There are 37 numeric variablesall_numVar <- all[, numericVars]
cor_numVar <- cor(all_numVar, use="pairwise.complete.obs") #correlations of all numeric variables
#sort on decreasing correlations with SalePrice
cor_sorted <- as.matrix(sort(cor_numVar[,'SalePrice'], decreasing = TRUE))
#select only high corelations
CorHigh <- names(which(apply(cor_sorted, 1, function(x) abs(x)>0.5)))
cor_numVar <- cor_numVar[CorHigh, CorHigh]
corrplot.mixed(cor_numVar, tl.col="black", tl.pos = "lt")
En el resto de esta sección, visualizaré la relación entre el Precio de Venta y los dos predictores con la mayor correlación con el Precio de Venta; la Calidad General y la Superficie Habitable “Por Encima del Grado” (esta es la proporción de la casa que no está en un sótano; enlace).
También queda claro que la multicolinealidad es un problema. Por ejemplo: la correlación entre GarageCars y GarageArea es muy alta (0,89), y ambas tienen correlaciones similares (altas) con SalePrice. Las otras seis variables con una correlación superior a 0,5 con Precio de venta son: -TotalBsmtSF: Total de pies cuadrados de superficie de sótano -1stFlrSF: Pies cuadrados de la primera planta -FullBath: Baños completos por encima del grado -TotRmsAbvGrd: Total de habitaciones sobre el grado (no incluye los baños) -YearBuilt: Fecha de construcción original -YearRemodAdd: Fecha de remodelación (la misma que la fecha de construcción si no hay remodelaciones o adiciones)
###Calidad general
La calidad general tiene la mayor correlación con el precio de venta entre las variables numéricas (0,79). Califica el material y el acabado general de la casa en una escala de 1 (muy pobre) a 10 (muy excelente).
ggplot(data=all[!is.na(all$SalePrice),], aes(x=factor(OverallQual), y=SalePrice))+
geom_boxplot(col='blue') + labs(x='Overall Quality') +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma)
La correlación positiva está ciertamente ahí, y parece ser una curva ligeramente ascendente. En cuanto a los valores atípicos, no veo ningún valor extremo. Si hay un candidato para sacar como valor atípico más adelante, parece ser la casa cara con grado 4.
###Above Grade (Ground) Living Area (square feet)
La variable numérica con la segunda mayor correlación con el precio de venta es Above Grade Living Area. Esto tiene mucho sentido; las casas grandes suelen ser más caras.
ggplot(data=all[!is.na(all$SalePrice),], aes(x=GrLivArea, y=SalePrice))+
geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +
geom_text_repel(aes(label = ifelse(all$GrLivArea[!is.na(all$SalePrice)]>4500, rownames(all), '')))## `geom_smooth()` using formula 'y ~ x'
Especialmente las dos casas con grandes superficies y bajos precios de venta parecen valores atípicos (casas 524 y 1299, ver etiquetas en el gráfico). No los eliminaré todavía, ya que tomar valores atípicos puede ser peligroso. Por ejemplo, una puntuación baja en la calidad general podría explicar un precio bajo. Sin embargo, como se puede ver a continuación, estas dos casas también obtienen la máxima puntuación en la calidad general. Por lo tanto, tendré en cuenta las casas 1299 y 524 como principales candidatas a ser eliminadas como valores atípicos.
all[c(524, 1299), c('SalePrice', 'GrLivArea', 'OverallQual')]## SalePrice GrLivArea OverallQual
## 524 184750 4676 10
## 1299 160000 5642 10#Datos perdidos, codificación de etiquetas y factorización de variables
##Integridad de los datos
En primer lugar, me gustaría ver qué variables contienen valores perdidos.
NAcol <- which(colSums(is.na(all)) > 0)
sort(colSums(sapply(all[NAcol], is.na)), decreasing = TRUE)## PoolQC MiscFeature Alley Fence SalePrice FireplaceQu
## 2909 2814 2721 2348 1459 1420
## LotFrontage GarageYrBlt GarageFinish GarageQual GarageCond GarageType
## 486 159 159 159 159 157
## BsmtCond BsmtExposure BsmtQual BsmtFinType2 BsmtFinType1 MasVnrType
## 82 82 81 80 79 24
## MasVnrArea MSZoning Utilities BsmtFullBath BsmtHalfBath Functional
## 23 4 2 2 2 2
## Exterior1st Exterior2nd BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF
## 1 1 1 1 1 1
## Electrical KitchenQual GarageCars GarageArea SaleType
## 1 1 1 1 1cat('There are', length(NAcol), 'columns with missing values')## There are 35 columns with missing valuesPor supuesto, los 1459 NAs en SalePrice coinciden perfectamente con el tamaño del conjunto de pruebas. Esto significa que tengo que arreglar los NAs en 34 variables predictoras.
##Imputar los datos que faltan {.tabset}
En esta sección, voy a corregir los 34 predictores que contienen valores perdidos. Iré a través de ellos trabajando hacia abajo desde la mayoría de los NA hasta que los haya arreglado todos. Si me encuentro con una variable que realmente forma un grupo con otras variables, también las trataré como un grupo. Por ejemplo, hay varias variables relacionadas con Piscina, Garaje y Sótano.
Como quiero mantener el documento lo más legible posible, he decidido utilizar la opción de “Tabs” que proporciona knitr. Puedes encontrar un breve análisis para cada (grupo de) variables bajo cada pestaña. No es necesario recorrer todas las secciones, y también se puede echar un vistazo a unas pocas pestañas. Si lo haces, creo que especialmente las secciones de Garaje y Sótano merecen la pena, ya que he sido cuidadoso a la hora de determinar qué casas realmente no tienen sótano o garaje.
Además de asegurarme de que se eliminan las NA, también he convertido las variables de carácter en enteros ordinales si hay una clara ordinalidad, o en factores si los niveles son categorías sin ordinalidad. Convertiré estos factores en numéricos más adelante utilizando one-hot encoding (usando la función model.matrix).
###Variables de la piscina
Calidad de la piscina y la variable PoolArea
El PoolQC es la variable con más NAs. La descripción es la siguiente:
PoolQC: Pool quality
Ex Excellent
Gd Good
TA Average/Typical
Fa Fair
NA No Pool
Por lo tanto, es obvio que tengo que asignar simplemente “No pool” a las NA. Además, el elevado número de NA tiene sentido, ya que normalmente sólo una pequeña proporción de casas tiene piscina.
all$PoolQC[is.na(all$PoolQC)] <- 'None'También está claro que puedo etiquetar esta variable ya que los valores son ordinales. Como hay múltiples variables que utilizan los mismos niveles de calidad, voy a crear un vector que puedo reutilizar más adelante.
Qualities <- c('None' = 0, 'Po' = 1, 'Fa' = 2, 'TA' = 3, 'Gd' = 4, 'Ex' = 5)all$PoolQC<-as.integer(revalue(all$PoolQC, Qualities))
table(all$PoolQC)##
## 0 2 4 5
## 2909 2 4 4Sin embargo, hay una segunda variable que se relaciona con las Piscinas. Se trata de la variable PoolArea (en pies cuadrados). Como puede ver a continuación, hay 3 casas sin PoolQC. En primer lugar, comprobé si había una relación clara entre PoolArea y PoolQC. Como no vi una relación clara (piscinas más grandes o más pequeñas con mejor PoolQC), voy a imputar los valores de PoolQC basándome en la Calidad General de las casas (que no es muy alta para esas 3 casas).
all[all$PoolArea>0 & all$PoolQC==0, c('PoolArea', 'PoolQC', 'OverallQual')]## PoolArea PoolQC OverallQual
## 2421 368 0 4
## 2504 444 0 6
## 2600 561 0 3all$PoolQC[2421] <- 2
all$PoolQC[2504] <- 3
all$PoolQC[2600] <- 2###Miscellaneous Feature
Miscellaneous feature no incluídos en otras categorías
Con Miscellaneous Feature, hay 2814 NAs. Como los valores no son ordinales, convertiré MiscFeature en un factor. Valores:
Elev Elevator
Gar2 2nd Garage (if not described in garage section)
Othr Other
Shed Shed (over 100 SF)
TenC Tennis Court
NA Noneall$MiscFeature[is.na(all$MiscFeature)] <- 'None'
all$MiscFeature <- as.factor(all$MiscFeature)
ggplot(all[!is.na(all$SalePrice),], aes(x=MiscFeature, y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue') +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..))## Warning: Ignoring unknown parameters: fun.y## No summary function supplied, defaulting to `mean_se()`
table(all$MiscFeature)##
## Gar2 None Othr Shed TenC
## 5 2814 4 95 1Al observar las frecuencias, la variable me parece irrelevante. Tener un cobertizo probablemente signifique “no tener garaje”, lo que explicaría el menor precio de venta del cobertizo. Además, aunque tiene mucho sentido que una casa con pista de tenis sea cara, sólo hay una casa con pista de tenis en el conjunto de entrenamiento.
###Callejón
Tipo de callejón de acceso a la propiedad
Con Alley, hay 2721 NAs. Como los valores no son ordinales, convertiré a Alley en un factor. Valores:
Grvl Gravel
Pave Paved
NA No alley accessall$Alley[is.na(all$Alley)] <- 'None'
all$Alley <- as.factor(all$Alley)
ggplot(all[!is.na(all$SalePrice),], aes(x=Alley, y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue')+
scale_y_continuous(breaks= seq(0, 200000, by=50000), labels = comma)## Warning: Ignoring unknown parameters: fun.y## No summary function supplied, defaulting to `mean_se()`
table(all$Alley)##
## Grvl None Pave
## 120 2721 78###Valla
Calidad de la valla
Con Fence, hay 2348 NAs. Los valores parecen ser ordinales. Valores:
GdPrv Good Privacy
MnPrv Minimum Privacy
GdWo Good Wood
MnWw Minimum Wood/Wire
NA No Fenceall$Fence[is.na(all$Fence)] <- 'None'
table(all$Fence)##
## GdPrv GdWo MnPrv MnWw None
## 118 112 329 12 2348all[!is.na(all$SalePrice),] %>% group_by(Fence) %>% summarise(median = median(SalePrice), counts=n())## # A tibble: 5 x 3
## Fence median counts
## <chr> <dbl> <int>
## 1 GdPrv 167500 59
## 2 GdWo 138750 54
## 3 MnPrv 137450 157
## 4 MnWw 130000 11
## 5 None 173000 1179Mi conclusión es que los valores no parecen ordinales (ninguna valla es mejor). Por lo tanto, voy a convertir Fence en un factor.
all$Fence <- as.factor(all$Fence)###Variables de la chimenea
Calidad de la chimenea, y Número de chimeneas
Dentro de la calidad de las chimeneas, hay 1420 NAs. El número de chimeneas está completo.
Fireplace quality
El número de NAs en FireplaceQu coincide con el número de casas con 0 chimeneas. Esto significa que puedo sustituir con seguridad los NAs en FireplaceQu por “sin chimenea”. Los valores son ordinales, y puedo utilizar el vector de Cualidades que ya he creado para la Cualidad Piscina. Valores:
Ex Excellent - Exceptional Masonry Fireplace
Gd Good - Masonry Fireplace in main level
TA Average - Prefabricated Fireplace in main living area or Masonry Fireplace in basement
Fa Fair - Prefabricated Fireplace in basement
Po Poor - Ben Franklin Stove
NA No Fireplaceall$FireplaceQu[is.na(all$FireplaceQu)] <- 'None'
all$FireplaceQu<-as.integer(revalue(all$FireplaceQu, Qualities))
table(all$FireplaceQu)##
## 0 1 2 3 4 5
## 1420 46 74 592 744 43Número de chimeneas
Chimeneas es una variable entera, y no hay valores perdidos.
table(all$Fireplaces)##
## 0 1 2 3 4
## 1420 1268 219 11 1sum(table(all$Fireplaces))## [1] 2919###Variables de la parcela
3 variables. Una con 1 NA, y 2 variables completas.
LotFrontage: Pies lineales de calle conectados a la propiedad
486 NAs. La imputación más razonable parece tomar la mediana por barrio.
ggplot(all[!is.na(all$LotFrontage),], aes(x=as.factor(Neighborhood), y=LotFrontage)) +
geom_bar(stat='summary', fun.y = "median", fill='blue') +
theme(axis.text.x = element_text(angle = 45, hjust = 1))## Warning: Ignoring unknown parameters: fun.y## No summary function supplied, defaulting to `mean_se()`
for (i in 1:nrow(all)){
if(is.na(all$LotFrontage[i])){
all$LotFrontage[i] <- as.integer(median(all$LotFrontage[all$Neighborhood==all$Neighborhood[i]], na.rm=TRUE))
}
}LotShape: Forma general de la propiedad
No NAs. Los valores parecen ordinales (Regular=best)
Reg Regular
IR1 Slightly irregular
IR2 Moderately Irregular
IR3 Irregularall$LotShape<-as.integer(revalue(all$LotShape, c('IR3'=0, 'IR2'=1, 'IR1'=2, 'Reg'=3)))
table(all$LotShape)##
## 0 1 2 3
## 16 76 968 1859sum(table(all$LotShape))## [1] 2919LotConfig: Configuración de la parcela
No NAs. The values seemed possibly ordinal to me, but the visualization does not show this. Therefore, I will convert the variable into a factor.
Inside Inside lot
Corner Corner lot
CulDSac Cul-de-sac
FR2 Frontage on 2 sides of property
FR3 Frontage on 3 sides of property
ggplot(all[!is.na(all$SalePrice),], aes(x=as.factor(LotConfig), y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue')+
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..))## Warning: Ignoring unknown parameters: fun.y## No summary function supplied, defaulting to `mean_se()`
all$LotConfig <- as.factor(all$LotConfig)
table(all$LotConfig)##
## Corner CulDSac FR2 FR3 Inside
## 511 176 85 14 2133sum(table(all$LotConfig))## [1] 2919###Garage variables
En total, hay 7 variables relacionadas con los garajes
Dos de ellas tienen un NA (GarageCars y GarageArea), una tiene 157 NAs (GarageType), 4 variables tienen 159 NAs.
En primer lugar, voy a sustituir los 159 valores que faltan de GarageYrBlt: Año en que se construyó el garaje con los valores de YearBuilt (esto es similar a YearRemodAdd, que también tiene como valor por defecto YearBuilt si no hay remodelaciones o adiciones).
all$GarageYrBlt[is.na(all$GarageYrBlt)] <- all$YearBuilt[is.na(all$GarageYrBlt)]Como las NAs significan ‘No Garage’ para las variables de carácter, ahora quiero averiguar de dónde provienen las diferencias entre la 157 NA GarageType y las otras 3 variables de carácter con 159 NAs.
#comprobar si las 157 NAs son las mismas observaciones entre las variables con 157/159 NAs
length(which(is.na(all$GarageType) & is.na(all$GarageFinish) & is.na(all$GarageCond) & is.na(all$GarageQual)))## [1] 157#Encontrar 2 NAs adicionales
kable(all[!is.na(all$GarageType) & is.na(all$GarageFinish), c('GarageCars', 'GarageArea', 'GarageType', 'GarageCond', 'GarageQual', 'GarageFinish')])| GarageCars | GarageArea | GarageType | GarageCond | GarageQual | GarageFinish | |
|---|---|---|---|---|---|---|
| 2127 | 1 | 360 | Detchd | NA | NA | NA |
| 2577 | NA | NA | Detchd | NA | NA | NA |
Los 157 NAs dentro de GarageType resultan ser todos NA en GarageCondition, GarageQuality y GarageFinish también. Las diferencias se encuentran en las casas 2127 y 2577. Como se puede ver, la casa 2127 parece tener un garaje y la casa 2577 no. Por lo tanto, debería haber 158 casas sin Garaje. Para arreglar la casa 2127, imputaré los valores más comunes (modos) para GarageCond, GarageQual y GarageFinish.
#Imputing modes.
all$GarageCond[2127] <- names(sort(-table(all$GarageCond)))[1]
all$GarageQual[2127] <- names(sort(-table(all$GarageQual)))[1]
all$GarageFinish[2127] <- names(sort(-table(all$GarageFinish)))[1]
#display "fixed" house
kable(all[2127, c('GarageYrBlt', 'GarageCars', 'GarageArea', 'GarageType', 'GarageCond', 'GarageQual', 'GarageFinish')])| GarageYrBlt | GarageCars | GarageArea | GarageType | GarageCond | GarageQual | GarageFinish | |
|---|---|---|---|---|---|---|---|
| 2127 | 1910 | 1 | 360 | Detchd | TA | TA | Unf |
GarageCars and GarageArea: Tamaño del garaje en capacidad de coches y Tamaño del garaje en cuadrados
Ambos tienen 1 NA. Como puede ver arriba, se trata de la casa 2577 para ambas variables. El problema probablemente se produjo porque el tipo de garaje de esta casa es “independiente”, mientras que todas las demás variables de garaje parecen indicar que esta casa no tiene garaje.
#fixing 3 values for house 2577
all$GarageCars[2577] <- 0
all$GarageArea[2577] <- 0
all$GarageType[2577] <- NA
#comprobar si los NAs de las variables de carácter son ahora todos 158
length(which(is.na(all$GarageType) & is.na(all$GarageFinish) & is.na(all$GarageCond) & is.na(all$GarageQual)))## [1] 158Ahora, las 4 variables de carácter relacionadas con el garaje tienen todas el mismo conjunto de 158 NAs, que corresponden a “Sin garaje”. Los arreglaré todos en el resto de esta sección
GarageType: Localización del garaje
Los valores no parecen ordinales, así que los convertiré en un factor.
2Types More than one type of garage
Attchd Attached to home
Basment Basement Garage
BuiltIn Built-In (Garage part of house - typically has room above garage)
CarPort Car Port
Detchd Detached from home
NA No Garageall$GarageType[is.na(all$GarageType)] <- 'No Garage'
all$GarageType <- as.factor(all$GarageType)
table(all$GarageType)##
## 2Types Attchd Basment BuiltIn CarPort Detchd No Garage
## 23 1723 36 186 15 778 158GarageFinish: Acabado interior del garaje
Los valores son ordinales.
Fin Finished
RFn Rough Finished
Unf Unfinished
NA No Garage all$GarageFinish[is.na(all$GarageFinish)] <- 'None'
Finish <- c('None'=0, 'Unf'=1, 'RFn'=2, 'Fin'=3)
all$GarageFinish<-as.integer(revalue(all$GarageFinish, Finish))
table(all$GarageFinish)##
## 0 1 2 3
## 158 1231 811 719GarageQual: Calidad del garaje
Otra variable que puede hacerse ordinal con el vector Cualidades.
Ex Excellent
Gd Good
TA Typical/Average
Fa Fair
Po Poor
NA No Garage
all$GarageQual[is.na(all$GarageQual)] <- 'None'
all$GarageQual<-as.integer(revalue(all$GarageQual, Qualities))
table(all$GarageQual)##
## 0 1 2 3 4 5
## 158 5 124 2605 24 3GarageCond: Condición del garaje
Otra variable que puede hacerse ordinal con el vector Cualidades.
Ex Excellent
Gd Good
TA Typical/Average
Fa Fair
Po Poor
NA No Garageall$GarageCond[is.na(all$GarageCond)] <- 'None'
all$GarageCond<-as.integer(revalue(all$GarageCond, Qualities))
table(all$GarageCond)##
## 0 1 2 3 4 5
## 158 14 74 2655 15 3###Variables del sotano
En total, hay 11 variables relacionadas con los sótanos de una casa
Cinco de ellos tienen 79-82 NAs, seis tienen uno o dos NAs.
#comprobar si las 79 NAs son las mismas observaciones entre las variables con 80+ NAs
length(which(is.na(all$BsmtQual) & is.na(all$BsmtCond) & is.na(all$BsmtExposure) & is.na(all$BsmtFinType1) & is.na(all$BsmtFinType2)))## [1] 79#Find the additional NAs; BsmtFinType1 is the one with 79 NAs
all[!is.na(all$BsmtFinType1) & (is.na(all$BsmtCond)|is.na(all$BsmtQual)|is.na(all$BsmtExposure)|is.na(all$BsmtFinType2)), c('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2')]## BsmtQual BsmtCond BsmtExposure BsmtFinType1 BsmtFinType2
## 333 Gd TA No GLQ <NA>
## 949 Gd TA <NA> Unf Unf
## 1488 Gd TA <NA> Unf Unf
## 2041 Gd <NA> Mn GLQ Rec
## 2186 TA <NA> No BLQ Unf
## 2218 <NA> Fa No Unf Unf
## 2219 <NA> TA No Unf Unf
## 2349 Gd TA <NA> Unf Unf
## 2525 TA <NA> Av ALQ UnfAsí que, en total, parece que hay 79 casas sin sótano, porque las variables de sótano de las otras casas con valores perdidos están todas completas al 80% (falta 1 de cada 5 valores). Voy a imputar los modos para arreglar esas 9 casas.
#Imputing modes.
all$BsmtFinType2[333] <- names(sort(-table(all$BsmtFinType2)))[1]
all$BsmtExposure[c(949, 1488, 2349)] <- names(sort(-table(all$BsmtExposure)))[1]
all$BsmtCond[c(2041, 2186, 2525)] <- names(sort(-table(all$BsmtCond)))[1]
all$BsmtQual[c(2218, 2219)] <- names(sort(-table(all$BsmtQual)))[1]Ahora que las 5 variables consideradas coinciden en 79 casas con “sin sótano”, voy a factorizarlas/codificarlas en caliente a continuación.
BsmtQual: Evalúa la altura del sótano
Una variable que puede convertirse en ordinal con el vector Cualidades.
Ex Excellent (100+ inches)
Gd Good (90-99 inches)
TA Typical (80-89 inches)
Fa Fair (70-79 inches)
Po Poor (<70 inches
NA No Basementall$BsmtQual[is.na(all$BsmtQual)] <- 'None'
all$BsmtQual<-as.integer(revalue(all$BsmtQual, Qualities))
table(all$BsmtQual)##
## 0 2 3 4 5
## 79 88 1285 1209 258BsmtCond: Evalúa el estado general del sótano
Una variable que puede convertirse en ordinal con el vector Cualidades.
Ex Excellent
Gd Good
TA Typical - slight dampness allowed
Fa Fair - dampness or some cracking or settling
Po Poor - Severe cracking, settling, or wetness
NA No Basementall$BsmtCond[is.na(all$BsmtCond)] <- 'None'
all$BsmtCond<-as.integer(revalue(all$BsmtCond, Qualities))
table(all$BsmtCond)##
## 0 1 2 3 4
## 79 5 104 2609 122BsmtExposure: Se refiere a los muros a nivel del jardín o de la calle
Una variable que puede convertirse en ordinal.
Gd Good Exposure
Av Average Exposure (split levels or foyers typically score average or above)
Mn Mimimum Exposure
No No Exposure
NA No Basementall$BsmtExposure[is.na(all$BsmtExposure)] <- 'None'
Exposure <- c('None'=0, 'No'=1, 'Mn'=2, 'Av'=3, 'Gd'=4)
all$BsmtExposure<-as.integer(revalue(all$BsmtExposure, Exposure))
table(all$BsmtExposure)##
## 0 1 2 3 4
## 79 1907 239 418 276BsmtFinType1: Calificación de la superficie del sótano terminado
Una variable que puede convertirse en ordinal.
GLQ Good Living Quarters
ALQ Average Living Quarters
BLQ Below Average Living Quarters
Rec Average Rec Room
LwQ Low Quality
Unf Unfinshed
NA No Basement
all$BsmtFinType1[is.na(all$BsmtFinType1)] <- 'None'
FinType <- c('None'=0, 'Unf'=1, 'LwQ'=2, 'Rec'=3, 'BLQ'=4, 'ALQ'=5, 'GLQ'=6)
all$BsmtFinType1<-as.integer(revalue(all$BsmtFinType1, FinType))
table(all$BsmtFinType1)##
## 0 1 2 3 4 5 6
## 79 851 154 288 269 429 849BsmtFinType2: Calificación de la superficie acabada del sótano (si hay varios tipos)
Una variable que puede convertirse en ordinal con el vector FinType.
GLQ Good Living Quarters
ALQ Average Living Quarters
BLQ Below Average Living Quarters
Rec Average Rec Room
LwQ Low Quality
Unf Unfinshed
NA No Basementall$BsmtFinType2[is.na(all$BsmtFinType2)] <- 'None'
FinType <- c('None'=0, 'Unf'=1, 'LwQ'=2, 'Rec'=3, 'BLQ'=4, 'ALQ'=5, 'GLQ'=6)
all$BsmtFinType2<-as.integer(revalue(all$BsmtFinType2, FinType))
table(all$BsmtFinType2)##
## 0 1 2 3 4 5 6
## 79 2494 87 105 68 52 34El resto de variabes del sótano con sólo unos pocos NAs
Ahora todavía hay que lidiar con esas 6 variables que tienen 1 o 2 NAs.
#mostrar los restantes NA. Utilizando BsmtQual como referencia para las 79 viviendas sin sótano acordadas anteriormente
all[(is.na(all$BsmtFullBath)|is.na(all$BsmtHalfBath)|is.na(all$BsmtFinSF1)|is.na(all$BsmtFinSF2)|is.na(all$BsmtUnfSF)|is.na(all$TotalBsmtSF)), c('BsmtQual', 'BsmtFullBath', 'BsmtHalfBath', 'BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF', 'TotalBsmtSF')]## BsmtQual BsmtFullBath BsmtHalfBath BsmtFinSF1 BsmtFinSF2 BsmtUnfSF
## 2121 0 NA NA NA NA NA
## 2189 0 NA NA 0 0 0
## TotalBsmtSF
## 2121 NA
## 2189 0Debería ser obvio que los restantes NAs se refieren todos a “no presente”. A continuación, estoy arreglando esas variables restantes.
BsmtFullBath: Baños completos en el sótano
Una variable entera.
all$BsmtFullBath[is.na(all$BsmtFullBath)] <-0
table(all$BsmtFullBath)##
## 0 1 2 3
## 1707 1172 38 2BsmtHalfBath: Medios baños en el sótano
Una variable entera.
all$BsmtHalfBath[is.na(all$BsmtHalfBath)] <-0
table(all$BsmtHalfBath)##
## 0 1 2
## 2744 171 4BsmtFinSF1: Tipo 1 pies cuadrados terminados
Una variable entera.
all$BsmtFinSF1[is.na(all$BsmtFinSF1)] <-0BsmtFinSF2: Tipo 2 pies cuadrados terminados
Una variable entera.
all$BsmtFinSF2[is.na(all$BsmtFinSF2)] <-0BsmtUnfSF: Pies cuadrados de superficie de sótano sin terminar
Una variable entera.
all$BsmtUnfSF[is.na(all$BsmtUnfSF)] <-0TotalBsmtSF: Total de pies cuadrados de superficie de sótano
Una variable entera.
all$TotalBsmtSF[is.na(all$TotalBsmtSF)] <-0###Variables de mampostería
Tipo de revestimiento de mampostería y área de revestimiento de mampostería
El tipo de revestimiento de mampostería tiene 24 NAs. El área de revestimiento de mampostería tiene 23 NAs. Si una casa tiene un área de revestimiento, también debería tener un tipo de revestimiento de mampostería. Arreglemos esto primero.
#compruebe si las 23 casas con superficie de chapa de madera NA son también NA en el tipo de chapa de madera
length(which(is.na(all$MasVnrType) & is.na(all$MasVnrArea)))## [1] 23#encontrar el que debe tener un MasVnrType
all[is.na(all$MasVnrType) & !is.na(all$MasVnrArea), c('MasVnrType', 'MasVnrArea')]## MasVnrType MasVnrArea
## 2611 <NA> 198#fijar este tipo de chapa imputando el modo
all$MasVnrType[2611] <- names(sort(-table(all$MasVnrType)))[2] #taking the 2nd value as the 1st is 'none'
all[2611, c('MasVnrType', 'MasVnrArea')]## MasVnrType MasVnrArea
## 2611 BrkFace 198Esto me deja con 23 casas que realmente no tienen mampostería.
Tipo de revestimiento de mampostería
Comprobará la ordinalidad a continuación.
BrkCmn Brick Common
BrkFace Brick Face
CBlock Cinder Block
None None
Stone Stoneall$MasVnrType[is.na(all$MasVnrType)] <- 'None'
all[!is.na(all$SalePrice),] %>% group_by(MasVnrType) %>% summarise(median = median(SalePrice), counts=n()) %>% arrange(median)## # A tibble: 4 x 3
## MasVnrType median counts
## <chr> <dbl> <int>
## 1 BrkCmn 139000 15
## 2 None 143125 872
## 3 BrkFace 181000 445
## 4 Stone 246839 128Parece que hay una diferencia significativa entre el “ladrillo común/ninguno” y los otros tipos. Supongo que las piedras simples y, por ejemplo, las casas de madera son simplemente más baratas. Haré la ordenación en consecuencia.
Masonry <- c('None'=0, 'BrkCmn'=0, 'BrkFace'=1, 'Stone'=2)
all$MasVnrType<-as.integer(revalue(all$MasVnrType, Masonry))
table(all$MasVnrType)##
## 0 1 2
## 1790 880 249MasVnrArea: Área de revestimiento de mampostería en pies cuadrados
An integer variable.
all$MasVnrArea[is.na(all$MasVnrArea)] <-0###MS Zoning
MSZoning: Identifica la clasificación zonal general de la venta
4 NAs. Los valores son categóricos.
A Agriculture
C Commercial
FV Floating Village Residential
I Industrial
RH Residential High Density
RL Residential Low Density
RP Residential Low Density Park
RM Residential Medium Density#imputing the mode
all$MSZoning[is.na(all$MSZoning)] <- names(sort(-table(all$MSZoning)))[1]
all$MSZoning <- as.factor(all$MSZoning)
table(all$MSZoning)##
## C (all) FV RH RL RM
## 25 139 26 2269 460sum(table(all$MSZoning))## [1] 2919###Kitchen variables
Kitchen quality and numer of Kitchens above grade
Kitchen quality has 1 NA. Number of Kitchens is complete.
Kitchen quality
1NA. Can be made ordinal with the qualities vector.
Ex Excellent
Gd Good
TA Typical/Average
Fa Fair
Po Poorall$KitchenQual[is.na(all$KitchenQual)] <- 'TA' #replace with most common value
all$KitchenQual<-as.integer(revalue(all$KitchenQual, Qualities))
table(all$KitchenQual)##
## 2 3 4 5
## 70 1493 1151 205sum(table(all$KitchenQual))## [1] 2919Number of Kitchens above grade
An integer variable with no NAs.
table(all$KitchenAbvGr)##
## 0 1 2 3
## 3 2785 129 2sum(table(all$KitchenAbvGr))## [1] 2919Please return to the 5.2 Tabs menu to select other (groups of) variables
###Utilities
Utilities: Type of utilities available
2 NAs. Ordinal as additional utilities is better.
AllPub All public Utilities (E,G,W,& S)
NoSewr Electricity, Gas, and Water (Septic Tank)
NoSeWa Electricity and Gas Only
ELO Electricity onlyHowever, the table below shows that only one house does not have all public utilities. This house is in the train set. Therefore, imputing ‘AllPub’ for the NAs means that all houses in the test set will have ‘AllPub’. This makes the variable useless for prediction. Consequently, I will get rid of it.
table(all$Utilities)##
## AllPub NoSeWa
## 2916 1kable(all[is.na(all$Utilities) | all$Utilities=='NoSeWa', 1:9])| MSSubClass | MSZoning | LotFrontage | LotArea | Street | Alley | LotShape | LandContour | Utilities | |
|---|---|---|---|---|---|---|---|---|---|
| 945 | 20 | RL | 82 | 14375 | Pave | None | 2 | Lvl | NoSeWa |
| 1916 | 30 | RL | 109 | 21780 | Grvl | None | 3 | Lvl | NA |
| 1946 | 20 | RL | 64 | 31220 | Pave | None | 2 | Bnk | NA |
all$Utilities <- NULLPlease return to the 5.2 Tabs menu to select other (groups of) variables
###Home functionality
Functional: Home functionality
1NA. Can be made ordinal (salvage only is worst, typical is best).
Typ Typical Functionality
Min1 Minor Deductions 1
Min2 Minor Deductions 2
Mod Moderate Deductions
Maj1 Major Deductions 1
Maj2 Major Deductions 2
Sev Severely Damaged
Sal Salvage only#impute mode for the 1 NA
all$Functional[is.na(all$Functional)] <- names(sort(-table(all$Functional)))[1]
all$Functional <- as.integer(revalue(all$Functional, c('Sal'=0, 'Sev'=1, 'Maj2'=2, 'Maj1'=3, 'Mod'=4, 'Min2'=5, 'Min1'=6, 'Typ'=7)))
table(all$Functional)##
## 1 2 3 4 5 6 7
## 2 9 19 35 70 65 2719sum(table(all$Functional))## [1] 2919###Exterior variables
Hay 4 variables exteriores
2 variables tienen 1 NA, 2 variables no tienen NAs.
Exterior1st: Revestimiento exterior de la casa
1 NA. Los valores son categóricos.
AsbShng Asbestos Shingles
AsphShn Asphalt Shingles
BrkComm Brick Common
BrkFace Brick Face
CBlock Cinder Block
CemntBd Cement Board
HdBoard Hard Board
ImStucc Imitation Stucco
MetalSd Metal Siding
Other Other
Plywood Plywood
PreCast PreCast
Stone Stone
Stucco Stucco
VinylSd Vinyl Siding
Wd Sdng Wood Siding
WdShing Wood Shingles#imputing mode
all$Exterior1st[is.na(all$Exterior1st)] <- names(sort(-table(all$Exterior1st)))[1]
all$Exterior1st <- as.factor(all$Exterior1st)
table(all$Exterior1st)##
## AsbShng AsphShn BrkComm BrkFace CBlock CemntBd HdBoard ImStucc MetalSd Plywood
## 44 2 6 87 2 126 442 1 450 221
## Stone Stucco VinylSd Wd Sdng WdShing
## 2 43 1026 411 56sum(table(all$Exterior1st))## [1] 2919Exterior2nd: Revestimiento exterior de la casa (si hay más de un material)
1 NA. Los valores son categóricos.
AsbShng Asbestos Shingles
AsphShn Asphalt Shingles
BrkComm Brick Common
BrkFace Brick Face
CBlock Cinder Block
CemntBd Cement Board
HdBoard Hard Board
ImStucc Imitation Stucco
MetalSd Metal Siding
Other Other
Plywood Plywood
PreCast PreCast
Stone Stone
Stucco Stucco
VinylSd Vinyl Siding
Wd Sdng Wood Siding
WdShing Wood Shingles#imputing mode
all$Exterior2nd[is.na(all$Exterior2nd)] <- names(sort(-table(all$Exterior2nd)))[1]
all$Exterior2nd <- as.factor(all$Exterior2nd)
table(all$Exterior2nd)##
## AsbShng AsphShn Brk Cmn BrkFace CBlock CmentBd HdBoard ImStucc MetalSd Other
## 38 4 22 47 3 126 406 15 447 1
## Plywood Stone Stucco VinylSd Wd Sdng Wd Shng
## 270 6 47 1015 391 81sum(table(all$Exterior2nd))## [1] 2919ExterQual: Evalúa la calidad del material en el exterior
No hay NAs. Se puede hacer ordinal utilizando el vector Cualidades.
Ex Excellent
Gd Good
TA Average/Typical
Fa Fair
Po Poor
all$ExterQual<-as.integer(revalue(all$ExterQual, Qualities))## The following `from` values were not present in `x`: None, Potable(all$ExterQual)##
## 2 3 4 5
## 35 1798 979 107sum(table(all$ExterQual))## [1] 2919ExterCond: Evalúa el estado actual del material en el exterior
No hay NAs. Se puede hacer ordinal utilizando el vector Cualidades.
Ex Excellent
Gd Good
TA Average/Typical
Fa Fair
Po Poorall$ExterCond<-as.integer(revalue(all$ExterCond, Qualities))## The following `from` values were not present in `x`: Nonetable(all$ExterCond)##
## 1 2 3 4 5
## 3 67 2538 299 12sum(table(all$ExterCond))## [1] 2919###Sistema eléctrico
Electrical: Sistema eléctrico
1 NA. Los valores son categóricos.
SBrkr Standard Circuit Breakers & Romex
FuseA Fuse Box over 60 AMP and all Romex wiring (Average)
FuseF 60 AMP Fuse Box and mostly Romex wiring (Fair)
FuseP 60 AMP Fuse Box and mostly knob & tube wiring (poor)
Mix Mixed#imputing mode
all$Electrical[is.na(all$Electrical)] <- names(sort(-table(all$Electrical)))[1]
all$Electrical <- as.factor(all$Electrical)
table(all$Electrical)##
## FuseA FuseF FuseP Mix SBrkr
## 188 50 8 1 2672sum(table(all$Electrical))## [1] 2919###Tipo y condición de venta
SaleType: Tipo de venta
1 NA. Los valores son categóricos..
WD Warranty Deed - Conventional
CWD Warranty Deed - Cash
VWD Warranty Deed - VA Loan
New Home just constructed and sold
COD Court Officer Deed/Estate
Con Contract 15% Down payment regular terms
ConLw Contract Low Down payment and low interest
ConLI Contract Low Interest
ConLD Contract Low Down
Oth Other#imputing mode
all$SaleType[is.na(all$SaleType)] <- names(sort(-table(all$SaleType)))[1]
all$SaleType <- as.factor(all$SaleType)
table(all$SaleType)##
## COD Con ConLD ConLI ConLw CWD New Oth WD
## 87 5 26 9 8 12 239 7 2526sum(table(all$SaleType))## [1] 2919SaleCondition: Condición de venta
No hay NAs. Los valores son categóricos.
Normal Normal Sale
Abnorml Abnormal Sale - trade, foreclosure, short sale
AdjLand Adjoining Land Purchase
Alloca Allocation - two linked properties with separate deeds, typically condo with a garage unit
Family Sale between family members
Partial Home was not completed when last assessed (associated with New Homes)all$SaleCondition <- as.factor(all$SaleCondition)
table(all$SaleCondition)##
## Abnorml AdjLand Alloca Family Normal Partial
## 190 12 24 46 2402 245sum(table(all$SaleCondition))## [1] 2919##Codificación de etiquetas/factorización de las restantes variables de caracteres {.tabset}
En este punto, me he asegurado de que todas las variables con NAs son atendidas. Sin embargo, todavía tengo que ocuparme de las restantes variables de caracteres que no tienen valores perdidos. Al igual que en la sección anterior, he creado pestañas para grupos de variables.
Charcol <- names(all[,sapply(all, is.character)])
Charcol## [1] "Street" "LandContour" "LandSlope" "Neighborhood" "Condition1"
## [6] "Condition2" "BldgType" "HouseStyle" "RoofStyle" "RoofMatl"
## [11] "Foundation" "Heating" "HeatingQC" "CentralAir" "PavedDrive"cat('There are', length(Charcol), 'remaining columns with character values')## There are 15 remaining columns with character values###Cimiento
Foundation: Tipo de cimiento
BrkTil Brick & Tile
CBlock Cinder Block
PConc Poured Contrete
Slab Slab
Stone Stone
Wood Wood#No ordinality, so converting into factors
all$Foundation <- as.factor(all$Foundation)
table(all$Foundation)##
## BrkTil CBlock PConc Slab Stone Wood
## 311 1235 1308 49 11 5sum(table(all$Foundation))## [1] 2919###Calefacción y aire acondicionado
Hay 2 variables de calefacción, y una que indica Airco Sí/No.
Heating: Tipo de calefacción
Floor Floor Furnace
GasA Gas forced warm air furnace
GasW Gas hot water or steam heat
Grav Gravity furnace
OthW Hot water or steam heat other than gas
Wall Wall furnace
#No ordinality, so converting into factors
all$Heating <- as.factor(all$Heating)
table(all$Heating)##
## Floor GasA GasW Grav OthW Wall
## 1 2874 27 9 2 6sum(table(all$Heating))## [1] 2919HeatingQC: Calidad y estado de la calefacción
Ex Excellent
Gd Good
TA Average/Typical
Fa Fair
Po Poor
#making the variable ordinal using the Qualities vector
all$HeatingQC<-as.integer(revalue(all$HeatingQC, Qualities))## The following `from` values were not present in `x`: Nonetable(all$HeatingQC)##
## 1 2 3 4 5
## 3 92 857 474 1493sum(table(all$HeatingQC))## [1] 2919CentralAir: Aire acondicionado central
N No
Y Yesall$CentralAir<-as.integer(revalue(all$CentralAir, c('N'=0, 'Y'=1)))
table(all$CentralAir)##
## 0 1
## 196 2723sum(table(all$CentralAir))## [1] 2919###Tejado
Hay dos variables que tienen que ver con el techo de las casas.
RoofStyle: Tipos de tejado
Flat Flat
Gable Gable
Gambrel Gabrel (Barn)
Hip Hip
Mansard Mansard
Shed Shed#No ordinality, so converting into factors
all$RoofStyle <- as.factor(all$RoofStyle)
table(all$RoofStyle)##
## Flat Gable Gambrel Hip Mansard Shed
## 20 2310 22 551 11 5sum(table(all$RoofStyle))## [1] 2919RoofMatl: Material del techo
ClyTile Clay or Tile
CompShg Standard (Composite) Shingle
Membran Membrane
Metal Metal
Roll Roll
Tar&Grv Gravel & Tar
WdShake Wood Shakes
WdShngl Wood Shingles#No ordinality, so converting into factors
all$RoofMatl <- as.factor(all$RoofMatl)
table(all$RoofMatl)##
## ClyTile CompShg Membran Metal Roll Tar&Grv WdShake WdShngl
## 1 2876 1 1 1 23 9 7sum(table(all$RoofMatl))## [1] 2919###Terreno
2 variables que especifican la planicidad y la inclinación de la propiedad.
LandContour: La planicie de la propiedad
Lvl Near Flat/Level
Bnk Banked - Quick and significant rise from street grade to building
HLS Hillside - Significant slope from side to side
Low Depression#No ordinality, so converting into factors
all$LandContour <- as.factor(all$LandContour)
table(all$LandContour)##
## Bnk HLS Low Lvl
## 117 120 60 2622sum(table(all$LandContour))## [1] 2919LandSlope: Pendiente de la propiedad
Gtl Gentle slope
Mod Moderate Slope
Sev Severe Slope#Ordinal, so label encoding
all$LandSlope<-as.integer(revalue(all$LandSlope, c('Sev'=0, 'Mod'=1, 'Gtl'=2)))
table(all$LandSlope)##
## 0 1 2
## 16 125 2778sum(table(all$LandSlope))## [1] 2919###Vivienda
2 variables que especifican el tipo y el estilo de la vivienda.
BldgType: Tipo de vivienda
1Fam Single-family Detached
2FmCon Two-family Conversion; originally built as one-family dwelling
Duplx Duplex
TwnhsE Townhouse End Unit
TwnhsI Townhouse Inside UnitEsto me parece ordinario (unifamiliar = mejor). Vamos a comprobarlo con la visualización.
ggplot(all[!is.na(all$SalePrice),], aes(x=as.factor(BldgType), y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue')+
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..))## Warning: Ignoring unknown parameters: fun.y## No summary function supplied, defaulting to `mean_se()`
Sin embargo, la visualización no muestra la ordinalidad.
#No ordinality, so converting into factors
all$BldgType <- as.factor(all$BldgType)
table(all$BldgType)##
## 1Fam 2fmCon Duplex Twnhs TwnhsE
## 2425 62 109 96 227sum(table(all$BldgType))## [1] 2919HouseStyle: Estilo de la vivienda
1Story One story
1.5Fin One and one-half story: 2nd level finished
1.5Unf One and one-half story: 2nd level unfinished
2Story Two story
2.5Fin Two and one-half story: 2nd level finished
2.5Unf Two and one-half story: 2nd level unfinished
SFoyer Split Foyer
SLvl Split Level#No ordinality, so converting into factors
all$HouseStyle <- as.factor(all$HouseStyle)
table(all$HouseStyle)##
## 1.5Fin 1.5Unf 1Story 2.5Fin 2.5Unf 2Story SFoyer SLvl
## 314 19 1471 8 24 872 83 128sum(table(all$HouseStyle))## [1] 2919###Vecindario y condiciones
3 variables que especifican la ubicación física y la proximidad de las “condiciones”.
Neighborhood: Ubicaciones físicas dentro de los límites de la ciudad de Ames
Blmngtn Bloomington Heights
Blueste Bluestem
BrDale Briardale
BrkSide Brookside
ClearCr Clear Creek
CollgCr College Creek
Crawfor Crawford
Edwards Edwards
Gilbert Gilbert
IDOTRR Iowa DOT and Rail Road
MeadowV Meadow Village
Mitchel Mitchell
Names North Ames
NoRidge Northridge
NPkVill Northpark Villa
NridgHt Northridge Heights
NWAmes Northwest Ames
OldTown Old Town
SWISU South & West of Iowa State University
Sawyer Sawyer
SawyerW Sawyer West
Somerst Somerset
StoneBr Stone Brook
Timber Timberland
Veenker Veenker#No ordinality, so converting into factors
all$Neighborhood <- as.factor(all$Neighborhood)
table(all$Neighborhood)##
## Blmngtn Blueste BrDale BrkSide ClearCr CollgCr Crawfor Edwards Gilbert IDOTRR
## 28 10 30 108 44 267 103 194 165 93
## MeadowV Mitchel NAmes NoRidge NPkVill NridgHt NWAmes OldTown Sawyer SawyerW
## 37 114 443 71 23 166 131 239 151 125
## Somerst StoneBr SWISU Timber Veenker
## 182 51 48 72 24sum(table(all$Neighborhood))## [1] 2919Condition1: Proximidad a diversas condiciones
Artery Adjacent to arterial street
Feedr Adjacent to feeder street
Norm Normal
RRNn Within 200' of North-South Railroad
RRAn Adjacent to North-South Railroad
PosN Near positive off-site feature--park, greenbelt, etc.
PosA Adjacent to postive off-site feature
RRNe Within 200' of East-West Railroad
RRAe Adjacent to East-West Railroad#No ordinality, so converting into factors
all$Condition1 <- as.factor(all$Condition1)
table(all$Condition1)##
## Artery Feedr Norm PosA PosN RRAe RRAn RRNe RRNn
## 92 164 2511 20 39 28 50 6 9sum(table(all$Condition1))## [1] 2919Condition2: Proximidad de varias condiciones (si hay más de una)
Artery Adjacent to arterial street
Feedr Adjacent to feeder street
Norm Normal
RRNn Within 200' of North-South Railroad
RRAn Adjacent to North-South Railroad
PosN Near positive off-site feature--park, greenbelt, etc.
PosA Adjacent to postive off-site feature
RRNe Within 200' of East-West Railroad
RRAe Adjacent to East-West Railroad#No ordinality, so converting into factors
all$Condition2 <- as.factor(all$Condition2)
table(all$Condition2)##
## Artery Feedr Norm PosA PosN RRAe RRAn RRNn
## 5 13 2889 4 4 1 1 2sum(table(all$Condition2))## [1] 2919###Pavimento de la calle y de la calzada
2 variables
Street: Tipo de acceso a la propiedad
Grvl Gravel
Pave Paved#Ordinal, so label encoding
all$Street<-as.integer(revalue(all$Street, c('Grvl'=0, 'Pave'=1)))
table(all$Street)##
## 0 1
## 12 2907sum(table(all$Street))## [1] 2919PavedDrive: Calzada pavimentada
Y Paved
P Partial Pavement
N Dirt/Gravel#Ordinal, so label encoding
all$PavedDrive<-as.integer(revalue(all$PavedDrive, c('N'=0, 'P'=1, 'Y'=2)))
table(all$PavedDrive)##
## 0 1 2
## 216 62 2641sum(table(all$PavedDrive))## [1] 2919##Cambio de algunas variables numéricas en factores
En este punto, todas las variables están completas (sin NAs), y todas las variables de carácter se convierten en etiquetas numéricas o en factores. Sin embargo, hay 3 variables que se registran como numéricas pero que en realidad deberían ser categóricas.
###Año y mes de venta
Aunque la oridinalidad dentro del YearBuilt (o remodelado) tiene sentido (las casas antiguas valen menos), estamos hablando de sólo 5 años de ventas. Estos años también incluyen una crisis económica. Por ejemplo: Los precios de venta en 2009 (después del colapso) es muy probable que sean mucho más bajos que en 2007. Convertiré YrSold en un factor antes de modelar, pero como necesito la versión numérica de YrSold para crear una variable Age, no lo estoy haciendo todavía.
El mes de venta también es una variable entera. Sin embargo, diciembre no es “mejor” que enero. Por lo tanto, volveré a convertir los valores de MoSold en factores.
str(all$YrSold)## int [1:2919] 2008 2007 2008 2006 2008 2009 2007 2009 2008 2008 ...str(all$MoSold)## int [1:2919] 2 5 9 2 12 10 8 11 4 1 ...all$MoSold <- as.factor(all$MoSold)Aunque posiblemente sea un poco menos pronunciado de lo que se esperaba, los efectos de las crisis bancarias que tuvieron lugar a finales de 2007 pueden apreciarse de hecho. Después de los precios medios más altos de 2007, los precios disminuyeron gradualmente. Sin embargo, la estacionalidad parece desempeñar un papel más importante, como puede verse a continuación.
ys <- ggplot(all[!is.na(all$SalePrice),], aes(x=as.factor(YrSold), y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue')+
scale_y_continuous(breaks= seq(0, 800000, by=25000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..)) +
coord_cartesian(ylim = c(0, 200000)) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice## Warning: Ignoring unknown parameters: fun.yms <- ggplot(all[!is.na(all$SalePrice),], aes(x=MoSold, y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue')+
scale_y_continuous(breaks= seq(0, 800000, by=25000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..)) +
coord_cartesian(ylim = c(0, 200000)) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice## Warning: Ignoring unknown parameters: fun.ygrid.arrange(ys, ms, widths=c(1,2))## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
###MSSubClass
MSSubClass: Identifica el tipo de vivienda implicada en la venta.
20 1-STORY 1946 & NEWER ALL STYLES
30 1-STORY 1945 & OLDER
40 1-STORY W/FINISHED ATTIC ALL AGES
45 1-1/2 STORY - UNFINISHED ALL AGES
50 1-1/2 STORY FINISHED ALL AGES
60 2-STORY 1946 & NEWER
70 2-STORY 1945 & OLDER
75 2-1/2 STORY ALL AGES
80 SPLIT OR MULTI-LEVEL
85 SPLIT FOYER
90 DUPLEX - ALL STYLES AND AGES
120 1-STORY PUD (Planned Unit Development) - 1946 & NEWER
150 1-1/2 STORY PUD - ALL AGES
160 2-STORY PUD - 1946 & NEWER
180 PUD - MULTILEVEL - INCL SPLIT LEV/FOYER
190 2 FAMILY CONVERSION - ALL STYLES AND AGESEstas clases están codificadas como números, pero realmente son categorías.
str(all$MSSubClass)## int [1:2919] 60 20 60 70 60 50 20 60 50 190 ...all$MSSubClass <- as.factor(all$MSSubClass)
#revalue for better readability
all$MSSubClass<-revalue(all$MSSubClass, c('20'='1 story 1946+', '30'='1 story 1945-', '40'='1 story unf attic', '45'='1,5 story unf', '50'='1,5 story fin', '60'='2 story 1946+', '70'='2 story 1945-', '75'='2,5 story all ages', '80'='split/multi level', '85'='split foyer', '90'='duplex all style/age', '120'='1 story PUD 1946+', '150'='1,5 story PUD all', '160'='2 story PUD 1946+', '180'='PUD multilevel', '190'='2 family conversion'))
str(all$MSSubClass)## Factor w/ 16 levels "1 story 1946+",..: 6 1 6 7 6 5 1 6 5 16 ...#Visualización de variables importantes
Por fin he llegado al punto en que todas las variables de carácter se han convertido en factores categóricos o se han codificado en números. Además, 3 variables numéricas se han convertido en factores, y he eliminado una variable (Utilities). Como puede ver a continuación, el número de variables numéricas es ahora de 56 (incluida la variable de respuesta), y las 23 variables restantes son categóricas.
numericVars <- which(sapply(all, is.numeric)) #index vector numeric variables
factorVars <- which(sapply(all, is.factor)) #index vector factor variables
cat('There are', length(numericVars), 'numeric variables, and', length(factorVars), 'categoric variables')## There are 56 numeric variables, and 23 categoric variables##Correlaciones de nuevo
all_numVar <- all[, numericVars]
cor_numVar <- cor(all_numVar, use="pairwise.complete.obs") #correlations of all numeric variables
#sort on decreasing correlations with SalePrice
cor_sorted <- as.matrix(sort(cor_numVar[,'SalePrice'], decreasing = TRUE))
#select only high corelations
CorHigh <- names(which(apply(cor_sorted, 1, function(x) abs(x)>0.5)))
cor_numVar <- cor_numVar[CorHigh, CorHigh]
corrplot.mixed(cor_numVar, tl.col="black", tl.pos = "lt", tl.cex = 0.7,cl.cex = .7, number.cex=.7)
##Encontrar la importancia de las variables con un Random Forest rápido
Aunque las correlaciones están dando una buena visión de las variables numéricas más importantes y la multicolinealidad entre esas variables, quería obtener una visión de las variables más importantes incluyendo las variables categóricas antes de pasar a la visualización.
Intenté obtener la importancia relativa de las variables con un modelo de regresión lineal rápido con la función calc.relimp del paquete , y también probé la función boruta del paquete boruta que separa las variables en grupos que son importantes o no. Sin embargo, estos métodos han llevado mucho tiempo. Como sólo quiero obtener una indicación de la importancia de las variables, finalmente decidí mantenerlo simple y utilizar un modelo rápido y sucio de Bosque Aleatorio con sólo 100 árboles. Esto también me sirve y no me lleva mucho tiempo, ya que puedo especificar un número (relativamente) pequeño de árboles.
set.seed(2018)
quick_RF <- randomForest(x=all[1:1460,-79], y=all$SalePrice[1:1460], ntree=100,importance=TRUE)
imp_RF <- importance(quick_RF)
imp_DF <- data.frame(Variables = row.names(imp_RF), MSE = imp_RF[,1])
imp_DF <- imp_DF[order(imp_DF$MSE, decreasing = TRUE),]
ggplot(imp_DF[1:20,], aes(x=reorder(Variables, MSE), y=MSE, fill=MSE)) + geom_bar(stat = 'identity') + labs(x = 'Variables', y= '% increase MSE if variable is randomly permuted') + coord_flip() + theme(legend.position="none")
Sólo 3 de esas variables más importantes son categóricas según la RF; Neighborhood, MSSubClass y GarageType.
1.0.1 Above Ground Living Area, y otras variables relacionadas con la superficie (en pies cuadrados)
As I have already visualized the relation between the Above Ground Living Area and SalePrice in my initial explorations, I will now just display the distribution itself. As there are more ‘square feet’ surface measurements in the Top 20, I am taking the opportunity to bundle them in this section. Note: GarageArea is taken care of in the Garage variables section.
I am also adding ‘Total Rooms Above Ground’ (TotRmsAbvGrd) as this variable is highly correlated with the Above Ground Living Area(0.81).
s1 <- ggplot(data= all, aes(x=GrLivArea)) +
geom_density() + labs(x='Square feet living area')
s2 <- ggplot(data=all, aes(x=as.factor(TotRmsAbvGrd))) +
geom_histogram(stat='count') + labs(x='Rooms above Ground')
s3 <- ggplot(data= all, aes(x=X1stFlrSF)) +
geom_density() + labs(x='Square feet first floor')
s4 <- ggplot(data= all, aes(x=X2ndFlrSF)) +
geom_density() + labs(x='Square feet second floor')
s5 <- ggplot(data= all, aes(x=TotalBsmtSF)) +
geom_density() + labs(x='Square feet basement')
s6 <- ggplot(data= all[all$LotArea<100000,], aes(x=LotArea)) +
geom_density() + labs(x='Square feet lot')
s7 <- ggplot(data= all, aes(x=LotFrontage)) +
geom_density() + labs(x='Linear feet lot frontage')
s8 <- ggplot(data= all, aes(x=LowQualFinSF)) +
geom_histogram() + labs(x='Low quality square feet 1st & 2nd')
layout <- matrix(c(1,2,5,3,4,8,6,7),4,2,byrow=TRUE)
multiplot(s1, s2, s3, s4, s5, s6, s7, s8, layout=layout)
Más adelante investigaré varias de estas variables en busca de valores atípicos. Para la visualización del lote, ya he sacado los lotes de más de 100.000 pies cuadrados (4 casas).
GrLivArea parecía ser sólo el total de pies cuadrados del primer y segundo piso. Sin embargo, en una versión posterior, descubrí que también hay una variable llamada LowQualFinSF: Pies cuadrados terminados de baja calidad (todos los pisos). Como se puede ver arriba (Low quality square feet 1st and 2nd) casi todas las casas no tienen nada de esto (sólo 40 casas tienen algo). Resulta que estos pies cuadrados sí están incluidos en el GrLivArea. La correlación entre esas 3 variables y GrLivArea es exactamente 1.
cor(all$GrLivArea, (all$X1stFlrSF + all$X2ndFlrSF + all$LowQualFinSF))## [1] 1head(all[all$LowQualFinSF>0, c('GrLivArea', 'X1stFlrSF', 'X2ndFlrSF', 'LowQualFinSF')])## GrLivArea X1stFlrSF X2ndFlrSF LowQualFinSF
## 52 1176 816 0 360
## 89 1526 1013 0 513
## 126 754 520 0 234
## 171 1382 854 0 528
## 186 3608 1518 1518 572
## 188 1656 808 704 144###La variable categórica más importante;Neighborhood
El primer gráfico muestra la mediana del precio de venta por barrio. La frecuencia (número de casas) de cada barrio en el conjunto de entrenamiento se muestra en las etiquetas.
El segundo gráfico muestra las frecuencias de todos los datos.
n1 <- ggplot(all[!is.na(all$SalePrice),], aes(x=Neighborhood, y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue') +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice
n2 <- ggplot(data=all, aes(x=Neighborhood)) +
geom_histogram(stat='count')+
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3)+
theme(axis.text.x = element_text(angle = 45, hjust = 1))
grid.arrange(n1, n2)## No summary function supplied, defaulting to `mean_se()`
###Calidad global y otras variables de calidad
Ya he visualizado la relación entre la calidad general y el precio de venta en mis exploraciones iniciales, pero quiero visualizar también la distribución de frecuencias. Como hay más medidas de calidad, aprovecho para agruparlas en esta sección.
q1 <- ggplot(data=all, aes(x=as.factor(OverallQual))) +
geom_histogram(stat='count')
q2 <- ggplot(data=all, aes(x=as.factor(ExterQual))) +
geom_histogram(stat='count')
q3 <- ggplot(data=all, aes(x=as.factor(BsmtQual))) +
geom_histogram(stat='count')
q4 <- ggplot(data=all, aes(x=as.factor(KitchenQual))) +
geom_histogram(stat='count')
q5 <- ggplot(data=all, aes(x=as.factor(GarageQual))) +
geom_histogram(stat='count')
q6 <- ggplot(data=all, aes(x=as.factor(FireplaceQu))) +
geom_histogram(stat='count')
q7 <- ggplot(data=all, aes(x=as.factor(PoolQC))) +
geom_histogram(stat='count')
layout <- matrix(c(1,2,8,3,4,8,5,6,7),3,3,byrow=TRUE)
multiplot(q1, q2, q3, q4, q5, q6, q7, layout=layout)
La Calidad Global es muy importante, y también más granular que las otras variables. La Calidad externa también es improtante, pero tiene una alta correlación con la Calidad general (0,73). La Calidad de la Cocina también parece una de las que hay que mantener, ya que todas las casas tienen cocina y hay una varianza con cierta sustancia. La Calidad del Garaje no parece distinguirse mucho, ya que la mayoría de los garajes tienen Q3. La Calidad de la Chimenea está en la lista de correlaciones altas, y en la lista de variables importantes. El PoolQC es sólo muy escaso (las 13 piscinas ni siquiera se puede ver en esta escala). Miraré de crear una variable “tiene piscina” más adelante.
###La segunda variable categórica más importante; MSSubClass
La primera visualización muestra la mediana del precio de venta por MSSubClass. La frecuencia (número de casas) de cada MSSubClass en el conjunto de entrenamiento se muestra en las etiquetas.
El histrograma muestra las frecuencias de todos los datos. La mayoría de las casas son relativamente nuevas y tienen una o dos plantas.
ms1 <- ggplot(all[!is.na(all$SalePrice),], aes(x=MSSubClass, y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue') +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice
ms2 <- ggplot(data=all, aes(x=MSSubClass)) +
geom_histogram(stat='count')+
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
grid.arrange(ms1, ms2)## No summary function supplied, defaulting to `mean_se()`
###Garaje variables
Varias variables de Garaje tienen una alta correlación con Precio de Venta, y también están en la lista de las 20 primeras del bosque aleatorio rápido. Sin embargo, hay multicolinealidad entre ellas y creo que 7 variables de garaje son demasiadas de todos modos. Creo que algo así como 3 variables deberían ser suficientes (posiblemente GarageCars, GarageType, y una medida de Calidad), pero antes de hacer cualquier selección estoy visualizando todas ellas en esta sección.
#correct error
all$GarageYrBlt[2593] <- 2007 #this must have been a typo. GarageYrBlt=2207, YearBuilt=2006, YearRemodAdd=2007.g1 <- ggplot(data=all[all$GarageCars !=0,], aes(x=GarageYrBlt)) +
geom_histogram()
g2 <- ggplot(data=all, aes(x=as.factor(GarageCars))) +
geom_histogram(stat='count')
g3 <- ggplot(data= all, aes(x=GarageArea)) +
geom_density()
g4 <- ggplot(data=all, aes(x=as.factor(GarageCond))) +
geom_histogram(stat='count')
g5 <- ggplot(data=all, aes(x=GarageType)) +
geom_histogram(stat='count')
g6 <- ggplot(data=all, aes(x=as.factor(GarageQual))) +
geom_histogram(stat='count')
g7 <- ggplot(data=all, aes(x=as.factor(GarageFinish))) +
geom_histogram(stat='count')
layout <- matrix(c(1,5,5,2,3,8,6,4,7),3,3,byrow=TRUE)
multiplot(g1, g2, g3, g4, g5, g6, g7, layout=layout)
Como ya se mencionó en la sección 4.2, GarageCars y GarageArea están altamente correlacionados. Aquí, GarageQual y GarageCond también parecen altamente correlacionados, y ambos están dominados por el nivel =3.
###Variables del sótano
Al igual que las variables del garaje, las múltiples variables del sótano son importantes en la matriz de correlaciones y en la lista de los 20 mejores predictores de la RF. Sin embargo, 11 variables de sótano parecen excesivas. Antes de decidir qué voy a hacer con ellas, a continuación visualizo 8 de ellas. Las 2 variables “Cuarto de baño” se tratan en la ingeniería de características (sección 7.1), y los “Pies cuadrados del sótano” ya se tratan en la sección 6.2.1.
b1 <- ggplot(data=all, aes(x=BsmtFinSF1)) +
geom_histogram() + labs(x='Type 1 finished square feet')
b2 <- ggplot(data=all, aes(x=BsmtFinSF2)) +
geom_histogram()+ labs(x='Type 2 finished square feet')
b3 <- ggplot(data=all, aes(x=BsmtUnfSF)) +
geom_histogram()+ labs(x='Unfinished square feet')
b4 <- ggplot(data=all, aes(x=as.factor(BsmtFinType1))) +
geom_histogram(stat='count')+ labs(x='Rating of Type 1 finished area')
b5 <- ggplot(data=all, aes(x=as.factor(BsmtFinType2))) +
geom_histogram(stat='count')+ labs(x='Rating of Type 2 finished area')
b6 <- ggplot(data=all, aes(x=as.factor(BsmtQual))) +
geom_histogram(stat='count')+ labs(x='Height of the basement')
b7 <- ggplot(data=all, aes(x=as.factor(BsmtCond))) +
geom_histogram(stat='count')+ labs(x='Rating of general condition')
b8 <- ggplot(data=all, aes(x=as.factor(BsmtExposure))) +
geom_histogram(stat='count')+ labs(x='Walkout or garden level walls')
layout <- matrix(c(1,2,3,4,5,9,6,7,8),3,3,byrow=TRUE)
multiplot(b1, b2, b3, b4, b5, b6, b7, b8, layout=layout)
Así que parecía como si la superficie total del sótano en pies cuadrados (TotalBsmtSF) se desglosa en áreas terminadas (2 si más de un tipo de acabado), y el área sin terminar. Hice una comprobación entre la correlación del total de esas 3 variables, y TotalBsmtSF. La correlación es exactamente 1, así que eso es algo bueno (sin errores o pequeñas discrepancias).
La calidad del sótano es un nombre de variable confuso, ya que resulta que califica específicamente la altura del sótano.
#Ingeniería de funciones
##Número total de baños
Hay 4 variables de baño. Por separado, estas variables no son muy importantes. Sin embargo, supongo que si las sumo en un solo predictor, es probable que éste se convierta en un predictor fuerte.
“Un medio baño, también conocido como tocador o baño de invitados, sólo tiene dos de los cuatro componentes principales del cuarto de baño, normalmente un inodoro y un lavabo”. En consecuencia, también contaré los medios baños como medios.
all$TotBathrooms <- all$FullBath + (all$HalfBath*0.5) + all$BsmtFullBath + (all$BsmtHalfBath*0.5)Como se puede ver en el primer gráfico, ahora parece haber una clara correlación (es de 0,63). La distribución de frecuencias de los Baños en todos los datos se muestra en el segundo gráfico.
tb1 <- ggplot(data=all[!is.na(all$SalePrice),], aes(x=as.factor(TotBathrooms), y=SalePrice))+
geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma)
tb2 <- ggplot(data=all, aes(x=as.factor(TotBathrooms))) +
geom_histogram(stat='count')
grid.arrange(tb1, tb2)## `geom_smooth()` using formula 'y ~ x'
##Añadir las variables “Edad de la casa”, “Remodelado (Sí/No)” y “Es nuevo”.
En total, hay 3 variables que son relevantes con respecto a la Edad de una casa; YearBlt, YearRemodAdd, y YearSold. YearRemodAdd es por defecto YearBuilt si no ha habido ninguna remodelación/adición. Utilizaré YearRemodeled y YearSold para determinar la edad. Sin embargo, como siempre quedarán partes de construcciones antiguas y puede que sólo se hayan renovado partes de la casa, también introduciré una variable Remodelado Sí/No. Esto debería verse como una especie de parámetro de penalización que indica que si la Edad se basa en una fecha de remodelación, probablemente valga menos que las casas que se construyeron desde cero en ese mismo año.
all$Remod <- ifelse(all$YearBuilt==all$YearRemodAdd, 0, 1) #0=No Remodeling, 1=Remodeling
all$Age <- as.numeric(all$YrSold)-all$YearRemodAddggplot(data=all[!is.na(all$SalePrice),], aes(x=Age, y=SalePrice))+
geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma)## `geom_smooth()` using formula 'y ~ x'
As expected, the graph shows a negative correlation with Age (old house are worth less).
cor(all$SalePrice[!is.na(all$SalePrice)], all$Age[!is.na(all$SalePrice)])## [1] -0.5090787ggplot(all[!is.na(all$SalePrice),], aes(x=as.factor(Remod), y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue') +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=6) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
theme_grey(base_size = 18) +
geom_hline(yintercept=163000, linetype="dashed") #dashed line is median SalePrice## Warning: Ignoring unknown parameters: fun.y## No summary function supplied, defaulting to `mean_se()`
Por último, voy a crear la variable IsNew a continuación. En total, hay
116 casas nuevas en el conjunto de datos.
all$IsNew <- ifelse(all$YrSold==all$YearBuilt, 1, 0)
table(all$IsNew)##
## 0 1
## 2803 116Estas 116 casas nuevas se distribuyen de forma bastante equitativa entre el conjunto de entrenamiento y el de prueba, y como se puede ver las casas nuevas valen considerablemente más en promedio.
ggplot(all[!is.na(all$SalePrice),], aes(x=as.factor(IsNew), y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue') +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=6) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
theme_grey(base_size = 18) +
geom_hline(yintercept=163000, linetype="dashed") #dashed line is median SalePrice## Warning: Ignoring unknown parameters: fun.y## No summary function supplied, defaulting to `mean_se()`
all$YrSold <- as.factor(all$YrSold) #the numeric version is now not needed anymore##Binning Neighborhood
nb1 <- ggplot(all[!is.na(all$SalePrice),], aes(x=reorder(Neighborhood, SalePrice, FUN=median), y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue') + labs(x='Neighborhood', y='Median SalePrice') +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice## Warning: Ignoring unknown parameters: fun.ynb2 <- ggplot(all[!is.na(all$SalePrice),], aes(x=reorder(Neighborhood, SalePrice, FUN=mean), y=SalePrice)) +
geom_bar(stat='summary', fun.y = "mean", fill='blue') + labs(x='Neighborhood', y="Mean SalePrice") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice## Warning: Ignoring unknown parameters: fun.ygrid.arrange(nb1, nb2)## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
Tanto la mediana como la media de los precios de venta coinciden en 3 barrios con precios de venta sustancialmente más altos. La separación de los 3 barrios relativamente pobres es menos clara, pero al menos ambos gráficos coinciden en los mismos 3 barrios pobres. Como no quiero “sobreponderar”, sólo estoy creando categorías para esos “extremos”.
all$NeighRich[all$Neighborhood %in% c('StoneBr', 'NridgHt', 'NoRidge')] <- 2
all$NeighRich[!all$Neighborhood %in% c('MeadowV', 'IDOTRR', 'BrDale', 'StoneBr', 'NridgHt', 'NoRidge')] <- 1
all$NeighRich[all$Neighborhood %in% c('MeadowV', 'IDOTRR', 'BrDale')] <- 0table(all$NeighRich)##
## 0 1 2
## 160 2471 288##Pies cuadrados totales
Como el espacio habitable total suele ser muy importante cuando la gente compra casas, estoy añadiendo un predictor que suma el espacio habitable por encima y por debajo del suelo.
all$TotalSqFeet <- all$GrLivArea + all$TotalBsmtSFggplot(data=all[!is.na(all$SalePrice),], aes(x=TotalSqFeet, y=SalePrice))+
geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +
geom_text_repel(aes(label = ifelse(all$GrLivArea[!is.na(all$SalePrice)]>4500, rownames(all), '')))## `geom_smooth()` using formula 'y ~ x'
Como era de esperar, la correlación con SalePrice es muy fuerte (0,78).
cor(all$SalePrice, all$TotalSqFeet, use= "pairwise.complete.obs")## [1] 0.7789588Los dos posibles valores atípicos parecen “salir” aún más que antes. Al eliminar estos dos valores atípicos, la correlación aumenta un 5%.
cor(all$SalePrice[-c(524, 1299)], all$TotalSqFeet[-c(524, 1299)], use= "pairwise.complete.obs")## [1] 0.829042##Consolidación de las variables del porche
A continuación, he enumerado las variables que parecen estar relacionadas con los porches.
WoodDeckSF: Área de la cubierta de madera en pies cuadrados
OpenPorchSF: Área de porche abierto en pies cuadrados
EnclosedPorch: Área de porche cerrado en pies cuadrados
3SsnPorch: Área de porche de tres estaciones en pies cuadrados
ScreenPorch: Superficie del porche en pies cuadrados
Por lo que sé, los porches son áreas protegidas fuera de la casa, y una cubierta de madera no está protegida. Por lo tanto, estoy dejando WoodDeckSF solo, y sólo se consolidan las 4 variables de porche.
all$TotalPorchSF <- all$OpenPorchSF + all$EnclosedPorch + all$X3SsnPorch + all$ScreenPorchAlthough adding up these Porch areas makes sense (there should not be any overlap between areas), the correlation with SalePrice is not very strong.
cor(all$SalePrice, all$TotalPorchSF, use= "pairwise.complete.obs")## [1] 0.1957389ggplot(data=all[!is.na(all$SalePrice),], aes(x=TotalPorchSF, y=SalePrice))+
geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma)## `geom_smooth()` using formula 'y ~ x'
#Preparación de los datos para la modelización
##Eliminación de variables altamente correlacionadas
En primer lugar, descarto una variable si dos variables están muy correlacionadas. Para encontrar estos pares correlacionados, he vuelto a utilizar la matriz de correlaciones (véase la sección 6.1). Por ejemplo: GarajeCoches y GarajeÁrea tienen una correlación de 0,89. De estas dos, estoy descartando la variable con la menor correlación con Precio de venta (que es Área de garaje con una correlación Precio de venta de 0,62. GarageCars tiene una correlación con SalePrice de 0,64).
dropVars <- c('YearRemodAdd', 'GarageYrBlt', 'GarageArea', 'GarageCond', 'TotalBsmtSF', 'TotalRmsAbvGrd', 'BsmtFinSF1')
all <- all[,!(names(all) %in% dropVars)]##Eliminación de valores atípicos
Por el momento, me limito a eliminar manualmente las dos casas realmente grandes con un precio de venta bajo. Sin embargo, tengo la intención de investigar esto más a fondo en una etapa posterior (posiblemente utilizando el paquete “outliers”).
all <- all[-c(524, 1299),]##Preprocesamiento de las variables predictoras
Antes de modelar necesito centrar y escalar los predictores “numéricos verdaderos” (es decir, las variables que han sido codificadas con etiquetas), y crear variables ficticias para los predictores categóricos. A continuación, divido el marco de datos en uno con todas las variables numéricas (verdaderas) y otro marco de datos con los factores (ordinales).
numericVarNames <- numericVarNames[!(numericVarNames %in% c('MSSubClass', 'MoSold', 'YrSold', 'SalePrice', 'OverallQual', 'OverallCond'))] #numericVarNames was created before having done anything
numericVarNames <- append(numericVarNames, c('Age', 'TotalPorchSF', 'TotBathrooms', 'TotalSqFeet'))
DFnumeric <- all[, names(all) %in% numericVarNames]
DFfactors <- all[, !(names(all) %in% numericVarNames)]
DFfactors <- DFfactors[, names(DFfactors) != 'SalePrice']
cat('There are', length(DFnumeric), 'numeric variables, and', length(DFfactors), 'factor variables')## There are 30 numeric variables, and 49 factor variables###Asimetría y normalización de los predictores numéricos
Asimetría La asimetría es una medida de la simetría de una distribución. Un conjunto de datos simétricos tendrá una asimetría igual a 0. Así, una distribución normal tendrá una asimetría de 0. La asimetría mide esencialmente el tamaño relativo de las dos colas. Como regla general, la asimetría debe estar entre -1 y 1. En este rango, los datos se consideran bastante simétricos. Para corregir la asimetría, estoy tomando el logaritmo para todos los predictores numéricos con una asimetría absoluta superior a 0,8 (en realidad: log+1, para evitar problemas de división por cero).
for(i in 1:ncol(DFnumeric)){
if (abs(skew(DFnumeric[,i]))>0.8){
DFnumeric[,i] <- log(DFnumeric[,i] +1)
}
}Normalizing the data
PreNum <- preProcess(DFnumeric, method=c("center", "scale"))
print(PreNum)## Created from 2917 samples and 30 variables
##
## Pre-processing:
## - centered (30)
## - ignored (0)
## - scaled (30)DFnorm <- predict(PreNum, DFnumeric)
dim(DFnorm)## [1] 2917 30###One hot encoding las variables categóricas
El último paso necesario para garantizar que todos los predictores se
conviertan en columnas numéricas (lo cual es necesario para la mayoría
de los algoritmos de aprendizaje automático) es “one-hot encode” las
variables categóricas. Esto significa, básicamente, que todos los
valores de los factores (no ordinales) reciben columnas separadas con 1
y 0 (1 significa básicamente Sí/Presente). Para hacer esta one-hot
encoding, estoy utilizando la función model.matrix().
DFdummies <- as.data.frame(model.matrix(~.-1, DFfactors))
dim(DFdummies)## [1] 2917 201###Eliminación de niveles con pocas o ninguna observación en el tren o la prueba
En versiones anteriores, trabajé con la función
Varianza cercana a cero de Caret. Aunque esto funciona,
también es una solución rápida y se pierde demasiada información. Por
ejemplo, al utilizar los valores predeterminados, todos los barrios con
menos de 146 casas se omiten como variables (codificadas con un punto)
(relación de frecuencia superior a 95/5). Por lo tanto, he tomado un
enfoque manual más cuidadoso en esta versión.
#check if some values are absent in the test set
ZerocolTest <- which(colSums(DFdummies[(nrow(all[!is.na(all$SalePrice),])+1):nrow(all),])==0)
colnames(DFdummies[ZerocolTest])## [1] "Condition2RRAe" "Condition2RRAn" "Condition2RRNn"
## [4] "HouseStyle2.5Fin" "RoofMatlMembran" "RoofMatlMetal"
## [7] "RoofMatlRoll" "Exterior1stImStucc" "Exterior1stStone"
## [10] "Exterior2ndOther" "HeatingOthW" "ElectricalMix"
## [13] "MiscFeatureTenC"DFdummies <- DFdummies[,-ZerocolTest] #removing predictors#check if some values are absent in the train set
ZerocolTrain <- which(colSums(DFdummies[1:nrow(all[!is.na(all$SalePrice),]),])==0)
colnames(DFdummies[ZerocolTrain])## [1] "MSSubClass1,5 story PUD all"DFdummies <- DFdummies[,-ZerocolTrain] #removing predictorTambién se eliminan las variables con menos de 10 “unos” en el conjunto de trenes.
fewOnes <- which(colSums(DFdummies[1:nrow(all[!is.na(all$SalePrice),]),])<10)
colnames(DFdummies[fewOnes])## [1] "MSSubClass1 story unf attic" "LotConfigFR3"
## [3] "NeighborhoodBlueste" "NeighborhoodNPkVill"
## [5] "Condition1PosA" "Condition1RRNe"
## [7] "Condition1RRNn" "Condition2Feedr"
## [9] "Condition2PosA" "Condition2PosN"
## [11] "RoofStyleMansard" "RoofStyleShed"
## [13] "RoofMatlWdShake" "RoofMatlWdShngl"
## [15] "Exterior1stAsphShn" "Exterior1stBrkComm"
## [17] "Exterior1stCBlock" "Exterior2ndAsphShn"
## [19] "Exterior2ndBrk Cmn" "Exterior2ndCBlock"
## [21] "Exterior2ndStone" "FoundationStone"
## [23] "FoundationWood" "HeatingGrav"
## [25] "HeatingWall" "ElectricalFuseP"
## [27] "GarageTypeCarPort" "MiscFeatureOthr"
## [29] "SaleTypeCon" "SaleTypeConLD"
## [31] "SaleTypeConLI" "SaleTypeConLw"
## [33] "SaleTypeCWD" "SaleTypeOth"
## [35] "SaleConditionAdjLand"DFdummies <- DFdummies[,-fewOnes] #removing predictors
dim(DFdummies)## [1] 2917 152En total, he eliminado 49 predictores codificados de un solo golpe con poca o ninguna varianza. Aunque pueda parecer un número significativo, en realidad es mucho menor que el número de predictores que se eliminaron utilizando la función de varianza cercana a cero de caret (utilizando sus umbrales por defecto).
combined <- cbind(DFnorm, DFdummies) #combining all (now numeric) predictors into one dataframe ##Tratamiento de la asimetría de la variable de respuesta
skew(all$SalePrice)## [1] 1.877427qqnorm(all$SalePrice)
qqline(all$SalePrice)
El sesgo de 1,87 indica un sesgo a la derecha demasiado alto, y el gráfico Q-Q muestra que los precios de venta tampoco se distribuyen normalmente. Para solucionar esto, estoy tomando el logaritmo de SalePrice.
all$SalePrice <- log(all$SalePrice) #default is the natural logarithm, "+1" is not necessary as there are no 0's
skew(all$SalePrice)## [1] 0.1213182Como puedes ver, el sesgo es ahora bastante bajo y el gráfico Q-Q también se ve mucho mejor.
qqnorm(all$SalePrice)
qqline(all$SalePrice)
##Composición de los conjuntos de entrenamiento y prueba
train1 <- combined[!is.na(all$SalePrice),]
test1 <- combined[is.na(all$SalePrice),]#Modelado
##Modelo de regresión Lasso
También he probado los modelos Ridge y Elastic Net, pero como el lasso da los mejores resultados de esos 3 modelos, sólo mantengo el modelo lasso en el documento.
La penalización de la red elástica está controlada por alfa, y cubre la brecha entre el lazo (alfa=1) y la cresta (alfa=0). El parámetro de ajuste lambda controla la fuerza general de la penalización. Se sabe que la penalización de cresta encoge los coeficientes de los predictores correlacionados entre sí, mientras que el lazo tiende a elegir uno de ellos y descartar los demás.
A continuación, utilizo la validación cruzada de caret para encontrar el mejor valor de lambda, que es el único hiperparámetro que debe ajustarse para el modelo lasso.
set.seed(27042018)
my_control <-trainControl(method="cv", number=5)
lassoGrid <- expand.grid(alpha = 1, lambda = seq(0.001,0.1,by = 0.0005))
lasso_mod <- train(x=train1, y=all$SalePrice[!is.na(all$SalePrice)], method='glmnet', trControl= my_control, tuneGrid=lassoGrid)
lasso_mod$bestTune## alpha lambda
## 1 1 0.001min(lasso_mod$results$RMSE)## [1] 0.1127325La documentación de la función caret `varImp’ dice: para glmboost y glmnet se utiliza el valor absoluto de los coeficientes correspondientes al modelo ajustado.
Aunque esto significa que no se almacena una clasificación real de las variables más importantes, me da la oportunidad de averiguar cuántas de las variables no se utilizan en el modelo (y por tanto tienen coeficiente 0).
lassoVarImp <- varImp(lasso_mod,scale=F)
lassoImportance <- lassoVarImp$importance
varsSelected <- length(which(lassoImportance$Overall!=0))
varsNotSelected <- length(which(lassoImportance$Overall==0))
cat('Lasso uses', varsSelected, 'variables in its model, and did not select', varsNotSelected, 'variables.')## Lasso uses 132 variables in its model, and did not select 50 variables.Así que el lazo hizo lo que se supone que tiene que hacer: parece haber tratado bien la multicolinealidad al no utilizar cerca del 45% de las variables disponibles en el modelo.
LassoPred <- predict(lasso_mod, test1)
predictions_lasso <- exp(LassoPred) #need to reverse the log to the real values
head(predictions_lasso)## 1461 1462 1463 1464 1465 1466
## 115562.1 162498.8 179170.1 197566.9 205357.5 168499.1##Modelo XGBoost
Inicialmente, sólo trabajé con el paquete XGBoost directamente. La razón principal para ello era que el paquete utiliza su propia estructura de datos eficiente (xgb.DMatrix). El paquete también proporciona una función de validación cruzada. Sin embargo, esta función de CV sólo determina el número óptimo de rondas, y no admite una búsqueda completa de hiperparámetros.
Aunque caret no parece utilizar la estructura de datos (rápida) del paquete xgb, finalmente decidí hacer el ajuste de hiperparámetros con él de todos modos, ya que al menos soporta una búsqueda de cuadrícula completa. Según tengo entendido, los principales parámetros que hay que ajustar para evitar el sobreajuste son max_depth, y min_child_weight (ver documentación de XGBoost). A continuación estoy configurando una rejilla que ajusta ambos parámetros, y también la eta (tasa de aprendizaje).
xgb_grid = expand.grid(
nrounds = 1000,
eta = c(0.1, 0.05, 0.01),
max_depth = c(2, 3, 4, 5, 6),
gamma = 0,
colsample_bytree=1,
min_child_weight=c(1, 2, 3, 4 ,5),
subsample=1
)El siguiente paso es dejar que caret encuentre los mejores valores de los hiperparámetros (utilizando una validación cruzada de 5 veces).
#xgb_caret <- train(x=train1, y=all$SalePrice[!is.na(all$SalePrice)], method='xgbTree', trControl= my_control, tuneGrid=xgb_grid)
#xgb_caret$bestTuneComo era de esperar, esto tomó bastante tiempo (localmente). Como quiero limitar el tiempo de ejecución en Kaggle, deshabilité el código, y sólo continúo con los resultados. Según caret, los parámetros de ‘bestTune’ son:
- Max_depth=3
- eta=0.05
- Min_child_weight=4
En el resto de esta sección, seguiré trabajando con el paquete xgboost directamente. A continuación, empiezo con la preparación de los datos en el formato recomendado.
label_train <- all$SalePrice[!is.na(all$SalePrice)]
# put our testing & training data into two seperates Dmatrixs objects
dtrain <- xgb.DMatrix(data = as.matrix(train1), label= label_train)
dtest <- xgb.DMatrix(data = as.matrix(test1))Además, tomo los valores mejor sintonizados de la validación cruzada del caret.
default_param<-list(
objective = "reg:linear",
booster = "gbtree",
eta=0.05, #default = 0.3
gamma=0,
max_depth=3, #default=6
min_child_weight=4, #default=1
subsample=1,
colsample_bytree=1
)El siguiente paso es hacer una validación cruzada para determinar el mejor número de rondas (para el conjunto de parámetros dado).
xgbcv <- xgb.cv( params = default_param, data = dtrain, nrounds = 500, nfold = 5, showsd = T, stratified = T, print_every_n = 40, early_stopping_rounds = 10, maximize = F)## [19:12:55] WARNING: amalgamation/../src/objective/regression_obj.cu:203: reg:linear is now deprecated in favor of reg:squarederror.
## [19:12:55] WARNING: amalgamation/../src/objective/regression_obj.cu:203: reg:linear is now deprecated in favor of reg:squarederror.
## [19:12:55] WARNING: amalgamation/../src/objective/regression_obj.cu:203: reg:linear is now deprecated in favor of reg:squarederror.
## [19:12:55] WARNING: amalgamation/../src/objective/regression_obj.cu:203: reg:linear is now deprecated in favor of reg:squarederror.
## [19:12:55] WARNING: amalgamation/../src/objective/regression_obj.cu:203: reg:linear is now deprecated in favor of reg:squarederror.
## [1] train-rmse:10.955589+0.004521 test-rmse:10.955560+0.019198
## Multiple eval metrics are present. Will use test_rmse for early stopping.
## Will train until test_rmse hasn't improved in 10 rounds.
##
## [41] train-rmse:1.428255+0.000799 test-rmse:1.429178+0.017698
## [81] train-rmse:0.219896+0.000924 test-rmse:0.232053+0.013013
## [121] train-rmse:0.102608+0.002015 test-rmse:0.130365+0.011149
## [161] train-rmse:0.090555+0.002014 test-rmse:0.123322+0.011224
## [201] train-rmse:0.084558+0.002009 test-rmse:0.120920+0.011448
## [241] train-rmse:0.079900+0.001708 test-rmse:0.119551+0.011576
## [281] train-rmse:0.076070+0.001421 test-rmse:0.118870+0.011577
## [321] train-rmse:0.073007+0.001213 test-rmse:0.118383+0.011633
## [361] train-rmse:0.070253+0.001015 test-rmse:0.117983+0.011803
## Stopping. Best iteration:
## [358] train-rmse:0.070455+0.001041 test-rmse:0.117967+0.011768Aunque ha costado un poco de trabajo, el ajuste de los hiperparámetros ha merecido la pena, ya que el RMSE validado de forma cruzada ha mejorado considerablemente (de 0,1225 sin el ajuste de caret, a 0,1162 en esta versión).
#train the model using the best iteration found by cross validation
xgb_mod <- xgb.train(data = dtrain, params=default_param, nrounds = 454)## [19:13:01] WARNING: amalgamation/../src/objective/regression_obj.cu:203: reg:linear is now deprecated in favor of reg:squarederror.XGBpred <- predict(xgb_mod, dtest)
predictions_XGB <- exp(XGBpred) #need to reverse the log to the real values
head(predictions_XGB)## [1] 116387.0 162307.1 186493.8 187440.4 187258.2 166241.1#view variable importance plot
library(Ckmeans.1d.dp) #required for ggplot clustering## Warning: package 'Ckmeans.1d.dp' was built under R version 4.1.3mat <- xgb.importance (feature_names = colnames(train1),model = xgb_mod)
xgb.ggplot.importance(importance_matrix = mat[1:20], rel_to_first = TRUE)
##Promedio de predicciones
Dado que los algoritmos de lazo y XGBoost son muy diferentes, promediar las predicciones probablemente mejore las puntuaciones. Como el modelo de lazo es mejor en cuanto a la puntuación de RMSE validada (0,1121 frente a 0,1162), estoy ponderando el modelo de lazo al doble.
sub_avg <- data.frame(Id = test_labels, SalePrice = (predictions_XGB+2*predictions_lasso)/3)
head(sub_avg)## Id SalePrice
## 1461 1461 115837.0
## 1462 1462 162434.9
## 1463 1463 181611.3
## 1464 1464 194191.4
## 1465 1465 199324.4
## 1466 1466 167746.4write.csv(sub_avg, file = 'average.csv', row.names = F)