Regressão Logística – Aula 1a

Regressão Logística – Aula 01 – parte a

A execução dessa primeira aula, com uma turma iniciante nos conceitos estatísticos e no uso do R, durou em torno de 6 horas. As duas partes (a e b) foram dividas em 3 horas e 3 horas (antes e depois do almoço)

Os Conceitos abordados foram:

  • Introdução a Modelos Lineares Generalizados
  • Modelo de Regressão Logística – LOGIT
    • Elementos do modelo LOGIT
    • Pressupostos do modelo LOGIT
    • Estimação do modelo LOGIT com o R
    • Interpretação do modelo LOGIT
    • Ajuste do modelo
    • Precisão da estimativa (estatística e gráfica)
    • Análise da contribuição dos previsores
    • Análise dos resíduos

Início da aula

A seguir, o ambiente de trabalho do R é preparado para a aula.
Uso normalmente os pacotes dplyr, tidyr, ggplot2 e readxl para facilitar a manipulação do banco de dados. Os outros pacotes broom,rms,pscl,car e pROC serão úteis no decorrer da aula (maioria na parte ‘b’).

if(!require(dplyr)){install.packages("dplyr")}
if(!require(tidyr)){install.packages("tidyr")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(broom)){install.packages("broom")}

Introdução aos modelos lineares generalizados

Quando pensamos em um Modelo Linear, pensamos em uma variável resposta (dependente) contínua $Y$ em relação a uma (ou mais) variável independente explicativa $X1$.

$latex\Y~=~beta0~+~beta1*X1~+~epsilon$

Olhando uma relação clássica e bem conhecida por exemplo é a relação entre hp (cavalos de potência, unidade: cavalos-vapor) e mpg (consumo de combustível, unidade: milhas por galão) do banco de dados mtcars.

ggplot(mtcars,aes(x=hp,y=mpg))+
  geom_point(color="red")+
  geom_smooth(method='lm',fill='yellow')+
  ylab("consumo (milhas por galão)")+xlab("cavalos de força (10 bhp)")+
  theme_bw()

E tal relação apresenta o seguinte resumo do seu modelo linear:

(mod <- lm(mpg~hp,data=mtcars))

Ou seja, a equação desse modelo linear seria:
$$
mpg = {30.1} -0.07*hp
$$

Com os valores de $\beta0$=30.098 e $\beta1$= -0.068

Mas um modelo linear não deve terminar (ou continuar) sem a avaliação dos resíduos dessa relação. Uma maneira muito rápida e eficiente é fazermos essa avaliação visualmente, observando os 4 gráficos a seguir, que levam em consideração os resíduos

par(mfrow=c(2,2))
plot(mod)

Podemos perceber, por exemplo, que se retirarmos o Maserati Bora, pois é um carro muito potente fabricado nos anos 1970, década na qual os carros apresentavam alto consumo antes que esse tipo de motor não fosse mais economicamente viável devido os altos valores que o petróleo veio alcançar.

A área em amarelo do gráfico refere-se a um intervalo de confiança de 0.95 (95%) para cada intervalo de valor da variável explicativa $X1$, no nosso caso: hp.
Essa variação faz com que a média de milhas percorridas por galão abastecido (mpg) caia conforme o valor do coeficiente calculado $\beta0$.

Conforme um tutorial publicado no site freakonometrics, podemos adaptar e construir uma função para visualizarmos como ocorre essa mudança da média conforme o nosso modelo linear entre duas variáveis:

plotar_curvas<-function(n=2,m=8,X,Y){
  df <- data.frame(X,Y)
  vX <- seq(min(X)-2,max(X)+2,length=n)
  vY <- seq(min(Y)-2,max(Y)+2,length=n)
  mat <- persp(vX,vY,matrix(0,n,n),zlim=c(0,.1),theta=-30,ticktype ="detailed", box = FALSE)
  reggig <- glm(Y~X,data=df,family=gaussian(link="identity"))
  x <- seq(min(X),max(X),length=501)

C <-trans3d(x,predict(reggig,newdata=data.frame(X=x),type="response"),rep(0,length(x)),mat)
  lines(C,lwd=2)
  sdgig <- sqrt(summary(reggig)$dispersion)
  x <- seq(min(X),max(X),length=501)
  y1=qnorm(.95,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig)
  C <- trans3d(x,y1,rep(0,length(x)),mat)
  lines(C,lty=2)
  y2 <- qnorm(.05,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig)
  C <- trans3d(x,y2,rep(0,length(x)),mat)
  lines(C,lty=2)
  C <- trans3d(c(x,rev(x)),c(y1,rev(y2)),rep(0,2*length(x)),mat)
  polygon(C,border=NA,col="yellow")
  C <- trans3d(X,Y,rep(0,length(X)),mat)
  points(C,pch=19,col="red")
  vX <- seq(min(X),max(X),length=m)
  mgig <- predict(reggig,newdata=data.frame(X=vX))
  sdgig <- sqrt(summary(reggig)$dispersion)
  for(j in m:1){
    stp <- 251
    x <- rep(vX[j],stp)
    y <- seq(min(min(Y)-15,qnorm(.05,predict(reggig,newdata=data.frame(X=vX[j]),type="response"), sdgig)),max(Y)+15,length=stp)
    z0 <- rep(0,stp)
    z <- dnorm(y, mgig[j], sdgig)
    C <- trans3d(c(x,x),c(y,rev(y)),c(z,z0),mat)
    polygon(C,border=NA,col="light blue",density=40)
    C <- trans3d(x,y,z0,mat)
    lines(C,lty=2)
    C <- trans3d(x,y,z,mat)
    lines(C,col="blue")}
}

