DEFINIÇÃO DO PROBLEMA DE NEGÓCIO:

Nosso objetivo é construir um modelo de Machine Learning que seja capaz de fazer previsões sobre a taxa média de ocupação de casas na região de Boston, EUA, por proprietários. A variável a ser prevista é um valor numérico que representa a mediana da taxa de ocupação das casas em Boston. Para cada casa temos diversas variáveis explanatórias. Sendo assim, podemos resolver este problema empregando Regressão Linear Simples ou Múltipla.

DEFININDO O DATASET:

Usaremos o Boston Housing Dataset, que é um conjunto de dados que tem a taxa média de ocupação das casas, juntamente com outras 13 variáveis que podem estar relacionadas aos preços das casas. Esses são os fatores como condições socioeconômicas, condições ambientais, instalações educacionais e alguns outros fatores semelhantes.

DICIONÁRIO DE DADOS (Traduzido para pt-br):

CRIM: Taxa de crime per capita por cidade.

ZN: Proporção de terrenos residenciais divididos em lotes com mais de 25.000 pés quadrados.

INDUS: Proporção de acres comerciais não comerciais por cidade.

CHAS: Variável fictícia Charles River (= 1 se o trecho limita o rio; 0 caso contrário).

NOX: Concentração de óxidos nítricos (partes por 10 milhões).

RM: Número médio de quartos por habitação.

AGE: Proporção de unidades ocupadas pelo proprietário construídas antes de 1940.

DIS: Distâncias ponderadas para cinco centros de emprego em Boston.

RAD: Índice de acessibilidade às rodovias radiais.

TAX: Taxa de imposto sobre a propriedade de valor total por US $ 10.000.

PTRATIO: Proporção aluno-professor por cidade.

BLACK: 1000 (Bk - 0,63) ^ 2 em que Bk é a proporção de negros por cidade.

LSTAT: % menor status da população.

MEDV: Valor médio das casas ocupadas pelos proprietários em US $ 1000.

LIBRARYS:
pacman::p_load(tidyverse,caTools,knitr,corrplot,car,olsrr,moments,MASS,caret,klaR)
CARREGANDO O DATASET:
df <- Boston
COLETANDO VARIÁVEIS:
glimpse(df)
## Rows: 506
## Columns: 14
## $ crim    <dbl> 0.00632, 0.02731, 0.02729, 0.03237, 0.06905, 0.02985, 0.08829…
## $ zn      <dbl> 18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.5, 12.5, 12.5, 12.5, 12.5, …
## $ indus   <dbl> 2.31, 7.07, 7.07, 2.18, 2.18, 2.18, 7.87, 7.87, 7.87, 7.87, 7…
## $ chas    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ nox     <dbl> 0.538, 0.469, 0.469, 0.458, 0.458, 0.458, 0.524, 0.524, 0.524…
## $ rm      <dbl> 6.575, 6.421, 7.185, 6.998, 7.147, 6.430, 6.012, 6.172, 5.631…
## $ age     <dbl> 65.2, 78.9, 61.1, 45.8, 54.2, 58.7, 66.6, 96.1, 100.0, 85.9, …
## $ dis     <dbl> 4.0900, 4.9671, 4.9671, 6.0622, 6.0622, 6.0622, 5.5605, 5.950…
## $ rad     <int> 1, 2, 2, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4…
## $ tax     <dbl> 296, 242, 242, 222, 222, 222, 311, 311, 311, 311, 311, 311, 3…
## $ ptratio <dbl> 15.3, 17.8, 17.8, 18.7, 18.7, 18.7, 15.2, 15.2, 15.2, 15.2, 1…
## $ black   <dbl> 396.90, 396.90, 392.83, 394.63, 396.90, 394.12, 395.60, 396.9…
## $ lstat   <dbl> 4.98, 9.14, 4.03, 2.94, 5.33, 5.21, 12.43, 19.15, 29.93, 17.1…
## $ medv    <dbl> 24.0, 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 27.1, 16.5, 18.9, 1…
dim(df)
## [1] 506  14

É possível notar que, as variáveis numéricas possuem escalas diferentes. Precisaremos aplicar uma técnica de padronização nessas variáveis.

O dataset também possui 506 observacoes e 14 variáveis.

PADRONIZANDO VARIÁVEIS NUMÉRICAS:
dados_padronizados <- scale(df[,1:13])

Seleciona todas as variaveis, exceto a a variavel target (MEDV).

