Como escolher o ponto de corte em um classificador binário: Uma abordagem pela Teoria da Decisão¶
Autor: Felipe Maia Polo
Site pessoal: https://felipemaiapolo.github.io/
Linkedin: https://www.linkedin.com/in/felipemaiapolo/
E-mail: felipemaiapolo@gmail.com
import numpy as np
import pandas as pd
from pandas import DataFrame
import random
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MaxAbsScaler
from sklearn import metrics
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from sklearn.metrics import precision_recall_curve
from sklearn.metrics import plot_roc_curve
from sklearn.metrics import roc_curve, auc
from sklearn.metrics import roc_auc_score
from sklearn.linear_model import LogisticRegression
import matplotlib
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import seaborn as sns
from IPython.display import Image
from IPython.core.display import HTML
Pipeline para se treinar um modelo de Machine Learning em trabalhar dados reais¶
Limpar dados e criar variáveis dummy a partir das variáveis qualitativas.
- Quando transformar uma variável categórica (com $k$ categorias) em $k$ variáveis binárias, lembrar de excluir uma das categorias;
- Se houver classes com muito poucos dados, podemos fundir as categorias;
- Se tivermos valores missing (olha que legal), podemos considerá-los como uma nova classe, a classe dos missings!
Aqui é importante ter um pipeline claro de limpeza e que pode ser replicado para quando novos dados chegarem.
Entender bem o que cada variável quer dizer. Qual é o meu $Y$ e quais variáveis são meu $X$, ou seja, quais fariam sentido colocar no meu modelo?
Dividir seus dados em treino (aprox. 80\%) e teste (aprox. 20\%), de preferência garantindo que a distribuição de $Y$ seja parecida nos dois conjuntos de dados. Teremos (X_train, y_train) e (X_test, y_test);
Verificar como a variável $Y$ está distribuída no conjunto de treino;
Das variáveis $X$, quais variáveis são quantitativas e quais são qualitativas (categóricas)?
Como as variáveis $X$ estão distribuídas no conjunto de treino? Temos outliers? As variáveis têm boa variabilidade?
Validação: se o modelo escolhido tiver hiperparâmetros que vocễ queira ajustar (e.g. ponto de corte, força da regularização etc), destinar parte dos dados de treino para validação, afim de escolher os melhores valores para os hiperparâmetros. Se você tiver a intenção de normalizar/padronizar suas variáveis quantitativas, é importante que os parâmetros da normalização/padronização não sejam calculados na base de validação!
Em relação a X_train, normalizar/padronizar as variáveis quantitativas.
Treinar modelo final em (X_train, y_train), ou seja, no conjunto de teste inteiro;
Normalização/padronização: a normalização/padronização deve ser feita com base nos parâmetros do conjunto de treino;
Avaliar performance do modelo em (X_test, y_test): acurácia, matriz de confusão, recall, precision etc;
Exemplo (SUS)¶
A base de dados utilizada foi retirada do Kaggle no link: https://www.kaggle.com/joniarroba/noshowappointments (está no meu Drive também https://drive.google.com/open?id=13MIRaMI1QDCePEmLW59Fkaosvvuq4akX). Esta é uma base brasileira de consultas médicas que foram agendadas no SUS de Vitória-ES. Cada observação da base de dados traz consigo dados do paciente, do agendamento e uma variável indicadora (0 ou 1) que nos diz se o paciente faltou ou não à consulta. Esta base é interessante, pois nos traz um caso real que nos permite entender melhor porque as pessoas faltam às consultas e nos ajuda a antecipar situações de falta para melhor gerenciamento do sistema como um todo. Apesar do nosso exemplo ser um caso específico da relação prestador de serviços e clientes, a metodologia que discutiremos a seguir pode ser utilizada em muitos outros casos. A descrição das variáveis de acordo com o Kaggle é:
PatientId: Id dos pacientes
Gender: gênero do paciente
DataMarcacaoConsulta: dia consulta
DataAgendamento: dia do agendamento
Age: idade do paciente
Neighbourhood: vizinhança onde a consulta foi marcada
Scholarship: bolsa família
Hipertension: tem ou não hipertensão
Diabetes: tem ou não diabetes
Alcoholism: tem ou não alcoolismo
Handcap: tem ou não debilidade física
SMSreceived: se foram ou não enviadas mensagens para o paciente, lembrando-o da consulta
No-show: se foi ou não na consulta
Vamos dar uma explorada/preparada na base para começarmos os trabalhos!
Conectando ao Google Drive:
from google.colab import drive
drive.mount('/content/drive')
Definindo caminho para os dados:
path = "/content/drive/My Drive/Data Science Course - Felipe/Materiais/dados/"
Abrindo os dados:
data=pd.read_csv(path+"medical.csv")
Visualizando dados:
data.head()
Vamos dar uma arrumada em algumas variáveis:
data.ScheduledDay = data.ScheduledDay.apply(np.datetime64) #Formatando para data
data.AppointmentDay = data.AppointmentDay.apply(np.datetime64) #Formatando para data
data['WaitTime'] = (data.AppointmentDay - data.ScheduledDay).dt.days+1 #Dias de espera
data['NoShow']= data['No-show'].apply(lambda x: 1 if x =="Yes" else 0) #1=não foi à consulta e 0=foi
data['Gender']= data['Gender'].apply(lambda x: 0 if x =="M" else 1) #1=Feminino e 0=Masculino
data['Dia'] = data.AppointmentDay.dt.weekday #0=segunda-feira e 6=domingo
Tirando algumas estatísticas descritivas:
data.describe()
Percebemos que há três coisas estranhas com nossa base: (i) há pessoas com idades negativas, (ii) há pessoas que tem a variável "Handcap" (debilidade física) como sendo maior que 1 (na teoria era para ser 0 ou 1) e (iii) pacientes com tempo de espera negativo em dias. Vamos contar quantos indivíduos carregam esses problemas para então transformar seus valores em "missing":
#Primeiramente para a idade:
sum(data['Age']<0)
#Handcap
sum(data['Handcap']>1)
#Tempo de espera
sum(data['WaitTime']<0)
Mandando valores estranhos para missing:
data['Age']=np.where(data['Age']<0,np.nan,data['Age'])
data['Handcap']=np.where(data['Handcap']>1,np.nan,data['Handcap'])
data['WaitTime']=np.where(data['WaitTime']<0,np.nan,data['WaitTime'])
Como queremos criar um modelo para prever quem irá comparecer ou não à consulta, não faz sentido mantermos aqueles pacientes que têm 'WaitTime'=0, pois não queremos criar um modelo para essas pessoas. Vamos manter somente aqueles que têm 'WaitTime'>0:
data=data[data['WaitTime']>0]
Transformando variáveis categóricas em dummies (variáveis binárias), começando pelo dia da semana:
d=pd.get_dummies(data["Dia"], dummy_na=False)
d=pd.DataFrame(d)
d.columns=["Seg", "Ter","Qua", "Qui","Sex","Sab"]
data=data.join(d)
Agora que já mexemos no que queríamos para deixar nossa base mais "trabalhável', vamos mexer com algumas tabelas que podem nos dar pistas do porquê as faltas. Nossa primeira hipótese é que pessoas com alguma debilidade física tendem a faltar mais. Vamos então olhar mais de perto a distribuição conjunta de faltas e debilidade física:
tab=pd.crosstab(data['Handcap'],data['NoShow'], normalize=True)
tab
É muito claro que a maioria das pessoas não tem debilidade alguma e isso acaba deixando a primeira linha da tabela muito densa, fazendo com que fique difícil avaliar a probabilidade de ir ou não à consulta condicionada ao fato de ter ou não debilidade física. Vamos então normalizar nas linhas, fazendo com que a soma por linha seja 1, para obtermos assim a probabilidade condicional que queremos:
tab=pd.crosstab(data['Handcap'],data['NoShow'], normalize= 'index')
tab
É possível ver uma diferença de mais ou menos 3,5p.p. na probabilidade de aparecer na consulta, o que invalida, a princípio, minha hipótese. Vamos agora checar como a variável gênero se relaciona com aparecer na consulta:
tab=pd.crosstab(data['Gender'],data['NoShow'], normalize= 'index')
tab
É possível ver que parece não haver nenhuma relação entre Gênero e presença na consulta. Vamos ver a variável de SMS:
tab=pd.crosstab(data['SMS_received'],data['NoShow'], normalize= 'index')
tab
Vemos que quem recebeu o SMS tem menor probabilidade de faltar na consulta, algo intuitivo. No entanto, essa diferença não é muito grande - vamos investigar o porquê:
tab=pd.crosstab(data['SMS_received'],data['WaitTime'], normalize='index')
sns.heatmap(tab, cmap="Greys")
plt.show()
tab
É possível ver que grande parte das pessoas que não receberam SMS marcou a consulta um ou dois dias antes da própria consulta. Possivelmente, o pessoal do SUS de Vitória considera que essas pessoas têm probabilidade muito baixa de faltar e por isso não enviam os avisos - isso pode, de fato, fazer com que subestimamos a importância do SMS.
Utilizaremos essa informação e criaremos uma variável nova que indica se a pessoa marcou a consulta 1 ou 2 dias antes, pois acreditamos que essa seja uma informação importante:
data['Recent']=None
data.loc[data['WaitTime']<=2, 'Recent']=1
data.loc[data['WaitTime']>2, 'Recent']=0
Vamos avaliar agora a relação entre a diferença de dias entre a consulta e o dia em que a consulta foi marcada:
tab=pd.crosstab(data['WaitTime'],data['NoShow'], normalize='index')
sns.heatmap(tab, cmap="Greys")
plt.show()
Algumas coisas são notáveis do gráfico acima: (i) quando a pessoa marca no dia, muito provavelmente ela vai (até porque pode ser que o registro de agendamento seja feito na hora), (ii) nos primeiros 20 dias, há uma clara tendência em diminuir a probabilidade de ir à consulta, (iii) depois fica difícil pegar um padrão, talvez porque o número de consultas diminua muito e isso impacte na confiabilidade da estimação. Vamos agora trabalhar com a variável de idade. Como a variável de idade tem muitos valores, vamos trabalhar com o Heatmap, para dar uma noção mais geral dos fatos de forma rápida:
tab=pd.crosstab(data['Age'],data['NoShow'], normalize= 'index')
sns.heatmap(tab, cmap="Greys")
plt.show()
No gráfico acima, é possível ver que a idade pode ter alguma influência na probabilidade de aparecer ou não na consulta. O que nos parece é que até certa idade, os anos estão relacionados positivamente com não aparecer na consulta e depois a relação passa a ser positiva - parece que temos uma relação quadrática. Na tentativa de pegar alguma relação quadrática entre idade e a probabilidade de comparecimento na consulta, vamos incluir idade 2 na equação:
data['Age2']=data['Age']**2
Até o momento avaliamos de forma descritiva o impacto de cada uma das variáveis na probabilidade de comparecer nas consultas. Agora vamos de fato treinar nosso modelo de Machine Learning:
X=data[['WaitTime', 'Recent', 'Gender', 'Age', 'Age2', 'Scholarship', 'Hipertension', 'Diabetes', 'Alcoholism', 'Handcap', 'SMS_received', 'Ter', 'Qua', 'Qui', 'Sex', 'Sab']]
y=data["NoShow"]
Percebi que há algumas linhas em que há valores faltantes. Isso porque mudamos alguns valores para NA anteriormente. Vamos exlcuir essas linhas da nossa análise:
indices = np.where(np.isnan(X))[0]
X=X.drop(X.index[indices])
y=y.drop(y.index[indices])
Dividindo nossos dados nos conjuntos de teste e treino:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, stratify=y, random_state=27)
Para o processo de validação, vamos separar parte do conjunto de treino para validação:
X_train2, X_val, y_train2, y_val = train_test_split(X_train, y_train, test_size=.2, stratify=y_train, random_state=27)
Vamos utilizar o Max-Abs Scaler para normalizar os dados de validação:
scaler_val = MaxAbsScaler()
scaler_val.fit(X_train2) #dando o fit na base de treino
#transformando
X_train2=scaler_val.transform(X_train2)
X_val=scaler_val.transform(X_val)
Desenhando curvas de Recall e Precisão:
### Processo de validação ###
logreg = LogisticRegression(solver='liblinear').fit(X_train2, y_train2)
#Prevendo probabilidade (segunda coluna para y=1)
y_proba=logreg.predict_proba(X_val)
#plotando Precisão-Recall
precision, recall, thresholds = precision_recall_curve(y_val, y_proba[:,1])
plt.title("Precisão-Recall vs Ponto de Corte")
plt.plot(thresholds, precision[: -1], "b--", label="Precisão")
plt.plot(thresholds, recall[: -1], color='r',label="Recall")
plt.ylabel("Precisão, Recall")
plt.xlabel("Ponto de Corte")
plt.legend(loc="lower left")
plt.ylim([0,1])
Vamos fixar três valores de $c$ e ver os seus reports na base de validação:
c1=.2
c2=.25
c3=.3
y_pred1 = 1*(y_proba[:,1]>c1)
y_pred2 = 1*(y_proba[:,1]>c2)
y_pred3 = 1*(y_proba[:,1]>c3)
report1=classification_report(y_val, y_pred1)
report2=classification_report(y_val, y_pred2)
report3=classification_report(y_val, y_pred3)
print("########### c1 ###########\n",report1)
print("########### c2 ###########\n",report2)
print("########### c3 ###########\n",report3)
Reescalando para modelo final:
scaler = MaxAbsScaler()
scaler.fit(X_train) #dando o fit na base de treino
#transformando
X_train=scaler_val.transform(X_train)
X_test=scaler_val.transform(X_test)
Perceba que é complicado definir um ponto de corte, pois sempre ganhamos de um lado mas perdemos de outro. Vamos tomar $c=.25$ e depois apresentaremos uma maneira mais inteligente de definir o ponto de corte:
c=.25
logreg = LogisticRegression(solver='liblinear').fit(X_train, y_train)
#Prevendo probabilidade (segunda coluna para y=1)
y_proba=logreg.predict_proba(X_test)
y_pred = 1*(y_proba[:,1]>c)
report=classification_report(y_test, y_pred)
print(report)
Teoria da decisão: como escolher um bom ponto de corte?¶
Como vimos no exemplo passado, não é tão óbvio dizer qual o melhor ponto de corte. Se escolhemos corte=.5 tendemos a maximizar a acurácia, no entanto, nem sempre é isso que queremos fazer. Pode ser, por exemplo, que termos Falsos Negativos seja muito pior do que tem Falsos Positivos, logo deveríamos subir o ponto de corte. Concorda? Isso ocorre por exemplo no caso em que queremos diagnosticar uma doença grave, como é o caso do câncer - dizer que uma pessoa não tem câncer quando ela tem é muito mais grave do que dizer que ela tem quando não é verdade, pois estamos colocando a vida dessa pessoa em xeque. Se 1 é ter câncer e 0 não ter, parece que deveríamos diminuir o ponto de corte para diminuir nossa taxa de Falsos Negativos, mas até que ponto? Devemos lembrar que também existe um custo em termos Falsos Positivos, pois submetemos pessoas saudáveis a procedimentos custosos, logo não podemos diminuir o ponto de corte indefinidamente. Vamos descobrir agora como fazer essa escolha.
\
Suponha que tenhamos o ponto de corte $c\in [0,1]$ e a seguinte matriz de confusão normalizada, ou seja, com probabilidades e não números absolutos:
\
| Modelo/Realidade | Classe 0 | Classe 1 |
|---|---|---|
| Classe 0 | $p_1(c)$ | $p_2(c)$ |
| Classe 1 | $p_3(c)$ | $p_4(c)$ |
\
Considere que a normalização foi feita com relação à matriz inteira e não dentro de cada uma das linhas, como também é comum. Logo $p_1(c)+p_2(c)+p_3(c)+p_4(c)=1$ para qualquer $c \in [0,1]$, sendo que todas as entradas da matriz devem ser não negativas. Como o esperado, os valores da matriz de confusão dependem do ponto de corte $c$. Suponha que além da matriz de confusão podemos definir uma matriz de perdas da seguinte maneira:
\
| Modelo/Realidade | Classe 0 | Classe 1 |
|---|---|---|
| Classe 0 | $l_1$ | $l_2$ |
| Classe 1 | $l_3$ | $l_4$ |
\
Na matriz de perdas temos as perdas (e.g. monetárias) que incorremos ao fazer uma escolha e ao contrá-la com a realidade. Por exemplo, $l_2$ é a perda que incorremos ao termos um caso de Falso Positivo. Aqui assumimos que as perdas são constantes, ou seja, são fixadas a priori e não dependem de $c$, por exemplo. Agora que temos as probabilidades e perdas para cada um dos cenários, é natural definirmos a seguinte função de perda esperada:
\begin{align} L^e(c)=p_1(c) \cdot l_1+p_2(c) \cdot l_2+p_3(c) \cdot l_3+p_4(c) \cdot l_4=\sum_{j=1}^4 p_j(c) \cdot l_j \end{align}Também é natural assumir que $l_1=l_4=0$, pois essas perdas são de situações em que não erramos!!! Logo:
\begin{align} L^e(c)=p_2(c)\cdot l_2+p_3(c)\cdot l_3 \end{align}Na prática, definimos as perdas e plotamos $L^e(c)$ como função de $c$. Quando vemos o gráfico, decidimos qual será o ponto de corte ótimo!
Problema¶
Suponha que consultamos um especialista que nos disse que a cada Falso Positivo temos um custo de R\$30 e que a cada Falso Negativo temos um custo de R\$75. Se falamos que o paciente não faltará na consulta e ele falta, o custo é de R\$75 (o custo de uma consulta). Se falamos que o paciente faltará e ele não falta, o custo é de R\$30 ('bagunça na agenda').
Escolha um ponto de corte de acordo com a base de validação a fim de minimizar a função de perda esperada;
Utilize esse ponto de corte para treinar um modelo final e classificar exemplos do base de teste;
logreg = LogisticRegression(solver='liblinear').fit(X_train2, y_train2)
y_proba=logreg.predict_proba(X_val)
Escolhendo o ponto de corte
cs=np.linspace(start = 0, stop = 1, num = 100)
L=[]
l2=75
l3=30
for c in cs:
y_pred = 1*(logreg.predict_proba(X_val)[:,1]>c)
conf=confusion_matrix(y_pred,y_val)
p2=conf[0,1]/np.sum(conf)
p3=conf[1,0]/np.sum(conf)
L.append(p2*l2+p3*l3)
plt.plot(cs,L)
plt.show()
O custo médio dos erros que nosso modelo comete é de aproximadamente R\$18,50. Potencialmente diminuímos os custos de R\$21,50 para R\$18,50 (14\%) se utilizamos o modelo com outro ponto de corte (e.g. c=0.5)!!!
c=cs[np.argmin(L)]
c
Treinando modelo e classificando:
logreg = LogisticRegression(solver='liblinear').fit(X_train, y_train)
y_proba=logreg.predict_proba(X_test)
y_pred = 1*(logreg.predict_proba(X_test)[:,1]>c)
report=classification_report(y_test, y_pred)
print(report)