No qual observaríamos a seguinte figura:

plotar_curvas(X = mtcars$hp, Y = mtcars$mpg)

Dessa maneira há a variação linear da média
$$
E(Y|X=x) = \beta0 +\beta1*x
$$

Com a variância constante, seguindo uma distribuição normal
$$
Var(Y|X=x)= \sigma^2
$$

Mas e se quisermos saber se a relação entre 'mpg' e uma outra variável, como se o carro possui câmbio automático ou manual (am):

ggplot(mtcars,aes(y=am,x=mpg))+
  geom_point(color="red")+
  geom_smooth(method='lm',fill='yellow')+
  theme_classic()+
  ylab("Manual - 0  ou Automático - 1")+
  xlab("Cavalos de potência (10 bhp)")
(mod_linear <-lm(am~mpg, data=mtcars) )

Mas observando os resíduos:

par(mfrow=c(2,2))
plot(mod_linear)

E observando como a variável resposta (am) se comportam em relação a uma variável indenpendente (mpg), podemos concluir que o modelo linear como foi aplicado não é indicado, pois corrompe os pressupostos para a construção de um modelo linear, pois:

  • A variável dpendente não segue uma distribuição normal
    • Ela é discreta (1 ou 0);
    • É bimodal (moda diferente da média e da mediana)
    • Ela segue uma distribuição binomial ($Y~Binomial(np, np(1-p)$)
  • Resíduos não mostram aleatoriedade
  • Heterocedasticidade dos dados
plotar_curvas(X = mtcars$mpg, Y = mtcars$am)

Nesse caso, o melhor é usarmos um Modelo Linear Generalizado (Generalized Linear Model – GLM), pois podem ser usados para inferirmos uma variável dependente (resposta $Y$ ) que não sigam uma distribuição normal

GLM consiste em três componentes:
* Componente estrutural (structural component): $\beta0 +\beta1X1… $
* Função de ligação (
link function): $g(\mu)$
* Distribuição da variável resposta (
response distribution*): Neste caso seria $Y \~ Binomial$

$$
g(\mu) = \beta0 + \beta0Xg(\mu) = \beta0 + \beta1X1…
$$

Outro detalhe importante é que ao invés de usar o método de mínimos quadrados (Ordinary Least Squares – OLS), GLM utiliza estimativa de verossimilhança (Maximum Likelihood Estimation – MLE), o que modifica como calculamos o ajuste de um modelo.

ggplot(mtcars,aes(x=mpg,y=am))+
  geom_point()+
  geom_smooth(method='lm',aes(color="Modelo Linear"),se=F)+
  geom_smooth(method = "glm", method.args = list(family = "binomial"), aes(color="GLM Binomial"),se=F )+
  scale_y_discrete(limits=c(0,1))+
  xlab("Consumo (milhas por galão)")+ylab("Manual - 0 / Automático - 1")+
  theme_bw()+  geom_jitter(width=1,height=0)+
  scale_color_manual(name="", values=c("blue","red") )

Portanto, verificar o ajuste de um modelo linear generalizado (GLM) não é mais possível fazer na mesma maneira que fazemos com os modelos lineares (LMs):
$R^2$ (para verificarmos ajuste) e plot(mod) (para verificarmos os resíduos).

Regressão Logística

Esse tipo de regressão é muito interessante, pois é muito útil para inferirmos valores em que as variáveis independentes (explicativas) têm e que levam a um ponto de ruptura para a nossa variável dependente (resposta), pois a nossa variável resposta só apresenta dois valores (1 ou 0). Podemos entender esse tipo de resposta para alguns exemplos, como:
* Compra ou não compra de um determinado item
* Morte ou sobreviência
* Falência ou sucesso de uma empresa
* Risco no Crédito ou não
* Infecção ou não
* Ruptura do material ou não
* Gol/cesta ou pra fora

Como $p$ (probabilidade de sucesso) tem valores entre 0 e 1 ($0 \ge p \ge 1$), precisamos adequar a nossa equação para a nossa variável resposta.

Portanto, para que $p \ge 0$:
$$
p = exp(\beta0+\beta1X1) = e^{\beta0+\beta1X1}
$$

E para que $p \le 1$:
$$
p = \frac{exp(\beta0+\beta1X1)}{1+exp(\beta0+\beta1X1)} = \frac{e^{\beta0+\beta1X1}}{1+e^{\beta0+\beta1X1}}
$$

sendo que $p$ é igual a $\hat{y}$. E quando nos referimos a $p$ (ou $\hat{y}$), estamos nos referindo a uma escala de probabilidade, que é diferente de:

  • Escala de Chance (Odd`s Scale)
    $$
    probit(\hat{y})= \frac{\hat{y}}{1-\hat{y}}=exp(\beta0+\beta1X1) = e^{\beta0+\beta1X1}
    $$
  • Escala de Log da Chance (Log-odd`s Scale)
    $$
    logit(\hat{y})=log(\frac{\hat{y}}{1-\hat{y}}) = \beta0+\beta1*X1
    $$

Cada uma dessas escalas (de probabilidade, PROBIT e LOGIT) tem suas vantagens e desvantagens, e nós faremos uma análise comparativa posteriormente. Mas para fazermos a escolha entre as escalas será necessário uma combinação de:
+ Conhecimento da distribuição da variável resposta
+ Considerações teóricas
+ Ajuste empírico aos dados

Mas um bom material disponível pode ajudar a entender a diferença entre probabilidade de um evento ($p$ ou $\hat{y}$) e chance ($odd$) de um evento.

Diferenças gráficas entre probabilidades, PROBIT e LOGIT

mod %
  mutate(y_hat = .fitted) %>%  # Valores de probabilidade para cada observacao
  mutate(odds = y_hat/(1-y_hat)) %>% # Valores de chance para cada observacao (PROBIT)
  mutate(log_odds =log(odds)) %>%  # Valores de log(chance) para cada observacao (LOGIT)
  select(am,mpg,y_hat,odds,log_odds)

print(ajuste)

Gráfico da relação entre $p$ e mpg. Perceba que o gráfico é igual a curva da relação GLM Binomial (e os valores variam de 0 a 1).

ggplot(ajuste, aes(x = mpg, y = y_hat)) +
  geom_point() + geom_line() +
  scale_y_continuous("Probabilidade de ser automático")

Gráfico da relação entre $odd$ e mpg, Perceba que a chance de ser automático aumenta conforme aumentamos o valor de mpg, em uma relação exponencial

ggplot(ajuste, aes(x = mpg, y = odds)) +
  geom_point() + geom_line() +
  scale_y_continuous("Chance (PROBIT) de ser automático")

Gráfico da relação entre $log(odd)$ e mpg. Perceba que a relação fica reta, mas com valores que não fazerm muito sentido para $log(odd)$.

ggplot(ajuste, aes(x = mpg, y = log_odds)) +
  geom_point() + geom_line() +
  scale_y_continuous("log(chance) (LOGIT) de ser automático")

Modelo LOGIT com o R

Seguiremos os seguintes passos:

1. Importar os dados;

2. Fazer uma análise exploratória dos dados

  • Verificar valores NA ou NaN
  • Verificar outliers
  • Verificar se as variáveis preditoras que utilizaremos não são altamente autocorrelacionadas

3. Construir dois grupos de dados (treino e teste)

4. Construir o modelo

5. Interpretar os coeficientes do modelo com os dados treino

6. Analisar a contribuição dos previsores

  • Sumário do modelo
  • Estatística Wald
  • Intervalo de Confiança dos coeficientes

7. Analisar o modelo treino utilizando seus resíduos

  • Análise Gráfica de Modelo Marginal (Marginal Model Plots)
  • Análise gráfica de Componentes + Resíduos e Resíduos Parciais (Component + Residual (Partial Residual) Plots)

8. Verificar o ajuste do modelo treino

  • $Pseudo R^{2}$
  • $p-value$ de um teste $\chi^2$ com um modelo nulo
  • Curva ROC (Receiver Operating Characteristic) e AUC (Area Under the Curve)
  • Tabela confusão
  • AIC (Akaike Information Criteria)

9. Analisar a precisão da estimativa com os próprios dados treino e com os dados teste

  • Tabela confusão

Todas esses passos serão realizados e postados na segunda parte da aula (parte b), onde cada um desses passos serão explicados detalhadamente.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s