GLMs no R: modelos para proporções e contagens

Modelos Lineares

GLM

Logística

Poisson

emmeans

Diagnósticos

Regressão logística (binomial) e Poisson, interpretação com emmeans, comparação de modelos e diagnósticos.

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 Setup

library(dplyr)
library(tidyr)
library(ggplot2)
library(car)         # Anova(type = 3)
library(emmeans)     # médias marginais estimadas
library(performance) # diagnósticos práticos

options(contrasts = c("contr.sum", "contr.poly"))

A ideia em 1 frase

Um GLM é o lm() quando a variável dependente não é bem descrita por erros normais (p.ex., 0/1, proporções, contagens).

1 GLM, Quando, Como e Porquê?

Se a VD é quantitativa (aprox. contínua) e os resíduos fazem sentido como “ruído normal”, começamos por lm().

Mas se a VD for:

Binária (0/1; sim/não; sucesso/insucesso) → GLM binomial (regressão logística).
Proporção (k sucessos em n tentativas) → GLM binomial (com cbind(sucessos, falhas)), ou binom.test(), ou prop.test(), ou teste de independência do qui-quadrado com chisq.test().
Contagem (0, 1, 2, 3, …) → GLM Poisson (ou alternativas se houver overdispersion, ver pacote glmmTMB).

Cábula de família + link

Binomial: family = binomial(link = "logit")
Poisson: family = poisson(link = "log")

Na prática, quase nunca precisas de mudar o link por defeito (mas podes, se houver motivo).

2 Comparação de GLMs

A lógica que já vimos nos modelos lineares continua a funcionar:

Temos um modelo compacto (mais simples).
Temos um modelo aumentado (mais parâmetros)
Avalias se os ganhos no desempenho justificam o aumento de complexidade.

A diferença é que, em GLMs, o erro não é medido com SSE (soma de quadrados), mas com log-likelihood / deviance.

3 Binomial: Regressão Logística

Vamos prever a probabilidade de sobrevivência (Survived) em função da classe (Pclass) e do sexo (Sex).

ds <- read.csv("../data/titanic.csv")

# tratar a classe como factor nominal
# (muitas bases vêm com Pclass como número, mas aqui é qualitativa)
ds$Pclass <- factor(paste0("class", ds$Pclass))
ds$Sex <- factor(ds$Sex)

head(ds)

X	PassengerId	Survived	Pclass	Name	Sex	Age	SibSp	Ticket	Fare	Cabin	Embarked
1	1	0	class3	Braund, Mr. Owen Harris	male	22	1	A/5 21171	7.2500		S
2	2	1	class1	Cumings, Mrs. John Bradley (Florence Briggs Thayer)	female	38	1	PC 17599	71.2833	C85	C
3	3	1	class3	Heikkinen, Miss. Laina	female	26	0	STON/O2. 3101282	7.9250		S
4	4	1	class1	Futrelle, Mrs. Jacques Heath (Lily May Peel)	female	35	1	113803	53.1000	C123	S
5	5	0	class3	Allen, Mr. William Henry	male	35	0	373450	8.0500		S
7	7	0	class1	McCarthy, Mr. Timothy J	male	54	0	17463	51.8625	E46	S

Ajustar o Modelo

m_log0 <- glm(Survived ~ Pclass + Sex,
              data = ds,
              family = binomial(link = "logit"))

m_log1 <- glm(Survived ~ Pclass * Sex,
              data = ds,
              family = binomial(link = "logit"))

Tabela da ANOVA

library(car)

Anova(m_log1, type = 3)

	LR Chisq	Df	Pr(>Chisq)
Pclass	92.34621	2	0e+00
Sex	223.13415	1	0e+00
Pclass:Sex	30.15567	2	3e-07

# Linha da interacção na tabela da ANOVA
anova(m_log0, m_log1, test = "Chisq")

Resid. Df	Resid. Dev	Df	Deviance	Pr(>Chi)
710	672.4307	NA	NA	NA
708	642.2750	2	30.15567	3e-07

O que é que estás a testar aqui?

Pclass: diferenças entre classes em média.
Sex: diferenças entre sexos em média.
Pclass:Sex: se a diferença entre sexos muda com a classe (interacção).

Interpretar Coeficientes

Os coeficientes do modelo estão em log-odds.

coef(m_log1) → log-odds.
exp(coef(m_log1)) → odds ratios.

coef(m_log1)

 (Intercept)      Pclass1      Pclass2         Sex1 Pclass1:Sex1 Pclass2:Sex1 
  0.28347600   1.15957992   0.06901481   1.57608030   0.28897003   0.49917712

