| 1.051 | 4.006 | 5.905 | 8.108 | 10.205 | 11.608 | 14.110 |
| 1.337 | 4.012 | 5.956 | 8.546 | 10.396 | 11.745 | 14.496 |
| 1.389 | 4.063 | 6.068 | 8.666 | 10.861 | 11.762 | 15.395 |
| 1.921 | 4.921 | 6.121 | 8.831 | 11.026 | 11.895 | 16.179 |
| 1.942 | 5.445 | 6.473 | 9.106 | 11.214 | 12.044 | 17.092 |
| 2.322 | 5.620 | 7.501 | 9.711 | 11.362 | 13.520 | 17.568 |
| 3.629 | 5.817 | 7.886 | 9.806 | 11.604 | 13.670 | 17.568 |
Modelos Lineares Generalizados
Apostila — Capítulo 4: Estimação
Nota ao leitor: Esta apostila resume o Capítulo 4 de An Introduction to Generalized Linear Models (Dobson & Barnett, 4ª ed., 2018). O capítulo desenvolve a teoria de estimação por máxima verossimilhança para MLGs, introduz o algoritmo de Newton–Raphson, apresenta o método de scoring e demonstra que os estimadores de MLGs são obtidos por mínimos quadrados iterativamente reponderados (IRLS). Todos os resultados são ilustrados com exemplos numéricos completos, incluindo tabelas e gráficos reproduzidos.
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 Introdução
O objetivo do capítulo é responder à pergunta: dado um MLG especificado, como estimar o vetor de parâmetros \(\boldsymbol{\beta}\)?
Em alguns casos especiais existem fórmulas fechadas (ex.: médias amostrais para Poisson, mínimos quadrados para Normal). Em geral, porém, as equações de máxima verossimilhança para MLGs não têm solução analítica e precisam ser resolvidas iterativamente.
A estrutura do capítulo é:
- Um exemplo motivador — distribuição de Weibull — para ilustrar o algoritmo de Newton–Raphson.
- A teoria geral de estimação para MLGs.
- Um exemplo completo de regressão de Poisson com ligação identidade.
2 Exemplo Motivador: Tempos de Falha de Vasos de Pressão
2.1 Os Dados
Os dados são tempos de falha (em horas) de vasos de pressão de fibra de Kevlar submetidos a 70% da carga máxima. \(N = 49\) observações.
2.2 Distribuição de Weibull
Um modelo amplamente utilizado para tempos de falha é a distribuição de Weibull:
\[ f(y;\,\lambda,\theta) = \frac{\lambda y^{\lambda-1}}{\theta^\lambda} \exp\!\left[-\left(\frac{y}{\theta}\right)^\lambda\right], \qquad y > 0 \tag{4.1} \]
onde \(\lambda\) é o parâmetro de forma e \(\theta\) é o parâmetro de escala.
Reescrevendo na forma da família exponencial:
\[ f(y;\theta) = \exp\!\bigl[\underbrace{y^\lambda}_{a(y)} \cdot \underbrace{(-\theta^{-\lambda})}_{b(\theta)} + \underbrace{(\log\lambda - \lambda\log\theta)}_{c(\theta)} + \underbrace{(\lambda-1)\log y}_{d(y)}\bigr] \]
Observação: Esta não está na forma canônica (pois \(a(y) = y^\lambda \neq y\) para \(\lambda \neq 1\)). Quando \(\lambda = 1\) obtemos a distribuição Exponencial, que é canônica. Mesmo assim, a estrutura da família exponencial facilita o cálculo da informação de Fisher.
Para este exemplo, fixamos \(\lambda = 2\) (compatível com o gráfico de probabilidade da Figura 4.2 do livro) e estimamos apenas \(\theta\).
2.3 Log-Verossimilhança e a Função Escore
Para \(\lambda = 2\) e \(N = 49\) observações independentes, a log-verossimilhança é:
\[ \ell(\theta;\,\mathbf{y}) = \sum_{i=1}^N \left[(\lambda-1)\log y_i + \log\lambda - \lambda\log\theta - \left(\frac{y_i}{\theta}\right)^\lambda\right] \]
A função escore (derivada em relação a \(\theta\)):
\[ U(\theta) = \frac{d\ell}{d\theta} = \sum_{i=1}^N \left[-\frac{\lambda}{\theta} + \frac{\lambda y_i^\lambda}{\theta^{\lambda+1}}\right] = -\frac{2N}{\theta} + \frac{2\sum y_i^2}{\theta^3} \tag{4.8} \]
A derivada do escore (necessária para Newton–Raphson):
\[ U'(\theta) = \frac{dU}{d\theta} = \frac{2N}{\theta^2} - \frac{6\sum y_i^2}{\theta^4} \tag{4.9} \]
3 O Algoritmo de Newton–Raphson
3.1 Ideia Geométrica
O objetivo é encontrar \(\hat\theta\) tal que \(U(\hat\theta) = 0\). O algoritmo de Newton–Raphson aproxima a função \(U\) por sua tangente no ponto atual e avança para onde a tangente cruza o eixo horizontal.
Derivação: Por Taylor de primeira ordem,
\[ U(\theta^{(m)}) \approx U(\theta^{(m-1)}) + U'(\theta^{(m-1)})\bigl(\theta^{(m)} - \theta^{(m-1)}\bigr) \]
Impondo \(U(\theta^{(m)}) = 0\) e isolando \(\theta^{(m)}\) chegamos a (4.7).
Intuição: Se a log-verossimilhança tem um único máximo, o algoritmo “desce” (ou “sobe”) pela tangente da função escore em direção ao zero. Quanto mais próximo do máximo, mais preciso se torna cada passo — a convergência costuma ser quadrática perto da solução.
3.2 Método de Scoring
Uma alternativa é substituir \(U'(\theta)\) por seu valor esperado \(E[U'(\theta)] = -I(\theta)\), onde \(I\) é a informação de Fisher. Isso dá o método de scoring:
\[ \theta^{(m)} = \theta^{(m-1)} + \frac{U(\theta^{(m-1)})}{I(\theta^{(m-1)})} \tag{4.11} \]
Para a distribuição de Weibull com \(\lambda = 2\), usando as expressões de \(b(\theta)\) e \(c(\theta)\) de (4.2) na fórmula da informação (equação 3.15 do Capítulo 3):
\[ I(\theta) = \frac{\lambda^2 N}{\theta^2} = \frac{4N}{\theta^2} \tag{4.10} \]
Vantagem do scoring: \(E[U']\) é tipicamente mais estável numericamente que \(U'\) em si — especialmente longe do máximo. Além disso, leva naturalmente à matriz de informação de Fisher para MLGs multidimensionais, que tem forma analítica conveniente.
3.3 Iterações para os Vasos de Pressão
Partindo de \(\theta^{(1)} = \bar{y} = 8805{,}7\) (média amostral como chute inicial):
lambda <- 2; N <- length(y); sum_y2 <- sum(y^2)
U_func <- function(th) -2*N/th + 2*sum_y2/th^3
Up_func <- function(th) 2*N/th^2 - 6*sum_y2/th^4
I_func <- function(th) lambda^2*N/th^2
theta_vals <- numeric(5)
theta_vals[1] <- mean(y)
for(m in 2:5) theta_vals[m] <- theta_vals[m-1] + U_func(theta_vals[m-1])/I_func(theta_vals[m-1])
iter_df <- data.frame(
Iteração = 1:4,
`θ(m)` = round(theta_vals[1:4], 1),
`U × 10⁶` = round(U_func(theta_vals[1:4]) * 1e6, 2),
`U′ × 10⁶` = round(Up_func(theta_vals[1:4]) * 1e6, 2),
`E(U′) × 10⁶` = round(-I_func(theta_vals[1:4]) * 1e6, 2),
`U/U′` = round(-U_func(theta_vals[1:4])/Up_func(theta_vals[1:4]), 3),
`U/E(U′)` = round( U_func(theta_vals[1:4])/I_func(theta_vals[1:4]), 3)
)
kable(iter_df,
caption = "Tabela 4.2 — Iterações de Newton–Raphson / Scoring para o parâmetro de escala θ da distribuição Weibull (λ = 2).",
digits = 3, align = "r",
col.names = c("Iteração", "θ(m)", "U×10⁶", "U′×10⁶", "E(U′)×10⁶", "U/U′", "U/E(U′)")) |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13) |>
column_spec(2, bold = TRUE)| Iteração | θ(m) | U×10⁶ | U′×10⁶ | E(U′)×10⁶ | U/U′ | U/E(U′) |
|---|---|---|---|---|---|---|
| 1 | 8805.7 | 2915.75 | -3.52 | -2.53 | 828.084 | 1153.510 |
| 2 | 9959.2 | -132.01 | -1.94 | -1.98 | -68.174 | -66.802 |
| 3 | 9892.4 | -0.45 | -2.00 | -2.00 | -0.226 | -0.226 |
| 4 | 9892.2 | 0.00 | -2.00 | -2.00 | 0.000 | 0.000 |
Leitura da tabela:
- A coluna \(\theta^{(m)}\) mostra a convergência rápida para \(\hat\theta \approx 9892\).
- \(U \times 10^6\) se aproxima de zero — confirmando que chegamos ao máximo.
- \(U'\) e \(E(U')\) têm valores muito similares em todas as iterações — justificando o uso do método de scoring em vez de Newton–Raphson puro.
- As colunas \(U/U'\) e \(U/E(U')\) (os passos de atualização) também são muito próximas.
Estimativa final e erro padrão:
Sendo \(\hat\theta = 9892{,}2\) e a informação de Fisher avaliada no máximo \(I = \lambda^2 N / \hat\theta^2 = 4 \times 49 / 9892{,}2^2 \approx 2{,}00 \times 10^{-6}\):
\[ \widehat{\text{ep}}(\hat\theta) = \frac{1}{\sqrt{I}} = \frac{1}{\sqrt{2{,}00 \times 10^{-6}}} \approx 707 \]
Intervalo de confiança de 95% aproximado (assumindo normalidade assintótica do EMV):
\[ \hat\theta \pm 1{,}96 \times \widehat{\text{ep}}(\hat\theta) = 9892 \pm 1{,}96 \times 707 \approx (8506,\; 11278) \]
Interpretação: A curvatura da log-verossimilhança no máximo determina a precisão da estimativa. Uma log-verossimilhança muito “plana” (curvatura pequena, \(I\) pequeno) indica grande incerteza — e portanto erro padrão grande. Uma curva “pontiaguda” (alta curvatura, \(I\) grande) indica alta precisão.
4 Teoria Geral de Estimação para MLGs
4.1 Configuração
Considere \(Y_1, \ldots, Y_N\) independentes satisfazendo as propriedades de um MLG. A log-verossimilhança total é:
\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^N \ell_i = \sum_{i=1}^N \bigl[y_i b(\theta_i) + c(\theta_i) + d(y_i)\bigr] \tag{4.13} \]
onde \(\mu_i = E(Y_i) = -c'(\theta_i)/b'(\theta_i)\) e \(g(\mu_i) = \mathbf{x}_i^T\boldsymbol{\beta} = \eta_i\).
4.2 Derivação do Escore para \(\beta_j\)
Para obter o EMV de \(\beta_j\) precisamos de \(\partial\ell/\partial\beta_j = 0\). Pela regra da cadeia:
\[ U_j = \frac{\partial\ell}{\partial\beta_j} = \sum_{i=1}^N \frac{\partial\ell_i}{\partial\theta_i} \cdot \frac{\partial\theta_i}{\partial\mu_i} \cdot \frac{\partial\mu_i}{\partial\beta_j} \tag{4.17} \]
Calculando cada fator separadamente:
Fator 1: \(\dfrac{\partial\ell_i}{\partial\theta_i} = y_i b'(\theta_i) + c'(\theta_i) = b'(\theta_i)(y_i - \mu_i)\)
(usamos que \(c'(\theta_i) = -b'(\theta_i)\mu_i\) da Eq. 4.14)
Fator 2: \(\dfrac{\partial\theta_i}{\partial\mu_i} = \dfrac{1}{\partial\mu_i/\partial\theta_i} = \dfrac{1}{b'(\theta_i)\,\text{var}(Y_i)}\)
(pois \(\partial\mu_i/\partial\theta_i = b'(\theta_i)\,\text{var}(Y_i)\), que vem da diferenciação de \(\mu_i = -c'/b'\))
Fator 3: \(\dfrac{\partial\mu_i}{\partial\beta_j} = \dfrac{\partial\mu_i}{\partial\eta_i} \cdot \dfrac{\partial\eta_i}{\partial\beta_j} = \dfrac{\partial\mu_i}{\partial\eta_i} \cdot x_{ij}\)
Combinando os três fatores:
\[ \boxed{U_j = \sum_{i=1}^N \frac{(y_i - \mu_i)}{\text{var}(Y_i)} \cdot \frac{\partial\mu_i}{\partial\eta_i} \cdot x_{ij}} \tag{4.18} \]
Leitura da fórmula: O escore \(U_j\) é uma soma ponderada de resíduos \((y_i - \mu_i)\). O peso de cada observação é \(\frac{1}{\text{var}(Y_i)} \cdot \frac{\partial\mu_i}{\partial\eta_i}\) — inversamente proporcional à variância (observações mais precisas têm mais peso) e proporcional à sensibilidade da média ao preditor linear.
4.3 Matriz de Informação \(\mathbf{I}\)
O elemento \((j,k)\) da matriz de informação de Fisher é \(I_{jk} = E[U_j U_k]\). Como os \(Y_i\) são independentes, os termos cruzados se anulam e:
\[ I_{jk} = \sum_{i=1}^N \frac{x_{ij}\,x_{ik}}{\text{var}(Y_i)} \left(\frac{\partial\mu_i}{\partial\eta_i}\right)^2 \tag{4.20} \]
Definindo os pesos:
\[ w_{ii} = \frac{1}{\text{var}(Y_i)} \left(\frac{\partial\mu_i}{\partial\eta_i}\right)^2 \tag{4.23} \]
e a matriz diagonal \(\mathbf{W} = \text{diag}(w_{11}, \ldots, w_{NN})\), a matriz de informação toma a forma matricial compacta:
\[ \mathbf{I} = \mathbf{X}^T \mathbf{W} \mathbf{X} \]
Por que os pesos \(w_{ii}\) têm esta forma? - O fator \(1/\text{var}(Y_i)\): observações com menor variância carregam mais informação. - O fator \((\partial\mu_i/\partial\eta_i)^2\): mede o quanto a média responde ao preditor linear. Se a ligação é a identidade, \(\partial\mu/\partial\eta = 1\) e todos os pesos são \(1/\text{var}(Y_i)\). Para a ligação log, \(\partial\mu/\partial\eta = \mu_i\), ponderando mais observações onde \(\mu_i\) é grande.
4.4 Mínimos Quadrados Iterativamente Reponderados (IRLS)
A atualização do método de scoring para o vetor \(\boldsymbol{\beta}\) é:
\[ \mathbf{b}^{(m)} = \mathbf{b}^{(m-1)} + \mathbf{I}^{-1}_{(m-1)}\,\mathbf{U}_{(m-1)} \tag{4.21} \]
Pré-multiplicando por \(\mathbf{I}_{(m-1)}\) e usando \(\mathbf{I} = \mathbf{X}^T\mathbf{W}\mathbf{X}\):
\[ \underbrace{(\mathbf{X}^T\mathbf{W}\mathbf{X})}_{\mathbf{I}}\mathbf{b}^{(m)} = \underbrace{(\mathbf{X}^T\mathbf{W}\mathbf{X})\mathbf{b}^{(m-1)} + \mathbf{U}_{(m-1)}}_{\mathbf{X}^T\mathbf{W}\mathbf{z}} \tag{4.25} \]
onde o vetor de respostas ajustadas \(\mathbf{z}\) tem elementos:
\[ z_i = \sum_{k=1}^p x_{ik}\,b_k^{(m-1)} + (y_i - \mu_i)\frac{\partial\eta_i}{\partial\mu_i} = \hat\eta_i + (y_i - \mu_i)\frac{\partial\eta_i}{\partial\mu_i} \tag{4.24} \]
Interpretação do vetor \(\mathbf{z}\): Em cada iteração, \(z_i\) é uma “resposta linearizada” — uma versão de \(y_i\) que incorpora a não-linearidade da função de ligação por meio do termo de correção \((y_i - \mu_i)\partial\eta_i/\partial\mu_i\). À medida que o algoritmo converge (\(y_i \approx \mu_i\)), \(z_i \to \hat\eta_i\).
4.4.1 Resumo do Algoritmo IRLS
ALGORITMO IRLS PARA MLGs
════════════════════════════════════════════════
Entrada: X, y, distribuição, função de ligação g
1. Inicializar b(0) (ex.: média da resposta)
2. REPITA até convergência:
a. Calcular μᵢ = g⁻¹(xᵢᵀ b)
b. Calcular pesos: wᵢ = (∂μᵢ/∂ηᵢ)² / var(Yᵢ)
c. Calcular respostas ajustadas: zᵢ = η̂ᵢ + (yᵢ - μᵢ)(∂ηᵢ/∂μᵢ)
d. Resolver: (XᵀWX) b(novo) = XᵀWz
e. b ← b(novo)
3. SAÍDA: β̂ = b(converge)
Var(β̂) ≈ (XᵀWX)⁻¹ = I⁻¹
════════════════════════════════════════════════5 Exemplo Completo: Regressão de Poisson com Ligação Identidade
5.1 Dados e Modelo
df_p <- data.frame(
i = 1:9,
y = c(2, 3, 6, 7, 8, 9, 10, 12, 15),
x = c(-1, -1, 0, 0, 0, 0, 1, 1, 1)
)
kable(t(df_p[,2:3]),
caption = "Tabela 4.3 — Dados para a regressão de Poisson.",
col.names = paste0("Obs. ", 1:9),
align = "r") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13)| Obs. 1 | Obs. 2 | Obs. 3 | Obs. 4 | Obs. 5 | Obs. 6 | Obs. 7 | Obs. 8 | Obs. 9 | |
|---|---|---|---|---|---|---|---|---|---|
| y | 2 | 3 | 6 | 7 | 8 | 9 | 10 | 12 | 15 |
| x | -1 | -1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 |
Modelo especificado:
\[ E(Y_i) = \mu_i = \beta_1 + \beta_2 x_i = \mathbf{x}_i^T\boldsymbol{\beta}; \qquad Y_i \sim \text{Poisson}(\mu_i) \]
Função de ligação: identidade, \(g(\mu_i) = \mu_i\).
Por que a ligação identidade aqui? No gráfico de dispersão (ver abaixo), a relação entre \(y\) e \(x\) parece linear e positiva, sugerindo que a ligação identidade (em vez da canônica log) é adequada. Além disso, a variabilidade aumenta com \(y\) — consistente com \(\text{var}(Y_i) = \mu_i\) (Poisson).
ggplot(df_p, aes(x = x, y = y)) +
geom_point(size = 4, color = "#4C72B0", alpha = 0.85) +
geom_smooth(method = "lm", se = FALSE, color = "#C44E52",
linetype = "dashed", linewidth = 0.9) +
scale_x_continuous(breaks = c(-1, 0, 1)) +
labs(x = "x", y = "y (contagem)",
title = "Regressão de Poisson — Ligação Identidade",
subtitle = "Reta de regressão linear sobreposta") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))`geom_smooth()` using formula = 'y ~ x'
5.2 Simplificação com Ligação Identidade
Para a ligação identidade \(g(\mu) = \mu\):
\[ \frac{\partial\mu_i}{\partial\eta_i} = 1 \quad \Rightarrow \quad w_{ii} = \frac{1}{\text{var}(Y_i)} = \frac{1}{\mu_i} = \frac{1}{\beta_1 + \beta_2 x_i} \]
A resposta ajustada simplifica para:
\[ z_i = \hat\eta_i + (y_i - \mu_i)\cdot 1 = \mu_i + y_i - \mu_i = y_i \]
Portanto \(\mathbf{z} = \mathbf{y}\) em cada iteração — os dados originais são sempre a “resposta ajustada”.
5.3 Iterações IRLS
Partindo de \(\mathbf{b}^{(1)} = (7,\; 5)^T\) (lidos aproximadamente do gráfico):
y_p <- df_p$y; x_p <- df_p$x; N_p <- length(y_p)
X <- cbind(1, x_p)
# IRLS manual
b <- c(7, 5)
iters <- list()
for(m in 1:4) {
mu <- X %*% b
W <- diag(as.vector(1/mu))
XtWX <- t(X) %*% W %*% X
XtWz <- t(X) %*% W %*% y_p # z = y sempre
b <- solve(XtWX, XtWz)
iters[[m]] <- round(b, 5)
}
iter_df2 <- data.frame(
m = 1:4,
b1_m = c(7, iters[[1]][1], iters[[2]][1], iters[[3]][1]),
b2_m = c(5, iters[[1]][2], iters[[2]][2], iters[[3]][2])
)
kable(iter_df2,
caption = "Tabela 4.4 — Aproximações sucessivas para os coeficientes da regressão de Poisson.",
col.names = c("m", "b₁(m)", "b₂(m)"),
digits = 5, align = "r") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13) |>
column_spec(2:3, bold = TRUE)| m | b₁(m) | b₂(m) |
|---|---|---|
| 1 | 7.00000 | 5.00000 |
| 2 | 7.45139 | 4.93750 |
| 3 | 7.45163 | 4.93531 |
| 4 | 7.45163 | 4.93530 |
O algoritmo converge já na 3ª iteração: \(\hat\beta_1 = 7{,}4516\) e \(\hat\beta_2 = 4{,}9353\).
5.4 Erros Padrão e Intervalo de Confiança
A matriz de variância-covariância estimada de \(\hat{\boldsymbol{\beta}}\) é \(\mathbf{I}^{-1} = (\mathbf{X}^T\mathbf{W}\mathbf{X})^{-1}\):
b_final <- iters[[3]]
mu_f <- X %*% b_final
W_f <- diag(as.vector(1/mu_f))
I_mat <- t(X) %*% W_f %*% X
Iinv <- solve(I_mat)
vcov_df <- data.frame(
` ` = c("β₁", "β₂"),
b1 = round(Iinv[,1], 4),
b2 = round(Iinv[,2], 4)
)
kable(vcov_df,
caption = "Matriz de variância-covariância estimada I⁻¹.",
col.names = c("", "β₁", "β₂"), align = "r") |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 13)| β₁ | β₂ | ||
|---|---|---|---|
| β₁ | 0.7817 | 0.4166 | |
| x_p | β₂ | 0.4166 | 1.1863 |
Erros padrão:
\[ \widehat{\text{ep}}(\hat\beta_1) = \sqrt{0{,}7817} \approx 0{,}884, \qquad \widehat{\text{ep}}(\hat\beta_2) = \sqrt{1{,}1863} \approx 1{,}089 \]
Intervalo de 95% para \(\beta_2\) (inclinação):
\[ \hat\beta_2 \pm 1{,}96 \times \widehat{\text{ep}}(\hat\beta_2) = 4{,}9353 \pm 1{,}96 \times 1{,}0892 \approx (2{,}80,\; 7{,}07) \]
df_p$mu_hat <- as.vector(X %*% b_final)
ggplot(df_p, aes(x = x)) +
geom_segment(aes(xend = x, y = y, yend = mu_hat),
color = "gray60", linetype = "dotted") +
geom_point(aes(y = y), size = 4, color = "#4C72B0",
shape = 16, alpha = 0.9) +
geom_point(aes(y = mu_hat), size = 4, color = "#C44E52",
shape = 17, alpha = 0.9) +
scale_x_continuous(breaks = c(-1, 0, 1)) +
labs(x = "x", y = "y",
title = "Observados (●) vs. Ajustados (▲)",
subtitle = expression(paste(hat(beta)[1], " = 7.45, ", hat(beta)[2], " = 4.94"))) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
5.5 Código R e Stata
O modelo acima pode ser ajustado diretamente com glm() no R:
library(dobson)
data(poisson)
res <- glm(y ~ x, family = poisson(link = "identity"), data = poisson)
summary(res)Ou em Stata:
.glm y x, family(poisson) link(identity) vce(eim)Ambos produzem \(\hat\beta_1 = 7{,}4516\,(0{,}8841)\) e \(\hat\beta_2 = 4{,}9353\,(1{,}0892)\).
6 Síntese: Pesos e Respostas Ajustadas por Distribuição
A tabela abaixo mostra como os elementos do IRLS se especializam para as três distribuições principais, com as duas funções de ligação mais comuns:
| Distribuição | Ligação | var(Y) | ∂μ/∂η | w | z |
|---|---|---|---|---|---|
| Normal | Identidade | σ² | 1 | 1/σ² | η+(y−μ) |
| Normal | Log | σ² | μ | μ²/σ² | η+(y−μ)/μ |
| Poisson | Identidade | μ | 1 | 1/μ | y |
| Poisson | Log | μ | μ | μ | η+(y−μ)/μ |
| Binomial | Logit | nπ(1−π) | nπ(1−π) | nπ(1−π) | η+(y−μ)/[nπ(1−π)] |
| Binomial | Log | nπ(1−π) | nπ | nπ/(1−π) | η+(y−μ)/(nπ) |
7 Demonstrações Didáticas
7.1 Derivação do Fator 2 da Regra da Cadeia
Precisamos calcular \(\partial\theta_i / \partial\mu_i\). Sabemos que \(\mu_i = -c'(\theta_i)/b'(\theta_i)\) (Eq. 4.14). Diferenciando \(\mu_i\) em relação a \(\theta_i\):
\[ \frac{\partial\mu_i}{\partial\theta_i} = \frac{-c''(\theta_i)b'(\theta_i) + c'(\theta_i)b''(\theta_i)}{[b'(\theta_i)]^2} \]
Usando a fórmula de \(\text{var}(Y_i)\) do Capítulo 3 (Eq. 3.12):
\[ \text{var}(Y_i) = \frac{b''(\theta_i)c'(\theta_i) - c''(\theta_i)b'(\theta_i)}{[b'(\theta_i)]^3} = \frac{\partial\mu_i/\partial\theta_i}{b'(\theta_i)/b'(\theta_i)} \cdot \frac{1}{b'(\theta_i)} = \frac{\partial\mu_i}{\partial\theta_i}\cdot\frac{1}{b'(\theta_i)} \]
Portanto \(\partial\mu_i/\partial\theta_i = b'(\theta_i)\,\text{var}(Y_i)\), e:
\[ \frac{\partial\theta_i}{\partial\mu_i} = \frac{1}{b'(\theta_i)\,\text{var}(Y_i)} \]
7.2 Simplificação do Escore (Equação 4.18)
Combinando os três fatores da regra da cadeia (Eq. 4.17):
\[ U_j = \sum_{i=1}^N \underbrace{b'(\theta_i)(y_i - \mu_i)}_{\partial\ell_i/\partial\theta_i} \cdot \underbrace{\frac{1}{b'(\theta_i)\,\text{var}(Y_i)}}_{\partial\theta_i/\partial\mu_i} \cdot \underbrace{\frac{\partial\mu_i}{\partial\eta_i}\,x_{ij}}_{\partial\mu_i/\partial\beta_j} \]
Os fatores \(b'(\theta_i)\) cancelam:
\[ U_j = \sum_{i=1}^N \frac{(y_i - \mu_i)}{\text{var}(Y_i)}\cdot\frac{\partial\mu_i}{\partial\eta_i}\cdot x_{ij} \quad \checkmark \]
7.3 Derivação de \(\mathbf{X}^T\mathbf{W}\mathbf{z}\) a partir de \(\mathbf{I}\mathbf{b}^{(m-1)} + \mathbf{U}\)
O lado direito de (4.22) é \(\mathbf{I}\mathbf{b}^{(m-1)} + \mathbf{U}\). O \(j\)-ésimo elemento é:
\[ \sum_{k=1}^p I_{jk}\,b_k^{(m-1)} + U_j = \sum_{i=1}^N w_{ii}\,x_{ij} \underbrace{\sum_{k=1}^p x_{ik}\,b_k^{(m-1)}}_{\hat\eta_i} + \sum_{i=1}^N \frac{(y_i-\mu_i)}{\text{var}(Y_i)}\frac{\partial\mu_i}{\partial\eta_i} x_{ij} \]
\[ = \sum_{i=1}^N w_{ii}\,x_{ij}\left[\hat\eta_i + \frac{(y_i-\mu_i)}{\text{var}(Y_i)}\frac{\partial\mu_i}{\partial\eta_i} \cdot \frac{1}{w_{ii}}\right] \]
Usando \(w_{ii} = (\partial\mu_i/\partial\eta_i)^2/\text{var}(Y_i)\), o termo entre colchetes é exatamente \(z_i\) (Eq. 4.24):
\[ = \sum_{i=1}^N w_{ii}\,x_{ij}\,z_i = (\mathbf{X}^T\mathbf{W}\mathbf{z})_j \quad \checkmark \]
8 Resumo do Capítulo 4
ESTIMAÇÃO EM MLGs — FLUXO GERAL
═══════════════════════════════════════════════════════════
PONTO DE PARTIDA: ℓ(β) = Σ[yᵢb(θᵢ) + c(θᵢ) + d(yᵢ)]
↓ regra da cadeia
ESCORE: Uⱼ = Σᵢ [(yᵢ−μᵢ)/var(Yᵢ)] × (∂μᵢ/∂ηᵢ) × xᵢⱼ
↓ valor esperado
INFORMAÇÃO: I = XᵀWX, com wᵢ = (∂μᵢ/∂ηᵢ)²/var(Yᵢ)
↓ método de scoring iterativo
IRLS: (XᵀWX) b(m) = XᵀWz
zᵢ = η̂ᵢ + (yᵢ−μᵢ)(∂ηᵢ/∂μᵢ)
↓ convergência
ESTIMATIVA: β̂, Var(β̂) ≈ I⁻¹ = (XᵀWX)⁻¹
ep(β̂ⱼ) = √[I⁻¹]ⱼⱼ
═══════════════════════════════════════════════════════════| Conceito | Fórmula / Resultado |
|---|---|
| Algoritmo de Newton–Raphson | \(\theta^{(m)} = \theta^{(m-1)} - U/U'\) |
| Método de scoring | \(\theta^{(m)} = \theta^{(m-1)} + U/I\) |
| Escore para \(\beta_j\) | \(U_j = \sum_i [(y_i-\mu_i)/\text{var}(Y_i)](\partial\mu_i/\partial\eta_i)x_{ij}\) |
| Pesos IRLS | \(w_i = (\partial\mu_i/\partial\eta_i)^2 / \text{var}(Y_i)\) |
| Resposta ajustada | \(z_i = \hat\eta_i + (y_i - \mu_i)(\partial\eta_i/\partial\mu_i)\) |
| Equação IRLS | \((\mathbf{X}^T\mathbf{W}\mathbf{X})\mathbf{b} = \mathbf{X}^T\mathbf{W}\mathbf{z}\) |
| Variância assintótica | \(\widehat{\text{Var}}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^T\mathbf{W}\mathbf{X})^{-1}\) |
| Erro padrão | \(\widehat{\text{ep}}(\hat\beta_j) = \sqrt{[(\mathbf{X}^T\mathbf{W}\mathbf{X})^{-1}]_{jj}}\) |
Próximo capítulo: O Capítulo 5 usa as estimativas e seus erros padrão para construir testes de hipóteses e intervalos de confiança para parâmetros de MLGs, introduzindo as estatísticas de desvio (deviance) e de Wald.
9 Exercícios
(4.1 - Dobson): Os dados abaixo mostram o número de casos de AIDS na Austrália por trimestre, de 1984 a 1988. Nesta fase inicial da epidemia, os casos pareciam crescer exponencialmente.
(a) Plotar \(y_i\) contra o período \(i = 1, \ldots, 20\).
(b) Verificar o modelo \(\log\lambda_i = \theta\log i\) plotando \(\log y_i\) vs. \(\log i\).
(c) Ajustar o MLG de Poisson com ligação log, \(\log\lambda_i = \beta_1 + \beta_2 x_i\) (com \(x_i = \log i\)), pelos primeiros princípios (IRLS manual).
(d) Repetir com glm() e comparar.
(4.2 - Dobson): Os dados abaixo são tempos até a morte (em semanas) após o diagnóstico e o \(\log_{10}\) da contagem inicial de células brancas (WBC) para 17 pacientes com leucemia.
(a) Plotar \(y_i\) vs. \(x_i\) — há tendência?
(b) Identificar a função de ligação adequada para \(E(Y_i) = \exp(\beta_1 + \beta_2 x_i)\).
(c) Mostrar que \(E(Y) = 1/\theta\) e \(\text{var}(Y) = 1/\theta^2\) para a distribuição Exponencial.
(d) Ajustar o modelo com glm().
(e) Comparar observados e ajustados; analisar resíduos padronizados.
(4.3 - Dobson): Seja \(Y_1, \ldots, Y_N\) uma amostra aleatória de \(N({\log\beta},\; \sigma^2)\), com \(\sigma^2\) conhecido.
(a) Encontrar o EMV de \(\beta\) pelos primeiros princípios.
(b) Verificar as Equações (4.18) e (4.25) para este caso.
- Referências
Dobson, A. J. & Barnett, A. G. (2018). An Introduction to Generalized Linear Models (4ª ed.). CRC Press / Chapman & Hall.
Andrews, D. F. & Herzberg, A. M. (1985). Data: A Collection of Problems from Many Fields for the Student and Research Worker. Springer.