library(tidyverse)
library(patchwork)
library(viridis)
library(GGally)
theme_set(theme_bw(base_size = 14))
set.seed(2025)Modelos Lineares Bayesianos
Regressão Múltipla, ANOVA e ANCOVA — Capítulo 7 (Kinas & Andrade)
1 Regressão Linear Múltipla
1.1 Formulação do Modelo
O modelo de regressão linear simples é generalizado de forma direta quando há \(k-1\) variáveis preditoras quantitativas \(x_j\), \(j = 1, \ldots, k-1\). O modelo para a \(i\)-ésima observação é:
\[ y_i = \beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i} + \cdots + \beta_{k-1} x_{k-1,i} + w_i \]
com \(w_i \overset{\text{perm}}{\sim} N(0, \sigma)\), isto é, os erros são permutáveis com distribuição normal, média zero e desvio padrão constante \(\sigma\) (homocedasticidade).
Na notação matricial compacta, com \(\mathbf{X}\) sendo a matriz de planejamento \(n \times k\) (incluindo a coluna de uns para o intercepto) e \(\boldsymbol{\beta} = (\beta_0, \beta_1, \ldots, \beta_{k-1})^\top\):
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{w}, \quad \mathbf{w} \sim N_n(\mathbf{0},\, \sigma^2 \mathbf{I}) \]
1.2 Distribuições Posteriores (Priori de Jeffreys)
Adotando a priori não-informativa de Jeffreys \(p(\beta_0, \ldots, \beta_{k-1}, \sigma^2) \propto 1/\sigma^2\) e com os dados \(D = \{(\mathbf{x}_i, y_i) : i = 1, \ldots, n\}\), as distribuições posteriores marginais são:
A posterior conjunta não tem forma fechada conveniente para simulação direta, mas a estrutura hierárquica \(\boldsymbol{\beta} \mid \sigma^2, D\) é gaussiana, o que sugere o seguinte algoritmo de simulação:
- Simular \(\tau^{(l)} \sim \text{Gamma}\!\left(\dfrac{n-k}{2},\; \dfrac{n-k}{2}S_e^2\right)\) para \(l = 1, \ldots, m\);
- Calcular \(\sigma^{2(l)} = 1/\tau^{(l)}\);
- Simular \(\boldsymbol{\beta}^{(l)} \sim N_k\!\left(\hat{\boldsymbol{\beta}},\; \sigma^{2(l)} V_\beta\right)\).
1.3 Padronização das Covariáveis
Antes de ajustar o modelo, é conveniente padronizar as variáveis preditoras:
\[ z_j = \frac{x_j - \bar{x}_j}{S_{x_j}} \]
Essa transformação produz dois benefícios importantes:
- Remove a correlação dos coeficientes, pois centraliza cada variável em torno de zero;
- Torna os coeficientes comparáveis: cada \(\beta_j\) expressa o efeito de um incremento de um desvio padrão em \(x_j\) sobre \(y\).
1.4 Exemplo: Captura de Tubarão-Martelo
Os dados de Kotas et al. (2005) contêm capturas de tubarão-martelo (Sphyrna spp.) em espinhel pelágico. Modelamos \(y = \log(\text{capkg})\) em função das preditoras padronizadas \(z_1 = [\log(\text{temp}) - \overline{\log(\text{temp})}]/S_{\log(\text{temp})}\) e \(z_2 = [\log(\text{nanz}) -
\overline{\log(\text{nanz})}]/S_{\log(\text{nanz})}\), onde capkg representa a captura em kg, temp a temperatura superficial do mar, e nanz é o número de anzóis (esforço de pesca).
## Dados do tubarao-martelo (Kotas et al., 2005)
if(!file.exists("martelo.csv")) {
download.file(
"http://leg.ufpr.br/~fernandomayer/data/martelo.csv", "martelo.csv")
}
da <- read_csv("martelo.csv")da |>
select(temp, nanz, capkg) |>
ggpairs(
lower = list(continuous = wrap("points", alpha = 0.6,
color = "steelblue4")),
diag = list(continuous = wrap("densityDiag", fill = "steelblue",
alpha = 0.4))) +
theme(strip.text = element_text(size = 11))
ggplot(da, aes(x = capkg)) +
geom_histogram(bins = 15, fill = "steelblue", alpha = 0.6,
color = "steelblue4") +
labs(x = "Captura (kg)", y = "Frequencia")
## Transformacoes e padronizacao
y <- log(da$capkg)
x1 <- log(da$temp)
x2 <- log(da$nanz)
z1 <- (x1 - mean(x1)) / sd(x1)
z2 <- (x2 - mean(x2)) / sd(x2)ggpairs(tibble(y, x1, x2),
lower = list(continuous = wrap("points", alpha = 0.6,
color = "steelblue4")),
diag = list(continuous = wrap("densityDiag", fill = "steelblue",
alpha = 0.4))) +
theme(strip.text = element_text(size = 11))
## Ajuste por minimos quadrados (estrutura para a posterior)
mod_lm <- lm(y ~ z1 + z2)
se <- summary(mod_lm)$sigma
est <- coef(mod_lm)
Vbeta <- summary(mod_lm)$cov.unscaled
n <- nrow(da)
## Simulacao da distribuicao posterior conjunta (m = 5000)
m <- 5000
prec <- rgamma(m, (n - 3) / 2, se^2 * (n - 3) / 2)
sig2 <- 1 / prec
beta_post <- t(sapply(seq_len(m), function(l)
MASS::mvrnorm(1, est, sig2[l] * Vbeta)))
colnames(beta_post) <- c("beta0", "beta1", "beta2")
bind_cols(beta_post, sigma2 = sig2) |>
as_tibble()# A tibble: 5,000 × 4
beta0 beta1 beta2 sigma2
<dbl> <dbl> <dbl> <dbl>
1 6.13 -0.218 1.28 0.471
2 6.09 -0.284 1.02 0.528
3 6.25 -0.0982 1.19 0.458
4 5.99 -0.0965 1.40 0.418
5 6.12 -0.0587 1.12 0.494
6 6.21 -0.371 1.45 0.550
7 6.27 -0.239 1.17 0.832
8 6.08 -0.296 1.22 0.571
9 5.92 -0.314 1.25 0.695
10 6.04 -0.283 1.28 0.786
# ℹ 4,990 more rows## Sumario posterior
sumario <- apply(
cbind(beta_post, sigma = sqrt(sig2)), 2,
function(x) c(
"2.5%" = quantile(x, 0.025),
"Mediana" = median(x),
"97.5%" = quantile(x, 0.975),
"Media" = mean(x),
"dp" = sd(x)
)
)
knitr::kable(
round(t(sumario), 3),
caption = paste(
"Sumario das distribuicoes posteriores marginais —",
"regressao multipla (dados simulados, estrutura de",
"Kotas et al., 2005)."
)
)| 2.5%.2.5% | Mediana | 97.5%.97.5% | Media | dp | |
|---|---|---|---|---|---|
| beta0 | 5.936 | 6.136 | 6.334 | 6.136 | 0.102 |
| beta1 | -0.432 | -0.228 | -0.022 | -0.228 | 0.103 |
| beta2 | 1.030 | 1.228 | 1.424 | 1.228 | 0.101 |
| sigma | 0.604 | 0.727 | 0.898 | 0.733 | 0.075 |
df_post <- as.data.frame(beta_post) |>
mutate(sigma = sqrt(sig2))
nomes <- c(
beta0 = expression(beta[0]),
beta1 = expression(beta[1]),
beta2 = expression(beta[2]),
sigma = expression(sigma)
)
plist <- lapply(names(nomes), function(par) {
ggplot(df_post, aes(x = .data[[par]])) +
geom_density(fill = "steelblue", alpha = 0.4,
color = "steelblue4") +
geom_vline(xintercept = 0, linetype = "dashed",
color = "tomato") +
labs(x = nomes[[par]], y = "Densidade") +
theme(axis.title.x = element_text(size = 13))
})
wrap_plots(plist, nrow = 1)
Interpretação: Com preditoras padronizadas, o impacto de um desvio padrão adicional em \(z_2\) (número de anzóis) sobre o logaritmo da captura é aproximadamente 5 vezes maior que o impacto equivalente de \(z_1\) (temperatura), pois \(\hat{b}_2 / |\hat{b}_1| \approx 5.4\).
2 Análise de Variância (ANOVA)
2.1 Formulação do Modelo
Quando todas as variáveis preditoras são categóricas, o modelo linear assume a estrutura:
\[ y_{g,i} = \mu_g + w_{g,i}, \quad g = 1, \ldots, G,\quad i = 1, \ldots, n_g \]
onde \(\mu_g\) é a média do \(g\)-ésimo nível (efeito fixo) e \(w_{g,i} \overset{\text{perm}}{\sim} N(0, \sigma)\). A hipótese de homocedasticidade (variância igual entre grupos) é central.
Para o modelo de um fator fixo, as estatísticas suficientes são:
\[ \bar{y}_g = \frac{1}{n_g}\sum_{i=1}^{n_g} y_{g,i},\quad SQ_g = \sum_{i=1}^{n_g}(y_{g,i} - \bar{y}_g)^2,\quad n = \sum_{g=1}^G n_g,\quad S_e^2 = \frac{\sum_{g=1}^G SQ_g}{n - G} \]
2.2 Distribuições Posteriores (Priori de Jeffreys)
O algoritmo de simulação é análogo ao da regressão:
- Simular \(\tau^{(l)} \sim \text{Gamma}\!\left(\dfrac{n-G}{2},\; \dfrac{n-G}{2}S_e^2\right)\);
- Calcular \(\sigma^{2(l)} = 1/\tau^{(l)}\);
- Simular \((\mu_1^{(l)}, \ldots, \mu_G^{(l)}) \sim N_G\!\left((\bar{y}_1, \ldots, \bar{y}_G),\; \sigma^{2(l)} V_\mu\right)\), com \(V_\mu = \text{diag}(1/n_1, \ldots, 1/n_G)\).
2.3 Conexão com a ANOVA Frequentista
Na abordagem frequentista, o foco é testar \(H_0: \mu_1 = \cdots = \mu_G\). A estatística de teste é a razão entre variâncias:
\[ F = \frac{\text{QMTr}}{\text{QMR}} = \frac{SQTr/(G-1)}{SQR/(n-G)} \overset{H_0}{\sim} F(G-1,\; n-G) \]
onde \(SQTr = SQT - SQR = \sum_{g,i}(y_{g,i} - \bar{\bar{y}})^2 - \sum_g SQ_g\).
Na abordagem bayesiana, o foco recai nas distribuições posteriores de cada \(\mu_g\) e nas probabilidades de interesse sobre contrastes, como \(P(\mu_a - \mu_b > 0 \mid D)\) — sem a rigidez de uma decisão binária.
2.4 Exemplo: CPUE do Tubarão-Martelo por Área de Pesca
O objetivo é comparar as médias de CPUE entre 3 áreas de pesca. A Captura por Unidade de Esforço (CPUE) é obtida pela razão entre a captura e uma medida padrão de esforço de pesca, e pode ser considerada como proporcional à abundância de uma população.
## Dados simulados com estrutura do Exemplo 7.3 (Kotas et al., 2005)
## Tres areas de pesca (9, 14, 15); resposta: log(CPUE)
## Cria area como fator
da$areaf <- as.factor(da$area)
## Calcula CPUE
da$cpue <- (da$capkg * 1000)/da$nanz
da |>
split(~ area) |>
map(\(x) summary(x$cpue))$`9`
Min. 1st Qu. Median Mean 3rd Qu. Max.
13.82 42.85 58.94 70.04 75.54 206.25
$`14`
Min. 1st Qu. Median Mean 3rd Qu. Max.
35.33 114.48 142.86 199.44 164.92 756.25
$`15`
Min. 1st Qu. Median Mean 3rd Qu. Max.
15.00 43.79 57.12 87.97 115.27 264.21 ## Calcula log(CPUE)
da$lcpue <- log(da$cpue)
p1 <- ggplot(da, aes(x = cpue)) +
geom_histogram(
bins = 15, fill = "steelblue", alpha = 0.6, color = "steelblue4") +
labs(x = "CPUE", y = "Frequência")
p2 <- ggplot(da, aes(x = lcpue)) +
geom_histogram(
bins = 15, fill = "steelblue", alpha = 0.6, color = "steelblue4") +
labs(x = "log(CPUE)", y = "Frequência")
p1 + p2
## Verifica distribuicao
p1 <- ggplot(da, aes(x = areaf, y = cpue, fill = areaf)) +
geom_boxplot(alpha = 0.6, show.legend = FALSE) +
scale_fill_viridis_d() +
labs(x = "Área", y = "CPUE")
p2 <- ggplot(da, aes(x = areaf, y = lcpue, fill = areaf)) +
geom_boxplot(alpha = 0.6, show.legend = FALSE) +
scale_fill_viridis_d() +
labs(x = "Área", y = "log(CPUE)")
p1 + p2
## Estatisticas suficientes
fnsq <- function(x) sum((x - mean(x))^2)
(ygbar <- with(da, tapply(lcpue, areaf, mean))) 9 14 15
4.085268 5.002832 4.215782 (sqg <- with(da, tapply(lcpue, area, fnsq))) 9 14 15
6.733349 7.782122 9.138331 (ng_obs <- with(da, tapply(lcpue, area, length))) 9 14 15
20 15 18 (nG <- sum(ng_obs))[1] 53(G <- length(ng_obs))[1] 3(se2G <- sum(sqg) / (nG - G))[1] 0.473076## Simulacao da posterior conjunta (m = 3000)
m <- 3000
(Vmu <- diag(1/ng_obs)) [,1] [,2] [,3]
[1,] 0.05 0.00000000 0.00000000
[2,] 0.00 0.06666667 0.00000000
[3,] 0.00 0.00000000 0.05555556tauG <- rgamma(m, (nG - G) / 2, se2G * (nG - G) / 2)
varG <- 1 / tauG
muGpost <- t(sapply(seq_len(m), function(l)
MASS::mvrnorm(1, as.numeric(ygbar), varG[l] * Vmu)))
colnames(muGpost) <- paste0("mu", names(ygbar))
muGpost <- as_tibble(muGpost)
## Posterior na escala original da CPUE
emuGpost <- exp(muGpost)
## Explora posterior
icr <- function(x, cr = 0.95) {
li <- (1 - cr)/2
ls <- cr + li
quant <- quantile(x, probs = c(li, ls, 0.5))
c(quant, media = mean(x), dp = sd(x))
}
apply(emuGpost, 2, icr) mu9 mu14 mu15
2.5% 43.601472 104.48008 48.01703
97.5% 81.323077 211.93654 92.67513
50% 59.594448 148.55403 67.43059
media 60.330391 150.84134 68.15799
dp 9.584031 27.47907 11.48647## Gráficos das posteriores
mu_post <- emuGpost |>
pivot_longer(everything(), names_to = "area",
values_to = "cpue") |>
mutate(area = factor(gsub("mu", "", area), levels = c(9, 14, 15)))
ggplot(mu_post, aes(x = cpue, color = area, fill = area)) +
geom_density(alpha = 0.3) +
scale_color_viridis_d() +
scale_fill_viridis_d() +
labs(x = expression(mu[CPUE]), y = "Densidade",
title = "Posteriores de CPUE",
color = "Area", fill = "Area") +
xlim(0, 400)
## Contrastes
mean(emuGpost$mu14 - emuGpost$mu9 > 0)[1] 0.9996667mean(emuGpost$mu14 - emuGpost$mu15 > 0)[1] 0.9996667## Probabilidade P(mu15 - mu9 > 0 | dados)
(prob15_9 <- mean(emuGpost$mu15 - emuGpost$mu9 > 0))[1] 0.705## Diferenca posterior
dif15_9 <- emuGpost$mu15 - emuGpost$mu9
mean(dif15_9 > 0)[1] 0.705mean(dif15_9)[1] 7.8276icr(dif15_9) 2.5% 97.5% 50% media dp
-20.922151 38.237868 7.830061 7.827600 14.905048 mu_dif <- tibble(dif = dif15_9)
ggplot(mu_dif, aes(x = dif)) +
geom_histogram(aes(y = after_stat(density)), bins = 40,
fill = "steelblue", alpha = 0.7,
color = "white") +
geom_density(color = "steelblue4", linewidth = 1) +
geom_vline(xintercept = 0, linetype = "dashed",
color = "tomato") +
labs(x = expression(mu[15]^{cpue} - mu[9]^{cpue}),
y = "Densidade", title = "Diferenca area 15 vs 9")
Leitura dos resultados: A área 14 destaca-se claramente das demais, com CPUE estocasticamente superior. A diferença entre as áreas 9 e 15 é mais sutil: \(P(\mu_{15} - \mu_9 > 0 \mid D) \approx 0.73\), mas o intervalo de credibilidade para essa diferença inclui zero, evidenciando incerteza substancial.
Outra forma de analisar as diferenças é por meio das distribuições acumuladas das posteriores.
p1 <- mu_post |>
group_by(area) |>
arrange(cpue, .by_group = TRUE) |>
mutate(cdf = seq(0, 1, length.out = n())) |>
ggplot(aes(x = cpue, y = cdf, color = area)) +
geom_line(linewidth = 1) +
scale_color_viridis_d() +
xlim(0, 300) +
labs(x = expression(mu[CPUE]), y = "Probabilidade acumulada",
color = "Área")
p2 <- mu_post |>
filter(area %in% c(9, 15)) |>
group_by(area) |>
arrange(cpue, .by_group = TRUE) |>
mutate(cdf = seq(0, 1, length.out = n())) |>
ggplot(aes(x = cpue, y = cdf, color = area)) +
geom_line(linewidth = 1) +
scale_color_viridis_d() +
xlim(0, 150) +
labs(x = expression(mu[CPUE]), y = "Probabilidade acumulada",
color = "Área")
p1 + p2
Como se compara à tradicional?
mod <- lm(lcpue ~ factor(area), data = da)
summary(mod)
Call:
lm(formula = lcpue ~ factor(area), data = da)
Residuals:
Min 1Q Median 3Q Max
-1.50773 -0.36092 -0.05852 0.17746 1.62554
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.0853 0.1538 26.563 < 2e-16 ***
factor(area)14 0.9176 0.2349 3.906 0.000282 ***
factor(area)15 0.1305 0.2235 0.584 0.561811
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6878 on 50 degrees of freedom
Multiple R-squared: 0.2536, Adjusted R-squared: 0.2238
F-statistic: 8.494 on 2 and 50 DF, p-value: 0.0006671ygbar 9 14 15
4.085268 5.002832 4.215782 coef(mod)[1] + coef(mod)[2](Intercept)
5.002832 coef(mod)[1] + coef(mod)[3](Intercept)
4.215782 anova(mod)Analysis of Variance Table
Response: lcpue
Df Sum Sq Mean Sq F value Pr(>F)
factor(area) 2 8.037 4.0185 8.4944 0.0006671 ***
Residuals 50 23.654 0.4731
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## SQ Total
sum(anova(mod)$`Sum Sq`)[1] 31.69079sum((da$lcpue - mean(da$lcpue))^2)[1] 31.69079## SQ Residual
sum(sqg)[1] 23.6538## SQ Tratamento
sum(anova(mod)$`Sum Sq`) - sum(sqg)[1] 8.036984sum((da$lcpue - mean(da$lcpue))^2) - sum(sqg)[1] 8.036984## QM Residuo
se2G[1] 0.4730763 Exercícios
Os exercícios abaixo são selecionados da lista do Capítulo 7 de Kinas & Andrade para complementar os tópicos abordados.
3.1 Exercício 7.3 — Relação Peso-Comprimento da Cabrinha
Os dados abaixo correspondem ao comprimento (cm) e peso (kg) de 25 exemplares de cabrinha (Prionotus punctatus):
comp_p <- c(10.2, 19.0, 22.3, 23.8, 24.9, 26.1, 26.9, 27.9,
28.5, 29.2, 30.1, 31.0, 31.7, 32.4, 32.9, 33.3,
34.1, 34.6, 35.5, 36.3, 37.0, 38.0, 39.0, 39.6, 43.0)
peso_p <- c(0.01, 0.09, 0.12, 0.18, 0.19, 0.23, 0.24, 0.31,
0.33, 0.32, 0.33, 0.40, 0.45, 0.49, 0.67, 0.48,
0.54, 0.63, 0.59, 0.65, 0.50, 0.81, 0.87, 0.71, 0.91)
prionotus <- data.frame(comp = comp_p, peso = peso_p)Pretende-se ajustar o modelo alométrico \(y = \beta_0 x^{\beta_1}\), onde \(y\) é o peso e \(x\) é o comprimento, e que na forma logaritmizada fica:
\[ \log(y) = \beta_0' + \beta_1' \log(x) + v,\quad v \sim N(0, \sigma) \]
com \(\beta_0' = \log(\beta_0)\) e \(\beta_1' = \beta_1\).
a) Construa diagramas de dispersão de \(y\) vs. \(x\) e de \(\log(y)\) vs. \(\log(x)\). Justifique a adoção do modelo linear na escala logarítmica.
b) Usando a priori não-informativa de Jeffreys, obtenha uma amostra de tamanho \(m = 3000\) da distribuição posterior conjunta de \((\beta_0', \beta_1', \sigma)\). Construa diagramas de dispersão de \(\beta_0'\) vs. \(\beta_1'\) e de \(\beta_0\) vs. \(\beta_1\) (escala original).
c) Represente as distribuições posteriores marginais de \(\beta_0\) e \(\beta_1\) com histogramas e calcule os ICr95%. Interprete o ICr95% para \(\beta_1\), considerando que peixes isométricos têm \(\beta_1 \approx 3\).
d) Construa a distribuição preditiva para o peso de um indivíduo com \(x_p = 20\) cm. Apresente o ICr99% e avalie a precisão da predição.
3.2 Exercício 7.4 — Regressão Múltipla para Fósforo no Solo
Os dados abaixo (Snedecor & Cochran, 1989, p. 335) descrevem fósforo inorgânico (\(X_1\), ppm), orgânico (\(X_2\), ppm) e disponível para as plantas (\(Y\), ppm) em amostras de solo:
X1 <- c(0.4, 0.4, 3.1, 0.6, 4.7, 1.7, 9.4, 10.1, 11.6,
12.6, 10.9, 23.1, 23.1, 21.6, 23.1, 1.9, 29.9)
X2 <- c(53, 23, 19, 34, 24, 65, 44, 31, 29,
58, 37, 46, 50, 44, 56, 36, 51)
Y <- c(64, 60, 71, 61, 54, 77, 81, 93, 93,
51, 76, 96, 77, 93, 95, 54, 99)
fosforo <- data.frame(Y = Y, X1 = X1, X2 = X2)a) Construa diagramas de dispersão de \(Y\) vs. \(X_1\) e \(Y\) vs. \(X_2\).
b) Ajuste um modelo bayesiano de regressão múltipla com variáveis explanatórias padronizadas e priori de Jeffreys. Obtenha \(m = 1000\) amostras da distribuição posterior.
c) Represente graficamente e construa um sumário das distribuições posteriores marginais de todos os parâmetros.
d) Obtenha 1000 amostras da distribuição preditiva para \(X_1 = 15\,\text{ppm}\) e \(X_2 = 40\,\text{ppm}\).
e) Apresente o sumário e o gráfico da distribuição preditiva posterior.
3.3 Exercício 7.5 — ANOVA para Comprimento de Tartarugas
Em um estudo sobre quatro praias de desova, pesquisadores mediram o comprimento curvilíneo da carapaça (CC, em cm):
praia <- c(1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4)
CC <- c(100.4, 96.7, 96.3, 96.2,
94.6, 92.8, 89.0,
92.2, 96.6, 95.3, 95.9, 93.2, 93.2,
95.3, 97.3, 99.3, 99.0, 96.4)
tartarugas <- data.frame(praia = factor(praia), CC = CC)Considere o modelo de fatores fixos: \(y_{g,i} = \mu_g + w_{g,i}\), com \(w_{g,i} \sim N(0, \sigma)\).
a) Calcule médias e variâncias por praia; represente com boxplots.
b) Usando a priori de Jeffreys, retire \(m = 3000\) amostras da distribuição posterior conjunta de \((\mu_1, \mu_2, \mu_3, \mu_4, \sigma)\).
c) Calcule sumários e represente as distribuições posteriores marginais de \(\mu_g\), \(g = 1, \ldots, 4\).
d) Calcule \(P(\mu_3 > \mu_2 \mid D)\) e outros contrastes que julgar útil.
e) Faça os gráficos das posteriores acumuladas.
4 Gabarito dos Exercícios
4.1 Gabarito — Exercício 7.3
p1 <- ggplot(prionotus, aes(x = comp, y = peso)) +
geom_point(color = "steelblue", size = 2.5) +
labs(x = "Comprimento (cm)", y = "Peso (kg)",
title = "Escala original")
p2 <- ggplot(prionotus, aes(x = log(comp), y = log(peso))) +
geom_point(color = "tomato", size = 2.5) +
geom_smooth(method = "lm", se = TRUE, color = "tomato4") +
labs(x = "log(Comprimento)", y = "log(Peso)",
title = "Escala logaritmica")
p1 + p2
Na escala original a relação é claramente não-linear e a variância aumenta com \(x\). Na escala logarítmica a relação é bem aproximada por uma reta e a dispersão é mais homogênea — o que justifica o modelo linear em \(\log(y)\) vs. \(\log(x)\).
## b) Amostra posterior conjunta
y_p <- log(prionotus$peso)
x_p <- log(prionotus$comp)
n_p <- length(y_p)
prio_lm <- lm(y_p ~ x_p)
se_p <- summary(prio_lm)$sigma
est_p <- coef(prio_lm)
Vb_p <- summary(prio_lm)$cov.unscaled
m_p <- 3000
prec_p <- rgamma(m_p, (n_p - 2) / 2, se_p^2 * (n_p - 2) / 2)
sig2_p <- 1 / prec_p
beta_p <- t(sapply(seq_len(m_p), function(l)
MASS::mvrnorm(1, est_p, sig2_p[l] * Vb_p)))
colnames(beta_p) <- c("b0l", "b1l")
## Conversao para escala original: beta0 = exp(b0'), beta1 = b1'
beta0_orig <- exp(beta_p[, "b0l"])
beta1_orig <- beta_p[, "b1l"]
df_b <- tibble(
b0l = beta_p[, "b0l"],
b1l = beta_p[, "b1l"],
beta0 = beta0_orig,
beta1 = beta1_orig
)p1 <- ggplot(df_b, aes(x = b0l, y = b1l)) +
geom_point(alpha = 0.1, color = "steelblue", size = 0.5) +
stat_density_2d(color = "steelblue4", linewidth = 0.5) +
labs(x = expression(beta[0]^"'"),
y = expression(beta[1]^"'"),
title = "Escala logaritmica")
p2 <- ggplot(df_b, aes(x = beta0, y = beta1)) +
geom_point(alpha = 0.1, color = "tomato", size = 0.5) +
stat_density_2d(color = "tomato4", linewidth = 0.5) +
labs(x = expression(beta[0]), y = expression(beta[1]),
title = "Escala original")
p1 + p2
## c) Marginais e ICr95%
icr_b0 <- quantile(beta0_orig, c(0.025, 0.975))
icr_b1 <- quantile(beta1_orig, c(0.025, 0.975))
round(icr_b0, 5) 2.5% 97.5%
0e+00 1e-05 round(icr_b1, 3) 2.5% 97.5%
3.002 3.390 p1 <- ggplot(df_b, aes(x = beta0)) +
geom_histogram(aes(y = after_stat(density)), bins = 40,
fill = "steelblue", alpha = 0.7,
color = "white") +
geom_vline(xintercept = icr_b0, linetype = "dashed",
color = "tomato") +
labs(x = expression(beta[0]), y = "Densidade",
title = expression(beta[0]))
p2 <- ggplot(df_b, aes(x = beta1)) +
geom_histogram(aes(y = after_stat(density)), bins = 40,
fill = "tomato", alpha = 0.7,
color = "white") +
geom_vline(xintercept = icr_b1, linetype = "dashed",
color = "steelblue") +
geom_vline(xintercept = 3, linetype = "solid",
color = "black", linewidth = 1) +
labs(x = expression(beta[1]), y = "Densidade",
title = expression(beta[1]))
p1 + p2
Interpretação de \(\beta_1\): O ICr95% para \(\beta_1\) é 3.0022 a 3.39. O valor de referência \(\beta_1 = 3\) (crescimento isométrico) está fora do ICr95%, sugerindo que o crescimento da espécie não é isométrico (peso aumenta de forma linear com o comprimento), o que pode ser visto pelo gráfico da relação peso x comprimento acima.
Neste caso, como o valor de referência e o limite inferior do intervalo de credibilidade são peraticamente iguais, outras medidas podem ser calculadas a partir da posterior. Por exemplo, podemos calcular \(P(\beta_1 < 3)\), que pode ajudar a chegar à uma conclusão:
mean(beta1_orig < 3)[1] 0.02366667## e) Distribuicao preditiva para xp = 20 cm
xp <- log(20)
yhat_p <- beta_p[, "b0l"] + beta_p[, "b1l"] * xp
ytil_log <- rnorm(m_p, yhat_p, sqrt(sig2_p))
ytil_orig <- exp(ytil_log)
icr99 <- quantile(ytil_orig, c(0.005, 0.995))
round(icr99, 3) 0.5% 99.5%
0.063 0.148 ggplot(data.frame(ytil = ytil_orig), aes(x = ytil)) +
geom_density(fill = "steelblue", alpha = 0.5,
color = "steelblue4") +
geom_vline(xintercept = icr99, linetype = "dashed",
color = "tomato") +
xlim(0, 0.5) +
labs(x = "Peso predito (kg)", y = "Densidade",
title = expression(
paste("Distribuicao preditiva para ",
x[p], " = 20 cm")
))
O ICr99% para o peso de um indivíduo com 20 cm é [0.063, 0.148] kg. A amplitude evidencia incerteza preditiva considerável, reflexo tanto da incerteza paramétrica quanto da variabilidade individual \(\sigma\).
4.2 Gabarito — Exercício 7.4
p1 <- ggplot(fosforo, aes(x = X1, y = Y)) +
geom_point(color = "steelblue", size = 2.5) +
geom_smooth(method = "lm", se = FALSE, color = "steelblue4") +
labs(x = expression(X[1]~(ppm)), y = "Y (ppm)",
title = expression(Y~vs.~X[1]))
p2 <- ggplot(fosforo, aes(x = X2, y = Y)) +
geom_point(color = "tomato", size = 2.5) +
geom_smooth(method = "lm", se = FALSE, color = "tomato4") +
labs(x = expression(X[2]~(ppm)), y = "Y (ppm)",
title = expression(Y~vs.~X[2]))
p1 + p2
Ambas as variáveis apresentam associação positiva com \(Y\), com \(X_1\) mostrando relação mais forte e \(X_2\) relação mais moderada.
## b) Regressao multipla bayesiana com preditoras padronizadas
z1_f <- (X1 - mean(X1)) / sd(X1)
z2_f <- (X2 - mean(X2)) / sd(X2)
n_f <- length(Y)
fosforo_lm <- lm(Y ~ z1_f + z2_f)
se_f <- summary(fosforo_lm)$sigma
est_f <- coef(fosforo_lm)
Vb_f <- summary(fosforo_lm)$cov.unscaled
m_f <- 1000
prec_f <- rgamma(m_f, (n_f - 3) / 2, se_f^2 * (n_f - 3) / 2)
sig2_f <- 1 / prec_f
beta_f <- t(sapply(seq_len(m_f), function(l)
MASS::mvrnorm(1, est_f, sig2_f[l] * Vb_f)))
colnames(beta_f) <- c("beta0", "beta1", "beta2")## c) Sumario e graficos
sumario_f <- apply(
cbind(beta_f, sigma = sqrt(sig2_f)), 2,
function(x) c(
"2.5%" = quantile(x, 0.025),
"Mediana" = median(x),
"97.5%" = quantile(x, 0.975),
"Media" = mean(x),
"dp" = sd(x)
)
)
knitr::kable(
round(t(sumario_f), 2),
caption = paste("Sumario posterior — regressao multipla",
"para fosforo disponivel.")
)| 2.5%.2.5% | Mediana | 97.5%.97.5% | Media | dp | |
|---|---|---|---|---|---|
| beta0 | 70.03 | 76.05 | 82.79 | 76.08 | 3.17 |
| beta1 | 5.29 | 12.66 | 19.34 | 12.56 | 3.55 |
| beta2 | -8.08 | -1.15 | 4.90 | -1.36 | 3.46 |
| sigma | 8.94 | 12.63 | 19.09 | 13.02 | 2.62 |
Exercicio 7.4c: posteriores marginais — regressao multipla para fosforo.
df_f <- as.data.frame(cbind(beta_f, sigma = sqrt(sig2_f)))
nomes_f <- c(
beta0 = expression(beta[0]),
beta1 = expression(beta[1]),
beta2 = expression(beta[2]),
sigma = expression(sigma)
)
plist_f <- lapply(names(nomes_f), function(par) {
ggplot(df_f, aes(x = .data[[par]])) +
geom_density(fill = "steelblue", alpha = 0.4,
color = "steelblue4") +
geom_vline(xintercept = 0, linetype = "dashed",
color = "tomato") +
labs(x = nomes_f[[par]], y = "Densidade")
})
wrap_plots(plist_f, nrow = 1)
## d) e e) Distribuicao preditiva para X1=15, X2=40
x1_new <- (15 - mean(X1)) / sd(X1)
x2_new <- (40 - mean(X2)) / sd(X2)
yhat_f <- beta_f[, "beta0"] + beta_f[, "beta1"] * x1_new +
beta_f[, "beta2"] * x2_new
ytil_f <- rnorm(m_f, yhat_f, sqrt(sig2_f))
icr95_f <- quantile(ytil_f, c(0.025, 0.975))
round(icr95_f, 1) 2.5% 97.5%
53.3 110.0 ggplot(data.frame(ytil = ytil_f), aes(x = ytil)) +
geom_density(fill = "steelblue", alpha = 0.5,
color = "steelblue4") +
geom_vline(xintercept = icr95_f, linetype = "dashed",
color = "tomato") +
labs(x = "Y predito (ppm)", y = "Densidade",
title = expression(
paste("Preditiva para ", X[1], "=15, ", X[2], "=40")
))
4.3 Gabarito — Exercício 7.5
## a) Sumario e boxplots
knitr::kable(
tartarugas |>
group_by(praia) |>
summarise(
n = n(),
media = mean(CC),
variancia = var(CC)
) |>
rename(Praia = praia,
`Media (cm)` = media,
`Variancia` = variancia) |>
mutate(across(where(is.numeric), \(x) round(x, 2))),
caption = "Estatisticas descritivas do CC por praia."
)| Praia | n | Media (cm) | Variancia |
|---|---|---|---|
| 1 | 4 | 97.40 | 4.05 |
| 2 | 3 | 92.13 | 8.17 |
| 3 | 6 | 94.40 | 3.12 |
| 4 | 5 | 97.46 | 2.89 |
Exercicio 7.5a: comprimento da carapaca de tartarugas por praia.
ggplot(tartarugas, aes(x = praia, y = CC, fill = praia)) +
geom_boxplot(alpha = 0.6, show.legend = FALSE) +
scale_fill_viridis_d(option = "C", end = 0.85) +
labs(x = "Praia", y = "CC (cm)",
title = "Comprimento da carapaca por praia")
## b) Posterior via simulacao (Jeffreys)
CC_vec <- tartarugas$CC
area_v <- tartarugas$praia
G_t <- nlevels(area_v)
n_t <- length(CC_vec)
ygbar_t <- tapply(CC_vec, area_v, mean)
sqg_t <- tapply(CC_vec, area_v,
function(x) sum((x - mean(x))^2))
ng_t <- tapply(CC_vec, area_v, length)
se2_t <- sum(sqg_t) / (n_t - G_t)
m_t <- 3000
Vmu_t <- diag(1 / as.numeric(ng_t))
tau_t <- rgamma(m_t, (n_t - G_t) / 2,
se2_t * (n_t - G_t) / 2)
var_t <- 1 / tau_t
muTpost <- t(sapply(seq_len(m_t), function(l)
MASS::mvrnorm(1, as.numeric(ygbar_t), var_t[l] * Vmu_t)))
colnames(muTpost) <- paste0("mu", 1:4)## c) Sumario e graficos
sumario_t <- apply(
cbind(muTpost, sigma = sqrt(var_t)), 2,
function(x) c(
"2.5%" = quantile(x, 0.025),
"Mediana" = median(x),
"97.5%" = quantile(x, 0.975),
"Media" = mean(x),
"dp" = sd(x)
)
)
knitr::kable(
round(t(sumario_t), 2),
caption = "Sumario posterior — ANOVA para comprimento de tartarugas."
)| 2.5%.2.5% | Mediana | 97.5%.97.5% | Media | dp | |
|---|---|---|---|---|---|
| mu1 | 95.28 | 97.41 | 99.67 | 97.42 | 1.09 |
| mu2 | 89.60 | 92.15 | 94.56 | 92.14 | 1.24 |
| mu3 | 92.70 | 94.40 | 96.16 | 94.39 | 0.88 |
| mu4 | 95.56 | 97.49 | 99.52 | 97.48 | 0.99 |
| sigma | 1.44 | 2.04 | 3.15 | 2.10 | 0.44 |
Exercicio 7.5c: posteriores marginais dos efeitos de praia.
df_mu_t <- as.data.frame(muTpost) |>
pivot_longer(everything(), names_to = "praia",
values_to = "mu") |>
mutate(praia = gsub("mu", "Praia ", praia))
ggplot(df_mu_t, aes(x = mu, color = praia, fill = praia)) +
geom_density(alpha = 0.3, linewidth = 0.9) +
scale_color_viridis_d(option = "C", end = 0.85) +
scale_fill_viridis_d(option = "C", end = 0.85) +
labs(x = expression(mu[g]~(cm)), y = "Densidade",
color = "Praia", fill = "Praia",
title = "Posteriores marginais dos efeitos de praia")
## d) P(mu3 > mu2 | D)
muTpost <- as_tibble(muTpost)
(prob_3_2 <- mean(muTpost$mu3 > muTpost$mu2))[1] 0.9326667\(P(\mu_3 > \mu_2 \mid D)\) = 0.933. A probabilidade é relativamente alta, indicando que tartarugas da praia 3 tendem a ser maiores que as da praia 2.
## CDF
mu_post <- muTpost |>
pivot_longer(everything(), names_to = "praia",
values_to = "CC") |>
mutate(praia = factor(gsub("mu", "", praia)))
mu_post |>
group_by(praia) |>
arrange(CC, .by_group = TRUE) |>
mutate(cdf = seq(0, 1, length.out = n())) |>
ggplot(aes(x = CC, y = cdf, color = praia)) +
geom_line(linewidth = 1) +
scale_color_viridis_d() +
labs(x = expression(mu[CC]), y = "Probabilidade acumulada",
color = "Praia",
title = "CDF posterior de CC por praia")