1 Introdução

Para realização deste trabalho, iremos utilizar um banco de dados no qual os dados mobiliários da cidade de São Paulo em diversas variáveis. Os dados possuem milhares de observações, porém selecionamos apenas algumas das variáveis (as qual nos interessa) para montar o modelo, inicialmente usamos as seguintes variáveis:

  • Logradouro:lugar onde está o imóvel (rua,hortos,passeios);
  • Bairro:bairro de São Paulo que se localiza o imóvel;
  • Cep: código de endereçameno postal;
  • Tipo de Imóvel:se referindo a apartamento,casas, etc;
  • Area útil:área em metros² do imóvel;
  • Banheiros:quantidade de banheiros no imóvel;
  • Vagas garagem:quantidade de vagas na garagem;
  • Preço venda:preço de venda do imóvel;
  • Preço aluguel:preço do aluguel do imóvel;


Além dessas variáveis, há diversas outras porém não muito úteis para tal análise. Pela quantidade e variedade dos dados, escolhermos fazer uma análise apenas dos apartamentos para alugar, desta forma, usaremos:

  1. A variável resposta: preço do aluguel
  2. As variáveis explicativas: área útil,banheiros,suites,quartos,vagas de garagem,preço do aluguel e Zonas(a qual criaremos utilizando o Cep e Bairros)

Os dados analisados estão disponíveis para download no Kaggle

2 Manipulação dos dados

Para realização deste trabalho, usamos os seguintes pacotes:

library(yaml)
library(rmarkdown)
library(dplyr)                                
library(car)                                
library(rstatix)                                
library(emmeans)
library(ggplot2)
library(knitr)
library(kableExtra)
library(htmltools)
library(dplyr)                                
library(tidyverse)                                
library(GGally)                                
library(corrr)
library(nortest)
library(readxl)
library(cowplot)
library(faraway)
library(psych)
library(olsrr)
library(equatiomatic)

2.1 Limpando os dados

Realizamos a limpeza dos dados principalmente usando o pacote Dplyr o qual ajudou muito a deixar mais claro as variáveis para o uso no modelo.

BaseOriginal<-read_excel("housing_sp_city1.xls")
BaseTop<-BaseOriginal %>% select(preco_aluguel,bairro,cep,tipo_imovel,area_util,banheiros,suites,quartos,vagas_garagem,periodicidade)
BaseTop<-BaseTop %>% filter(periodicidade=="MONTHLY",tipo_imovel=="Apartamento")
kable(head(BaseTop,10), col.names = c("Preço Aluguel", "Bairro", "Cep", "Tipo de imóvel","Área(m²)","Banheiros","Suítes","Quartos","Vagas Garagem","Periodicidade"))
Preço Aluguel Bairro Cep Tipo de imóvel Área(m²) Banheiros Suítes Quartos Vagas Garagem Periodicidade
12299 Indianópolis 4081010 Apartamento 240 3 2 4 3 MONTHLY
1941 Cidade São Francisco 5351025 Apartamento 72 2 1 3 2 MONTHLY
0 Parque das Paineiras 3692060 Apartamento 35 1 0 2 0 MONTHLY
NA Perdizes 5018011 Apartamento 123 3 1 3 2 MONTHLY
0 Artur Alvim 3692060 Apartamento 42 1 0 2 0 MONTHLY
0 Jardim Umarizal 5756240 Apartamento 84 1 0 2 1 MONTHLY
2485 Pinheiros 3178200 Apartamento 30 1 1 1 0 MONTHLY
3290 Penha 3632020 Apartamento 154 5 3 3 3 MONTHLY
1014 Vila Maria 2127001 Apartamento 64 1 0 2 1 MONTHLY
5950 Vila Romana 5051030 Apartamento 180 5 2 4 3 MONTHLY

