Revendo o Dataset

Olá! Seja bem vindo mais uma vez ao dataset de nascimentos e casamentos do IBGE. A base de dados é composta de 4536 observações e 6 variáveis (estado, ano, faixa etária da mulher, quantidade de nascimentos de bebês do sexo masculino, quantidade de nascimentos de bebês do sexo feminino e quantidade de casamentos. A “key” de cada observação é um conjunto de estado, ano da observação e faixa etária da mulher que casou naquele período ou faixa etária da mãe que teve bebês naquele período. Vamos fazer uma breve revisão, antes de partir para o modelo preditivo.

library(tidyverse)
library(readxl)
casamentos_nascimentosbr<-read_xlsx("C:\\Users\\lferreira\\Downloads\\live university\\Projetos\\casamentos_nascimentosbr.xlsx")

glimpse(casamentos_nascimentosbr)
## Observations: 4,536
## Variables: 6
## $ estado         <chr> "Acre", "Acre", "Acre", "Acre", "Acre", "Acre",...
## $ ano            <dbl> 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003,...
## $ classe         <chr> "15 a 19 anos", "20 a 24 anos", "25 a 29 anos",...
## $ nasc_masc      <dbl> 3049, 3692, 2143, 1211, 566, 257, 58, 16, 16, 1...
## $ nasc_fem       <dbl> 2886, 3267, 2059, 1139, 554, 212, 50, 18, 18, 1...
## $ qtd_casamentos <dbl> 410, 685, 493, 340, 184, 129, 66, 52, 27, 19, 2...
summary(casamentos_nascimentosbr)
##     estado               ano          classe            nasc_masc      
##  Length:4536        Min.   :2003   Length:4536        Min.   :    0.0  
##  Class :character   1st Qu.:2006   Class :character   1st Qu.:   34.0  
##  Mode  :character   Median :2010   Mode  :character   Median :  496.5  
##                     Mean   :2010                      Mean   : 4809.2  
##                     3rd Qu.:2013                      3rd Qu.: 5316.0  
##                     Max.   :2016                      Max.   :92918.0  
##     nasc_fem     qtd_casamentos 
##  Min.   :    0   Min.   :    0  
##  1st Qu.:   35   1st Qu.:  140  
##  Median :  475   Median :  683  
##  Mean   : 4572   Mean   : 2976  
##  3rd Qu.: 5054   3rd Qu.: 2924  
##  Max.   :87869   Max.   :76407

Lembrando que o período das observações vai entre 2003 a 2016. Essa é uma base livre de NA’s!

Vamos nessa análise elaborar um modelo preditivo que possa para cada estado prever o número de nascimentos, dependendo de variáveis, como faixa etária da mulher e quantidade de casamentos.

Vamos criar dois modelos de aprendizagem supervisionada a Regressão Múltipla e Random Forest e decidir qual modelo tem melhor capacidade para prever a quantidade de nascimentos, através de uma das mais famosas medidas de avaliação de desempenho, o SMAPE!(Symmetric Mean Absolute Percentage Error).

Para simplificar, vamos agrupar a quantidade de nascimentos de homens e mulheres em uma variável só!

casamentos_nascimentosbr_ml<-casamentos_nascimentosbr%>%
  group_by(estado,classe,ano)%>%
  summarise(nasc_masc=sum(nasc_masc),nasc_fem=sum(nasc_fem),qtd_casamentos=sum(qtd_casamentos))%>%
  mutate(qtd_nascimentos=nasc_masc+nasc_fem)
#Arrumando as variáveis
 casamentos_nascimentosbr_ml<-casamentos_nascimentosbr_ml[,c(1,2,3,6,7)]
 casamentos_nascimentosbr_ml$estado<-as.factor(casamentos_nascimentosbr_ml$estado)
 casamentos_nascimentosbr_ml$ano<-as.integer(casamentos_nascimentosbr_ml$ano)
  #trocando os valores de classe para numero ordinal
 casamentos_nascimentosbr_ml$classe<-factor(casamentos_nascimentosbr_ml$classe,levels=c(
    "15 a 19 anos","20 a 24 anos","25 a 29 anos","30 a 34 anos",
   "35 a 39 anos","40 a 44 anos","Menos de 15 anos","45 a 49 anos","50 anos ou mais"
 ))

Não sei se repararam, mas fiz algo interessante, já que a faixa “menos de 15 anos” possui índice de correlação parecido com as faixas acima de 45 anos, ou seja, sem muito padrão. Para explicar a variável target (Visto em EDA IBGE), coloquei uma próxima da outra no argumento levels. Isso vai ser crucial para o nosso modelo.

Correlação entre quantidade de casamentos e quantidade de nascimentos

casamentos_nascimentosbr_ml%>%
    ggplot(aes(x=qtd_casamentos,y=qtd_nascimentos,col=estado))+
    geom_point()+
    geom_smooth(method ="lm",se=FALSE)+
    labs(title = "Correlação entre casamentos e nascimentos",subtitle = "Observações em estado,ano e classe",
         caption = "Fonte: IBGE - Estatísticas do Registro Civil (2003-2016)",x="Qtd de Casamentos",
         y="Qtd de Nascimentos")+
    scale_y_continuous(label=scales::comma)+
    scale_x_continuous(label=scales::comma)+
    theme_light()

library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## 
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
## 
##     nasa
casamentos_nascimentosbr_ml%>%
  ggpairs(3:5)

Aparentemente sem multicolinearidade entre as variáveis.

Vamos agrupar o dataset em estado. Gosto muito quando vou elaborar vários modelos de uma vez só usar o purrr e o tidyr package! Principalmente por duas grandes e maravilhosas funções, a nest() e a map() para trabalhar com listas!

Vamos agrupar e dar uma olhada na correlação entre quantidade de nascimentos e casamentos para cada estado?

library(purrr)

casamentos_nascimentosbr_ml_nest<-casamentos_nascimentosbr_ml%>%
  group_by(estado)%>%nest()

casamentos_nascimentosbr_ml_nest%>%
  mutate(cor_pearson = map_dbl(.x=data,.f= ~cor(.$qtd_casamentos,.$qtd_nascimentos)))%>%
  ggplot(aes(y=cor_pearson))+
  geom_boxplot()+
  labs(x="",y="Correlação de Pearson",title = "BoxPlot - Correlação de Pearson",
       caption = "Fonte: IBGE - Estatísticas do Registro Civil (2003-2016)")+
  theme_dark()+
  theme(
    axis.text.x = element_blank()
  )

Uau! Que poder! Correlação excelente de um modo geral, mas vemos um outlier em nosso boxplot. Qual estado será? Vamos descobrir.

casamentos_nascimentosbr_ml_nest%>%
  mutate(cor_pearson = map_dbl(.x=data,.f= ~cor(.$qtd_casamentos,.$qtd_nascimentos)))%>%
  filter(cor_pearson<.8)
## # A tibble: 1 x 3
##   estado data               cor_pearson
##   <fct>  <list>                   <dbl>
## 1 Amapa  <tibble [126 x 4]>       0.712

Ainda assim uma correlação de 0.71 é excelente! Bem, vamos dar início à modelagem.

Modelando com Regressão Linear Múltipla

library(rsample) #Criando a partição
set.seed(42) 
casamentos_nascimentosbr_ml_split<-casamentos_nascimentosbr_ml_nest%>%
  mutate(split = map(.x=data,.f=~initial_split(.x,prop=0.75))) #vamos fazer uma separação entre treinamento e teste com uma relação 75/25

casamentos_nascimentosbr_ml_split<-casamentos_nascimentosbr_ml_split%>%
  mutate(train_CV = map(.x=split,.f=~training(.x)),test = map(.x=split,.f=~testing(.x)))

head(casamentos_nascimentosbr_ml_split)
## # A tibble: 6 x 5
##   estado   data              split          train_CV        test           
##   <fct>    <list>            <list>         <list>          <list>         
## 1 Acre     <tibble [126 x 4~ <split [95/31~ <tibble [95 x ~ <tibble [31 x ~
## 2 Alagoas  <tibble [126 x 4~ <split [95/31~ <tibble [95 x ~ <tibble [31 x ~
## 3 Amapa    <tibble [126 x 4~ <split [95/31~ <tibble [95 x ~ <tibble [31 x ~
## 4 Amazonas <tibble [126 x 4~ <split [95/31~ <tibble [95 x ~ <tibble [31 x ~
## 5 Bahia    <tibble [126 x 4~ <split [95/31~ <tibble [95 x ~ <tibble [31 x ~
## 6 Ceara    <tibble [126 x 4~ <split [95/31~ <tibble [95 x ~ <tibble [31 x ~

Beleza, agora na amostra de treinamento train_CV vamos fazer uma validação cruzada criando 5 partições e para cada partição vamos fazer a mesma divisão proporcional entre treinamento e teste.

casamentos_nascimentosbr_ml_split<-casamentos_nascimentosbr_ml_split%>%
  mutate(CV_split = map(.x=train_CV,.f=~vfold_cv(.x,v=5)))

#vamos agora dar unnest e trabalhar com um dataset da validação cruzada a parte para facilitar
# a modelagem
casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_split%>%
 unnest(CV_split)
head(casamentos_nascimentosbr_ml_CV)
## # A tibble: 6 x 3
##   estado  splits          id   
##   <fct>   <list>          <chr>
## 1 Acre    <split [76/19]> Fold1
## 2 Acre    <split [76/19]> Fold2
## 3 Acre    <split [76/19]> Fold3
## 4 Acre    <split [76/19]> Fold4
## 5 Acre    <split [76/19]> Fold5
## 6 Alagoas <split [76/19]> Fold1
#criando uma amostra de treino e uma de teste para cada partição da validação cruzada
casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_CV%>%
  mutate(
    train = map(splits, ~training(.x)), 
    validate = map(splits, ~testing(.x))
  )

head(casamentos_nascimentosbr_ml_CV)
## # A tibble: 6 x 5
##   estado  splits          id    train             validate         
##   <fct>   <list>          <chr> <list>            <list>           
## 1 Acre    <split [76/19]> Fold1 <tibble [76 x 4]> <tibble [19 x 4]>
## 2 Acre    <split [76/19]> Fold2 <tibble [76 x 4]> <tibble [19 x 4]>
## 3 Acre    <split [76/19]> Fold3 <tibble [76 x 4]> <tibble [19 x 4]>
## 4 Acre    <split [76/19]> Fold4 <tibble [76 x 4]> <tibble [19 x 4]>
## 5 Acre    <split [76/19]> Fold5 <tibble [76 x 4]> <tibble [19 x 4]>
## 6 Alagoas <split [76/19]> Fold1 <tibble [76 x 4]> <tibble [19 x 4]>

Beleza, tudo certo. Partições criadas, amostras de treinamento e testes prontos. Agora, vamos para a criação do modelo de regressão!

casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_CV%>%
  mutate(model = map(.x=train,.f=~lm(qtd_nascimentos~.,data=.x)))
head(casamentos_nascimentosbr_ml_CV)
## # A tibble: 6 x 6
##   estado  splits          id    train             validate          model  
##   <fct>   <list>          <chr> <list>            <list>            <list> 
## 1 Acre    <split [76/19]> Fold1 <tibble [76 x 4]> <tibble [19 x 4]> <S3: l~
## 2 Acre    <split [76/19]> Fold2 <tibble [76 x 4]> <tibble [19 x 4]> <S3: l~
## 3 Acre    <split [76/19]> Fold3 <tibble [76 x 4]> <tibble [19 x 4]> <S3: l~
## 4 Acre    <split [76/19]> Fold4 <tibble [76 x 4]> <tibble [19 x 4]> <S3: l~
## 5 Acre    <split [76/19]> Fold5 <tibble [76 x 4]> <tibble [19 x 4]> <S3: l~
## 6 Alagoas <split [76/19]> Fold1 <tibble [76 x 4]> <tibble [19 x 4]> <S3: l~

Vamos ver como ficou nosso R² ajustado, para cada modelo:

library(broom)
casamentos_nascimentosbr_ml_CV%>%
  mutate(
    rsquad = map_dbl(.x=model,.f=~glance(.x)[["adj.r.squared"]])
  )%>%
  ggplot(aes(y=rsquad))+
  geom_boxplot()+
  labs(x="",y="R² Ajustado",title = "BoxPlot - R² Ajustado distribuição por modelo",
       caption = "Fonte: IBGE - Estatísticas do Registro Civil (2003-2016)")+
  theme_dark()+
  theme(
    axis.text.x = element_blank()
  )

Assombroso! Parece que as variáveis do dataset estão conseguindo explicar de forma excelente a variação no fator quantidade de nascimentos. Mas será que o modelo vai ser capaz de prever bem nas amostras de teste da validação cruzada?

Mas antes disso, será que podemos melhorar ainda mais o nosso modelo? Gosto sempre quando trabalho com regressão, pois tento aperfeiçoar o modelo, minimizando o valor do critério de Akaike com a stepAIC() função do MASS package.

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_CV%>%
  mutate(model_otimo = map(.x=train,.f=~(lm(qtd_nascimentos~.,data=.x)%>%
  stepAIC(trace=0))
  ))

casamentos_nascimentosbr_ml_CV%>%
  mutate(
    rsquad = map_dbl(.x=model_otimo,.f=~glance(.x)[["adj.r.squared"]])
  )%>%
  ggplot(aes(y=rsquad))+
  geom_boxplot()+
  labs(x="",y="R² Ajustado",title = "BoxPlot - R² Ajustado distribuição por modelo otimo",
       caption = "Fonte: IBGE - Estatísticas do Registro Civil (2003-2016)")+
  theme_dark()+
  theme(
    axis.text.x = element_blank()
  )

Como o modelo já tinha um nível muito bom, não tivemos uma diferença significativa entre eles.

Mas vamos seguir com ambos os modelos para validação com o SMAPE!

casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_CV%>%
  mutate(
    
    validate_actual = map(validate, ~.x$qtd_nascimentos),
    
    validate_predicted = map2(.x = model_otimo, .y = validate, ~predict(.x, .y))
    
  )

head(casamentos_nascimentosbr_ml_CV)
## # A tibble: 6 x 9
##   estado splits id    train validate model model_otimo validate_actual
##   <fct>  <list> <chr> <lis> <list>   <lis> <list>      <list>         
## 1 Acre   <spli~ Fold1 <tib~ <tibble~ <S3:~ <S3: lm>    <dbl [19]>     
## 2 Acre   <spli~ Fold2 <tib~ <tibble~ <S3:~ <S3: lm>    <dbl [19]>     
## 3 Acre   <spli~ Fold3 <tib~ <tibble~ <S3:~ <S3: lm>    <dbl [19]>     
## 4 Acre   <spli~ Fold4 <tib~ <tibble~ <S3:~ <S3: lm>    <dbl [19]>     
## 5 Acre   <spli~ Fold5 <tib~ <tibble~ <S3:~ <S3: lm>    <dbl [19]>     
## 6 Alago~ <spli~ Fold1 <tib~ <tibble~ <S3:~ <S3: lm>    <dbl [19]>     
## # ... with 1 more variable: validate_predicted <list>
library(Metrics)
casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_CV%>%
  mutate(validate_smape = map2_dbl(
    .x=validate_actual,.y= validate_predicted,.f= ~smape(actual = .x, predicted = .y)))

casamentos_nascimentosbr_SMAPE<-casamentos_nascimentosbr_ml_CV%>%
  group_by(estado)%>%
  summarise(SMAPE = mean(validate_smape))

library(knitr)

kable(casamentos_nascimentosbr_SMAPE)
estado SMAPE
Acre 0.3887492
Alagoas 0.5712747
Amapa 0.2981771
Amazonas 0.2281959
Bahia 0.3931730
Ceara 0.3372935
Distrito Federal 0.5070146
Espirito Santo 0.4344515
Goias 0.4022064
Maranhao 0.5761041
Mato Grosso 0.5559720
Mato Grosso do Sul 0.2549416
Minas Gerais 0.4231878
Para 0.4897178
Paraiba 0.5985835
Parana 0.3032747
Pernambuco 0.5397053
Piaui 0.5579125
Rio de Janeiro 0.3641744
Rio Grande do Norte 0.5528111
Rio Grande do Sul 0.3324656
Rondonia 0.4961022
Roraima 0.3739956
Santa Catarina 0.4287751
Sao Paulo 0.3742262
Sergipe 0.4243968
Tocantins 0.6402341

Ops! Tivemos um overfitting, parece que na amostra de treinamento, os fatores conseguem explicar muito bem a variação do fator quantidade de nascimentos, mas para previsão fora da amostra o desempenho ficou muito abaixo do esperado!

Vamos tentar então simplificar o modelo, vamos usar somente as variaveis qtd_casamentos e classe, como variáveis explicativas dessa vez!

casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_CV%>%
  mutate(model2 = map(.x=train,.f=~lm(qtd_nascimentos~qtd_casamentos+classe,data=.x)))

casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_CV%>%
  mutate(
    
    
    
    validate_predicted2 = map2(.x = model2, .y = validate, ~predict(.x, .y))
    
  )


casamentos_nascimentosbr_ml_CV<-casamentos_nascimentosbr_ml_CV%>%
  mutate(validate_smape2 = map2_dbl(
    .x=validate_actual,.y= validate_predicted2,.f= ~smape(actual = .x, predicted = .y)))

casamentos_nascimentosbr_SMAPE2<-casamentos_nascimentosbr_ml_CV%>%
  group_by(estado)%>%
  summarise(SMAPE2 = mean(validate_smape2))

kable(casamentos_nascimentosbr_SMAPE2)
estado SMAPE2
Acre 0.1595873
Alagoas 0.2918275
Amapa 0.1024682
Amazonas 0.2281959
Bahia 0.3982810
Ceara 0.2408246
Distrito Federal 0.3484189
Espirito Santo 0.4613920
Goias 0.3935993
Maranhao 0.2459672
Mato Grosso 0.4792189
Mato Grosso do Sul 0.2568969
Minas Gerais 0.3032223
Para 0.3121873
Paraiba 0.3066841
Parana 0.3465375
Pernambuco 0.1713558
Piaui 0.4822571
Rio de Janeiro 0.3196625
Rio Grande do Norte 0.3531862
Rio Grande do Sul 0.4519925
Rondonia 0.3262205
Roraima 0.2561849
Santa Catarina 0.4305370
Sao Paulo 0.3216776
Sergipe 0.3051049
Tocantins 0.3521932 
analise_SMAPE<-full_join(casamentos_nascimentosbr_SMAPE,casamentos_nascimentosbr_SMAPE2,by="estado")

analise_SMAPE%>%gather(Modelo,Valor,2:3)%>%
  ggplot(aes(x=estado,y=Valor,fill=Modelo))+
  geom_col(position = "dodge")+
  geom_hline(yintercept = .2,linetype = "dashed",col = "black",size=1)+
  coord_flip()

Parece que o modelo 2 mais simples, teve um efeito preditivo melhor.

Vamos finalizar fazendo calculo na amostra de teste que separamos no inicio da modelagem.

casamentos_nascimentosbr_ml_reg<-full_join(
 casamentos_nascimentosbr_ml_split,
  analise_SMAPE,
  by="estado")

#facilitar a leitura
casamentos_nascimentosbr_ml_reg<-casamentos_nascimentosbr_ml_reg[,c(1,4,5,7,8)]
 
head(casamentos_nascimentosbr_ml_reg)
## # A tibble: 6 x 5
##   estado   train_CV          test              SMAPE SMAPE2
##   <fct>    <list>            <list>            <dbl>  <dbl>
## 1 Acre     <tibble [95 x 4]> <tibble [31 x 4]> 0.389  0.160
## 2 Alagoas  <tibble [95 x 4]> <tibble [31 x 4]> 0.571  0.292
## 3 Amapa    <tibble [95 x 4]> <tibble [31 x 4]> 0.298  0.102
## 4 Amazonas <tibble [95 x 4]> <tibble [31 x 4]> 0.228  0.228
## 5 Bahia    <tibble [95 x 4]> <tibble [31 x 4]> 0.393  0.398
## 6 Ceara    <tibble [95 x 4]> <tibble [31 x 4]> 0.337  0.241
#criando o modelo melhor para cada estado
casamentos_nascimentosbr_ml_reg<-casamentos_nascimentosbr_ml_reg%>%
  mutate(best_model = if_else(
    
  SMAPE<SMAPE2,"MODEL1","MODEL2"  
    
  ))


casamentos_nascimentosbr_ml_reg<-casamentos_nascimentosbr_ml_reg%>%
  mutate(model =if_else(best_model=="MODEL1",
   map(.x=train_CV,.f=~(lm(
   qtd_nascimentos~.,data=.x)%>%      
   stepAIC(trace=0))),
   map(.x=train_CV,
  .f=~lm(qtd_nascimentos~qtd_casamentos+classe,data=.x))))
                        
casamentos_nascimentosbr_ml_reg<-casamentos_nascimentosbr_ml_reg%>%
  mutate(
    
    validate_actual = map(test, ~.x$qtd_nascimentos),
    
    validate_predicted = map2(.x = model, .y = test, ~predict(.x, .y))
    
  )                        
                        
                        
casamentos_nascimentosbr_ml_reg<-casamentos_nascimentosbr_ml_reg%>%
  mutate(SMAPE_final = map2_dbl(
    .x=validate_actual,.y= validate_predicted,.f= ~smape(actual = .x, predicted = .y)))                       
                        
                        
 ggplot(casamentos_nascimentosbr_ml_reg,aes(x=estado,y=SMAPE_final))+
   geom_col()+
   geom_hline(yintercept = .2,linetype="dashed",col="red")+
   coord_flip()

  casamentos_nascimentosbr_ml_reg%>%
   count(SMAPE_final<=.2)%>%
    kable()
SMAPE_final <= 0.2 n
FALSE 23
TRUE 4

Bem, amenizamos um pouco a situação, mas ainda temos 23 modelos com um SMAPE maior do que 20%. Inaceitável!

Será que conseguimos algo melhor com o Random Forest?

Modelando com Random Forest (Floresta Aleatória)

Vamos manter a mesma partição da validação cruzada usada na regressão. No nosso modelo, para facilitar o processamento vamos deixar o número de árvores em 100. Acredito que já seja um valor razoável para criação de um bom modelo. Para o parâmetro mtry, vamos inserir todas as possibilidades de um a três fatores, usando a função crossing() assim vamos saber qual o parâmetro mtry terá o melhor desempenho.

library(ranger)
  casamentos_nascimentosbr_ml_CV2<-casamentos_nascimentosbr_ml_CV[,1:5]
  
  casamentos_nascimentosbr_ml_tune<-casamentos_nascimentosbr_ml_CV2%>%
    crossing(mtry=1:3)

  #random forest para cada combinação de estado, fold e mtry
  casamentos_nascimentosbr_ml_tune<-casamentos_nascimentosbr_ml_tune%>%
    mutate(
      model = map2(.x=train,.y=mtry,.f= ~ranger(formula=qtd_nascimentos~.,data= .x,
      num.trees=100,seed = 42,mtry = .y))
   
    )
  
  casamentos_nascimentosbr_ml_tune<-casamentos_nascimentosbr_ml_tune%>%
    mutate(
      validate_predicted = map2(.x = model, .y = validate, ~predict(.x, .y)$predictions))
  
  casamentos_nascimentosbr_ml_tune<-casamentos_nascimentosbr_ml_tune%>%
    mutate(
      validate_actual = map(.x=validate,.f= ~.x$qtd_nascimentos))
 
  casamentos_nascimentosbr_ml_tune<-casamentos_nascimentosbr_ml_tune%>%
    mutate(validate_smape = map2_dbl(.x=validate_actual,
   .y=validate_predicted, .f= ~smape(actual = .x, predicted = .y)))
  
  casamentos_nascimentosbr_ml_rf<-casamentos_nascimentosbr_ml_tune%>%
    group_by(estado,mtry)%>%
    summarise(SMAPE = mean(validate_smape))
  
  
  kable(casamentos_nascimentosbr_ml_rf)
estado mtry SMAPE
Acre 1 0.3589207
Acre 2 0.1748306
Acre 3 0.1395432
Alagoas 1 0.3611932
Alagoas 2 0.1563820
Alagoas 3 0.1311870
Amapa 1 0.4470668
Amapa 2 0.1810180
Amapa 3 0.1005982
Amazonas 1 0.2980661
Amazonas 2 0.1835755
Amazonas 3 0.1430832
Bahia 1 0.3672615
Bahia 2 0.1383488
Bahia 3 0.0892367
Ceara 1 0.3687329
Ceara 2 0.1077914
Ceara 3 0.0720162
Distrito Federal 1 0.3744352
Distrito Federal 2 0.1673408
Distrito Federal 3 0.1406324
Espirito Santo 1 0.4180209
Espirito Santo 2 0.1943504
Espirito Santo 3 0.1405170
Goias 1 0.4384444
Goias 2 0.1527381
Goias 3 0.0843674
Maranhao 1 0.3243044
Maranhao 2 0.2044007
Maranhao 3 0.1757804
Mato Grosso 1 0.4739290
Mato Grosso 2 0.2018041
Mato Grosso 3 0.1653786
Mato Grosso do Sul 1 0.3662247
Mato Grosso do Sul 2 0.1573716
Mato Grosso do Sul 3 0.1266351
Minas Gerais 1 0.3915629
Minas Gerais 2 0.1470130
Minas Gerais 3 0.0876439
Para 1 0.3490172
Para 2 0.1415865
Para 3 0.0867624
Paraiba 1 0.3787655
Paraiba 2 0.1564813
Paraiba 3 0.1243608
Parana 1 0.3713845
Parana 2 0.1646090
Parana 3 0.1115468
Pernambuco 1 0.3895686
Pernambuco 2 0.1673348
Pernambuco 3 0.0910034
Piaui 1 0.3815494
Piaui 2 0.1489806
Piaui 3 0.1238502
Rio de Janeiro 1 0.4123803
Rio de Janeiro 2 0.1659959
Rio de Janeiro 3 0.0811324
Rio Grande do Norte 1 0.3659522
Rio Grande do Norte 2 0.1654939
Rio Grande do Norte 3 0.0978601
Rio Grande do Sul 1 0.4407174
Rio Grande do Sul 2 0.2342263
Rio Grande do Sul 3 0.0786478
Rondonia 1 0.4109574
Rondonia 2 0.2701019
Rondonia 3 0.2176122
Roraima 1 0.3912241
Roraima 2 0.2455535
Roraima 3 0.2264044
Santa Catarina 1 0.4084320
Santa Catarina 2 0.2382095
Santa Catarina 3 0.1536827
Sao Paulo 1 0.4013597
Sao Paulo 2 0.2035089
Sao Paulo 3 0.0833817
Sergipe 1 0.3623651
Sergipe 2 0.1579764
Sergipe 3 0.0964302
Tocantins 1 0.4216540
Tocantins 2 0.1806331
Tocantins 3 0.1585769 
   casamentos_nascimentosbr_ml_rf2<-casamentos_nascimentosbr_ml_tune%>%
    group_by(mtry)%>%
    summarise(SMAPE = mean(validate_smape))
   
   kable(casamentos_nascimentosbr_ml_rf2)
mtry SMAPE
1 0.3879070
2 0.1780614
3 0.1232545

Hum, parece que o parâmetro mtry(3) foi o melhor entre os demais. Vamos continuar com a amostra de teste final.

casamentos_nascimentos_br_rf<-casamentos_nascimentosbr_ml_split[,c(1,4,5)]
  
  casamentos_nascimentos_br_rf<-casamentos_nascimentos_br_rf%>%
   mutate(
     
     model = map(.x=train_CV,.f= ~ranger(formula=qtd_nascimentos~.,data= .x,
                                               num.trees=100,seed = 42,mtry = 3))
     
   )
       
  casamentos_nascimentos_br_rf<-casamentos_nascimentos_br_rf%>%
    mutate(
      validate_predicted = map2(.x = model, .y = test, ~predict(.x, .y)$predictions))
  casamentos_nascimentos_br_rf<-casamentos_nascimentos_br_rf%>%
    mutate(
      validate_actual = map(.x=test,.f= ~.x$qtd_nascimentos))     
 
    
  casamentos_nascimentos_br_rf<-casamentos_nascimentos_br_rf%>%
    mutate(validate_smape = map2_dbl(.x=validate_actual,
                                    .y=validate_predicted, .f= ~smape(actual = .x, predicted = .y)))

#comparando os MAPE de Regressão multipla e random forest
  
  smape_ml_models<-full_join(casamentos_nascimentosbr_ml_reg[,c(1,10)],
  casamentos_nascimentos_br_rf[,c(1,7)],by="estado")
    colnames(smape_ml_models)<-c("estado","smape_regressao","smape_random_forest")
  smape_ml_models2<-smape_ml_models%>%
    gather(modelo,smape,2:3)
  ggplot(smape_ml_models2,aes(x=estado,y=smape,fill=modelo))+
    geom_col(position = "dodge")+
    geom_hline(yintercept = .2,linetype="dashed",col="red",size=1)+
    coord_flip()

  kable(smape_ml_models)
estado smape_regressao smape_random_forest
Acre 0.1719814 0.1007378
Alagoas 0.2613108 0.1562917
Amapa 0.1384032 0.1333460
Amazonas 0.2802289 0.1339450
Bahia 0.7389492 0.1326222
Ceara 0.2744310 0.0823935
Distrito Federal 0.4457913 0.1407420
Espirito Santo 0.2810086 0.0641352
Goias 0.4321118 0.0691178
Maranhao 0.2156510 0.1449370
Mato Grosso 0.3324441 0.1373863
Mato Grosso do Sul 0.2798702 0.0994778
Minas Gerais 0.4489222 0.1416210
Para 0.3266229 0.0853473
Paraiba 0.2707645 0.1102267
Parana 0.3571254 0.1230340
Pernambuco 0.1320879 0.0910613
Piaui 0.4885375 0.1200653
Rio de Janeiro 0.4097590 0.0680517
Rio Grande do Norte 0.4838143 0.0891332
Rio Grande do Sul 0.3658414 0.1058016
Rondonia 0.4347397 0.1958627
Roraima 0.1405856 0.1474623
Santa Catarina 0.2438924 0.0857926
Sao Paulo 0.6021793 0.2003700
Sergipe 0.2844346 0.0601751
Tocantins 0.3809025 0.1497633 
  smape_ml_models%>%
    count(smape_regressao<smape_random_forest)
## # A tibble: 2 x 2
##   `smape_regressao < smape_random_forest`     n
##   <lgl>                                   <int>
## 1 FALSE                                      26
## 2 TRUE                                        1
  smape_ml_models2%>%
    filter(modelo == "smape_regressao")%>%
    count(smape<=.2)
## # A tibble: 2 x 2
##   `smape <= 0.2`     n
##   <lgl>          <int>
## 1 FALSE             23
## 2 TRUE               4
  smape_ml_models2%>%
    filter(modelo == "smape_random_forest")%>%
    count(smape<=.2)
## # A tibble: 2 x 2
##   `smape <= 0.2`     n
##   <lgl>          <int>
## 1 FALSE              1
## 2 TRUE              26
  mean(smape_ml_models$smape_random_forest)
## [1] 0.1173667

Excelente, exceto pelo estado de Roraima, o modelo de Random Forest foi superior! Mesmo assim a diferença do estado de Roraima foi muito pequena e podemos dizer que tivemos para todos os estados com o modelo de Random Forest um SMAPE menor do que 21% com média de 11,7%!

Como fizemos uma análise exploratória de dados no artigo anterior, descobrimos que para faixas etárias “menos de 15 anos” e para classes acima de 45 anos possuem comportamento muito parecido. Não temos tanto padrão só com número de casamentos, o que faz todo sentido, correto? Precisaríamos de mais variáveis para tentar explicar o número de nascimentos para essas classes. Poderíamos até agrupar essas classes em uma só, para nossa modelagem (ex: Faixas etárias fora do intervalo 15-44). Mas achei melhor deixá-las próximas no argumento levels(). Sem esse conhecimento de “ordem”, se colocássemos o argumento classe, somente com ordem de idade, obteríamos um modelo muito pior, com média do smape em torno de 50%! Por este motivo, sempre acho válido analisar bem o dataset antes para descobrir insights valiosos para modelagem!

Final: Previsão para 2017

Vamos agora importar o dataset de casamentos_nascimentosbr_2017 para incrementar nossa base de dados atual, e ver o quanto nosso modelo acerta para o ano de 2017!

casamentos_nascimentosbr_2017<-read_xlsx("C:\\Users\\lferreira\\Downloads\\live university\\casamentos_nascimentosbr_2017.xlsx")

casamentos_nascimentosbr_2017<- casamentos_nascimentosbr_2017%>%
  mutate(
    estado = as.factor(estado),
    classe = factor(classe,levels = c( "15 a 19 anos","20 a 24 anos","25 a 29 anos","30 a 34 anos",
   "35 a 39 anos","40 a 44 anos","Menos de 15 anos","45 a 49 anos","50 anos ou mais"
 )),
    ano = as.integer(ano)  
  )

casamentos_nascimentosbr_ml2<-bind_rows(casamentos_nascimentosbr_ml,casamentos_nascimentosbr_2017)

casamentos_nascimentosbr_ml2_nest<-casamentos_nascimentosbr_ml2%>%
  group_by(estado)%>%nest()

casamentos_nascimentosbr_ml2_nest<-casamentos_nascimentosbr_ml2_nest%>%
  mutate(train = map(.x=data,.f=~.x%>%
        filter(ano<2017)),
        validate = map(.x=data,.f=~.x%>%
        filter(ano>=2017))
                     )

casamentos_nascimentosbr_ml2_nest<-casamentos_nascimentosbr_ml2_nest%>%
  mutate( model = map(.x=train,.f=~ranger(formula=qtd_nascimentos~.,data= .x,
                                               num.trees=100,seed = 42,mtry = 3)))

casamentos_nascimentosbr_ml2_nest<-casamentos_nascimentosbr_ml2_nest%>%
  mutate( validate_predicted = map2(.x = model, .y = validate, ~predict(.x, .y)$predictions),
      validate_actual = map(.x=validate,.f= ~.x$qtd_nascimentos))

casamentos_nascimentosbr_ml2_nest<-casamentos_nascimentosbr_ml2_nest%>%
    mutate(validate_smape = map2_dbl(.x=validate_actual,
                                    .y=validate_predicted, .f= ~smape(actual = .x, predicted = .y)))

casamentos_nascimentosbr_ml2_result<-casamentos_nascimentosbr_ml2_nest[,c(1,8)]

 casamentos_nascimentosbr_ml2_result%>%
  kable()
estado validate_smape
Acre 0.3075230
Alagoas 0.3329709
Amapa 0.3062975
Amazonas 0.3999786
Bahia 0.2596776
Ceara 0.3059167
Distrito Federal 0.2922559
Espirito Santo 0.3011266
Goias 0.3084185
Maranhao 0.3512116
Mato Grosso 0.3725225
Mato Grosso do Sul 0.3492735
Minas Gerais 0.2480678
Para 0.4205127
Paraiba 0.2303846
Parana 0.3050680
Pernambuco 0.3033716
Piaui 0.2682199
Rio de Janeiro 0.2337834
Rio Grande do Norte 0.3042838
Rio Grande do Sul 0.2649118
Rondonia 0.2856358
Roraima 0.3543977
Santa Catarina 0.2950892
Sao Paulo 0.2328368
Sergipe 0.3018810
Tocantins 0.3072274 
analise_2017<-casamentos_nascimentosbr_ml2_nest[,c(1,6,7)]%>%
  unnest(validate_actual,validate_predicted)%>%
  group_by(estado)%>%
  summarise(Real = sum(validate_actual),Previsto = sum(validate_predicted))%>%
  gather(Tipo,Valor,2:3)

  ggplot(analise_2017,aes(x=estado,y=Valor,fill=Tipo))+
  geom_col(position = "dodge")+
  labs(x="Estado",y="Qtd_Nascimentos",title = "Real x Previsto 2017",subtitle= "Usando o
       modelo de Random Forest, mtry = 3, num.trees = 100",
       caption = "Fonte: IBGE - Estatísticas do Registro Civil (2003-2016)")+
    scale_y_continuous(label=scales::comma)+
  coord_flip()

  casamentos_nascimentosbr_ml2_nest[,c(1,6,7)]%>%
  unnest(validate_actual,validate_predicted)%>%
  group_by(estado)%>%
  summarise(Real = sum(validate_actual),Previsto = sum(validate_predicted))%>%
    kable()
estado Real Previsto
Acre 15485 17146.217
Alagoas 49761 49214.526
Amapa 14190 16230.002
Amazonas 71508 81577.706
Bahia 201243 207638.049
Ceara 125243 131107.304
Distrito Federal 44613 45322.897
Espirito Santo 55999 52887.074
Goias 93016 88091.148
Maranhao 107889 114149.652
Mato Grosso 52258 57247.472
Mato Grosso do Sul 44518 45418.918
Minas Gerais 260630 257913.451
Para 125666 147683.171
Paraiba 56024 55034.058
Parana 156208 154528.146
Pernambuco 129168 131462.686
Piaui 46815 44701.125
Rio de Janeiro 216315 214493.788
Rio Grande do Norte 45926 48791.979
Rio Grande do Sul 140185 141701.828
Rondonia 27543 26990.190
Roraima 9875 9611.608
Santa Catarina 98136 95526.083
Sao Paulo 611097 602944.658
Sergipe 32971 33450.128
Tocantins 24275 23472.081

Bem próximo da realidade não acha?

Chegamos ao fim! Não criamos o melhor modelo do mundo, mas foi bastante eficiente!

Acredito que possamos explorar mais esses dados na busca de variáveis que ajudam a explicar melhor a quantidade de nascimentos, principalmente para classes fora do arco 15-44 anos.

Mas foi bacana não é?

Obrigado por ter tido paciência para chegar até aqui! Foi uma longa jornada!

Abraços! Sucesso e até a próxima!