Modelo Preditivo: Valor de Casa

1 Introdução

O Kaggle é uma plataforma gratuita que tem como objetivo inicial promover desafios entre profissionais da área de análise, porém atualmente se tornou uma grande comunidade de Data Sciente.

Para demonstrar o trabalho de modelagem de dados, foi coletado um dataset fornecido no campeonato “House Sales in King Country, USA” da plataforma Kaggle

# Carregando os pacotes necessários

pacman::p_load("ggplot2",'tidyverse',"plyr","cowplot","plotly","factoextra","corrplot","Boruta","lubridate","data.table", "ggthemes", "ggrepel", "GGally", "gridExtra", "grid", "caret","lubridate", "hydroGOF", "miscTools", "MASS","reshape2", "randomForest", "pls", "xgboost", "mltools")

#Carregando o data set
df <- read.table("kc_house_data.csv", header = TRUE, sep =',' )

df <- df[!df$bedrooms==33,]

1.1 Objetivo

Temos como objetivo principal deste projeto, realizar Bussiness Analyst no conjunto de dados. E como segundo objetivo desenvolver um modelo de regressão que possa prever o valor de um imóvel com base nas variáveis apresentadas no banco de dados.

1.2 Metodologia

Primeiramente coletamos o dataset fornecido pelo Kaggles. Ao coletar tais informações percebemos que o conjunto de dados é composto por 21 variáveis (tamanho da casa, quantidade de quartos, banheiros, entre outros) e o valor vendido da casa entre o período de Maio de 2014 a Maio de 2015, variável dependendo, ou seja, a variável que queremos que o nosso modelo nos apresente.

2 Explorando Dataset

Ao analisar o conjunto de dados, podemos observar que ele é composto por 21.613 linhas e 21 colunas, sendo cada coluna uma característica da casa vendida. Abaixo podemos observar cada coluna e o seu significado:

• Price: Preço da Casa
• Id: uma identificação individual única para cada casa vendida;
• date: Data realizada de venda da casa;
• Bedrooms: Quantidade de quarta que a casa possui;
• Bathrooms: Quantidade de banheiros que a casa possui;
• sqft_living: Metragem quadrados da casa;
• sqlt_lot: Metragem quadrada do lote;
• floors: Quantidade de andares que a casa possui;
• Waterfront: Se a casa possui vista para o mar;
• View: Visualização da casa;
• Condition: Condição geral da casa;
• Grade: Classificação geral de atributos de unidade residual;
• sqft_above: Metragem quadrada além do porão;
• sqft_basement: metragem quadrada do basamento;
• yr_built: Ano de construção;
• yr renovated: Ano em que a casa foi reformada;
• zip code: Cep da casa;
• lat: Latitude;
• long: Longitude;
• sqlt_living15: Área da sala de estar em 2015;
• sqlt_lot15: Outra metragem de área;

2.1 Regressão Linear simples:

Regressão Linear é um nome dado pela econometria ou estatística, a uma equação que se estima a condicional (valor esperado) de uma variável Y (Valor da casa), recedendo o valor de uma outra variável, valor x (metragem do terreno).

Em geral a regressão linear tem como objetivo, traçar uma reta através dos dados em um diagrama de dispersão. Vejamos o gráfico abaixo:

ggplot(df, aes(x= sqft_living, y =price/100000)) +
  geom_point(color = 'lightblue') + 
  geom_smooth(method = lm, size=1.2, color = "red") +
  theme_classic() +
    labs(title = "Preço de Casa X Área da casa",
       x = "Sqft",
       y = "Preço")

Como podemos observar acima, uma linear simples, ou seja apenas uma únia variável (área do terreno) não consegue representar efetivamente o valor de uma casa. É preciso acrescentar outras variáveis que agregam valor para o preço da casa, mas antes de realizar o Feature Selection iremos analisar as outras variáveis.

No gráfico abaixo iremos observar a correlação entre as variáveis:

corH <- cor(df[,c(3:21)])

corrplot::corrplot(corH, type="full", method = "circle", main="Correlation")

par(mfrow=c(3,3))
for(i in 3:17){
  plot(df[,i], df$price, main=names(df[i]), ylab=names(df$price), xlab="", col='steelblue')
}

par(mfrow=c(1,1))

