Modelos lineares e variáveis binárias

Embora os modelos lineares nos dêem a base para entendermos os fundamentos da análise de regressão, estes modelos não são os mais comuns nas Ciências Sociais. Se nossas observações se dão sobre a ocorrência ou nao de um golpe de estado, a entrada ou não de mulheres no mercado de trabalho, a evasão ou não de um aluno do ensino médio etc. estamos lidando com variáveis que assumem apenas dois estados, isto é, variáveis dicotômicas. Neste caso um modelo linear não é o mais adequado para estimarmos o efeito de possíveis fatores explicativos destas observações.

Para entender o por quê disso vamos seguir o exemplo de Fox(1997:439) e analisar a votação do plebiscito que restaurou a democracia no Chile em 1988. Utilizaremos um modelo linear de regressão para mostrar os problemas que surgem com a adoçao deste tipo de modelo e como a solução destes problemas nos levam à adoção de modelos como logit e probit.

Os dados foram obtidos por meio de um Survey nacional conduzido pela FLACSO/Chile envolvendo 2700 eleitores chilenos selecionados aleatoriamente. Abaixo vemos um resumo dos dados:

kable(summary(Chile)) %>% kable_classic(full_width = T, position = "center", )
region population sex age education income statusquo vote
C :600 Min. : 3750 F:1379 Min. :18.00 P :1107 Min. : 2500 Min. :-1.80301 A :187
M :100 1st Qu.: 25000 M:1321 1st Qu.:26.00 PS : 462 1st Qu.: 7500 1st Qu.:-1.00223 N :889
N :322 Median :175000 NA Median :36.00 S :1120 Median : 15000 Median :-0.04558 U :588
S :718 Mean :152222 NA Mean :38.55 NA’s: 11 Mean : 33876 Mean : 0.00000 Y :868
SA:960 3rd Qu.:250000 NA 3rd Qu.:49.00 NA 3rd Qu.: 35000 3rd Qu.: 0.96857 NA’s:168
NA Max. :250000 NA Max. :70.00 NA Max. :200000 Max. : 2.04859 NA
NA NA NA NA’s :1 NA NA’s :98 NA’s :17 NA

Vamos observar apenas os eleitores que expressaram uma preferência (Si (“Y”), a favor da continuidade do regime militar ou No (“N”), a favor do retorno a um governo civil). Para isso vamos criar uma nova variável chamada vote.dico, com os valores 1 para o voto “Y” e 0 para “N”. Esse tipo de codificação é o mais utiizado na representação de variáveis dicotômicas. Em seguida vamos reproduzir a Figura 15.1 de Fox(1997). Nesta figura procuramos entender a relação do voto com o apoio ao status quo, uma escala formada por respostas a questões sobre o apoio a medidas sociais, econômicas e políticas do governo de Pinochet e onde valores altos significam maior apoio.

Chile$voto.dico <-  ifelse(Chile$vote == "Y",1, ifelse(Chile$vote == "N",0,NA))
kable(table(Chile$voto.dico)) %>% kable_classic(full_width = F, position = "center", )
Var1 Freq
0 889
1 868 
g1 <- ggplot(Chile, aes(statusquo, voto.dico)) + geom_point()
g1

g1 <- ggplot(Chile, aes(statusquo, voto.dico)) + geom_jitter(width = 0.1, height = 0.1)
g1

Na primeria figura os pontos se sobrepuseram dificultando a interpretação. Para lidar com isso utilizei a função geom_jitter do ggplot (equivalente ao “jitter” do pacote base) para espalhar um pouco os dados e observar melhor sua distribuição. Essa função afeta apenas o gráfico e não os dados.

Costumamos pensar uma regressão linear como a média condicional. Mas como vemos nos gráficos a variável \(Y\) assume os valores zero e um. Qual o sentido de se falar em média? Na população a média condicional é a proporção de 1’s entre os indivíduos que compartilham a variável \(x\), a quantidade de eleitores que votam “Sí” em um certo grupo. Ou seja:

\[\pi_{i} \equiv Pr(Y_{i}) \equiv Pr(Y_{i} = 1|X = x_{i})\]

e

\[E(Y|x_{i}) = \pi_{i}(1) + (1 - \pi_{i})(0) = \pi_{i}\] (1)

Como podemos calcular \(E(Y|x_{i})\)? Se \(X\) é discreta (qualitativa) podemos contar a proporção de 1’s em cada categoria de \(X\). Se, por exemplo, quisermos saber a média de votos no “Sí” em cada região do Chile, podemos observar a tabela cruzada abaixo onde as proporções nos são dadas diretamente.

kable(prop.table(xtabs(~Chile$voto.dico + Chile$region),2)) %>% kable_classic(full_width = F, position = "center", )
C M N S SA
0 0.546875 0.3214286 0.4303797 0.4376278 0.5837563
1 0.453125 0.6785714 0.5696203 0.5623722 0.4162437