Com a tabela acima, temos uma visão geral do banco, porém uma série de mudanças vão acontecer agora para ajudar a deixar o banco mais limpo possível para usarmos os dados no ajuste, um adendo é que, por causa de alguns Ceps mal formatados ou faltantes, fizemos um trabalho manual pelo excel para corrigir os números e também achar Ceps em brancos (usamos os logradouros e bairros), algumas das manipulações dos dados apresentamos abaixo: 
BaseTop=BaseTop %>% filter(preco_aluguel!=0)
BaseTop = BaseTop %>% mutate_at(vars(cep,area_util,banheiros,suites,quartos,vagas_garagem,preco_aluguel),parse_number)#passando os valores <chr> para numérico
BaseTop=BaseTop %>% mutate(Zonas=case_when(cep>=1000000 & cep<=1599000 ~ "Centro",
                                      cep>=2000000 & cep<=2999000 ~ "Norte",
                                      cep>=3000000 & cep<=3999000 ~ "Leste",
                                      cep>=8000000 & cep<=8499000 ~ "Leste",
                                      cep>=4000000 & cep<=4999000 ~ "Sul",
                                      cep>=5000000 & cep<=5999000 ~ "Oeste"))#Categorizando os ceps pelas Zonas da cidade de São paulo
BaseTop=BaseTop %>% select(preco_aluguel,area_util,banheiros,suites,quartos,vagas_garagem,Zonas)#Retirando as variáveis agora desnecessárias para o modelo
BaseTop=na.omit(BaseTop)#Omitindo os dados NA
BaseTop=BaseTop %>% mutate(Z1=case_when(Zonas=="Centro" ~ 1,
                                    Zonas!="Centro" ~ 0),
                       Z2=case_when(Zonas=="Leste" ~ 1,
                                    Zonas!="Leste" ~ 0),
                       Z3=case_when(Zonas=="Sul" ~ 1,
                                    Zonas!="Sul" ~ 0),
                       Z4=case_when(Zonas=="Oeste" ~ 1,
                                    Zonas!="Oeste" ~ 0),
                       Z5=case_when(Zonas=="Norte" ~ 1,
                                    Zonas!="Norte" ~ 0))#Separar as Zonas em dummy

2.2 Visualizando o Banco Final


Por fim, após essa limpeza nos dados, temos nossa Base final já tratada abaixo:

kable(head(BaseTop,20), col.names = c("Preço Aluguel","Área(m²)","Banheiros","Suítes","Quartos","Vagas Garagem","Zonas","Z1","Z2","Z3","Z4","Z5")) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Preço Aluguel Área(m²) Banheiros Suítes Quartos Vagas Garagem Zonas Z1 Z2 Z3 Z4 Z5
12299 240 3 2 4 3 Sul 0 0 1 0 0
1941 72 2 1 3 2 Oeste 0 0 0 1 0
2485 30 1 1 1 0 Leste 0 1 0 0 0
3290 154 5 3 3 3 Leste 0 1 0 0 0
1014 64 1 0 2 1 Norte 0 0 0 0 1
5950 180 5 2 4 3 Oeste 0 0 0 1 0
1470 40 1 0 1 1 Sul 0 0 1 0 0
1666 45 1 1 1 1 Oeste 0 0 0 1 0
2765 128 3 1 4 2 Oeste 0 0 0 1 0
3360 38 1 0 1 1 Sul 0 0 1 0 0
3080 160 5 3 3 2 Oeste 0 0 0 1 0
8001 183 5 3 3 2 Sul 0 0 1 0 0
1547 35 2 1 1 1 Oeste 0 0 0 1 0
5040 200 3 3 3 3 Oeste 0 0 0 1 0
1456 50 1 0 2 1 Norte 0 0 0 0 1
1866 40 1 0 1 1 Oeste 0 0 0 1 0
7349 210 0 3 3 4 Oeste 0 0 0 1 0
1400 55 2 0 2 1 Oeste 0 0 0 1 0
10651 230 5 3 3 4 Sul 0 0 1 0 0
2233 84 2 1 2 2 Leste 0 1 0 0 0


O banco contém 4447 observações, as quais vamos agora usar para fazer algumas análises descritivas.

3 Resultados

Agora, vamos começar apresentar os resultados obtidos, através de uma análise da base dos dados, por fim ajustando o modelo.

