Trabalho Regressão

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\()= Valor do aluguel; property tax (R\))= IPTU da casa; fire insurance (R$)= Seguro de incendio da casa; Total= Total pago ao mês pela casa;

Biblioteca utilizada

#install.packages("readr")
#install.packages("plotly")
library(readr)
library(visdat)
library(dplyr)
library(ggplot2)
library(plotly)
library(dplyr)
library(corrplot)
library(randomForest)
library(caret)
#install.packages("faraway")
library(faraway)
#install.packages("lmtest")
library(lmtest)
#install.packages("MVar.pt")
library(MVar.pt)
#install.packages(nortest)
library(nortest)
#install.packages("fpp")
library(fpp)
#install.packages("leaps")
library(leaps)

Amostra dos dados

df<- read_csv("houses_to_rent_v2.csv")

## 
## -- Column specification -------------------
## cols(
##   city = col_character(),
##   area = col_double(),
##   rooms = col_double(),
##   bathroom = col_double(),
##   `parking spaces` = col_double(),
##   floor = col_character(),
##   animal = col_character(),
##   furniture = col_character(),
##   `hoa (R$)` = col_double(),
##   `rent amount (R$)` = col_double(),
##   `property tax (R$)` = col_double(),
##   `fire insurance (R$)` = col_double(),
##   `total (R$)` = col_double()
## )

head(df)

## # A tibble: 6 x 13
##   city   area rooms bathroom `parking spaces` floor animal furniture `hoa (R$)`
##   <chr> <dbl> <dbl>    <dbl>            <dbl> <chr> <chr>  <chr>          <dbl>
## 1 São ~    70     2        1                1 7     acept  furnished       2065
## 2 São ~   320     4        4                0 20    acept  not furn~       1200
## 3 Port~    80     1        1                1 6     acept  not furn~       1000
## 4 Port~    51     2        1                0 2     acept  not furn~        270
## 5 São ~    25     1        1                0 1     not a~ not furn~          0
## 6 São ~   376     3        3                7 -     acept  not furn~          0
## # ... with 4 more variables: `rent amount (R$)` <dbl>, `property tax
## #   (R$)` <dbl>, `fire insurance (R$)` <dbl>, `total (R$)` <dbl>

Metadados

glimpse(df)

## Rows: 10,692
## Columns: 13
## $ city                  <chr> "São Paulo", "São Paulo", "Porto Alegre", "Po...
## $ area                  <dbl> 70, 320, 80, 51, 25, 376, 72, 213, 152, 35, 2...
## $ rooms                 <dbl> 2, 4, 1, 2, 1, 3, 2, 4, 2, 1, 1, 1, 1, 1, 2, ...
## $ bathroom              <dbl> 1, 4, 1, 1, 1, 3, 1, 4, 2, 1, 1, 1, 1, 1, 2, ...
## $ `parking spaces`      <dbl> 1, 0, 1, 0, 0, 7, 0, 4, 1, 0, 0, 1, 0, 1, 2, ...
## $ floor                 <chr> "7", "20", "6", "2", "1", "-", "7", "4", "3",...
## $ animal                <chr> "acept", "acept", "acept", "acept", "not acep...
## $ furniture             <chr> "furnished", "not furnished", "not furnished"...
## $ `hoa (R$)`            <dbl> 2065, 1200, 1000, 270, 0, 0, 740, 2254, 1000,...
## $ `rent amount (R$)`    <dbl> 3300, 4960, 2800, 1112, 800, 8000, 1900, 3223...
## $ `property tax (R$)`   <dbl> 211, 1750, 0, 22, 25, 834, 85, 1735, 250, 35,...
## $ `fire insurance (R$)` <dbl> 42, 63, 41, 17, 11, 121, 25, 41, 191, 30, 27,...
## $ `total (R$)`          <dbl> 5618, 7973, 3841, 1421, 836, 8955, 2750, 7253...

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 x 13
##   cidade  area quartos banheiros n_estacionamento andar animal Furniture
##   <chr>  <dbl>   <dbl>     <dbl>            <dbl> <chr> <chr>  <chr>    
## 1 São P~    70       2         1                1 7     acept  furnished
## 2 São P~   320       4         4                0 20    acept  not furn~
## 3 Porto~    80       1         1                1 6     acept  not furn~
## 4 Porto~    51       2         1                0 2     acept  not furn~
## 5 São P~    25       1         1                0 1     not a~ not furn~
## 6 São P~   376       3         3                7 -     acept  not furn~
## # ... with 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