y <- df[,14]

Seleciona somente a variável target (MEDV).

df_p <- cbind(dados_padronizados,y);df_p <- as.data.frame(df_p)

Unifica o dataset com padronização.

MATRIZ DE CORRELAÇÃO:
num_vars <- sapply(df_p,is.numeric) 
corr <- cor(df_p[num_vars]) 
corrplot(corr, method = "color")

Avaliando a matrix de correlacao, é possivel ver que há grande correlacao entre as variáveis:

Ou seja, quanto mais comercial for a área, maior a sua concentracao de oxido nitrico. O que faz todo sentido, pois áreas comerciais costumam, ser mais poluentes em comparação à areas não comercias.

p1 <- ggplot(df, aes(indus,nox))+
  geom_point()+
  labs(x = "Proporcao de acres comerciais", y = "Concentracao de oxido nitrico",title = "Analise de correlação - indus, nox");p1

Ou seja, acres comerciais, tendem a pagar mais impostos.

p2 <- ggplot(df, aes(indus,tax))+
  geom_point()+
  labs(x = "Proporcao de acres comerciais",y = "Taxa de impostos sobre a propriedade", title = "Análise de correlacao - indus, tax");p2

p3 <- ggplot(df, aes(age, nox))+
  geom_point()+
  labs(y = "Concentracao de oxido nitrico",x = "Proporcao de unidades contruidas antes de 1940", 
       title = "Análise de correlacao - nox, age");p3

p4 <- ggplot(df, aes(age,dis))+
  geom_point()+
  labs(x = "Proporcao de unidades contruidas antes de 1940",y = "distâncias ponderadas para cinco centros de emprego em Boston", 
       title = "Análise de correlacao - age, dis");p4

DIS E NOX -0.77 - Distâncias ponderadas para cinco centros de emprego em Boston possui forte relacionamento negativo com a concentracao de oxido nitrico. Ou seja, a medida que se distancia dos cinco centros de emprego em Boston, a concentracao de oxido nitrico tende a cair.

p5 <- ggplot(df, aes(dis,nox))+
  geom_point()+
  labs(x = "concentracao de oxido nitrico",y = "distâncias ponderadas para cinco centros de emprego em Boston", 
       title = "Análise de correlacao - age, dis");p5

Teremos que cogitar a hipotese de dropar algumas dessas variáveis antes da construção do modelo afim de, evitar problemas com multicolinearidade quando construirmos o modelo preditivo. Isso evita que nosso modelo se torne tendencioso, pois há correlacao entre variáveis que não envolvem a variável target (MEDV)

Além dessas variáveis multicolineares no dataset, podemos identificar as variáveis que mais se relaciona à variavel Target, que são: lstat -0,74 (% menor status da população) e rm 0,7 (número médio de quartos por habitação)

p6 <- ggplot(df) +
  aes(x = lstat, y = medv, size = rm) +
  geom_point(colour = "#0c4c8a", alpha = 0.6)+
  labs(title = "Análise de correlação - medv, lstats, rm");p6

Aqui podemos identifcar que, a media que o lstat aumenta, a taxa média aumenta, assim como a medida que a quantidade de quartos diminui, a taxa média também diminui.

DADOS DE TREINO E DE TESTE:
sample <- sample.split(df_p$crim, SplitRatio = 0.8) 
treino <- subset(df_p, sample == TRUE) # Sample de treino
teste <- subset(df_p, sample == FALSE) # Sample de teste
DROPANDO VARIÁVEIS MULTICOLINEARES:

Utilizaremos o algoritmo stepAIC para indicaçao do melhor modelo, levando-se em conta o fator da não multicolinearidade, também coletaremos o RSE do modelo (residual stard error). Isso será usado na análise do STEPAIC.

