Stone Data Challenge - Previsão de TPV
Por Max Mitteldorf
Índice
Introdução
Neste pequeno html buscarei fazer um resumo sobre o tratamento que conferi às bases de dados do desafio. Comentarei tanto sobre as bases de clientes e faturamento da Stone, quanto das bases adicionais que julguei interessante adicionar à análise.
Limpeza dos dados dos clientes da Stone
O tratamento da base de dados cadastrais sobre os clientes foi meu primeiro passo. Não descreverei todas as modificações (que estão disponíveis no script “stone_data_cleaning.R”), mas buscarei tratar dos aspectos mais importantes.
Meu primeiro passo foi importar todos os pacotes que utilizaria no processo:
library(tidyverse)
library(skimr)
library(stringi)
library(data.table)
library(lubridate)
Em seguida, importei as bases de dados referentes aos clientes e faturamentos:
### Importing dataframes
df_cad = read_csv("data/raw/cadastrais.csv")
df_data = read_csv("data/raw/tpv-mensais-treinamento.csv")
Antes de prosseguir, juntei as bases de dados de cadastro e de faturamento:
#### Merging dataframes
df_cad = distinct(df_cad)
df_total = inner_join(df_cad, df_data, by="id")
Em seguida, podemos investigar os tipos das variáveis e, também, dar uma olhada nas primeiras linhas da base de dados:
#### Merging dataframes
options(dplyr.width = Inf)
head(df_total, 2)## # A tibble: 2 x 14
## id StoneCreatedDate StoneFirstTransactionDate MCC MacroClassificacao
## <dbl> <dttm> <dbl> <dbl> <chr>
## 1 76217 2019-12-13 00:00:00 20191220 7531 Serviços
## 2 76217 2019-12-13 00:00:00 20191220 7531 Serviços
## segmento sub_segmento
## <chr> <chr>
## 1 Autopeças e Serviços Automotivos Oficinas Automotivas
## 2 Autopeças e Serviços Automotivos Oficinas Automotivas
## persona porte TPVEstimate tipo_documento Estado
## <chr> <chr> <dbl> <chr> <chr>
## 1 SMB - Pequeno Porte e Ticket Medio 10k-25k 15000 PJ SP
## 2 SMB - Pequeno Porte e Ticket Medio 10k-25k 15000 PJ SP
## mes_referencia TPV_mensal
## <dbl> <dbl>
## 1 20191231 500
## 2 20200131 6601
O primeiro problema que podemos ver na base é a formatação errada das datas:
#### Merging dataframes
head(df_total$mes_referencia, 2)## [1] 20191231 20200131Resolver tal problema não é difícil:
df_total$mes_referencia = ymd(df_total$mes_referencia)
df_total$StoneFirstTransactionDate = ymd(df_total$StoneFirstTransactionDate)
head(df_total$mes_referencia, 2)## [1] "2019-12-31" "2020-01-31"
O próximo grande problema está na coluna “Estado”, que contém mais de 61 valores únicos:
### Investigating why there are 61 unique values for "Estado"
unique(df_total$Estado) # both fullnames and initials are appearing for each state## [1] "SP" "MT" "RJ"
## [4] "MG" "RS" "GO"
## [7] "MA" "DF" "SC"
## [10] "PR" "CE" "PE"
## [13] "Minas Gerais" "PA" "BA"
## [16] NA "PB" "ES"
## [19] "AM" "Rio Grande do Norte" "Santa Catarina"
## [22] "RR" "SE" "PI"
## [25] "RO" "MS" "RN"
## [28] "Pará" "TO" "AC"
## [31] "São Paulo" "Rio de Janeiro" "Bahia"
## [34] "AL" "Rio Grande do Sul" "Goiás"
## [37] "Paraná" "Ceará" "Espírito Santo"
## [40] "Pernambuco" "AP" "Roraima"
## [43] "Rondônia" "Maranhão" "Mato Grosso"
## [46] "Paraíba" "Sao Paulo" "Mato Grosso do Sul"
## [49] "Sergipe" "Amazonas" "Tocantins"
## [52] "Distrito Federal" "Acre" "Amapá"
## [55] "sc" "Piauí" "Alagoas"
## [58] "Rj" "Parana" "Goias"
## [61] "Sc" "Sp"Para corrigir isso utilizarei um dataframe extra com siglas e nomes de estados e, também, manipularei o padrão das strings:
df_siglas = read_csv("data/raw/sigla_estados.csv") # Dataframe w/ state icons
## Fixing "Estado" variable
df_siglas = mutate(df_siglas, NOME = paste0(toupper(NOME),"$")) %>% # Converting to lowercase
mutate(NOME = stri_trans_general(NOME,"Latin-ASCII"))# removing accents
## Creating named list with replacements
siglas = df_siglas$SIGLA
names(siglas) = df_siglas$NOME
df_total = mutate(df_total, Estado = toupper(Estado)) %>% # Converting to uppercase
mutate(Estado = stri_trans_general(Estado,"Latin-ASCII")) # removing accents
## Replacing strings
df_total$Estado = str_replace_all(df_total$Estado, siglas)
unique(df_total$Estado)## [1] "SP" "MT" "RJ" "MG" "RS" "GO" "MA" "DF" "SC" "PR" "CE" "PE" "PA" "BA" NA
## [16] "PB" "ES" "AM" "RN" "RR" "SE" "PI" "RO" "MS" "TO" "AC" "AL" "AP"
O próximo problema pode ser encontrado na coluna porte, que contém em cada célula uma string com um range numérico:
## Changing porte datatype to number
head(df_total$porte, 3)## [1] "10k-25k" "10k-25k" "10k-25k"Dessa forma, podemos manipular a string para tornar o range numérico e calcular sua média:
## Changing porte datatype to number
df_total = mutate(df_total, porte=str_replace_all(porte, "2.5k", "2500"))
df_total = mutate(df_total, porte=str_replace_all(porte, "k", "000"))
df_total = mutate(df_total, porte=str_replace_all(porte, "\\+", ""))
df_total = mutate(df_total, porte=str_replace_all(porte, "-", "+"))
eval_parse = function (x) {eval(parse(text = x))}
df_total$porte = sapply(df_total$porte, eval_parse)/2
head(df_total$porte, 3)## 10000+25000 10000+25000 10000+25000
## 17500 17500 17500
Para encerrar a parte de limpeza de dados, vou descrever o maior problema que encontrei na base. Durante minhas manipulações, descobri que o identificador único de clientes possuía entradas repetidas para as mesmas datas. Isso ocorre porque alguns clientes tiverem algumas variáveis alteradas durante o tempo. Dessa forma, haviam clientes com mais de uma Macro Classificação, mais de um segmento, mais de uma estimativa de TPV, entre outros problemas:
### The dataframe has different entries for the same id's
distinct_values = df_total %>%
as_tibble() %>%
group_by(id) %>%
summarise(n_mcc = n_distinct(MCC), n_macro=n_distinct(MacroClassificacao),
n_segmento=n_distinct(segmento), n_subsegmento=n_distinct(sub_segmento),
n_persona=n_distinct(persona), n_porte=n_distinct(porte),
n_tpvestimate = n_distinct(TPVEstimate), n_documento=n_distinct(tipo_documento),
n_estado=n_distinct(Estado)) Vamos observar alguns ids com algumas dessas variáveis que se alteram:
## These ids have different caracteristics
ids_errados = distinct_values %>%
filter(n_mcc>1 | n_macro>1 | n_segmento>1 | n_subsegmento>1 | n_persona>1 |
n_porte>1 | n_tpvestimate>1 | n_documento>1 | n_estado>1)
head(ids_errados, 4)## # A tibble: 4 x 10
## id n_mcc n_macro n_segmento n_subsegmento n_persona n_porte n_tpvestimate
## <dbl> <int> <int> <int> <int> <int> <int> <int>
## 1 21 2 2 2 2 2 1 2
## 2 29 1 1 1 1 1 2 2
## 3 30 1 1 1 1 1 1 2
## 4 32 1 1 1 1 2 2 2
## n_documento n_estado
## <int> <int>
## 1 1 1
## 2 1 1
## 3 1 1
## 4 1 1
Agora, podemos definir alguma estratégia para remover esses valores problemáticos. Meu critério de decisão será priorizar entradas com valor válido para estado, maior estimativa de TPV e maior porte. Julgo que esse critério selecionará as entradas mais recentes para cada id.
Esse é o código que utilizei para resolver isso. Não comentarei ele em detalhes pois acredito ser algo um pouco cansativo:
## Defining a strategy to keep the most recent entries for each id
linhas_certas = tibble()
contador = 1
for (id_problematico in ids_errados$id)
{
df_i = df_total %>%
as_tibble() %>% x
filter(id==id_problematico) %>%
arrange(Estado, desc(TPVEstimate), desc(porte)) ## Prioritizing entries with a valid state variable
## and with the highest TPV estimate and size
df_i = slice_head(df_i, n=n_distinct(df_i$mes_referencia))
linhas_certas = bind_rows(linhas_certas, df_i)
print(contador)
contador = contador + 1
}
## Eliminating wrong id's from dataframe
df_correct = df_total %>%
as_tibble() %>%
filter(!(id %in% ids_errados$id))
df_correct = bind_rows(df_correct, linhas_certas) %>%
arrange(id, mes_referencia)O último passo é apenas eliminar os ids errados selecionados e salvar a base de dados:
## Eliminating wrong id's from dataframe
df_correct = df_total %>%
as_tibble() %>%
filter(!(id %in% ids_errados$id))
df_correct = bind_rows(df_correct, linhas_certas) %>%
arrange(id, mes_referencia)
#### Saving cleaned dataframe
write_csv(df_correct, "Data/Clean/stone_data_cleaned.csv")Dados adicionais do ipeadata
Decidi trazer algumas variáveis (a nível de estado e de país) que acredito úteis para medir o nível de atividade econômica e de preços. São elas:
- Nível país:
- Produção total de caixas de papelão no mês
- Variação da produção de caixas de papelão no mês
- IPCA mensal
- Taxa de câmbio mensal
- Nível estado:
- Número de demissões no mês (dados do CAGED)
- Variação da produção de caixas de papelão
- Índice do varejo no estado
Todas essas variáveis estavam bem formatadas, motivo pelo qual decidi omitir o código que une tais variáveis com as da Stone.
Análise Exploratória
A partir de agora, buscarei conduzir uma análise explorátória da base de dados unificada, com dados da Stone e do ipeadata. Nesta etapa, focarei meus comentários na análise dos resultados e nos trechos de código mais importantes. Omitirei os códigos para os gráficos porque eles ocupariam um espaço muito grande.
TPV Mensal
Primeiramente, podemos carregar as bases e, também, a fonte que decidi usar nos gráficos:
Começarei pelo TPV Mensal, nossa variável explicativa e, portanto, a mais importante. Vamos observar sua distribuição: 
É possível observar que a distribuição dos TPVs mensais possui uma grande assimetria positiva, indicando que a maior parte das empresas e clientes Stone concentram-se em faixas menores de faturamento.
Vamos investigar mais um pouco a distribuição do TPV, ver como ela se distribui pelos estados e, agrupar a base por estado para calcular a média do TPV de cada um:
### Number of clients by state
df %>% group_by(Estado) %>%
summarise(TPV_mensal = mean(TPV_mensal)) %>%
head(n=10)## # A tibble: 10 x 2
## Estado TPV_mensal
## <chr> <dbl>
## 1 AC 13778.
## 2 AL 21871.
## 3 AM 16293.
## 4 AP 21823.
## 5 BA 19383.
## 6 CE 22902.
## 7 DF 27776.
## 8 ES 16858.
## 9 GO 19706.
## 10 MA 17032.Agora já podemos desenhar o gráfico de barras: 
Podemos ver que Roraima é o estado com a maior média de TPV mensal. Isso pode parecer surpreendente, a princípio, mas devemos nos atentar também à quantidade de clientes nesse estado para verificar se isso não é efeito de algum outlier. Do outro lado, o estado com menor faturamento médio é o Acre. Apesar desses valores não esperados, vemos uma tendência geral dos estados mais ricos terem uma média de faturamento maior.
Agora, vamos investigar o número de clientes por estado para confirmar a suspeita sobre Roraima:
### Number of clients by state
states = df %>%
as_tibble() %>%
count(Estado) %>%
drop_na() %>%
rename("n_clientes"=n) %>%
arrange(n_clientes)
head(states, n=5)## # A tibble: 5 x 2
## Estado n_clientes
## <chr> <int>
## 1 AC 7437
## 2 AP 7654
## 3 RR 11167
## 4 PI 11348
## 5 TO 12087
Podemos ver que nossa suspeita sobre Roraima se confirmou. Muito provavelmente existe algum outlier com faturamento muito alto em Roraima. Além disso, vemos que São Paulo tem mais que o dobro de clientes que o segundo colocado (Rio de Janeiro). Vemos que a Stone tem uma presença maior na região Sudeste, seguida pela região Sul.
Agora, vamos investigar como varia a média do TPV para cada Macro Classificação de setores:
Antes de fazer o gráfico, precisamos agrupar o dataframe por Macro Classificação:
## Agrupando
df_est = df %>% group_by(MacroClassificacao) %>%
summarise(TPV_mensal = mean(TPV_mensal)) %>%
drop_na() %>%
arrange(TPV_mensal)
## Ordem
limits_mac = df_est$MacroClassificacao
Podemos ver que, dentre as áreas de Macro Classificação, os postos de gasolina apresentam um faturameto muito superior em relação aos demais grupos.
Vamos investigar um pouco mais a fundo a média do TPV entre diversos setores. Dessa vez vamos agrupar a base pela variável segmento:
## Agrupando
df_seg = df %>% group_by(segmento) %>%
summarise(TPV_mensal = mean(TPV_mensal)) %>%
drop_na() %>%
arrange(TPV_mensal)
## Ordem
limits_seg = df_seg$segmento
Podemos ver que os segmentos de logística e de postos de gasolina são os com maior média de faturamento. Serviços de beleza e companhias aparecem por último. Apesar disso, o valor de companhias aéreas parece ser anormalmente baixo.
Variáveis do ipeadata
Agora, vamos analisar um pouco as variáveis do ipeadata. A melhor forma de começar é construindo um heatmap das correlações entre elas e o TPV mensal:
Observando-se a matriz, vemos que as correlações entre as variáveis do ipeadata e o TPV mensal não são muito altas. Apesar disso, vemos uma correlação ligeiramente positiva com o ipca e o índice de varejo, e uma correlação lugeiramente negativa com o dólar. Esses efeitos estão dentro do esperado, uma vez que o nível de câmbio e da inflação tem um efeito direto no faturamento. Devemos nos manter atentos porque há diversos segmentos de clientes que podem ter diferentes interações com essas variáveis, que acabam mascaradas quando agregadas.
Ainda, quero ressaltar que será adicionado lag temporal e médias móveis a todas essas variáveis posteriormente. Uma vez que esperamos uma forte característica autoregressiva no modelo.
Como acabei de comentar, vamos ver o impacto dessas variáveis no TPV de forma mais detalhada. Primeiramente, vamos começar vendo o impacto do nível do ipca de acordo com as Macro Classificações de clientes: 
Observando os gráficos, podemos ver uma ligeira tendência positiva (embora pequena) entre inflação(IPCA) e o TPV. Além disso, vemos que o faturamento de postos de gasolina é o mais sensível à inflação. Provavelmente, por ser um setor monopolizado, o ajuste dos preços ocorre mais rapidamente.
Vamos, agora, construir o mesmo gráfico para a taxa de câmbio e ver se nossas expectativas estavam corretas: 
Estávamos corretos! É possível ver uma clara relação negativa entre o preço do dólar em reais e o TPV. O setor de postos tem seu faturamento muito prejudicado com um dólar mais apreciado. Por outro lado, vemos que o setor de supermercados e farmácias tem um aumento em seu faturamento (provavelmente por serem setores essenciais). Varejo e entretenimento também são prejudicados pelo aumento do dólar.
O que aprendemos com isso? Que analisar a correlação, com todos os setores agrupados, esconde o poder preditivo de muitas variáveis. Dessa forma, nosso “feeling econômico” de que essas variáveis eram importantes estava correto.
Vamos observar, agora, a relação entre o TPV e o índice de nível do varejo. É de se esperar que tal relação seja mais forte para o macro setor de varejo. Vamos ver:
Mais uma vez, vemos que nossa intuição econômica estava correta. O índice de varejo aparenta ser uma boa proxy para a performance econômica deste setor e de outros e, certamente, contribuirá com o poder preditivo do modelo.
Por último, vamos visualizar a relação entre o número de demissões contabilizado pelo CAGED e o TPV Mensal: 
Neste caso, a relação não fica tão clara ou interpretável. Vemos que o faturamento de postos de gasolina claramente aumenta em meses com mais demissões, algo muito interessante.
Para finalizar, gostaria de dizer, mais uma vez, que todas as variáveis que trouxe do ipea podem ter uma importância especial para um grupo específico de clientes. Dessa forma, é válido deixá-las no modelo. Não adicionei mais variáveis pois pretendo adicionar vários lags temporais e médias móveis para todas as variáveis listadas. Dessa forma, espero que o número de variáveis regressoras aumente de forma substancial na etapa de feature engineering.
Feature Engineering
Minha estratégia de feature engineering foi bastante convencional. Meu método foi aplicar uma série de lags e médias móveis à todas as variáveis, a fim de evitar data leakage entre os períodos.
Vale comentar que a partir desse momento continuei o desafio com Python.
Primeiramente, podemos importar os pacotes e os dataframes utilizados no script. Nesse markdown utilizarei apenas as primeiras 1000 linhas do dataframe por ser somente uma demonstração dos métodos aplicados:
from stone_data_challenge.config import data_dir
from sklearn.impute import SimpleImputer
import datetime as dt
import numpy as np
import datetime
import pandas as pd
# # Reading dataframe
pd.set_option('display.max_columns', 10)
pd.set_option('display.width', 1000)
df = pd.read_csv(f"{data_dir}/data/clean/spine_stone_ipea.csv")
df = df.head(100000)
df = df.replace("nan", np.NaN)
df = df.drop_duplicates()
df.head(5)## id StoneCreatedDate StoneFirstTransactionDate MCC MacroClassificacao ... caixas_papelao_varpct ipca_pct cambio_dol demissoes ind_varejo
## 0 1 2018-12-12T00:00:00Z 2018-12-13 5942 Varejo ... -14.092461 0.15 3.882136 142853.0 149.6
## 1 1 2018-12-12T00:00:00Z 2018-12-13 5942 Varejo ... 8.555862 0.32 3.741055 140224.0 123.5
## 2 1 2018-12-12T00:00:00Z 2018-12-13 5942 Varejo ... -3.899083 0.43 3.723025 138702.0 114.4
## 3 1 2018-12-12T00:00:00Z 2018-12-13 5942 Varejo ... 4.087390 0.75 3.845884 139766.0 125.0
## 4 1 2018-12-12T00:00:00Z 2018-12-13 5942 Varejo ... 1.185331 0.57 3.895557 140085.0 122.3
##
## [5 rows x 21 columns]
Primeiramente, vamos checar os tipos das variáveis para ver se podemos armazená-las em algum tipo de dado mais eficiente:
df.dtypes## id int64
## StoneCreatedDate object
## StoneFirstTransactionDate object
## MCC int64
## MacroClassificacao object
## segmento object
## sub_segmento object
## persona object
## porte int64
## TPVEstimate float64
## tipo_documento object
## Estado object
## mes_referencia object
## TPV_mensal float64
## mes object
## caixas_papelao int64
## caixas_papelao_varpct float64
## ipca_pct float64
## cambio_dol float64
## demissoes float64
## ind_varejo float64
## dtype: object
Vamos agora corrigir os tipos de algumas variáveis:
# # Fixing types
df[["MCC", "MacroClassificacao", "StoneCreatedDate", "StoneFirstTransactionDate", "segmento",
"sub_segmento", "persona", "tipo_documento", "Estado", "mes_referencia", "id"]] = df[["MCC", "MacroClassificacao", "StoneCreatedDate", "StoneFirstTransactionDate", "segmento", "sub_segmento", "persona", "tipo_documento", "Estado", "mes_referencia", "id"]].astype(str)
df.dtypes## id object
## StoneCreatedDate object
## StoneFirstTransactionDate object
## MCC object
## MacroClassificacao object
## segmento object
## sub_segmento object
## persona object
## porte int64
## TPVEstimate float64
## tipo_documento object
## Estado object
## mes_referencia object
## TPV_mensal float64
## mes object
## caixas_papelao int64
## caixas_papelao_varpct float64
## ipca_pct float64
## cambio_dol float64
## demissoes float64
## ind_varejo float64
## dtype: object
Agora podemos consertar o formato das datas:
# # Fixing date formatting
df["StoneCreatedDate"] = pd.to_datetime(df["StoneCreatedDate"]).dt.tz_localize(None)
df["StoneFirstTransactionDate"] = pd.to_datetime(df["StoneFirstTransactionDate"])
df["mes_referencia"] = pd.to_datetime(df["mes_referencia"])
df[["StoneFirstTransactionDate", "StoneCreatedDate"]].head(5)## StoneFirstTransactionDate StoneCreatedDate
## 0 2018-12-13 2018-12-12
## 1 2018-12-13 2018-12-12
## 2 2018-12-13 2018-12-12
## 3 2018-12-13 2018-12-12
## 4 2018-12-13 2018-12-12
Algo que podemos fazer para o modelo é transformar as variáveis sobre a data de cadastro do cliente e sua primeira transação em números ordinais. Observe:
# # Transforming date variables in cardinal numbers
initial_date = dt.datetime(2020, 7, 31)
df["DaysSinceCreation"] = -(df["StoneCreatedDate"] - initial_date).dt.days
df["DaysSinceFirstTrans"] = -(df["StoneFirstTransactionDate"] - initial_date).dt.days
df[["DaysSinceCreation", "DaysSinceFirstTrans"]]## DaysSinceCreation DaysSinceFirstTrans
## 0 597 596.0
## 1 597 596.0
## 2 597 596.0
## 3 597 596.0
## 4 597 596.0
## ... ... ...
## 99990 242 241.0
## 99992 242 241.0
## 99994 242 241.0
## 99996 242 241.0
## 99998 242 241.0
##
## [77817 rows x 2 columns]
Para que nosso modelo possa identificar os efeitos de sazonalidade, é interessante dividir a coluna “mes_referencia” em duas variáveis diferentes: uma para mês e outra para ano. Apliquei o mesmo tratamento para “StoneCreatedDate” e “StoneFirstTransactionDate”.
# # Creating separate year and month variables for date imputing
df["YearCreated"] = df["StoneCreatedDate"].dt.year
df["MonthCreated"] = df["StoneCreatedDate"].dt.month
df["YearFirstTrans"] = df["StoneFirstTransactionDate"].dt.year
df["MonthFirstTrans"] = df["StoneFirstTransactionDate"].dt.month
df["AnoReferencia"] = df["mes_referencia"].dt.year
df["MesReferencia"] = df["mes_referencia"].dt.month
df[["YearCreated", "MonthCreated", "YearFirstTrans", "MonthFirstTrans", "AnoReferencia", "MesReferencia"]]## YearCreated MonthCreated YearFirstTrans MonthFirstTrans AnoReferencia MesReferencia
## 0 2018 12 2018.0 12.0 2018 12
## 1 2018 12 2018.0 12.0 2019 1
## 2 2018 12 2018.0 12.0 2019 2
## 3 2018 12 2018.0 12.0 2019 3
## 4 2018 12 2018.0 12.0 2019 4
## ... ... ... ... ... ... ...
## 99990 2019 12 2019.0 12.0 2019 6
## 99992 2019 12 2019.0 12.0 2019 7
## 99994 2019 12 2019.0 12.0 2019 8
## 99996 2019 12 2019.0 12.0 2019 9
## 99998 2019 12 2019.0 12.0 2019 10
##
## [77817 rows x 6 columns]
Criadas essas variáveis, agora podemos remover as variáveis de data não ordinais, com excessão de “mes_referencia” que será utilizada como índice mais adiante:
# # # Removing unwanted variables
df.drop(["mes", "StoneFirstTransactionDate", "StoneCreatedDate"], axis="columns", inplace=True)
Vamos checar a porcentagem de valores faltantes em cada coluna no dataframe:
# # Checking for missing data
df.replace('nan', np.nan, inplace=True)
percent_missing = df.isna().sum() * 100 / len(df)
missing_value_df = pd.DataFrame({"col_name": df.columns, "pct_missing": percent_missing})
missing_value_df## col_name pct_missing
## id id 0.000000
## MCC MCC 0.000000
## MacroClassificacao MacroClassificacao 0.109231
## segmento segmento 0.109231
## sub_segmento sub_segmento 0.109231
## persona persona 0.000000
## porte porte 0.000000
## TPVEstimate TPVEstimate 0.000000
## tipo_documento tipo_documento 0.000000
## Estado Estado 1.453410
## mes_referencia mes_referencia 0.000000
## TPV_mensal TPV_mensal 0.000000
## caixas_papelao caixas_papelao 0.000000
## caixas_papelao_varpct caixas_papelao_varpct 0.000000
## ipca_pct ipca_pct 0.000000
## cambio_dol cambio_dol 0.000000
## demissoes demissoes 1.453410
## ind_varejo ind_varejo 1.453410
## DaysSinceCreation DaysSinceCreation 0.000000
## DaysSinceFirstTrans DaysSinceFirstTrans 1.904468
## YearCreated YearCreated 0.000000
## MonthCreated MonthCreated 0.000000
## YearFirstTrans YearFirstTrans 1.904468
## MonthFirstTrans MonthFirstTrans 1.904468
## AnoReferencia AnoReferencia 0.000000
## MesReferencia MesReferencia 0.000000
Primeiramente, podemos imputar alguns valores categóricos como “Desconhecido” para as colunas referentes ao setor dos clientes:
# # Checking for missing data
df["MacroClassificacao"] = df["MacroClassificacao"].fillna("Desconhecido")
df["segmento"] = df["segmento"].fillna("Desconhecido")
df["sub_segmento"] = df["sub_segmento"].fillna("Desconhecido")
Da análise exploratória, sabemos que São Paulo é, de longe, o estado com maior número de clientes Stone. Dessa forma podemos imputar os valores faltantes com São Paulo:
# # Imputing missing categorical data for state
df["Estado"] = df["Estado"].fillna("SP")
Por último, sobraram os valores faltantes para as colunas “ind_varejo”, “demissoes”, “DaysSinceFirstTrans”, “YearFirstTrans”e “MonthFirstTrans”. Para imputar os valores nulos, utilizarei o imputador do scikit-learn e substituirei os valores faltantes pela mediana:
# Creating column transformer
missing_vars = ["ind_varejo", "demissoes", "DaysSinceFirstTrans", "YearFirstTrans", "MonthFirstTrans"]
imputer = SimpleImputer(strategy="median")
df[missing_vars] = imputer.fit_transform(df[missing_vars]) # Imputing data
Agora, podemos conferir que não há mais valores nulos:
# # Checking for missing data
percent_missing = df.isnull().sum() * 100 / len(df)
missing_value_df = pd.DataFrame({"col_name": df.columns, "pct_missing": percent_missing})
missing_value_df## col_name pct_missing
## id id 0.0
## MCC MCC 0.0
## MacroClassificacao MacroClassificacao 0.0
## segmento segmento 0.0
## sub_segmento sub_segmento 0.0
## persona persona 0.0
## porte porte 0.0
## TPVEstimate TPVEstimate 0.0
## tipo_documento tipo_documento 0.0
## Estado Estado 0.0
## mes_referencia mes_referencia 0.0
## TPV_mensal TPV_mensal 0.0
## caixas_papelao caixas_papelao 0.0
## caixas_papelao_varpct caixas_papelao_varpct 0.0
## ipca_pct ipca_pct 0.0
## cambio_dol cambio_dol 0.0
## demissoes demissoes 0.0
## ind_varejo ind_varejo 0.0
## DaysSinceCreation DaysSinceCreation 0.0
## DaysSinceFirstTrans DaysSinceFirstTrans 0.0
## YearCreated YearCreated 0.0
## MonthCreated MonthCreated 0.0
## YearFirstTrans YearFirstTrans 0.0
## MonthFirstTrans MonthFirstTrans 0.0
## AnoReferencia AnoReferencia 0.0
## MesReferencia MesReferencia 0.0Chegamos em um momento muito importante. Pela forte característica autoregressiva de nosso problema, devemos gerar uma série de variáveis com defasagem temporal de 1 e mais períodos. Também é interessante gerar médias móveis para distintas janelas de tempo.
Antes de começar a adicionar variáveis, devemos adicionar 5 meses à base de dados, para que possamos utilizar essas variáveis com lag e médias móveis na previsão final do desafio. Esse é o código que utilizei para isso:
df = df.set_index(["mes_referencia"])
# # We need to add 5 months to the dataset before shifting the columns
df_agosto = df.groupby("id").tail(1)
df_agosto.index = df_agosto.index + datetime.timedelta(days=30)
df_agosto.loc[:, "TPV_mensal"] = np.nan## /Users/maxmitteldorf/Library/r-miniconda/envs/r-reticulate/lib/python3.6/site-packages/pandas/core/indexing.py:1763: SettingWithCopyWarning:
## A value is trying to be set on a copy of a slice from a DataFrame.
## Try using .loc[row_indexer,col_indexer] = value instead
##
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
## isetter(loc, value)df_agosto.loc[:, "caixas_papelao"] = np.nan
df_agosto.loc[:, "caixas_papelao_varpct"] = np.nan
df_agosto.loc[:, "ipca_pct"] = np.nan
df_agosto.loc[:, "cambio_dol"] = np.nan
df_agosto.loc[:, "demissoes"] = np.nan
df_agosto.loc[:, "ind_varejo"] = np.nan
df_agosto[["DaysSinceCreation", "DaysSinceFirstTrans"]] = df_agosto[["DaysSinceCreation", "DaysSinceFirstTrans"]] + 30## /Users/maxmitteldorf/Library/r-miniconda/envs/r-reticulate/lib/python3.6/site-packages/pandas/core/frame.py:3069: SettingWithCopyWarning:
## A value is trying to be set on a copy of a slice from a DataFrame.
## Try using .loc[row_indexer,col_indexer] = value instead
##
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
## self[k1] = value[k2]df_agosto["MesReferencia"] = df_agosto["MesReferencia"] + 1## /Users/maxmitteldorf/Library/r-miniconda/envs/r-reticulate/bin/python:1: SettingWithCopyWarning:
## A value is trying to be set on a copy of a slice from a DataFrame.
## Try using .loc[row_indexer,col_indexer] = value instead
##
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copydf_setembro = df_agosto[:]
df_setembro.index = df_setembro.index + datetime.timedelta(days=30)
df_setembro[["DaysSinceCreation", "DaysSinceFirstTrans"]] = df_setembro[["DaysSinceCreation", "DaysSinceFirstTrans"]] + 30
df_setembro["MesReferencia"] = df_setembro["MesReferencia"] + 1
df_outubro = df_setembro[:]
df_outubro.index = df_outubro.index + datetime.timedelta(days=30)
df_outubro[["DaysSinceCreation", "DaysSinceFirstTrans"]] = df_outubro[["DaysSinceCreation", "DaysSinceFirstTrans"]] + 30
df_outubro["MesReferencia"] = df_outubro["MesReferencia"] + 1
df_novembro = df_outubro[:]
df_novembro.index = df_novembro.index + datetime.timedelta(days=30)
df_novembro[["DaysSinceCreation", "DaysSinceFirstTrans"]] = df_novembro[["DaysSinceCreation", "DaysSinceFirstTrans"]] + 30
df_novembro["MesReferencia"] = df_novembro["MesReferencia"] + 1
df_dezembro = df_novembro[:]
df_dezembro.index = df_dezembro.index + datetime.timedelta(days=30)
df_dezembro[["DaysSinceCreation", "DaysSinceFirstTrans"]] = df_dezembro[["DaysSinceCreation", "DaysSinceFirstTrans"]] + 30
df_dezembro["MesReferencia"] = df_dezembro["MesReferencia"] + 1
df = pd.concat([df, df_agosto,df_setembro,df_outubro,df_novembro,df_dezembro], axis=0)
df = df.sort_values(["id", "mes_referencia"])
Agora, podemos conferir que temos mais 5 meses na base de dados:
df.tail(8)## id MCC MacroClassificacao segmento sub_segmento ... MonthCreated YearFirstTrans MonthFirstTrans AnoReferencia MesReferencia
## mes_referencia ...
## 2020-05-31 999 5411 Supermercado/Farmácia Supermercados Supermercados ... 3 2020.0 3.0 2020 5
## 2020-06-30 999 5411 Supermercado/Farmácia Supermercados Supermercados ... 3 2020.0 3.0 2020 6
## 2020-07-31 999 5411 Supermercado/Farmácia Supermercados Supermercados ... 3 2020.0 3.0 2020 7
## 2020-08-30 999 5411 Supermercado/Farmácia Supermercados Supermercados ... 3 2020.0 3.0 2020 12
## 2020-09-29 999 5411 Supermercado/Farmácia Supermercados Supermercados ... 3 2020.0 3.0 2020 12
## 2020-10-29 999 5411 Supermercado/Farmácia Supermercados Supermercados ... 3 2020.0 3.0 2020 12
## 2020-11-28 999 5411 Supermercado/Farmácia Supermercados Supermercados ... 3 2020.0 3.0 2020 12
## 2020-12-28 999 5411 Supermercado/Farmácia Supermercados Supermercados ... 3 2020.0 3.0 2020 12
##
## [8 rows x 25 columns]
Vamos começar a adicionar as médias móveis e variáveis com lag. Vou exemplificar o tratamento que conferi ao TPV mensal. Repeti exatamente o mesmo procedimento para as outras variáveis:
# Lag TPV
df["TPV_lag_1"] = df.groupby("id")["TPV_mensal"].shift(1)
df["TPV_lag_2"] = df.groupby("id")["TPV_mensal"].shift(2)
df["TPV_lag_3"] = df.groupby("id")["TPV_mensal"].shift(3)
df["TPV_lag_4"] = df.groupby("id")["TPV_mensal"].shift(4)
df["TPV_lag_5"] = df.groupby("id")["TPV_mensal"].shift(5)
df["TPV_lag_6"] = df.groupby("id")["TPV_mensal"].shift(6)
df["TPV_lag_12"] = df.groupby("id")["TPV_mensal"].shift(12)
# Creating moving averages for TPV
df["TPV_MA_2"] = df.groupby("id")["TPV_mensal"].transform(lambda x: x.rolling(2, 2).mean())
df["TPV_MA_3"] = df.groupby("id")["TPV_mensal"].transform(lambda x: x.rolling(3, 2).mean())
df["TPV_MA_4"] = df.groupby("id")["TPV_mensal"].transform(lambda x: x.rolling(4, 2).mean())
df["TPV_MA_5"] = df.groupby("id")["TPV_mensal"].transform(lambda x: x.rolling(5, 2).mean())
df["TPV_MA_6"] = df.groupby("id")["TPV_mensal"].transform(lambda x: x.rolling(6, 2).mean())
df["TPV_MA_12"] = df.groupby("id")["TPV_mensal"].transform(lambda x: x.rolling(12, 2).mean())
df["TPV_MA_18"] = df.groupby("id")["TPV_mensal"].transform(lambda x: x.rolling(18, 2).mean())
Vamos verificar que essas transformações funcionaram:
df[["TPV_mensal","TPV_lag_1","TPV_lag_3","TPV_lag_3","TPV_lag_4","TPV_MA_2","TPV_MA_3","TPV_MA_4"]].tail(10)## TPV_mensal TPV_lag_1 TPV_lag_3 TPV_lag_3 TPV_lag_4 TPV_MA_2 TPV_MA_3 TPV_MA_4
## mes_referencia
## 2020-03-31 3585.26 NaN NaN NaN NaN NaN NaN NaN
## 2020-04-30 6040.58 3585.26 NaN NaN NaN 4812.920 4812.920000 4812.920000
## 2020-05-31 6224.67 6040.58 NaN NaN NaN 6132.625 5283.503333 5283.503333
## 2020-06-30 9100.48 6224.67 3585.26 3585.26 NaN 7662.575 7121.910000 6237.747500
## 2020-07-31 11630.05 9100.48 6040.58 6040.58 3585.26 10365.265 8985.066667 8248.945000
## 2020-08-30 NaN 11630.05 6224.67 6224.67 6040.58 NaN 10365.265000 8985.066667
## 2020-09-29 NaN NaN 9100.48 9100.48 6224.67 NaN NaN 10365.265000
## 2020-10-29 NaN NaN 11630.05 11630.05 9100.48 NaN NaN NaN
## 2020-11-28 NaN NaN NaN NaN 11630.05 NaN NaN NaN
## 2020-12-28 NaN NaN NaN NaN NaN NaN NaN NaN
Por último, após aplicar tal procedimento em todas as variáveis, podemos salvar o dataframe final para a modelagem:
# Writing dataframe
df = df.reset_index()
df.to_csv(f"{data_dir}/feat_eng/model_data_v5.csv", index=False)Modelagem
Diante de um problema com um componente temporal, e que exige a previsão de 5 lags para frente, optei por construir 5 modelos independentes (um para cada lag temporal).
Decidi por utilizar o LightGBM para fazer as previsões por sua rapidez e robustez. Diante de um grande número de hiperparâmetros, também optei por utilizar o pacote Optuna para otimizá-los.
Mostrarei o código para o modelo que prevê um lag temporal (t+1). O código para os modelos (t+2) em diante são muito similares e estão disponíveis no repositório.
Primeiramente, vamos importar os pacotes necessários:
import pandas as pd
import numpy as np
from sklearn.model_selection._split import _BaseKFold, indexable, _num_samples
from sklearn.utils.validation import _deprecate_positional_args
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import mean_absolute_error
import lightgbm as lgb
import optuna
import pickle
Em seguida, as configurações do script devem ser definidas: * TUNING - Indica se a otimização dos parâmetros será realizada * TUNING_SIZE - Indica o tamanho do dataframe para a otimização dos parâmetros * TRAINING - Indica se o modelo será treinado ou não * TRAINING_SIZE - Indica o tamanho do dataframe para o treino do modelo * FINNAL_PREDICTION - Indica se queremos treinar o modelo em todos os dados e exportar a previsão final do desafio
Agora, podemos começar lendo o dataframe e removendo as colunas sem lag temporal para evitarmos problemas de data leakage:
# # Reading dataframe
print("Reading Dataframe")
df = pd.read_csv(f"{data_dir}/data/feat_eng/model_data_v5.csv")
print("Dataframe succesfully readed")
# # Dropping columns without lag e MA2
df.drop(["caixas_papelao", "caixas_papelao_varpct", "ipca_pct", "cambio_dol", "demissoes", "ind_varejo",
"TPV_MA_2", "ipca_MA_2", "caixas_papelao_MA_2", "caixas_papelao_varpct_MA_2",
"demissoes_MA_2", "ind_varejo_MA_2", "cambio_dol_MA_2"], axis="columns", inplace=True)
# # Renaming TPV columns
df.rename(columns={"TPV_mensal": "TPV_t1"}, inplace=True)
Em seguida, usei o LabelEncoder() para transformar as variáveis discretas em numéricas. Utilizei o LabelEncoder() no lugar do OneHotEcoder() porque modelos de árvore como o LightGBM lidam bem com tal formato de dados.
# # Using LabelEncoder for categorical variables
# # Lightgbm can handle label encoded variables very efficiently
le = LabelEncoder()
df["MCC"] = le.fit_transform(df["MCC"])
df["MacroClassificacao"] = le.fit_transform(df["MacroClassificacao"])
df["segmento"] = le.fit_transform(df["segmento"])
df["sub_segmento"] = le.fit_transform(df["sub_segmento"])
df["persona"] = le.fit_transform(df["persona"])
df["Estado"] = le.fit_transform(df["Estado"])
df["tipo_documento"] = le.fit_transform(df["tipo_documento"])
Sendo assim, também podemos separar os meses da previsão final (agosto até dezembro) da base de dados e separá-los para um momento posterior:
# # Saving dataset for prediction
df["mes_referencia"] = pd.to_datetime(df["mes_referencia"])
df_prediction = df.groupby("id").tail(5).set_index("mes_referencia")
df_prediction_august = df_prediction.groupby("id").head(1).drop("TPV_t1", axis="columns")
df_prediction_august["ind"] = df_prediction_august["id"]
df_prediction_august = df_prediction_august.set_index("ind")
df = df[df["mes_referencia"] < "2020-08-30"]
Agora, como desejamos poder avaliar as métricas de performance do modelo, devemos remover os próximos 5 meses de dados para deixá-los como base de validação final. Também tive que fazer algumas manipulações para remover os dados de algumas variáveis regressoras que não deveriam estar disponíveis no momento da execução do modelo.
df_val = df.groupby("id").tail(5).groupby("id").head(1)
df_val = df_val[pd.notnull(df_val["TPV_t1"])].set_index("mes_referencia")
df_train = df.drop(df.groupby("id").tail(5).index, axis=0).set_index("mes_referencia")
df_train = df_train[pd.notnull(df_train["TPV_t1"])]
df_train = df_train.head(TRAINING_SIZE)
df = df.set_index("mes_referencia")
df = df[pd.notnull(df["TPV_t1"])]
# # Saving column names for later
colnames = df.columns
# # Saving categorical vars names
x_colnames = list(df.drop("TPV_t1", axis="columns").columns.values)
categorical_vars = ["id", "MCC", "MacroClassificacao", "segmento", "sub_segmento", "persona", "tipo_documento", "Estado"]
Agora, gostaria de falar um pouco sobre minha estratégia para separar as bases de treino e teste. Estamos diante de um problema de séries temporais, por isso, a princípio, poderíamos utilizar o TimeSerisSplit() para treinar nosso modelo em vários ranges de tempo sem que ocorra time leakage. A questão é que nosso problema é de séries temporais com grupos (nesse caso o id dos clientes). Dessa forma, o TimeSeriesSplit() poderia misturar clientes com ids diferentes, o que não queremos de jeito alguma.
É aí que entra o GroupTimeSeriesSplit(), uma estratégia de divisão que respeita tanto a temporalidade, quanto o grupo das observações. Essa funcionalidade ainda não está completamente implementada no scikit-learn, mas está em vias de ser. Dessa forma, conseguimos acessar o código para ela nesse link.
class GroupTimeSeriesSplit(_BaseKFold):
"""Time Series cross-validator variant with non-overlapping groups.
Provides train/test indices to split time series data samples
that are observed at fixed time intervals according to a
third-party provided group.
In each split, test indices must be higher than before, and thus shuffling
in cross validator is inappropriate.
This cross-validation object is a variation of :class:`KFold`.
In the kth split, it returns first k folds as train set and the
(k+1)th fold as test set.
The same group will not appear in two different folds (the number of
distinct groups has to be at least equal to the number of folds).
Note that unlike standard cross-validation methods, successive
training sets are supersets of those that come before them.
Read more in the :ref:`User Guide <cross_validation>`.
Parameters
----------
n_splits : int, default=5
Number of splits. Must be at least 2.
max_train_size : int, default=None
Maximum size for a single training set.
Examples
--------
>>> import numpy as np
>>> from sklearn.model_selection import GroupTimeSeriesSplit
>>> groups = np.array(['a', 'a', 'a', 'a', 'a', 'a',\
'b', 'b', 'b', 'b', 'b',\
'c', 'c', 'c', 'c',\
'd', 'd', 'd'])
>>> gtss = GroupTimeSeriesSplit(n_splits=3)
>>> for train_idx, test_idx in gtss.split(groups, groups=groups):
... print("TRAIN:", train_idx, "TEST:", test_idx)
... print("TRAIN GROUP:", groups[train_idx],\
"TEST GROUP:", groups[test_idx])
TRAIN: [0, 1, 2, 3, 4, 5] TEST: [6, 7, 8, 9, 10]
TRAIN GROUP: ['a' 'a' 'a' 'a' 'a' 'a']\
TEST GROUP: ['b' 'b' 'b' 'b' 'b']
TRAIN: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] TEST: [11, 12, 13, 14]
TRAIN GROUP: ['a' 'a' 'a' 'a' 'a' 'a' 'b' 'b' 'b' 'b' 'b']\
TEST GROUP: ['c' 'c' 'c' 'c']
TRAIN: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]\
TEST: [15, 16, 17]
TRAIN GROUP: ['a' 'a' 'a' 'a' 'a' 'a' 'b' 'b' 'b' 'b' 'b' 'c' 'c' 'c' 'c']\
TEST GROUP: ['d' 'd' 'd']
"""
@_deprecate_positional_args
def __init__(self,
n_splits=5,
*,
max_train_size=None
):
super().__init__(n_splits, shuffle=False, random_state=None)
self.max_train_size = max_train_size
def split(self, X, y=None, groups=None):
"""Generate indices to split data into training and test set.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training data, where n_samples is the number of samples
and n_features is the number of features.
y : array-like of shape (n_samples,)
Always ignored, exists for compatibility.
groups : array-like of shape (n_samples,)
Group labels for the samples used while splitting the dataset into
train/test set.
Yields
------
train : ndarray
The training set indices for that split.
test : ndarray
The testing set indices for that split.
"""
if groups is None:
raise ValueError(
"The 'groups' parameter should not be None")
X, y, groups = indexable(X, y, groups)
n_samples = _num_samples(X)
n_splits = self.n_splits
n_folds = n_splits + 1
group_dict = {}
u, ind = np.unique(groups, return_index=True)
unique_groups = u[np.argsort(ind)]
n_samples = _num_samples(X)
n_groups = _num_samples(unique_groups)
for idx in np.arange(n_samples):
if (groups[idx] in group_dict):
group_dict[groups[idx]].append(idx)
else:
group_dict[groups[idx]] = [idx]
if n_folds > n_groups:
raise ValueError(
("Cannot have number of folds={0} greater than"
" the number of groups={1}").format(n_folds,
n_groups))
group_test_size = n_groups // n_folds
group_test_starts = range(n_groups - n_splits * group_test_size,
n_groups, group_test_size)
for group_test_start in group_test_starts:
train_array = []
test_array = []
for train_group_idx in unique_groups[:group_test_start]:
train_array_tmp = group_dict[train_group_idx]
train_array = np.sort(np.unique(
np.concatenate((train_array,
train_array_tmp)),
axis=None), axis=None)
train_end = train_array.size
if self.max_train_size and self.max_train_size < train_end:
train_array = train_array[train_end -
self.max_train_size:train_end]
for test_group_idx in unique_groups[group_test_start:
group_test_start +
group_test_size]:
test_array_tmp = group_dict[test_group_idx]
test_array = np.sort(np.unique(
np.concatenate((test_array,
test_array_tmp)),
axis=None), axis=None)
yield [int(i) for i in train_array], [int(i) for i in test_array]
A partir daqui, podemos definir uma função de otimização para o Optuna e realizar o tuning dos parâmetros do LightGBM em uma porção do dataframe. Após a otimização, os melhores parâmetros são salvados em um dicionário através de um pickle.
if TUNING:
df_train_tuning = df_train.head(TUNING_SIZE)
# # Objective Function
def objective(trial, cv_fold_func=np.average):
params = {
'objective': 'regression',
'metric': 'mae',
'boosting': 'rf',
'lambda_l1': trial.suggest_loguniform('lambda_l1', 1e-8, 10.0),
'lambda_l2': trial.suggest_loguniform('lambda_l2', 1e-8, 10.0),
'num_leaves': trial.suggest_int('num_leaves', 2, 200),
'feature_fraction': trial.suggest_uniform('feature_fraction', 0.4, 1.0),
'bagging_fraction': trial.suggest_uniform('bagging_fraction', 0.4, 1.0),
'bagging_freq': trial.suggest_int('bagging_freq', 1, 7),
'min_child_samples': trial.suggest_int('min_child_samples', 5, 100),
'max_depth': trial.suggest_int('max_depth', 5, 15),
'learning_rate': trial.suggest_uniform('learning_rate', 0.01, 0.10)}
# fitting and returning RMSE scores
mae_list = []
for train_idx, test_idx in GroupTimeSeriesSplit().split(df_train_tuning, groups=df_train_tuning["id"]):
train_x, test_x = df_train_tuning.drop("TPV_t1", axis="columns").values[train_idx], df_train_tuning.drop("TPV_t1", axis="columns").values[test_idx]
train_y, test_y = df_train_tuning["TPV_t1"][train_idx], df_train_tuning["TPV_t1"][test_idx]
train_data = lgb.Dataset(train_x, label=train_y, categorical_feature=categorical_vars, feature_name=x_colnames)
model = lgb.train(params, train_data)
pred = model.predict(test_x)
mae = mean_absolute_error(pred, test_y)
mae_list.append(mae)
print("Trial done: MAE values on folds: {}".format(mae_list))
return cv_fold_func(mae_list)
# # Optuna results
study = optuna.create_study(direction="minimize")
study.optimize(objective, n_trials=80)
print("Number of finished trials:", len(study.trials))
print("Best trial:", study.best_trial.params)
best_params = study.best_trial.params
# Building parameter dictionary
fixed_params = {'objective': 'regression',
'metric': 'mae',
'boosting': 'gbdt',
'num_threads': 4}
best_params = {**best_params, **fixed_params}
# # Saving best params
filename = f"{data_dir}/data/model/parameters/best_params_t1.pickle"
outfile = open(filename, "wb")
pickle.dump(best_params, outfile)
outfile.close()
else:
# # Reading optimal parameters
infile = open(f"{data_dir}/data/model/parameters/best_params_t1.pickle", "rb")
best_params = pickle.load(infile)
Agora, com uma boa seleção de parâmetros, podemos partir para o treino do LightGBM com o objetivo de avaliar a performance nos dados que deixamos separados. Depois disso, salvamos em pickles o erro absoluto médio do modelo, a importância das variáveis, e as previsões na base de avaliação:
if TRAINING:
# # Training model
print("Splitting train and test sets...")
for train_idx, test_idx in GroupTimeSeriesSplit().split(df_train, groups=df_train["id"]):
train_x, test_x = df_train.drop("TPV_t1", axis="columns").values[train_idx], df_train.drop("TPV_t1", axis="columns").values[test_idx]
train_y, test_y = df_train["TPV_t1"][train_idx], df_train["TPV_t1"][test_idx]
print("Beginning training...")
train_data = lgb.Dataset(train_x, label=train_y, categorical_feature=categorical_vars, feature_name=x_colnames)
test_data = lgb.Dataset(test_x, label=test_y, categorical_feature=categorical_vars, feature_name=x_colnames)
lgbm_model = lgb.train(best_params, train_data, 250, valid_sets=test_data, early_stopping_rounds=50, verbose_eval=10)
# # Predicting on Test Data
# Creating validation sets for each month
df_val_t1_y = df_val["TPV_t1"]
df_val_t1_x = df_val.drop(["TPV_t1"], axis="columns")
y_pred_lgb_test = lgbm_model.predict(df_val_t1_x, num_iteration=lgbm_model.best_iteration)
# # Model performance on test data
MAE_lgb_test = mean_absolute_error(y_pred_lgb_test, df_val_t1_y)
# # Feature importances
imp = pd.DataFrame(lgbm_model.feature_importance()).T
imp.columns = x_colnames
imp = imp.T
imp.columns = ["feat_importance"]
imp = imp.sort_values("feat_importance")
# Saving Results
imp.to_csv(f"{data_dir}/data/model/performance/importances_t1.csv")
# Saving Results
filename = f"{data_dir}/data/model/performance/MAE_t1.pickle"
outfile = open(filename, "wb")
pickle.dump(MAE_lgb_test, outfile)
outfile.close()
# # Saving test prediction
# Joining test set and predictions
df_y_pred = pd.DataFrame(y_pred_lgb_test, columns=["TPV_pred_t1"])
df_val_t1_x = df_val_t1_x.reset_index()
df_val_t1_y = pd.DataFrame(df_val_t1_y)
df_val_t1_y = df_val_t1_y.reset_index().drop("mes_referencia", axis="columns")
df_pred_t1 = pd.concat([df_val_t1_x, df_val_t1_y, df_y_pred], axis="columns").set_index("mes_referencia")
df_pred_t1 = df_pred_t1[["id", "TPV_t1", "TPV_pred_t1"]]
# Saving dataframe
df_pred_t1.to_csv(f"{data_dir}/data/performance/test_pred_t1.csv")
Por último após termos decidido que estamos satisfeitos com o modelo podemos treinar o LightGBM em todos os dados e realizar a previsão finnal do desafio para o mês de Agosto:
if FINAL_PREDICTION:
print("Training on all data...")
# # Final Prediction
print("Splitting train and test sets...")
for train_idx, test_idx in GroupTimeSeriesSplit().split(df, groups=df["id"]):
train_x, test_x = df.drop("TPV_t1", axis="columns").values[train_idx], df.drop("TPV_t1", axis="columns").values[test_idx]
train_y, test_y = df["TPV_t1"][train_idx], df["TPV_t1"][test_idx]
print("Beginning training...")
train_data = lgb.Dataset(train_x, label=train_y, categorical_feature=categorical_vars, feature_name=x_colnames)
test_data = lgb.Dataset(test_x, label=test_y, categorical_feature=categorical_vars, feature_name=x_colnames)
lgbm_model = lgb.train(best_params, train_data, 250, valid_sets=test_data, early_stopping_rounds=50, verbose_eval=10)
# # Predicting on Test Data
# Creating validation sets for each month
y_pred_lgb_final = lgbm_model.predict(df_prediction_august, num_iteration=lgbm_model.best_iteration)
# # Exporting final results
y_pred_lgb_final = pd.DataFrame(y_pred_lgb_final, columns=["TPV Agosto"])
y_pred_lgb_final.to_csv(f"{data_dir}/data/predictions/prediction_august.csv")Métricas de Performance
Esta seção tem por objetivo explorar um pouco os resultados da performance do modelo na base de teste. Pretendo explorar como varia a performance do modelo de acordo com o setor ou região do cliente. Também desejo comprovar a superioridade dos resultados de meu modelo em relação à média histórica do TPV.
Vamos começar importando os pacotes, a fonte de texto que iremos utilizar para os gráficos, e os dataframes com as métricas de performance:
library(tidyverse)
library(showtext)
library(Metrics)
library(grid)
library(ggcharts)
library(mdthemes)
## Importing font
font_add(family = "Sharon", regular = "visualization/fonts/Sharon.ttf")
showtext_auto()
### Reading dataframes
df_perf = read_csv("data/model/performance/df_performance.csv") %>%
select(-1)
df_tpvs = read_csv("data/clean/stone_data_cleaned.csv")
Primeiramente, vou juntar a base de dados cadastrais e de TPVs com a de performance. Tive que realizar algumas pivotagens e manipulações para isso, como pode ser observado no código abaixo:
### Getting last five rows for each id
df_tpvs = df_tpvs %>%
arrange(id, mes_referencia) %>%
group_by(id) %>%
slice_tail(n=5) %>%
select(-mes_referencia, -TPV_mensal) %>%
unique()
### Pivoting dataframes for merge
df_tpv_real = df_perf %>%
pivot_longer(starts_with("TPV_t"), values_to="TPV_mensal", names_to="lag") %>%
select(id, TPV_mensal, lag)
df_tpv_pred = df_perf %>%
pivot_longer(starts_with("TPV_pred_t"), values_to="TPV_pred", names_to="lag") %>%
select(TPV_pred, TPV_historic_mean)
### Merging dataframes
df_comp = bind_cols(df_tpv_pred, df_tpv_real) %>%
select(id, everything()) %>%
mutate(lag=as.numeric(str_remove(lag, "TPV_t")))
rm(df_tpv_real)
rm(df_tpv_pred)
### Transforming df_comp to long format to plot graphs
df_comp_long = df_comp %>%
pivot_longer(-c(id, lag), values_to="valor", names_to="tipo")
Antes de continuar, gostaria de mostrar as primeiras linhas da base que utilizarei para os gráficos. Quero ressaltar que não há entradas para alguns clientes pois tive que selecionar apenas clientes com, no mínimo, 6 meses de histórico para treinar o modelo
head(df_comp, n=15)## # A tibble: 15 x 5
## id TPV_pred TPV_historic_mean TPV_mensal lag
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 11408. 11570. 12207. 1
## 2 1 19065. 11570. 10459. 2
## 3 1 27956. 11570. 17708. 3
## 4 1 19995. 11570. 12749. 4
## 5 1 20597. 11570. 17156. 5
## 6 3 2726. 2267. 2643. 1
## 7 3 2531. 2267. 4430. 2
## 8 3 2830. 2267. 4552. 3
## 9 3 2680. 2267. 4005. 4
## 10 3 2822. 2267. 3444. 5
## 11 4 38850. 52883. 38428. 1
## 12 4 41464. 52883. 21209. 2
## 13 4 49477. 52883. 26268. 3
## 14 4 65904. 52883. 28057. 4
## 15 4 45860. 52883. 31950. 5
Para começar, algo muito interessante que podemos visualizar, é a performance do modelo e da média histórica do TPV para cada cliente. Vou mostrar esses resultados para os primeiros clientes:
Podemos ver que, na maior parte dos casos, a reta de previsões de nosso modelo está mais próxima do valor real do que a média móvel, o que é um sinal muito promissor.
Agora, vamos juntar à nossa base os dados cadastrais dos clientes para investigarmos a performance do modelo em diferentes setores econômicos:
### Merging performance dataframe with df_comp
df_complete = inner_join(df_comp, df_tpvs, by="id") %>%
select(-lag)
### Grouping by id and getting the mean MAE for each client
df_complete = df_complete %>%
group_by(id) %>%
summarise(id, MAE_model = mae(TPV_mensal, TPV_pred), MAE_mean=mae(TPV_historic_mean, TPV_mensal),
MacroClassificacao, Estado, segmento) %>% unique()
Agora que temos todos os dados juntos, vamos investigar a performance de nosso modelo de acordo com a Macro Classificação:
Podemos observar que nosso modelo possui um erro absoluto médio bastante inferior ao da média histórica na base de teste, o que é um ótimo sinal. Também podemos ver que o setor de postos de gasolina é o que o modelo tem mais dificuldades de prever. Já o setor de serviços foi onde o modelo teve sua melhor performance.
Vamos detalhar mais a performance do modelo e investigar seu comportamento de acordo com a variável “segmento”:
Nessa visão menos agrupada, podemos ver que a performance do modelo é bem elevada no setor de Petshops e Veterinárias, Serviços de Beleza e Estética, e Autopeças e Serviços Autimotivos. Mais uma vez, o setor de Postos de Gasolina demonstra ser difícil de prever para o modelo, sendo superado pelo de Logística e Mobilidade. Muito provavelmente, a baixa performance em relação a esses setores tem relação com o preço dos combustíveis, uma variável que pretendo incluir no modelo caso eu prossiga no desafio.
Agora que já vimos a performance do modelo de acordo com o segmento econômico, vamos observar a performance de acordo com o estado da federação para ver se o modelo tem dificuldades particulares em algum estado:
É possível observar que nosso modelo apresenta uma variabilidade muito menor em relação ao erro absoluto médio quando comparado com a performance da média histórica, um sinal muito promissor. Mesmo assim, vemos que o modelo possui uma considerável diferença de performance quando comparamos estados como o Espírito Santo e Piauí. Uma variável que poderíamos adicionar ao modelo para tentar captar essas diferenças é o número de casos de covid no último mês, porque sabemos que a pandemia não ocorre de forma homogênea no território nacional e, certamente, tem um impacto nas atividades econômicas.
Para finalizar, vamos comparar os gráficos de dispersão entre as previsões de nosso modelo com os valores reais do TPV e entre a média histórica do TPV e os valores reais: 
Agora vamos juntar as duas linhas de regressão em um só gráfico para ver qual possui a maior inclinação:
Podemos ver que a inclinação da reta do modelo (verde) é, claramente superior. Demonstrando que, pelo menos, para a base de teste nosso modelo claramente supera a performance da média histórica do TPV.
Considerações
Gostaria de dizer que foi um grande prazer participar deste desafio. Independentemente de qualquer resultado ou métrica de performance, aprendi bastante sobre problemas de séries temporais e alguns métodos específicos a eles. Espero que o conteúdo tenha sido interessante e prazeroso para vocês, assim como foi para mim.