Para começarmos a fazer qualquer tipo de análise, precisamos antes de tudo validar se as variáveis do nosso dataset estão formatadas com tipo de dados certo. Para ficar visualmente mais atrativo, decidimos utilizar o pacote visdat.

vis_dat(df)

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 x 13
##    cidade  area quartos banheiros n_estacionamento andar animal Furniture
##    <fct>  <dbl>   <dbl>     <int>            <dbl> <chr> <fct>  <fct>    
##  1 São P~    70       2         1                1 7     acept  furnished
##  2 São P~   320       4         4                0 20    acept  not furn~
##  3 Porto~    80       1         1                1 6     acept  not furn~
##  4 Porto~    51       2         1                0 2     acept  not furn~
##  5 São P~    25       1         1                0 1     not a~ not furn~
##  6 São P~   376       3         3                7 -     acept  not furn~
##  7 Rio d~    72       2         1                0 7     acept  not furn~
##  8 São P~   213       4         4                4 4     acept  not furn~
##  9 São P~   152       2         2                1 3     acept  furnished
## 10 Rio d~    35       1         1                0 2     acept  furnished
## # ... with 10,682 more rows, and 5 more variables: condominio <dbl>,
## #   valor_aluguel <dbl>, IPTU <dbl>, taxa_incendio <dbl>, total_pago <dbl>

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)

## $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

## 
## $n_estacionamento

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)") 

Fica claro que a relação do valor do 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

Media_diff= media_mobiliado- media_n_mobiliado

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.001

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

Media_diff= media_c_animais- media_s_animais

Plotando o bloxplot para analisar a distribuição dos dados em torno da sua mediana não foi possível encontrar um deslocamento discrepante que possa justificar uma significancia da diferença encontrada entre os valores.

ggplot(data=casas_baratas, aes_string(x=casas_baratas$animal, y=casas_baratas$valor_aluguel, group=casas_baratas$animal)) + 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")

Avaliando agora o teste de permutação aleatória é possível notar que a diferença entre os valores do aluguel é significativa, mostrando que o acaso tem uma grande chance de não ter causado essa diferença.

#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

Media_diff= media_mobiliado- media_n_mobiliado

#tamanho da amostra de mobiliados (3894)
length(casas_caras$valor_aluguel[casas_caras$animal==1])

## [1] 4422

#Tamanho da amostra de não mobiliados (1275)
length(casas_caras$valor_aluguel[casas_caras$animal==0])

## [1] 1101