exp(coef(m_log1))

 (Intercept)      Pclass1      Pclass2         Sex1 Pclass1:Sex1 Pclass2:Sex1 
    1.327737     3.188594     1.071452     4.835963     1.335052     1.647365

Interpretação rápida

Se um coeficiente for 0.7, então exp(0.7) ≈ 2.01: as odds são ~2x maiores (não é o mesmo que “probabilidade 2x maior”).

`emmeans`: probabilidades e comparações

O grande truque é pedir type = "response", para obteres probabilidades (em vez de log-odds).

emm_cells <- emmeans(m_log1, ~ Pclass * Sex, type = "response")

emm_pclass <- emmeans(m_log1, ~ Pclass, type = "response")
emm_sex <- emmeans(m_log1, ~ Sex, type = "response")

emm_cells

 Pclass Sex     prob     SE  df asymp.LCL asymp.UCL
 class1 female 0.965 0.0200 Inf    0.8963     0.989
 class2 female 0.919 0.0317 Inf    0.8310     0.963
 class3 female 0.461 0.0494 Inf    0.3667     0.558
 class1 male   0.396 0.0487 Inf    0.3056     0.494
 class2 male   0.152 0.0360 Inf    0.0935     0.236
 class3 male   0.150 0.0225 Inf    0.1113     0.200

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale

emm_pclass

 Pclass  prob     SE  df asymp.LCL asymp.UCL
 class1 0.809 0.0480 Inf     0.697     0.886
 class2 0.587 0.0618 Inf     0.463     0.701
 class3 0.280 0.0267 Inf     0.231     0.335

Results are averaged over the levels of: Sex 
Confidence level used: 0.95 
Intervals are back-transformed from the logit scale

emm_sex

 Sex     prob     SE  df asymp.LCL asymp.UCL
 female 0.865 0.0292 Inf     0.797     0.913
 male   0.215 0.0219 Inf     0.176     0.261

Results are averaged over the levels of: Pclass 
Confidence level used: 0.95 
Intervals are back-transformed from the logit scale

Pairwise comparisons (ex.: classes, corrigindo para múltiplas comparações conforme o default do emmeans):

pairs(emm_pclass)

 contrast        odds.ratio   SE  df null z.ratio p.value
 class1 / class2       2.98 1.20 Inf    1   2.713  0.0183
 class1 / class3      10.89 3.68 Inf    1   7.066 <0.0001
 class2 / class3       3.66 1.05 Inf    1   4.515 <0.0001

Results are averaged over the levels of: Sex 
P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale

pairs(emmeans(m_log1, ~ Pclass | Sex, type = "response"))

Sex = female:
 contrast        odds.ratio     SE  df null z.ratio p.value
 class1 / class2       2.41  1.750 Inf    1   1.213  0.4453
 class1 / class3      31.99 19.800 Inf    1   5.588 <0.0001
 class2 / class3      13.26  6.230 Inf    1   5.501 <0.0001

Sex = male:
 contrast        odds.ratio     SE  df null z.ratio p.value
 class1 / class2       3.67  1.270 Inf    1   3.756  0.0005
 class1 / class3       3.71  0.998 Inf    1   4.874 <0.0001
 class2 / class3       1.01  0.334 Inf    1   0.031  0.9995

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale

pairs(emmeans(m_log1, ~ Sex | Pclass, type = "response"))

Pclass = class1:
 contrast      odds.ratio    SE  df null z.ratio p.value
 female / male      41.68 25.90 Inf    1   6.000 <0.0001

Pclass = class2:
 contrast      odds.ratio    SE  df null z.ratio p.value
 female / male      63.47 32.40 Inf    1   8.141 <0.0001

Pclass = class3:
 contrast      odds.ratio    SE  df null z.ratio p.value
 female / male       4.83  1.28 Inf    1   5.938 <0.0001

Tests are performed on the log odds ratio scale

Um gráfico de interacção (sem drama)

emm_df <- as.data.frame(emm_cells)

# Nota: em GLMs binomiais, a coluna do emmeans com probabilidades chama-se tipicamente `prob`
# e os intervalos assimptóticos `asymp.LCL` / `asymp.UCL`.

ggplot(emm_df, aes(x = Pclass, y = prob, color = Sex, group = Sex)) +
geom_point(size = 2) +
geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.15) +
geom_line() +
theme_classic() +
labs(x = "Classe", y = "P(Sobreviver)")

Diagnósticos

Em logística, o que costuma matar análises não é “não-normalidade dos resíduos” (isso já estava assumido), mas sim:

