Modelos Lineares Generalizados

Quando a resposta é uma probabilidade ou uma contagem

Modelos Lineares
GLM
Logística
Poisson
Inferência
Uma introdução intuitiva a GLMs binomiais e Poisson, mantendo a lógica de modelos, contrastes e comparações.
Autores
Afiliações

João O. Santos

ISPA

Cristina Mendonça

ISPA; WJCR

Vitória Melita

FP-UL

Mariona Pascual Peñas

UIB

João Raposo

ISPA

Marta Barros

FP-UL

0 A Mesma Pergunta, Uma Resposta Diferente

No capítulo de comparação de modelos começámos com uma ideia muito simples:

\[ \text{Dados}=\text{Modelo}+\text{Erro}. \]

Depois vimos que uma média, um teste-t, uma regressão e uma ANOVA podiam ser formas diferentes de construir e comparar modelos lineares.

Agora imaginem duas novas perguntas:

A probabilidade de um organismo sobreviver é diferente entre duas condições?

O número esperado de organismos num quadrado depende do habitat?

As perguntas continuam a pedir diferenças, declives, contrastes e interacções. O problema está na resposta.

Sobrevivência só pode ser 0 ou 1. Uma probabilidade só pode ficar entre 0 e 1. Uma contagem pode ser 0, 1, 2, 3, … mas não -2.7.

Se usarmos uma recta como se nada tivesse mudado, o modelo pode fazer previsões impossíveis. E há uma segunda dificuldade ainda mais importante: nestes dados, a quantidade de variabilidade que esperamos costuma depender da própria média.

Os modelos lineares generalizados, ou GLMs, deixam-nos conservar a parte do modelo que representa a pergunta científica e mudar a forma como representamos a resposta e o erro.

A frase a reter é:

não estamos a abandonar o modelo linear; estamos a pôr o seu preditor linear numa escala onde ele faça sentido para a resposta que temos.

Primeiro Exemplo: Sobreviveu ou Não?

Suponham que codificamos:

  • morreu = 0;
  • sobreviveu = 1.

Se 20 de 100 organismos sobreviverem, a média de 0 e 1 é 0.20. Portanto, a média desta variável binária é também uma proporção.

Até aqui tudo parece familiar.

Agora suponham que ajustamos simplesmente uma recta. Nada impede essa recta de prever -0.10 num extremo ou 1.08 no outro. O problema não é estético: esses valores não são probabilidades.

Também sabemos que uma variável 0/1 não tem a mesma variabilidade em todo o lado. Quando a probabilidade está muito perto de 0 ou 1, há menos variabilidade possível do que quando está perto de 0.5.

Precisamos então de duas coisas:

  1. uma maneira de respeitar os valores possíveis da resposta;
  2. uma maneira de representar a variabilidade própria desse tipo de dados.

É isso que a família e a função de ligação acrescentam.

O Que Fica Igual?

Antes da terminologia nova, vale a pena olhar para a parte que não mudou.

O preditor linear continua a ser:

\[ \eta = \beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_1X_2+\ldots \]

Aqui continuam a viver:

  • interceptos;
  • declives;
  • contrastes;
  • factores;
  • interacções;
  • modelos compactos e aumentados.

Se já perceberam como uma coluna -1/+1 entra num modelo linear, já perceberam uma grande parte de como entra numa regressão logística ou Poisson.

A diferença é que () deixa de ser necessariamente a média da resposta na escala em que a medimos.

2 Respostas Binárias

Binomial: Manter Probabilidades Entre 0 e 1

Na regressão logística, a ligação habitual é o logit:

\[ \log\left(\frac{p}{1-p}\right)=\eta. \]

Se esta expressão parece hostil, não comecem pelo logaritmo. Comecem pela função inversa:

\[ p=\frac{1}{1+e^{-\eta}}. \]

() pode ir de menos infinito a mais infinito. Depois da transformação, a previsão (p) fica sempre entre 0 e 1.

É esse o trabalho essencial da ligação logística.

NotaDois Exemplos da Psicologia

Memória. Num ensaio de reconhecimento, a resposta pode ser codificada como correcta = 1 e incorrecta = 0. Um GLM binomial pode perguntar se a probabilidade de acertar muda com o atraso, a condição experimental ou a idade.

Controlo inibitório. Numa tarefa go/no-go, podemos modelar se cada resposta foi correcta ou incorrecta e perguntar se a probabilidade de erro muda com a fadiga ou com o tempo na tarefa.