No caso de uma variável dependente contínua, como é o caso do apoio ao status quo, essa estratégia não funcionaria, pois não temos categorias naturais onde observar a quantidade de “Sis”. Com uma amostra grande podemos utilizar alguma técnica de regressão não paramétrica, como regressão local, para “discretizar” nossa variável observando a proporção de 1’s para diferentes faixas de valores de \(X\).

A figura abaixo mostra uma linha de regressão local (“loess”) sobreposta ao gráfico da intenção de votos pela escala de atitude diante do status quo. Vemos que nos níveis mais baixos de apoio ao status quo a proporção de “Sí” fica perto de 0 e que nos níveis mais altos de apoio ela fica perto de 1 formando um “S, o que se ajusta bem aos dados”:

g1 + geom_smooth(method = "loess")
## `geom_smooth()` using formula 'y ~ x'

Embora o ajuste da “regressão local” seja bom, nem sempre é possível fazer isso e, mais importante, é difícil interpretar o que essa regresão nos diz sobre a relação entre o apoio ao status quo e o voto no “si”, por exemplo, o quanto a mudança em uma unidade de apoio ao SQ impacta a probabilidade do eleitor votar “si”. O que queremos é um modelo que explique nossa incerteza sobre \(Y\) (voto “Sí”) em função de \(X\) (“apoio ao Status Quo”).

Para construir esse modelo, vamos inicialmente tentar uma distribuição gaussiana (como a normal) para \(Y\) onde a relação entre a média de \(Y\) e \(X\) é linear nos parâmetros, isto é:

\[Y_{i} = \alpha + \beta X_{i} + \epsilon_{i}\] (2)

onde \(\epsilon_{i} \sim N(0,\sigma_{\epsilon}^2)\), \(\epsilon_{i} \perp \epsilon_{j}\) para todo \(i \neq j\) e, quando \(X\) é aleatório, \(X \perp \epsilon\).

Usando as equações (1) e (2) temos:

\[E(Y_{i}) = \alpha + \beta X_{i} = \pi_{i}\]

Por conta da identidade \(E(Y_{i}) = \pi_{i}\) o modelo linear aplicado a uma variável dependente dicotômica é chamada de modelo linear de probabilidade. O que acontece se aplicamos esse modelo aos dados do plebiscito?

chile.mlp <- lm(voto.dico ~ statusquo, data = Chile)
par(mfrow = c(2,2))
plot(chile.mlp)

Compare essa figura com a figura abaixo com os gráficos de diagnóstico de uma regressão que atende às premissas de um modelo linear:

X <- rnorm(100,10,2)
Y = X + rnorm(100)
ok.lm <- lm(Y ~ X)
par(mfrow = c(2,2))
plot(ok.lm)

O gráfico no canto superio esquerdo das figuras mostra a relação entre os resíduos \(\epsilon_{i}\) e os valores ajustados \(\hat{Y} = \alpha + \beta X\). Se a relação fosse linear (como é o caso da simulação) não esperaríamos encontrar nenhum padrão discernível. Já os resíduos da regressão do voto no status quo seguem um padrão bem discernível: se \(Y_{i} = 1\), o que ocorre com probabilidade \(\pi_{i}\) temos:

\[\epsilon_{i} = 1 - E(Y_{i}) = 1 - (\alpha + \beta X_{i}) = 1 - \pi\]

Se $Y_{i} = 0 $, o que ocorre com probabilidade \(1 - \pi_{i}\) temos:

\[\epsilon_{i} = 0 - E(Y_{i}) = 0 - (\alpha + \beta X_{i}) = - \pi\]

O gráfico no canto superior direito das figuras é um Q-Q Plot e indica se os resíduos seguem uma distribuição normal. Se os resíduos caem sobre a linha pontilhada a dsitribuição se aproxima da normal. Vemos claramente que no caso dos resíduos da regressão do voto no status quo a distribuição do erro não é a normal.

O gráfico no canto inferior esquerdo da figura mostra se os resíduos se espalham homogeneamente por diferentes faixas de valores de \(E(Y_{i})\). Se houver homogeneidade isso quer dizer que a variância é constante, como exige a premissa de homoscedasticidade. No caso da regresssão do voto no satus quo, por usarmos um modelo linear assumimos que \(E(\epsilon) = 0\) e combinando o que vimos acima a variância de \(\epsilon_{i}\) é:

\[V(\epsilon_{i}) = \pi_{i}(1 - \pi_{i})^2 + (1 - \pi_{i})(- \pi_{i})^2 = \pi_{i}(1-\pi_{i})\]

Isto é, a variância dos resíduos da regressão de uma variável dicotômica depende de \(\pi_{i}\) que, por sua vez, varia com \(X_{i}\) e, portanto, não pode ser constante, o gráfico abaixo mstra a curva da variância calculada acima:

curve(x*(1-x))

O problema mais sério com o modelo linear de probabilidade aparece no Q-Q Plot da regressão do voto no status quo. Vemos que poucos valores se encontram próximos a \(E(\epsilon_{i}) = 0\) e que a maior parte dos valores se distanciam rapidamente de 0. Isso quer dizer que para valores extremos de \(X\) os valores de \(\pi\) podem estar fora do intervalo [0,1], o que não faz sentido quando interpretamos probabilidades.

Se uma grande parte dos valores de \(X\) se concentram próximos a 0.5, o modelo linear de probabilidade vai apresentar \(\pi\) dentro do intervalo, mas esse não é o caso do nosso exemplo o que pode ser visto quando sobrepomos a linha de mínimos quadrados no gráfico da relação entre voto e status quo.

ggplot(Chile, aes(statusquo, voto.dico)) + geom_point() + geom_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

Modelos logit e probit

Se os modelos lineares podem gerar valores fora do intervalo [0,1] a solução seria encontrar algum modelo que resolva esses problemas. Uma alternativa é a função retangular do tipo:

\[\begin{equation} \pi = \begin{cases} 0 & : 0 > \alpha + \beta X \\ \alpha + \beta X & : 0 \leq \alpha + \beta X \\ 0 & : \alpha + \beta X > 1 \end{cases} \end{equation}\]

Neste caso teríamos a seguinte curva:

curve(punif(x), -4, 4, col = "blue")

Embora a distribuição retangular esteja dentro do intervalo [0,1] ela tem alguns incovenientes, o principal é que a inclinação é muito abrupta, enquanto que, como vimos acima, no caso da curva não paramétrica, a relação entre \(\pi\) e \(X\) pode ser mais suave. Então além de uma distribuição que fique no intervalo [0,1], também gostarámos que essa distribuição fosse suave. As distribuições de probabilidade acumuladas são candidatas naturais e nosso modelo seria:

\[\pi_{i} = P(\eta_i) = P(\alpha + \beta X_{i})\]

Onde \(P(.)\) é uma função de distribuição acumulada (fda) e \(\alpha\) e \(\beta\) são parâmetros desta função a serem estimados. Se a função for estritamente crescente podemos inverte-la obtendo:

\[ P^{-1}(\pi_{i}) = \alpha + \beta X_{i} \]

Dese modo podemos ter um modelo linear de uma transformação de \(\pi\) (\(P^{-1}(.)\)) ou um modelo não linear de de \(\pi\) (\(P(.)\)).

Duas transformações frequentemente escolhidas são a função de distribuição acumulada da normal padronizada:

\[\Phi(z) = \frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{z}exp\Bigl(-\frac{1}{2}Z^{2}\Bigl)dZ\]

ou, mais usualmente, a distribuição logística:

\[\Lambda(z) = \frac{1}{1 + e^{-z}}\]

Graficamente temos:

curve(pnorm(x),-4,4, col = 2)
curve(plogis(x,0,0.6), -4, 4, add = T, col = "green")

Se substituirmos \(z\) por \(\alpha + \beta X_{i}\) na distribuição normal padronizada obtemos o modelo probit

\[\pi = \Phi(\alpha + \beta X_{i}) = \frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\alpha + \beta X_{i}}exp\Bigl(-\frac{1}{2}(\alpha + \beta X_{i})^{2}\Bigl)d(\alpha + \beta X_{i})\]

e fazendo a mesma substituição na distribuição logística temos o modelo

\[\pi = \Lambda(\alpha + \beta X_{i}) = \frac{1}{1 + e^{-(\alpha + \beta X_{i})}}\]

Avaliar o probit para diversos valores de \(z_{i} = \alpha + \beta X_{i}\) exige o cálculo de uma integral. A função logística é mais fácil de ser calculada, o que faz com que seja mais utilizada.

Outra propriedade útil destas distribuições é que elas são simétricas e monotônicas (uma mudança de escala mantém a ordem entre os valores) o que permite ‘reverter’ a função.

A inversa da função logística pode ser interpretada como o logaritmo da chance de P(\(Y_{i} = 1\)). Rearranjando a função logística temos

\[\begin{align} \pi &= \frac{1}{1 + e^{-(\alpha + \beta X_{i})}} \\ \pi(1 + e^{-(\alpha + \beta X_{i})}) &= 1 \\ \pi + \pi e^{-(\alpha + \beta X_{i})} &= 1 \\ \pi e^{-(\alpha + \beta X_{i})} &= 1 -\pi \\ e^{-(\alpha + \beta X_{i})} &= \frac{1 -\pi}{\pi} \\ \frac{1}{e^{\alpha + \beta X_{i}}} &= \frac{1 -\pi}{\pi} \\ e^{\alpha + \beta X_{i}} &= \frac{\pi}{1-\pi} \end{align}\]

