Quase-verossimilhança e Modelo Binomial Negativo com aplicações em R
Sobre esta apostila. O modelo de Poisson assume que média e variância são iguais. Na prática, dados de contagem quase sempre violam essa suposição — fenômeno chamado superdispersão. Esta apostila apresenta o diagnóstico de superdispersão e duas soluções práticas: o quasi-Poisson e a Binomial Negativa, aplicadas ao banco de dados
quinedo pacote MASS.
quine)O banco quine contém dados de 146 crianças australianas do ensino fundamental. Registra o número de dias de ausência escolar (Days) e quatro variáveis explicativas categóricas.
| Variável | Tipo | Descrição |
|---|---|---|
Days |
Resposta (contagem) | Dias ausentes no ano letivo |
Eth |
Fator (A/N) | Etnia: A = aborígine, N = não-aborígine |
Sex |
Fator (F/M) | Sexo: feminino / masculino |
Age |
Fator (F0/F1/F2/F3) | Faixa etária / série escolar |
Lrn |
Fator (AL/SL) | Ritmo de aprendizagem: AL = médio, SL = lento |
library(MASS) # para o banco quine e para glm.nb()
library(ggplot2) # para os gráficos
data(quine)
# Visão geral do banco
str(quine)'data.frame': 146 obs. of 5 variables:
$ Eth : Factor w/ 2 levels "A","N": 1 1 1 1 1 1 1 1 1 1 ...
$ Sex : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
$ Age : Factor w/ 4 levels "F0","F1","F2",..: 1 1 1 1 1 1 1 1 2 2 ...
$ Lrn : Factor w/ 2 levels "AL","SL": 2 2 2 1 1 1 1 1 2 2 ...
$ Days: int 2 11 14 5 5 13 20 22 6 6 ...# Estatísticas da variável resposta
cat("=== Estatísticas de Days (dias ausentes) ===\n\n")=== Estatísticas de Days (dias ausentes) ===cat(sprintf("N = %d observações\n", nrow(quine)))N = 146 observaçõescat(sprintf("Mínimo = %d\n", min(quine$Days)))Mínimo = 0cat(sprintf("Mediana = %.1f\n", median(quine$Days)))Mediana = 11.0cat(sprintf("Média = %.3f\n", mean(quine$Days)))Média = 16.459cat(sprintf("Máximo = %d\n", max(quine$Days)))Máximo = 81cat(sprintf("Variância = %.3f\n", var(quine$Days)))Variância = 264.167cat(sprintf("\nVar / Média = %.2f\n", var(quine$Days) / mean(quine$Days)))
Var / Média = 16.05cat("(Poisson exige Var/Média ≈ 1)\n")(Poisson exige Var/Média ≈ 1)# Frequências por grupo
cat("\n=== Frequências por variável ===\n\n")
=== Frequências por variável ===cat("Etnia:\n"); print(table(quine$Eth))Etnia:
A N
69 77 cat("\nSexo:\n"); print(table(quine$Sex))
Sexo:
F M
80 66 cat("\nSérie:\n"); print(table(quine$Age))
Série:
F0 F1 F2 F3
27 46 40 33 cat("\nRitmo:\n"); print(table(quine$Lrn))
Ritmo:
AL SL
83 63 ggplot(quine, aes(x = Days)) +
geom_histogram(binwidth = 5, fill = "#4C72B0",
color = "white", alpha = 0.85) +
geom_vline(xintercept = mean(quine$Days),
color = "#C44E52", linewidth = 1.2, linetype = "dashed") +
annotate("text", x = mean(quine$Days) + 3, y = 28,
label = paste0("média = ", round(mean(quine$Days), 1)),
color = "#C44E52", size = 4) +
labs(title = "Distribuição dos dias ausentes — banco quine",
subtitle = "Cauda longa e variância muito maior que a média indicam superdispersão",
x = "Dias ausentes (Days)", y = "Frequência") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
ggplot(quine, aes(x = Age, y = Days, fill = Lrn)) +
geom_boxplot(alpha = 0.7, outlier.shape = 16) +
facet_wrap(~ Eth + Sex, labeller = label_both) +
scale_fill_manual(values = c("#4C72B0","#C44E52"),
name = "Ritmo\n(AL=médio,SL=lento)") +
labs(title = "Dias ausentes por série, ritmo, etnia e sexo",
subtitle = "Alta variabilidade dentro dos grupos confirma superdispersão",
x = "Série (Age)", y = "Dias ausentes") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
O que os gráficos mostram: a variância/média \(\approx 20\) (muito acima de 1) e a cauda longa à direita indicam superdispersão forte. O modelo de Poisson padrão não é adequado para esses dados.
\[\text{var}(Y_i) = E(Y_i) = \lambda_i\]
Na superdispersão, a variância real é maior que a média:
\[\text{var}(Y_i) = \phi\,\lambda_i, \qquad \phi > 1\]
Consequência prática: ignorar a superdispersão produz erros padrão subestimados por um fator \(\sqrt{\hat\phi}\), tornando os testes \(z\) e \(p\)-valores incorretamente otimistas.
\[\hat\phi = \frac{X^2}{N - p} = \frac{\sum(y_i - \hat\mu_i)^2/\hat\mu_i}{N - p}\]
ou equivalentemente via deviance: \(\hat\phi = D/(N-p)\).
Quasi-Poisson: mantém os coeficientes do Poisson e multiplica os erros padrão por \(\sqrt{\hat\phi}\). Não especifica uma distribuição completa — apenas a relação \(\text{var}(Y) = \phi\,\mu\).
Binomial Negativa: especifica uma distribuição completa com parâmetro extra \(\theta\):
\[\text{var}(Y_i) = \mu_i + \frac{\mu_i^2}{\theta}\]
Quanto menor \(\theta\), maior a superdispersão. Quando \(\theta \to \infty\), converge para Poisson.
Note que os erros-padrão são todos muito próximos de zero. Na prática, ao observarmos dados reais, este detalhe já acende um alerta para possível superdirpersão, o que consequentemente pode estar levando os valores-p erroneamente a serem extremamente significativos.
# Ajuste do modelo de Poisson padrão
mod_pois <- glm(Days ~ Eth + Sex + Age + Lrn,
family = poisson(link = "log"),
data = quine)
summary(mod_pois)
Call:
glm(formula = Days ~ Eth + Sex + Age + Lrn, family = poisson(link = "log"),
data = quine)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.71538 0.06468 41.980 < 2e-16 ***
EthN -0.53360 0.04188 -12.740 < 2e-16 ***
SexM 0.16160 0.04253 3.799 0.000145 ***
AgeF1 -0.33390 0.07009 -4.764 1.90e-06 ***
AgeF2 0.25783 0.06242 4.131 3.62e-05 ***
AgeF3 0.42769 0.06769 6.319 2.64e-10 ***
LrnSL 0.34894 0.05204 6.705 2.02e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 2073.5 on 145 degrees of freedom
Residual deviance: 1696.7 on 139 degrees of freedom
AIC: 2299.2
Number of Fisher Scoring iterations: 5# Estatísticas de dispersão
X2 <- sum(residuals(mod_pois, type = "pearson")^2)
gl <- df.residual(mod_pois)
phi <- X2 / gl
cat("=== Diagnóstico de superdispersão ===\n\n")=== Diagnóstico de superdispersão ===cat(sprintf("X² de Pearson = %.2f\n", X2))X² de Pearson = 1830.19cat(sprintf("Graus de liberdade = %d\n", gl))Graus de liberdade = 139cat(sprintf("phi_hat = X²/gl = %.2f\n\n", phi))phi_hat = X²/gl = 13.17cat(sprintf("Os erros padrão do Poisson estão subestimados\n"))Os erros padrão do Poisson estão subestimadoscat(sprintf("por um fator de sqrt(%.2f) = %.2f\n", phi, sqrt(phi)))por um fator de sqrt(13.17) = 3.63\(\hat\phi = 13,17 \gg 1\). Isso significa que os erros padrão do Poisson estão subestimados por um fator de \(\sqrt{\hat\phi} \approx 3,6\). Os valores-\(p\) estão incorretos — precisamos corrigir.
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
# Gráfico 1: resíduos de Pearson vs. valores ajustados
plot(log(fitted(mod_pois)),
residuals(mod_pois, type = "pearson"),
xlab = "log(valores ajustados)",
ylab = "Resíduos de Pearson",
main = "Resíduos vs. Ajustados\n(Poisson)",
pch = 16, col = "#4C72B088")
abline(h = c(-2, 0, 2), lty = c(2, 1, 2),
col = c("red","gray50","red"))
# Gráfico 2: QQ-plot dos resíduos de deviance
qqnorm(residuals(mod_pois, type = "deviance"),
main = "QQ-plot — Resíduos Poisson",
pch = 16, col = "#4C72B088")
qqline(residuals(mod_pois, type = "deviance"),
col = "red", lwd = 1.5)
par(mfrow = c(1, 1))Esquerda (padrão em leque): os resíduos se espalham cada vez mais à medida que os valores ajustados crescem. Isso é a assinatura visual da superdispersão — a variância real cresce mais rápido do que o modelo Poisson prevê.
Direita (QQ-plot): os pontos nas extremidades saem muito da reta diagonal, especialmente no extremo superior. Isso indica caudas mais pesadas do que a distribuição normal, consistente com superdispersão forte.
O quasi-Poisson é a solução mais simples: muda apenas family = quasipoisson. Os coeficientes são idênticos ao Poisson, mas os erros padrão são corrigidos.
# Modelo quasi-Poisson — mudar apenas a family
mod_quasi <- glm(Days ~ Eth + Sex + Age + Lrn,
family = quasipoisson(link = "log"),
data = quine)
summary(mod_quasi)
Call:
glm(formula = Days ~ Eth + Sex + Age + Lrn, family = quasipoisson(link = "log"),
data = quine)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.7154 0.2347 11.569 < 2e-16 ***
EthN -0.5336 0.1520 -3.511 0.000602 ***
SexM 0.1616 0.1543 1.047 0.296914
AgeF1 -0.3339 0.2543 -1.313 0.191413
AgeF2 0.2578 0.2265 1.138 0.256938
AgeF3 0.4277 0.2456 1.741 0.083831 .
LrnSL 0.3489 0.1888 1.848 0.066760 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 13.16691)
Null deviance: 2073.5 on 145 degrees of freedom
Residual deviance: 1696.7 on 139 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5# Extrair erros padrão dos dois modelos
ep_pois <- coef(summary(mod_pois))[, "Std. Error"]
ep_quasi <- coef(summary(mod_quasi))[, "Std. Error"]
cat("=== Erros padrão: Poisson vs. Quasi-Poisson ===\n\n")=== Erros padrão: Poisson vs. Quasi-Poisson ===cat(sprintf("%-20s %12s %12s %12s\n",
"Parâmetro", "EP Poisson", "EP Quasi", "Razão"))Parâmetro EP Poisson EP Quasi Razãocat(strrep("-", 58), "\n")---------------------------------------------------------- for (nm in names(ep_pois)) {
cat(sprintf("%-20s %12.4f %12.4f %12.4f\n",
nm, ep_pois[nm], ep_quasi[nm],
ep_quasi[nm]/ep_pois[nm]))
}(Intercept) 0.0647 0.2347 3.6286
EthN 0.0419 0.1520 3.6286
SexM 0.0425 0.1543 3.6286
AgeF1 0.0701 0.2543 3.6286
AgeF2 0.0624 0.2265 3.6286
AgeF3 0.0677 0.2456 3.6286
LrnSL 0.0520 0.1888 3.6286cat(sprintf("\nRazão média = %.4f | sqrt(phi) = %.4f\n",
mean(ep_quasi/ep_pois), sqrt(phi)))
Razão média = 3.6286 | sqrt(phi) = 3.6286# Modelo Binomial Negativa — usar glm.nb() do pacote MASS
mod_nb <- glm.nb(Days ~ Eth + Sex + Age + Lrn,
data = quine)
summary(mod_nb)
Call:
glm.nb(formula = Days ~ Eth + Sex + Age + Lrn, data = quine,
init.theta = 1.274892646, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.89458 0.22842 12.672 < 2e-16 ***
EthN -0.56937 0.15333 -3.713 0.000205 ***
SexM 0.08232 0.15992 0.515 0.606710
AgeF1 -0.44843 0.23975 -1.870 0.061425 .
AgeF2 0.08808 0.23619 0.373 0.709211
AgeF3 0.35690 0.24832 1.437 0.150651
LrnSL 0.29211 0.18647 1.566 0.117236
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(1.2749) family taken to be 1)
Null deviance: 195.29 on 145 degrees of freedom
Residual deviance: 167.95 on 139 degrees of freedom
AIC: 1109.2
Number of Fisher Scoring iterations: 1
Theta: 1.275
Std. Err.: 0.161
2 x log-likelihood: -1093.151 # Parâmetro de dispersão theta
theta <- mod_nb$theta
cat(sprintf("Parâmetro theta = %.4f\n", theta))Parâmetro theta = 1.2749cat(sprintf("var(Y) = mu + mu²/%.4f\n\n", theta))var(Y) = mu + mu²/1.2749cat("Exemplos de variância prevista pela BN:\n")Exemplos de variância prevista pela BN:for (mu_ex in c(5, 10, 20, 50)) {
cat(sprintf(" mu = %2d => var = %.1f (Poisson teria var = %d)\n",
mu_ex, mu_ex + mu_ex^2/theta, mu_ex))
} mu = 5 => var = 24.6 (Poisson teria var = 5)
mu = 10 => var = 88.4 (Poisson teria var = 10)
mu = 20 => var = 333.8 (Poisson teria var = 20)
mu = 50 => var = 2010.9 (Poisson teria var = 50)# Reunir coeficientes dos três modelos em uma tabela simples
b_p <- coef(mod_pois)
b_q <- coef(mod_quasi)
b_nb <- coef(mod_nb)
ep_p <- coef(summary(mod_pois))[, 2]
ep_q <- coef(summary(mod_quasi))[, 2]
ep_nb <- coef(summary(mod_nb))[, 2]
cat("=== Coeficientes e erros padrão — três modelos ===\n\n")=== Coeficientes e erros padrão — três modelos ===cat(sprintf("%-18s %8s %6s %8s %6s %8s %6s\n",
"", "Pois b", "EP", "Quasi b", "EP", "BN b", "EP")) Pois b EP Quasi b EP BN b EPcat(strrep("-", 72), "\n")------------------------------------------------------------------------ for (nm in names(b_p)) {
cat(sprintf("%-18s %8.4f %6.4f %8.4f %6.4f %8.4f %6.4f\n",
nm,
b_p[nm], ep_p[nm],
b_q[nm], ep_q[nm],
b_nb[nm], ep_nb[nm]))
}(Intercept) 2.7154 0.0647 2.7154 0.2347 2.8946 0.2284
EthN -0.5336 0.0419 -0.5336 0.1520 -0.5694 0.1533
SexM 0.1616 0.0425 0.1616 0.1543 0.0823 0.1599
AgeF1 -0.3339 0.0701 -0.3339 0.2543 -0.4484 0.2397
AgeF2 0.2578 0.0624 0.2578 0.2265 0.0881 0.2362
AgeF3 0.4277 0.0677 0.4277 0.2456 0.3569 0.2483
LrnSL 0.3489 0.0520 0.3489 0.1888 0.2921 0.1865cat("=== Estatísticas de ajuste ===\n\n")=== Estatísticas de ajuste ===phi_quasi <- summary(mod_quasi)$dispersion
X2_nb <- sum(residuals(mod_nb, type = "pearson")^2)
phi_nb <- X2_nb / df.residual(mod_nb)
cat(sprintf("Poisson: D = %7.2f gl = %d phi_hat = %6.2f AIC = %8.1f\n",
deviance(mod_pois), df.residual(mod_pois), phi, AIC(mod_pois)))Poisson: D = 1696.71 gl = 139 phi_hat = 13.17 AIC = 2299.2cat(sprintf("Quasi-Poisson: D = %7.2f gl = %d phi_hat = %6.2f AIC = (N/A)\n",
deviance(mod_quasi), df.residual(mod_quasi), phi_quasi))Quasi-Poisson: D = 1696.71 gl = 139 phi_hat = 13.17 AIC = (N/A)cat(sprintf("Binomial Neg.: D = %7.2f gl = %d phi_res = %6.2f AIC = %8.1f\n",
deviance(mod_nb), df.residual(mod_nb), phi_nb, AIC(mod_nb)))Binomial Neg.: D = 167.95 gl = 139 phi_res = 0.99 AIC = 1109.2cat(sprintf("\ntheta (BN) = %.4f\n", theta))
theta (BN) = 1.2749cat(sprintf("AIC(BN) < AIC(Poisson): %.1f vs %.1f — BN é melhor\n",
AIC(mod_nb), AIC(mod_pois)))AIC(BN) < AIC(Poisson): 1109.2 vs 2299.2 — BN é melhordf_ep <- data.frame(
param = rep(names(ep_p), 3),
modelo = rep(c("Poisson","Quasi-Poisson","Bin. Negativa"), each=length(ep_p)),
EP = c(ep_p, ep_q, ep_nb)
)
df_ep$param <- factor(df_ep$param, levels = names(ep_p))
df_ep$modelo <- factor(df_ep$modelo,
levels = c("Poisson","Quasi-Poisson","Bin. Negativa"))
ggplot(df_ep, aes(x = param, y = EP, fill = modelo)) +
geom_col(position = "dodge", alpha = 0.85, width = 0.7) +
scale_fill_manual(values = c("#C44E52","#4C72B0","#2CA02C"),
name = "Modelo") +
labs(title = "Erros padrão por modelo",
subtitle = "Poisson (vermelho) subestima os EP — quasi e BN corrigem",
x = NULL, y = "Erro padrão") +
theme_minimal(base_size = 11) +
theme(axis.text.x = element_text(angle = 40, hjust = 1),
plot.title = element_text(face = "bold"),
legend.position = "top")
cat("=== Critério de escolha ===\n\n")=== Critério de escolha ===cat(sprintf("phi_hat (Poisson) = %.2f => superdispersão FORTE\n", phi))phi_hat (Poisson) = 13.17 => superdispersão FORTEcat(sprintf("phi_res (BN) = %.2f => BN capturou bem a dispersão\n", phi_nb))phi_res (BN) = 0.99 => BN capturou bem a dispersãocat(sprintf("AIC(Poisson) = %.1f\n", AIC(mod_pois)))AIC(Poisson) = 2299.2cat(sprintf("AIC(BN) = %.1f (menor = melhor)\n\n", AIC(mod_nb)))AIC(BN) = 1109.2 (menor = melhor)# Teste de razão de verossimilhança BN vs Poisson
lrt_stat <- 2*(as.numeric(logLik(mod_nb)) - as.numeric(logLik(mod_pois)))
p_lrt <- pchisq(lrt_stat, df = 1, lower.tail = FALSE)
cat(sprintf("Teste LRT (BN vs Poisson): stat = %.2f, p = %.2e\n", lrt_stat, p_lrt))Teste LRT (BN vs Poisson): stat = 1192.03, p = 3.29e-261cat("=> BN é SIGNIFICATIVAMENTE melhor que Poisson\n\n")=> BN é SIGNIFICATIVAMENTE melhor que Poissoncat("MODELO ESCOLHIDO: Binomial Negativa\n")MODELO ESCOLHIDO: Binomial Negativacat("(mas quasi-Poisson também é adequado quando AIC não é necessário)\n")(mas quasi-Poisson também é adequado quando AIC não é necessário)cf <- coef(mod_nb)
ep_b <- coef(summary(mod_nb))[, 2]
irr <- exp(cf)
ic_inf <- exp(cf - 1.96 * ep_b)
ic_sup <- exp(cf + 1.96 * ep_b)
z_val <- cf / ep_b
p_val <- 2 * pnorm(abs(z_val), lower.tail = FALSE)
cat("=== Tabela de resultados — Binomial Negativa ===\n\n")=== Tabela de resultados — Binomial Negativa ===cat(sprintf("%-18s %8s %8s %8s %8s %8s %8s\n",
"Parâmetro", "b", "EP", "z", "p-valor",
"IRR", "IC95%"))Parâmetro b EP z p-valor IRR IC95%cat(strrep("-", 80), "\n")-------------------------------------------------------------------------------- for (nm in names(cf)) {
sig <- ifelse(p_val[nm] < 0.001, "***",
ifelse(p_val[nm] < 0.01, "**",
ifelse(p_val[nm] < 0.05, "*",
ifelse(p_val[nm] < 0.1, ".", ""))))
cat(sprintf("%-18s %8.4f %8.4f %8.3f %8.4f %8.4f [%.3f,%.3f] %s\n",
nm, cf[nm], ep_b[nm], z_val[nm], p_val[nm],
irr[nm], ic_inf[nm], ic_sup[nm], sig))
}(Intercept) 2.8946 0.2284 12.672 0.0000 18.0759 [11.552,28.284] ***
EthN -0.5694 0.1533 -3.713 0.0002 0.5659 [0.419,0.764] ***
SexM 0.0823 0.1599 0.515 0.6067 1.0858 [0.794,1.486]
AgeF1 -0.4484 0.2397 -1.870 0.0614 0.6386 [0.399,1.022] .
AgeF2 0.0881 0.2362 0.373 0.7092 1.0921 [0.687,1.735]
AgeF3 0.3569 0.2483 1.437 0.1507 1.4289 [0.878,2.325]
LrnSL 0.2921 0.1865 1.566 0.1172 1.3392 [0.929,1.930] cat("\nSignif.: *** p<0.001 ** p<0.01 * p<0.05 . p<0.1\n")
Signif.: *** p<0.001 ** p<0.01 * p<0.05 . p<0.1\[IRR = \dfrac{\text{taxa esperada do grupo 1}}{\text{taxa esperada do grupo referência}}\]
A RC é usada no modelo Binomial (resposta binária) e a IRR nos modelos Poisson, Binomial Negativa (BN) e Quasi-Poisson (QP), em que as respostas são contagens.
O IRR na Poisson, BN e QP compara taxas médias de eventos. O RC na Binomial compara as chances (odds) de sucesso.
\(IRR \rightarrow e^{\hat{\beta}} = \mu_{1}/\mu_{2}\)
\(RC \rightarrow e^{\hat{\beta}} = odds_{1}/odds_{2}\)
Resumo: RC compara chances de um evento binário; IRR compara taxas médias de contagem. A fórmula \(e^{\hat{\beta}}\) é a mesma, mas os denominadores são completamente diferentes.
Somente os coeficientes \(\beta_{0}\) (intercepto) \(\beta_{1}\) (etnia não aborígine) foram significativos.
Normalmente, rodaríamos o modelo novamente apenas com a covariável que apresenta significância estatística. Todavia, para fins didáticos, faremos as análises de todas covariáveis.
cat("=== INTERPRETAÇÃO DOS COEFICIENTES (modelo BN) ===\n\n")=== INTERPRETAÇÃO DOS COEFICIENTES (modelo BN) ===cat("INTERCEPTO:\n")INTERCEPTO:cat(sprintf(" Criança de referência: etnia N, sexo F, série F0, ritmo AL\n")) Criança de referência: etnia N, sexo F, série F0, ritmo ALcat(sprintf(" Número esperado de faltas = exp(%.4f) = %.1f dias\n\n",
cf["(Intercept)"], exp(cf["(Intercept)"]))) Número esperado de faltas = exp(2.8946) = 18.1 diascat("ETNIA (EthN = não-aborígine vs. aborígine):\n")ETNIA (EthN = não-aborígine vs. aborígine):cat(sprintf(" b = %.4f => IRR = exp(%.4f) = %.4f\n",
cf["EthN"], cf["EthN"], exp(cf["EthN"]))) b = -0.5694 => IRR = exp(-0.5694) = 0.5659cat(sprintf(" Não-aborígines faltam %.0f%% %s que aborígines\n",
abs(exp(cf["EthN"])-1)*100,
ifelse(exp(cf["EthN"])>1,"mais","menos"))) Não-aborígines faltam 43% menos que aboríginescat(sprintf(" p = %.4f %s\n\n", p_val["EthN"],
ifelse(p_val["EthN"]<0.05,"(significativo)","(não significativo)"))) p = 0.0002 (significativo)cat("SEXO (SexM = masculino vs. feminino):\n")SEXO (SexM = masculino vs. feminino):cat(sprintf(" b = %.4f => IRR = exp(%.4f) = %.4f\n",
cf["SexM"], cf["SexM"], exp(cf["SexM"]))) b = 0.0823 => IRR = exp(0.0823) = 1.0858cat(sprintf(" Meninos faltam %.0f%% %s que meninas\n",
abs(exp(cf["SexM"])-1)*100,
ifelse(exp(cf["SexM"])>1,"mais","menos"))) Meninos faltam 9% mais que meninascat(sprintf(" p = %.4f %s\n\n", p_val["SexM"],
ifelse(p_val["SexM"]<0.05,"(significativo)","(não significativo)"))) p = 0.6067 (não significativo)cat("SÉRIE (referência = F0):\n")SÉRIE (referência = F0):series <- c("AgeF1","AgeF2","AgeF3")
for (s in series) {
if (s %in% names(cf)) {
cat(sprintf(" %s: b = %.4f => IRR = %.4f (%.0f%% %s) p = %.4f\n",
s, cf[s], exp(cf[s]),
abs(exp(cf[s])-1)*100,
ifelse(exp(cf[s])>1,"mais","menos"),
p_val[s]))
}
} AgeF1: b = -0.4484 => IRR = 0.6386 (36% menos) p = 0.0614
AgeF2: b = 0.0881 => IRR = 1.0921 (9% mais) p = 0.7092
AgeF3: b = 0.3569 => IRR = 1.4289 (43% mais) p = 0.1507cat("\nRITMO DE APRENDIZAGEM (LrnSL = lento vs. AL = médio):\n")
RITMO DE APRENDIZAGEM (LrnSL = lento vs. AL = médio):cat(sprintf(" b = %.4f => IRR = exp(%.4f) = %.4f\n",
cf["LrnSL"], cf["LrnSL"], exp(cf["LrnSL"]))) b = 0.2921 => IRR = exp(0.2921) = 1.3392cat(sprintf(" Crianças com ritmo lento faltam %.0f%% %s\n",
abs(exp(cf["LrnSL"])-1)*100,
ifelse(exp(cf["LrnSL"])>1,"mais","menos"))) Crianças com ritmo lento faltam 34% maiscat(sprintf(" p = %.4f %s\n", p_val["LrnSL"],
ifelse(p_val["LrnSL"]<0.05,"(significativo)","(não significativo)"))) p = 0.1172 (não significativo)par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
# Gráfico 1: resíduos vs. ajustados
plot(log(fitted(mod_nb)),
residuals(mod_nb, type = "pearson"),
xlab = "log(valores ajustados)",
ylab = "Resíduos de Pearson",
main = "Resíduos vs. Ajustados\n(Binomial Negativa)",
pch = 16, col = "#2CA02C88")
abline(h = c(-2, 0, 2), lty = c(2, 1, 2),
col = c("red","gray50","red"))
# Gráfico 2: QQ-plot
qqnorm(residuals(mod_nb, type = "deviance"),
main = "QQ-plot — Resíduos BN",
pch = 16, col = "#2CA02C88")
qqline(residuals(mod_nb, type = "deviance"),
col = "red", lwd = 1.5)
par(mfrow = c(1, 1))Note que o padrão em leque (resíduos espalhando cada vez mais à medida que \(\log (\hat{\mu})\) cresce — é exatamente a assinatura visual da superdispersão) some e a distribuição fica mais simétrica com a BN.
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
# Poisson: resíduos vs. ajustados
plot(log(fitted(mod_pois)),
residuals(mod_pois, type = "pearson"),
xlab = "log(ajustados)", ylab = "Resíduo Pearson",
main = "Poisson: resíduos vs. ajustados",
pch = 16, col = "#C44E5288")
abline(h = c(-2,0,2), lty=c(2,1,2), col=c("red","gray50","red"))
# Poisson: QQ
qqnorm(residuals(mod_pois, type = "deviance"),
main = "Poisson: QQ-plot", pch = 16, col = "#C44E5288")
qqline(residuals(mod_pois, type = "deviance"),
col = "red", lwd = 1.5)
# BN: resíduos vs. ajustados
plot(log(fitted(mod_nb)),
residuals(mod_nb, type = "pearson"),
xlab = "log(ajustados)", ylab = "Resíduo Pearson",
main = "Bin. Negativa: resíduos vs. ajustados",
pch = 16, col = "#2CA02C88")
abline(h = c(-2,0,2), lty=c(2,1,2), col=c("red","gray50","red"))
# BN: QQ
qqnorm(residuals(mod_nb, type = "deviance"),
main = "Bin. Negativa: QQ-plot", pch = 16, col = "#2CA02C88")
qqline(residuals(mod_nb, type = "deviance"),
col = "red", lwd = 1.5)
par(mfrow = c(1, 1))No script da apostila foi usado Pearson no gráfico resíduos vs. ajustados porque é o mais intuitivo para mostrar o padrão em leque da superdispersão — a relação direta com \(V(\hat{\mu_{i}})\) torna o leque visualmente óbvio. Mas o correto seria usar os dois, como foi feito no QQ-plot (que usou deviance). O ideal para as análises seria:
Resíduos vs. log(ajustados): usar deviance (diagnóstico geral).
QQ-plot: usar deviance (verificar normalidade).
Detecção do padrão em leque: mostrar Pearson como complemento.
cat("=== Verificação da dispersão residual ===\n\n")=== Verificação da dispersão residual ===phi_res_nb <- sum(residuals(mod_nb, type="pearson")^2) / df.residual(mod_nb)
cat(sprintf("Poisson: phi_hat = %.2f (muito longe de 1 — RUIM)\n", phi))Poisson: phi_hat = 13.17 (muito longe de 1 — RUIM)cat(sprintf("BN: phi_res = %.2f (próximo de 1 — BOM)\n\n", phi_res_nb))BN: phi_res = 0.99 (próximo de 1 — BOM)cat("Interpretar:\n")Interpretar:cat(" phi_res ≈ 1 significa que a BN capturou toda\n") phi_res ≈ 1 significa que a BN capturou todacat(" a superdispersão dos dados.\n") a superdispersão dos dados.cat(sprintf(" phi_res = %.2f => pequena superdispersão residual\n",
phi_res_nb)) phi_res = 0.99 => pequena superdispersão residualcat(sprintf(" (fator de %.2f — muito melhor que %.2f do Poisson)\n",
phi_res_nb, phi)) (fator de 0.99 — muito melhor que 13.17 do Poisson)quine$res_nb <- residuals(mod_nb, type = "pearson")
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
boxplot(res_nb ~ Eth, data = quine,
col = c("#4C72B0","#C44E52"),
main = "Resíduos por Etnia",
xlab = "Eth", ylab = "Resíduo Pearson")
abline(h = 0, lty = 2, col = "gray50")
boxplot(res_nb ~ Sex, data = quine,
col = c("#4C72B0","#C44E52"),
main = "Resíduos por Sexo",
xlab = "Sex", ylab = "")
abline(h = 0, lty = 2, col = "gray50")
boxplot(res_nb ~ Age, data = quine,
col = c("#4C72B0","#C44E52","#55A868","#FF7F0E"),
main = "Resíduos por Série",
xlab = "Age", ylab = "Resíduo Pearson")
abline(h = 0, lty = 2, col = "gray50")
boxplot(res_nb ~ Lrn, data = quine,
col = c("#4C72B0","#C44E52"),
main = "Resíduos por Ritmo",
xlab = "Lrn", ylab = "")
abline(h = 0, lty = 2, col = "gray50")
par(mfrow = c(1, 1))Leitura dos boxplots de resíduos:
- As medianas de todos os grupos estão próximas de zero — sem viés sistemático por grupo.
- A dispersão é razoavelmente homogênea entre os grupos — o modelo captura bem as diferenças.
- Alguns valores extremos existem, mas são esperados em dados com variabilidade alta.
cat("=== SÍNTESE DOS RESULTADOS ===\n\n")=== SÍNTESE DOS RESULTADOS ===cat("1. DIAGNÓSTICO:\n")1. DIAGNÓSTICO:cat(sprintf(" Var(Days)/Média(Days) = %.1f/%.1f = %.1f >> 1\n",
var(quine$Days), mean(quine$Days),
var(quine$Days)/mean(quine$Days))) Var(Days)/Média(Days) = 264.2/16.5 = 16.1 >> 1cat(sprintf(" phi_hat (Poisson) = %.2f — superdispersão FORTE\n\n",phi)) phi_hat (Poisson) = 13.17 — superdispersão FORTEcat("2. MODELO ESCOLHIDO: Binomial Negativa (AIC = %.1f)\n\n",AIC(mod_nb))2. MODELO ESCOLHIDO: Binomial Negativa (AIC = %.1f)
1109.151cat("3. RESULTADOS SUBSTANTIVOS (modelo BN):\n\n")3. RESULTADOS SUBSTANTIVOS (modelo BN):pars <- c("EthN","SexM","AgeF1","AgeF2","AgeF3","LrnSL")
pars <- pars[pars %in% names(cf)]
for (nm in pars) {
dir <- ifelse(exp(cf[nm])>1,"mais","menos")
sig <- ifelse(p_val[nm]<0.05," (*)"," ")
cat(sprintf(" %-8s: IRR = %.3f => %.0f%% %s faltas %s\n",
nm, exp(cf[nm]), abs(exp(cf[nm])-1)*100, dir, sig))
} EthN : IRR = 0.566 => 43% menos faltas (*)
SexM : IRR = 1.086 => 9% mais faltas
AgeF1 : IRR = 0.639 => 36% menos faltas
AgeF2 : IRR = 1.092 => 9% mais faltas
AgeF3 : IRR = 1.429 => 43% mais faltas
LrnSL : IRR = 1.339 => 34% mais faltas cat("\n4. QUALIDADE DO AJUSTE:\n")
4. QUALIDADE DO AJUSTE:cat(sprintf(" phi_res (BN) = %.2f — dispersão controlada\n", phi_res_nb)) phi_res (BN) = 0.99 — dispersão controladacat(" QQ-plot e resíduos vs. ajustados: sem padrão sistemático\n") QQ-plot e resíduos vs. ajustados: sem padrão sistemáticocat(" Resíduos por grupo: medianas próximas de zero\n\n") Resíduos por grupo: medianas próximas de zerocat("5. COMPARAÇÃO:\n")5. COMPARAÇÃO:cat(sprintf(" AIC Poisson = %.1f\n", AIC(mod_pois))) AIC Poisson = 2299.2cat(sprintf(" AIC BN = %.1f (-%.1f => BN muito melhor)\n",
AIC(mod_nb), AIC(mod_pois)-AIC(mod_nb))) AIC BN = 1109.2 (-1190.0 => BN muito melhor)cat("=== FLUXO RECOMENDADO PARA DADOS DE CONTAGEM ===\n\n")=== FLUXO RECOMENDADO PARA DADOS DE CONTAGEM ===cat(" 1. Calcule Variância/Média da resposta\n") 1. Calcule Variância/Média da respostacat(" - Se ≈ 1: Poisson padrão pode ser adequado\n") - Se ≈ 1: Poisson padrão pode ser adequadocat(" - Se >> 1: há superdispersão -> siga os passos\n\n") - Se >> 1: há superdispersão -> siga os passoscat(" 2. Ajuste o Poisson e calcule phi = X2/gl\n") 2. Ajuste o Poisson e calcule phi = X2/glcat(" - phi < 1.5: Poisson OK\n") - phi < 1.5: Poisson OKcat(" - phi 1.5–5: quasi-Poisson\n") - phi 1.5–5: quasi-Poissoncat(" - phi > 5: Binomial Negativa\n\n") - phi > 5: Binomial Negativacat(" 3. Ajuste quasi-Poisson e/ou BN\n\n") 3. Ajuste quasi-Poisson e/ou BNcat(" 4. Compare AIC (apenas BN tem AIC) e phi residual\n\n") 4. Compare AIC (apenas BN tem AIC) e phi residualcat(" 5. Analise resíduos: padrão em leque? QQ-plot?\n\n") 5. Analise resíduos: padrão em leque? QQ-plot?cat(" 6. Interprete os IRR = exp(coeficiente)\n") 6. Interprete os IRR = exp(coeficiente)cat(" - IRR > 1: mais eventos (fator > 100%)\n") - IRR > 1: mais eventos (fator > 100%)cat(" - IRR < 1: menos eventos (fator < 100%)\n") - IRR < 1: menos eventos (fator < 100%)cat(" - (IRR - 1) × 100% = variação percentual\n") - (IRR - 1) × 100% = variação percentualVenables, W. N. & Ripley, B. D. (2002). Modern Applied Statistics with S (4ª ed.). Springer. [banco de dados quine]
Dobson, A. J. & Barnett, A. G. (2018). An Introduction to Generalized Linear Models (4ª ed.). CRC Press.
Cameron, A. C. & Trivedi, P. K. (2013). Regression Analysis of Count Data (2ª ed.). Cambridge University Press.