for(i in c(4,5,8,9,10,11,12)){
  boxplot(df[,1]~df[,i], xlab='', main=names(df[i]), col=c("blue"))
}

df1 <- df

df1$floors <- mapvalues(df$floors,
                        from = c("1","1.5","2","2.5","3","3.5"),
                        to =c("Casa com 1 andar",
                              "Casa com 1 andar e 1/2",
                              "Casa com 2 andares",
                              "Casa com 2 andares e 1/2",
                              "Casa com 3 andares",
                              "Casa com 3 andares e 1/2"
                              ))


ggplot(df1, aes(x= sqft_living/1000, y =price/100000)) +
  geom_point(color = 'yellow') + 
  geom_smooth(span = 0.5, size=1.2, color = "red") +
  theme_classic() +
    labs(title = "Preço de Casa X Área da casa",
       x = "Metragem quadrada da casa",
       y = "Preço")+
  facet_wrap(~floors)

Agora iremos comparar a quantidade de quartos:

#df$floors <- factor(as.character(df$floors),
#                    levels = c("1","1.5","2","2.5","3","3.5"))

ggplot( df1, aes(x = as.factor(bedrooms), y = price/100000, color = floors)) +
geom_point() + 
    scale_colour_brewer(palette="Greens", name="Andar") +
  theme_classic() +
    labs(title = "Preço de Casa X Quantidade de quartos",
       x = "Quantidade de Quartos",
       y = "Preço") + 
  ggplot2::annotate("text", x=1, y=20, label= "40.95", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=2, y=18, label= "31.76", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=3, y=38, label= "40.14", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=4, y=43, label= "46.62", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=5, y=50, label= "63.54", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=6, y=75, label= "78.66", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=7, y=80, label= "82.55", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=8, y=40, label= "95.12", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=9, y=37, label= "110.51", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=10, y=18, label= "89.40", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=11, y=15, label= "81.93", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=12, y=15, label= "52.00", cex= 3, colour = "black")

#df %>%
#  dplyr::group_by(bedrooms) %>%
#  dplyr::summarise(media = round(mean(price/10000), digits = 2))

O gráfico acima nos mostra que as casa com maiores preços possuem entre 04 a 06 quartos. Os valores apresentado no gráfico é a média de preço da casa em relação aos quartos.

Porém em relação a quantidade de banheiros, conforme podemos demonstrar no gráfico abaixo,

df1 <- df1 %>%
  separate(bathrooms, into = c("bathrooms","not"), sep = "\\.")

df1 <- df1[!df1$bathrooms==0,]



ggplot( df1, aes(x = as.factor(bathrooms), y = price/100000, color = floors)) +
geom_point() + 
    scale_colour_brewer(palette="Greens", name="Andar") +
  theme_classic() +
    labs(title = "Preço de Casa X Quantidade de banheiros",
       x = "Quantidade de banheiros",
       y = "Preço") + 
  ggplot2::annotate("text", x=1, y=38, label= "39.74", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=2, y=35, label= "54.42", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=3, y=50, label= "88.50", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=4, y=43, label= "139.88", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=5, y=65, label= "195.88", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=6, y=60, label= "272.97", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=7, y=75, label= "366.75", cex= 3, colour = "black")+ 
  ggplot2::annotate("text", x=8, y=80, label= "499", cex= 3, colour = "black")

#Removendo "T000000" da coluna date

df1$date <- gsub("T000000","",df1$date)


# Excluindo a clununa id
df1$id <- NULL

# Separando a coluna em ano mês e dia

df1 <- df1 %>%
  mutate(ano = str_sub(date,end = 4),
         mes = str_sub(str_sub(df1$date, start = 5),end = 2),
         dia = str_sub(df1$date, start = 7))

df1$not <- NULL

#Criando uma coluna de mês com os nomes dos meses

df1$month <- df1$mes

df1$month <- mapvalues(df1$month,
                       from = c("01","02","03","04","05","06","07","08","09","10","11","12"),
                       to = c("Janeiro","Feveveiro","Março","Abril","Maio","Junho","Julho","Agosto","Setembro","Outubro","Novembro","Dezembro"))

df1$month <- factor(as.character(df1$month),
                   levels = c("Janeiro","Feveveiro","Março","Abril","Maio","Junho","Julho","Agosto","Setembro","Outubro","Novembro","Dezembro"))