\(\pi/(1-\pi)\) é a chance de P(\(Y_{i} = 1\)). Tirando o logaritmo dos dois lados obtemos:

\[log\bigl(\frac{\pi}{1-\pi}\bigl) = \alpha + \beta X_{i}\]

Também conhecida como logit. Com essas duas fórmulas podemos interpretar o coeficiente \(\beta\) como sendo o aumento no logit \(log(\frac{\pi}{1-\pi})\) dado o incremento de uma unidade de \(X\).

Se invertemos o logit temos \(\text{logit}^{-1} = \text{exp}^{\alpha+\beta X}/(1+\text{exp}^{\alpha+\beta X}) = \pi\)

Podemos derivar uma interpretação mais direta observando que a relação entre \(X\) e \(\pi\) é quase linear em torno de \(\pi = 0.5\) e que a inclinação neste ponto é a derivada \(\beta\text{exp}^{\alpha+\beta X}/(1+\text{exp}^{\alpha+\beta X})^2\) quando \(\alpha + \beta X = 0\). Daí temos que \(\beta\text{exp}^0/(1+\text{exp}^{0})^2 = \beta\Big(\frac{1}{2^2}\Big) = \beta/4\). Esse valor pode ser interretado como o efeito máximo na \(P(Y = 1)\) dada a mudança em uma unidade de \(X\).

Os coeficientes \(\alpha\) e \(\beta\) são desconhecidos. Para estima-los utilizamos a máxima verossimilhança. Primeiro definimos o componente estocástico, nossa incerteza sobre \(Y\). No exemplo do plebiscito cada decisão “Si” ou “No” segue uma distribuição Bernoulli de probabilidade:

\[P(Y|\pi_{i}) = \pi^{y_{i}}(1-\pi)^{1 - y_{i}}\]

Ou seja um eleitor pode ou não votar sim. Se juntarmos todos os eleitores e suas probabilidades, partindo da premissa que todas as decisões são independentes:

\[P(Y|\pi) = \prod_{i=1}^{n}\pi^{y_{i}}(1-\pi)^{1- y_{i}}\]

Usando \(\pi = \frac{1}{1 + e^{-x_{i}\beta}}\), onde \(\beta\) é um vetor com todos os parâmetros \(\beta_{i}\) e a constante \(\alpha\), obtemos nosso componente sistemático e substituimos na distribuição acima:

\[\begin{align} P(Y|\pi) &= \prod_{i=1}^{n}\Big(\frac{1}{1 + e^{-x_{i}\beta}}\Big)^{y_{i}}\Big(1-\frac{1}{1 +e^{-x_{i}\beta}}\Big)^{1- y_{i}}\\ P(Y|\pi) &= \prod_{i=1}^{n}\Big(1 + e^{-x_{i}\beta}\Big)^{-y_{i}}\Big(1 + e^{x_{i}\beta}\Big)^{-(1- y_{i})} \end{align}\]

Como a verossimilhança é proporcional à essa probabilidade e tirando o logaritmo da equação acima para nos livrar do produto, obtendo o log da verossimilhança (log-likelihood):

\[\text{log}\mathcal{L}(\tilde\beta \mid y) = \sum_{i = 1}^{n}-y \text{ ln}[1 + e^{-x_{i}\beta}]-(1-y) \text{ ln}[1 + e^{x_{i}\beta}]\]

O log da verossimilhança é diferenciável e, igualando a derivativa a zero, podemos obter o valor de \(\beta\) que maximiza a função. Não há uma solução analítica para isso e, como veremos mais a fundo na próxima aula, usamos soluções numéricas para estimar \(\beta\).

Usando o R:

## glm(formula = voto.dico ~ statusquo, family = binomial(link = "logit"), 
##     data = Chile)
##             coef.est coef.se
## (Intercept) 0.22     0.10   
## statusquo   3.21     0.14   
## ---
##   n = 1754, k = 2
##   residual deviance = 752.6, null deviance = 2431.3 (difference = 1678.7)

Calculando a probabilidade quando X = 0, usando a inversa do logit \(\text{logit}^{-1}(0.21)\)

## [1] 0.5523079

Calculando a mudança na probabilidade quando aumentamos X em uma unidade (primeiras diferenças) ou \(\text{logit}^{-1}(0.21 + 3.2) - \text{logit}^{-1}(0.21)\):

## [1] 0.4157077
## lm(formula = voto.dico ~ statusquo, data = Chile)
##             coef.est coef.se
## (Intercept) 0.49     0.01   
## statusquo   0.39     0.01   
## ---
## n = 1754, k = 2
## residual sd = 0.26, R-Squared = 0.73

Logit na prática

(GELMAN & HILL, Chap. 5. poços em Bangladesh)