Há, porém, uma pergunta diferente. Num teste de reconhecimento, duas pessoas podem ter a mesma percentagem de acertos e tendências muito diferentes para responder «sim, já vi». Quando queremos separar sensibilidade de critério de resposta, entramos na Teoria da Detecção de Sinal.

Os Contrastes Continuam Iguais

Suponham duas condições codificadas como -1 e +1:

Condição (X)
A -1
B +1

O preditor linear é:

\[ \eta=\beta_0+\beta_1X. \]

Substituímos os dois valores exactamente como fizemos no modelo linear.

Para A:

\[ \eta_A=\beta_0-\beta_1. \]

Para B:

\[ \eta_B=\beta_0+\beta_1. \]

Logo,

\[ \eta_B-\eta_A=2\beta_1. \]

A estrutura do contraste não mudou. O que mudou foi a escala em que essa diferença vive.

Num modelo Gaussiano com ligação identidade, esta diferença podia estar directamente em gramas, milissegundos ou outra unidade da resposta.

Na regressão logística, ela está em log-odds. Depois traduzimos as duas previsões para probabilidades para responder à pergunta substantiva.

É por isso que vale a pena separar duas tarefas:

  1. o contraste diz que comparação estamos a fazer;
  2. a ligação diz em que escala o modelo a representa.

Odds Sem Misticismo

As odds são apenas outra maneira de escrever uma probabilidade:

\[ odds=\frac{p}{1-p}. \]

Se (p=0.5), as odds são 1: há uma possibilidade de sucesso por cada possibilidade de insucesso.

Se (p=0.8),

\[ odds=\frac{0.8}{0.2}=4. \]

Há quatro sucessos esperados por cada insucesso, na linguagem das odds.

Um coeficiente logístico é uma mudança em log-odds por uma unidade do preditor. Se o exponenciarmos, obtemos a razão de odds correspondente a essa mudança de uma unidade:

$ OR=e^. $

A codificação importa. No contraste -1/+1 que acabámos de usar, passar do grupo A (-1) para o grupo B (+1) é uma mudança de duas unidades. Portanto, a razão de odds B/A é

$ OR_{B/A}=e^{2_1}. $

Isto é exactamente paralelo ao facto de a diferença B - A ser (2_1) na escala do preditor linear.

Uma razão de odds de 2 não quer dizer «mais 2 pontos percentuais» nem «o dobro da probabilidade». Por isso, mesmo quando reportamos odds ratios, costuma ser útil mostrar também probabilidades previstas.

O logit transforma o intervalo aberto (0<p<1) em toda a recta real. Além disso, uma diferença aditiva em log-odds corresponde a uma razão multiplicativa de odds.

Existem outras ligações binomiais, como probit e complementary log-log. A logística é muito usada porque a transformação é conveniente, simétrica e produz odds ratios com uma interpretação estabelecida. Não é a única ligação matematicamente possível.

Interacções Continuam a Ser Produtos

Este ponto merece ficar no caminho principal.

Suponham dois factores A e B codificados como -1/+1:

\[ \eta = \beta_0+\beta_AA+\beta_BB+\beta_{AB}AB. \]

A coluna da interacção continua a ser o produto das outras duas.

A B (AB)
-1 -1 +1
+1 -1 -1
-1 +1 -1
+1 +1 +1

E (_{AB}) continua a representar uma diferença de diferenças na escala do preditor linear.

Com esta codificação,

\[ 4\beta_{AB} \]

corresponde à diferença de diferenças em ().

A novidade é importante: numa regressão logística, essa diferença de diferenças é em log-odds, não necessariamente em probabilidades.

Duas linhas que parecem ter uma interacção forte na escala da probabilidade podem ter outro aspecto na escala logit, e vice-versa. Portanto, depois de testar a interacção, olhem sempre para as previsões na escala da resposta.

3 Contagens e Taxas

Poisson: Quando a Resposta É uma Contagem

Agora mudemos a pergunta.

Quantos organismos esperamos encontrar num quadrado, e esse número depende do habitat?

Uma contagem não pode ser negativa. Por isso, um modelo Poisson usa habitualmente a ligação log:

\[ \log(\mu)=\eta. \]

Voltando à escala da resposta:

\[ \mu=e^\eta. \]

Como (e^) é sempre positivo, a contagem média prevista também é.

A Interpretação Torna-se Multiplicativa

Suponham um preditor quantitativo:

\[ \eta=\beta_0+\beta_1X. \]