especificação errada (falta de preditores importantes / não-linearidades),
separação (quase) completa,
overdispersion em dados agregados,
pontos influentes.

performance::check_model(m_log1)

4 Binomial com Proporções

Se não tens uma VD 0/1 por linha, mas sim contagens de sucessos e falhas por condição, ajusta um binomial assim:

# Exemplo didáctico (simulado) com 3 grupos
dat_prop <- data.frame(grupo = factor(c("A", "B", "C")),
                       sucessos = c(12, 20, 28),
                       total = c(40, 40, 40))

dat_prop$falhas <- dat_prop$total - dat_prop$sucessos

dat_prop

grupo	sucessos	total	falhas
A	12	40	28
B	20	40	20
C	28	40	12

m_prop <- glm(cbind(sucessos, falhas) ~ grupo,
             data = dat_prop,
             family = binomial(link = "logit"))

Anova(m_prop, type = 3)

	LR Chisq	Df	Pr(>Chisq)
grupo	13.16526	2	0.0013842

emmeans(m_prop, ~ grupo, type = "response")

 grupo prob     SE  df asymp.LCL asymp.UCL
 A      0.3 0.0725 Inf     0.179     0.457
 B      0.5 0.0791 Inf     0.350     0.650
 C      0.7 0.0725 Inf     0.543     0.821

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale

pairs(emmeans(m_prop, ~ grupo, type = "response"))

 contrast odds.ratio     SE  df null z.ratio p.value
 A / B         0.429 0.2010 Inf    1  -1.810  0.1661
 A / C         0.184 0.0896 Inf    1  -3.473  0.0015
 B / C         0.429 0.2010 Inf    1  -1.810  0.1661

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale

Porquê cbind(sucessos, falhas)?

Porque assim o modelo sabe quantas tentativas existem em cada linha, e não trata “0.70” e “7/10” como a mesma coisa.

5 Poisson: modelos para contagens

Quando a VD é uma contagem (número de eventos), um ponto de partida natural é Poisson com link log.

Vamos usar um exemplo real do curso: mortes de manatins (ManateeDeaths) em função do número de barcos (Powerboats).

ds <- read.csv("../data/manatees.csv")
head(ds)

Year	ManateeDeaths	Powerboats
1982	13	447
1983	21	460
1984	24	481
1985	16	498
1986	24	513
1987	20	512

Ajustar Poisson

m_pois0 <- glm(ManateeDeaths ~ 1,
              data = ds,
              family = poisson(link = "log"))

m_pois1 <- glm(ManateeDeaths ~ Powerboats,
              data = ds,
              family = poisson(link = "log"))

anova(m_pois0, m_pois1, test = "Chisq")

Resid. Df	Resid. Dev	Df	Deviance	Pr(>Chi)
34	543.05050	NA	NA	NA
33	66.97316	1	476.0773	0

summary(m_pois1)


Call:
glm(formula = ManateeDeaths ~ Powerboats, family = poisson(link = "log"), 
    data = ds)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 1.6319121  0.1278815   12.76   <2e-16 ***
Powerboats  0.0030391  0.0001503   20.22   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 543.051  on 34  degrees of freedom
Residual deviance:  66.973  on 33  degrees of freedom
AIC: 272

Number of Fisher Scoring iterations: 4

Interpretar: multiplicadores (IRR)

Em Poisson, exp(beta) é um multiplicador na taxa esperada.

exp(coef(m_pois1))

(Intercept)  Powerboats 
   5.113643    1.003044

Se exp(beta1) = 1.10, isso sugere ~10% mais eventos por cada unidade de X (dependendo da escala de X).

Diagnóstico-chave: overdispersion

Poisson assume (muito aproximadamente) que variância ≈ média. Se a variância for muito maior, os erros-padrão ficam optimistas e os p-values ficam “bons demais”.

performance::check_overdispersion(m_pois1)

# Overdispersion test

       dispersion ratio =  1.983
  Pearson's Chi-Squared = 65.429
                p-value =  0.001

Se houver overdispersion: 3 caminhos comuns

Quasi-Poisson (ajusta os SE, mas não dá AIC “a sério”):

m_qpois <- glm(ManateeDeaths ~ Powerboats,
              data = ds,
              family = quasipoisson(link = "log"))
summary(m_qpois)


Call:
glm(formula = ManateeDeaths ~ Powerboats, family = quasipoisson(link = "log"), 
    data = ds)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.6319121  0.1800673   9.063 1.79e-10 ***
Powerboats  0.0030391  0.0002116  14.361 9.44e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 1.982688)

    Null deviance: 543.051  on 34  degrees of freedom