ggplot(df1, aes(x = month)) +
  geom_bar(stat="count", fill = c("#abd4ff","#abd4ff","#abd4ff",
                                  "#b3db00","#b3db00","#b3db00",
                                  "#bd0000","#bd0000","#bd0000",
                                  "#b56700","#b56700","#b56700"), color = "black", size = 0.1)+
  theme_classic() +
  geom_text(aes(label = paste0(round(prop.table(..count..)*100,2),"%"), y= ..count.. ), stat= "count", vjust=3, color="black",size=3) +
  labs(title = "Quantidade de casa vendida por mês",x = "Mês",y = "Quantidade de Casa Vendidas")+
theme(axis.text.x = element_text(angle = 45, hjust = 1))+ 
  ggplot2::annotate("text", x=2, y=200, label= "Inverno", cex= 6, colour = "black") + 
  ggplot2::annotate("text", x=5, y=200, label= "Primavera", cex= 6, colour = "black") + 
  ggplot2::annotate("text", x=8, y=200, label= "Verão", cex= 6, colour = "black") + 
  ggplot2::annotate("text", x=11, y=200, label= "Outono", cex= 6, colour = "black") 

3 Similaridade

Para averiguar esse resultado, consideremos uma análise multivariada de modo a analisar disposição de locação em cada cidade em relação às variáveis preço de aluguel, quantidade de quartos, quantidade de camas, banheiros, quantidade de hóspedes e avaliação do imóvel. Em função de problemas de escala, vamos considerar a padronização em scale. O objetivo aqui é verificar se há imóvel que destacam-se em relação a demais locações e, caso existam, quais são estes. Vejamos:

df2 <- df
df2$bathrooms <- as.numeric(df2$bathrooms)

pca <- prcomp(df2[,c(3,4,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21)], scale.=T)
            
summary(pca)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.2459 1.5538 1.3500 1.17092 1.09809 1.00190 0.95175
## Proportion of Variance 0.2802 0.1341 0.1013 0.07617 0.06699 0.05577 0.05032
## Cumulative Proportion  0.2802 0.4144 0.5156 0.59179 0.65878 0.71454 0.76487
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.91427 0.80327 0.74039 0.71449 0.64729 0.57471 0.53951
## Proportion of Variance 0.04644 0.03585 0.03045 0.02836 0.02328 0.01835 0.01617
## Cumulative Proportion  0.81130 0.84715 0.87760 0.90597 0.92924 0.94759 0.96376
##                           PC15    PC16    PC17      PC18
## Standard deviation     0.51044 0.45469 0.43009 2.958e-15
## Proportion of Variance 0.01447 0.01149 0.01028 0.000e+00
## Cumulative Proportion  0.97824 0.98972 1.00000 1.000e+00

df2$waterfront <- mapvalues(df2$waterfront,
                           from = c("0","1"),
                           to = c("Não possui vista para o Mar", "Possui vista para o Mar"))

group = as.factor(df2$waterfront)

fviz_pca_ind(pca,
             col.ind = group, 
             palette = c("#00AFBB","#FC4E07" ),
             addEllipses = F,
             ellipse.type = F,
             legend.title = "Vista para o mar")

No No contexto multivariado, temos que as duas variáveis latentes criadas resumem os 41,30% dos dados do dataset. Note que com isso, podemos notar um comportamento bem similar das locações de casas em King country - EUA que tendem a estarem próximos no gráfico acima.

3.1 Feature Selection

A seleção de Atributos, ou o Feature Selection tem como objetivo a simpleficação do modelo, para facilitar sua interpretação, dredução do tempo de treinamento do modelo e melhoria da generalização do modelo, evitando overfitting. Utilizaremos a tecnicas de feaure selection para automitizar a seleção de variáveis com maior potencialpara variáveis preditoras. Sendo uma espécia de filtro, que remove do seu dataset as variáveis qwue não serão úteis para a criação do modelo preditivo. Tem como principal objetivo a criação de um modelo preditoco com a maior precisão possível e que seja generalizável. As técnicas de Feaure Selection basicamente calculam o nível de signifiância de cada variável e eliminam aquelas com significância mais baixa.

