Análise Preditiva da Taxa de Utilização Casual das Bicicletas Alugadas
Laion Lima Boaventura
25 de outubro de 2019
Escopo
A Nova Mobilidade tem mudado a maneira como nos deslocamos pelas cidades. Estamos falando de serviços de transporte A Nova Mobilidade tem mudado a maneira como nos deslocamos pelas cidades. Estamos falando de serviços de transporte sob demanda (como Uber, 99 e Cabify), veículos elétricos e aplicativos de planejamento de rotas, mas, também, das bicicletas. Novas empresas apostam na bicicleta como meio para atender às necessidades de transporte nas áreas urbanas e os programas de compartilhamento, públicos ou privados, mudam a maneira como enxergamos as bikes: de opção de lazer a meio de transporte e até mesmo instrumento de trabalho.
Os programas de bicicleta compartilhada evoluíram e hoje são meio de transporte e até instrumento de trabalho. Atualmente, existem mais de 500 programas de compartilhamento de bicicletas em todo o mundo, com serviço para os cidadãos locais. Buscando mais adeptos, alguns programas permitem o chamado “uso casual”. Esses usuários, podem utilizar as bicicletas sem a necessidade de se registrar à empresa. Neste sentido, é de grande interesse prever a taxa de utilização casual dessas bicicletas, uma vez que o aumento de demanda, poderia refletir em uma baixa qualidade no atendimento dos usuários registrados. De acordo com especialistas da própria empresa, nos períodos em que o percentual de clientes casuais supera 20% se faz necessário a aquisição de novas bicicletas a serem ofertadas, evitando assim filas e descontentamento dos clientes registrados.
Neste case, utilizei dados reais de uma empresa privada da Inglaterra que forcene o serviço de aluguel de bicicleta em diferentes cidades. Para uma previsão de qualidade, conceitos aprendizagem estatística de máquina, implementados à diferentes modelos, bem como noções descritivas e exploratória dos dados, são aplicados.
Toda parte computacional, da implementação à apresentação foi feita no Software R, e os códigos estão disponíveis para visualização no corpo do próprio trabalho.
Análise dos Dados
library(readr)
dados = read_csv("dados1.csv")
library(lubridate)o banco de dados original contava com as seguintes variáveis:
- instante (índice de registro do usuário)
- dteday (data em ano, mês e dia)
- season (estação do ano)
- holiday (indicador se o dia era feriado o não)
- weekday (indicador se o dia era final de semana ou não)
- workday (indicador se o dia dia útil ou não)
- weathersit (indicador da condição climática)
- temp (temperatura real em graus Celsius)
- atemp (sensação térmica em graus Celsius)
- hum (umidade do ar)
- windspeed (velocidade do vento)
- casual (contagem dos usuários casuais)
- registred (contagem dos usuários registrados)
- cnt (contagem das bicicletas requeridas = casual + registred)
E, as seguintes tomadas de decisão foram realizadas para tratamento dos dados:
# Retirei a Variável instante (já que atuava como um índice particular de cada medição)
dados = dados[,-1]
# Transformei a Variável season em fator
dados$season = as.factor(dados$season)
# Transformei a Variável yr em fator
dados$yr = as.factor(dados$yr)
# Criei a variável mnth a partir da Transformação da Variável dteday em meses
dados$mnth= month(as.POSIXlt(dados$dteday, format="%Y/%m/%d"), label = TRUE, abbr = FALSE)
# Transformei a Variável holiday em fator
dados$holiday = as.factor(dados$holiday)
# Criei a variável mnth a partir da Transformação da Variável dteday em dias da semana
dados$weekday= as.factor(weekdays(dados$dteday))
# Retirei A variavel dteda, já que tinha extraído toda informação necessária
dados = dados[,-1]
# Transformei a Variável workingday em fator
dados$workingday = as.factor(dados$workingday)
# Transformei a Variável weathersit em fator
dados$weathersit = as.factor(dados$weathersit)
# Transformei as Variáveis casual em taxa (%)
dados$casual = dados$casual/dados$cnt
## Como a variável registered transmite informações complementares a variável casual, retirei-a
dados = dados[,-13]Resultando num banco de dados com 731 observações e 13 variáveis, sendo 12 explicativas. Podemos observar o resumo do banco da seguinte forma:
library(skimr)
skim_with(numeric = list(hist = NULL))
skim(dados)## Skim summary statistics
## n obs: 731
## n variables: 13
##
## -- Variable type:factor -------------------------------------------------------------------------------------------------
## variable missing complete n n_unique
## holiday 0 731 731 2
## mnth 0 731 731 12
## season 0 731 731 4
## weathersit 0 731 731 3
## weekday 0 731 731 7
## workingday 0 731 731 2
## yr 0 731 731 2
## top_counts ordered
## 0: 710, 1: 21, NA: 0 FALSE
## jan: 62, mar: 62, mai: 62, jul: 62 TRUE
## 3: 188, 2: 184, 1: 181, 4: 178 FALSE
## 1: 463, 2: 247, 3: 21, NA: 0 FALSE
## dom: 105, sáb: 105, seg: 105, qua: 104 FALSE
## 1: 500, 0: 231, NA: 0 FALSE
## 1: 366, 0: 365, NA: 0 FALSE
##
## -- Variable type:numeric ------------------------------------------------------------------------------------------------
## variable missing complete n mean sd p0 p25 p50
## atemp 0 731 731 0.47 0.16 0.079 0.34 0.49
## casual 0 731 731 0.18 0.11 0.013 0.099 0.15
## cnt 0 731 731 4504.35 1937.21 22 3152 4548
## hum 0 731 731 0.63 0.14 0 0.52 0.63
## temp 0 731 731 0.5 0.18 0.059 0.34 0.5
## windspeed 0 731 731 0.19 0.077 0.022 0.13 0.18
## p75 p100
## 0.61 0.84
## 0.21 0.51
## 5956 8714
## 0.73 0.97
## 0.66 0.86
## 0.23 0.51Análise univariada
Para uma melhor exploração do banco, é interessante realizar uma análise descritiva univariada dos dados disponíveis. Tentando assim, entender previamente o comportamento de cada variável mediante uma outra variável de referência. Para isso, escolhi estudar a relação das variáveis explicativas com a variável cnt. Isto porque, como ela é derivada do número de usuários casuais, e espera-se que seu comportamento em relação as demais variáveis, reflita o comportamento das demais variáveis a propria variável casual. E como a variável casual sofreu alterações, passando a ser taxa, utilizar cnt garante uma melhor interpretabilidade dos resultados gráficos.
Variáveis Categóricas:
Aqui, construí os Gráficos Boxplot, a fim de identificar diferenças no comportamento do interesse do usuário, conforme os níveis das variáveis categóricas explicativas presentes no banco.
boxplot(dados$cnt~dados$mnth, main = 'Contagem de Bicicletas vs Meses ', col='lightblue', xaxt = "n")
text(x = 1:length(levels(dados$mnth)),
labels = levels(dados$mnth),
xpd = NA,
srt = 35,
cex = 0.8)
boxplot(dados$cnt~dados$season, ylim= c(0,9800), main = 'Contagem de Bicicletas vs Estacoes ', col=c('lightblue','tomato','orange','green'))
legend( 'topleft',pch=19,c("Primavera ","Verao ","Outono","Inverno"),col=c('lightblue','tomato','orange','green'),bty="n",cex=0.8)
boxplot(dados$cnt~dados$yr, main = 'Contagem de Bicicletas nos anos de 2011 e 2012', col=c('lightblue','tomato'))
legend( 'topleft',pch=19,c("2011","2012"),col=c('lightblue','tomato'),bty="n",cex=0.8)
boxplot(dados$cnt~dados$holiday,ylim = c(0,9800), main = 'Contagem de Bicicletas vs Feriado', col=c('lightblue','tomato'))
legend( 'topleft',pch=19,c("Nao Feriado","Feriado"),col=c('lightblue','tomato'),bty="n",cex=0.8)
f = as.numeric(dados$weekday)
f = as.factor(f)
f = ordered(interaction(f),levels = c("1","5",
"7","2","3","6","4"))
boxplot(dados$cnt~f, xaxt = "n", main = 'Contagem de Bicicletas vs Dias da Semana', col='lightblue', ylim = c(0,12000))
text(x = 1:length(levels(f)),
labels = c("domingo","segunda-feira",
"terca-feira","quarta-feira","quinta-feira","sexta-feira","sabado"),
xpd = NA,
srt = 35,
cex = 1)
boxplot(dados$cnt~dados$workingday,ylim = c(0,9800), main = 'Contagem de Bicicletas vs Dia Util', col=c('lightblue','tomato'))
legend( 'topleft',pch=19,c("Feriado ou Final de Semana","Dia Util"),col=c('lightblue','tomato'),bty="n",cex=0.8)
boxplot(dados$cnt~dados$weathersit, ylim= c(0,12000), main = 'Contagem de Bicicletas vs Condicoes Climaticas', col=c('lightblue','tomato','orange'))
legend( 'topleft',pch=19,c("[Ensolarado ; Parcialmente Nublado]","]Parcialmente Nublado ; Nublado]","]Nublado ; Chovoso]"),col=c('lightblue','tomato','orange'),bty="n",cex=0.8)
Pelos gráficos ilustrados podemos analisar o período de férias de verão na europa apresenta um maior volume de alugueis que nos demais períodos. Podemos concluir também que existe uma aparente diferença no volume dos alugueis de acordo com as condições climáticas do dia; como esperado, em dias chuvosos o número de bicicletas alugadas é menor que em dias com baixa probabilidade de chuva. Outra análise é que, de fato, houve um crescimento geral no número de bicicletas alugadas entre os anos de 2011 e 2012, evidenciando uma crescente na demanda por esse tipo de produto.
Variáveis Contínuas:
Aqui, analisei a relação das variáveis explicativas pela matriz de correlação. O que permitiu a observação acerca da relação linear entre as covariáveis explicativas e a variável resposta, bem como a relação linear entre as próprias variáveis explicativas.
Assim elipses ilustradas no gráfico abaixo indicam tanto o sentido quanto a itensidade e a magnitude dessas relações lineares. Como pode-se observar, a única relação que merece um certo tipo de atenção é a que ocorre entre as variáveis temp e atemp. A alta relação linear entre elas é forte indicador de colinearidade, o que torna a matriz hessiana não inversível, e consequentemente, prejudica o uso de classificadores lineares fudamentados em conceitos verossimilhancistas.
library(corrplot)
corrplot.mixed(cor(dados[,c(8:13)]),upper = "ellipse")
Modelagem
Nesta etapa do trabalho, mostrarei o passo a passo do desenvolvimento da modelagem, desde a criação de funções auxiliares, até a apresentação do resultado sobre a perfonmance de cada modelo, passando pela seleção e escolha mediante processo de turning em alguns modelos.
Antes disso, é importante verificar o comportamento univariado da variável de interesse (a variável resposta), ou seja, analisar se a ocorrência de interesse (a taxa de clientes usuais superior à 20%) é um caso raro o suficiente para que a variável resposta seja considerada como desbalanceada.
dados$casual = ifelse(dados$casual >0.2 ,1,0)
dados$casual = as.factor(dados$casual)
boxplot(dados$cnt~dados$casual,ylim = c(0,9800), main = 'Contagem de Bicicletas vs Taxa Casual', col=c('lightblue','tomato'))
legend( 'topleft',pch=19,c("de 0 a 20% de usuarios casuais","Mais de 20% de usuarios casuais"),col=c('lightblue','tomato'),bty="n",cex=0.8)
Como observado pelo boxplot ilustrado acima, a ocorrência de dias em que a taxa de clientes usuais é superior à 20% não pode ser considerado um caso raro, e assim, não faz-se necessário utilizar uma métrica de desempenho mais robusta para análise dos resultados. Desta forma, nos passos a seguir, referentes a seleção dos modelos, utiliza-se a acurácia para análise da performance preditiva de cada um deles.
Funções Auxiliares
Aqui, apresento funções auxiliares para aplicação dos modelos. A função medida, além da acurácia, trás a formulação de outras métricas de desempenho. Isto, possibilita a quem de interesse for, repetir este trabalho alterando medida de desenpenho utilizada.
varResp = 'casual'
Medidas = function(mc){
vp = mc[2,2]
vn = mc[1,1]
fp = mc[2,1]
fn = mc[1,2]
acc = (vp + vn)/sum(mc)
sen = vp/(vp + fn) #Recall/Sensibilidade
esp = vn/(vn + fp) #Especificidade
vpp = vp/(vp + fp) #Precision/Valor preditivo positivo
vpn = vn/(vn + fn) #Valor preditivo negativo
mcc = (vp*vn - fp*fn)/(sqrt((vp + fp)*(vp + fn)*(vn + fp)*(vn + fn))) #Coeficiente de correlacao de Matthews
mcc = ifelse(is.na(mcc), 0, mcc)
f1 = 2*sen*vpp/(sen + vpp)
return(list(Acc = acc, Sens = sen, Espec = esp, VPP = vpp, VPN = vpn, MCC = mcc, F1Score = f1))
}
AuxSaida = function(pred, teste, formula, alg){
ind = which(colnames(teste) == as.character(formula)[2])
mc = table(pred, teste[, ind][[1]])
rs = Medidas(mc)
rs[['Classificador']] = alg
return(do.call(data.frame, rs))
}
AuxSaida.bayes = function(pred, teste, formula, alg){
ind = which(colnames(teste) == as.character(formula)[2])
mc = table(pred, teste[, ind])
rs = Medidas(mc)
rs[['Classificador']] = alg
return(do.call(data.frame, rs))
}
formula = as.formula(paste(varResp, '~ .'))
FormFac = function(formula){
aux = as.character(formula)[2]
return(as.formula(paste('factor(',aux, ')', '~ .')))
}Funções de Data Splitting
Aqui, apresento funções referentes ao processo de data splitting. Neste trabalho, fiz a validação baseando-se na validação cruzada k-fold repetida. Para um k=10, e 10 repetições.
Contudo, assim como na seção anterior, trago outras alternativas a quem for de interesse a replicar esse trabalho. Desta forma, além da própria validação k-fold, sem repetições, uma outra possibilidade é o uso do Holdout, implementado através da função denominada Holdout, em que o parâmetro p, indica o percentual utilizado na amostra de treino.
###### DATA SPLITTING
Holdout = function(dados, p = 0.7){
n_treino = ceiling(dim(dados)[1]*p)
ind = c(rep('treino', n_treino), rep('teste', dim(dados)[1] - n_treino) )
ind = sample(ind)
treino = dados[ind == 'treino', ]
teste = dados[ind == 'teste' , ]
return(list(treino = treino, teste = teste))
}
ValidCruzada = function(Bayesiano = 'n', cls, dados, formula, kfolds = 10, ...){
library(dismo)
library(bnlearn)
dados.discretizados = discretize(dados,breaks = 3)
if (Bayesiano == 's') {dados = dados.discretizados}
id = kfold(1:nrow(dados), kfolds)
vcs = data.frame()
for(i in 1:kfolds){
treino = dados[id != i, ]
teste = dados[id == i, ]
kcls = cls(treino, teste, formula)
vcs = rbind(vcs, kcls)
}
vcs = vcs[, c(ncol(vcs), 1:(ncol(vcs) - 1))]
avg = aggregate(. ~ Classificador, data = vcs, FUN = mean)
std = aggregate(. ~ Classificador, data = vcs, FUN = sd)
return(list(Media = avg, Desvio = avg, Modelos = vcs))
}
ValCruzRep = function(Bayesiano = 'n',cls, dados, formula, kfolds = 10, reps = 10,...){
x = data.frame()
for(i in 1:reps) x = rbind(x, ValidCruzada(Bayesiano = Bayesiano,cls, dados, formula, kfolds)$Media)
avg = aggregate(. ~ Classificador, data = x, FUN = mean)
std = aggregate(. ~ Classificador, data = x, FUN = sd)
return(list(Media = avg, Desvio = std, Modelos = x))
}Funções dos Modelos
Nesta seção, apresento as funções de cada um dos modelos estudados: Regressão Logística, KNN, LDA, Árvore de Decisão, Naive Bayes e KDB. Em que, para os modelos KNN e KDB apresento as funções acerca do processo de turning sobre seus parâmetros, o número de vizinhos mais próximos e o número de dependências entre as covariáveis explicativas, respectivamente.
########## Log?stica
RegLog = function(treino, teste, formula){
cls = glm(formula, data = treino, family = 'binomial')
pred = predict(cls, newdata = teste, type = 'response')
pred = ifelse(pred > 0.5, 1, 0)
pred = as.numeric(pred)
pred = as.factor(pred)
return(pred)
}
medRegLog = function(treino, teste, formula){
alg = 'Regressao Logistica'
pred = RegLog(treino, teste, formula)
return(AuxSaida(pred, teste, formula, alg))
}
########## KNN
KNN = function(treino, teste, formula, K = 7){
library(kknn)
cls = kknn(FormFac(formula), train = treino, test = teste, k = K)
pred = cls$fitted.values
pred = as.numeric(pred) - 1
return(pred)
}
medKNN = function(treino, teste, formula, K = 7){
alg = 'KNN'
pred = KNN(treino, teste, formula, K)
return(AuxSaida(pred, teste, formula, alg))
}
# Turning de K vizinhos do KNN, via Validação Cruzada K=10 folds
KNNmelhorK = function(dados, formula, inicio = 1, fim = 31, metrica = 'Acc', ...){
knns = data.frame()
s = seq(inicio, fim, by = 2)
for(k in s) knns = rbind(knns, ValidCruzada(Bayesiano = 'n',medKNN, dados, formula, K = k, ...)$Media)
ordem = order(knns[metrica])
melhorK = s[tail(ordem, 1)]
ord = knns[ordem, ]
plot(s, knns[[metrica]], pch = 17, xlab = 'Valores de K (Vizinhos)', main = paste('Melhor K - metrica:', metrica),
ylab = metrica, xaxt = "n", col ='blue')
axis(1, at = s)
segments(x0 = melhorK, y0 = 0, x1 = melhorK, y1 = as.numeric(tail(ord[metrica], 1)),
lty = 'dotted', col = 'gray50', lwd = 2)
return(melhorK)
}
########## LDA
LDA = function(treino, teste, formula){
library(MASS)
cls = lda(formula, data = treino)
pred = predict(cls, newdata = teste)$class
pred = as.numeric(pred) - 1
return(pred)
}
medLDA = function(treino, teste, formula){
alg = 'Analise Discriminante Linear'
pred = LDA(treino, teste, formula)
return(AuxSaida(pred, teste, formula, alg))
}
########## ARVORE
Arvore = function(treino, teste, formula){
library(rpart)
cls = rpart(formula, data = treino, method = 'class')
pred = predict(cls, newdata = teste, type = 'class')
pred = as.numeric(pred) - 1
return(pred)
}
medArvore = function(treino, teste, formula){
alg = 'Arvore de Decisao'
pred = Arvore(treino, teste, formula)
return(AuxSaida(pred, teste, formula, alg))
}
########## Naive Bayes
NB = function(treino, teste,varResp, formula){
library(bnclassify)
cls = bnc('nb', varResp, treino, smooth = 1)
pred = predict(cls, teste)
return(pred)
}
medNB = function(treino, teste,formula){
alg = 'Naive Bayes'
pred = NB(treino, teste,varResp, formula)
return(AuxSaida.bayes(pred, teste, formula, alg))
}
########## KDB
KDB = function(treino, teste,varResp, formula,K=2){
library(bnclassify)
cls = kdb(varResp,treino, k=10, kdbk = K, epsilon = 0.01, smooth = 0,cache_reset = NULL)
cls = lp(cls, treino, smooth = 1)
pred = predict(cls,teste,prob = FALSE)
return(pred)
}
medKDB = function(treino, teste,formula,K=2){
alg = paste('KDB')
pred = KDB(treino, teste,varResp, formula,K=K)
return(AuxSaida.bayes(pred, teste, formula, alg))
}
# Turning de K dependentes do KDB, via Validação Cruzada K=10 folds
KDBmelhorK = function(dados, formula, inicio = 1, fim = (dim(dados)[2]-2), metrica = 'Acc', ...){
knns = data.frame()
s = seq(inicio, fim, by = 1)
for(k in s) knns = rbind(knns, ValidCruzada(Bayesiano = 's',cls = medKDB,dados, formula, kfolds = 10,K=k)$Media)
ordem = order(knns[metrica])
melhorK = s[tail(ordem, 1)]
ord = knns[ordem, ]
plot(s, knns[[metrica]], pch = 17, xlab = 'K- Dependentes', main = paste('Melhor K - metrica:', metrica),
ylab = metrica, xaxt = "n", col = 'blue')
axis(1, at = s)
segments(x0 = melhorK, y0 = 0, x1 = melhorK, y1 = as.numeric(tail(ord[metrica], 1)),
lty = 'dotted', col = 'gray50', lwd = 2)
return(melhorK)
}Turning dos KNN e KDB
Aqui, apliquei as funções de tunagem apresentadas acima, a fim de escolher o melhor modelo entre os KNN, variando o número de vizinhos mais próximos a serem considerados; bem como o modelo KDB, variando o número de k-dependentes considerados. Lembre que esses parâmetros representam a complexidades desses modelos, e portanto merece atenção quanto ao tradeoff viés-variância. Isto evita possíveis situações de over(under)fitting.
#TURNING KNN
KNNmelhorK(dados, formula, 1 , (dim(dados)[1]-1), 'Acc', kfolds = 10)
## [1] 391#TURNING KDB
KDBmelhorK(dados, formula, 1 ,(dim(dados)[2]-2), 'Acc', kfolds = 10)
## [1] 1Desempenho dos Modelos
Por fim, aplicando a função ValCruzRep, ilustrada anteriormente, construímos, com base na acurácia, a tabela de desempenho dos modelos estudados. Observe que o desempenho em todos os modelos foi muito bom, o que evidencia um tratamento adequado dos dados.
Além disso, ao analisar a tabela abaixo, destaco que, entre os modelos estudados, o de melhor performance preditiva foi a Regressão Logística, o que indica a possibilidade de traçar uma fronteira de decisão linear, facilitando a interpretabilidade dos resultados.
tabela.desempenho = matrix(NA,6,2)
colnames(tabela.desempenho) = c('Modelo', 'Acuracia')
saida = ValCruzRep(Bayesiano = 'n',cls = medRegLog,dados, formula, kfolds = 10)
tabela.desempenho[1,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 'n',cls = medLDA,dados, formula, kfolds = 10)
tabela.desempenho[2,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 'n',cls = medKNN,dados, formula, kfolds = 10,K= 85)
tabela.desempenho[3,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 'n',cls = medArvore,dados, formula, kfolds = 10)
tabela.desempenho[4,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 's',cls = medNB,dados, formula, kfolds = 10)
tabela.desempenho[5,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
saida = ValCruzRep(Bayesiano = 's',cls = medKDB,dados, formula, kfolds = 10, K= 8)
tabela.desempenho[6,]=c(levels(saida$Media$Classificador),round(saida$Media[[2]],3))
tabela.desempenho= tabela.desempenho[order(tabela.desempenho[,2]),]
library('formattable')
tabela.desempenho = data.frame(tabela.desempenho)
formattable(tabela.desempenho, list(Acuracia =color_tile("white", "orange")))| Modelo | Acuracia |
|---|---|
| Naive Bayes | 0.904 |
| Analise Discriminante Linear | 0.905 |
| KNN | 0.916 |
| KDB | 0.919 |
| Arvore de Decisao | 0.922 |
| Regressao Logistica | 0.936 |