Aumentar (X) uma unidade acrescenta (_1) ao log da média.

Na escala original, isso significa multiplicar a média esperada por:

\[ e^{\beta_1}. \]

Se

\[ e^{\beta_1}=2, \]

a contagem esperada duplica.

Se

$ e^{_1}=0.8, $

a contagem esperada passa para 80% do valor anterior.

Uma palavra importante aqui é factor multiplicativo. Em inglês, esta interpretação é frequentemente descrita como multiplicative change. Quando comparamos duas médias previstas, podemos também falar numa mean ratio, uma razão entre médias esperadas:

$ =e^{_1}. $

Não é, em geral, uma rate of change: essa expressão sugere uma velocidade de mudança ou uma derivada. (e^{_1}) diz por que factor a média esperada é multiplicada quando (X) aumenta uma unidade. Só quando a resposta representa uma taxa por unidade de exposição, como veremos no offset, faz sentido falar numa rate ratio, ou razão de taxas.

O efeito também se vê melhor do que se lê. Se (e^{_1}=1.5), cada passo de uma unidade em (X) multiplica a média anterior por 1.5:

~{r} #| echo: false #| fig-alt: “Curva crescente da média esperada num modelo com ligação log em que cada aumento de uma unidade no preditor multiplica a média por 1.5.”

x_log <- 0:8 mu_log <- 2 * 1.5^x_log

ggplot2::ggplot( data.frame(x = x_log, media = mu_log), ggplot2::aes(x = x, y = media) ) + ggplot2::geom_line(linewidth = 0.9) + ggplot2::geom_point(size = 2) + ggplot2::theme_classic() + ggplot2::labs( x = “X”, y = “Média esperada” ) ~

Na escala do logaritmo, os passos são aditivos. Na escala da resposta, tornam-se multiplicativos. É esta tradução que produz a curva.

Não é uma nova lógica de regressão. É o mesmo declive, traduzido através de uma ligação diferente.

Família Também Quer Dizer Variância

Esta é a parte que convém não perder atrás das palavras «binomial» e «Poisson».

Num modelo Poisson simples:

\[ Var(Y\mid X)=E(Y\mid X)=\mu. \]

Isto significa que, quando a média aumenta, a variabilidade esperada também aumenta.

Uma zona onde esperamos 100 organismos pode variar muito mais, em números absolutos, do que uma zona onde esperamos 2.

No binomial, a variância depende da probabilidade e do número de tentativas.

Por isso, escolher uma família não decide apenas que valores o modelo pode prever. Decide também quanto erro o modelo espera encontrar em diferentes partes da escala.

E isso afecta erros-padrão, intervalos de confiança e testes.

Sobredispersão: Primeiro Perguntar Porquê

Num Poisson simples, a variância esperada está presa à média. Dizemos que há sobredispersão quando os dados apresentam mais variabilidade do que essa estrutura prevê.

Mas «há sobredispersão» ainda não é um diagnóstico substantivo.

Pode acontecer porque:

  • falta um preditor importante;
  • a relação média–preditor está mal especificada;
  • há observações repetidas da mesma unidade;
  • locais ou organismos diferem entre si de uma forma não representada;
  • a exposição varia;
  • há uma estrutura de contagens mais variável do que o Poisson;
  • os zeros podem nascer de um processo adicional.

Essas histórias pedem soluções diferentes.

Se a causa é dependência dentro de tanques, trocar Poisson por uma distribuição mais flexível não cria tanques independentes.

Se a causa é heterogeneidade de contagens que cresce mais depressa do que a média, uma binomial negativa pode ser uma representação mais adequada.

Se a causa é uma variável explicativa em falta, talvez o primeiro passo seja melhorar o modelo da média.

AvisoNão Tratem Sobredispersão Como um Botão

A sequência «Poisson → detectei sobredispersão → binomial negativa» pode ser uma boa decisão em alguns problemas. Não é uma regra universal.

O diagnóstico útil pergunta que processo poderia produzir a variabilidade extra. A mudança de modelo deve responder a essa história.

Alternativas Comuns

Binomial negativa: acrescenta um parâmetro de dispersão e permite uma relação média–variância mais flexível.

Quasi-Poisson: mantém a estrutura da média do Poisson e introduz um parâmetro de escala para a variância. Como não define a mesma verosimilhança completa, não entra da mesma forma em comparações por AIC.

Efeito aleatório: representa heterogeneidade ou dependência partilhada por uma unidade de agrupamento.

