O modelo logit é utilizado para medir probabilidade da existência de uma determinada características em um unidade \(i\), representado por \(Y_i=1\), utilizando informações \(X\):
\(P(Y_i = 1 | X) = \frac{e^{\beta_0 + \sum_i^{n} \beta_i X_i} }{1 + e^{\beta_0 + \sum_i^{n} \beta_i X_i}}\)
Dessa forma, esse modelo tem como objetivo determinar a probabilidade da existência ou não de terminada características a partir de variáveis explicativas.
A definição de quais informações fazem parte do \(X\) pode seguir um processo sistemático, útil, principalmente, há um grande número de variáveis candidatas a explicativas.
O método Stepwise permite a escolha do conjunto regressor de forma iterativa, admitindo a variação foward, partindo o modelo naive são incluída variáveis para o completo (todas as variáveis), ou backward, partindo do modelo completo são excluídas variáveis para um modelo mais parcimonioso.
Considerando o modelo Stepwise forward, segundo Izenman(2008), o processo de decisão entre inclusão/exclusão de variáveis se dá por meio da comparação da razão F com o valor 2 ou 4:
\(F = \frac{\frac{SS0 - SS1}{df0 - df1}}{\frac{SS1}{df1}}\)
sendo SS0 a soma de resíduos e df0 graus de liberdade do modelo de comparação, SS1 a soma de resíduos e df0 graus de liberdade do modelo mais completo.
Como exemplo, será utilizada a base de dados que identificam se um e-mail é spam com base conteúdo lexical do corpo do texto. Mais informações da base, tais como descrição e download podem ser acessadas no link.
Dessa forma, a variável \( Y\) identifica se a mensagem é spam ou não será modelada utilizando o modelo logit com seleção stepwise forward de \(X\), formado por 57 variáveis que determinam a intensidade de algumas palavras no corpo da mensagem.
Para isso, as bibliotecas utilizadas foram:
from scipy import stats
from scipy.optimize import fmin
import numpy as np
import math as math
import copy
A primeira função definida é estabelece o cálculo da log-verossimilhança do modelo logit:
def logitEst(betas):
mLike=0;
z = np.matrix(np.zeros((len(y),1),dtype = float))
for i in range(n):
for j in range(np.shape(xs)[1]):
z[i,0] = z[i,0] + xs[i,j] * betas[j]
for i in range(n):
try:
prob = (math.exp(z[i,0])) / (math.exp(z[i,0]) + 1)
except OverflowError:
prob = float(.9999999999)
logprob = math.log(stats.binom.pmf(y[i],1,prob))
mLike = mLike + logprob
return(-mLike)
Essa função é importante, pois será utilizada para o processo de otimização utilizando o pacote scipy.optimize. O chute inicial para o vetor de \(\beta\) do modelo é de 0.5 para todos.
def estimacaoRegBin(y,xs):
betas = np.transpose(np.matrix(np.ones(( 1 ,np.shape(xs)[1])) / 2 ))
betas = np.matrix(fmin(logitEst, betas ,xtol=1e-4))
return(betas)
Em seguida a base de dados baixa é importada por meio da seguinte programação:
martiz_resultado = np.matrix(np.empty((58,3),float))
## Importando a base
results = []
with open('C:\\Users\\Igor\\Documents\\Doutorado\\Métodos Computacionais\\Estrela4\\base.txt') as inputfile:
for line in inputfile:
results.append(line.strip().split(','))
base = np.array(results)
b = np.empty((4601,58),float)
for i in range(0,4601):
for j in range(0,58):
b[i,j] = float(base[i,j])
y = b[:,57]
n = len(y)
Em seguida, são definidas as variáveis iniciais para o processo stepwise, sendo vecvars o vetor que armazena as variáveis utilizadas pelo modelo.
vecvars = []
medida_ajuste = float(0)
RSS0 = 10**6
df0 = n
Por fim, o código a seguir realiza as chamadas das funções para o processo sequencial de inclusão de variáveis e armazena os resultados.
for j in range(-1,np.shape(b)[1]-1): # scanning variables
if j == -1:
xs = np.matrix(np.ones((len(y),1),dtype = float))
vecvars_temp = copy.deepcopy(vecvars)
else:
vecvars_temp = copy.deepcopy(vecvars)
vecvars_temp.append(j)
nvars = len(vecvars_temp)
xs = np.matrix(np.ones((len(y),nvars+1),dtype = float))
for i in vecvars_temp :
xs[:,i +1] = np.transpose(np.matrix(b[:,i]))
resultado = estimacaoRegBin(y,xs)
z = xs * np.transpose(resultado)
martiz_probz = np.matrix(np.empty((4601,3),float))
for i in range(0,n):
try:
prob = (math.exp(z[i,0])) / (math.exp(z[i,0]) + 1)
martiz_probz[i,0] = prob
except OverflowError:
prob = float(.9999999999)
martiz_probz[i,0] = prob
if (prob >= 0.5) :
z[i,0] = 1
else:
z[i,0] = 0
diag = 0
sum_erro = 0
for i in range(n):
if z[i,0] == float(y[i]):
diag = diag + 1
sum_erro = sum_erro + (z[i,0] - y[i])**2
df1 = n - np.shape(resultado)[1]
## Utilizando aproximação normal
Fvalue = ((RSS0 - sum_erro)/(df0 - df1))/(sum_erro/df1)
print("rss0= ",RSS0," e rss1=",sum_erro," com valor F=",Fvalue)
print("vf= ",vecvars)
print("vt= ",vecvars_temp)
print("df0=",df0," e df1=",df1)
## testando o quantidade de valores corretos
medida_ajuste_temp = float(diag / n )
if medida_ajuste_temp > medida_ajuste:
medida_ajuste = copy.deepcopy(medida_ajuste_temp)
## testando o valor F
if Fvalue > 2:
RSS0 = copy.deepcopy(sum_erro)
df0 = copy.deepcopy(df1)
vecvars= copy.deepcopy(vecvars_temp)
#armazenando SS
martiz_resultado[j+1,0] = sum_erro
martiz_resultado[j+1,1] = Fvalue
martiz_resultado[j+1,2] = medida_ajuste_temp
As variáveis que estão foram selecionadas para compor o modelo estão armazenadas no objeto vecvars.
A imagem a seguir representa o processo, sendo cada barra uma variável e seu tamanho a proporções da diagonal (acurácia) da matriz de confusão. As palavras representadas pelas barras cinzas não entram no modelo. Ao todo são 40 barras vermelhas, isto é, são´incluídas 40 variáveis das 57, sendo a acurácia final do modelo de 92.37%.
Referências
Izenman, A. (2008). Modern Multivariate Statistical Techniques. Assessment (Vol. 65, p. 733). Springer New York. https://doi.org/10.1007/978-0-387-78189-1
Lichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science.