Aplicação de uma regressão logística em Python

     Em análise discriminante, o caráter não-métrico de uma variável dependente dicotômica é acomodado fazendo-se previsões de pertinência a grupo baseadas em escores Z discriminantes. Isso requer o cálculo de escores de corte e a designação de observações a grupos. A regressão logística se aproxima desta tarefa de uma maneira, mas semelhante à encontrada na regressão múltipla (veja nosso post sobre regressão múltipla no Python). Ela difere da regressão múltipla, contudo, no sentido de que ela prevê diretamente a probabilidade de um evento ocorrer. O valor previsto deve ser limitado, de modo a recair no intervalo entre zero e um.

     Enquanto a regressão múltipla emprega o método dos mínimos quadrados, o qual minimiza a soma de diferenças quadradas entre os valores reais e os previstos, a natureza não-linear da transformação logística demanda que um outro procedimento, o de máxima verossimilhança, seja usado de forma iterativa para encontrar as estimativas mais prováveis para os coeficientes.

     Indo para a prática, vamos implementar uma regressão logística em Python, utilizando uma base fictícia de clientes de uma financeira. Nosso objetivo é determinar a probabilidade de default de um indivíduo. Este exemplo foi tirado de minha aula de pós graduação, resolvido em SPSS, mas adaptado aqui para o Python.

     Começamos importando as bibliotecas necessárias, carregando os dados e, em seguida, olhando o dataset:

import matplotlib.pyplot as plt 
import pandas as pd  
import statsmodels.api as sm  
import pylab as pl  
import numpy as np  
from patsy import dmatrix  

#read the data in  
df = pd.read_csv('C:/Python27/bankloan.csv')  

#take a look at the dataset  
print 'Dataset sample'  
print df.head()  
print

Como saída desta primeira etapa temos:

   
     Aqui visualizamos as 5 primeiras linhas (enumeradas de 0 a 4) das 8 variáveis que temos disponíveis (idade, educação, t_emprego, t_endereco, renda, divida, divida_cc, outras_div) e nossa variável dependente default.

     Prosseguimos com uma análise descritiva da base, a visualização do desvio de cada variável e um exemplo de relação entre a idade da pessoa e a quantidade de defaults.

#summarize the data 
print '------------------------------------' 
print 'Data Summary' 
print df.describe() 
print 


#take a look at the standard deviation of each column 

print '------------------------------------' 

print 'Std Deviation of each column' 

print df.std() 

print

     Vamos agora plotar o histograma de cada uma de nossas variáveis:

#plot all of the columns
df.hist()
pl.show()

     Agora um ponto de atenção: como nossa variável idade é do tipo categórica (1 indica o menor nível e 5 o maior), utilizaremos variáveis dummies para cara classificação. Assim, teremos 5 novos atributos, um para cada nível.

#dummify educacao
dummy_educacao = pd.get_dummies(df['educacao'],prefix='educacao')

     Agora selecionamos os atributos que serão usados no processo de modelagem, retirando a variável educação e acrescentando as suas dummies, excluindo-se a primeira, para evitar os efeitos de correlação:

cols_to_keep = ['default','idade','t_emprego','t_endereco','renda','divida','divida_cc','outras_div']
data = df[cols_to_keep].join(dummy_educacao.ix[:,'educacao_2':])

#manually add the intercept
data['intercept'] = 1.0
Nossa próxima etapa é rodar a regressão:

#performing the regression
train_cols = data.columns[1:]
mLogit = sm.Logit(data['default'],data[train_cols])

#fit the model
mLogit_res = mLogit.fit()
print
print '------------------------------------'
print 'Regression Summary'
print '------------------------------------'
print mLogit_res.summary()

     Avaliando nossa saída, vemos que as variáveis realmente significativas para elaborarmos uma boa previsão de risco de default são: t_emprego, t_endereço e divida_cc, pois apresentam p valor inferior a 0.05 (nosso nível de significância aceitável), ou seja, rejeitamos a hipótese nula.

     Desta forma, devemos rodar novamente a regressão, selecionando apenas estas variáveis. Assim, reescrevemos:

cols_to_keep = ['default','t_emprego','t_endereco','divida_cc',]
data = df[cols_to_keep]

     E temos nossa nova saída:

     Note que todas as variáveis estão dentro do nosso nível de significância aceitável. Agora, vamos plotar o diagnóstico dos resíduos, a distância de cook e a curva ROC.

Q = mLogit.predict(mLogit_res.params,linear=True)

## Calculate sensitivity and specificity
S = np.linspace(min(Q), max(Q), 100)
Sens = np.zeros(len(S), dtype=np.float64)
Spec = np.zeros(len(S), dtype=np.float64)
i1 = np.flatnonzero(data['default'] == 1)
i0 = np.flatnonzero(data['default'] == 0)
for i,s in enumerate(S):
    Sens[i] = np.mean(Q[i1] > s)
    Spec[i] = np.mean(Q[i0] <= s)

#calculate the AUC using the trapezoidal rule
auc = 0.
for i in range(1,len(Sens)):
    auc += (Spec[i]-Spec[i-1]) * (Sens[i] + Sens[i-1])/2

#plot the ROC curve
plt.clf()
plt.grid(True)
plt.plot([0,1], [0,1], '-', color='black', lw=2)
plt.plot(1-Spec, Sens, '-', color='blue', lw=2)
plt.xlabel("1 - Specificity", size=12)
plt.ylabel("Sensitivity", size=12)
plt.title("AUC=%.2f" % auc)
plt.show()

#parameters estimated
print
print '------------------------------------'
print 'Parameters'
print mLogit_res.params

#Analise de residuos
plt.figure()
plt.scatter(Q,mLogit_res.resid_pearson)
plt.plot([0.0,1.0],[0.0,0.0],'k-')
plt.title('Residual Dependence Plot')
plt.ylabel('Pearson Residuals')
plt.xlabel('Fitted Values')
plt.show()


plt.figure()
plt.scatter(Q,df['default'])
plt.title('Cooks Distance')
plt.ylabel('Fitted values')
plt.xlabel('Measured values')
plt.show()

     Geralmente, a sensibilidade e a especificidade são características difíceis de conciliar, isto é, é complicado aumentar a sensibilidade e a especificidade de um teste ao mesmo tempo. As curvas ROC (receiver operator characteristic curve) são uma forma de representar a relação, normalmente antagónica, entre a sensibilidade e a especificidade de um teste diagnóstico quantitativo, ao longo de um contínuo de valores de "cutoff point". A interpretação é: quanto mais a curva se aproxima do canto superior esquerdo do diagrama, maior é a exatidão da regressão. Nesse caso, de uma escala de 0 a 1 possível, conseguimos 0,83.

     Finalizamos nosso estudo montando a equação da regressão e calculando a probabilidade de um evento ocorrer:

#parameters estimated
print
print '------------------------------------'
print 'Parameters'
print mLogit_res.params

Com nossos parâmetros estimados, montamos a equação:

Assim, como exemplo, calculando a probabilidade de default de um indivíduo com t_emprego = 3, t_endereço = 5 e divida_cc = 0,70:

Finalizamos aqui nosso post. Um abraço e até o próximo!

Bibliografia de Referência
HAIR, J. F., ANDERSON, R. E. TATHAM, R. L., BLACK, W. C. Análise Multivariada de Dados. 5. ed. Porto Alegre: Bookman, 2006.