Esse documento mostra a aplicação de modelos de regressão para a precificação de imóveis em Curitiba. Os dados são referêntes a imoveis de dois bairros de Curitiba: Cristo Rei e Rebouças. Os dados foram retirados do site http://www.imovelweb.com.br/ usando recursos para Web Scraping do pacote XML.

Carregando Pacotes

#-----------------------------------------------------------------------
# Load packages.

library(lattice)
library(latticeExtra)
library(asbio) # press().
library(car)   # residualPlots() além de outras.

Carregando os Dados

Os dados foram salvos em arquivos CSV. Eles contém informações que indenticam o imóvel (id), o tamanho (tam, \(m^2\)), o número de quartos (quar) e banheiros (banh) e o preço de venda do imóvel (valor, R$).

#-----------------------------------------------------------------------
# Dados de inventário florestal.

cr <- read.csv2("data/cristo-rei.csv",
                header = TRUE, stringsAsFactors = FALSE)
rb <- read.csv2("data/reboucas.csv",
                header = TRUE, stringsAsFactors = FALSE)
str(cr)
## 'data.frame':    861 obs. of  5 variables:
##  $ id   : chr  "aviso-2927564861" "aviso-2924957329" "aviso-2925955953" "aviso-2927461774" ...
##  $ tam  : int  133 32 32 75 32 0 51 131 69 60 ...
##  $ quar : int  3 1 1 3 1 2 2 2 2 2 ...
##  $ banh : int  3 1 1 2 1 1 2 1 2 1 ...
##  $ valor: num  1100000 208492 149000 650000 260000 ...
da <- rbind(cbind(bair = "CRei", cr),
            cbind(bair = "Rebo", rb))
str(da)
## 'data.frame':    1582 obs. of  6 variables:
##  $ bair : Factor w/ 2 levels "CRei","Rebo": 1 1 1 1 1 1 1 1 1 1 ...
##  $ id   : chr  "aviso-2927564861" "aviso-2924957329" "aviso-2925955953" "aviso-2927461774" ...
##  $ tam  : int  133 32 32 75 32 0 51 131 69 60 ...
##  $ quar : int  3 1 1 3 1 2 2 2 2 2 ...
##  $ banh : int  3 1 1 2 1 1 2 1 2 1 ...
##  $ valor: num  1100000 208492 149000 650000 260000 ...
da <- subset(da, tam > 10)

#-----------------------------------------------------------------------
# Visualização.

xyplot(valor ~ tam, groups = bair,
       data = da,
       auto.key = TRUE,
       type = c("p", "smooth", "g"),
       xlab = expression("Tamanho do imóvel" ~ (m^2)),
       ylab = "Valor de venda (R$)")

xyplot(valor ~ tam, groups = bair,
       data = da,
       auto.key = TRUE,
       scales = list(x = list(log = 10), y = list(log = 10)),
       type = c("p", "smooth", "g"),
       xlab = expression("Tamanho do imóvel" ~ (log ~ m^2)),
       ylab = "Valor de venda (R$)",
       xscale.components = xscale.components.log10.3,
       yscale.components = yscale.components.log10.3)

Na escala log-log os dados apresentam mais compatibilidade com as suposições do modelo linear usual: relação linear, variância homogênea, erros normais. Sendo assim, a análise será aplicada aos dados transformados.

Ajuste do Modelo de Precificação

#-----------------------------------------------------------------------

# Dividir o preço por 1000 e criar o log das variáveis.
da <- within(da, {
    valor <- valor/1000
    lvalor <- log10(valor)
    ltam <- log10(tam)
})

str(da)
## 'data.frame':    1568 obs. of  8 variables:
##  $ bair  : Factor w/ 2 levels "CRei","Rebo": 1 1 1 1 1 1 1 1 1 1 ...
##  $ id    : chr  "aviso-2927564861" "aviso-2924957329" "aviso-2925955953" "aviso-2927461774" ...
##  $ tam   : int  133 32 32 75 32 51 131 69 60 92 ...
##  $ quar  : int  3 1 1 3 1 2 2 2 2 3 ...
##  $ banh  : int  3 1 1 2 1 2 1 2 1 2 ...
##  $ valor : num  1100 208 149 650 260 ...
##  $ ltam  : num  2.12 1.51 1.51 1.88 1.51 ...
##  $ lvalor: num  3.04 2.32 2.17 2.81 2.41 ...
splom(da[-(1:2)], groups = da$bair, cex = 0.2)

# Modelo com as variáveis na escala original.
m0 <- lm(valor ~ bair * tam, data = da)

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

A análise feita acima apenas confirma que na escala original os pressupostos do modelo não são verificados.

# Modelo com as variáveis na escala log.
m0 <- lm(lvalor ~ bair * ltam, data = da)

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Modelo com inclusão de quartos e banheiros.
m1 <- lm(lvalor ~ bair * ltam + quar + banh, data = da)

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m1)

layout(1)

O modelo que contém as preditoras quartos e banheiros (m1) exibe um padrão mais compatível com os pressupostos do que o modelo mais simples (m0). Quando um modelo tem muitas variáveis preditoras, é mais difícil reconhecer se desvios de pressupostos são causados pela má especificação de algumas delas. Nesse caso, gráfico com os resíduos parciais são úteis.

# Resíduos parciais com teste para inclusão de termo quadrático.
residualPlots(m1)