Residual deviance:  66.973  on 33  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

Negative Binomial (muito comum em contagens):

library(MASS)

m_nb <- glm.nb(ManateeDeaths ~ Powerboats, data = ds)
summary(m_nb)


Call:
glm.nb(formula = ManateeDeaths ~ Powerboats, data = ds, init.theta = 63.1074526, 
    link = log)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 1.5519755  0.1624194   9.555   <2e-16 ***
Powerboats  0.0031387  0.0001973  15.906   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(63.1075) family taken to be 1)

    Null deviance: 298.374  on 34  degrees of freedom
Residual deviance:  32.973  on 33  degrees of freedom
AIC: 261.84

Number of Fisher Scoring iterations: 1

              Theta:  63.1 
          Std. Err.:  30.2 

 2 x log-likelihood:  -255.836

Rever o modelo (muitas vezes, overdispersion é sinal de variável omitida ou estrutura de agrupamento → GLMM).

6 Offset: Taxas em Vez de Contagens

Às vezes a pergunta não é “quantos eventos”, mas sim “qual é a taxa por exposição” (tempo, população, km percorridos, etc.).

Um offset(log(exposicao)) diz ao modelo: “a contagem esperada deve crescer proporcionalmente à exposição”.

# Exemplo didáctico (simulado)
set.seed(2)

dat_rate <- data.frame(exposicao = c(10, 10, 20, 20), grupo = factor(c("A", "B", "A", "B")))

# taxa verdadeira: grupo B tem taxa 2x
lambda <- ifelse(dat_rate$grupo == "B", 0.6, 0.3) * dat_rate$exposicao

dat_rate$y <- rpois(nrow(dat_rate), lambda)

dat_rate

exposicao	grupo	y
10	A	1
10	B	7
20	A	6
20	B	8

m_rate <- glm(y ~ grupo + offset(log(exposicao)),
              data = dat_rate,
              family = poisson(link = "log"))

summary(m_rate)


Call:
glm(formula = y ~ grupo + offset(log(exposicao)), family = poisson(link = "log"), 
    data = dat_rate)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.0742     0.2289  -4.694 2.68e-06 ***
grupo1       -0.3811     0.2289  -1.665   0.0959 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 5.4383  on 3  degrees of freedom
Residual deviance: 2.4615  on 2  degrees of freedom
AIC: 19.865

Number of Fisher Scoring iterations: 4

emmeans(m_rate, ~ grupo, type = "response")

 grupo rate   SE  df asymp.LCL asymp.UCL
 A     3.30 1.25 Inf      1.57      6.92
 B     7.07 1.83 Inf      4.26     11.73

Results are averaged over the levels of: exposicao 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale

7 Mini-checklist antes de escrever resultados

A VD é mesmo binária/proporção/contagem? (Se for 0–100 contínuo, talvez não.)
As variáveis categóricas estão como factor()?
Fizeste pelo menos uma comparação de modelos (compacto vs aumentado)?
Usaste emmeans(..., type = "response") para comunicar probabilidades/taxas?
Verificaste overdispersion em Poisson?
Fizeste um check_model() para apanhar problemas óbvios?

8 Exercícios

Titanic: acrescentar idade

No Titanic, refaz o modelo com Age e compara:

Survived ~ Pclass * Sex
Survived ~ Pclass * Sex + Age

Dica: remove NAs de Age.

Manatees: comparar Poisson vs quasi-Poisson

Ajusta:

glm(..., family = poisson)
glm(..., family = quasipoisson)

Compara as conclusões (coeficientes semelhantes, SE/p-values diferentes).

`emmeans` em logística

No modelo Titanic, pede:

emmeans(m_log1, list(pairwise ~ Pclass | Sex), type = "response")

E resume em 2 frases o padrão de diferenças.

0 Setup

1 GLM, Quando, Como e Porquê?

2 Comparação de GLMs

3 Binomial: Regressão Logística

Ajustar o Modelo

Tabela da ANOVA

Interpretar Coeficientes

emmeans: probabilidades e comparações

Um gráfico de interacção (sem drama)

Diagnósticos

4 Binomial com Proporções

5 Poisson: modelos para contagens

Ajustar Poisson

Interpretar: multiplicadores (IRR)

Diagnóstico-chave: overdispersion

Se houver overdispersion: 3 caminhos comuns

6 Offset: Taxas em Vez de Contagens

7 Mini-checklist antes de escrever resultados

8 Exercícios

Titanic: acrescentar idade

Manatees: comparar Poisson vs quasi-Poisson

emmeans em logística

`emmeans`: probabilidades e comparações

`emmeans` em logística