x <- c(-1, -1, 0, 0, 0, 0, 1, 1, 1)
y <- c( 2, 3, 6, 7, 8, 9, 10, 12, 15)Modelos Lineares Generalizados
Apostila — Exemplos 5.6.3 e 6.7.1: Deviance de Poisson e Polinômios Fracionários
Nota ao leitor. Esta apostila implementa passo a passo dois exemplos de An Introduction to Generalized Linear Models (Dobson & Barnett, 4ª ed., 2018): o Exemplo 5.6.3 (deviance para o modelo de Poisson, com comparação de modelos aninhados) e o Exemplo 6.7.1 (associação não-linear entre comprimento de título e número de autores, via modelos linear, quadrático e polinômio fracionário). Todo o código R é exibido e comentado.
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 Exemplo 5.6.3 — Deviance para o Modelo de Poisson
1.1 Contexto e Dados
Este exemplo (Tabela 5.1 do livro) retoma os dados artificiais da Seção 4.4, onde uma reta foi ajustada a respostas de Poisson. Há 9 observações com um preditor categórico \(x \in \{-1, 0, 1\}\) representando três grupos.
kable(
data.frame(Observacao = 1:9, x = x, y = y),
caption = "**Tabela 5.1** — Dados do Exemplo 5.6.3 (Tabela 4.3 do livro).",
col.names = c("Observacao", "x", "y"),
align = "c"
) |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13)| Observacao | x | y |
|---|---|---|
| 1 | -1 | 2 |
| 2 | -1 | 3 |
| 3 | 0 | 6 |
| 4 | 0 | 7 |
| 5 | 0 | 8 |
| 6 | 0 | 9 |
| 7 | 1 | 10 |
| 8 | 1 | 12 |
| 9 | 1 | 15 |
1.2 O Modelo Ajustado
O modelo especificado é:
\[ E(Y_i) = \mu_i = \beta_1 + \beta_2\, x_i, \qquad Y_i \sim \text{Poisson}(\mu_i) \]
A ligação identidade é usada aqui (não a canônica logarítmica), portanto o preditor linear é diretamente a média. Como a ligação não-canônica pode causar problemas de convergência, é necessário fornecer um chute inicial para o IRLS.
mod_poisson <- glm(
y ~ x,
family = poisson(link = "identity"), # ligacao identidade
start = c(7, 5) # chute inicial necessario
)
summary(mod_poisson)
Call:
glm(formula = y ~ x, family = poisson(link = "identity"), start = c(7,
5))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.4516 0.8841 8.428 < 2e-16 ***
x 4.9353 1.0892 4.531 5.86e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 18.4206 on 8 degrees of freedom
Residual deviance: 1.8947 on 7 degrees of freedom
AIC: 40.008
Number of Fisher Scoring iterations: 3b1 <- coef(mod_poisson)[1] # intercepto beta_1
b2 <- coef(mod_poisson)[2] # inclinacao beta_2
y_hat <- fitted(mod_poisson) # valores ajustadosAs estimativas \(\hat\beta_1 = 7.45163\) e \(\hat\beta_2 = 4.9353\) coincidem com os valores do livro (\(b_1 = 7{,}45163\), \(b_2 = 4{,}93530\)).
1.3 Cálculo Manual da Deviance
Para o modelo de Poisson, a deviance é dada pela expressão (livro, p. 92):
\[ D = 2\sum_{i=1}^{N}\left[y_i\log\frac{y_i}{\hat{y}_i} - (y_i - \hat{y}_i)\right] \]
Como \(\sum y_i = \sum \hat{y}_i\) para a maioria dos modelos de Poisson, os termos lineares se cancelam e a expressão simplifica para \(D = 2\sum o_i \log(o_i/e_i)\).
# Contribuicao de cada observacao
contrib <- y * log(y / y_hat) - (y - y_hat)
tab51 <- data.frame(
i = 1:9,
x = x,
y = y,
y_hat = round(y_hat, 5),
contrib = round(contrib, 5)
)
D_manual <- 2 * sum(contrib)kable(
tab51,
caption = "**Tabela 5.1 reproduzida** — Contribuicoes para a deviance do modelo de Poisson.",
col.names = c("i", "xi", "yi", "y_hat_i", "yi * log(yi/y_hat_i) - (yi - y_hat_i)"),
align = "c"
) |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) |>
footnote(
general = sprintf("D = 2 x %.5f = %.4f (livro: D = 1.8947)",
sum(contrib), D_manual),
general_title = "", footnote_as_chunk = TRUE
)| i | xi | yi | y_hat_i | yi * log(yi/y_hat_i) - (yi - y_hat_i) |
|---|---|---|---|---|
| 1 | -1 | 2 | 2.51633 | 0.05702 |
| 2 | -1 | 3 | 2.51633 | 0.04376 |
| 3 | 0 | 6 | 7.45163 | 0.15159 |
| 4 | 0 | 7 | 7.45163 | 0.01397 |
| 5 | 0 | 8 | 7.45163 | 0.01970 |
| 6 | 0 | 9 | 7.45163 | 0.15076 |
| 7 | 1 | 10 | 12.38693 | 0.24636 |
| 8 | 1 | 12 | 12.38693 | 0.00611 |
| 9 | 1 | 15 | 12.38693 | 0.25805 |
| D = 2 x 0.94733 = 1.8947 (livro: D = 1.8947) |
Verificacao:
deviance(mod_poisson)= 1.8947 — identico ao valor manual.
1.4 Inferência: Bondade de Ajuste
Sob o modelo correto, \(D \sim \chi^2(N - p)\). Aqui \(N = 9\) e \(p = 2\), logo \(D \sim \chi^2(7)\).
N <- length(y)
p <- 2
gl <- N - p # 7
p_valor <- pchisq(D_manual, df = gl, lower.tail = FALSE)tab_inf <- data.frame(
Estatistica = c("Deviance D", "Graus de liberdade (N-p)", "E[chi2(7)]",
"D / gl", "p-valor"),
Valor = c(
sprintf("%.4f", D_manual),
as.character(gl),
as.character(gl),
sprintf("%.4f", D_manual / gl),
sprintf("%.4f", p_valor)
)
)
kable(tab_inf,
caption = "**Inferência sobre a bondade de ajuste** do modelo de Poisson.",
col.names = c("Estatistica", "Valor"), align = "lr") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13)| Estatistica | Valor |
|---|---|
| Deviance D | 1.8947 |
| Graus de liberdade (N-p) | 7 |
| E[chi2(7)] | 7 |
| D / gl | 0.2707 |
| p-valor | 0.9654 |
Interpretacao: \(D = 1{,}89 \ll \text{gl} = 7\) (esperado \(= 7\)); \(D/\text{gl} = 0{,}27\); \(p = 0{,}97 \gg 0{,}05\). O modelo linear ajusta muito bem — não há evidência de falta de ajuste.
df_chi <- data.frame(d = seq(0.1, 22, by = 0.05))
df_chi$dens <- dchisq(df_chi$d, df = gl)
crit_5 <- qchisq(0.95, gl)
ggplot(df_chi, aes(x = d, y = dens)) +
geom_line(color = "#4C72B0", linewidth = 1.1) +
geom_area(data = subset(df_chi, d >= crit_5),
fill = "#C44E52", alpha = 0.25) +
geom_vline(xintercept = D_manual, color = "#2CA02C",
linetype = "dashed", linewidth = 1.2) +
annotate("text", x = D_manual + 0.4, y = 0.07,
label = sprintf("D = %.2f\np = %.4f", D_manual, p_valor),
color = "#2CA02C", size = 3.8, hjust = 0) +
annotate("text", x = crit_5 + 0.4, y = 0.012,
label = sprintf("critico\n%.2f", crit_5),
color = "#C44E52", size = 3.5, hjust = 0) +
labs(title = "Bondade de ajuste — D vs. chi2(7)",
subtitle = "Area: regiao critica (5%) | Tracejado: D observado",
x = "d", y = "Densidade") +
tema
1.5 Resíduos Padronizados de Pearson
Para Poisson com \(\text{var}(Y_i) = \mu_i\), o resíduo padronizado de Pearson é:
\[ r_i = \frac{y_i - \hat{\mu}_i}{\sqrt{\hat{\mu}_i}} \]
r_pearson <- (y - y_hat) / sqrt(y_hat)
tab_res <- data.frame(
i = 1:9,
y = y,
y_hat = round(y_hat, 4),
residuo = round(y - y_hat, 4),
r_p = round(r_pearson, 4)
)
kable(tab_res,
caption = "**Residuos padronizados de Pearson** — Exemplo 5.6.3.",
col.names = c("i", "yi", "y_hat_i", "yi - y_hat_i", "ri (Pearson)"),
align = "c") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) |>
footnote(
general = sprintf("X2 de Pearson = sum(ri^2) = %.4f (similar a D = %.4f)",
sum(r_pearson^2), D_manual),
general_title = "", footnote_as_chunk = TRUE
)| i | yi | y_hat_i | yi - y_hat_i | ri (Pearson) |
|---|---|---|---|---|
| 1 | 2 | 2.5163 | -0.5163 | -0.3255 |
| 2 | 3 | 2.5163 | 0.4837 | 0.3049 |
| 3 | 6 | 7.4516 | -1.4516 | -0.5318 |
| 4 | 7 | 7.4516 | -0.4516 | -0.1654 |
| 5 | 8 | 7.4516 | 0.5484 | 0.2009 |
| 6 | 9 | 7.4516 | 1.5484 | 0.5672 |
| 7 | 10 | 12.3869 | -2.3869 | -0.6782 |
| 8 | 12 | 12.3869 | -0.3869 | -0.1099 |
| 9 | 15 | 12.3869 | 2.6131 | 0.7425 |
| X2 de Pearson = sum(ri^2) = 1.8944 (similar a D = 1.8947) |
df_diag <- data.frame(y = y, y_hat = y_hat, x = x, r = r_pearson)
p1 <- ggplot(df_diag, aes(x = y_hat, y = y)) +
geom_point(color = "#4C72B0", size = 3) +
geom_abline(intercept = 0, slope = 1,
linetype = "dashed", color = "gray50") +
labs(title = "Observados vs. Ajustados",
subtitle = "Linha: y = y_hat (ajuste perfeito)",
x = "Valores ajustados", y = "Valores observados") +
tema
p2 <- ggplot(df_diag, aes(x = x, y = r)) +
geom_point(color = "#4C72B0", size = 3) +
geom_hline(yintercept = c(-2, 0, 2),
linetype = c("dashed","solid","dashed"),
color = c("#C44E52","gray50","#C44E52")) +
scale_x_continuous(breaks = c(-1, 0, 1)) +
labs(title = "Residuos padronizados vs. x",
subtitle = "Linhas vermelhas: +/- 2",
x = "x", y = "Residuo padronizado") +
tema
p1 + p2
Diagnostico: Todos os residuos estao dentro de \([-2, +2]\); nenhuma observacao e discrepante. O grafico de observados vs. ajustados nao mostra padrao sistematico.
1.6 Comparação de Modelos Aninhados via ANOVA
1.6.1 As Duas Perguntas
Ao comparar modelos aninhados via deviance, ha duas perguntas distintas, cada uma com sua propria comparacao:
tab_pq <- data.frame(
Pergunta = c("A — Bondade de ajuste", "B — Necessidade do preditor"),
Objetivo = c("M1 descreve bem os dados?",
"O preditor x e necessario?"),
Comparacao = c("M1 vs. M_sat (saturado)",
"M0 (nulo) vs. M1 (completo)"),
Estatistica = c("D(M1) ~ chi2(N-p)", "deltaD = D(M0)-D(M1) ~ chi2(1)"),
`p bom` = c("GRANDE (D pequeno = bom ajuste)",
"PEQUENO (deltaD grande = x importante)")
)
kable(tab_pq,
caption = "**As duas perguntas** ao comparar modelos de Poisson aninhados.",
col.names = c("Pergunta","Objetivo","Comparacao","Estatistica","p-valor bom e"),
align = "l") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = TRUE, font_size = 12) |>
column_spec(1, bold = TRUE)| Pergunta | Objetivo | Comparacao | Estatistica | p-valor bom e |
|---|---|---|---|---|
| A — Bondade de ajuste | M1 descreve bem os dados? | M1 vs. M_sat (saturado) | D(M1) ~ chi2(N-p) | GRANDE (D pequeno = bom ajuste) |
| B — Necessidade do preditor | O preditor x e necessario? | M0 (nulo) vs. M1 (completo) | deltaD = D(M0)-D(M1) ~ chi2(1) | PEQUENO (deltaD grande = x importante) |
1.6.2 Ajuste dos Três Modelos
# M0: modelo nulo — E(Y_i) = mu (media constante)
mod_nulo <- glm(y ~ 1,
family = poisson(link = "identity"),
start = c(mean(y)))
# M1: modelo completo — E(Y_i) = b1 + b2*x (ja ajustado)
# M_sat: modelo saturado — um parametro por observacao
# Deviance = 0 por definicao; serve de referencia para D(M0) e D(M1)
obs_id <- factor(seq_along(y))
mod_sat <- glm(y ~ obs_id,
family = poisson(link = "identity"),
start = c(mean(y), rep(0, length(y) - 1)))tab_mod <- data.frame(
Modelo = c("M0 (nulo)", "M1 (completo)", "M_sat (saturado)"),
Formula = c("E(Yi) = mu", "E(Yi) = b1 + b2*x", "E(Yi) = mu_i"),
Parametros = c(1, 2, 9),
gl = c(df.residual(mod_nulo), df.residual(mod_poisson), df.residual(mod_sat)),
Deviance = c(round(deviance(mod_nulo), 4),
round(deviance(mod_poisson), 4),
round(deviance(mod_sat), 4))
)
kable(tab_mod,
caption = "**Os tres modelos** — hierarquia M0 c M1 c M_sat.",
col.names = c("Modelo","Formula","Parametros","gl","Deviance"),
align = "lllrr") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13) |>
row_spec(3, background = "#f5f5f5")| Modelo | Formula | Parametros | gl | Deviance |
|---|---|---|---|---|
| M0 (nulo) | E(Yi) = mu | 1 | 8 | 18.4206 |
| M1 (completo) | E(Yi) = b1 + b2*x | 2 | 7 | 1.8947 |
| M_sat (saturado) | E(Yi) = mu_i | 9 | 0 | 0.0000 |
1.6.3 Pergunta A — Bondade de Ajuste: M1 vs. M_sat
Esta e a comparacao do Exemplo 5.6.3 (ja calculada na Secao 1.3). A deviance reportada pelo summary() e exatamente \(D(M1) = 2[\ell(M_{\text{sat}}) - \ell(M_1)]\).
D_M1 <- deviance(mod_poisson)
gl_M1 <- df.residual(mod_poisson) # 7
p_M1 <- pchisq(D_M1, df = gl_M1, lower.tail = FALSE)
# Via anova()
av_sat <- anova(mod_poisson, mod_sat, test = "Chisq")kable(av_sat,
caption = "**Pergunta A** — anova(mod_poisson, mod_sat): M1 vs. M_sat.",
digits = 4) |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12)| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 7 | 1.8947 | NA | NA | NA |
| 0 | 0.0000 | 7 | 1.8947 | 0.9654 |
Resultado A: \(D(M1) = 1.8947\), \(\text{gl} = 7\), \(p = 0.9654\). p-valor grande \(\Rightarrow\) M1 nao e significativamente pior que o saturado — bom ajuste.
1.6.4 Pergunta B — Necessidade do Preditor: M0 vs. M1
D_M0 <- deviance(mod_nulo)
gl_M0 <- df.residual(mod_nulo) # 8
delta_D <- D_M0 - D_M1
gl_delta <- gl_M0 - gl_M1 # 1
p_delta <- pchisq(delta_D, df = gl_delta, lower.tail = FALSE)
# Via anova()
av_nulo <- anova(mod_nulo, mod_poisson, test = "Chisq")kable(av_nulo,
caption = "**Pergunta B** — anova(mod_nulo, mod_poisson): M0 vs. M1.",
digits = 4) |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12)| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 8 | 18.4206 | NA | NA | NA |
| 7 | 1.8947 | 1 | 16.526 | 0 |
Resultado B: \(\Delta D = D(M0) - D(M1) = 18.4206 - 1.8947 = 16.526\), \(\text{gl} = 1\), \(p = 4.80e-05\). p-valor muito pequeno \(\Rightarrow\) o preditor \(x\) e indispensavel.
1.6.5 Tabela ANOVA Unificada
tab_anova <- data.frame(
Fonte = c("M0: nulo",
"Melhoria B: efeito de x (M0->M1)",
"M1: completo, residuo (Pergunta A)",
"M_sat: saturado"),
Par = c(1, 1, 2, 9),
gl = c(gl_M0, gl_delta, gl_M1, 0),
Deviance = c(round(D_M0, 4), round(delta_D, 4), round(D_M1, 4), 0),
Ref = c(sprintf("chi2(%d)", gl_M0), sprintf("chi2(%d)", gl_delta),
sprintf("chi2(%d)", gl_M1), "— ref."),
p = c(sprintf("%.4f", pchisq(D_M0, gl_M0, lower.tail=FALSE)),
sprintf("%.2e", p_delta),
sprintf("%.4f", p_M1),
"—"),
Q = c("—", "B", "A", "—")
)
kable(tab_anova,
caption = "**Tabela ANOVA de deviances** — Exemplo 5.6.3. Coluna Q indica qual pergunta cada linha responde.",
col.names = c("Fonte","Par.","gl","Deviance","Distrib. H0","p-valor","Q"),
align = "lrrrrrc") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) |>
row_spec(2, background = "#FFF3CD") |>
row_spec(3, background = "#D4EDDA")| Fonte | Par. | gl | Deviance | Distrib. H0 | p-valor | Q |
|---|---|---|---|---|---|---|
| M0: nulo | 1 | 8 | 18.4206 | chi2(8) | 0.0183 | — |
| Melhoria B: efeito de x (M0->M1) | 1 | 1 | 16.5260 | chi2(1) | 4.80e-05 | B |
| M1: completo, residuo (Pergunta A) | 2 | 7 | 1.8947 | chi2(7) | 0.9654 | A |
| M_sat: saturado | 9 | 0 | 0.0000 | — ref. | — | — |
1.6.6 Gráfico Comparativo das Duas Perguntas
x_grid <- seq(0.01, 25, by = 0.05)
# Painel A
df_A <- data.frame(d = x_grid, dens = dchisq(x_grid, df = gl_M1))
crit_A <- qchisq(0.95, gl_M1)
gA <- ggplot(df_A, aes(x = d, y = dens)) +
geom_line(color = "#4C72B0", linewidth = 1.1) +
geom_area(data = subset(df_A, d >= crit_A), fill = "#4C72B0", alpha = 0.20) +
geom_vline(xintercept = D_M1, color = "#4C72B0",
linetype = "dashed", linewidth = 1.2) +
geom_vline(xintercept = crit_A, color = "#4C72B0",
linetype = "dotted", linewidth = 0.9) +
annotate("text", x = D_M1 + 0.3, y = 0.075,
label = sprintf("D(M1)=%.2f\np=%.4f", D_M1, p_M1),
color = "#4C72B0", size = 3.5, hjust = 0) +
coord_cartesian(xlim = c(0, 20), ylim = c(0, 0.16)) +
labs(title = "A) Bondade de ajuste: M1 vs. M_sat",
subtitle = sprintf("D(M1)=%.4f ~ chi2(%d) | p=%.4f (GRANDE = bom)", D_M1, gl_M1, p_M1),
x = "d", y = "Densidade") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", color = "#4C72B0"))
# Painel B
df_B <- data.frame(d = x_grid, dens = dchisq(x_grid, df = gl_delta))
crit_B <- qchisq(0.95, gl_delta)
gB <- ggplot(df_B, aes(x = d, y = dens)) +
geom_line(color = "#C44E52", linewidth = 1.1) +
geom_area(data = subset(df_B, d >= crit_B), fill = "#C44E52", alpha = 0.20) +
geom_vline(xintercept = delta_D, color = "#C44E52",
linetype = "dashed", linewidth = 1.2) +
geom_vline(xintercept = crit_B, color = "#C44E52",
linetype = "dotted", linewidth = 0.9) +
annotate("text", x = delta_D + 0.3, y = 0.22,
label = sprintf("deltaD=%.2f\np=%.2e", delta_D, p_delta),
color = "#C44E52", size = 3.5, hjust = 0) +
coord_cartesian(xlim = c(0, 20), ylim = c(0, 0.50)) +
labs(title = "B) Necessidade do preditor: M0 vs. M1",
subtitle = sprintf("deltaD=%.4f ~ chi2(%d) | p=%.2e (PEQUENO = x necessario)", delta_D, gl_delta, p_delta),
x = "d", y = "Densidade") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", color = "#C44E52"))
gA + gB +
plot_annotation(
title = "ANOVA de deviances — Modelo Poisson (Exemplo 5.6.3)",
subtitle = "Tracejado: deviance observada | Pontilhado: valor critico chi2 (95%) | Area: regiao critica"
)
1.6.7 Interpretação Final
O modelo linear M1 e ao mesmo tempo adequado e necessario:
tab_interp <- data.frame(
Comparacao = c("A: M1 vs. M_sat", "B: M0 vs. M1"),
Estatistica = c(sprintf("D(M1) = %.4f ~ chi2(%d)", D_M1, gl_M1),
sprintf("deltaD = %.4f ~ chi2(%d)", delta_D, gl_delta)),
p_valor = c(sprintf("%.4f", p_M1), sprintf("%.2e", p_delta)),
Interpretacao = c(
"p GRANDE: M1 ajusta bem os dados. A relacao linear em x e adequada.",
"p PEQUENO: x melhora muito o ajuste. O preditor e necessario."
)
)
kable(tab_interp,
caption = "**Resumo das duas comparacoes** — Exemplo 5.6.3.",
col.names = c("Comparacao","Estatistica","p-valor","Interpretacao"),
align = "llll") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = TRUE, font_size = 12) |>
column_spec(3, bold = TRUE) |>
row_spec(1, background = "#D4EDDA") |>
row_spec(2, background = "#FFF3CD")| Comparacao | Estatistica | p-valor | Interpretacao |
|---|---|---|---|
| A: M1 vs. M_sat | D(M1) = 1.8947 ~ chi2(7) | 0.9654 | p GRANDE: M1 ajusta bem os dados. A relacao linear em x e adequada. |
| B: M0 vs. M1 | deltaD = 16.5260 ~ chi2(1) | 4.80e-05 | p PEQUENO: x melhora muito o ajuste. O preditor e necessario. |
Regra pratica para interpretar p-valores em deviances: - Bondade de ajuste (vs. saturado): p grande e bom — o modelo nao e pior que o saturado. - Teste do preditor (vs. nulo): p pequeno e bom — o preditor melhora significativamente o ajuste.
2 Exemplo 6.7.1 — PLOS Medicine: Comprimento de Titulo vs. Numero de Autores
2.1 Contexto e Dados Simulados
O livro analisa 878 artigos da revista PLOS Medicine (2011–2015). A variavel resposta e o comprimento do titulo (em caracteres) e a variavel explicativa e o numero de autores (truncado em 30). Como os dados nao estao disponiveis publicamente, simulamos uma amostra com os parametros exatos reportados nas Tabelas 6.19 e 6.20.
set.seed(2018)
N_plos <- 878
# Distribuicao do numero de autores (log-normal, truncada em [1, 30])
autores <- pmin(pmax(round(rlnorm(N_plos, meanlog = 2.2, sdlog = 0.7)), 1), 30)
# Comprimento do titulo gerado pelo modelo quadratico do livro + ruido
titulo <- 81.40 + 6.07 * autores - 0.15 * autores^2 + rnorm(N_plos, sd = 30)
titulo <- pmax(titulo, 20)
dados_plos <- data.frame(authors = autores, nchar = titulo,
x_scaled = autores / 30)desc <- data.frame(
Variavel = c("Numero de autores", "Comprimento do titulo (char)"),
Min = c(min(autores), round(min(titulo), 1)),
Med = c(median(autores), round(median(titulo), 1)),
Max = c(max(autores), round(max(titulo), 1)),
Media = c(round(mean(autores), 1), round(mean(titulo), 1))
)
kable(desc,
caption = "**Estatisticas descritivas** — dados simulados do Exemplo 6.7.1 (N = 878).",
col.names = c("Variavel","Min","Mediana","Max","Media"),
align = "lrrrr") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13)| Variavel | Min | Mediana | Max | Media |
|---|---|---|---|---|
| Numero de autores | 1.0 | 9.0 | 30.0 | 11.0 |
| Comprimento do titulo (char) | 21.2 | 123.3 | 253.1 | 121.7 |
2.2 Modelo Linear
\[ E(Y_i) = \beta_0 + \beta_1 x_i \tag{6.15} \]
mod_lin <- lm(nchar ~ authors, data = dados_plos)
summary(mod_lin)
Call:
lm(formula = nchar ~ authors, data = dados_plos)
Residuals:
Min 1Q Median 3Q Max
-89.24 -20.40 -0.85 21.46 117.47
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 102.6258 1.9272 53.25 <2e-16 ***
authors 1.7366 0.1462 11.88 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 31.52 on 876 degrees of freedom
Multiple R-squared: 0.1388, Adjusted R-squared: 0.1378
F-statistic: 141.1 on 1 and 876 DF, p-value: < 2.2e-16cf_lin <- coef(summary(mod_lin))
kable(round(cf_lin, 4),
caption = "**Modelo linear** — estimativas e erros padrao.",
align = "r") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 102.6258 | 1.9272 | 53.2499 | 0 |
| authors | 1.7366 | 0.1462 | 11.8800 | 0 |
Interpretacao: Cada autor adicional esta associado, em media, a \(1.74\) caracteres a mais no titulo. O modelo assume que essa associacao e constante para qualquer numero de autores.
2.3 Modelo Quadratico
\[ E(Y_i) = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 \tag{6.16} \]
mod_quad <- lm(nchar ~ authors + I(authors^2), data = dados_plos)
summary(mod_quad)
Call:
lm(formula = nchar ~ authors + I(authors^2), data = dados_plos)
Residuals:
Min 1Q Median 3Q Max
-81.941 -18.595 0.566 20.635 108.919
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 77.25578 3.18831 24.231 <2e-16 ***
authors 6.63196 0.52234 12.697 <2e-16 ***
I(authors^2) -0.16368 0.01684 -9.722 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 29.96 on 875 degrees of freedom
Multiple R-squared: 0.2227, Adjusted R-squared: 0.2209
F-statistic: 125.4 on 2 and 875 DF, p-value: < 2.2e-16cf_quad <- coef(summary(mod_quad))
kable(round(cf_quad, 4),
caption = "**Modelo quadratico** — estimativas e erros padrao.",
align = "r") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 77.2558 | 3.1883 | 24.2309 | 0 |
| authors | 6.6320 | 0.5223 | 12.6967 | 0 |
| I(authors^2) | -0.1637 | 0.0168 | -9.7223 | 0 |
b0q <- coef(mod_quad)[1]; b1q <- coef(mod_quad)[2]; b2q <- coef(mod_quad)[3]
pico <- -b1q / (2 * b2q)Interpretacao: O crescimento inicial e mais rapido (\(\hat\beta_1 = 6.63\)) mas desacelera com mais autores (\(\hat\beta_2 = -0.164 < 0\)). O pico estimado esta em \(x^* = -\hat\beta_1 / (2\hat\beta_2) = 20.3\) autores.
2.3.1 Tabela 6.19 Reproduzida
tab619 <- data.frame(
Parametro = c("b0", "b1", "b2"),
Lin_est = c(round(coef(mod_lin)[1], 2), round(coef(mod_lin)[2], 2), NA),
Lin_ep = c(round(cf_lin[1, 2], 2), round(cf_lin[2, 2], 2), NA),
Quad_est = c(round(b0q, 2), round(b1q, 2), round(b2q, 3)),
Quad_ep = c(round(cf_quad[1, 2], 2), round(cf_quad[2, 2], 2), round(cf_quad[3, 2], 3))
)
kable(tab619, na_label = "—",
caption = "**Tabela 6.19 reproduzida** — Estimativas e erros padrao para os modelos linear (6.15) e quadratico (6.16).",
col.names = c("Param.",
"Estim. linear", "EP linear",
"Estim. quadrat.", "EP quadrat."),
align = "lrrrr") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13)| Param. | Estim. linear | EP linear | Estim. quadrat. | EP quadrat. | |
|---|---|---|---|---|---|
| (Intercept) | b0 | 102.63 | 1.93 | 77.260 | 3.190 |
| authors | b1 | 1.74 | 0.15 | 6.630 | 0.520 |
| b2 | NA | NA | -0.164 | 0.017 |
2.3.2 Tabela 6.20 — ANOVA: Linear vs. Quadratico
av <- anova(mod_lin, mod_quad)
tab620 <- data.frame(
Modelo = c("Linear","Quadratico","Diferenca"),
gl_res = c(df.residual(mod_lin), df.residual(mod_quad), 1),
SQR_res = c(round(sum(resid(mod_lin)^2)),
round(sum(resid(mod_quad)^2)),
round(av$`Sum of Sq`[2])),
F_stat = c(NA, NA, round(av$F[2], 3)),
p_valor = c(NA, NA, format.pval(av$`Pr(>F)`[2], digits = 3))
)
kable(tab620, na_label = "—",
caption = "**Tabela 6.20 reproduzida** — ANOVA para os modelos linear e quadratico.",
col.names = c("Modelo","gl residual","SQR residual","F","p-valor"),
align = "lrrrr") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13) |>
row_spec(3, bold = TRUE, background = "#D4EDDA")| Modelo | gl residual | SQR residual | F | p-valor |
|---|---|---|---|---|
| Linear | 876 | 870127 | NA | NA |
| Quadratico | 875 | 785295 | NA | NA |
| Diferenca | 1 | 84832 | 94.523 | <2e-16 |
Conclusao: O modelo quadratico e significativamente melhor (\(F = 94.5\), \(p < 0{,}001\)). O termo \(x^2\) e necessario.
2.4 Polinomios Fracionarios — Tabela 6.21
O livro propoe polinomios fracionarios com candidatos \(p \in \{-2, -1, -0{,}5, 0, 0{,}5, 1, 2, 3\}\) (com \(p = 0\) correspondendo a \(\log x\)), aplicados a \(x_{\text{scaled}} = \text{autores}/30\).
p_vals <- c(-2, -1, -0.5, 0, 0.5, 1, 2, 3)
# Calcular SQR para cada potencia p
calc_sqr <- function(p, dados) {
xp <- if (p == 0) log(dados$x_scaled) else dados$x_scaled^p
sum(resid(lm(nchar ~ xp, data = dados))^2)
}
sqr_vals <- sapply(p_vals, calc_sqr, dados = dados_plos)
p_melhor <- p_vals[which.min(sqr_vals)]sqr_ref <- c(924091, 868424, 846761, 852395, 890891, 943484, 1024999, 1065298)
tab621 <- data.frame(
p = p_vals,
SQR_sim = round(sqr_vals),
SQR_livro = sqr_ref,
Melhor = ifelse(p_vals == p_melhor, "*", "")
)
kable(tab621,
caption = "**Tabela 6.21 reproduzida** — SQR para polinomios fracionarios (x escalado = autores/30). (*) indica o melhor ajuste.",
col.names = c("p","SQR (simulado)","SQR (livro)",""),
align = "rrrr") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13) |>
row_spec(which(p_vals == p_melhor), bold = TRUE, background = "#D4EDDA") |>
footnote(general = "Livro usa dados reais; SQR simulados diferem um pouco dos originais.",
general_title = "", footnote_as_chunk = TRUE)| p | SQR (simulado) | SQR (livro) | |
|---|---|---|---|
| -2.0 | 930317 | 924091 | |
| -1.0 | 837534 | 868424 | |
| -0.5 | 809478 | 846761 | * |
| 0.0 | 810826 | 852395 | |
| 0.5 | 835604 | 890891 | |
| 1.0 | 870127 | 943484 | |
| 2.0 | 929975 | 1024999 | |
| 3.0 | 964058 | 1065298 | |
| Livro usa dados reais; SQR simulados diferem um pouco dos originais. |
Melhor transformacao: tanto na simulacao quanto no livro, \(p = -0{,}5\) (raiz quadrada reciproca) produz o menor SQR. A curva correspondente mostra crescimento inicial muito rapido e estabilizacao para muitos autores — diferente do comportamento simetrico do quadratico.
2.5 Selecao Automatica com mfp
O pacote mfp busca automaticamente a melhor potencia \(p\) entre os candidatos via criterio de maxima verossimilhanca.
library(mfp)
mod_mfp <- mfp(
nchar ~ fp(x_scaled, df = 2, select = 1),
data = dados_plos,
family = gaussian
)
print(mod_mfp)Call:
mfp(formula = nchar ~ fp(x_scaled, df = 2, select = 1), data = dados_plos,
family = gaussian)
Deviance table:
Resid. Dev
Null model 1010316
Linear model 870127
Final model 809478.2
Fractional polynomials:
df.initial select alpha df.final power1 power2
x_scaled 2 1 0.05 2 -0.5 .
Transformations of covariates:
formula
x_scaled I(x_scaled^-0.5)
Coefficients:
Intercept x_scaled.1
163.18 -21.16
Degrees of Freedom: 877 Total (i.e. Null); 876 Residual
Null Deviance: 1010000
Residual Deviance: 809500 AIC: 8491 sqr_mfp <- sum(resid(mod_mfp)^2)
fp_power <- mod_mfp$fptransform[[1]]$power
cat(sprintf("\nSQR do modelo mfp = %.0f\n", sqr_mfp))
SQR do modelo mfp = 809478cat(sprintf("Potencia selecionada: p = %s\n", paste(fp_power, collapse = ", ")))Potencia selecionada: p = 2.6 Grafico com os Tres Ajustes
x_grade <- seq(1, 30, by = 0.5)
df_curvas <- data.frame(
x = rep(x_grade, 3),
y = c(
predict(mod_lin, newdata = data.frame(authors = x_grade)),
predict(mod_quad, newdata = data.frame(authors = x_grade)),
predict(mod_mfp, newdata = data.frame(x_scaled = x_grade / 30))
),
Modelo = rep(c(
sprintf("Linear (SQR=%.0f)", sum(resid(mod_lin)^2)),
sprintf("Quadratico (SQR=%.0f)", sum(resid(mod_quad)^2)),
sprintf("Frac. pol. p=%s (SQR=%.0f)", paste(fp_power, collapse=","), sqr_mfp)
), each = length(x_grade))
)
medias_obs <- aggregate(nchar ~ authors, data = dados_plos, FUN = mean)
ggplot() +
geom_point(data = medias_obs, aes(x = authors, y = nchar),
color = "gray35", size = 2.5, alpha = 0.8) +
geom_line(data = df_curvas,
aes(x = x, y = y, color = Modelo, linetype = Modelo),
linewidth = 1.1) +
scale_color_manual(values = c("#4C72B0","#C44E52","#55A868"), name = "Modelo") +
scale_linetype_manual(values = c("dashed","solid","dotdash"), name = "Modelo") +
labs(title = "Exemplo 6.7.1 — PLOS Medicine",
subtitle = "Pontos: medias observadas | Curvas: modelos ajustados",
x = "Numero de autores",
y = "Comprimento medio do titulo (caracteres)") +
tema +
theme(legend.position = "bottom",
legend.direction = "vertical")
2.7 Interpretacao Final e Comparacao dos Tres Modelos
tab_final <- data.frame(
Modelo = c("Linear","Quadratico","Frac. polinomial (mfp)"),
SQR = c(round(sum(resid(mod_lin)^2)),
round(sum(resid(mod_quad)^2)),
round(sqr_mfp)),
gl_res = c(df.residual(mod_lin), df.residual(mod_quad), df.residual(mod_mfp)),
Forma = c("Reta (crescimento constante)",
"Parabola (pico e queda)",
"Curva assimetrica (rapido inicio, estabilizacao)"),
Interpretabilidade = c("Alta","Media","Baixa")
)
kable(tab_final,
caption = "**Comparacao dos tres modelos** — Exemplo 6.7.1.",
col.names = c("Modelo","SQR","gl res.","Forma da curva","Interpretabilidade"),
align = "lrrlc") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = TRUE, font_size = 12) |>
column_spec(1, bold = TRUE) |>
row_spec(3, background = "#D4EDDA")| Modelo | SQR | gl res. | Forma da curva | Interpretabilidade |
|---|---|---|---|---|
| Linear | 870127 | 876 | Reta (crescimento constante) | Alta |
| Quadratico | 785295 | 875 | Parabola (pico e queda) | Media |
| Frac. polinomial (mfp) | 809478 | 876 | Curva assimetrica (rapido inicio, estabilizacao) | Baixa |
O modelo quadratico ja e uma melhoria substancial sobre o linear (Tabela 6.20: \(F = 94.5\), \(p < 0{,}001\)). O polinomio fracionario com \(p = -0{,}5\) ajusta ainda melhor ao capturar a assimetria real da relacao: crescimento muito rapido para poucos autores e estabilizacao para muitos. Na pratica, o modelo quadratico e preferivel pela sua interpretabilidade mais simples, a nao ser que o melhor ajuste do polinomio fracionario seja substantivamente importante.
3 Resumo Geral
tab_res <- data.frame(
Exemplo = c("5.6.3 — Bondade de ajuste",
"5.6.3 — Teste do preditor",
"6.7.1 — Linear vs. Quadratico",
"6.7.1 — Melhor polinomio fracionario"),
Funcao = c("deviance(mod_poisson) + pchisq()",
"anova(mod_nulo, mod_poisson, test='Chisq')",
"anova(mod_lin, mod_quad)",
"mfp(nchar ~ fp(x_scaled, df=2))"),
Resultado = c(
sprintf("D=%.4f, gl=7, p=%.4f (bom ajuste)", D_M1, p_M1),
sprintf("deltaD=%.4f, gl=1, p=%.2e (x necessario)", delta_D, p_delta),
sprintf("F=%.1f, p<0.001 (quadratico melhor)", av$F[2]),
sprintf("p=%s, SQR=%.0f (melhor ajuste)", paste(fp_power,collapse=","), sqr_mfp)
)
)
kable(tab_res,
caption = "**Resumo geral** dos resultados dos dois exemplos.",
col.names = c("Exemplo / Comparacao","Funcao R usada","Resultado"),
align = "lll") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = TRUE, font_size = 12) |>
column_spec(1, bold = TRUE, width = "30%")| Exemplo / Comparacao | Funcao R usada | Resultado |
|---|---|---|
| 5.6.3 — Bondade de ajuste | deviance(mod_poisson) + pchisq() | D=1.8947, gl=7, p=0.9654 (bom ajuste) |
| 5.6.3 — Teste do preditor | anova(mod_nulo, mod_poisson, test='Chisq') | deltaD=16.5260, gl=1, p=4.80e-05 (x necessario) |
| 6.7.1 — Linear vs. Quadratico | anova(mod_lin, mod_quad) | F=94.5, p<0.001 (quadratico melhor) |
| 6.7.1 — Melhor polinomio fracionario | mfp(nchar ~ fp(x_scaled, df=2)) | p=, SQR=809478 (melhor ajuste) |
3.1 Referencias
Dobson, A. J. & Barnett, A. G. (2018). An Introduction to Generalized Linear Models (4ª ed.). CRC Press / Chapman & Hall.
Royston, P., Altman, D. G. & Sauerbrei, W. (1999). Dichotomizing continuous predictors in multiple regression. Statistics in Medicine, 18(8), 945–974.