Modelos com processo extra de zeros: representam explicitamente uma estrutura adicional para produzir zeros. Devem ser usados quando essa história faz sentido para o processo, não apenas porque o histograma tem muitos zeros.

Contagem ou Taxa?

Imaginem dois transectos.

Num deles observámos durante 1 hora. No outro, durante 10 horas.

Se ambos tiverem 12 avistamentos, têm a mesma contagem, mas claramente não a mesma taxa de avistamentos.

Num modelo Poisson com ligação log podemos incluir uma exposição (E_i):

\[ \log(\mu_i) = \eta_i+\log(E_i). \]

O termo ((E_i)) entra como offset, com coeficiente fixado em 1.

Em R, a estrutura seria:

count ~ habitat + offset(log(hours))

Sem o offset perguntamos por contagens totais.

Com o offset perguntamos por uma taxa por unidade de exposição.

Isto não é uma correcção técnica invisível. Muda a pergunta científica.

4 Inferência e Diagnóstico

Comparar Modelos Outra Vez

No modelo linear, comparámos modelos perguntando quanto diminuía o erro quadrático.

Num GLM continuamos a poder construir:

  • um modelo compacto, que impõe a restrição de (H_0);
  • um modelo aumentado, que liberta o parâmetro ou conjunto de parâmetros que representa a alternativa.

A medida de ajuste já não é simplesmente o SSE. A família dá-nos uma verosimilhança: uma medida de quão compatíveis são os dados observados com um conjunto de parâmetros dentro daquele modelo probabilístico.

A deviance traduz essa informação para uma escala conveniente. Para modelos aninhados, uma redução de deviance significa que o modelo aumentado se ajusta melhor.

Podemos resumir a comparação como:

\[ \Delta D=D_{compacto}-D_{aumentado}. \]

Quanto maior a redução, mais o modelo aumentado ganhou ao libertar a restrição.

Em muitos GLMs, sob condições regulares, comparamos (D) com uma distribuição qui-quadrado.

O valor-p mantém exactamente a interpretação que introduzimos em Do t ao Valor-p. Aqui mudou a estatística de teste — usamos a melhoria de verosimilhança em vez de um t ou um F Gaussiano — mas não mudou a lógica da inferência.

Se (L()) é a verosimilhança e (()=L()), o teste de razão de verosimilhanças para dois modelos aninhados usa

\[ 2\left[ \ell(\hat\theta_A)-\ell(\hat\theta_C) \right]. \]

Sob condições regulares, esta quantidade tem aproximadamente uma distribuição (^2), com graus de liberdade iguais à diferença no número de parâmetros.

A deviance é definida por comparação com um modelo saturado e, para esta comparação aninhada, a diferença de deviances corresponde à mesma estatística de razão de verosimilhanças.

A aproximação qui-quadrado tem condições; não é uma licença para comparar quaisquer dois modelos de qualquer forma.

Como o Computador Ajusta um GLM?

Esta parte é importante para quem quer perceber a computação, mas não para interpretar o primeiro modelo.

No modelo linear Gaussiano, os mínimos quadrados levam a uma solução fechada para os coeficientes quando a matriz do modelo tem posto completo.

Num GLM, em geral, não existe uma fórmula única que resolva directamente (). O computador escolhe os coeficientes que maximizam a verosimilhança.

Uma estratégia clássica é iteratively reweighted least squares (IRLS), intimamente relacionada com Newton–Raphson/Fisher scoring. Em cada iteração, o problema é aproximado por um problema de mínimos quadrados ponderados; os pesos e a resposta de trabalho são actualizados até o algoritmo convergir.

Esta maquinaria explica como glm() encontra as estimativas. Não muda o que um contraste -1/+1, um declive ou uma interacção significam.

Contrastes e Tipo III Continuam Cá

A parte (X) continua a ser um modelo linear. Portanto, a discussão sobre contrastes e Tipo III continua relevante.

Se A e B estão codificados com contrastes centrados e o modelo contém A × B, podemos continuar a perguntar:

  • diferença média associada a A, mantendo os restantes termos no modelo;
  • diferença média associada a B;
  • diferença de diferenças A × B.

A diferença é que estes efeitos são definidos primeiro na escala da ligação.

Por isso, uma boa análise não termina na tabela omnibus. Depois do teste, voltamos à escala que interessa:

  • probabilidades;
  • diferenças de probabilidades;
  • contagens esperadas;
  • razões entre médias esperadas;
  • razões de taxas, quando existe uma exposição que define uma taxa;
  • gráficos de previsões.