modelo1 <- lm(y~.,treino);summary(modelo1)
## 
## Call:
## lm(formula = y ~ ., data = treino)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.6563  -2.6534  -0.4251   1.4473  24.2140 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.5229     0.2398  93.933  < 2e-16 ***
## crim         -0.9917     0.2937  -3.376 0.000809 ***
## zn            1.0623     0.3590   2.960 0.003269 ** 
## indus         0.2823     0.5160   0.547 0.584709    
## chas          0.7841     0.2495   3.143 0.001799 ** 
## nox          -2.0972     0.5069  -4.137 4.31e-05 ***
## rm            2.2984     0.3328   6.906 2.03e-11 ***
## age           0.3201     0.4114   0.778 0.437047    
## dis          -3.2281     0.4843  -6.666 9.00e-11 ***
## rad           2.6765     0.6455   4.147 4.14e-05 ***
## tax          -2.0815     0.7322  -2.843 0.004709 ** 
## ptratio      -2.1063     0.3304  -6.374 5.20e-10 ***
## black         0.7821     0.2726   2.869 0.004344 ** 
## lstat        -4.4262     0.4102 -10.789  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.797 on 390 degrees of freedom
## Multiple R-squared:  0.7402, Adjusted R-squared:  0.7316 
## F-statistic: 85.48 on 13 and 390 DF,  p-value: < 2.2e-16
step(modelo1, direction = 'both',scale = 4.732^2) 
## Start:  AIC=24.81
## y ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + 
##     ptratio + black + lstat
## 
##           Df Sum of Sq     RSS      Cp
## - indus    1      6.88  8981.8  23.118
## - age      1     13.93  8988.8  23.433
## <none>                  8974.9  24.811
## - tax      1    185.96  9160.8  31.116
## - black    1    189.40  9164.3  31.269
## - zn       1    201.56  9176.5  31.813
## - chas     1    227.35  9202.2  32.964
## - crim     1    262.32  9237.2  34.526
## - nox      1    393.86  9368.7  40.400
## - rad      1    395.69  9370.6  40.482
## - ptratio  1    935.04  9909.9  64.569
## - dis      1   1022.55  9997.4  68.477
## - rm       1   1097.58 10072.5  71.828
## - lstat    1   2678.87 11653.8 142.447
## 
## Step:  AIC=23.12
## y ~ crim + zn + chas + nox + rm + age + dis + rad + tax + ptratio + 
##     black + lstat
## 
##           Df Sum of Sq     RSS      Cp
## - age      1     13.00  8994.8  21.699
## <none>                  8981.8  23.118
## + indus    1      6.88  8974.9  24.811
## - black    1    188.46  9170.2  29.535
## - zn       1    195.11  9176.9  29.832
## - tax      1    203.24  9185.0  30.195
## - chas     1    237.59  9219.4  31.729
## - crim     1    265.91  9247.7  32.994
## - nox      1    392.62  9374.4  38.652
## - rad      1    402.52  9384.3  39.094
## - ptratio  1    929.08  9910.8  62.610
## - rm       1   1090.79 10072.6  69.832
## - dis      1   1132.00 10113.8  71.673
## - lstat    1   2677.00 11658.8 140.671
## 
## Step:  AIC=21.7
## y ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
##           Df Sum of Sq     RSS      Cp
## <none>                  8994.8  21.699
## + age      1     13.00  8981.8  23.118
## + indus    1      5.96  8988.8  23.433
## - zn       1    185.39  9180.2  27.978
## - black    1    193.84  9188.6  28.356
## - tax      1    197.82  9192.6  28.534
## - chas     1    241.10  9235.9  30.466
## - crim     1    265.56  9260.3  31.559
## - nox      1    382.50  9377.3  36.781
## - rad      1    392.64  9387.4  37.234
## - ptratio  1    917.52  9912.3  60.675
## - rm       1   1204.93 10199.7  73.510
## - dis      1   1299.29 10294.1  77.724
## - lstat    1   2872.71 11867.5 147.992
## 
## Call:
## lm(formula = y ~ crim + zn + chas + nox + rm + dis + rad + tax + 
##     ptratio + black + lstat, data = treino)
## 
## Coefficients:
## (Intercept)         crim           zn         chas          nox           rm  
##     22.5098      -0.9971       1.0004       0.8032      -1.9307       2.3415  
##         dis          rad          tax      ptratio        black        lstat  
##     -3.3906       2.5250      -1.8511      -2.0657       0.7901      -4.3029

Resultado StepAIC: lm(formula = y ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat, data = treino).

O algoritmo step, nos indica que devemos dropar apenas as variáveis INDUS,AGE.

Além do que nos indica o STEPAIC, vamos listas as variaveis de maior importancia para cogitarmos dropar mais algumas delas para o modelo.

varImp(modelo1)
##            Overall
## crim     3.3762233
## zn       2.9595376
## indus    0.5469745
## chas     3.1431249
## nox      4.1370215
## rm       6.9061327
## age      0.7779889
## dis      6.6659200
## rad      4.1466440
## tax      2.8426636
## ptratio  6.3743079
## black    2.8688161
## lstat   10.7893078