3.1 Análise descritiva

Para nossa análise, frizamos os variáveis criadas para as Zonas (Z1,Z2,Z3,Z4,Z5),para transformação de variáveis qualitativas em dummys, para que possamos fazer as correlações necessárias, também vamos tirar a variável Zonas, pois não será mais utilizada.

BaseTop=BaseTop %>% select(preco_aluguel,area_util,banheiros,suites,quartos,vagas_garagem,Zonas,Z1,Z2,Z3,Z4,Z5)


Vamos plotar dois gráficos de correlação para vermos quais as variáveis mais correlacionadas, para entendermos um pouco melhor como ficará o modelo.

BaseTop %>%
  correlate(quiet = TRUE) %>%
  rplot(shape=20,colors=c("red","green"),print_cor = TRUE)#Uma mera impressão inicial das correlações usando pacote Corrr

Agora vamos visualizar essa correlação novamente, porém de uma forma diferente, que passa a mesma informação mais de forma mais simples:

ggcorr(BaseTop,label=T,label_size = 3, label_round = 2, label_alpha = TRUE,
       hjust = 0.72, size = 4, color = "grey50", layout.exp = 1)#usando pacote GGally temos uma melhor forma de visualizar as correlações, nessa visualização em modo triangulo
## Warning in ggcorr(BaseTop, label = T, label_size = 3, label_round = 2,
## label_alpha = TRUE, : data in column(s) 'Zonas' are not numeric and were
## ignored


A partir dos gráficos, vamos analisar as variáveis mais correlacionadas.
Vemos que , area, suites, quartos e vagas garagem se relacionam bem com preço de aluguel, todavia também esperariamos uma boa correlação de algumas das Zonas, vamos ver o que podemos fazer.
Primeiramente vamos plotar o histograma com os preços dos imóveis com área maior que 157m²(valor da média das áreas), e acima que esse valor, para termos uma base.

h1=BaseTop %>% filter(area_util>=157)
h2=BaseTop %>% filter(area_util<157)
plot_grid(ggplot(h1, aes(x = preco_aluguel)) +
            geom_histogram(color = "white",bins = 45) +
            coord_cartesian(xlim = c(0,25000), ylim = c(0, 500))+
            labs(title = "Histograma dos preços com área maior que 157m²")+
            xlab("Preço dos alugueis em R$") + 
            ylab("Frequência"),
          ggplot(h2, aes(x = preco_aluguel)) +
            geom_histogram(color = "white",bins = 45) +
            coord_cartesian(xlim = c(0,25000), ylim = c(0, 500))+
            labs(title = "Histograma dos preços com área menor que 157m²")+
            xlab("Preço dos alugueis em R$") + 
            ylab("Frequência"),ncol=1)

Como visto no gráfico, as casas com maiores áreas tem preços maiores, podemos afirmar portanto que o aluguel está relacionado a área do imóvel.

Podemos também mostrar as outras variáveis significavelmente relacionadas com o custo (banheiros, suítes, quartos, vagas de garagem). O exemplo abaixo é do banheiro, que também, através da média, fizemos dois diagramas para comparação.


h1=BaseTop %>% filter(suites<=2)
h2=BaseTop %>% filter(suites>2)
plot_grid(ggplot(h2, aes(y = preco_aluguel)) +
            geom_boxplot(show.legend = T, alpha = 1) +
            theme_classic(base_size = 10) +
            coord_cartesian(xlim = c(-2,2), ylim = c(0, 25000))+
            xlab("Boxplot imóveis com mais de 2 suítes") + 
            ylab("Custo do aluguel"),
          ggplot(h1, aes(y = preco_aluguel)) +
            geom_boxplot(show.legend = T, alpha = 12) +
            coord_cartesian(xlim = c(-2,2), ylim = c(0, 25000))+
            theme_classic(base_size = 10) +
            xlab("Boxplot imóveis com 2 suítes ou menos") + 
            ylab("Custo aluguel"))


