Estimando regressões logísticas no R (com razão de chance)

Uma das coisas que dá dor de cabeça a usuários do Stata que migram para o R é a estimação de modelos estatísticos. Como o Stata já oferece, de forma simples, diversos modelos e opções para alterar as suas especificações, a migração para o R pode ser frustrante neste aspecto: muitas coisas simples no Stata, como incluir erros-padrão robustos ou efeitos fixos, geralmente demandam mais linhas de código (e chamada a vários pacotes) no R.

Neste post, vou iniciar um pequeno guia para estimar regressões logísticas no R. O procedimento básico é bastante simples, como se poderá ver, mas alguns detalhes, como converter os coeficientes para razão de chance, envolvem algumas manhas adicionais.

Estimando modelos logísticos: explorando dados do Titanic

Antes de começar, precisamos de alguns dados para rodar um modelo logístico. Para isto, vou utilizar dados sobre os sobreviventes do Titanic, onde a variável de interesse (que usaremos como variável dependente nos modelos) é binária (ter sobrevivido, 1, ou não, 0). Esta base pode ser baixada aqui (meus tutoriais de como carregar dados no R, aqui e aqui), mas também podemos armazená-la num objeto diretamente no R via:

# Salva a base de dados sobre os sobreviventes do Titanic
titanic <- read.table("http://web.univ-ubs.fr/lmam/blanche/data-titanic.txt", sep = ",", header = T)

Feito isto, podemos explorar um pouco a base. Em particular, ela possui 891 observações e 12 variáveis, que são:

names(titanic)

##  [1] "PassengerId" "Survived"    "Pclass"      "Name"        "Sex"        
##  [6] "Age"         "SibSp"       "Parch"       "Ticket"      "Fare"       
## [11] "Cabin"       "Embarked"

O que mais nos interessa aqui é a variável binária Survived, que assume o valor de 1 para passageiros sobreviventes, e 0 para os não sobreviventes. Outras variáveis podem nos ajudar a explicá-la, como o sexo dos passageiros. Com alguns gráficos simples, podemos ter uma noção da relação entre estas duas variáveis:

# Primeiro, carregamos o pacote ggplot2 (para gerar os graficos)
library(ggplot2)

# Geramos um grafico de barras para visualizar o numero de sobreviventes e nao sobreviventes
ggplot(titanic, aes(x = as.factor(Survived))) + geom_bar() + theme_light() + 
  labs(title = "Titanic", x = "Sobreviventes", y = "Frequência")

# Por fim, cruzamos a variavel sexo com a variavel sobreviventes
ggplot(titanic, aes(x = as.factor(Survived), fill = Sex)) + geom_bar() + theme_light() + 
  labs(title = "Titanic", x = "Sobreviventes", y = "Frequência")

O que podemos concluir? Que a maior parte dos passageiros do Titanic, infelizmente, não sobreviveu ao acidente, e que, dos que sobreviveram, a maioria é do sexo feminino. Estimando uma regressão logística, podemos investigar esta relação mais detidamente.

Estimando modelos logísticos: a função glm do R

Para estimar essa regressão logística, uso a função glm, que já vem no R, especificamente no pacote stats. A função aceita vários argumentos, mas os principais que usaremos são três: a fórmula, o link (o que usaremos é o logístico, mas a função também aceita outros) e os dados usados (no caso, da base dos sobreviventes do Titanic, que já carregamos). Exemplo:

# Estima um modelo logistico com a funcao glm
modelo <- glm(Survived ~ Sex, family = "binomial", data = titanic)

Explicando o código acima, fizemos o seguinte. Primeiro, passamos a fórmula que queremos estimar, com a variável Survived (sobreviventes) como dependente (Y) sendo explicada (o ~ indica que tudo após isto é variável independente) pela variávei Sex (sexo dos passageiros). Segundo, explicitamos o link (family = "binomial"). Por fim, passamos os nossos dados (data = titanic).

Para visualizarmos o resultado deste modelo, apenas precisamos usar a função summary:

# Acessa os resultados do modelo logistico
summary(modelo)

## 
## Call:
## glm(formula = Survived ~ Sex, family = "binomial", data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6462  -0.6471  -0.6471   0.7725   1.8256  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.0566     0.1290   8.191 2.58e-16 ***
## Sexmale      -2.5137     0.1672 -15.036  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.7  on 890  degrees of freedom
## Residual deviance:  917.8  on 889  degrees of freedom
## AIC: 921.8
## 
## Number of Fisher Scoring iterations: 4

Como já havíamos visto, a probabilidade de passageiros homens terem sobrevivido é menor do que a de passageiras mulheres (aqui há um bom tutorial de como interpretar os resultados de modelos logísticos). Mas, da forma como é reportado, este coeficiente (Sexmale, que retorna a probabilidade, em log odds, de um homem ter sobrevivido em relação a uma mulher) é difícil de ser interpretado. Uma forma mais simples é por meio de razão de chance.

Tendo salvo o modelo, podemos calcular as razões de chance de forma fácil usando a função exp, que serve para exponenciar um número:

# Converte os coeficientes do modelo de log odds para razao de chance
exp(modelo$coefficients)

## (Intercept)     Sexmale 
##  2.87654321  0.08096732

Agora, podemos interpretar os resultado de forma mais direta: para cada mulher que se salvou no Titanic, 0.08 homens se salvaram (ou, de forma mais intuitiva, 100 mulheres para cada 8 homens; o artigo da Wikipedia sobre razão de chance é um bom guia para interpretar essas razões de chance).