domingo, 28 de outubro de 2012

Regressão Logística em R

Olá pessoal, tudo certo? Já faz um tempo que não posto nada, mas a partir deste mês pretendo escrever com maior frequência. No post de hoje vamos falar um pouco sobre regressão logística, uma técnica bastante popular e sútil do ponto de vista teórico.

Em Janeiro de 1986 o mundo inteiro acompanhou o desastre no lançamento do ônibus espacial Challenger, que causou a morte dos sete tripulantes. Posteriormente, uma falha dos anéis de vedação (no tanque externo de combustível sólido) foi apontada como a causa do desastre. Como exemplo vamos usar os dados de Temperatura x Falha do anél O-ring, que estavam disponíveis antes da tragédia.

Introdução

A regressão logística pode ser usada quando desejamos modelar uma variável categórica em função de uma ou mais variáveis preditoras (numéricas, categóricas ou ambas). Na maioria dos casos, temos uma variável resposta binária e estamos interessado em quantificar as 'chances' de obter uma resposta positiva. Repare que não podemos simplesmente usar uma regressão linear pelos seguintes motivos:

  • A linha da regressão linear define o valor médio da variável resposta, podendo assumir valores entre -Inf e +Inf. No caso de uma variável binomial, o valor esperado está restrito entre [0, 1].
  • Em segundo lugar, a regressão linear assume que a variância do termo de erro é constante e independe dos valores assumidos pelos preditores. Isso não ocorre quando a variável resposta é binomial.
  • Por fim, os resíduos não podem ser normalmente distribuídos, já que a variável resposta assume apenas dois valores possíveis.
Uma possível solução é pensar no conceito de odds (chances). Seja Y a variável resposta, que pode assumir os valores 0 ou 1, e Xn o vetor representando as n variáveis regressoras. A idéia é transformar Y em uma variável numérica, e em seguida aplicar uma regressão linear. Para isso, seja p a probabilidade de Y assumir o valor 1, podemos então definir:

odds (p) = p/(1 - p)

Repare que odds (Y) está limitada entre [0, +Inf), o que ainda nos impede de aplicar uma regressão linear. Podemos então definir a função logit:

logit (p) = log(odds(p)) = log (p/(1 - p))

Agora estamos no intervalo (-Inf, +Inf) e podemos assumir um modelo linear nos preditores:

log (p/(1 - p)) = βo + β1 x1 + β2 x2 + ... + βn xn

Calculando a função inversa de logit (p), obtemos a função logística:

p (t) = 1/(1 + e-t)

ou

p = 1/(1 + e-(βo + β1 x1 + β2 x2 + ... + βn xn))

A equação acima nos diz a probabilidade de se obter Y = 1 dado Xn. Observe que o modelo é linear entre logit (p) e Xn, porém, a relação entre p e Xn é não-linear. Por exemplo, ao se aumentar o preditor x1 em uma unidade (mantendo os demais constantes), odds (p) sofrerá uma variação dada por eβ1. Agora vamos ver como implementar os conceitos acima em R.

R e GLM's

Em R não existe uma função específica para ajustar um modelo de regressão logística, e o motivo é simples: regressão logística é apenas um caso de modelo linear generalizado, ou GLM em inglês. Nesse tipo de modelo especificamos apenas a distribuição do erro e a função de link. O código abaixo vai deixar tudo claro, então baixe o dataset e salve no seu diretório de trabalho.

    dados.nasa = read.csv('nasa_disaster.csv')
    modelo = glm(failure ~ temperature, data = dados.nasa, family = binomial(link = 'logit'))
    summary(modelo)
   

Como nossa variável de interesse é binária (1 = houve falha, 0 = não houve falha), especificamos uma distribuição binomial e usamos a função logística como link (que é o default nesse caso, veja help(glm)). Pelo output vemos que a temperatura possui efeito significativo na possibilidade de falhas (p value ~ 0.04). Observando os coeficientes do modelo, podemos escrever:

log (p/(1 - p)) = 10.87535 - 0.17132 * T

ou

p = 1/(1 + e-(10.87535 - 0.17132 * T))

Onde p é a probabilidade de falha e T a temperatura. Graficamente:

    temp = dados.nasa$temperature
    p = modelo$fitted.values
    xlb = expression(paste('Temperatura [', degree,'F]'))
    titulo = 'NASA - Desastre Challenger\nProb falha x Temperatura'
    plot(temp, p,  col = 'navy', pch = 4, cex = 1.5, ylab = 'p', xlab = xlb, main = titulo)
    grid(10, 10, col = '#CCCCCC')
    

O gráfico deixa claro a influência da temperatura sobre a possibilidade de falha. Como comentado na introdução, ao variarmos a temperatura em 1 grau, as chances de falha variam em e-0.17132 = 0.8425515. Ou seja, ao aumentarmos a temperatura em 1 grau, as chances de falha são reduzidas em 0.8425515. Esse valor é chamado de odds ratio, e é uma constante característica do modelo.

Podemos facilmente prever as probabilidades de falha para as temperaturas que não estão no dataset original:

    new.temp = c(54, 55, 58, 59, 60, 61, 62, 64, 65, 71, 74, 77)
    new.data = data.frame(temperature = new.temp)
    probs = predict(modelo, new.data, type = 'response')
    points(new.temp, probs, pch = 3, col  ='red', cex = 1.5)
    legenda = c('modelo', 'previsto')
    pontos.tipo = c(4, 3)
    cores = c('navy', 'red')
    legend('topright', legenda, col = cores, pch = pontos.tipo)
    

A temperatura dos anéis de vedação no dia do lançamento era de 31°F, valor bem distante dos valores usados no nosso modelo (dados de testes anteriores da NASA). Podemos dizer que uma previsão a essa temperatura provavelmente não seria precisa, mas observando o gráfico, fica claro que um lançamento a 31°F provavelmente vai resultar em uma falha nos anéis de vedação. É incrível como uma análise simples poderia ter ao menos alertado os Engenheiros.

Voltando ao código, a função predict é bem simples de ser usada: fornecemos o modelo, um data frame com as novas temperaturas e especificamos o tipo como 'response', o que permite obter as probabilidades diretamente, ao invés dos log-odds. Lembrando que essa é uma função genérica, portanto seu uso varia de acordo com o tipo de modelo (veja help(predict.glm)).

O estudo de GLM's é vasto, essa foi apenas uma introdução prática! Para mais exemplos de uma olhada em help(glm), e fique ligado nos posts. Era isso, até a próxima!

3 comentários:

  1. Ótimo, mais um belo blog tratando de R! Bem legal mesmo vou fazer a ligação lá no notasr.wordpress.com
    Sigamos na blogagem R!!

    Abraço!

    ResponderExcluir
  2. Opa, valeu! Tambem vou dar um folllow no obiomassa
    abs

    ResponderExcluir
  3. Olá pessoal...possuo uma variável dependente categórica, e duas independentes(uma quantitativa e outra categórica) o teste que uso para analisar é regressão logística? Preciso transformar a variável independente categórica em quantitativa ?

    ResponderExcluir