Fizemos o mesmo para as outras variáveis citadas acima e apresentaram resumidamente a mesma performance de correlação, portanto vamos entrar agora na variável Zonas, a qual dividimos cada um em variáveis Dummys porém o gráfico abaixo mostrando a correlação do custo com a área, por zona, já demonstra bem uma interação entre as variáveis

ggplot(BaseTop, aes(y = preco_aluguel, x = area_util
                  ,color=Zonas )) +
  geom_point(size=3,alpha=0.9)+
  theme_classic(base_size = 18) +
  theme(legend.position = "top") +
  xlab("Áreas úteis dos imóveis em m²") + 
  ylab("Custo do aluguel do imóvel")


Vemos os dados de 3 Regiões (Z1,Z3,Z4) acima que as outras, sugerindo assim uma ideia que quanto mais alto o aluguel, possivelmente será um imóvel em uma dessas regiões.

É interessante também mostrarmos o Box-plot por Zonas relacionado ao custo.

# Edita níveis do fator.
BaseTop$Zonas <- factor(BaseTop$Zonas)
boxplot(BaseTop$preco_aluguel ~ BaseTop$Zonas,
        data = BaseTop,
        xlab = "Zonas",
        ylab = "Preço de aluguel dos Apartamentos (R$)")


Infelizmente a base contou com poucos dados das Zonas do Norte e Leste, sendo assim, por mais que sugira uma leve diferença para maiores custos nas outras Zonas(principalmente o Centro como maiores preços), não podemos indicar com certeza esse fato.

4 Modelo Proposto

4.1 Breve Noção do Modelo Inicial

Vamos começar a análise de regressão linear múltipla dos nossos dados, para tentar prever o valor do aluguel baseado nas variáveis explicativas no conjunto de dados. Iniciamos essa análise com um modelo completo, contendo todas as variáveis presentes no nosso banco de dados (preço do aluguel,área útil, banheiros, suites, quartos, vagas garagem e Variáveis dummys das Zonas).

dadossemZonas=BaseTop %>% select( preco_aluguel, area_util, banheiros, suites, quartos, vagas_garagem,    Z1,    Z2,    Z3,    Z4,    Z5)
mod <- lm(preco_aluguel ~ .,dadossemZonas)

O modelo obtido foi:

\[ Y= 47.34 +23.45X_1+77.34X_2+517.77X_3-624.09X_4+635.60X_5+2305.44X_6-291.08X_7+1724.78X_8+219.83X_9 \]

Agora vamos continuar o ajuste no modelo, após a análise de variâncias do modelo, concluimos algumas variáveis não sigificativas, juntamente, utilizamos o método stepwise, vale lembrar que é importante criar um modelo com o mínimo de parâmetros possíveis a serem estimados e que explique bem o comportamento da variável resposta, portanto:

modelonovo= step(mod,direction = "both")

Desta forma, o novo modelo ajustado é composto sem as variáveis Z2 e Z5 iniciáis, que foram excluídos após o stepwise e também apresentaram sinais tanto no teste F quanto no teste T, e também retiramos a variável Banheiros, pois ela nao é tão significativa para o modelo, sendo que o R² ajustado quase inalterável após sua retirada.

\[Y= 288.16 +23.61X_1+561.14X_2-597.15X_3+642.62X_4+2093.70X_5-514.02X_6+1515.61X_7\]


O modelo inicial com todas as variáveis tinha um coeficiente de determinação ajustado igual ao modelo final até então proposto acima, reforçando a ideia que as variáveis não contribuiam em nada no nosso modelo, vamos agora testar a presença de multicolinearidade.

pairs.panels(dadossemZonas)

vif(modelonovo) 
##     area_util        suites       quartos vagas_garagem            Z1 
##      3.227298      2.649867      2.113322      3.038011      1.656742 
##            Z2            Z3 
##      1.195955      1.472215

Testamos para ver se há uma alta correlação entre as variáveis, usando a função pairs.panels, no qual geralmente uma correlação superior a 0.9 é danosa ao modelo, e também usando a função vif, que é uma medida da proporção em que a variância de um coeficiente de regressão é inflacionada pela presença de outra variável explicativa. O recomendável é um VIF < 10. Ambos estão ok com nosso modelo.

