GLMs no R: Proporções e Contagens na Prática

Modelos Lineares

GLM

Logística

Poisson

emmeans

Diagnósticos

Fluxos práticos no R para respostas binárias, proporções, contagens e taxas.

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 Preparação

Este é o tutorial prático. A explicação de família, função de ligação, logit, deviance, sobredispersão e offsets está em Modelos Lineares Generalizados. Leia esse capítulo quando precisar de perceber por que o modelo funciona. Aqui vamos praticar como passar dos dados para um GLM e voltar à escala da resposta.

Os exemplos seguem sempre a mesma sequência:

olhar para a resposta e para os preditores;
escolher a família a partir do tipo de resposta;
ajustar o modelo que representa a pergunta;
obter uma tabela de testes dos efeitos;
voltar à escala da resposta com previsões e comparações;
desenhar o padrão e verificar problemas específicos do GLM.

Se fórmulas, factores, contrastes ou interacções ainda forem pouco familiares, reveja primeiro modelos lineares.

# Run once in the console, not in an analysis script.
install.packages(c("ggplot2", "car", "emmeans", "performance"))

library(ggplot2)
library(car)
library(emmeans)
library(performance)

1 Escolher a Família

Resposta	Primeiro modelo	Exemplo
0/1	`glm(..., family = binomial)`	sobrevivência no Titanic
sucessos em várias tentativas	`glm(cbind(sucessos, falhas) ~ ..., family = binomial)`	proporções por grupo
contagem 0, 1, 2, …	`glm(..., family = poisson)`	mortes de manatins
contagem por exposição	Poisson + `offset(log(exposicao))`	taxa por unidade de exposição

Esta tabela é uma cábula para começar. A teoria de GLMs explica quando a relação média–variância, a função de ligação ou a dependência exigem outra escolha.

2 Binomial: Regressão Logística

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

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

# Confirm the variables entering the model before recoding them.
stopifnot(all(c("Survived", "Pclass", "Sex") %in% names(ds)))

# Treat class as a nominal factor; this example does not use age or fare.
ds$Pclass <- factor(paste0("class", ds$Pclass))
ds$Sex <- factor(ds$Sex)
ds <- ds[complete.cases(ds[, c("Survived", "Pclass", "Sex")]), , drop = FALSE]

str(ds[, c("Survived", "Pclass", "Sex")])

'data.frame':   714 obs. of  3 variables:
 $ Survived: int  0 1 1 1 0 0 0 1 1 1 ...
 $ Pclass  : Factor w/ 3 levels "class1","class2",..: 3 1 3 1 3 1 3 3 2 3 ...
 $ Sex     : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 1 1 1 ...

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

Ver os Dados

Como Survived está codificado 0/1, a média em cada célula é a proporção observada de sobreviventes.

sobrevivencia_observada <- aggregate(
  Survived ~ Pclass + Sex,
  data = ds,
  FUN = mean
)

sobrevivencia_observada

Pclass	Sex	Survived
class1	female	0.9647059
class2	female	0.9189189
class3	female	0.4607843
class1	male	0.3960396
class2	male	0.1515152
class3	male	0.1501976

ggplot(
  sobrevivencia_observada,
  aes(x = Pclass, y = Survived, shape = Sex, group = Sex)
) +
  geom_point(size = 2.5) +
  geom_line() +
  theme_classic() +
  labs(
    x = "Classe",
    y = "Proporção observada de sobreviventes",
    shape = "Sexo"
  )

Antes de olhar para valores-p, já sabemos que o modelo tem de representar diferenças entre classes, entre sexos e a possibilidade de a diferença entre sexos mudar de classe para classe.

Ajustar o Modelo

Ajustamos directamente o modelo factorial que corresponde à pergunta. Como vamos pedir testes Tipo III, definimos contrastes de soma nos factores.

m_log <- glm(
  Survived ~ Pclass * Sex,
  data = ds,
  family = binomial(link = "logit"),
  contrasts = list(
    Pclass = "contr.sum",
    Sex = "contr.sum"
  )
)

Testar Todos os Efeitos

Para um glm, car::Anova() pode produzir testes qui-quadrado de razão de verosimilhanças. Esta é a tabela que queremos no fluxo principal:

car::Anova(m_log, 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

A tabela testa Pclass, Sex e Pclass:Sex no mesmo modelo. Não precisamos de ajustar primeiro um modelo só com efeitos principais apenas para obter o teste da interacção.

Os contrastes importam para testes Tipo III. Aqui ficam especificados no próprio modelo, sem alterar as opções globais da sessão. A razão conceptual para Tipo III e a relação com modelos restritos pertence ao capítulo teórico de comparação de modelos.

Que Perguntas Correspondem aos Termos?

Pclass: diferenças entre classes, resumidas sobre os sexos segundo o modelo e a codificação.
Sex: diferenças entre sexos, resumidas sobre as classes.
Pclass:Sex: se a diferença entre sexos varia com a classe.

Se a interacção for importante, estas médias globais deixam de contar a história toda. Nesse caso interessa sobretudo comparar classes dentro de cada sexo e sexos dentro de cada classe. Para a geometria e interpretação de interacções entre factores, veja interacção entre dois factores.

Reconhecer Coeficientes em Log-Odds

summary(m_log) apresenta coeficientes na escala logit. Para os exponenciar:

summary(m_log)


Call:
glm(formula = Survived ~ Pclass * Sex, family = binomial(link = "logit"), 
    data = ds, contrasts = list(Pclass = "contr.sum", Sex = "contr.sum"))

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.28348    0.14112   2.009   0.0446 *  
Pclass1       1.15958    0.22831   5.079 3.79e-07 ***
Pclass2       0.06901    0.20390   0.338   0.7350    
Sex1          1.57608    0.14112  11.169  < 2e-16 ***
Pclass1:Sex1  0.28897    0.22831   1.266   0.2056    
Pclass2:Sex1  0.49918    0.20390   2.448   0.0144 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 964.52  on 713  degrees of freedom
Residual deviance: 642.28  on 708  degrees of freedom
AIC: 654.28

Number of Fisher Scoring iterations: 5

exp(coef(m_log))

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

A exponencial de um coeficiente é uma razão de odds para a comparação representada por esse coeficiente e pela codificação usada. Para comunicar o padrão factorial, o percurso principal abaixo usa probabilidades previstas e comparações de emmeans.

`emmeans`: Probabilidades e Comparações

Para comunicar o padrão é normalmente mais fácil voltar à escala da probabilidade. emmeans(..., type = "response") apresenta as médias marginais como probabilidades previstas. Os contrastes continuam a ser calculados na escala linear do modelo. Quando também os pedimos com type = "response", numa logística aparecem como odds ratios. Assim evitamos obrigar o leitor a interpretar log-odds quando a probabilidade ou a razão de odds responde melhor à pergunta.

emm_log <- emmeans(
  m_log,
  list(
    pairwise ~ Pclass,
    ~ Sex,
    pairwise ~ Pclass | Sex,
    pairwise ~ Sex | Pclass
  ),
  type = "response"
)

emm_log

$`emmeans of 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 

$`pairwise differences of Pclass`
 1               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 

$`emmeans of 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 

$`emmeans of Pclass | Sex`
Sex = female:
 Pclass  prob     SE  df asymp.LCL asymp.UCL
 class1 0.965 0.0200 Inf    0.8963     0.989
 class2 0.919 0.0317 Inf    0.8310     0.963
 class3 0.461 0.0494 Inf    0.3667     0.558

Sex = male:
 Pclass  prob     SE  df asymp.LCL asymp.UCL
 class1 0.396 0.0487 Inf    0.3056     0.494
 class2 0.152 0.0360 Inf    0.0935     0.236
 class3 0.150 0.0225 Inf    0.1113     0.200

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

$`pairwise differences of Pclass | Sex`
Sex = female:
 2               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:
 2               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 

$`emmeans of Sex | Pclass`
Pclass = class1:
 Sex     prob     SE  df asymp.LCL asymp.UCL
 female 0.965 0.0200 Inf    0.8963     0.989
 male   0.396 0.0487 Inf    0.3056     0.494

Pclass = class2:
 Sex     prob     SE  df asymp.LCL asymp.UCL
 female 0.919 0.0317 Inf    0.8310     0.963
 male   0.152 0.0360 Inf    0.0935     0.236

Pclass = class3:
 Sex     prob     SE  df asymp.LCL asymp.UCL
 female 0.461 0.0494 Inf    0.3667     0.558
 male   0.150 0.0225 Inf    0.1113     0.200

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

$`pairwise differences of Sex | Pclass`
Pclass = class1:
 2             odds.ratio    SE  df null z.ratio p.value
 female / male      41.68 25.90 Inf    1   6.000 <0.0001

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

Pclass = class3:
 2             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

Ao escrever os resultados, diga se está a comunicar probabilidades previstas, diferenças na escala logit ou odds ratios. São maneiras diferentes de olhar para o mesmo modelo, mas não são números intercambiáveis.

Ver a Interacção

emm_cells <- emmeans(m_log, ~ Pclass * Sex, type = "response")
emm_df <- as.data.frame(emm_cells)

# The probability column is usually `prob`; intervals are `asymp.LCL` and `asymp.UCL`.
ggplot(emm_df, aes(x = Pclass, y = prob, colour = 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)")

O gráfico ajuda a ver se as diferenças entre sexos parecem semelhantes entre classes e mostra a incerteza das probabilidades previstas. Para a geometria e interpretação de interacções entre factores, veja interacção entre dois factores. Não transforma, por si só, esta associação observacional numa explicação causal.

Diagnósticos

Numa regressão logística não esperamos resíduos normais. Por isso não vale a pena importar mecanicamente a checklist do lm(). Os problemas relevantes incluem, entre outros:

forma funcional inadequada para preditores quantitativos;
separação completa ou quase completa;
dispersão adicional em dados binomiais agregados;
observações muito influentes;
dependência que o modelo não representou.

performance::check_model(m_log)

Use estes diagnósticos como pistas sobre o modelo, não como uma tabela de «passou/falhou». Veja diagnóstico em GLMs e, para retomar o papel do termo de erro na inferência, a secção teórica sobre inferência. Para dados que partilham unidades, veja dependência entre observações.

3 Binomial com Proporções

Agora cada linha resume várias tentativas. Imagine três grupos com 40 organismos cada e números diferentes de sucessos.

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$proporcao <- dat_prop$sucessos / dat_prop$total

dat_prop

grupo	sucessos	total	falhas	proporcao
A	12	40	28	0.3
B	20	40	20	0.5
C	28	40	12	0.7

Ver as Proporções

ggplot(dat_prop, aes(x = grupo, y = proporcao)) +
  geom_point(size = 3) +
  ylim(0, 1) +
  theme_classic() +
  labs(x = "Grupo", y = "Proporção observada de sucessos")

Ajustar e Testar

cbind(sucessos, falhas) informa o modelo de quantas tentativas estão por trás de cada proporção.

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

car::Anova(m_prop, type = 3)

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

Uma proporção de .70 baseada em 7/10 e outra baseada em 700/1000 têm o mesmo valor observado, mas não a mesma precisão. É por isso que não analisamos apenas a coluna proporcao.

Probabilidades e Comparações

emm_prop <- emmeans(m_prop, ~ grupo, type = "response")
emm_prop

 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

summary(pairs(emm_prop), type = "response")

contrast	odds.ratio	SE	df	null	z.ratio	p.value
A / B	0.4285714	0.2005822	Inf	1	-1.810368	0.1661361
A / C	0.1836735	0.0896235	Inf	1	-3.472888	0.0014938
B / C	0.4285714	0.2005822	Inf	1	-1.810368	0.1661361

emm_prop_df <- as.data.frame(emm_prop)

ggplot(emm_prop_df, aes(x = grupo, y = prob)) +
  geom_point(size = 2.5) +
  geom_errorbar(
    aes(ymin = asymp.LCL, ymax = asymp.UCL),
    width = 0.12
  ) +
  ylim(0, 1) +
  theme_classic() +
  labs(x = "Grupo", y = "Probabilidade prevista")

E os Testes Clássicos de Proporções?

binom.test() e prop.test() continuam a responder a perguntas simples sobre proporções. O GLM torna-se especialmente útil quando queremos acrescentar covariáveis, interacções ou manter a mesma linguagem de modelação usada no resto do curso.

4 Poisson: Modelos para Contagens

Vamos usar as mortes anuais de manatins (ManateeDeaths) e o número de barcos a motor registados (Powerboats). A pergunta é descritiva: a contagem esperada de mortes varia com o número de barcos registados? Estes dados observacionais não identificam, sozinhos, um efeito causal.

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

Ver os Dados

ggplot(ds, aes(x = Powerboats, y = ManateeDeaths)) +
  geom_point(size = 2) +
  theme_classic() +
  labs(
    x = "Barcos a motor registados",
    y = "Mortes anuais de manatins"
  )

Ajustar e Testar

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

summary(m_pois)


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

car::Anova(m_pois, type = 3)

	LR Chisq	Df	Pr(>Chisq)
Powerboats	476.0773	1	0

`emmeans()` Sem Grupos e o Intercepto

Como Powerboats é quantitativo, este exemplo torna mais clara a diferença entre uma estimativa sem grupos e o intercepto.

Primeiro, uma única estimativa na grelha de referência:

emmeans(m_pois, ~ 1, type = "response")

 1       rate   SE  df asymp.LCL asymp.UCL
 overall 50.9 1.29 Inf      48.4      53.5

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

Atenção ao ~ 1

~ 1 pede uma única média marginal estimada. Não é uma sintaxe especial para «dar o intercepto». Com covariáveis quantitativas, emmeans usa por defeito os valores da grelha de referência, normalmente as médias das covariáveis. Para obter a previsão correspondente ao intercepto, fixe explicitamente essas covariáveis em zero — ou recentre-as primeiro num zero cientificamente útil.

Por defeito, a covariável quantitativa é avaliada no seu valor médio. O resultado é, portanto, a contagem esperada de mortes no número médio de barcos observado, não a exponencial do intercepto.

Para pedir exactamente a previsão correspondente ao intercepto, fixamos Powerboats = 0:

exp(coef(m_pois)["(Intercept)"])

(Intercept) 
   5.113643

emmeans(
  m_pois,
  ~ 1,
  at = list(Powerboats = 0),
  type = "response"
)

 1       rate    SE  df asymp.LCL asymp.UCL
 overall 5.11 0.654 Inf      3.98      6.57

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

As duas estimativas coincidem. Na escala logarítmica do modelo Poisson, o intercepto é a previsão quando Powerboats = 0. Com type = "response", emmeans aplica a exponencial e devolve a contagem esperada nessa condição.

Neste conjunto de dados, porém, zero barcos está muito fora do intervalo observado, que começa perto de 450. A igualdade matemática continua correcta, mas a previsão não é uma boa quantidade para comunicar. Recentremos, por exemplo, em 700 barcos:

ds$Powerboats_700 <- ds$Powerboats - 700

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

exp(coef(m_pois_700)["(Intercept)"])

(Intercept) 
    42.9189

emmeans(
  m_pois_700,
  ~ 1,
  at = list(Powerboats_700 = 0),
  type = "response"
)

 1       rate  SE  df asymp.LCL asymp.UCL
 overall 42.9 1.3 Inf      40.4      45.6

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

Agora zero significa 700 barcos e a estimativa pode ser lida directamente como a contagem esperada de mortes nesse ponto. O princípio é o mesmo na regressão logística: fixar todas as covariáveis quantitativas em zero dá o intercepto na escala logit. type = "response" transforma-o numa probabilidade prevista.

A tabela de car::Anova() dá o qui-quadrado do efeito de Powerboats. Neste modelo simples há apenas um preditor, mas o mesmo comando estende-se a modelos com vários termos e interacções.

A exponencial do declive põe a associação numa escala multiplicativa:

exp(coef(m_pois))

(Intercept)  Powerboats 
   5.113643    1.003044

Como este modelo não tem offset, a exponencial do declive é uma razão de contagens esperadas por unidade adicional de Powerboats. Ainda não é uma razão de taxas.

Desenhar a Previsão

grelha_pois <- data.frame(
  Powerboats = seq(
    min(ds$Powerboats),
    max(ds$Powerboats),
    length.out = 100
  )
)

grelha_pois$previsto <- predict(
  m_pois,
  newdata = grelha_pois,
  type = "response"
)

ggplot(ds, aes(x = Powerboats, y = ManateeDeaths)) +
  geom_point(alpha = 0.7) +
  geom_line(
    data = grelha_pois,
    aes(y = previsto),
    linewidth = 0.9
  ) +
  theme_classic() +
  labs(
    x = "Barcos a motor registados",
    y = "Mortes anuais de manatins"
  )

Verificar Sobredispersão

performance::check_overdispersion(m_pois)

# Overdispersion test

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

E olhe para o painel de diagnóstico numa figura suficientemente grande:

performance::check_model(m_pois)

Se houver sobredispersão, não responda apenas «trocar de distribuição». Pergunte se falta estrutura ao modelo: um preditor importante, dependência entre observações, heterogeneidade entre unidades ou uma relação média–variância que Poisson não representa bem. As alternativas estão explicadas em dispersão nos GLMs.

5 Offset: Taxas em Vez de Contagens

Às vezes a pergunta não é «quantos eventos ocorreram?», mas «qual é a taxa por unidade de exposição?». A exposição pode ser tempo, população, quilómetros percorridos ou outra oportunidade de observar o evento.

Antes de usar um offset, confirme que a exposição é positiva e que a unidade tem significado substantivo. A derivação teórica do offset mostra por que log(exposicao) entra com coeficiente fixo em 1. offset(log(exposicao)) diz ao modelo que a contagem esperada deve crescer proporcionalmente à exposição. O offset fixa o coeficiente da exposição em 1, o que equivale a modelar a taxa em vez da contagem bruta.

# Simulate a small teaching example reproducibly.
set.seed(2)

dat_rate <- data.frame(exposicao = c(10, 10, 20, 20), grupo = factor(c("A", "B", "A", "B")))
stopifnot(all(dat_rate$exposicao > 0))

# Set group B's rate to twice group A's rate.
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

dat_rate$taxa_observada <- dat_rate$y / dat_rate$exposicao
dat_rate

exposicao	grupo	y	taxa_observada
10	A	1	0.1
10	B	7	0.7
20	A	6	0.3
20	B	8	0.4

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

summary(m_rate)


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

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

car::Anova(m_rate, type = 3)

	LR Chisq	Df	Pr(>Chisq)
grupo	2.976856	1	0.0844632

# Set the reference exposure explicitly to one unit.
emm_rate_unidade <- emmeans(
  m_rate, ~ grupo,
  offset = log(1),
  type = "response"
)
emm_rate_unidade

 grupo  rate     SE  df asymp.LCL asymp.UCL
 A     0.233 0.0882 Inf     0.111     0.489
 B     0.500 0.1290 Inf     0.301     0.829

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

summary(pairs(emm_rate_unidade), type = "response")

contrast	ratio	SE	df	null	z.ratio	p.value
A / B	0.4666667	0.2136105	Inf	1	-1.665018	0.0959091

offset = log(1) fixa a exposição de referência em uma unidade. As previsões podem então ser lidas como taxas esperadas por unidade de exposição. Uma comparação por pares na escala da resposta é uma razão de taxas. Se a exposição for pessoa-tempo e o evento for uma incidência, essa razão pode ser descrita como incidence-rate ratio (IRR). Diga sempre qual é a unidade: pessoa-ano, quilómetro, dia de amostragem, etc.

6 Antes de Escrever os Resultados

Antes de fechar a análise, confirme:

se a resposta é realmente binária, uma proporção ou uma contagem;
se as variáveis categóricas foram representadas como factor() quando isso corresponde à pergunta;
que termos foram testados na tabela Tipo III e qual é a pergunta de cada um;
em que escala vai comunicar as estimativas: probabilidade, odds ratio, contagem esperada ou taxa;
no caso duma taxa, qual é a exposição e a sua unidade;
se há sobredispersão, separação, dependência ou outro padrão que torne o modelo demasiado simples. Quando a mesma unidade contribui com várias observações, a dependência deve ser modelada. Veja medidas repetidas e modelos mistos.

A checklist serve para recordar perguntas. Não substitui olhar para os dados, o delineamento e as estimativas.

7 Exercícios

Titanic: Acrescentar Idade

Remova os valores em falta de Age e ajuste:

glm(
  Survived ~ Pclass * Sex + Age,
  data = ...,
  family = binomial,
  contrasts = list(Pclass = "contr.sum", Sex = "contr.sum")
)

Depois use car::Anova(..., type = 3). Que efeitos têm evidência no modelo? Faça previsões com emmeans e explique por que o teste de Age e as comparações entre classes/sexos respondem a perguntas diferentes.

`emmeans` em Logística

No modelo do Titanic, peça:

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

Identifique quais resultados são probabilidades previstas e quais são odds ratios. Distinga efeitos globais da tabela Tipo III de comparações condicionais entre níveis.

Proporções

Altere sucessos para c(12, 20, 22). Volte a ajustar m_prop, peça car::Anova(m_prop, type = 3) e compare o gráfico das probabilidades previstas. O que mudou na magnitude das diferenças e na evidência global para grupo?

Manatins

Ajuste o Poisson, obtenha car::Anova(m_pois, type = 3), desenhe as contagens previstas e corra check_overdispersion(m_pois). Escreva quatro linhas que separem: associação estimada, teste do efeito, adequação da dispersão e limites causais.

Solução: Titanic com Idade

Removemos em falta nas quatro variáveis usadas e ajustamos a idade como covariável quantitativa:

ds_titanic <- read.csv("../data/titanic.csv")
ds_titanic$Pclass <- factor(paste0("class", ds_titanic$Pclass))
ds_titanic$Sex <- factor(ds_titanic$Sex)
ds_age <- ds_titanic[complete.cases(
  ds_titanic[, c("Survived", "Pclass", "Sex", "Age")]
), ]

m_log_age <- glm(
  Survived ~ Pclass * Sex + Age,
  data = ds_age,
  family = binomial(link = "logit"),
  contrasts = list(Pclass = "contr.sum", Sex = "contr.sum")
)

car::Anova(m_log_age, type = 3)

emmeans(m_log_age, ~ Pclass * Sex, type = "response")

A tabela Tipo III testa os termos Pclass, Sex, Pclass:Sex e Age no mesmo modelo. Age testa uma associação com a sobrevivência condicionada pelos outros termos; as comparações entre classes ou sexos descrevem diferenças de probabilidade prevista condicionadas na idade de referência da grelha.

Solução: emmeans no Modelo Logístico

emm_log <- emmeans(
  m_log,
  list(
    pairwise ~ Pclass,
    ~ Sex,
    pairwise ~ Pclass | Sex,
    pairwise ~ Sex | Pclass
  ),
  type = "response"
)

emm_log

As linhas de médias estimadas são probabilidades previstas. As linhas de contrast, pedidas com type = "response", são razões de odds. Os resultados sem | resumem comparações sobre o outro factor; os que têm | são comparações condicionais dentro de cada sexo ou classe. A tabela Tipo III testa os efeitos globais do modelo, não substitui estas comparações específicas.

Solução: Proporções com o Grupo C Alterado

dat_prop$sucessos <- c(12, 20, 22)
dat_prop$falhas <- dat_prop$total - dat_prop$sucessos
dat_prop$proporcao <- dat_prop$sucessos / dat_prop$total

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

car::Anova(m_prop, type = 3)

emm_prop <- emmeans(m_prop, ~ grupo, type = "response")
emm_prop

emm_prop_df <- as.data.frame(emm_prop)
ggplot(emm_prop_df, aes(grupo, prob)) +
  geom_point(size = 2.5) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.12) +
  coord_cartesian(ylim = c(0, 1)) +
  theme_classic() +
  labs(x = "Grupo", y = "Probabilidade prevista")

A proporção de C fica mais próxima da de B do que no exemplo original. Compare as probabilidades previstas e a tabela Tipo III para descrever a mudança na magnitude das diferenças e na evidência para grupo, sem tratar o limiar do valor-p como a única conclusão.

Solução: Poisson para as Mortes de Manatins

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

summary(m_pois)
car::Anova(m_pois, type = 3)
performance::check_overdispersion(m_pois)

grelha_pois <- data.frame(
  Powerboats = seq(min(ds$Powerboats), max(ds$Powerboats), length.out = 100)
)
grelha_pois$previsto <- predict(m_pois, newdata = grelha_pois, type = "response")

ggplot(ds, aes(Powerboats, ManateeDeaths)) +
  geom_point(alpha = 0.7) +
  geom_line(data = grelha_pois, aes(y = previsto), linewidth = 0.9) +
  theme_classic() +
  labs(x = "Barcos a motor registados", y = "Mortes anuais de manatins")

Um relato curto deve (1) indicar a direcção e a escala multiplicativa da associação estimada, usando o coeficiente ou exp(coef(m_pois)); (2) distinguir o teste do efeito na tabela Tipo III; (3) dizer o que check_overdispersion() sugere sobre a adequação da variância Poisson; e (4) lembrar que estes dados observacionais não identificam por si só um efeito causal dos barcos nas mortes.

0 Preparação

1 Escolher a Família

2 Binomial: Regressão Logística

Ver os Dados

Ajustar o Modelo

Testar Todos os Efeitos

emmeans: Probabilidades e Comparações

Ver a Interacção

Diagnósticos

3 Binomial com Proporções

Ver as Proporções

Ajustar e Testar

Probabilidades e Comparações

4 Poisson: Modelos para Contagens

Ver os Dados

Ajustar e Testar

emmeans() Sem Grupos e o Intercepto

Desenhar a Previsão

Verificar Sobredispersão

5 Offset: Taxas em Vez de Contagens

6 Antes de Escrever os Resultados

7 Exercícios

Titanic: Acrescentar Idade

emmeans em Logística

Proporções

Manatins

`emmeans`: Probabilidades e Comparações

`emmeans()` Sem Grupos e o Intercepto

`emmeans` em Logística