df$Renovated <- ifelse(df$yr_renovated == 0, 0, 1)
df$Basement <- ifelse(df$sqft_basement == 0, 0, 1)
df$HouseAge <- year(Sys.time()) - df$yr_built
df$Date <- as.POSIXct(str_replace_all(df$date, "T000000", ""), format = "%Y%m%d")
head(df[,3:25])

O número de casas construídas tem aumentado constantemente com o tempo com algumas irregularidades, principalmente em tempos de dificuldades financeiras (1930, 2008-12)

options(repr.plot.width = 10, repr.plot.height = 5)

a1 <- ggplot(df, aes(yr_built, price/100000)) +
geom_point( color = "#08519c") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
theme_classic() +
  labs( y= "Price") +
theme(text = element_text(face = "bold"),
     axis.text.x = element_blank(),
     axis.title.x = element_blank())

a2 <- ggplot(df, aes(yr_built, price/10000)) +
geom_smooth(se = FALSE, colour = "dodgerblue3") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
theme_classic() +
labs( x= "Ano de Construção",
      y = "Price") + 
theme(text = element_text(face = "bold"))

output1 <- grid.arrange(a1, a2)

Existem 2 casas particulares que são totalmente caras até os anos 40, mas o preço parece ser constante até os anos 70 e a partir daí geralmente uma tendência ascendente. Em geral, Price mostra muita assimetria correta, como é o caso de entidades não físicas do mundo real.

Relação entre insumos discretos e preço da casa

discreteVars <- c("bedrooms", "bathrooms", "floors", "waterfront", "view", "condition", "grade", "Renovated", "Basement")


options(repr.plot.width = 14, repr.plot.height = 10)

p <- ggpairs(df[,c("price", discreteVars)]) +
theme_minimal() +
theme(text = element_text(face = "bold"),
     axis.text.x = element_text(angle = 90),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank())
p

3.2 Modelo Preditivo

distribution <- as.data.frame(t(sapply(df[,3:21], quantile)))
distribution$Mean <- sapply(df[,3:21], mean)
distribution$SD <- sapply(df[,3:21], sd)
round(distribution, 2)

### Lets create discrete variables
discreteVars <- c("bedrooms", "floors", "waterfront", "view", "condition", "grade","Renovated", "Basement")
df[,discreteVars] <- lapply(df[,discreteVars], as.factor)
t(sapply(df, class))

##      id        date        price     bedrooms bathrooms sqft_living sqft_lot 
## [1,] "numeric" "character" "numeric" "factor" "numeric" "integer"   "integer"
##      floors   waterfront view     condition grade    sqft_above sqft_basement
## [1,] "factor" "factor"   "factor" "factor"  "factor" "integer"  "integer"    
##      yr_built  yr_renovated zipcode   lat       long      sqft_living15
## [1,] "integer" "integer"    "integer" "numeric" "numeric" "integer"    
##      sqft_lot15 Renovated Basement HouseAge  Date       
## [1,] "integer"  "factor"  "factor" "integer" Character,2

Vamos padronizar as entradas quantitativas e transformar a saída usando a transformação de log.

df1 <- df
df$price <- log(df$price)

preprocessParams <- preProcess(df[,!(colnames(df) %in% c(discreteVars, "Price"))], 
                               method=c("center", "scale"))
df[,!(colnames(df) %in% c(discreteVars, "Price"))] <- 
predict(preprocessParams, df[,!(colnames(df) %in% c(discreteVars, "Price"))])
head(df[,-c(1:2)])

Agora estamos prontos para fazer a primeira rodada do processo de construção de modelo.

options(repr.plot.width = 7, repr.plot.height = 5)
set.seed(123)
boruta.train <- Boruta(price~., data = df[,3:23], doTrace = 2, maxRuns = 200, getImp = getImpRfZ, ntree = 500)

plotImpHistory(boruta.train)

Verde e azul representam, respectivamente, recursos confirmados e sombreados

boruta_results <- attStats(boruta.train)
boruta_results <- rownames_to_column(boruta_results)
boruta_results <- arrange(boruta_results, desc(normHits), desc(medianImp))
borutaInputs <- getSelectedAttributes(boruta.train, withTentative = FALSE)
boruta_results

De acordo com o algoritmo de Boruta, a idade da casa é a característica mais importante na determinação do preço da casa.

modelling_df <- df[,c(borutaInputs, "price")]