4.2 Ajustando Melhor o Modelo e Analisando

Realizamos a análise do ajuste do modelo verificando o gráfico de resíduos,normalidade,homosedasciedade e distância de cook. Verificamos a existência de observações não usuais (outliers, pontos de alavancagem, ponto de influência) nos dados do nosso modelo pelo método da Distância de Cook.

cooksd<-cooks.distance(modelonovo)
sample_size <- nrow(dadossemZonas)

maximo=max(cooksd,1.1)
plot(cooksd,pch="*",ylab="Distância de Cook",
     ylim=c(0,maximo))
abline(h=1,col="orange")
abline(h=4/sample_size,col="red")

Se constata uma grande quantidade de outliers no primeiro critério. Não se torna viável eliminar todos para um novo ajuste.

residuos=residuals(modelonovo)
par(mfrow=c(1,3))
qqnorm(residuos,xlab="Quantis Teóricos", ylab="Quantis Empíricios")
qqline(residuos)

plot(fitted(modelonovo),residuos,ylab="resíduos",
     xlab="Valores ajustados")
abline(h=0,col="red",lty=2)

plot(residuos,ylab="resíduos",xlab="Ordem")
abline(h=0,col="red",lty=2)

Verificamos o seguinte com os gráficos:

  • ** A suposição de homoscedasticidade não é verdadeira.
  • ** A suposição de independência dos erros também não parece verdadeiro.
  • ** A suposição de normalidade dos erros não é verdadeira, os resíduos não seguem uma distribuição normal, vamos realizar algumas transformações para tentar corrigir o problema

Realizamos diversos tipos de transformações com variáveis explicativas e respostas, como log,SQRT,1/X. A mais útil para nosso modelo, e que concertou um pouco o problema da normalidade, além de aumentar nosso coeficiente de determinação foi a SQRT das variáveis explicativas e resposta. Vemos abaixo um breve resumo de como ficou nosso modelo.

summary(modeloraiz)
## 
## Call:
## lm(formula = preco_aluguel ~ ., data = D1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.509  -8.549  -1.479   7.551  61.784 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    16.13342    1.10941  14.542  < 2e-16 ***
## area_util       4.32681    0.09586  45.135  < 2e-16 ***
## suites          6.73833    0.47762  14.108  < 2e-16 ***
## quartos       -13.23938    0.84794 -15.614  < 2e-16 ***
## vagas_garagem   8.01355    0.62779  12.765  < 2e-16 ***
## Z1             12.61725    0.55867  22.584  < 2e-16 ***
## Z2             -3.99177    0.78916  -5.058  4.4e-07 ***
## Z3              9.69030    0.50004  19.379  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.13 on 4439 degrees of freedom
## Multiple R-squared:  0.7051, Adjusted R-squared:  0.7046 
## F-statistic:  1516 on 7 and 4439 DF,  p-value: < 2.2e-16
residuos=residuals(modeloraiz)
par(mfrow=c(1,3))
qqnorm(residuos,xlab="Quantis Teóricos", ylab="Quantis Empíricios")
qqline(residuos)

plot(fitted(modeloraiz),residuos,ylab="resíduos",
     xlab="Valores ajustados")
abline(h=0,col="red",lty=2)

plot(residuos,ylab="resíduos",xlab="Ordem")
abline(h=0,col="red",lty=2)

5 Conclusão

Após toda essa análise (descritiva,ajustes e transformações), temos um modelo proposto que compreendemos o custo do aluguel dos apartamentos de São Paulo em relação a diversas variáveis. Além disso, obtivemos uma boa melhora com a transformação dos dados ao fim do ajuste, nosso modelo é capaz de explicar cerca de 70% da variável resposta pelas outras variáveis explicativas.
O modelo final é:

\[Y= 16.133 +4.327X_1+6.738X_2+13.239X_3+8.014X_4+12.617X_5-3.992X_6+9.690X_7\]

6 Referências Bibliográficas