Modelagem HousePrices
Arthur Machado
04/10/2019
Objetivos
Neste relatório tem como principal objetivo realizar um modelo de rede neural e comparar com um modelo de regressão linear, tal que consiga melhor predizer sob a natureza dos dados. O critério de seleção é o erro quadrático médio, no qual o modelo que tiver menor valor dessa medida melhor.
Banco de dados
É uma base do site KAGGLE, a respeito de valores imobiliários de residências, neste relatório foi utilizado apenas variáveis de caráter numérico, sendo deixado pra próximas aplicações as variáveis qualitativas.
library(keras)
library(tidyverse)
library(MASS)
library(gridExtra)
train <- read_csv("/opt/datasets/houseprices/train.csv", na = "")set.seed(123)
train_num <- train %>%
select_if(is.numeric) %>%
mutate_at(vars(ends_with("Area")), function(x){log(x+1)}) %>%
mutate_at(vars(ends_with("Sold")), function(x){log(x+1)}) %>%
mutate_at(vars(ends_with("Price")), function(x){log(x+1)}) %>%
mutate_at(-1, scale) %>%
mutate(Flag = sample(c("Test", "Train"),
n(),
replace = TRUE,
prob = c(.2,.8)))Modelo de regressão linear
Nesta seção será aplicada uma modelagem mais difundida da literatura que é do regressão linear, no qual nessa metodologia, precisa de uma hipótese distribucional grande, que é que os dados são normalmente distribuídos. Como estamos manipulando a base de dados e aplicando o log na variável resposta, então essa opção pode ser bem viável.
Definindo as bases de treinamento e de teste do modelo.
#Base de treinamento
train_all <- train_num %>%
filter(Flag=="Train") %>%
dplyr::select(-Id, -Flag, -TotalBsmtSF) %>%
as.data.frame()
test<- train_num %>%
filter(Flag=="Test") %>%
dplyr::select(-Id, -Flag, - TotalBsmtSF) %>%
as.data.frame()Definindo o modelo de regressão pela função lm().
m1<- lm(SalePrice ~. , data=train_all )
summary(m1)##
## Call:
## lm(formula = SalePrice ~ ., data = train_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6330 -0.1415 0.0107 0.1700 1.0743
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0176800 0.0092536 1.911 0.056303 .
## MSSubClass -0.0023876 0.0126311 -0.189 0.850106
## LotArea 0.1110315 0.0127796 8.688 < 2e-16 ***
## OverallQual 0.2497075 0.0168339 14.834 < 2e-16 ***
## OverallCond 0.1484838 0.0120697 12.302 < 2e-16 ***
## YearBuilt 0.2403111 0.0185719 12.939 < 2e-16 ***
## YearRemodAdd 0.0494122 0.0138891 3.558 0.000389 ***
## BsmtFinSF1 0.2084667 0.0228521 9.122 < 2e-16 ***
## BsmtFinSF2 0.0550703 0.0118088 4.664 3.47e-06 ***
## BsmtUnfSF 0.1218484 0.0197204 6.179 8.96e-10 ***
## `1stFlrSF` 0.0232742 0.0367197 0.634 0.526316
## `2ndFlrSF` 0.0337338 0.0365463 0.923 0.356179
## LowQualFinSF -0.0054130 0.0103170 -0.525 0.599915
## GrLivArea 0.2881089 0.0450141 6.400 2.25e-10 ***
## BsmtFullBath 0.0347822 0.0140477 2.476 0.013430 *
## BsmtHalfBath -0.0094616 0.0095180 -0.994 0.320395
## FullBath 0.0145609 0.0158446 0.919 0.358299
## HalfBath 0.0060257 0.0135785 0.444 0.657295
## BedroomAbvGr -0.0300429 0.0145373 -2.067 0.038995 *
## KitchenAbvGr -0.0505453 0.0116770 -4.329 1.63e-05 ***
## TotRmsAbvGrd 0.0457145 0.0204405 2.236 0.025512 *
## Fireplaces 0.0535861 0.0117203 4.572 5.35e-06 ***
## GarageCars 0.0944362 0.0177755 5.313 1.30e-07 ***
## GarageArea 0.0001765 0.0148534 0.012 0.990523
## WoodDeckSF 0.0196851 0.0102996 1.911 0.056221 .
## OpenPorchSF 0.0059326 0.0101156 0.586 0.557667
## EnclosedPorch 0.0213352 0.0104984 2.032 0.042360 *
## `3SsnPorch` 0.0092104 0.0104338 0.883 0.377561
## ScreenPorch 0.0321203 0.0099129 3.240 0.001229 **
## PoolArea -0.0010062 0.0103547 -0.097 0.922608
## MiscVal -0.0299971 0.0148930 -2.014 0.044223 *
## MoSold 0.0030615 0.0094576 0.324 0.746217
## YrSold -0.0129617 0.0094918 -1.366 0.172341
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3172 on 1147 degrees of freedom
## Multiple R-squared: 0.8998, Adjusted R-squared: 0.897
## F-statistic: 321.8 on 32 and 1147 DF, p-value: < 2.2e-16Podemos perceber que muitas dessas 33 variáveis escolhidas para o modelo, muitas não são significantes para o modelo, com p-valores ácima de um nível de significância \(\alpha\). Com isso, podemos aplicar um critério de seleção por AIC, no qual as variáveis que tiverem maior AIC vão ser descartadas. Então, podemos utilizar a função stepAIC().
m2<- stepAIC(m1)summary(m2)##
## Call:
## lm(formula = SalePrice ~ LotArea + OverallQual + OverallCond +
## YearBuilt + YearRemodAdd + BsmtFinSF1 + BsmtFinSF2 + BsmtUnfSF +
## GrLivArea + BsmtFullBath + BedroomAbvGr + KitchenAbvGr +
## TotRmsAbvGrd + Fireplaces + GarageCars + WoodDeckSF + EnclosedPorch +
## ScreenPorch + MiscVal, data = train_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4910 -0.1466 0.0122 0.1759 1.0732
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.017160 0.009220 1.861 0.062976 .
## LotArea 0.111533 0.011079 10.068 < 2e-16 ***
## OverallQual 0.251755 0.016569 15.194 < 2e-16 ***
## OverallCond 0.147930 0.011857 12.476 < 2e-16 ***
## YearBuilt 0.246918 0.016774 14.720 < 2e-16 ***
## YearRemodAdd 0.050731 0.013558 3.742 0.000192 ***
## BsmtFinSF1 0.201303 0.016892 11.917 < 2e-16 ***
## BsmtFinSF2 0.051712 0.010482 4.934 9.25e-07 ***
## BsmtUnfSF 0.117191 0.013991 8.376 < 2e-16 ***
## GrLivArea 0.332502 0.021113 15.749 < 2e-16 ***
## BsmtFullBath 0.035550 0.012910 2.754 0.005984 **
## BedroomAbvGr -0.030944 0.014059 -2.201 0.027936 *
## KitchenAbvGr -0.052051 0.010409 -5.001 6.60e-07 ***
## TotRmsAbvGrd 0.051440 0.019705 2.611 0.009155 **
## Fireplaces 0.053948 0.011513 4.686 3.12e-06 ***
## GarageCars 0.095102 0.013023 7.303 5.23e-13 ***
## WoodDeckSF 0.017799 0.010155 1.753 0.079904 .
## EnclosedPorch 0.020676 0.010321 2.003 0.045380 *
## ScreenPorch 0.030839 0.009813 3.143 0.001716 **
## MiscVal -0.028776 0.014762 -1.949 0.051503 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3164 on 1160 degrees of freedom
## Multiple R-squared: 0.8991, Adjusted R-squared: 0.8975
## F-statistic: 544.2 on 19 and 1160 DF, p-value: < 2.2e-16Após isso, podemos fazer uma análise de resíduo para esse modelo normal, afim de estudar a evidência de pontos discrepantes e de certas suposições do modelo.
X <- model.matrix(m2)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
h <- diag(H)
lms <- summary(m2)
s <- lms$sigma
res<-m2 %>% residuals()
res_padronize<- res/(s*sqrt(1-h))
pred<- m2 %>% predict()Com os códigos abaixo, podemos ver no histograma que os resíduos apresentam-se centrados na média, porém é visto que existe a presença de alguns outliers. No boxplot a evidência de outliers é notável, com diversas observações fora das delimitações do boxplot.
p1<- ggplot() + geom_histogram(aes(x = res))
p2<- ggplot() + geom_boxplot(aes(x = res, y=res))
gridExtra::grid.arrange(p1,p2,nrow=1,ncol=2)
Seguindo o protocolo, podemos agora fazer dois gráficos, um é o \(\mathrm{resíduo} \times \mathrm{Ordem}\) e outro é o de \(\mathrm{resíduo} \times \mathrm{Ajuste}\), o primeira é mais pra estudar se os resíduos estão distribuídos aleatoriamente, sem nenhum padrão, e o segundo é a respeitos de pontos discrepantes no modelo, tal que, se algumas das obsevarções ultrapassarm os limites de -2 e 2, então essa obsevação pode ser tratada com um ponto discrepantes para o modelo, sendo a presença dessa variável prejudical para o modelo.
Podemos perceber, que os resíduos estão bem distribuídos, porém com algumaws observações fora das regiões de 2 e -2.
order<-1:1180
p3<- ggplot() + geom_point(aes(y = res, x=order)) +
geom_hline(yintercept = 0, color = "red") +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red")
p4<- ggplot() + geom_point(aes(y = res, x = pred)) +
geom_hline(yintercept = 0, color = "red") +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red")
grid.arrange(p3,p4,nrow=2,ncol=1)
Utilizando esta função pra a obtenção do gráfico de envelope, podemos observar que alguns dados apresentam dentro das regiões de confiança, porém apresentando uma certa não-linearidade, com alguns pontos discrepantes. Esta função foi obtida do site do professor Gilberto A.PAULA da Universidade Estadual de São Paulo (USP).
envel_norm<- function(fit.model){
par(mfrow=c(1,1))
X <- model.matrix(fit.model)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
h <- diag(H)
si <- lm.influence(fit.model)$sigma
r <- resid(fit.model)
tsi <- r/(si*sqrt(1-h))
#
ident <- diag(n)
epsilon <- matrix(0,n,100)
e <- matrix(0,n,100)
e1 <- numeric(n)
e2 <- numeric(n)
#
for(i in 1:100){
epsilon[,i] <- rnorm(n,0,1)
e[,i] <- (ident - H)%*%epsilon[,i]
u <- diag(ident - H)
e[,i] <- e[,i]/sqrt(u)
e[,i] <- sort(e[,i]) }
#
for(i in 1:n){
eo <- sort(e[i,])
e1[i] <- (eo[2]+eo[3])/2
e2[i] <- (eo[97]+eo[98])/2 }
#
med <- apply(e,1,mean)
faixa <- range(tsi,e1,e2)
#
par(pty="s")
qqnorm(tsi,xlab="Percentil da N(0,1)",
ylab="Residuo Studentizado", ylim=faixa, pch=16, main="")
par(new=TRUE)
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(med,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=2, main="")
}
envel_norm(m2)
E por final, foi obtido um erro quadrático medio de 0.316434.
mean_absolute_error<-lms$sigma
mean_absolute_error## [1] 0.316434Modelo de rede neural
Após a aplicação do modelo de regressão linear pela seção anterior, agora podemos fazer um modelo de regressão por rede neural, no qual a última camada de vez de ser o número de classficadores é apenas o valor predito do preço da residência.
#Variáveis preditoras na base de treino
x_train <- train_num %>%
filter(Flag == "Train") %>%
dplyr::select(-Id, -Flag, -SalePrice) %>%
as.matrix()
#Variável resposta na base de treino
y_train <- train_num %>%
filter(Flag == "Train") %>%
dplyr::select(SalePrice) %>%
as.matrix()
#Variáveis preditoras na base de teste
x_test <- train_num %>%
filter(Flag == "Test") %>%
dplyr::select(-Id, -Flag, -SalePrice) %>%
as.matrix()
#Variável resposta na base de teste
y_test<- train_num %>%
filter(Flag == "Test") %>%
dplyr::select(SalePrice) %>%
as.matrix()Definindo o modelo de rede naural abaixo, notamos que ele apresenta uma estrutura bem mais simples em relação aos modelos de rede neural para detecção de imagem. Vale salientar que nesses modelos encobrem casos que não necessariamente precisa de inúmeras camadas de neurônios, levando a implementação de uma forma mais simples.
O modelo introduzido é constituído de 3 camadas, a primeira camada tem 32 neurônios, a segunda o número de neurônios caí pra metade e por final, a camada obrigatoriamente precisa de apenas 1 neurônio, que no caso é justamente a previsão de preco de venda da casa. O monitoramente que foi feito para o modelo é do erro quadrático médio e não da taxa de acerto introduzida nos modelos de classificação anteriormente.
model <- keras_model_sequential() %>%
layer_dense(units = 32, activation = "relu",
input_shape = ncol(train_num) - 3, kernel_initializer='normal') %>%
layer_dense(units =16, activation = "relu", kernel_initializer='normal') %>%
layer_dense(units = 1)
model %>% compile(
optimizer = optimizer_adam(),
loss = 'mse',
metrics = list("mean_absolute_error")
)Como as operações são mais “básicas” computacionalmente falando, então foi utilizado no código abaixo 40 iterações no modelo, uma condição de parada no algoritmo, tal que, após as realizações dos resultados, se com isso for atingido o esperado, o programa vai parar com até 20 iterações.
history <- model %>% fit(
x_train, y_train,
epochs = 40,
validation_split = 0.2,
verbose = 1,
shuffle = TRUE,
callbacks = list(
callback_early_stopping(monitor = "val_mean_absolute_error", min_delta = 0.01,
patience = 20,
verbose = 0,
mode = "auto",
restore_best_weights = TRUE)
))
plot(history)
Por final, foi obtido um valor menor que o modelo de regressão linear usual. Porém, com uma complexidade bem maior de entendimento em relação ao modelos lineares comuns.
evaluate(model, x_test, y_test)## $loss
## [1] 0.2034168
##
## $mean_absolute_error
## [1] 0.2644975