É aí que vemos o tamanho e a direcção substantiva do efeito.

Diagnóstico: Qual Parte da História Falhou?

Num GLM, «o modelo não ajusta» é demasiado vago.

Tentem localizar o problema.

Média: falta curvatura, um preditor ou uma interacção?

Variância/família: os dados variam como a família prevê?

Ligação: a forma que liga os preditores à média parece plausível?

Dependência: várias linhas vêm do mesmo animal, local, tanque ou participante?

Exposição: algumas unidades tiveram mais oportunidade de produzir eventos?

Influência: algumas observações estão a determinar quase sozinhas o resultado?

Binomial: há separação ou previsões extremas sustentadas por poucos casos?

Contagens: a sobredispersão ou os zeros sugerem um processo que o modelo não representa?

Cada pergunta aponta para uma modificação diferente. O objectivo do diagnóstico não é fazer desaparecer uma mensagem do software. É descobrir que afirmação do modelo deixou de ser plausível.

5 Escolher e Ler o Modelo

Um Mapa Para Escolher

Se a resposta é… Um primeiro modelo possível O que estamos a modelar
contínua aproximadamente Gaussiana modelo linear média na escala original
0/1 binomial + logit probabilidade, via log-odds
sucessos em (n) tentativas binomial + logit proporção/probabilidade
contagem Poisson + log contagem média
contagem por exposição Poisson + log + offset taxa
contagem com variância extra plausível talvez binomial negativa média com dispersão adicional

A tabela é um ponto de partida, não uma árvore de decisão automática.

Um Fluxo de Leitura

Quando encontrarem um GLM num artigo, ou quando tiverem de escolher um:

  1. Qual é a pergunta de investigação?
  2. O que representa uma linha?
  3. Que valores pode assumir a resposta?
  4. Como esperamos que a variabilidade mude com a média?
  5. Que contraste, declive ou interacção representa a hipótese?
  6. Que ligação transforma a escala da resposta numa escala onde esses efeitos se possam combinar?
  7. Há exposição diferente entre observações?
  8. Há dependência ou agrupamento?
  9. O que dizem as previsões quando voltamos à escala da resposta?
  10. Se o diagnóstico falhar, que parte da história precisa realmente de mudar?

Se conseguirem responder a estas perguntas, já estão a interpretar o modelo em vez de apenas chamar glm().

6 Síntese

Os GLMs parecem trazer muita terminologia nova, mas a espinha dorsal é a mesma.

A pergunta científica continua a ser representada por:

\[ \eta=X\beta. \]

É aí que continuam os contrastes -1/+1, os declives, os produtos que criam interacções e os modelos compactos e aumentados.

A família diz que tipo de resposta e de variabilidade estamos preparados para observar.

A ligação cria uma escala onde o preditor linear pode trabalhar sem produzir probabilidades ou médias impossíveis.

A comparação de modelos continua a perguntar o que ganhámos ao libertar uma restrição. A interpretação do valor-p é a mesma que já aprendemos no capítulo anterior.

O que precisamos de reaprender não é «uma nova estatística para cada caso». Precisamos de aprender a perguntar:

Que parte da história dos dados mudou, e que peça do modelo deve mudar com ela?

7 Exercícios

  1. Uma resposta 0/1 tem média 0.30. Explique por que essa média também pode ser lida como uma proporção.
  2. Dois grupos estão codificados -1/+1 num modelo logístico. Derive (_A), (_B) e mostre por que a diferença é (2_1).
  3. Num 2 × 2 com contrastes -1/+1, construa a coluna produto A × B e explique por que a interacção continua a ser uma diferença de diferenças.
  4. Explique por que uma diferença constante em log-odds não implica uma diferença constante em pontos percentuais.
  5. Num Poisson, (e^). Escreva a interpretação na escala da contagem esperada.
  6. Dois locais tiveram 12 avistamentos; um foi observado 2 horas e outro 12. Explique a diferença entre modelar contagem e modelar taxa.
  7. Um Poisson apresenta forte sobredispersão. Dê três processos diferentes que a poderiam produzir e diga que mudança de modelo investigaria em cada caso.
  8. Explique por que uma comparação por diferença de deviance continua a ser uma comparação entre um modelo compacto e um aumentado, apesar de já não usarmos a soma de quadrados como medida de ajuste.
  9. Diga por que um efeito aleatório e uma binomial negativa não são duas maneiras intercambiáveis de «corrigir sobredispersão».

8 Recursos