Capítulo 13 - GLM
GLM é a sigla de Modelos Lineares Generalizados. A regressão logística pode ser considerada como um dos tipos de modelos lineares generalizados. Vamos nos concentrar em entender a regressão logística, como ela se conecta com GLM e, rapidamente, ver outros modelos de GLM que existem, como Poisson e probit. Até o momento, variáveis resposta binárias foram modeladas com regressão linear.
13.1 PLM
Até o momento, variáveis respostas binárias foram modeladas com regressão linear.
\[ y_i \sim N(\alpha + \beta x_i, \sigma^2) \]
Nessa parametrização do modelo de regressão, vemos que a resposta é modelada como uma variável Gaussiana e, portanto, é uma variável contínua (em vez de binária). Além disso, não está limitada, possuindo suporte entre −∞ e +∞. Esses modelos são chamados também de Modelos de Probabilidade Linear (MPL, ou LPM na sigla em inglês). Explicaremos mais à frente porque são chamados desse modo, após apresentarmos a regressão logística e probit. E então poderemos comparar modelos logísticos ou probit e MPL.
13.2 Logística
Há várias formas de apresentar e/ou justificar a regressão logística. Qual delas você irá usar para pensar esse tipo de modelo depende do seu problema de pesquisa e das suas preferências sobre o que funciona melhor para você.
13.2.1 Logística como melhoria em relação ao MPL
Como vimos, o problema da regressão linear para dados binários é que considera que a variável resposta possui distribuição Gaussiana. Faz muito mais sentido modelar dados binários como seguindo uma distribuição de Bernoulli, com parâmetro \(p_i\). Poderíamos, portanto, tentar reescrever um modelo para \(y\) binário da seguinte forma:
\[ y_i \sim Ber(p_i) \] Em que \(p_i\) é a probabilidade de sucesso para a unidade \(i\). Dessa forma, contudo, não incluí nenhum preditor para estimar o \(p_i\), e gostaríamos de fazê-lo. Posso então escrever algo como: \[ p_i(x_i) = \alpha + \beta x_i \]
O problema dessa formulação é que probabilidades devem estar entre 0 e 1, e nada garante que \(\alpha + \beta x_i\) seja um número entre 0 e 1. Então, uma saída é tentar achar uma função \(f\) que transforme \(\alpha + \beta x_i\) em números entre 0 e 1, ou seja, \(0 \le f(\alpha + \beta x_i) \le 1\). E uma função que tem essa propriedade é a logística padrão (também chamada de sigmoide), dada por:
\[ f(x) = \frac{1}{1 + \exp(-x)} \]
Então, nossa relação entre a probabilidade e os preditores fica:
\[ p_i = \frac{1}{1 + \exp(-(\alpha + \beta x_i))} \]
Conectando com nossa variável resposta, temos:
\[ y_i \sim Ber(\frac{1}{1 + \exp(-(\alpha + \beta x_i))}) \] A probabilidade de sucesso, condicional ao preditor, fica então: \[\begin{align} Pr(y_i=1|x) = \frac{1}{1 + \exp(-(\alpha + \beta x_i))} \tag{13.1} \end{align}\]
Assim, a modelagem de variável binária por meio da regressão linear pode ser pensada como um modelo em que a probabilidade \(p_i\) é modelada sem a função logística, o que significa que podemos acabar com previsões de probabilidades negativas ou maiores que 1, o que não faz sentido. Essa é uma razão para usar modelos logísticos em vez de lineares.
13.2.2 Logito
Às vezes a logística é designada como a função inversa da logito, em que a logito é dada por: \[ logito(p) = log(\frac{p}{1-p}) \] Para entender isso, vamos lembrar que, se \(f(x) = x + 2\) é uma função, sua inversa \(f^{-1}(x)\), se existir, pode ser descoberta pelo algoritmo em que chamamos \(f(x)\) de \(y\), trocamos \(y\) por \(x\) e resolvemos para \(y\): \[\begin{align} y = x + 2 \\ x = y + 2 \\ y = x - 2 \\ f^{-1}(x) = x - 2 \end{align}\]
Para uma função \(f(x) = log(x+2)\), a inversa é:
\[\begin{align} y = log(x + 2) \\ x = log(y + 2) \\ \exp(x) = \exp(log(y + 2)) \\ \exp(x) = y + 2 \\ y = \exp(x) + 2 \\ f^{-1}(x) = \exp(x) + 2 \end{align}\]
Então, se \(f(X) = log(x)\), \(f^{-1}(x) = \exp(x)\). Disso se segue que:
\[\begin{align} f(x) = log(\frac{x}{1-x}) \\ y = log(\frac{x}{1-x}) \\ x = log(\frac{y}{1-y}) \\ \exp(x) = \exp(log(\frac{y}{1-y})) \\ \exp(x) = \frac{y}{1-y} \\ \exp(x)(1-y) = y \\ \exp(x)- y\exp(x) = y \\ \exp(x) = y + y\exp(x) \\ \exp(x) = y(1 + \exp(x)) \\ \frac{\exp(x)}{(1 + \exp(x))} = y \\ f^{-1}(x) = \frac{\exp(x)}{1 + \exp(x)} \\ f^{-1}(x) = \frac{\frac{\exp(x)}{\exp(x)}}{\frac{1 + \exp(x)}{\exp(x)}} \\ f^{-1}(x) = \frac{1}{\exp(-x) + 1}\\ f^{-1}(x) = \frac{1}{1 + \exp(-x)}\\ \end{align}\]
No R, podemos acessar as duas funções por meio de:
E podemos modelar os dados de uma logística com nosso exemplo do Latinobarômetro usando apenas um preditor, ideologia, para simplificar.
library(here)
library(data.table)
library(tidyverse)
library(sjlabelled) # pra remover labelled variables
library(haven)
library(janitor)
library(lubridate)
library(knitr)
library(broom)
## dados
# https://www.latinobarometro.org/latContents.jsp
lat_bar23 <- sjlabelled::read_spss(here("Dados", "Latinobarometro_2023_Eng_Spss_v1_0.sav"),
drop.labels = TRUE) %>%
mutate(S17 = as.Date(as.character(S17), "%Y%m%d")) %>%
clean_names()## Invalid date string (length=9): 09 032 23# get_label(lat_bar23)
lat_bar23 <- lat_bar23 %>%
mutate(data_base = as.Date(paste(diareal, mesreal, "2023", sep="-"), "%d-%m-%Y"),
idade = year(as.period(interval(s17,data_base))),
econ_12_meses = ifelse(p6stgbs %in% c(1,2), "better",
ifelse(p6stgbs == 8, NA, "other")),
econ_12_meses = relevel(as.factor(econ_12_meses), ref = "other"),
aprovacao_presidente = ifelse(p15stgbs == 0, NA, p15stgbs),
ideologia = ifelse(p16st %in% c(97, 98, 99), NA, p16st),
votaria_governo = ifelse(perpart == 4, NA,
ifelse(perpart == 1, 1, 0)),
genero = factor(sexo, labels = c("homem", "mulher")),
evangelico = ifelse(s1 %in% c(0,98), NA,
ifelse(s1 %in% c(2,3,4,5), 1, 0))) # não considera adventista, testemunha Jeová, Mórmon## Warning: There was 1 warning in `.fun()`.
## ℹ In argument: `idade = year(as.period(interval(s17, data_base)))`.
## Caused by warning in `class(xx) <- cl`:
## ! Definindo class(x) para múltiplas strings ("POSIXct", "POSIXt", ...); resultado não será mais um objeto S4br_latbar_23 <- lat_bar23 %>%
mutate(idenpa = remove_all_labels(idenpa)) %>% # haven_labelled problems
filter(idenpa == 76) %>% ## seelciona brasil
filter(!is.na(votaria_governo) & !is.na(evangelico) & !is.na(ideologia) & !is.na(econ_12_meses))
reg_logistica <- glm(votaria_governo ~ ideologia, data=br_latbar_23,
family=binomial(link= "logit"))
reg_logistica %>%
tidy() %>%
kable()| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 3.1399152 | 0.3377743 | 9.295898 | 0 |
| ideologia | -0.2850991 | 0.0385935 | -7.387224 | 0 |
# plotando # usando base R
# código adaptado de Regression and Other Stories (Gelman et. al.)
library(arm)
n <- nrow(br_latbar_23)
ideologia_jitt <- br_latbar_23$ideologia + runif(n, -.2, .2)
vote_jitt <- br_latbar_23$votaria_governo + ifelse(br_latbar_23$votaria_governo==0, runif(n, .005, .05), runif(n, -.05, -.005))
curve(invlogit(reg_logistica$coef[1] + reg_logistica$coef[2]*x), from = -5,to=24, ylim=c(0,1), xlim=c(-5, 24), xaxt="n", xaxs="i",
ylab="Pr (voto governo)", xlab="Ideologia", lwd=.5, yaxs="i")
curve(invlogit(reg_logistica$coef[1] + reg_logistica$coef[2]*x), 1, 11, lwd=3, add=TRUE)
axis(1, 1:11)
mtext("(left)", side=1, line=1.7, at=1, adj=.5)
mtext("(right)", 1, 1.7, at=11, adj=.5)
points(ideologia_jitt, vote_jitt, pch=20, cex=.1)
O gráfico apresenta o modelo ajustado da logística, juntamente com os dados (jittered). Cada ponto na curva representa \(Pr(Y=1|X=x)\). A linha mais escura é o modelo ajustado para os pontos da amostra. Podemos observar que há mais dados de pessoas que votariam com o governo entre os indivíduos de esquerda, com mais observações votando contra o governo apenas entre as pessoas mais centristas, o que faz sentido. Também observamos muitas pessoas de direita votando com o governo, indicando que entre as pessoas de direita, o dado não diferencia tanto quem vota a favor ou contra o governo, sugerindo que precisamos considerar outras variáveis.
Como a logística é curva, o efeito preditivo de \(x\) sobre a probabilidade \(y=1\) não é constante (ao contrário de regressões lineares). No caso do nosso gráfico, passar de \(1\) para \(2\) na ideologia tem um efeito significativamente menor (de 95% para 93%) do que passar de \(6\) para \(7\) (de 81% para 76%), por exemplo. Para calcular o efeito preditivo nesses casos, fazemos:
newdata <- data.frame(ideologia = 1:11)
previsao <- predict(reg_logistica, newdata =newdata, type = "response")
print(round(previsao, 2))## 1 2 3 4 5 6 7 8 9 10 11
## 0.95 0.93 0.91 0.88 0.85 0.81 0.76 0.70 0.64 0.57 0.50A interpretação do coeficiente é um pouco difícil na logística, mas um caminho é pensar que o efeito é máximo no centro da curva, em que \(\alpha + \beta x = 0\). Nesse caso, a inclinação da curva (a derivada) é dada por \(\beta/4\). Ou seja, podemos dividir o coeficiente estimado por quatro para ter uma ideia do máximo impacto preditivo. Em nosso caso, com um coeficiente aproximado de \(-0,29\), temos que no máximo a mudança em uma unidade diminui a probabilidade em no máximo 7,25%. E nos dados, esse efeito máximo é alcançado quando \(3.14 + -0.29*x = 0\), ou seja, aproximadamente \(x = 11\).
13.2.3 Odds Ratio - Razão de chances
Outra forma de interpretar a regressão logística é utilizando a parametrização do logito. O coeficiente assume a forma de razão de chances, onde uma chance é \(p/(1-p)\), e razão de chances é a divisão de duas chances (o que quer que isso signifique). Uma vantagem de trabalhar com razão de chances é que a interpretação do coeficiente se torna linear, em vez de não-linear, conforme mostra a equação abaixo.
\[ log(\frac{Pr(Y=1|X=x)}{Pr(Y=0|X=x)}) = \alpha + \beta x \]
Eu acho complicado entender o que é uma razão de chances, então preferimos trabalhar com probabilidade. Mas existe essa outra parametrização.
13.3 Logística como variável latente
Uma outra forma de interpretar ou justificar a logística é com a formulação de variável latente. Imagine que existe uma variável latente (não observada), \(z\), que é contínua e reflete a propensão a votar no governo. Porém, nós só observamos os valores \(1\) ou \(0\), de tal modo que podemos escrever:
\[\begin{align} y_i = 1, se & & z_i > 0 \\ y_i = 0, c.c. \\ z_i = \alpha + \beta x_i + e_i, e_i \sim logistica \end{align}\]
A distribuição logística é dada por: \[ Pr(e_i < x) = \frac{\exp(x)}{1 + \exp(x)} \]
Portanto, temos que:
\[\begin{align} Pr(y_i = 1) = Pr(z_i > 0) = \\ Pr(e_i + \alpha + \beta x_i > 0) = \\ Pr(e_i > -\alpha - \beta x_i) \\ \textrm{Pela regra do complemento:} \\ Pr(e_i > -\alpha - \beta x_i) = 1 - Pr(e_i < -\alpha - \beta x_i) = \\ 1 - \frac{\exp(- \alpha - \beta x_i)}{1 + \exp(- \alpha - \beta x_i)} = \\ \frac{1 + \exp(- \alpha - \beta x_i) - \exp(- \alpha - \beta x_i)}{1 + \exp(- \alpha - \beta x_i)} = \\ \frac{1}{1 + \exp(- \alpha - \beta x_i)} \tag{13.2} \end{align}\]
Vimos assim que as parametrizações de (13.1) e (13.2) são equivalentes.
13.4 Probit
Se, por outro lado, supusermos que o termo de erro tem distribuição normal, em vez de distribuição logística, teremos um modelo probit.
\[ z_i = \alpha + \beta x_i + e_i, e_i \sim N(0, \sigma^2) \] Quando \(\sigma = 1.6\), essa formulação é aproximadamente a regressão logística. Isso significa que o coeficiente da probit é aproximadamente o da logística/1.6 e vice-versa, isto é, podemos multiplicar o coeficiente da probit por 1.6. Vamos checar no R.
reg_logistica <- glm(votaria_governo ~ ideologia, data=br_latbar_23,
family=binomial(link= "logit"))
reg_probit <- glm(votaria_governo ~ ideologia, data=br_latbar_23,
family=binomial(link= "probit"))
library(stargazer)
stargazer(reg_logistica, reg_probit, type = "html")| Dependent variable: | ||
| votaria_governo | ||
| logistic | probit | |
| (1) | (2) | |
| ideologia | -0.285*** | -0.165*** |
| (0.039) | (0.021) | |
| Constant | 3.140*** | 1.836*** |
| (0.338) | (0.178) | |
| Observations | 472 | 472 |
| Log Likelihood | -239.376 | -239.128 |
| Akaike Inf. Crit. | 482.753 | 482.256 |
| Note: | p<0.1; p<0.05; p<0.01 | |
13.5 Logística como GLM
Uma outra forma de pensar a logística é como um modelo linear generalizado. Nós vimos que o modelo linear precisou entrar em uma função \(f\) que transformasse os intervalos entre \(0\) e \(1\). Podemos considerar o modelo de regressão linear como um caso particular, em que temos a função identidade \(f(x) = x\) e a distribuição da resposta é modelada como Gaussiana. Na logística,a função de ligação é logística, e a resposta modelada como Bernoulli (e para \(n\) dados, binomial).
Similarmente, se tenho dados discretos \(\{1, 2, 3, ..., n\}\), posso modelar meus dados com uma Poisson e com uma função \(f\) específica, que chamamos de função de ligação.
13.6 Logística com múltiplos preditores
reg_logistica1 <- glm(votaria_governo ~ ideologia + evangelico, data=br_latbar_23,
family=binomial(link= "logit"))
reg_logistica1 %>%
tidy() %>%
kable()| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 3.3035245 | 0.3498086 | 9.443805 | 0.0000000 |
| ideologia | -0.2780476 | 0.0391497 | -7.102157 | 0.0000000 |
| evangelico | -0.7387838 | 0.2432515 | -3.037120 | 0.0023885 |
Para facilitar a comparação de coeficientes, pode ser mais fácil padronizar as variáveis. No caso, vamos padronizar apenas ideologia (precisamos supor que é uma variável contínua, medida com valores discretos).
br_latbar_23 <- br_latbar_23 %>%
mutate(ideologia_pad = (ideologia - mean(ideologia))/sd(ideologia))
reg_logistica2 <- glm(votaria_governo ~ ideologia_pad + evangelico, data=br_latbar_23,
family=binomial(link= "logit"))
reg_logistica2 %>%
tidy() %>%
kable()| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 1.4449650 | 0.1497791 | 9.647304 | 0.0000000 |
| ideologia_pad | -0.9975056 | 0.1404511 | -7.102157 | 0.0000000 |
| evangelico | -0.7387838 | 0.2432515 | -3.037120 | 0.0023885 |
Usando a regra de dividir por quatro, temos que: ideologia tem efeito máximo de redzir a probabilidade em 25%, e evangélico de diminuir em 18%.
Plotando graficamente os dois preditores, temos:
jitter_binary <- function(a, jitt=0.05){
ifelse(a==0, runif(length(a), 0, jitt), runif(length(a), 1 - jitt, 1))
}
br_latbar_23$votaria_governo_jitter <- jitter_binary(br_latbar_23$votaria_governo)
ideologia_jitt <- br_latbar_23$ideologia_pad + runif(n, -.05, .05)
plot(ideologia_jitt, br_latbar_23$votaria_governo_jitter, xlim=c(0,max(ideologia_jitt)))
curve(invlogit(cbind(1, x, 0) %*% coef(reg_logistica2)), add=TRUE)
curve(invlogit(cbind(1, x, 1.0) %*% coef(reg_logistica2)), add=TRUE)
13.7 Ajuste do modelo
Nós vimos que podemos usar o erro quadrático médio com regressão linear para comparar ajustes do modelo. Porém, para variáveis binárias, há formas melhores de comparar a capacidade preditiva do modelo. A mais tradicional é o log da verossimilhança, dado pela soma para todos as observações de: \((ylog(p) + (1-y)log(1-p))\). No limite, um modelo perfeito, que classifica corretmente todas as observações (isto é, coloca probabilidade \(1\) de sucesso e \(0\) de fracasso) teria um log da verossimilhança igual a zero. Portanto, quanto mais próximo de zero, melhor o ajuste do modelo. E um modelo “burro”, que previsse probabilidade de 50% para cada ponto, teria um total de \(n \cdot log(.5) = n \cdot -0.693\). Um modelo um pouco melhor seria um que previsse a média, isto é, se 60% dos casos são sucesso, pode prever para todas as observações \(p=.6\). Esses modelos simples servem como comparação. No R:
y_hat <- predict(reg_logistica, type="response")
y_obs <- br_latbar_23$votaria_governo
log_likelihood <- sum(y_obs*log(y_hat) + (1-y_obs)*log(1 - y_hat))
print(log_likelihood)[1] -239.3764
[1] -274.8595
[1] -327.1655
| Dependent variable: | ||
| votaria_governo | ||
| (1) | (2) | |
| ideologia | -0.285*** | -0.278*** |
| (0.039) | (0.039) | |
| evangelico | -0.739*** | |
| (0.243) | ||
| Constant | 3.140*** | 3.304*** |
| (0.338) | (0.350) | |
| Observations | 472 | 472 |
| Log Likelihood | -239.376 | -234.816 |
| Akaike Inf. Crit. | 482.753 | 475.632 |
| Note: | p<0.1; p<0.05; p<0.01 | |
Para mais detalhes sobre regressão logística, nossa referência favorita é o livro de Gelman, Hill e Vehtati, Regression and other Stories, de onde tiramos alguns dos códigos do presente capítulo. Capítulos 13 e 14 contém bastante material a um nível que exige um nível de matemática similar ao que temos usado em nosso livro.