set.seed(100)
p <- 0.8
train_df <- modelling_df[sample(1:nrow(modelling_df), p*nrow(modelling_df)),]

test_df <- modelling_df[-sample(1:nrow(modelling_df), p*nrow(modelling_df)),]

linear <- lm(price~., train_df)

fit_lin <- lm(price~1, data = train_df)

big <- formula(lm(price~.,data = train_df))

stepwise <- stepAIC(fit_lin, steps = floor(0.6*nrow(train_df)), scope = big,direction="both",trace = FALSE)

comp <- min(5,ncol(train_df))

PLS <- plsr(price~., data = train_df, comps = comp)

rF <- randomForest(price~., data = train_df, keep.forest = TRUE, keep.inbag = TRUE)

model_imp <- as.data.frame(rF$importance, stringsAsFactors = FALSE)

model_imp$IncNodePurity <- round(model_imp$IncNodePurity, 2)

model_imp <- rownames_to_column(model_imp)

model_imp <- arrange(model_imp, desc(IncNodePurity))

model_imp

3.3 Previsões

train_model_data <- data.frame(matrix(nrow = nrow(train_df), ncol = 5))

colnames(train_model_data) <- c("Actual", "Linear", "Stepwise", "randomForest", "XGBoost")

train_model_data$Actual <- train_df$price

train_model_data[,2] <- predict(linear, train_df)

## Warning in predict.lm(linear, train_df): prediction from a rank-deficient fit
## may be misleading

train_model_data[,"Stepwise"] <- predict(stepwise, train_df)

train_model_data[,"randomForest"] <- predict(rF, train_df)

#train_model_data[,"XGBoost"] <- predict(xgb_tune, test_df1)
train_model_data <- sapply(train_model_data, exp)

#########################################################################################

test_model_data <- data.frame(matrix(nrow = nrow(test_df), ncol = 5))

colnames(test_model_data) <- c("Actual", "Linear", "Stepwise", "randomForest")

test_model_data$Actual <- test_df$price

test_model_data[,"Linear"] <- predict(linear, test_df)

## Warning in predict.lm(linear, test_df): prediction from a rank-deficient fit may
## be misleading

test_model_data[,"Stepwise"] <- predict(stepwise, test_df)

test_model_data[,"randomForest"] <- predict(rF, test_df)

#test_model_data[,"XGBoost"] <- predict(xgb_tune, test_df1)
test_model_data <- sapply(test_model_data, exp)

err_train <- data.frame(matrix(nrow = 5, ncol = 4))

colnames(err_train) <- c("Actual", "Linear", "Stepwise", "randomForest")

rownames(err_train) <- c("RMSE","MAPE","MAE","maxError","rSquare")

err_test <- data.frame(matrix(nrow = 5, ncol = 4))

colnames(err_test) <- c("Actual", "Linear", "Stepwise", "randomForest")

rownames(err_test) <- c("RMSE","MAPE","MAE","maxError","rSquare")
errcalc <- function(Actual, Predicted){
  err <- c()
  res <- abs((Actual-Predicted)*100/Actual)
  perc_error <- res
  err[1] <- rmse(Actual, Predicted, na.rm=TRUE)
  err[2] <- sum(res)/length(res)
  err[3] <- sum(abs((Actual-Predicted)/length(res)))
  err[4] <- max(abs(Actual-Predicted))
  err[5] <- rSquared(Actual, Actual - Predicted)
  return(err)
}

err_train[,1] <- errcalc(train_model_data[,1],train_model_data[,2])
err_train[,2] <- errcalc(train_model_data[,1],train_model_data[,3])
err_train[,3] <- errcalc(train_model_data[,1],train_model_data[,4])

err_train <- as.data.frame(apply(err_train, 2, function(x) round(x, 2)), stringsAsFactors = FALSE)
#########################################################################################
                                 
err_test[,1] <- errcalc(test_model_data[,1],test_model_data[,2])
err_test[,2] <- errcalc(test_model_data[,1],test_model_data[,3])
err_test[,3] <- errcalc(test_model_data[,1],test_model_data[,4])
err_test <- as.data.frame(apply(err_test, 2, function(x) round(x, 2)), stringsAsFactors = FALSE)

err_train

err_test

---

E-mail: consultoriaterra@hotmail.com

site: www.rodolfoterra.com

GitHub: https://github.com/rodolffoterra