##            Test stat Pr(>|t|)
## bair              NA       NA
## ltam           3.767    0.000
## quar           1.827    0.068
## banh          -0.441    0.659
## Tukey test     3.326    0.001
m2 <- lm(lvalor ~ bair * (ltam + I(ltam^2)) + quar + banh, data = da)
anova(m1, m2)
## Analysis of Variance Table
## 
## Model 1: lvalor ~ bair * ltam + quar + banh
## Model 2: lvalor ~ bair * (ltam + I(ltam^2)) + quar + banh
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   1562 16.806                                  
## 2   1560 16.646  2   0.15997 7.4959 0.0005756 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas e medidas de ajuste.
summary(m2)
## 
## Call:
## lm(formula = lvalor ~ bair * (ltam + I(ltam^2)) + quar + banh, 
##     data = da)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41785 -0.05539  0.00363  0.06208  0.38741 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.615985   0.217597   7.426 1.83e-13 ***
## bairRebo           -0.035777   0.261605  -0.137   0.8912    
## ltam                0.145113   0.230634   0.629   0.5293    
## I(ltam^2)           0.183939   0.060240   3.053   0.0023 ** 
## quar               -0.009307   0.005903  -1.577   0.1151    
## banh                0.041622   0.004210   9.887  < 2e-16 ***
## bairRebo:ltam       0.180858   0.281844   0.642   0.5212    
## bairRebo:I(ltam^2) -0.067261   0.075199  -0.894   0.3712    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1033 on 1560 degrees of freedom
## Multiple R-squared:  0.8146, Adjusted R-squared:  0.8138 
## F-statistic: 979.5 on 7 and 1560 DF,  p-value: < 2.2e-16
# Distribuição de frequência de quarto e banheiro.
xtabs(~quar + banh, data = da)
##     banh
## quar   1   2   3   4   5   6
##    1 483  31   0   0   0   0
##    2 140 261  36   1   0   0
##    3  70 255 134  57  20   0
##    4   4   5  18  23  26   4

O teste para inclusão do termos quadrático foi significativo (m2). O teste é realizado ao exibir os resíduos parciais.

#-----------------------------------------------------------------------
# Predição.

pred <- with(da, expand.grid(bair = levels(bair),
                             quar = 2,
                             banh = 2,
                             ltam = seq(min(ltam), max(ltam),
                                        length.out = 20)))

fit <- predict(m2, newdata = pred, interval = "confidence")
pred <- cbind(pred, fit)
str(pred)
## 'data.frame':    40 obs. of  7 variables:
##  $ bair: Factor w/ 2 levels "CRei","Rebo": 1 2 1 2 1 2 1 2 1 2 ...
##  $ quar: num  2 2 2 2 2 2 2 2 2 2 ...
##  $ banh: num  2 2 2 2 2 2 2 2 2 2 ...
##  $ ltam: num  1.32 1.32 1.4 1.4 1.48 ...
##  $ fit : num  2.19 2.28 2.24 2.33 2.3 ...
##  $ lwr : num  2.15 2.25 2.21 2.31 2.27 ...
##  $ upr : num  2.24 2.31 2.28 2.35 2.32 ...
source(paste0("https://raw.githubusercontent.com/",
              "walmes/wzRfun/master/R/panel.cbH.R"))

p1 <- xyplot(fit ~ ltam, data = pred, type = "l",
             ly = pred$lwr, uy = pred$upr,
             xlab = expression("Tamanho do imóvel" ~ (log ~ m^2)),
             ylab = "Valor de venda (R$)",
             cty = "bands",
             groups = bair, alpha = 0.5,
             prepanel = prepanel.cbH,
             panel.groups = panel.cbH,
             panel = panel.superpose)
p1

xyplot(valor ~ tam, groups = bair,
       data = subset(da, quar == 2 & banh == 2),
       auto.key = TRUE,
       scales = list(x = list(log = 10), y = list(log = 10)),
       xlab = expression("Tamanho do imóvel" ~ (m^2)),
       ylab = "Valor de venda (R$)",
       xscale.components = xscale.components.logpower,
       yscale.components = yscale.components.logpower) +
    as.layer(p1)

Informações da Sessão

## Updated in 2016-10-10.
## 
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] tcltk     stats     graphics  grDevices utils     datasets 
## [7] base     
## 
## other attached packages:
## [1] car_2.1-2           asbio_1.3-4         latticeExtra_0.6-28
## [4] RColorBrewer_1.1-2  lattice_0.20-33     rmarkdown_1.0      
## [7] knitr_1.14         
## 
## loaded via a namespace (and not attached):
##  [1] pixmap_0.4-11        Rcpp_0.12.7          magrittr_1.5        
##  [4] splines_3.3.1        MASS_7.3-45          scatterplot3d_0.3-37
##  [7] minqa_1.2.4          stringr_1.0.0        tools_3.3.1         
## [10] parallel_3.3.1       nnet_7.3-12          pbkrtest_0.4-6      
## [13] grid_3.3.1           nlme_3.1-128         mgcv_1.8-13         
## [16] quantreg_5.26        multcompView_0.1-7   plotrix_3.6-2       
## [19] MatrixModels_0.4-1   htmltools_0.3.5      lme4_1.1-12         
## [22] yaml_2.1.13          digest_0.6.10        Matrix_1.2-6        
## [25] nloptr_1.0.4         formatR_1.4          deSolve_1.13        
## [28] evaluate_0.9         stringi_1.1.1        methods_3.3.1       
## [31] SparseM_1.7          mvtnorm_1.0-5