Vamos cogitar criar modelos selecionando somente algumas variáveis.

AJUSTES DE DATASET PARA MODELOS OTIMIZADOS:

TREINO:

df_treino_step <- treino %>% 
  dplyr::select(crim,zn,chas,rm,rad,ptratio,black,lstat,y)
df_treino_step2 <- treino %>% 
  dplyr::select(crim,zn,chas,rm,rad,ptratio,black,lstat,y) %>% 
  mutate(rad2 = rad^2,
         rm2 = rm^2,
        lstat2 = lstat^2,
        ptratio2 = ptratio^2)
df_treino_varimp1 <- treino %>% 
  dplyr::select(lstat,ptratio,dis,rm,y)
df_treino_varimp2 <- treino %>% 
  dplyr::select(lstat,ptratio,dis,rm,rad,y)

TESTE:

df_teste_step <- teste %>% 
  dplyr::select(crim,zn,chas,rm,rad,ptratio,black,lstat,y)
df_teste_step2 <- teste %>% 
  dplyr::select(crim,zn,chas,rm,rad,ptratio,black,lstat,y) %>% 
  mutate(rad2 = rad^2,
         rm2 = rm^2,
         lstat2 = lstat^2,
         ptratio2 = ptratio^2)
df_teste_varimp1 <- teste %>% 
  dplyr::select(lstat,ptratio,dis,rm,y)
df_teste_varimp2 <- teste %>% 
  dplyr::select(lstat,ptratio,dis,rm,rad,y)
CRIAÇÃO DE MODELOS OTIMIZADOS:
modelo1 <- lm(y~.,df_treino_step, method = "lm");summary(modelo1) 
## Warning in lm(y ~ ., df_treino_step, method = "lm"): method = 'lm' is not
## supported. Using 'qr'
## 
## Call:
## lm(formula = y ~ ., data = df_treino_step, method = "lm")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.5821  -2.8779  -0.8458   1.5775  27.3969 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.4535     0.2561  87.672  < 2e-16 ***
## crim         -0.7375     0.3126  -2.359 0.018812 *  
## zn           -0.1960     0.2994  -0.655 0.513031    
## chas          0.8858     0.2628   3.371 0.000824 ***
## rm            2.8928     0.3381   8.556 2.62e-16 ***
## rad           1.0048     0.3788   2.653 0.008309 ** 
## ptratio      -2.1109     0.3275  -6.446 3.36e-10 ***
## black         0.8225     0.2916   2.821 0.005030 ** 
## lstat        -4.3193     0.3835 -11.262  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.142 on 395 degrees of freedom
## Multiple R-squared:  0.6977, Adjusted R-squared:  0.6916 
## F-statistic:   114 on 8 and 395 DF,  p-value: < 2.2e-16
modelo2 <- lm(y~.,df_treino_varimp1,method="lm");summary(modelo2) 
## Warning in lm(y ~ ., df_treino_varimp1, method = "lm"): method = 'lm' is not
## supported. Using 'qr'
## 
## Call:
## lm(formula = y ~ ., data = df_treino_varimp1, method = "lm")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.352  -2.781  -0.560   1.876  25.985 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.4920     0.2568  87.602  < 2e-16 ***
## lstat        -5.3700     0.3740 -14.356  < 2e-16 ***
## ptratio      -2.1383     0.2914  -7.339 1.22e-12 ***
## dis          -1.5608     0.3089  -5.053 6.62e-07 ***
## rm            2.6874     0.3337   8.054 9.33e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.159 on 399 degrees of freedom
## Multiple R-squared:  0.6926, Adjusted R-squared:  0.6896 
## F-statistic: 224.8 on 4 and 399 DF,  p-value: < 2.2e-16
modelo3 <- lm(y~., df_treino_varimp2, method="lm");summary(modelo3) 
## Warning in lm(y ~ ., df_treino_varimp2, method = "lm"): method = 'lm' is not
## supported. Using 'qr'
## 
## Call:
## lm(formula = y ~ ., data = df_treino_varimp2, method = "lm")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.205  -2.779  -0.549   1.849  26.212 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.4994     0.2573  87.440  < 2e-16 ***
## lstat        -5.3132     0.3878 -13.700  < 2e-16 ***
## ptratio      -2.0702     0.3159  -6.553 1.75e-10 ***
## dis          -1.6235     0.3287  -4.939 1.16e-06 ***
## rm            2.7193     0.3388   8.027 1.13e-14 ***
## rad          -0.1929     0.3437  -0.561    0.575    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.163 on 398 degrees of freedom
## Multiple R-squared:  0.6929, Adjusted R-squared:  0.689 
## F-statistic: 179.6 on 5 and 398 DF,  p-value: < 2.2e-16
modelo4 <- lm(y~.,df_treino_step2, method = "lm");summary(modelo4)
## Warning in lm(y ~ ., df_treino_step2, method = "lm"): method = 'lm' is not
## supported. Using 'qr'
## 
## Call:
## lm(formula = y ~ ., data = df_treino_step2, method = "lm")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -28.5565  -2.5106  -0.6092   1.7301  27.2510 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.1732     0.7426  27.165  < 2e-16 ***
## crim         -1.0674     0.2692  -3.966 8.70e-05 ***
## zn           -0.4189     0.2712  -1.545 0.123249    
## chas          0.8045     0.2247   3.580 0.000387 ***
## rm            1.9951     0.2982   6.690 7.75e-11 ***
## rad           1.2004     0.7325   1.639 0.102068    
## ptratio      -1.0367     0.3613  -2.870 0.004331 ** 
## black         0.7644     0.2482   3.080 0.002215 ** 
## lstat        -6.0513     0.4193 -14.430  < 2e-16 ***
## rad2         -0.3510     0.6776  -0.518 0.604696    
## rm2           0.9471     0.1310   7.232 2.52e-12 ***
## lstat2        1.2754     0.2107   6.053 3.33e-09 ***
## ptratio2      0.4840     0.2426   1.995 0.046690 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.371 on 391 degrees of freedom
## Multiple R-squared:  0.7838, Adjusted R-squared:  0.7771 
## F-statistic: 118.1 on 12 and 391 DF,  p-value: < 2.2e-16
AVALIAÇÃO DO MODELOS 1-4:
par(mfrow = c(2,2))
plot(modelo1, which = c(1:4), pch = 20)

