# Contagens trimestrais — Tabela 4.5 do livro
y <- c(
1, 6, 16, 23, # 1984: T1, T2, T3, T4
27, 39, 31, 30, # 1985
43, 51, 63, 70, # 1986
88, 97, 91, 104, # 1987
110, 113, 149, 159 # 1988
)
i_seq <- 1:20
x <- log(i_seq) # covariavel: log do periodo
ano <- rep(1984:1988, each = 4)
trim <- rep(1:4, times = 5)
dados <- data.frame(periodo = i_seq, ano = ano,
trimestre = trim, y = y, x = round(x, 4))Modelos Lineares Generalizados
Apostila — Regressão de Poisson com Ligação Log: Casos de AIDS na Austrália
Sobre anexo. Este documento implementa passo a passo a regressão de Poisson com ligação log aplicada aos dados de casos de AIDS na Austrália (1984–1988), conforme o Exercício 4.1 de Dobson & Barnett (2018). O objetivo é apresentar de forma didática o contexto do problema, a justificativa do modelo, o ajuste via
glm(), a interpretação completa da saída e a avaliação da qualidade do ajuste.
Este material foi produzido em Quarto Markdown, com apoio de ferramentas de inteligência artificial na organização e síntese do conteúdo, tendo sido integralmente revisado e validado pelo autor.
1 Contexto do Problema
1.1 Os Dados
Os dados (Tabela 4.5 do livro) registram o número de casos de AIDS diagnosticados na Austrália em cada trimestre entre 1984 e 1988, na fase inicial da epidemia, quando os casos pareciam crescer de forma exponencial.
Variáveis:
- \(y_i\): número de casos de AIDS no trimestre \(i\) — variável resposta (contagem)
- \(i\): índice do trimestre (\(i = 1, \ldots, 20\)) — variável explicativa
- \(x_i = \log(i)\): transformação logarítmica do período — covariável usada no modelo
Estamos interessados em verificar se o número de casos de AIDS aumentou ao longo do tempo. Como os dados foram registrados em períodos consecutivos, usamos o período i=1,2,…,20 como variável explicativa para representar a passagem do tempo.
O período é usado como variável explicativa porque ele representa o tempo. Já o logaritmo do período é usado porque o crescimento dos casos parece não ser linear na escala original, mas pode ser aproximadamente linear na escala logarítmica.
| Ano | T1 | T2 | T3 | T4 | Total anual | |
|---|---|---|---|---|---|---|
| 1984 | 1984 | 1 | 6 | 16 | 23 | 46 |
| 1985 | 1985 | 27 | 39 | 31 | 30 | 127 |
| 1986 | 1986 | 43 | 51 | 63 | 70 | 227 |
| 1987 | 1987 | 88 | 97 | 91 | 104 | 380 |
| 1988 | 1988 | 110 | 113 | 149 | 159 | 531 |
| Estatistica | Valor |
|---|---|
| N (trimestres) | 20.0 |
| Total de casos | 1311.0 |
| Minimo | 1.0 |
| Maximo | 159.0 |
| Media | 65.6 |
2 Por que Usar a Distribuição de Poisson?
A distribuição de Poisson é adequada quando a variável resposta satisfaz três condições:
- 1. É uma contagem de eventos (inteiro \(\geq 0\))
- \(y_i\) = número de casos de AIDS no trimestre \(i\) — por definição um inteiro não-negativo.
- 2. Os eventos ocorrem de forma independente
- Cada diagnóstico não influencia diretamente a ocorrência de outro (na fase inicial da epidemia, sem grandes efeitos de rede).
- 3. A média e a variância são da mesma ordem de magnitude
- A propriedade fundamental da Poisson é \(E(Y_i) = \text{var}(Y_i) = \lambda_i\). Verificamos empiricamente:
| Ano | Media | Variancia | Variancia / Media | |
|---|---|---|---|---|
| 1984 | 1984 | 11.50 | 97.67 | 8.49 |
| 1985 | 1985 | 31.75 | 26.25 | 0.83 |
| 1986 | 1986 | 56.75 | 145.58 | 2.57 |
| 1987 | 1987 | 95.00 | 50.00 | 0.53 |
| 1988 | 1988 | 132.75 | 620.25 | 4.67 |
| Razao proxima de 1 sustenta o uso da Poisson. Valores maiores indicam possivel superdispersao. |
Observacao: As razoes variancia/media variam bastante entre anos (grupos de apenas 4 observacoes). Isso e esperado — com grupos tao pequenos, a estimativa de variancia e instavel. A verificacao formal de superdispersao sera feita apos o ajuste do modelo (Secao 7).
3 Por que Usar a Ligação Log?
O modelo especifica \(g(\lambda_i) = \log(\lambda_i)\). Ha três motivos para esta escolha:
- Motivo 1 — Positividade garantida
- O numero de casos é sempre positivo: \(\lambda_i > 0\). Com a ligação log, \(\lambda_i = e^{\eta_i} > 0\) qualquer que seja o valor do preditor linear \(\eta_i\). A ligação identidade (\(\lambda_i = b_1 + b_2 x_i\)) poderia gerar valores negativos.
- Motivo 2 — Ligacao canonica da Poisson
- Do Capítulo 3, o parâmetro natural da Poisson e \(\log(\lambda)\). Usar a ligação canônica simplifica o algoritmo IRLS e garante as melhores propriedades assintóticas dos estimadores.
- Motivo 3 — Interpretabilidade multiplicativa
- Com a ligação log, os coeficientes tem interpretação multiplicativa — \(e^{b_j}\) e o fator pelo qual a contagem média e multiplicada por unidade de aumento em \(x_j\). Isso é natural para dados com crescimento exponencial.
O modelo proposto pelo livro usa \(x_i = \log(i)\):
\[ \log(\lambda_i) = \beta_1 + \beta_2 \log(i) \quad \Longleftrightarrow \quad \lambda_i = e^{\beta_1} \cdot i^{\beta_2} \]
Este é um modelo de potência: os casos crescem proporcionalmente a \(i\) elevado ao expoente \(\beta_2\).
4 Analise Exploratoria
4.1 Graficos de dispersão
g1 <- ggplot(dados, aes(x = periodo, y = y)) +
geom_line(color = "blue4", linewidth = 0.9) +
geom_point(color = "blue4", size = 3) +
scale_x_continuous(breaks = c(1,5,9,13,17),
labels = c("1984-T1","1985-T1","1986-T1","1987-T1","1988-T1")) +
labs(title = "Serie temporal — escala original",
subtitle = "Crescimento claramente nao-linear",
x = "Periodo", y = "Numero de casos (y)") +
tema
g2 <- ggplot(dados, aes(x = x, y = log(y))) +
geom_point(color = "red4", size = 3) +
geom_smooth(method = "lm", se = TRUE,
color = "red4", fill = "red4", alpha = 0.15) +
labs(title = "Escala log-log — verificacao do modelo",
subtitle = "log(y) vs. log(i): relacao linear confirma modelo de potencia",
x = "x = log(periodo)", y = "log(casos)") +
tema
g1 + g2 +
plot_annotation(title = "Analise exploratoria — AIDS Australia (1984-1988)",
subtitle = "Esquerda: escala original | Direita: escala log-log")
4.2 Verificacao numerica
# Regressao linear de log(y) em x = log(i)
lm_explorat <- lm(log(y) ~ x, data = dados)
cf_exp <- coef(lm_explorat)
r2_exp <- summary(lm_explorat)$r.squared
print(r2_exp)[1] 0.9528825A regressão linear de \(\log(y)\) em \(x = \log(i)\) fornece \(R^2 = 0.9529\) — a relação log-log e quase perfeitamente linear, confirmando que o modelo de potência e adequado para estes dados.
5 Ajuste do Modelo via glm()
5.1 Especificação
\[ Y_i \sim \text{Poisson}(\lambda_i), \qquad \log(\lambda_i) = \beta_1 + \beta_2 x_i, \qquad x_i = \log(i) \]
modelo <- glm(
y ~ x, # formula: y em funcao de x = log(i)
family = poisson(link = "log"), # distribuicao Poisson, ligacao log
data = dados
)5.2 Saida completa do summary()
summary(modelo)
Call:
glm(formula = y ~ x, family = poisson(link = "log"), data = dados)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.99601 0.16971 5.869 4.39e-09 ***
x 1.32661 0.06463 20.525 < 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: 677.264 on 19 degrees of freedom
Residual deviance: 21.757 on 18 degrees of freedom
AIC: 138.06
Number of Fisher Scoring iterations: 46 Interpretacao da Saida do summary()
6.1 Os coeficientes
A saida do summary() mostra uma tabela de coeficientes com quatro colunas: estimativa, erro padrao, estatistica \(z\) e p-valor. Veja o que cada uma significa:
| Parametro | Estimativa | Erro Padrao | z (Wald) | p-valor | exp(b) | |
|---|---|---|---|---|---|---|
| (Intercept) | b1 (Intercept) | 0.99601 | 0.16971 | 5.869 | 0.0000 | 2.7075 |
| x | b2 (x = log i) | 1.32661 | 0.06463 | 20.525 | 1.29e-93 | 3.7682 |
6.1.1 Interpretação do intercepto \(b_1\)
Quando \(x = \log(i) = 0\), ou seja, quando \(i = 1\) (primeiro trimestre de 1984):
\[ \log(\hat\lambda_1) = b_1 \quad \Rightarrow \quad \hat\lambda_1 = e^{b_1} = e^{0.99601} = 2.7075 \]
O modelo estima em média \(2.71\) casos no primeiro trimestre. O valor observado foi \(y_1 = 1\) — a proximidade confirma um bom ponto de partida.
6.1.2 Interpretação da inclinação \(b_2\)
O modelo na forma multiplicativa é \(\lambda_i = e^{b_1} \cdot i^{b_2}\). O coeficiente \(b_2\) e o expoente da potência: os casos crescem como \(i\) elevado a \(b_2 = 1.32661\). Para cada aumento de 1 unidade em log(i), o logaritmo da média esperada de casos aumenta aproximadamente 1,32661 unidades
Interpretação prática: quando o periodo dobra (\(i \to 2i\)), os casos são multiplicados por:
\[ 2^{b_2} = 2^{1.32661} = 2.5081 \]
Ou seja, dobrar o tempo é associado a um aumento de aproximadamente 150.8% no número esperado de casos.
6.1.3 Estatística \(z\) de Wald
A estatística \(z\) testa \(H_0: \beta_j = 0\):
\[ z = \frac{\hat\beta_j}{\text{EP}(\hat\beta_j)} \;\overset{\text{assint.}}{\sim}\; N(0,1) \]
- \(|z| > 1{,}96 \Rightarrow p < 0{,}05\) (significativo ao nível de 5%)
- Para \(b_2\): \(z = 20.525\), \(|z| \gg 1{,}96\), \(p = 1.29e-93\) — b2 é fortemente significativo.
6.2 Deviance nula e deviance residual
A saida do summary() reporta duas deviances que tem papeis distintos:
| Nome | Formula | Valor | gl | Significado |
|---|---|---|---|---|
| Null deviance (deviance nula) | D(M0) = 2[l(M_sat) - l(M0)] | 677.2640 | 19 | Quao mal M0 (sem x) descreve os dados. Quanto maior, mais x pode ajudar. |
| Residual deviance (deviance residual) | D(M1) = 2[l(M_sat) - l(M1)] | 21.7573 | 18 | Quao mal M1 (com x) descreve os dados. Quanto menor, melhor o ajuste. |
| Diferenca deltaD | D(M0) - D(M1) | 655.5068 | 1 | Melhoria ao incluir x. Testa H0: b2 = 0. |
A covariavel \(x = \log(i)\) explica:
\[ \frac{\Delta D}{D(M_0)} = \frac{655.51}{677.26} = 96.8\% \text{ da deviance nula} \]
7 Qualidade e Bondade de Ajuste
7.1 Criterio 1 — Deviance residual vs. \(\chi^2(\text{gl})\)
Sob o modelo correto, \(D(M_1) \sim \chi^2(N - p)\). Portanto \(E[D(M_1)] = N - p = \text{gl}\).
| Item | Valor |
|---|---|
| Deviance residual D(M1) | 21.7573 |
| Graus de liberdade (gl = N-p) | 18 |
| Valor esperado E[chi2(gl)] | 18 |
| D(M1) / gl (regra pratica) | 1.2087 |
| p-valor P(chi2(gl) >= D) | 0.2429 |
Regra pratica: \(D/\text{gl} \approx 1\) indica bom ajuste. Aqui \(D/\text{gl} = 1.2087\) — proximo de 1, sem sinal de falta de ajuste.
Teste formal: \(p = 0.2429 > 0{,}05\) — nao ha evidencia de falta de ajuste.
7.2 Criterio 2 — Necessidade do preditor (\(\Delta D\))
mod_nulo <- glm(y ~ 1, family = poisson(link = "log"), data = dados)
av <- anova(mod_nulo, modelo, test = "Chisq")| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 19 | 677.2640 | NA | NA | NA |
| 18 | 21.7573 | 1 | 655.5068 | 0 |
\(\Delta D = 655.5068\), gl \(= 1\), \(p = 1.42e-144\): forte evidência de que o preditor \(x = \log(i)\) é indispensável.
7.3 Criterio 3 — Superdispersao
Para a Poisson, a variância teórica é \(\text{var}(Y_i) = \lambda_i\). O parametro de dispersao é:
\[ \hat\phi = \frac{D(M_1)}{\text{gl}} = \frac{21.7573}{18} = 1.2087 \]
| Item | Resultado |
|---|---|
| phi_hat = D/gl | 1.2087 |
| Interpretacao | phi ~ 1: sem superdispersao relevante |
| Acao recomendada | Manter family = poisson (adequado) |
7.4 Criterio 4 — Valores ajustados e residuos
lambda_hat <- fitted(modelo)
residuo_pad <- residuals(modelo, type = "pearson")
residuo_dev <- residuals(modelo, type = "deviance")
tab_res <- data.frame(
periodo = i_seq,
ano_trim = paste0(ano, "-T", trim),
y = y,
lambda_hat = round(lambda_hat, 2),
res_bruto = round(y - lambda_hat, 2),
res_pearson= round(residuo_pad, 3),
res_dev = round(residuo_dev, 3)
)
kable(tab_res,
caption = "**Tabela de ajuste** — observados, ajustados e residuos.",
col.names = c("i","Periodo","y","lambda_hat",
"y - lambda_hat","Pearson ri","Deviance di"),
align = "c") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 11) |>
row_spec(which(abs(residuo_pad) > 2), background = "#FFF3CD") |>
footnote(
general = sprintf(
"X2 de Pearson = sum(ri^2) = %.4f | Deviance D = %.4f (devem ser proximos)",
sum(residuo_pad^2), D_res),
general_title = "", footnote_as_chunk = TRUE
)| i | Periodo | y | lambda_hat | y - lambda_hat | Pearson ri | Deviance di |
|---|---|---|---|---|---|---|
| 1 | 1984-T1 | 1 | 2.71 | -1.71 | -1.038 | -1.193 |
| 2 | 1984-T2 | 6 | 6.79 | -0.79 | -0.303 | -0.309 |
| 3 | 1984-T3 | 16 | 11.63 | 4.37 | 1.282 | 1.212 |
| 4 | 1984-T4 | 23 | 17.03 | 5.97 | 1.446 | 1.372 |
| 5 | 1985-T1 | 27 | 22.90 | 4.10 | 0.857 | 0.833 |
| 6 | 1985-T2 | 39 | 29.17 | 9.83 | 1.821 | 1.731 |
| 7 | 1985-T3 | 31 | 35.78 | -4.78 | -0.799 | -0.818 |
| 8 | 1985-T4 | 30 | 42.72 | -12.72 | -1.946 | -2.056 |
| 9 | 1986-T1 | 43 | 49.94 | -6.94 | -0.982 | -1.006 |
| 10 | 1986-T2 | 51 | 57.44 | -6.44 | -0.849 | -0.866 |
| 11 | 1986-T3 | 63 | 65.18 | -2.18 | -0.269 | -0.271 |
| 12 | 1986-T4 | 70 | 73.15 | -3.15 | -0.368 | -0.371 |
| 13 | 1987-T1 | 88 | 81.34 | 6.66 | 0.739 | 0.729 |
| 14 | 1987-T2 | 97 | 89.75 | 7.25 | 0.765 | 0.755 |
| 15 | 1987-T3 | 91 | 98.36 | -7.36 | -0.742 | -0.751 |
| 16 | 1987-T4 | 104 | 107.14 | -3.14 | -0.304 | -0.305 |
| 17 | 1988-T1 | 110 | 116.11 | -6.11 | -0.567 | -0.572 |
| 18 | 1988-T2 | 113 | 125.27 | -12.27 | -1.096 | -1.115 |
| 19 | 1988-T3 | 149 | 134.57 | 14.43 | 1.244 | 1.223 |
| 20 | 1988-T4 | 159 | 144.05 | 14.95 | 1.246 | 1.225 |
| X2 de Pearson = sum(ri^2) = 21.6679 | Deviance D = 21.7573 (devem ser proximos) |
7.5 Criterio 5 — Graficos de diagnostico
df_diag <- data.frame(
i = i_seq,
y = y,
lambda_hat = lambda_hat,
res_pearson = residuo_pad,
res_dev = residuo_dev
)
# A: observados vs. ajustados no tempo
gA <- ggplot(df_diag, aes(x = i)) +
geom_point(aes(y = y), color = "blue4", size = 3, alpha = 0.85) +
geom_line(aes(y = lambda_hat), color = "red4", linewidth = 1.1) +
scale_x_continuous(breaks = c(1,5,9,13,17),
labels = c("1984-T1","1985-T1","1986-T1","1987-T1","1988-T1")) +
labs(title = "A) Observados vs. Ajustados",
subtitle = "Pontos: y obs. | Curva: lambda ajustado",
x = "Periodo", y = "Casos") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
# B: residuos de Pearson vs. periodo
gB <- ggplot(df_diag, aes(x = i, y = res_pearson)) +
geom_point(color = "blue4", size = 3, alpha = 0.85) +
geom_hline(yintercept = c(-2,0,2),
linetype = c("dashed","solid","dashed"),
color = c("red4","gray50","red4"),
linewidth= c(0.8, 0.7, 0.8)) +
scale_x_continuous(breaks = c(1,5,9,13,17),
labels = c("1984-T1","1985-T1","1986-T1","1987-T1","1988-T1")) +
labs(title = "B) Residuos de Pearson vs. Periodo",
subtitle = "Sem padrao sistematico = bom ajuste",
x = "Periodo", y = "Residuo de Pearson") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
# C: residuos vs. valores ajustados
gC <- ggplot(df_diag, aes(x = lambda_hat, y = res_pearson)) +
geom_point(color = "blue4", size = 3, alpha = 0.85) +
geom_hline(yintercept = c(-2,0,2),
linetype = c("dashed","solid","dashed"),
color = c("red4","gray50","red4"),
linewidth= c(0.8, 0.7, 0.8)) +
labs(title = "C) Residuos de Pearson vs. Ajustados",
subtitle = "Padrao em leque = heterocedasticidade",
x = "Valores ajustados (lambda_hat)", y = "Residuo de Pearson") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
# D: QQ-plot dos residuos de deviance
df_qq <- data.frame(
teorico = qnorm(ppoints(length(residuo_dev))),
observado = sort(residuo_dev)
)
gD <- ggplot(df_qq, aes(x = teorico, y = observado)) +
geom_point(color = "blue4", size = 3, alpha = 0.85) +
geom_abline(intercept = 0, slope = 1,
color = "red4", linetype = "dashed", linewidth = 0.9) +
labs(title = "D) QQ-plot dos residuos de deviance",
subtitle = "Pontos proximos a reta = residuos ~ N(0,1)",
x = "Quantis teoricos N(0,1)", y = "Residuos de deviance") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
(gA + gB) / (gC + gD) +
plot_annotation(
title = "Diagnosticos do modelo Poisson-log — AIDS Australia",
subtitle = "A: ajuste no tempo | B: residuos no tempo | C: residuos vs mu | D: normalidade"
)
Leitura dos graficos:
- A) A curva ajustada acompanha bem os pontos — o modelo captura o crescimento exponencial.
- B) Sem padrao sistematico no tempo e sem residuos alem de \(\pm 2\) — bom sinal.
- C) Sem padrao em leque — nao ha evidencia de heterocedasticidade nem superdispersao.
- D) Pontos proximos a diagonal — residuos de deviance aproximadamente \(N(0,1)\).
7.6 Criterio 6 — AIC
aic_nulo <- AIC(mod_nulo)
aic_model <- AIC(modelo)| Modelo | AIC | Delta AIC vs. nulo |
|---|---|---|
| M0: nulo (sem preditor) | 791.5620 | 0.0 |
| M1: com x = log(i) | 138.0552 | -653.5 |
\(\Delta\text{AIC} = -653.5\): o modelo com \(x\) e muito superior ao modelo nulo em termos de equilibrio ajuste/parcimonia.
8 Resumo dos Criterios de Qualidade
| Criterio | Valor | Diagnostico |
|---|---|---|
| 1. D(M1)/gl (regra pratica) | 1.2087 | Bom (proximo de 1) |
| 2. p-valor bondade de ajuste | 0.2429 | Bom (p > 0.05) |
| 3. Necessidade de x (deltaD) | 1.42e-144 | x fortemente necessario |
| 4. Superdispersao (phi_hat) | 1.2087 | Sem superdispersao |
| 5. Residuos fora de +/-2 | 0 de 20 | Nenhum (ideal) |
| 6. AIC (vs. modelo nulo) | 138.1 (delta = -653.5) | M1 muito melhor que M0 |
9 Interpretacao Final
9.1 Modelo ajustado
\[ \log(\hat\lambda_i) = 0.99601 + 1.32661 \log(i) \quad \Longleftrightarrow \quad \hat\lambda_i = 2.7075 \times i^{1.32661} \]
9.2 Resumo dos resultados
| Item | Resultado |
|---|---|
| Intercepto b1 | 0.99601 (EP = 0.16971, p = 0.0000) |
| Inclinacao b2 | 1.32661 (EP = 0.06463, p = 1.29e-93) |
| Numero esperado de casos em i=1 | exp(b1) = exp(0.99601) = 2.7075 |
| Fator de multiplicacao ao dobrar i | 2^b2 = 2^1.32661 = 2.5081 (+150.8%) |
| Percentual da deviance nula explicado por x | 96.8% |
| Deviance residual D(M1) | 21.7573 (gl = 18) |
| D(M1)/gl (bondade de ajuste) | 1.2087 (proximo de 1 = bom) |
| p-valor (D(M1) ~ chi2(gl)) | 0.2429 (> 0.05 = sem falta de ajuste) |
| Superdispersao phi_hat | 1.2087 (< 1.5 = sem superdispersao) |
| AIC | 138.0552 (vs. 791.5620 do modelo nulo) |
9.3 Conclusão
O modelo de Poisson com ligacao log e covariavel \(x_i = \log(i)\) descreve muito bem o crescimento dos casos de AIDS na Australia no periodo 1984–1988:
- O crescimento e consistente com um modelo de potência \(\hat\lambda_i = 2.7075 \times i^{1.32661}\)
- O preditor \(x = \log(i)\) e fortemente necessário (\(\Delta D = 655.51\), \(p \approx 0\)) e explica 96.8% da deviance nula
- O modelo ajusta bem os dados (\(D/\text{gl} = 1.2087\), \(p = 0.2429\)) sem evidência de falta de ajuste
- Nâo há superdispersão relevante (\(\hat\phi = 1.2087\)), justificando o uso da Poisson padrão
9.4 Referencias
Dobson, A. J. & Barnett, A. G. (2018). An Introduction to Generalized Linear Models (4ª ed.). CRC Press / Chapman & Hall.