Objetivo Geral Este trabalho tem como objetivo fazer uma análise exploratoria dos dados coletados no kaggle e criar um modelo que possa prever o valor do aluguel de casas nas regiões de São Paulo, Rio de Janeiro, Porto Alegre, Campinas e Belo Horizonte.

Infomormações Gerais Os dados coletados são referentes ao dataset disponibilizado no kaggle por Rubens Junior (https://www.kaggle.com/rubenssjr/brasilian-houses-to-rent), possuem informações retiradas através de Webcrawler em diversos sites com ofertas de alugueis de imoveis. Os dados contém 10962 observações e 13 vaariaveis. Possuem as seguintes colunas:

city= Cidade; area= Tamanho da casa em m2; rooms= Número de quartos; Bathroom= Número de banheiros; Parking spaces= Número de vagas em estacionamento; Floor= Número de andar; Animal= Se aceita animais ou não; Furniture= Se o local é mobiliado ou não; Hoa= Valor do condominio; rent amount (R)=Valordoaluguel; propertytax(R)= IPTU da casa; fire insurance (R$)= Seguro de incendio da casa; Total= Total pago ao mês pela casa;

Bibliotecas utilizadas

library(readr)
library(visdat)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(corrplot)
## corrplot 0.92 loaded
library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:dplyr':
## 
##     combine
library(caret)
## Carregando pacotes exigidos: lattice
library(faraway)
## 
## Attaching package: 'faraway'
## The following object is masked from 'package:lattice':
## 
##     melanoma
library(lmtest)
## Carregando pacotes exigidos: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(MVar.pt)
library(nortest)
library(fpp)
## Carregando pacotes exigidos: forecast
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Carregando pacotes exigidos: fma
## 
## Attaching package: 'fma'
## The following objects are masked from 'package:faraway':
## 
##     airpass, eggs, ozone, wheat
## Carregando pacotes exigidos: expsmooth
## Carregando pacotes exigidos: tseries
library(leaps)

Amostra dos dados

df<- read_csv("C:\\Users\\ruanr\\Downloads\\houses_to_rent_v2.csv")
## Rows: 10692 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): city, floor, animal, furniture
## dbl (9): area, rooms, bathroom, parking spaces, hoa (R$), rent amount (R$), ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(df)
## # A tibble: 6 × 13
##   city    area rooms bathroom `parking spaces` floor animal furniture `hoa (R$)`
##   <chr>  <dbl> <dbl>    <dbl>            <dbl> <chr> <chr>  <chr>          <dbl>
## 1 São P…    70     2        1                1 7     acept  furnished       2065
## 2 São P…   320     4        4                0 20    acept  not furn…       1200
## 3 Porto…    80     1        1                1 6     acept  not furn…       1000
## 4 Porto…    51     2        1                0 2     acept  not furn…        270
## 5 São P…    25     1        1                0 1     not a… not furn…          0
## 6 São P…   376     3        3                7 -     acept  not furn…          0
## # ℹ 4 more variables: `rent amount (R$)` <dbl>, `property tax (R$)` <dbl>,
## #   `fire insurance (R$)` <dbl>, `total (R$)` <dbl>

Mudança de nomes de variaveis Com o objetivo de deixar mais simples a manipulação dos dados e o entendimento, decidimos alterar o nome das variaveis.

names(df)= c("cidade","area","quartos","banheiros","n_estacionamento", "andar","animal","Furniture","condominio","valor_aluguel","IPTU","taxa_incendio","total_pago")
head(df)
## # A tibble: 6 × 13
##   cidade        area quartos banheiros n_estacionamento andar animal   Furniture
##   <chr>        <dbl>   <dbl>     <dbl>            <dbl> <chr> <chr>    <chr>    
## 1 São Paulo       70       2         1                1 7     acept    furnished
## 2 São Paulo      320       4         4                0 20    acept    not furn…
## 3 Porto Alegre    80       1         1                1 6     acept    not furn…
## 4 Porto Alegre    51       2         1                0 2     acept    not furn…
## 5 São Paulo       25       1         1                0 1     not ace… not furn…
## 6 São Paulo      376       3         3                7 -     acept    not furn…
## # ℹ 5 more variables: condominio <dbl>, valor_aluguel <dbl>, IPTU <dbl>,
## #   taxa_incendio <dbl>, total_pago <dbl>

Precisamos mudar os tipos das variáveis cidade, animal e furniture para o tipo fator. Embora realmente esses dados sejam strigs, elas representam uma variável de classificação.

df$cidade= as.factor(df$cidade)
df$animal= as.factor(df$animal)
df$Furniture= as.factor(df$Furniture)
df$banheiros=as.integer(df$banheiros)
df
## # A tibble: 10,692 × 13
##    cidade         area quartos banheiros n_estacionamento andar animal Furniture
##    <fct>         <dbl>   <dbl>     <int>            <dbl> <chr> <fct>  <fct>    
##  1 São Paulo        70       2         1                1 7     acept  furnished
##  2 São Paulo       320       4         4                0 20    acept  not furn…
##  3 Porto Alegre     80       1         1                1 6     acept  not furn…
##  4 Porto Alegre     51       2         1                0 2     acept  not furn…
##  5 São Paulo        25       1         1                0 1     not a… not furn…
##  6 São Paulo       376       3         3                7 -     acept  not furn…
##  7 Rio de Janei…    72       2         1                0 7     acept  not furn…
##  8 São Paulo       213       4         4                4 4     acept  not furn…
##  9 São Paulo       152       2         2                1 3     acept  furnished
## 10 Rio de Janei…    35       1         1                0 2     acept  furnished
## # ℹ 10,682 more rows
## # ℹ 5 more variables: condominio <dbl>, valor_aluguel <dbl>, IPTU <dbl>,
## #   taxa_incendio <dbl>, total_pago <dbl>

Tratando Dados Inicialmente iremos avaliar se o dataset possuem algum valor faltante.Abaixo é possível perceber que não possuimos nenhum dado faltante, não precisando fazer qualquer tipo de tratamento.

colSums(is.na(df))
##           cidade             area          quartos        banheiros 
##                0                0                0                0 
## n_estacionamento            andar           animal        Furniture 
##                0                0                0                0 
##       condominio    valor_aluguel             IPTU    taxa_incendio 
##                0                0                0                0 
##       total_pago 
##                0
vis_dat(df)

Análise exploratória Qual cidade possue a maior quantidade de casas para alugar ? A primeira pergunta que faremos aos dados é como as observações estão distribuidas em relação as cidades, quais cidades possuem mais observações ?

g=ggplot(data=df, aes(x=cidade)) + geom_bar(fill="lightblue",)+ coord_flip()+
  theme_classic() + labs(title = "Número de casas para alugar por cidades") + ylab("numero de casas")
ggplotly(g)

O gráfico de barras acima mostra que a grande maioria dos número de alugueis observados foi da cidade de São Paulo, essa má distribuição dos dados talvez possa fazer o modelo criar um pequeno viés quando a cidade é são Paulo.

Como esta distribuido o valor do aluguel e qual cidade é mais cara para se morar ? Para proteger esta análise de possíveis valores extremos decidimos utilizar algumas medidas de localização e variabilidades que sejam menos influenciadas por valores extremos, tal como a média aparada e o desvio absoluto mediano (MAD).

#Media aparada
aggregate(df[,"valor_aluguel"], list(local=df$cidade),trim=0.1,mean)
##            local valor_aluguel
## 1 Belo Horizonte      2956.688
## 2       Campinas      1862.445
## 3   Porto Alegre      1879.297
## 4 Rio de Janeiro      2671.480
## 5      São Paulo      4055.008
#MAD
aggregate(df[,"valor_aluguel"], list(local=df$cidade), mad)
##            local valor_aluguel
## 1 Belo Horizonte      1779.120
## 2       Campinas       978.516
## 3   Porto Alegre       963.690
## 4 Rio de Janeiro      1497.426
## 5      São Paulo      2490.768
aggregate(df[,"valor_aluguel"], list(local=df$cidade), max)
##            local valor_aluguel
## 1 Belo Horizonte         15000
## 2       Campinas         15000
## 3   Porto Alegre         19000
## 4 Rio de Janeiro         15000
## 5      São Paulo         45000
aggregate(df[,"valor_aluguel"], list(local=df$cidade), min)
##            local valor_aluguel
## 1 Belo Horizonte           450
## 2       Campinas           500
## 3   Porto Alegre           500
## 4 Rio de Janeiro           500
## 5      São Paulo           500

Aparando a média em 10% a cidade de São Paulo possuí em média os alugueis mais caros entre os locais analisados. O desvio absoluto mediano mostra que para todas as cidades existe uma grande variabilidade nos valores de alugueis, temos valores de R$ 500,00 a R$ 45000,00 em São Paulo por exemplo, mostrando que os dados estão bem distribuidos.

Quais variáveis possuem correlação com o valor do aluguel ?┬ Para ter um entedimento melhor sobre quais as variáveis podem influênciar as variações do valor do aluguel, abaixo segue uma gráfico que mostra a correlação entre as variaveis.

colnumericas=sapply(df, is.numeric)
correlacao= cor(df[,colnumericas])
corrplot(correlacao, method = "color",addCoef.col =T)

É possivel notar uma correlação positiva média a forte entre as variaveis taxa de incendio, número de quartos, número de banheiros e número de vagas de estacionamento com a váriavel valor do aluguel. (Outra correlação muito forte é com o total pago, mas essa variavel não ajudará nada nosso modelo, pois seu valor depende do preço do aluguel) Excluiremos a variável total_pago pois ela não trás nenhuma informação útil para o modelo

df$total_pago=NULL

Analise exploratoria multivariadas das variáveis que tiveram alta ou média correlação Devido a grande quantidade de variáveis para analisar, selecionamos apenas as features que tiveram uma alta ou média correlação (Taxa de incêndio, numero de banheiros, numero de estacionamento e numero de quartos) para enteder como elas se comportam em relação ao valor do aluguel.

Como o valor do aluguel é distribuido em relação as variáveis quartos, banheiros e estacionamento ? Para avaliar essa distribuição tendo em vista que as 3 variáveis (quartos, Banheiros e Estacionamentos) são do tipo fator, a melhor forma de analisar sua relação com o valor do aluguel é com um boxplot. Abaixo segue os gráficos gerados:

lista=c("banheiros","quartos","n_estacionamento")
legenda=c("Valor do aluguel pela quantidade de banheiros",
          "Valor do alguel pela quantidade de quartos",
          "Valor do aluguel pela quantidade de n_estacionamento")

plot.boxes<-function(X,legenda){
  ggplot(df,aes_string(x=X, y="valor_aluguel",group=X))+
      geom_boxplot() + labs(title = legenda)+
    theme_classic()
        
  
}

Map(plot.boxes,lista,legenda)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## $banheiros

## 
## $quartos

## 
## $n_estacionamento

Com o Boxplot fica claro o crescimento da mediana a medida que o número de quartos,estacionamento e banheiros aumentam. Além disso, começa-se a enchegar os valores extremos no valor do aluguel ,esses outlines eram esperados pois dependendo da localização de uma casa nas cidades os valores serão totalmente diferentes mesmo que elas possuam a mesma quantidade de quartos, banheiros ou estacionamento.

Utilizando um gráfico de violino é possível entender como a distribuição dos dados do valor do aluguel se comportam a medida que temos mudanças nas variaveis categóricas. É possível identificar que para todas as variáveis a medida que seus números são maiores(quartos, Banheiros,Vagas de estacionamento) mais esparço e menos entendível fica a distribuição dos dados (muito provavelmente devido a falta de dados extremos).

lista=c("banheiros","quartos","n_estacionamento")
legenda=c("Valor do aluguel pela quantidade de banheiros",
          "Valor do alguel pela quantidade de quartos",
          "Valor do aluguel pela quantidade de n_estacionamento")
plot.boxes<-function(X,legenda){
  ggplot(df,aes_string(x=X, y="valor_aluguel",group=X))+
      geom_violin() + labs(title = legenda)
        
  
}

Map(plot.boxes,lista,legenda)
## $banheiros

## 
## $quartos
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

## 
## $n_estacionamento
## Warning: Groups with fewer than two data points have been dropped.

Como o valor do aluguel se relaciona com a taxa de incêndio ? Abaixo segue um gráfico de dispersão entre a variavel valor aluguel e taxa de incêndio, devido a grande quantidade de pontos, decidimos utilizar uma visualização de compartimentação Hexagonal.

#Geral
ggplot(df, aes(x=taxa_incendio,y=valor_aluguel)) + theme_bw() + stat_binhex(colour="blue") + scale_fill_gradient(low="white", high="black")+
labs(x="Taxa de incêndio",y="Valor do aluguel", title = "Valor do Aluguel vs Taxa de incêncio (Geral)")

aluguel e a taxa de incêndio tem um comportamento linear onde a medida que a taxa de incêndio aumenta o valor do aluguel tende a aumentar também.

Abaixo temos uma análise estratificada por cidade, é possível identificar que para todas as cidades do dataset possuimos um comportamento linear com uma correlação forte e com padrões semelhantes.

#Por cidade
ggplot(df, aes(x=taxa_incendio,y=valor_aluguel)) + theme_bw() + stat_binhex(colour="blue") + scale_fill_gradient(low="white", high="black")+
labs(x="Taxa de incêndio",y="Valor do aluguel", title = "Valor do Aluguel vs Taxa de incêncio (Por cidade)") + facet_wrap(df$cidade)

Teste de hipóteses Serão feitos dois testes de hipóteses para entender se essas questões tem influência no valor do aluguel. Para tal, os dados serão divididos em duas classes, para alugueis abaixo de 2500 mil e acima deste valor, esta separação ajudará a entender se a hipótese é significativa ou não para casas mais baratas e mais caras. Abaixo temos o valor médio aparado de alugueis para os dois grupos:

#Criando a nova variavel categorica
df$tipo="Barata"
df$tipo[df$valor_aluguel>2500]="Caro"

#Calculando a média para os dois grupos
aggregate(df[,"valor_aluguel"], list(local=df$tipo),trim=0.1,mean)
##    local valor_aluguel
## 1 Barata      1538.130
## 2   Caro      5533.138

Com base nessa divisão as casas baratas tem um valor médio de R$ 1538,130 e as casas mais caras tem um valor médio de R$ 5533,138.

Hipotese 1: O valor do aluguel aumenta caso a casa já venha mobiliada ? Para ficar mais fácil de trabalhar substituiresmos os valores da variável Furniture para 1 caso seja mobiliada e 0 caso não seja mobiliada.

muda = function(x){
  if (x=="furnished"){
    return (1)
  }else{
    return(0)
  }
}
df$Furniture= sapply(df$Furniture, muda)

Para casas Baratas Analisando os valores abaixo é possível notar que em média quando uma casa é mobiliada fica R$ 313,835 mais caro, resta analisarmos se essa diferença é estatisticamente significativa ou se pode ocorrer pelo acaso. Para isso será utilizado o teste de permutação aleatória, este teste é útil pois não precisamos nos preocupar com a distribuição dos dados, tamanho da amostra e se os dados são igualmente distribuidos.

casas_baratas=df[df$tipo=="Barata",]
media_mobiliado=mean(casas_baratas$valor_aluguel[casas_baratas$Furniture==1], trim = 0.1)
media_n_mobiliado=mean(casas_baratas$valor_aluguel[casas_baratas$Furniture==0], trim = 0.1)
c(media_n_mobiliado, media_mobiliado)
## [1] 1486.600 1800.431
Media_diff= media_mobiliado- media_n_mobiliado

Visualizando os dados em Boxplot é possível perceber uma certa significancia pois o terceiro quartil das casas não mobiliadas são práticamente o primeiro quartil das casas mobiliadas, mostrando um aumento nos preços.

ggplot(data=casas_baratas, aes_string(x=casas_baratas$Furniture, y=casas_baratas$valor_aluguel, group=casas_baratas$Furniture)) + geom_boxplot() + theme_classic() + labs(x="Mobiliado ou não", y="Valor do aluguel", title = "Distribuição do aluguel de casas baratas pelo fator de ser mobiliado")

Abaixo segue o código utilizado para o teste de permutação, onde foi considerado a hipótese Ho como verdadeira (Não há diferênça entre as médias). Com base nisso retiramos 1000 vezes de forma aleatoria dois grupos (n1 e n2) e calculamos a média para ver a distribuição da diferença desses valores e comparar com a diferença calculada sem a permutação para encontrar as chances de o acaso ser responsável por essa variação no valor do aluguel.

#tamanho da amostra de mobiliados (819)
length(casas_baratas$valor_aluguel[casas_baratas$Furniture==1])
## [1] 819
#Tamanho da amostra de não mobiliados (4350)
length(casas_baratas$valor_aluguel[casas_baratas$Furniture==0])
## [1] 4350
#Tamanho total de 5169

#Criando funçãod e permutação considerado hipotese h0
perm_fun = function(x,n1,n2){
  n= n1 + n2
  idx_b= sample(1:n, n1)
  idx_a= setdiff(1:n, idx_b)
  mean_diff= mean(x[idx_b], trim = 0.1)- mean(x[idx_a], trim = 0.1)
  return (mean_diff)
}
perm_diffs= rep(0,1000)
for(i in 1:1000)
  perm_diffs[i]= perm_fun(casas_baratas$valor_aluguel, 819,4350)
hist(perm_diffs, xlim=c(-100,500),nclass = 80)
abline(v= (Media_diff))

#valor p
p=mean(perm_diffs>Media_diff)
p
## [1] 0

O histograma acima mostra a frequência da diferença entre as médias encontradas geradas de forma aleatória e a linha constante representa o valor da diferença das médias real que é de R$ 313.815.É possível perceber que é bem improvável que o acaso seja responsável por essa diferença. O p-value associado é 0, logo rejeitamos a hipótese nula de que a diferença pode ser explicada pelo acaso.

Casas Caras Para analisar a diferença do valor do aluguel em casas caras foi utilizado os mesmos procedimentos descritos acima, então para não ficar repetitivo só colocarei o código utilizado.

casas_cara=df[df$tipo=="Caro",]
media_mobiliado=mean(casas_cara$valor_aluguel[casas_cara$Furniture==1], trim = 0.1)
media_n_mobiliado=mean(casas_cara$valor_aluguel[casas_cara$Furniture==0], trim = 0.1)
c(media_n_mobiliado, media_mobiliado)
## [1] 5416.686 5780.968
ggplot(data=casas_cara, aes_string(x=casas_cara$Furniture, y=casas_cara$valor_aluguel, group=casas_cara$Furniture)) + geom_boxplot() + theme_classic() + labs(x="Mobiliado ou não", y="Valor do aluguel", title = "Distribuição do aluguel de casas Caras pelo fator de ser mobiliado")

#tamanho da amostra de mobiliados (1787)
length(casas_cara$valor_aluguel[casas_cara$Furniture==1])
## [1] 1787
#Tamanho da amostra de não mobliados (3736)
length(casas_cara$valor_aluguel[casas_cara$Furniture==0])
## [1] 3736
#Tamanho total de 5523

#Criando funçãod e permutação considerado hipotese h0
perm_fun = function(x,n1,n2){
  n= n1 + n2
  idx_b= sample(1:n, n1)
  idx_a= setdiff(1:n, idx_b)
  mean_diff= mean(x[idx_b], trim = 0.1)- mean(x[idx_a], trim = 0.1)
  return (mean_diff)
}

perm_diffs= rep(0,1000)
for(i in 1:1000)
  perm_diffs[i]= perm_fun(casas_cara$valor_aluguel, 1787,3736)
hist(perm_diffs, xlim=c(-100,500),nclass = 80)
abline(v= (Media_diff))

#valor p
p=mean(perm_diffs>Media_diff)
p
## [1] 0.002

O teste de permutação aleatória gerou um valor de p de 0.001, logo rejeitamos a hipótese nula de que a diferença entre as médias pode ter ocorrido pelo acaso.

Hipotese 2: O valor do aluguel aumenta caso a casa Permita animais ? Para facilitar as análises feitas mudaremos o valor da variável animal, onde caso se aceite animais o valor será 1, caso não aceite o valor será 0, abaixo utilizamos uma função da familia Apply do R para converter os valores categóricos em 0 e 1.

muda = function(x){
  if (x=="acept"){
    return (1)
  }else{
    return(0)
  }
}
df$animal= sapply(df$animal, muda)

Casas Baratas Analisando a média aparada para as casas que aceitam ou não animais se nota uma diferença de R$ 116,241. Agora basta saber se essa diferença pode ser gerada aleatoriamente ou não.

casas_baratas=df[df$tipo=="Barata",]
media_c_animais=mean(casas_baratas$valor_aluguel[casas_baratas$animal==1], trim = 0.1)
media_s_animais=mean(casas_baratas$valor_aluguel[casas_baratas$animal==0], trim = 0.1)
c(media_c_animais, media_s_animais)
## [1] 1566.734 1450.495
#tamanho da amostra de mobiliados (3894)
length(casas_baratas$valor_aluguel[casas_baratas$animal==1])
## [1] 3894
#Tamanho da amostra de não mobiliados (1275)
length(casas_baratas$valor_aluguel[casas_baratas$animal==0])
## [1] 1275
#Tamanho total de 5169

#Criando funçãod e permutação considerado hipotese h0
perm_fun = function(x,n1,n2){
  n= n1 + n2
  idx_b= sample(1:n, n1)
  idx_a= setdiff(1:n, idx_b)
  mean_diff= mean(x[idx_b], trim = 0.1)- mean(x[idx_a], trim = 0.1)
  return (mean_diff)
}

perm_diffs= rep(0,1000)
for(i in 1:1000)
  perm_diffs[i]= perm_fun(casas_baratas$valor_aluguel, 3894,1275)
hist(perm_diffs, xlim=c(-100,500),nclass = 80)
abline(v= (Media_diff))

#valor p
p=mean(perm_diffs>Media_diff)
p
## [1] 0

Casas Caras

casas_caras=df[df$tipo=="Caro",]
media_mobiliado=mean(casas_caras$valor_aluguel[casas_caras$animal==1], trim = 0.1)
media_n_mobiliado=mean(casas_caras$valor_aluguel[casas_caras$animal==0], trim = 0.1)
c(media_n_mobiliado, media_mobiliado)
## [1] 5194.056 5619.596
ggplot(df,aes(y=valor_aluguel))+
      geom_boxplot() + labs(title = "Distribuição dos valores de alugueis")+
    theme_classic()

n=as.data.frame(boxplot.stats(df$valor_aluguel)$out)
dim(n)
## [1] 715   1

Com um p-value de 0 reijeitamos a hipotese nula de que essa diferença possa ter ocorrido devido ao acaso, logo a diferença é estatisticamente significativa.

Criação do modelo de regressão Todos os passos a seguir são referentes a etapa de criação do modelo de regressão, previsão e análise dos pressupostos.

Outlines Para criar o modelo e evitar possíveis problemas com heteroscedasticidades nos resíduos foi decidido retirar os valores extremos da feature valor_aluguel, a qual decidimos fazer as previsões. Abaixo segue um boxplot mostrando a distribuição dos dados e os outlines que essa variável possui.

ggplot(df,aes(y=valor_aluguel))+
      geom_boxplot() + labs(title = "Distribuição dos valores de alugueis")+
    theme_classic()

Foi encontrado um total de 715 valores extremos que correspondem a apoximadamente 6% da coluna valor_aluguel. Devido a isso foi decidido remover essas linhas tendo em vista que ela não iria comprometer o modelo de regressão linear criado, foi utilizado como base para corte o intervalo inter-quartil.

n=as.data.frame(boxplot.stats(df$valor_aluguel)$out)
dim(n)
## [1] 715   1
q1= quantile(df$valor_aluguel,0.25)
q3= quantile(df$valor_aluguel,0.75)
iqr= IQR(df$valor_aluguel)
lim_superior= q1 + 1.5*(iqr)
lim_inferior= q3 - 1.5*(iqr)
df=df[(df$valor_aluguel < lim_superior) & (df$valor_aluguel> lim_inferior),]

ggplot(df,aes(y=valor_aluguel))+
      geom_boxplot() + labs(title = "Distribuição dos valores de alugueis")+
    theme_classic()

Normalizando dados Olhando a base de dados é possível notar variações grandes de escala entre as features do dataset e isso pode acabar afetando o graú de importancia que cada variável tem na regressão linear. Para evitar esse problema decidimos normalizar todas as variáveis do tipo númericas onde os valores irão variar entre -1 e 1.

colnumericas=sapply(df, is.numeric)
df[,colnumericas]=NormData(df[,colnumericas], type = 2)

Separando dados de treino e dados de Teste Antes de fazer a seleção de features e a criação do modelo, criaremos uma função que separará o nosso conjunto de dados em um dataset de treino e outro de teste, 70% dos dados serão utilizados para treinar o modelo e 30% para avaliar a eficácia do mesmo.

df$tipo=NULL
teste_treino <- function(df){
split<- createDataPartition(df$cidade,p=0.7, list=F)
return (split)
}
split= teste_treino(df)
dados_treino=df[split,]
dados_teste=df[-split,]

Feature Select Para garantir um modelo otimizado é necessário escolher as variáveis preditoras que serão utilizadas, essa etapa se faz importante para identificar o número mínimo de features que garantem uma acurácia alta. Pensando nisso utilizaremos será utilizado as seguintes métricas: -> A métrica AIC(Akaike’s Information Criteria), que penaliza a adição de termos a um modelo, onde através de uma regressão passo a passo serão incluidas e excluidas variáveis até que se encontre um modelo que diminua o AIC. -> CP de masllow onde ele compara a precisão e o vício do modelo completo a modelos com um subconjunto de preditores.

Modelo com base no critério de Akaike

completo= lm(valor_aluguel~., data = dados_treino)
modelo1=step(completo, data=dados_treino, direction="backward", trace=F)
summary(modelo1)
## 
## Call:
## lm(formula = valor_aluguel ~ cidade + area + banheiros + n_estacionamento + 
##     andar + animal + Furniture + IPTU + taxa_incendio, data = dados_treino)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.82446 -0.06827 -0.01406  0.04260  1.68147 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -0.297863   0.006851 -43.475  < 2e-16 ***
## cidadeCampinas        0.065067   0.008051   8.082 7.59e-16 ***
## cidadePorto Alegre   -0.129377   0.007398 -17.487  < 2e-16 ***
## cidadeRio de Janeiro  0.068124   0.007283   9.354  < 2e-16 ***
## cidadeSão Paulo       0.119217   0.006068  19.646  < 2e-16 ***
## area                 -0.039373   0.005412  -7.276 3.87e-13 ***
## banheiros             0.013773   0.002962   4.651 3.38e-06 ***
## n_estacionamento     -0.012279   0.002595  -4.732 2.27e-06 ***
## andar1                0.270288   0.007072  38.220  < 2e-16 ***
## andar10               0.294152   0.010731  27.410  < 2e-16 ***
## andar11               0.303509   0.012182  24.914  < 2e-16 ***
## andar12               0.308585   0.012893  23.935  < 2e-16 ***
## andar13               0.335333   0.014697  22.817  < 2e-16 ***
## andar14               0.342009   0.015095  22.657  < 2e-16 ***
## andar15               0.299344   0.017075  17.531  < 2e-16 ***
## andar16               0.300242   0.018860  15.920  < 2e-16 ***
## andar17               0.356547   0.020573  17.331  < 2e-16 ***
## andar18               0.287721   0.024067  11.955  < 2e-16 ***
## andar19               0.327393   0.029933  10.937  < 2e-16 ***
## andar2                0.285513   0.007239  39.439  < 2e-16 ***
## andar20               0.392798   0.038211  10.280  < 2e-16 ***
## andar21               0.349087   0.032855  10.625  < 2e-16 ***
## andar22               0.333828   0.045197   7.386 1.71e-13 ***
## andar23               0.294655   0.043019   6.849 8.13e-12 ***
## andar24               0.316501   0.043096   7.344 2.34e-13 ***
## andar25               0.281019   0.036959   7.604 3.31e-14 ***
## andar26               0.250324   0.041264   6.066 1.39e-09 ***
## andar27               0.424422   0.100553   4.221 2.47e-05 ***
## andar28               0.299446   0.142072   2.108  0.03510 *  
## andar3                0.288712   0.007426  38.881  < 2e-16 ***
## andar301              0.327287   0.142064   2.304  0.02127 *  
## andar32               0.288721   0.142021   2.033  0.04210 *  
## andar4                0.305208   0.007937  38.455  < 2e-16 ***
## andar5                0.302900   0.008720  34.737  < 2e-16 ***
## andar51               0.267712   0.142088   1.884  0.05959 .  
## andar6                0.301428   0.008975  33.584  < 2e-16 ***
## andar7                0.301501   0.009357  32.221  < 2e-16 ***
## andar8                0.309776   0.009357  33.106  < 2e-16 ***
## andar9                0.332922   0.010886  30.582  < 2e-16 ***
## animal               -0.006008   0.001846  -3.255  0.00114 ** 
## Furniture             0.013407   0.001921   6.979 3.29e-12 ***
## IPTU                  0.189890   0.020115   9.440  < 2e-16 ***
## taxa_incendio         0.988811   0.002606 379.378  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1419 on 6193 degrees of freedom
## Multiple R-squared:  0.9801, Adjusted R-squared:   0.98 
## F-statistic:  7280 on 42 and 6193 DF,  p-value: < 2.2e-16

A regressão encontrada se mostrou significativa pelo teste T para todos parâmetros, e conseguiu explicar boa parte da variação do valor do aluguel obtendo um R2 ajustado de 97,76%. #### Critério de cp Masllow

split=teste_treino(df)
dados_treino= df[split,]
dados_teste= df[-split,]
rs=summary(regsubsets(valor_aluguel~., nbest=2,data = dados_treino))
#rs$outmat
n_parametros= as.numeric(rownames(rs$which)) +1
#n_parametros
cp= rs$cp
#s2=rs$rss
R2_ajustado= rs$adjr2 
plot(n_parametros,cp,xlab = "Numero de parametros", ylab = "Estatistica Cp")

#plot(n_parametros,s2,xlab="N parametros",ylab = "S2")
plot(n_parametros,R2_ajustado)

#plot(n_parametros,s2,xlab="N parametros",ylab = "S2")
#n_variaveis=n_parametros-1
#cbind(n_variaveis,n_parametros,cp,R2_ajustado,s2)
modelo2=lm(valor_aluguel~taxa_incendio + cidade +Furniture + area,dados_treino)
summary(modelo2)
## 
## Call:
## lm(formula = valor_aluguel ~ taxa_incendio + cidade + Furniture + 
##     area, data = dados_treino)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.95578 -0.07720  0.00213  0.08058  1.98821 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -0.076005   0.006994 -10.866  < 2e-16 ***
## taxa_incendio         0.961263   0.002613 367.932  < 2e-16 ***
## cidadeCampinas        0.058169   0.010643   5.465 4.80e-08 ***
## cidadePorto Alegre   -0.119361   0.009720 -12.280  < 2e-16 ***
## cidadeRio de Janeiro  0.120855   0.009338  12.943  < 2e-16 ***
## cidadeSão Paulo       0.128655   0.007794  16.506  < 2e-16 ***
## Furniture             0.036487   0.002435  14.982  < 2e-16 ***
## area                 -0.012420   0.002208  -5.625 1.94e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.189 on 6228 degrees of freedom
## Multiple R-squared:  0.9643, Adjusted R-squared:  0.9643 
## F-statistic: 2.404e+04 on 7 and 6228 DF,  p-value: < 2.2e-16

A regressão encontrada se mostrou significativa pelo teste T para todos parâmetros, e conseguiu explicar boa parte da variação do valor do aluguel obtendo um R2 ajustado de 96,45%.

Escolhendo modelo Linear Serão avaliados dois modelos de regressão linear, modelo1 obtido utilizando o método de Akaike e o modelo2 que foi utilizado o críterio de masllow.

Para o modelo 1 foi decidido continuar com todas as 11 variáveis inicialmente colocadas no modelo, já o modelo2 foi decidido ficar apenas com apenas 4 variáveis.

Embora o modelo2 e o modelo1 tenham uma precisão muito próxima, decidimos testar a eficiência dos modelos para dados que eles não tiveram contato. O melhor modelo será aquele que trará melhores valores para as métricas R2 e RMSE. Para diminuir a influência do acaso foi utilizado o método bootstrap, onde serão testados exaustivamente diversos dados de teste e treino para coletar as métricas R2 e RMSE.

#plot(n_parametros,s2,xlab="N parametros",ylab = "S2")
r2_modelo1=c()
RMSE_modelo1= c()
r2_modelo2= c()
RMSE_modelo2= c()


set.seed(0)
for (i in 1:20){
split= teste_treino(df)
dados_treino=df[split,]
dados_teste=df[-split,]
modelo1<- lm(valor_aluguel ~ cidade + quartos + banheiros + n_estacionamento + animal + Furniture + taxa_incendio,data=dados_treino)


#Com o modelo criado será utilizado os dados de testes para avaliar a precisão do modelo criado


previsao<- predict(modelo1,dados_teste)
#Dataframe com valor real e valor previsto
score1= as.data.frame(cbind(dados_teste$valor_aluguel,previsao))
names(score1)=c("Valor_real","Valor_previsto")
head(score1)


#Calculando o R2 para os modelo
r2_modelo1 = rbind(r2_modelo1,R2(score1$Valor_previsto, score1$Valor_real))
RMSE_modelo1= rbind(RMSE_modelo1,RMSE(score1$Valor_previsto, score1$Valor_real))

split= teste_treino(df)
dados_treino=df[split,]
dados_teste=df[-split,]

modelo2 = lm(formula = valor_aluguel ~ taxa_incendio + cidade +Furniture + area, data = dados_treino, 
    na.action = na.omit)
#Com o modelo criado será utilizado os dados de testes para avaliar a precisão do modelo criado


previsao<- predict(modelo2,dados_teste)
#Dataframe com valor real e valor previsto
score2= as.data.frame(cbind(dados_teste$valor_aluguel,previsao))
names(score2)=c("Valor_real","Valor_previsto")
head(score2)


#Calculando o R2 para os modelo
r2_modelo2 = rbind(r2_modelo2,R2(score2$Valor_previsto, score2$Valor_real))
RMSE_modelo2= rbind(RMSE_modelo2,RMSE(score2$Valor_previsto, score2$Valor_real))

}

#####Comparando os dois modelos a partir do R2

c(mean(r2_modelo1),mean(r2_modelo2))
## [1] 0.9634122 0.9621905

Avaliando o R2 temos que o segundo modelo foi um pouco melhor, embora a diferênça encontrada não é muito significativa para decidir qual modelo é melhor.

Comparando os dois modelos a partir do RMSE

c(mean(RMSE_modelo1), mean(RMSE_modelo2))
## [1] 0.1906002 0.1935386

Utilizando o Erro médio quadratico o modelo 1 também teve o maior desemprenho, embora a diferênça também não seja significativa. Como os resultados encontrados foram praticamente iguais decidimos utilizar o critério da parcimônia que visa minimizar escolher o modelo mais simples caso a diferença dos resultados encontrados sejam próximos, logo será escolhido o modelo2 que utilizam apenas 4 variáveis para explicar o valor do aluguel.

Avaliando as premissas de regressão linear Linearidade dos coeficientes e variáveis. Esta primeira hipotese é verdadeira, o software R garante que o modelo criado tenha coeficientes lineares.

par(par(mfrow=c(1,2)))
plot(valor_aluguel~taxa_incendio, main="Valor Aluguel Vs Taxa de Incendio",df)

plot(valor_aluguel~area, main="Valor Aluguel vs área",df)

Normalidade dos resíduos Uma das principais premissas para um modelo de regressão linear é que os residuos sigam uma distribuição normal, inicialmente plotamos o histograma para avaliar gráficamente se os dados são normais. É possível notar que embora ele tenha toda a forma de sino caracteristica desse tipo de distribuição, as caudas são longas decorrentes valores considerados extremos (Embora tenhamos retirado valores extremos da coluna valo_aluguel, as outras variáveis também possuiam valores outlines)

par(mfrow=c(1,2))
hist(modelo2$residuals, main="Histograma dos Residuos", breaks = 50)
qqnorm(modelo2$residuals, main = "QQ-Norm")
qqline(modelo2$residuals)

O gráfico qq normal deixa claro que os resíduos do modelo não seguem uma distribuição normal, mas para se afimar com certeza foi feito o teste de shapiro wilk para as primeiras 5000 observaçõs (limite máximo que o teste suporta).

shapiro.test(modelo1$residuals[1:5000])
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo1$residuals[1:5000]
## W = 0.87234, p-value < 2.2e-16

Utilizando um nivel de confiança de 95% o teste de shapiro gerou um p-value de praticamente 0, indicando que podemos rejeitar a hipótese nula de que os dados seguem uma distribuição normal.

Multicolinearidade das variáveis independentes O primeiro indício da baixa multicolinearidade entre as variáveis regressoras se pode observar analisando os coeficientes do modelo que deram todos significativos para o teste t. Porém para se ter mais certeza será utilizado o FIV (Fator de inflação de variância) para avaliar a relação entre as variáveis explicativas

vif(modelo2)
##        taxa_incendio       cidadeCampinas   cidadePorto Alegre 
##             1.790567             1.614867             1.847294 
## cidadeRio de Janeiro      cidadeSão Paulo            Furniture 
##             1.984109             2.692913             1.084470 
##                 area 
##             1.690016

As variáveis explicativas estão com um fator de inflação baixos, como o valor do FIV para as features deram no máximo 2.69, pela classificação do método podemos afirmar que existe uma baixa multicolinearidade entre as features.

Homoscedásticidade dos resíduos Analisando o plot dos residuos do modelo não é possível avaliar bem o gráfico devido a quantidade alta de pontos, mesmo assim da para notar um desajuste nos resíduos, parece existir mais pontos subestimados do que superstimados gerando mais pontos extremos positivos que negativos e isso pode ser um indício de heteroscedasticidade.

plot(modelo2$residuals, main="Resíduos x Valores Ajustados", xlab="Valores",ylab="Resíduos")

Pelo teste de Breusch Pagan, onde temos a hipotese Ho indicando que os resíduos são homoscedásticos e H1 indicando que são heterocedásticos.

bptest(modelo2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 167.14, df = 7, p-value < 2.2e-16

Ausência de autocorrelação dos residuos Para uma análise mais formal foi utilizado o teste estátistico de Durbin Watson onde a Hipótese nula Ho indica a ausência de autocorrelação e a hipótese alternativa H1 a presença de autocorrelação. Para um nivel de significancia de 95% foi encontrado um p-value acima de 0.05, logo não podemos rejeitar a hipótese nula de que os resíduos não possúem autocorrelação.

dwtest(modelo2)
## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 1.9766, p-value = 0.178
## alternative hypothesis: true autocorrelation is greater than 0

Análise dos resultados encontrados pela previsão feita Como foi dito anteriormente, foi utilizado uma parte dos dados para criar o modelo e outra para testar a eficiência da regressão quando ela tenta prever valores novos. Abaixo segue uma pequena análise dos resíduos gerados quando o modelo tenta prever valores de alugueis que ele não teve contato.

residuo<- function(real,previsto){
  return(previsto-real)
}
residu=as.data.frame(mapply( residuo,score2$Valor_real,score2$Valor_previsto))
names(residu)="residuo"

par(mfrow=c(1,2))
hist(residu$residuo, main="Histograma dos Residuos", breaks = 50)
qqnorm(residu$residuo, main = "QQ-Norm")
qqline(residu$residuo)

O histograma mostra que o modelo tem a maioria das suas previsões próxima a zero o que indica uma alta precisão gerando um baixo erro, porém ele tende a superestimar os valores dos alugueis em alguns casos gerando um leve aumento aparente no lado direito do histograma. O gráfico qq normal mostra que tem uma grande chance de os resíduos não serem normais semelhante a análise feita dos erros gerados no modelo mostrados durante os pressupostos. Acredito que esse problema se dá devido a falta de uma feature muito importante que é a localização da casa nas cidades causando heterocedasticidade e a não normalidade.