#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

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 importância 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 para que variem 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 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 + quartos + banheiros + 
##     n_estacionamento + andar + animal + Furniture + taxa_incendio, 
##     data = dados_treino)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76445 -0.07207 -0.01719  0.04295  2.02784 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -0.316336   0.007022 -45.049  < 2e-16 ***
## cidadeCampinas        0.068809   0.008575   8.025 1.21e-15 ***
## cidadePorto Alegre   -0.126145   0.007899 -15.970  < 2e-16 ***
## cidadeRio de Janeiro  0.073454   0.007738   9.493  < 2e-16 ***
## cidadeSão Paulo       0.132324   0.006469  20.454  < 2e-16 ***
## area                 -0.004337   0.001840  -2.357 0.018463 *  
## quartos              -0.004369   0.002907  -1.503 0.132858    
## banheiros             0.020477   0.003065   6.680 2.60e-11 ***
## n_estacionamento     -0.012156   0.002648  -4.590 4.52e-06 ***
## andar1                0.283003   0.007277  38.888  < 2e-16 ***
## andar10               0.324118   0.011418  28.387  < 2e-16 ***
## andar11               0.325274   0.012405  26.222  < 2e-16 ***
## andar12               0.319562   0.013890  23.006  < 2e-16 ***
## andar13               0.348669   0.014288  24.403  < 2e-16 ***
## andar14               0.333804   0.015696  21.267  < 2e-16 ***
## andar15               0.317827   0.018744  16.957  < 2e-16 ***
## andar16               0.317426   0.020756  15.293  < 2e-16 ***
## andar17               0.335537   0.021724  15.446  < 2e-16 ***
## andar18               0.312278   0.023912  13.059  < 2e-16 ***
## andar19               0.400177   0.028773  13.908  < 2e-16 ***
## andar2                0.299852   0.007676  39.064  < 2e-16 ***
## andar20               0.458636   0.035747  12.830  < 2e-16 ***
## andar21               0.360596   0.035749  10.087  < 2e-16 ***
## andar22               0.355671   0.053370   6.664 2.89e-11 ***
## andar23               0.281288   0.043659   6.443 1.26e-10 ***
## andar24               0.292725   0.047790   6.125 9.61e-10 ***
## andar25               0.303784   0.037925   8.010 1.36e-15 ***
## andar26               0.258763   0.040533   6.384 1.85e-10 ***
## andar27               0.482343   0.087034   5.542 3.11e-08 ***
## andar28               0.378274   0.106492   3.552 0.000385 ***
## andar3                0.306092   0.007746  39.515  < 2e-16 ***
## andar301              0.352029   0.150513   2.339 0.019375 *  
## andar4                0.317422   0.008296  38.261  < 2e-16 ***
## andar5                0.310275   0.009171  33.831  < 2e-16 ***
## andar51               0.283248   0.150549   1.881 0.059960 .  
## andar6                0.308680   0.009352  33.006  < 2e-16 ***
## andar7                0.311947   0.009933  31.404  < 2e-16 ***
## andar8                0.317810   0.010040  31.654  < 2e-16 ***
## andar9                0.340501   0.011399  29.871  < 2e-16 ***
## animal               -0.008246   0.001974  -4.176 3.00e-05 ***
## Furniture             0.015442   0.002057   7.508 6.84e-14 ***
## taxa_incendio         0.983046   0.002743 358.364  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1504 on 6194 degrees of freedom
## Multiple R-squared:  0.9773, Adjusted R-squared:  0.9771 
## F-statistic:  6491 on 41 and 6194 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)

#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.91025 -0.07639  0.00283  0.07781  1.98817 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -0.083287   0.006944 -11.994  < 2e-16 ***
## taxa_incendio         0.962790   0.002607 369.376  < 2e-16 ***
## cidadeCampinas        0.064008   0.010577   6.052 1.52e-09 ***
## cidadePorto Alegre   -0.104237   0.009651 -10.800  < 2e-16 ***
## cidadeRio de Janeiro  0.128755   0.009265  13.897  < 2e-16 ***
## cidadeSão Paulo       0.136966   0.007749  17.676  < 2e-16 ***
## Furniture             0.037231   0.002414  15.424  < 2e-16 ***
## area                 -0.010789   0.002230  -4.839 1.34e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1878 on 6228 degrees of freedom
## Multiple R-squared:  0.9651, Adjusted R-squared:  0.9651 
## F-statistic: 2.462e+04 on 7 and 6228 DF,  p-value: < 2.2e-16

Para escolher o número de variáveis do modelo foi utilizado o gráfico do R2_ajustado onde com 4 variáveis o R2 não tem mais ganhos aparentes. 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.

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 primeiro 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 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 resíduos 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

O teste gerou um p-value menor que 0.05, logo para um nível de confiança de 95% podemos rejeitar a hipótese nula de que os resíduos sao homoscedásticos.

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. Acreditamos 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.

Plotando valores reais Vs Previstos

Em seguida foi plotado os valores reais e previstos em um gráfico para visualizar a diferença entre os valores. É possível notar que os valores ficaram bem proximos um dos outros, mostrando que a análise foi bastante eficaz. Foi mostrado apenas os 200 primeiros valores previstos por motivos de melhor visualização.

score2$index=c(1:2670)
ggplotly(ggplot(data=score2, aes(x=index,y=Valor_real, col="Real"))+ geom_line() + labs(title = "Valores Previstos e reais", x="Observações", y="Valores") +geom_line(data=score2, aes(x=index, y=Valor_previsto, col="Previsto")) +theme_classic() + xlim(c(0,200)) + scale_color_manual(values = c("black","red")) + guides(color=guide_legend(title=""))
)