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.
#-----------------------------------------------------------------------
# Load packages.
library(lattice)
library(latticeExtra)
library(asbio) # press().
library(car) # residualPlots() além de outras.
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.
#-----------------------------------------------------------------------
# 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)
## 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