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).
)