Grupo Bimbo Inventory Demand¶
06 de março, 2020
1. Descrição geral do problema¶
O Grupo Bimbo, se esforça para atender a demanda diária dos consumidores por produtos frescos de panificação nas prateleiras de mais de 1 milhão de lojas ao longo das suas 45.000 lojas em todo o México.
Atualmente, os cálculos diários de estoque são realizados por funcionários de vendas de entregas diretas, que devem, sozinhos, prever a necessidade de estoque dos produtos e demanda com base em suas experiências pessoais em cada loja. Como alguns pães têm uma vida útil de uma semana, a margem aceitável para o erro é pequena.
Objetivo: neste projeto de aprendizado de máquina, vamos desenvolver um modelo para prever com precisão a demanda de estoque com base nos dados históricos de vendas. Isso fará com que os consumidores dos mais de 100 produtos de panificação não fiquem olhando para as prateleiras vazias, além de reduzir o valor gasto com reembolsos para os proprietários de lojas com produtos excedentes impróprios para venda.
Para a construção desse projeto, utilizaremos a linguagem R e os datasets disponíveis no Kaggle em:
2. Carregando dados¶
2.1 Importando bibliotecas necessárias¶
Vamos começar nosso projeto importanto todas as bilbiotecas necessárias para a realização das fases iniciais de exploração e transformação dos dados (Data Munging).
# Caso não possua uma das bibliotecas importadas abaixo, a instale com um dos comandos a seguir:
install.packages(c(
'data.table',
'bigreadr',
'dplyr',
'ggplot2',
'fasttime',
'lubridate',
'corrplot',
'anomalize',
'stringr'
))
# Definindo a oculatação de warnings.
options(warn = -1)
# Importando bibliotecas.
library(data.table)
library(bigreadr)
library(dplyr)
library(ggplot2)
library(fasttime)
library(lubridate)
library(corrplot)
library(anomalize)
library(stringr)
2.2 Carregando dados do dataset cliente_tabla¶
# Importando dataset.
client <- fread('/content/datasets/cliente_tabla.csv')
# Verificando as primeiras linhas do dataset.
head(client)
2.3 Carregando dados do dataset producto_tabla¶
# Importando dataset.
product <- fread('/content/datasets/producto_tabla.csv')
# Verificando as primeiras linhas do dataset.
head(product)
2.4 Carregando dados do dataset town_state¶
# Importando dataset.
town <- fread('/content/datasets/town_state.csv')
# Verificando as primeiras linhas do dataset.
head(town)
2.5 Carregando dados de treino¶
# Importando dataset.
train <- fread('/content/datasets/train.csv')
# Verificando as primeiras linhas do dataset.
head(train)
2.6 Carregando dados de teste¶
# Importando dataset.
test <- fread('/content/datasets/test.csv')
# Verificando as primeiras linhas do dataset.
head(test)
3. Data Munging - Eliminando inconsistências nos datasets¶
A documentação nos alerta para a existência de alguns problemas que devem ser tratados dentro do conjunto de dados, como por exemplo registros duplicados. Por isso, iremos fazer uma breve exploração dos datasets para eliminar todas as inconsistências que possuam.
3.1 Dataset client¶
# Visualizando as primeiras 10 linhas do dataset.
head(client, 10)
Podemos observar que existem IDs que se repetem e nomes de clientes desconhecidos (denominados como "SIN NOMBRE") no conjunto de dados que precisarão ser tratados.
Vamos contabilizar o número de registros duplicados no dataset.
# Verificando linhas duplicadas no dataset.
table(duplicated(client))
Não encontramos nenhum registro duplicado o que nos leva a crer que a variável NombreCliente apresenta strings com diferentes tamanhos para cada ID duplicado. Vamos confirmar esta teoria.
Ao listar as primeira linhas do conjunto de dados, vimos que o ID 4 se repete duas vezes. Com base nisso, vamos capturar e analisar o valor da variável NombreCliente associado a este ID em cada observação.
# Definindo o número do ID que deve ser capturado.
id <- 4
# Capturando linhas que contenham o ID especificado.
client[client$Cliente_ID == id,]
# Capturando nomes associados a cada um dos registros duplicados.
fName <- client[client$Cliente_ID == id,][1, 'NombreCliente']
sName <- client[client$Cliente_ID == id,][2, 'NombreCliente']
# Definindo o número de caracteres de cada um dos nomes dos registros duplicados.
nchar(fName)
nchar(sName)
A partir deste resultado podemos confirmar que há uma diferença entre os valores das variáveis NombreCliente de cada registro duplicado. Isto provavelmente ocorre devido a diferença do número de espaços existentes em cada nome.
Iremos contabilizar o número de registros duplicados a partir da variável Cliente_ID.
# Verificando número de IDs duplicados no dataset.
table(duplicated(client$Cliente_ID))
# Removendo registros com número de ID duplicado.
client <- client[!duplicated(client$Cliente_ID),]
Definiremos o número de registros sem o nome do cliente.
# Verificando número de registros sem o nome do cliente.
nrow(client[client$NombreCliente == 'SIN NOMBRE', ])
Há 356 observações sem o nome do cliente.
# Verificando se existem valores nulos no dataset.
anyNA(client)
Não há valores nulos no dataset.
# Verificando o tipo de dado das variáveis do dataset.
glimpse(client)
O conjunto de dados contém registros de 930.500 clientes.
3.2 Dataset product¶
# Visualizando as primeiras 10 linhas do dataset.
head(product, 10)
Vamos contabilizar o número de registros duplicados no dataset.
# Verificando linhas duplicadas no dataset.
table(duplicated(product))
Nenhuma observação duplicada foi encontrada.
# Verificando número de IDs duplicados no dataset.
table(duplicated(product$Producto_ID))
Nenhum número de ID duplicado foi encontrado.
Porém, o produto com o ID 0 não possui o nome do produto (possui o valor "NO IDENTIFICADO"). Iremos verificar se isto só ocorre neste registro.
# Capturando a substring especificada.
pattern <- "NO IDENTIFICADO"
# Definindo linhas que não identifiquem o nome do produto.
rows <- grep(pattern, product[, 'NombreProducto'], value = F)
# Visualizando linhas que não contenham o nome do produto.
product[rows, ]
Concluímos que apenas o ID 0 não apresenta um nome de produto. Talvez vez isto ocorra porque um único ID pode estar sendo utilizado para identificar produtos que ainda não foram devidamente registrados no conjunto de dados.
# Verificando se existem valores nulos no dataset.
anyNA(product)
Não há valores nulos no dataset.
# Verificando o tipo de dado das variáveis do dataset.
glimpse(product)
O conjunto de dados contém registros de 2.592 produtos.
Repare que a variável NombreProducto contém outras informações além do nome do produto. Parece que a string segue o seguinte padrão:
| NombreProducto | Nome do produto | Número de peças | Peso | Sigla do fabricante | ID do produto |
Veja que este padrão não está presente em todos os valores da variável, mas predomina em grande parte dos dados. Bom, não precisamos de todas estas informações para a análise que iremos fazer e por isso iremos extair apenas o nome, o peso e a sigla do fabricante de cada produto.
## Extraindo a unidade de medida de massa do produto.
# Extraindo a substring com as informações brutas para uma variável temporária.
tmp <- str_extract(product$NombreProducto, "([0-9 ] |[0-9])+(G|g|Kg|kg|ml)")
# Criando uma varíavel para armazenar o número associado ao peso do produto.
product$Npeso <- as.integer(str_extract(tmp, "[0-9]+"))
# Criando uma varíavel para armazenar a unidade de medida do peso do produto.
product$UniPeso <- tolower(str_extract(tmp, "[A-z]+"))
# Criando uma variável para armazenar a sigla referente ao fabricante.
product$Productor <- toupper(str_extract(
str_extract(
product$NombreProducto,
"( [A-Z]+[a-z ]+[A-Z]+ [A-Z ]+ [0-9]+$| [A-Z ]+[A-Z ]+ [0-9]+$)"),
"( [A-Z]+[a-z ]+[A-Z]+ [A-Z ]+ | [A-Z ]+[A-Z ]+ )"
))
# Extraindo o nome do produto.
product$NombreProducto <- str_extract(product$NombreProducto, "[A-z ]+")
# Visualizando dataset após a extração das informações desejadas.
head(product)
# Verificando se existem valores nulos em cada variável do dataset.
sapply(product, function(v){
table(is.na(v))
})
Como resultado final, verificamos que não foi possível determinar o peso de 51 produtos e nem o sigla de 1 fabricante.
3.3 Dataset town¶
# Visualizando as primeiras 10 linhas do dataset.
head(town, 10)
Vamos verificar se há registros ou IDs de agência duplicados no dataset.
# Verificando linhas duplicadas no dataset.
table(duplicated(town))
Nenhum registro duplicado foi encontrado.
# Verificando número de IDs duplicados no dataset.
table(duplicated(town$Agencia_ID))
Nenhum registro contém um número de ID duplicado.
# Verificando se existem valores nulos no dataset.
anyNA(town)
Não há valores nulos no dataset.
# Verificando o tipo de dado das variáveis do dataset.
glimpse(town)
O conjunto de dados contém registros de 790 cidades e seus respectivos estados.
3.4 Dataset train¶
# Visualizando as primeiras 10 linhas do dataset.
head(train, 10)
# Verificando linhas duplicadas no dataset.
table(duplicated(train))
Não há registros duplicados no dataset.
# Verificando se existem valores nulos no dataset.
anyNA(train)
Não há valores nulos no dataset.
# Verificando o tipo de dado das variáveis do dataset.
glimpse(train)
O conjunto de dados de treino contém 74.180.464 registros e 11 colunas.
3.5 Dataset test¶
# Visualizando as primeiras 10 linhas do dataset.
head(test, 10)
# Verificando linhas duplicadas no dataset.
table(duplicated(test))
Não há registros duplicados no dataset.
# Verificando se existem valores nulos no dataset.
anyNA(test)
Não há valores nulos no dataset.
# Verificando o tipo de dado das variáveis do dataset.
glimpse(test)
O conjunto de dados de treino contém 6.999.251 registros e 7 colunas.
4. Análise exploratória dos dados¶
4.1 Visão geral¶
Segundo a documentação referente ao projeto, cada linha dos dados de treinamento contém um registro de vendas com as seguintes variáveis:
| Variável | Descrição |
|---|---|
| Semana | é o número da semana (de quinta a quarta-feira); |
| Agencia_ID | é o ID do depósito de vendas; |
| Canal_ID | é o ID do canal de vendas; |
| Ruta_SAK | é o ID da rota (várias rotas = depósito de vendas); |
| Cliente_ID | é o ID do cliente; |
| NombreCliente | é o nome do cliente; |
| Producto_ID | é o ID do produto; |
| NombreProducto | é o nome do produto; |
| Venta_uni_hoy | é o número de unidades vendidas na semana; |
| Venta_hoy | é o valor de venda na semana (unidade monetária: pesos); |
| Dev_uni_proxima | é o número de unidades retornadas na próxima semana; |
| Dev_proxima | é o valor retornado na próxima semana e (unidade monetária: pesos) e; |
| Demanda_uni_equil (Target) | é a variável a ser prevista, define a demanda ajustada. |
Nesta etapa vamos buscar entender a disposição e as características dos dados dentro do dataset de treino além de extrair insigths que possam auxiliar no processo de criação do modelo preditivo.
# Verificando o tipo de dado das variáveis do dataset.
glimpse(train)
O conjunto de dados de treino contém 74.180.464 registros e 11 colunas. Todas as variáveis apresentam o tipo de dado numérico.
# Verificando o número de valores únicos para cada variável do dataset.
t <- sapply(train, function(c) {
length(unique(c))
})
# Exibindo os resultados.
print(t)
Estes resultados são interessantes. Destacamos que podemos observar a existência de registros de 7 semanas; 880604 IDs de clientes e; 1799 IDs de Produtos diferentes no dataset.
4.2 Analisando cada variável separadamente¶
4.2.1 Criando funções auxiliares¶
Criaremos algumas funções para padronizar as plotagens de gráficos que efetuaremos.
# Definindo uma função para criar gráficos de barra.
barPlot <- function(col, data) {
data %>%
mutate_at(c(var = col), as.factor) %>%
group_by(var) %>%
summarise(absFreq = n()) %>%
ggplot(aes(x = var, y = absFreq)) +
geom_bar(stat = 'identity', alpha = 0.75, fill = '#086788') +
ylab('Frequency') +
xlab(col) +
labs(title = paste('Bar plot for variable:', col)) +
theme_bw()
}
# Definindo uma função para criar gráficos de boxplot.
boxPlot <- function(col, data) {
data %>%
group_by_at(col) %>%
summarise(absFreq = n()) %>%
ggplot(aes(x = absFreq)) +
geom_boxplot(fill = '#566E3D', color = '#373D20', alpha = 0.8) +
theme_bw() +
theme(axis.text.y = element_blank()) +
xlab(paste(col, 'frequency')) +
labs(title = paste('Boxplot for variable:', col))
}
4.2.2 Variável Semana¶
# Definindo o nome da variável a ser analisada.
col <- 'Semana'
# Criando um gráfico de barras para a variável especificada.
barPlot(col, train)
O gráfico nos permite observar há existênica de uma distribuição aproximadamente uniforme dos registros entre as semanas.
# Contabilizando o número de registros por semana.
table(Semana = train$Semana)
4.2.3 Variável Agencia_ID¶
# Criando um boxplot com as frequências com que os ID das agências aparecem no conjunto de dados.
# Definindo o nome da variável a ser analisada.
col <- 'Agencia_ID'
# Criando um gráfico boxplot para a variável especificada.
boxPlot(col, train)
Vemos que há outliers nas frequências com que os IDs das agências aparecem no conjunto de dados. Vamos determinar a cidade e o estado destas agências com que a frequência do ID foi discrepante.
# Determinando o endereço, o estado e ordenando de forma decrescente a frequência dos IDs das agências no dataset.
AgencyIdFreq <- train %>%
select(Agencia_ID) %>%
group_by(Agencia_ID) %>%
summarise(absFreq = n()) %>%
arrange(desc(absFreq)) %>%
inner_join(town, by = 'Agencia_ID')
# Visualizando os 5 IDs de agências mais frequentes.
head(AgencyIdFreq)
# Extraindo outliers do dataset.
AgencyIdStates <- AgencyIdFreq %>%
anomalize(absFreq, method = 'iqr', alpha = 0.10) %>%
filter(anomaly == 'Yes') %>%
select(-c(absFreq_l1, absFreq_l2, anomaly))
# Visualizando as informações dos outliers.
AgencyIdStates
Parece haver um grupo de estados mais frequente dentro das informações destas agências discrepantes. Para facilitar esta análise, iremos criar um gráfico de barras.
# Criando um gráfico de barras para os estados das agências que apresentaram uma frequência discrepante.
# Definindo o nome da variável a ser analisada.
col <- 'State'
# Criando um gráfico de barras para a variável especificada.
barPlot(col, AgencyIdStates)
# Determinando a proporção dos estados identificados entre os IDs das agências que apresentaram uma frequência discrepante.
prop.table(table(AgencyIdStates$State))
Concluímos que o Estado do México apresenta 6 das agências mais recorrentes no conjunto de dados e que o estado de Jalisco é o que possui a agência mais frequente.
4.2.4 Variável Canal_ID¶
# Plotando um gráfico de barras para o conjunto de dados da variável.
# Definindo o nome da variável a ser analisada.
col <- 'Canal_ID'
# Criando um gráfico de barras para a variável especificada.
barPlot(col, train)
Observamos que o canal de vendas com ID 1 é o mais frequente. Determinaremos a proporção da frequência de cada um destes canais dentro do conjunto de dados.
# Determinando a propoção de cada um dos canais.
train %>%
mutate(Canal_ID = as.factor(Canal_ID)) %>%
group_by(Canal_ID) %>%
summarise(Prop = round(n() / length(train$Canal_ID) * 100, digits = 3))
Concluímos que aproximadamente 91% dos registros do conjunto de dados está associado ao Canal com ID = 1.
4.2.5 Variável Ruta_SAK¶
# Criando um boxplot com as frequências com que as rotas usadas aparecem no conjunto de dados.
# Definindo o nome da variável a ser analisada.
col <- 'Ruta_SAK'
# Criando um gráfico de boxplot para a variável especificada.
boxPlot(col, train)
Há um grande número de outliers dentro desta variável. Pode ser interessante extraí-los e análisá-los separadamente.
# Determinando em ordem decrescente as rotas mais frequentes dentro do conjunto de dados.
routeFreq <- train %>%
group_by(Ruta_SAK) %>%
summarise(absFreq = n()) %>%
arrange(desc(absFreq))
# Visualizando as primeiras linhas do dataset.
head(routeFreq)
# Extraindo outliers do dataset.
routeOutFreq <- routeFreq %>%
anomalize(absFreq, method = 'iqr', alpha = 0.10) %>%
filter(anomaly == 'Yes') %>%
select(-c(absFreq_l1, absFreq_l2, anomaly))
# Determinando número de outliers.
nrow(routeOutFreq)
Bom, detectamos a existência de 605 rotas com frequências discrepantes dentro do conjunto de dados.
# Determinando a rota mais frequente dentro do conjunto de dados.
routeOutFreq[routeOutFreq$absFreq == max(routeOutFreq$absFreq), ]
Constatamos que a rota 1201 é a mais recorrente no dataset.
# Determinando a proporção de rotas com frequências discrepantes.
length(routeOutFreq$Ruta_SAK) / length(unique(train$Ruta_SAK)) * 100
# Determinando a proporção de registros associados as rotas discrepantes.
sum(routeOutFreq$absFreq) / length(train$Ruta_SAK) * 100
Por fim, concluímos que aproximadamente 16.8% das rotas apresentam frequências discrepantes e são responsáveis por 86.38% das entregas.
4.2.6 Variável Cliente_ID¶
# Criando um boxplot com as frequências com que os IDs dos clientes que estão presentes no dataset aparecem.
# Definindo o nome da variável a ser analisada.
col <- 'Cliente_ID'
# Criando um gráfico de boxplot para a variável especificada.
boxPlot(col, train)
Veja que interessante, há um cliente que apresenta uma frequência extremamente alta em relação aos demais e acaba distorcendo o gráfico. Vamos identificar o nome deste cliente junto com os dos demais outliers.
# Determinando em ordem decrescente os clientes mais frequentes dentro do conjunto de dados.
clientFreq <- train %>%
select(Cliente_ID) %>%
group_by(Cliente_ID) %>%
summarise(absFreq = n()) %>%
arrange(desc(absFreq)) %>%
inner_join(client, by = 'Cliente_ID')
# Visualizando as primeiras linhas do dataset.
head(clientFreq)
O cliente mais discrepante possui uma frequência que é aproximadamente 23.4 vezes maior do que a do cliente que ocupa a segunda posição nos dando uma breve noção do quão distante este primeiro colocado está dos demais.
# Extraindo outliers do dataset.
clientOutFreq <- clientFreq %>%
anomalize(absFreq, method = 'iqr', alpha = 0.0415) %>%
filter(anomaly == 'Yes') %>%
select(-c(absFreq_l1, absFreq_l2, anomaly))
# Determinando número de outliers.
nrow(clientOutFreq)
Verificamos a existência de 1622 IDs de clientes com frequências discrepantes dentro do conjunto de dados.
# Determinando o cliente mais frequente dentro do conjunto de dados.
mostFrequentClient <- clientOutFreq[clientOutFreq$absFreq == max(clientOutFreq$absFreq), ]
# Visualizando o cliente mais frequente dentro do conjunto de dados.
mostFrequentClient
Identificamos que o nome do cliente que possuí o ID mais recorrente dentro do dataset é "Puebla Remision".
# Determinando a proporção de registros que contém o ID do Cliente com a frequência mais discrepante.
mostFrequentClient$absFreq / length(train$Cliente_ID) * 100
O cliente "Puebla Remision" está associado a aproximadamente 0.167% dos registros do conjunto de dados.
# Determinando a proporção de registros que contém os IDs dos Clientes com as frequências mais discrepantes.
sum(clientOutFreq$absFreq) / length(train$Cliente_ID) * 100
Todos os registros associados a IDs de clientes que possuem uma frequência discrepante correspondem a aproximadamente 1.37% do total. Isto nos indica que a maior parte dos dados que estamos manipulando estão relacionados a muitos clientes que efetuam compras com uma frequência que não foge do padrão.
4.2.7 Variável Producto_ID¶
# Criando um boxplot com as frequências com que os IDs dos produtos que estão presentes no dataset aparecem.
# Definindo o nome da variável a ser analisada.
col <- 'Producto_ID'
# Criando um gráfico de boxplot para a variável especificada.
boxPlot(col, train)
Existem muitos outliers entre as frequências dos produtos. Isto indica que há um subconjunto de itens que fogem do padrão de venda do resto dos demais produtos.
# Determinando em ordem decrescente os produtos mais frequentes dentro do conjunto de dados.
productFreq <- train %>%
select(Producto_ID) %>%
group_by(Producto_ID) %>%
summarise(absFreq = n()) %>%
arrange(desc(absFreq)) %>%
inner_join(product, by = 'Producto_ID')
# Visualizando as primeiras linhas do dataset.
head(productFreq)
# Extraindo outliers do dataset.
productOutFreq <- productFreq %>%
anomalize(absFreq, method = 'iqr', alpha = 0.1) %>%
filter(anomaly == 'Yes') %>%
select(-c(absFreq_l1, absFreq_l2, anomaly))
# Determinando número de outliers.
nrow(productOutFreq)
Observamos que existem 333 frequências de produtos discreprantes.
# Determinando o produto mais frequente dentro do conjunto de dados.
mostFrequentProduct <- productOutFreq[productOutFreq$absFreq == max(productOutFreq$absFreq), ]
# Visualizando o produto mais frequente dentro do conjunto de dados.
mostFrequentProduct
Detectamos que o produto que mais apresenta registros de vendas dentro do nosso conjunto de dados é denominado "Mantecadas Vainilla".
# Determinando a proporção de registros que contém os IDs dos produtos com as frequências discrepantes.
sum(productOutFreq$absFreq) / length(train$Producto_ID) * 100
Os produtos que apresentam uma frequência discrepante dentro do conjunto de dados são responsáveis por aproximadamente 96.76% dos registros.
# Determinando a sigla do fabricante mais recorrente dentro do conjunto de dados que contém os IDs dos produtos
# com as frequências discrepantes.
manufacturersOut <- productOutFreq %>%
group_by(Productor) %>%
summarise(absFreq = n()) %>%
arrange(desc(absFreq))
# Visualizando as primeiras linhas do dataset.
head(manufacturersOut)
Concluímos que o fabricante identificado pela sigla BIM é o mais recorrente dentro dos produtos que apresentam uma frequência discrepante.
4.2.8 Variável Venta_uni_hoy¶
# Verificando a distribuição dos dados.
summary(train$Venta_uni_hoy)
# Verificando a frequência com que cada número de unidades aparece no dataset.
t <- train %>%
group_by(Venta_uni_hoy) %>%
summarise(absFreq = n())
# Ordenando dados das frequências em ordem decrescente.
t <- t[order(t$absFreq, decreasing = T), ]
# Visualizando o número de unidades mais frequentes no dataset.
head(t, 10)
Concluímos que as vendas com 2 unidades são as mais frequentes dentro do conjunto de dados e que o número das 10 unidades mais frequentes varia entre 1 e 10.
4.2.9 Variável Venta_hoy¶
# Verificando a distribuição dos dados.
summary(train$Venta_hoy)
# Definindo o valor total de vendas por semana.
train %>%
group_by(Semana) %>%
summarise(total_Venta_hoy = sum(Venta_hoy))
# Definindo o valor total de vendas por cliente.
train %>%
group_by(Cliente_ID ) %>%
summarise(total_Venta_hoy = sum(Venta_hoy)) %>%
arrange(desc(total_Venta_hoy)) %>%
head(10)
4.2.10 Variável Dev_uni_proxima¶
# Verificando a distribuição dos dados.
summary(train$Dev_uni_proxima)
# Definindo o número total de unidades retornadas na próxima semana.
train %>%
group_by(Semana) %>%
summarise(total_Dev_uni_proxima = sum(Dev_uni_proxima))
# Definindo o número total de unidades retornadas na próxima semana por cliente.
train %>%
group_by(Cliente_ID) %>%
summarise(total_Dev_uni_proxima = sum(Dev_uni_proxima)) %>%
arrange(desc(total_Dev_uni_proxima)) %>%
head()
4.2.11 Variável Dev_proxima¶
# Verificando a distribuição dos dados.
summary(train$Dev_proxima)
# Definindo o valor total retornado na próxima semana.
train %>%
group_by(Semana) %>%
summarise(total_Dev_proxima = sum(Dev_proxima))
# Definindo o valor total retornado na próxima semana por cliente.
train %>%
group_by(Cliente_ID) %>%
summarise(total_Dev_proxima = sum(Dev_proxima)) %>%
arrange(desc(total_Dev_proxima)) %>%
head()
4.2.12 Variável Demanda_uni_equil¶
# Criando um boxplot para visualizar a distribuição dos dados da variável Demanda_uni_equil.
# Definindo o nome da variável a ser analisada.
col <- 'Demanda_uni_equil'
# Criando um gráfico de boxplot para a variável especificada.
train %>%
ggplot(aes(x = Demanda_uni_equil)) +
geom_boxplot(fill = '#566E3D', color = '#373D20', alpha = 0.8) +
theme_bw() +
theme(axis.text.y = element_blank()) +
xlab('Values') +
labs(title = 'Boxplot for variable: Demanda uni equil')
Vemos que há uma distorção muito grande dentro dos valores da variável a ser prevista devido a presença de outliers muito discrepantes. Deveremos tratar este problema para conseguirmos aumentar a performance dos modelos preditivos que criarmos.
4.2.13 Variável Town¶
Agora, iremos verificar quantos e quais são os estados e as agências por cidades presentes no conjunto de dados.
# Contabilizando a frequência de agências por cidade dentro do dataset.
t <- town %>%
group_by(Town) %>%
summarise(absFreq = n())
# Ordenando resultados.
t <- t[order(t$absFreq, decreasing = T),]
# Visualizando as primeiras 10 linhas da tabela.
head(t, 10)
# Determinando o número de agências por cidade atendidas.
nrow(t)
Detectamos 260 agências por cidades atendidas.
Criaremos um boxplot para verificar se existem outliers dentro das frequências de agências por cidade.
# Definindo o nome da variável a ser analisada.
col <- 'Town'
# Criando um gráfico de boxplot para a variável especificada.
boxPlot(col, town)
Concluímos que a agência 2013 AG. MEGA NAUCALPAN apresenta um frequência absoluta que saí do pradrão das demais presentes no conjunto de dados.
4.2.14 Variável State¶
# Contabilizando a frequência de cada estado dentro do dataset.
t <- town %>%
group_by(State) %>%
summarise(absFreq = n())
# Ordenando resultados.
t <- t[order(t$absFreq, decreasing = T),]
# Visualizando as primeiras 10 linhas da tabela.
head(t, 10)
# Determinando o número de estados atendidos.
nrow(t)
Detectamos que 33 estados são atendidos.
# Definindo o nome da variável a ser analisada.
col <- 'State'
# Criando um gráfico de boxplot para a variável especificada.
boxPlot(col, town)
Concluímos que o ESTADO DE MÉXICO apresenta um frequência absoluta que saí do pradrão dos demais presentes no conjunto de dados.
5. Análise Preditiva¶
5.1 Importando bibliotecas necessárias¶
Importaremos todas as bilbiotecas necessárias para a realização dos processos de modelagem preditiva.
# Caso não possua uma das bibliotecas importadas abaixo, a instale com um dos comandos a seguir:
install.packages(c(
'Metrics',
'xgboost',
'randomForest',
'caret'
))
# Importando bibliotecas.
library(Metrics)
library(xgboost)
library(randomForest)
library(caret)
5.2 Feature Selection¶
Observe que as variáveis Venta_uni_hoy, Venta_hoy, Dev_uni_proxima e Dev_proxima não estão presentes no conjunto de dados de teste e por isso iremos excluí-las do conjunto de dados de treino.
# Selecionando as variáveis que serão utilizadas na fase de modelagem preditiva dentro do dataset de treino.
train <- train %>% select(Semana, Agencia_ID, Canal_ID, Ruta_SAK, Cliente_ID, Producto_ID, Demanda_uni_equil)
5.3 Feature Engineering I - Transformando variável Target¶
Agora, iremos retornar ao problema da distorção dos valores do conjunto de dados da variável alvo. Criaremos mais uma vez um boxplot para visualizar a distribuição dos dados bem como um gráfico de densidade.
# Criando um boxplot para visualizar a distribuição dos dados da variável Demanda_uni_equil.
train %>%
ggplot(aes(x = Demanda_uni_equil)) +
geom_boxplot(fill = '#566E3D', color = '#373D20', alpha = 0.8) +
theme_bw() +
theme(axis.text.y = element_blank()) +
xlab('Values') +
labs(title = 'Boxplot for variable: Demanda uni equil')
# Criando um gráfico de densidade para visualizar a distribuição dos dados da variável Demanda_uni_equil.
train %>%
ggplot(aes(x = Demanda_uni_equil)) +
geom_density(fill = '#A6EBC9') +
theme_bw() +
labs(title = 'Density graph for variable: Demanda uni equil') +
xlab('Demanda uni equil')
Muito bem, para contornar este problema da distorção nos dados usaremos a função de transformação log1p (ou log(x + 1)) para diminuir a irregularidade nos dados e tornar os padrões que esta variável possuí mais visíveis. Também utilizaremos a função expm1 (ou exp(x) - 1) para realizar o processo inverso de transformação dos resultados obtidos ao se aplicar a função log1p.
Ou seja, transformaremos a variável a ser prevista para realizar a execução da análise preditiva e no final converteremos os resultados gerados para a sua escala original.
Para mais informações sobre como a função log atua sobre dados distorcidos consulte este link. Para entender melhor as funções log1p e expm1 consulte este link.
Destacamos que o principal motivo de utilizarmos a função log1p é há existência de valores nulos dentro dos dados da variável o que inviabiliza o uso da função log pois o log(0) é um valor indefinido.
# Verificando a existência de valores nulos dentro do conjunto de dados da variável Demanda_uni_equil.
prop.table(table(train$Demanda_uni_equil == 0))
Detectamos que aproximadamente 1.8% dos dados da variável a ser prevista são iguais a 0.
Vamos aplicar a função log1p ao conjunto de dados.
# Calcula o logaritmo natural de cada valor da variável Demanda_uni_equil acrescido de 1 unidade, ou seja,
# log(Demanda_uni_equil + 1).
train$Demanda_uni_equil <- log1p(train$Demanda_uni_equil)
Criaremos mais uma vez um gráfico de boxplot e um gráfico de densidade para visualizar o efeito da transformação aplicada sobre os dados.
# Criando um boxplot para visualizar a distribuição dos dados da variável Demanda_uni_equil transformada.
train %>%
ggplot(aes(x = Demanda_uni_equil)) +
geom_boxplot(fill = '#566E3D', color = '#373D20', alpha = 0.8) +
theme_bw() +
theme(axis.text.y = element_blank()) +
xlab('log(Demanda uni equil + 1)') +
labs(title = 'Boxplot for variable: Demanda uni equil')
# Criando um gráfico de densidade para visualizar a distribuição dos dados da variável Demanda_uni_equil transformada.
train %>%
ggplot(aes(x = Demanda_uni_equil)) +
geom_density(fill = '#A6EBC9') +
theme_bw() +
labs(title = 'Density graph for variable: Demanda uni equil') +
xlab('log(Demanda uni equil + 1)')
Concluímos que a aplicação da função log1p diminuiu a distorção causada pelos valores discrepantes dentro do conjunto de dados da variável a ser prevista e isso nos ajudará a alçancar valores de acurácia melhores nos modelos que criarmos.
5.4 Unindo dados de treino e de teste em um mesmo dataset¶
Nosso objetivo nesta etapa é criar um único dataset contendo tanto os dados de treino quanto os dados de teste. Mas, antes de executarmos esta ação devemos observar que há variáveis excluivas para cada um destes conjuntos de dados.
O dataset de teste possui a variável id que não está contida nos dados de treino e por isso a criaremos para este conjunto de dados. O mesmo processo será efetuado para a variável Demanda_uni_equil no dataset de teste.
Para distinguir os registros que pertencem ao conjunto de dados de treino dos que pertencem aos dados de teste criaremos uma variável binária denominada toTest (0: são dados de treino; 1: são dados de teste).
# Criando variável ID para o conjunto de dados de treino com um valor auxiliar.
train$id <- 0
# Criando a variável para indicar se o registro pertence ou não ao dados de treino ou aos dados de teste.
train$toTest <- 0
# Criando variável Demanda_uni_equil para o conjunto de dados de teste com um valor auxiliar.
test$Demanda_uni_equil <- 0
# Criando a variável para indicar se o registro pertence ou não ao dados de treino ou aos dados de teste.
test$toTest <- 1
Para começar esta junção entre os datasets, iremos capturar os registros de apenas uma das semanas presentes no dataset de treino e unir com os dados de teste.
Os registros dos dados de treino que não irão ser manipulados nesta etapa serão utilizados na etapa a seguir de Feature Engennier.
# Unindo os registros do dataset de treino em que a variável Semana é igual a 9 com todos os registros de teste.
data <- rbind(train[Semana == 9], test)
Como já mesclamos todos os registros do dataset de teste, podemos liberar memória excluindo a variável test.
# Removendo o dataset test.
rm(test)
5.4.1 Feature Engennier II - Criando novas variáveis preditoras¶
Os registros de treino das semanas que não foram usados na etapa anterior serão agrupados ao novo dataset que estamos gerando para criar novas variáveis.
# Determinando a média da demanda ajustada de clientes por produto e a quantidade de registros de clientes por produto.
train[Semana <= 8][ , .(meanClientProd = mean(Demanda_uni_equil), countClientProd = .N),
by = .(Producto_ID, Cliente_ID)] %>%
merge(data, all.y = TRUE, by = c("Producto_ID", "Cliente_ID")) -> data
# Determinando a média da demanda ajustada por produto e a quantidade de registros por produto.
train[Semana <= 8][ , .(meanProd = mean(Demanda_uni_equil), countProd = .N),
by = .(Producto_ID)] %>%
merge(data, all.y = TRUE, by = c("Producto_ID")) -> data
# Determinando a média da demanda ajustada por cliente e a quantidade de registros por cliente.
train[Semana <= 8][ , .(meanClient = mean(Demanda_uni_equil), countCliente = .N),
by = .(Cliente_ID)] %>%
merge(data, all.y = TRUE, by = c("Cliente_ID")) -> data
# Visualizando as primeiras linhas do dataset.
head(data)
Observe que a partir da execução desta fase conseguimos eliminar o problema da existência de produtos dentro dos dados de teste que não estão presentes dentro dos dados de treino, pois passamos a avaliar os valores médios e quantidades de cada variável por agrupamentos.
Agora também podemos eliminar a variável train.
# Removendo o dataset train.
rm(train)
5.5 Feature Engineering III - Transformando variáveis preditoras¶
Nesta etapa vamos escalar os valores das variáveis preditoras entre 0 e 1.
# Definindo método de pré-processamento.
params <- preProcess(data[, !c('id', 'Demanda_uni_equil', 'toTest')], method = 'range')
# Transformando os dados.
data <- predict(params, data)
# Visualizando as primeiras linhas do dataset.
head(data)
5.6 Segmentando dados de treino e de teste¶
Iremos separar os dados de treino e de teste do conjunto de dados que criamos nas etapas anteriores.
# Extraindo registros de treino.
train <- data %>%
filter(toTest == 0) %>%
select(-c(id, toTest))
# Visualizando as primeiras linhas do dataset.
head(train)
# Extraindo registros de teste.
test <- data %>%
filter(toTest == 1) %>%
select(-c(Demanda_uni_equil, toTest))
# Visualizando as primeiras linhas do dataset.
head(test)
Agora podemos remover a variável data.
# Removendo o dataset data.
rm(data)
5.7 Criando função para gerar modelos com diferentes valores de parametrização baseados no algoritmo XGboost¶
Bom, optamos por utilizar o algoritmo XGboost para a criação do nosso modelo preditivo por apresentar uma boa performance para gerar os scores da métrica de avaliação a ser utilizada e por ser consideravelmente mais rápido quando comparado a outros algoritmos.
Como não sabemos quais valores utilizar para sua configuração, criaremos uma função que gere diferentes modelos com diferentes ajustes e selecionaremos aquele que obtiver o melhor desempenho para os dados de teste.
Iremos avaliar a performance dos modelos a serem criados com base nos dados de treino e por isso deveremos ter atenção ao overfitting quando formos selecionar qual deverá ser utilizado para fazer as previsões dos dados de teste.
Adotando esta estratégia conseguiremos extrair o melhor que este algoritmo pode nos oferecer.
# Definindo uma função para gerar diferentes modelos com diferentes valores de parametrização baseados no algoritmo XGboost.
getBetterXGboostParameters <- function(data, label, maxDepth = 13, nEta = 0.2, nRounds = 86, subsample = 0.85,
colsample = 0.7, statusPrint = F) {
# Criando o dataframe para salvar os resultados dos modelos.
featuresXGboost <- data.frame()
# Define uma varíavel auxiliar para permitir o acompanhamento do progresso na avaliação dos modelos criados.
count <- 0
# Define o número total de modelos a serem criados.
total <- length(maxDepth) * length(nEta) * length(nRounds) * length(subsample) * length(colsample)
# Convertendo os dados das variáveis do dataset para o tipo DMatrix (uma matriz densa).
dTrain <- xgb.DMatrix(
data = data, # Define as variáveis preditoras.
label = label # Define a variável a ser prevista.
)
for(m in maxDepth) {
for(e in nEta) {
for(r in nRounds) {
for(s in subsample) {
for(c in colsample) {
# Define um seed para permitir que o mesmo resultado do experimento seja reproduzível.
set.seed(100)
# Criando o modelo baseado no algoritmo XGboost.
model_xgb <- xgb.train(
params = list(
objective = "reg:linear", # Define que o modelo deve ser baseado em uma regressão
# logistica linear.
booster = "gbtree", # Definindo o booster a ser utilizado.
eta = e, # Define a taxa de aprendizado do modelo.
max_depth = m, # Define o tamanho máximo da árvore.
subsample = s, # Define a proporção de subamostra da instância de treinamento.
colsample_bytree = c # Define a proporção da subamostra de colunas ao construir cada árvore.
),
data = dTrain, # Define as variáveis preditoras e a variável a ser prevista.
feval = rmsle, # Define a função de avaliação a ser utilizada.
nrounds = r, # Define o número de iterações que o algoritmo deve executar.
verbose = F, # Define a exibição da queda da taxa de erro durante o treinamento.
maximize = FALSE, # Define que a pontuação da avaliação deve ser minimizada.
nthread = 16 # Define o número de threads que devem ser usadas.
# Quanto maior for esse número, mais rápido será o treinamento.
)
# Realizando as previsões com o modelo baseado no algoritmo XGboost.
pred <- predict(model_xgb, data)
# Armazena os parâmetros utilizados para criação do modelo e o score da métrica RMSLE obtido
# em um dataframe.
featuresXGboost <- rbind(featuresXGboost, data.frame(
maxDepth = m,
eta = e,
nRounds = r,
s = s,
c = c,
rmsle = rmsle(label, pred)
))
# Incrementa o número de modelos avaliados.
count <- count + 1
# Imprime a porcetagem de progresso do treinamento e o melhor score da métrica RMSLE já alcançado.
print(paste(100 * count / total, '%, best rmsle: ', min(featuresXGboost$rmsle)))
# Salvando dataframe com os resultados gerados em um arquivo .csv.
write.csv(
x = featuresXGboost, # Determinando o conjunto de dados a ser
# salvo.
file = "/content/outputs/featuresXGboost.csv", # Define o nome com o qual o conjunto de
# dados deve ser salvo.
row.names = FALSE # Indica que o nome das linhas não deve ser
# gravado no arquivo a ser salvo.
)
}
}
}
}
}
# Retorna o dataframe com os resultados obtidos pelo treinamento de cada modelo.
featuresXGboost
}
5.8 Criando modelo XGboost¶
O algoritmo XGboost tem a capacidade de lidar com valores NA e por isso não vamos transformar os dados para tratar estas aparições.
Dito isto, podemos criar nosso modelo.
# Gerando diferentes modelos baseados no algoritmo XGboost e determinando os scores para a métrica RMSLE de cada um.
featuresXGboost <- getBetterXGboostParameters(
data = as.matrix(train %>% select(- Demanda_uni_equil)),
label = train$Demanda_uni_equil,
maxDepth = 12:14,
nEta = 0.2,
nRounds = 85:87,
subsample = 0.85,
colsample = 0.7,
statusPrint = F
)
# Salvando dataframe com os resultados gerados em um arquivo .csv.
fwrite(featuresXGboost, '/content/outputs/featuresXGboost.csv')
Caso deseje pular a execução do bloco de código anterior, basta carregar os resultados já processados que estão salvos no arquivo CSV abaixo:
# Carregando dataframe com os resultados obtidos para cada modelo XGboost criado.
featuresXGboost <- fread('/content/outputs/featuresXGboost.csv')
Imprimiremos os registros dos modelos criados.
# Visualizando dataframe com os resultados obtidos no treinamento.
featuresXGboost
Após utilizar cada uma das configurações acima para realizar as previsões dos dados de teste, observamos que aquela que obteve o melhor resultado está descrita na linha 5.
Os modelos registrados após está linha apresentam desempenhos inferiores pois começam a apresentar overfitting.
# Visualizando a melhor configuração para realizar as previsões dos dados de teste.
bestXGboost <- featuresXGboost[5, ]
bestXGboost
# Convertendo os dados das variáveis do dataset para o tipo DMatrix (uma matriz densa).
dTrain <- xgb.DMatrix(
data = as.matrix(train %>% select(- Demanda_uni_equil)), # Define as variáveis preditoras.
label = train$Demanda_uni_equil # Define a variável a ser prevista.
)
# Define um seed para permitir que o mesmo resultado do experimento seja reproduzível.
set.seed(100)
# Criando o modelo baseado no algoritmo XGboost.
model_xgb <- xgb.train(
params = list(
objective = "reg:linear", # Define que o modelo deve ser baseado em uma regressão logistica linear.
booster = "gbtree", # Definindo o booster a ser utilizado.
eta = bestXGboost$eta, # Define a taxa de aprendizado do modelo.
max_depth = bestXGboost$maxDepth, # Define o tamanho máximo da árvore.
subsample = bestXGboost$s, # Define a proporção de subamostra da instância de treinamento.
colsample_bytree = bestXGboost$c # Define a proporção da subamostra de colunas ao construir cada árvore.
),
data = dTrain, # Define as variáveis preditoras e a variável a ser prevista.
feval = rmsle, # Define a função de avaliação a ser utilizada.
nrounds = bestXGboost$nRounds, # Define o número de iterações que o algoritmo deve executar.
verbose = T, # Define a exibição da queda da taxa de erro durante o treinamento.
print_every_n = 5, # Define o número de iterações que devem ocorrer para que a impressão
# da mensagem de avaliação seja efetuada.
maximize = FALSE, # Define que a pontuação da avaliação deve ser minimizada.
nthread = 16 # Define o número de threads que devem ser usadas. Quanto maior for esse
# número, mais rápido será o treinamento.
)
# Realizando as previsões com o modelo baseado no algoritmo XGboost.
pred <- predict(model_xgb, as.matrix(test %>% select(- id)))
# Convertendo os resultados previsto para a escala original da variável alvo (exp(Demanda_uni_equil) - 1).
pred <- expm1(pred)
# Transformando qualquer previsão negativa em um valor nulo.
pred[pred < 0] <- 0
# Salvando os resultados gerados em um arquivo CSV.
write.csv(
x = data.frame(id = as.integer(test$id), Demanda_uni_equil = pred), # Determinando o conjunto de dados a ser
# salvo.
file = "/content/outputs/results.csv", # Define o nome com o qual o conjunto de
# dados deve ser salvo.
row.names = FALSE # Indica que o nome das linhas não deve
# ser gravado no arquivo a ser salvo.
)
Por fim, obtemos as seguintes pontuações com a métrica RMSLE no Kaggle ao submeter as previsões geradas pelo nosso modelo para os dados de teste:
- Pontuação privada: 0.46821
- Pontuação pública: 0.45352
Entre em contato comigo!¶
Caso tenha alguma dúvida, sugestão ou apenas queira trocar uma ideia sobre este projeto, não hesite em entrar em contato comigo!