plot(modelo2, which = c(1:4), pch = 20)

plot(modelo3, which = c(1:4), pch = 20)

plot(modelo4, which = c(1:4), pch = 20)

INTERPRETANDO OS GRÁFICOS:

Gráfico 1: Temos os resíduos em função dos valores estimados. Aqui observamos a independência e a homocedasticidade, se os resíduos se distribuem de maneira razoavelmente aleatória e com mesma amplitude em torno do zero.

Gráfico 2: Podemos avaliar a normalidade dos resíduos. A linha diagonal pontilhada representa a distribuição normal teórica, e os pontos a distribuição dos resíduos observada. Espera-se que não exista grande fuga dos pontos em relação à reta.

Gráfico 3: Pode ser avaliado da mesma maneira que o primeiro, observando a aleatoriedade e amplitude, desta vez dos resíduos padronizados.

Gráfico 4: E o último gráfico permite visualizar as Distâncias de Cook das observações, uma medida de influência que pode indicar a presença de outliers quando possui valor maior do que 1.

Avaliando os 4 modelos criados, o que parece mais promissor é o modelo 4, cujo possui 0,79 em R2, e razoavel NORMAL Q-Q. Também é possível notar que no grafico 1 (a esquerda), os residuals parecem nao possuir tendencias e se apresentam razoalvemnte de forma aleatoria.

Também podemos comparar as distribuiçoes dos residuals ao redor do zero. Desta análise podemos verificar que o modelo 4 também é o que mais possui residuals próximos a zero e a uma distribuição gaussiana.

par(mfrow = c(2,2))
hist(modelo1$residuals, labels = TRUE, ylim = c(0,250))
hist(modelo2$residuals, labels = TRUE, ylim = c(0,250))
hist(modelo3$residuals,labels = TRUE, ylim = c(0,250))
hist(modelo4$residuals,labels = TRUE, ylim = c(0,250))

FAZENDO PREVISÕES:
previsao4 <- predict(modelo4, df_teste_step2)

Podemos ainda avaliar a acuracia do modelo, comparando as médias de previsoes vs a média dos valores observados:

mean(previsao4)
## [1] 22.88649
mean(df_teste_step2$y) 
## [1] 22.93039

Comparado ao nossos dados de teste a média da previsão do modelo ficou bem próxima às observaçoes presentes no dataset de teste. O que é muito bom.