library(tidyverse)
library(rjags)
library(coda)
library(patchwork)
set.seed(2026)
theme_set(theme_bw(base_size = 14))Análise Bayesiana de Decisão
1 Apresentação
2 Parte I – O exemplo: maturação sexual do peixe-galo
2.1 Contexto biológico e a decisão de manejo
O peixe-galo é um recurso pesqueiro de interesse econômico capturado predominantemente com redes de emalhe. O tamanho de malha define quais peixes a rede retém:
- Malhas pequenas selecionam indivíduos antes da maturidade reprodutiva (prejudicando o estoque)
- Malhas grandes deixam escapar peixes que poderiam ser comercializados
Estabelece-se como objetivo de manejo capturar preferencialmente fêmeas que já atingiram maturidade sexual, isto é, peixes com comprimento médio superior ao LT-50, o comprimento em que cerca de 50% das fêmeas estão sexualmente maturas.
A decisão a tomar é qual rede usar entre três alternativas com tamanhos crescentes de malha:
- \(\delta_1 \equiv\) rede A
- \(\delta_2 \equiv\) rede B
- \(\delta_3 \equiv\) rede C
Como o LT-50 é desconhecido, a escolha precisa levar em conta a incerteza sobre seu valor.
2.2 Os dados
Os dados disponíveis correspondem a 17 fêmeas amostradas em quatro classes de comprimento (cm) na costa sul do Brasil:
x <- c(12.5, 22.5, 32.5, 62.5)
n <- c(3, 5, 4, 5)
y <- c(0, 1, 3, 5)
galo <- tibble(
classe = c("[10, 15)", "[20, 25)", "[30, 35)", "[60, 65)"),
ponto_medio = x,
total = n,
maturas = y,
proporcao = y / n
)
galo# A tibble: 4 × 5
classe ponto_medio total maturas proporcao
<chr> <dbl> <dbl> <dbl> <dbl>
1 [10, 15) 12.5 3 0 0
2 [20, 25) 22.5 5 1 0.2
3 [30, 35) 32.5 4 3 0.75
4 [60, 65) 62.5 5 5 1 galo |>
ggplot(aes(ponto_medio, proporcao)) +
geom_point(size = 3) +
scale_x_continuous(limits = c(10, 70)) +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "Comprimento (cm)",
y = "Proporção de fêmeas maturas")
2.3 O modelo logístico
Sendo \(p_i\) a probabilidade de uma fêmea da classe \(i\) estar matura, assume-se que
\[ y_i \mid p_i \sim \text{Binomial}(n_i, p_i), \qquad i = 1, \ldots, k, \]
com a estrutura logística (centrada em \(\bar{x} = 32{,}5\) cm para facilitar a interpretação):
\[ \text{logit}(p_i) = \log\!\left(\frac{p_i}{1 - p_i}\right) = \beta_0 + \beta_1 (x_i - \bar{x}). \]
Aqui \(\beta_0\) é o logito da probabilidade de maturidade no comprimento médio e \(\beta_1\) é o incremento esperado no logito por cm adicional. O parâmetro de interesse para o manejo é
\[ \text{LT-50} = -\frac{\beta_0}{\beta_1} + \bar{x}. \]
2.3.1 Análise convencional via glm()
A abordagem convencional para esses dados consiste em ajustar uma regressão logística por máxima verossimilhança via glm().
## Centraliza
xbar <- mean(x)
xc <- x - xbar
## Ajuste do modelo
mglm <- glm(cbind(y, n - y) ~ xc, family = binomial(link = "logit"))
summary(mglm)
Call:
glm(formula = cbind(y, n - y) ~ xc, family = binomial(link = "logit"))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1630 1.1404 1.020 0.3078
xc 0.2669 0.1458 1.831 0.0671 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 14.00544 on 3 degrees of freedom
Residual deviance: 0.10686 on 2 degrees of freedom
AIC: 7.6181
Number of Fisher Scoring iterations: 6A estimativa pontual do LT-50 segue diretamente da definição \(\widehat{LT\text{-}50} = -\hat{\beta}_0/\hat{\beta}_1 + \bar{x}\):
beta_hat <- coef(mglm)
(LT50_hat <- unname(-beta_hat[1]/beta_hat[2] + xbar))[1] 28.14272A saída do glm() fornece estimativas pontuais e erros padrão para \(\beta_0\) e \(\beta_1\) individualmente, mas não entrega diretamente uma medida de incerteza para o LT-50, que é uma função não-linear desses parâmetros.
Diferentemente da abordagem bayesiana, em que a posterior conjunta de \((\beta_0, \beta_1)\) se propaga naturalmente para qualquer função dela, inclusive o LT-50, na abordagem convencional o LT-50 sai como estimativa pontual apenas.
Para construir um intervalo de confiança assintótico para o LT-50, pode-se recorrer ao método delta.
NotaMétodo delta para a variância do LT-50
Seja \(g(\boldsymbol{\beta}) = -\beta_0/\beta_1 + \bar{x}\). Pelas propriedades assintóticas dos estimadores de máxima verossimilhança,
\[ \hat{\boldsymbol{\beta}} \;\sim\; N(\boldsymbol{\beta},\, V) \quad \text{(assintoticamente)}, \qquad V = \text{Cov}(\hat{\boldsymbol{\beta}}), \]
com \(V\) estimada por vcov(mglm). Uma expansão de Taylor de primeira ordem de \(g\) em torno de \(\hat{\boldsymbol{\beta}}\) leva a
\[ \text{Var}\!\left(g(\hat{\boldsymbol{\beta}})\right) \;\approx\; \nabla g(\hat{\boldsymbol{\beta}})^{\!\top} \, V \, \nabla g(\hat{\boldsymbol{\beta}}), \]
com gradiente
\[ \nabla g(\boldsymbol{\beta}) = \left(\dfrac{\partial g}{\partial \beta_0},\, \dfrac{\partial g}{\partial \beta_1}\right)^{\!\top} = \left(-\dfrac{1}{\beta_1},\, \dfrac{\beta_0}{\beta_1^{2}}\right)^{\!\top}. \]
Daí
\[ \text{Var}(\widehat{LT\text{-}50}) \;\approx\; \frac{V_{00}}{\beta_1^{2}} + \frac{\beta_0^{2}\, V_{11}}{\beta_1^{4}} - \frac{2\,\beta_0\, V_{01}}{\beta_1^{3}}, \]
de onde o erro padrão é \(\text{EP} = \sqrt{\widehat{\text{Var}}}\) e o intervalo de confiança assintótico de 95% é \(\widehat{LT\text{-}50} \pm 1{,}96 \cdot \text{EP}\).
Aplicando o método delta:
## Matriz V
(V_beta <- vcov(mglm)) (Intercept) xc
(Intercept) 1.3006233 0.11662573
xc 0.1166257 0.02125125## Gradiente
(grad_g <- c(-1 / beta_hat[2], beta_hat[1] / beta_hat[2]^2)) xc (Intercept)
-3.746618 16.325046 ## Variancia assintótica
(var_LT50 <- as.numeric(t(grad_g) %*% V_beta %*% grad_g))[1] 9.654122(ep_LT50 <- sqrt(var_LT50))[1] 3.107108LT50_freq <- tibble(
estimativa = LT50_hat,
EP = ep_LT50,
IC_inf = LT50_hat - 1.96 * ep_LT50,
IC_sup = LT50_hat + 1.96 * ep_LT50
)
LT50_freq# A tibble: 1 × 4
estimativa EP IC_inf IC_sup
<dbl> <dbl> <dbl> <dbl>
1 28.1 3.11 22.1 34.22.3.2 Especificação bayesiana
A função de verossimilhança conjunta dos \(k = 4\) pontos amostrais é
\[ L(\beta_0, \beta_1 \mid \mathbf{y}) \propto p_i^{\,y_i} (1 - p_i)^{\,n_i - y_i}, \quad p_i = \frac{e^{\beta_0 + \beta_1(x_i - \bar{x})}} {1 + e^{\beta_0 + \beta_1(x_i - \bar{x})}}. \]
Adotam-se prioris independentes e pouco informativas (em função da precisão \(\tau\))
\[ \beta_0 \sim \text{N}(0, \tau = 0{,}001), \qquad \beta_1 \sim \text{N}(0, \tau = 0{,}001). \]
Pelo Teorema de Bayes, a posterior conjunta é
\[ p(\beta_0, \beta_1 \mid \mathbf{y}) \;\propto\; L(\beta_0, \beta_1 \mid \mathbf{y})\, p(\beta_0)\, p(\beta_1). \]
Esta posterior não tem forma analítica fechada, o que motiva o uso de métodos MCMC para se obter uma amostra dela.
2.3.3 Implementação em JAGS
O modelo é escrito na linguagem do JAGS e ajustado via pacote rjags. Note que o LT-50 não entra no bloco do modelo: por ser função determinística de \((\beta_0, \beta_1)\), é mais seguro calculá-lo depouis, a partir das amostras posteriores. Dessa forma, evita-se divisão por zero quando a cadeia visita valores de \(\beta_1\) muito próximos de zero.
## Modelo
mod_string <- "model {
for (i in 1:N) {
logit(p[i]) <- beta0 + beta1 * xc[i]
y[i] ~ dbin(p[i], n[i])
}
beta0 ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
## Dados
xbar <- mean(x)
xc <- x - xbar
dados_jags <- list(
y = y,
n = n,
xc = xc,
N = length(y)
)
## Calibra modelo
modelo <- jags.model(
file = textConnection(mod_string),
data = dados_jags,
inits = list(beta0 = 1, beta1 = 0.5),
n.chains = 1,
n.adapt = 1000,
quiet = TRUE
)
## Burn-in
update(modelo, n.iter = 5000, progress.bar = "none")
## Amostra final
amostra <- coda.samples(
model = modelo,
variable.names = c("beta0", "beta1"),
n.iter = 45000,
thin = 9, # 5000 amostras apos thinning
progress.bar = "none"
)A amostra é convertida em tibble e o LT-50 é calculado como \(-\beta_0/\beta_1 + \bar{x}\) para cada par simulado:
posterior <- as_tibble(as.matrix(amostra)) |>
mutate(LT50 = -beta0 / beta1 + xbar)
posterior# A tibble: 5,000 × 3
beta0 beta1 LT50
<dbl> <dbl> <dbl>
1 0.257 0.0790 29.2
2 2.58 0.395 26.0
3 1.16 0.197 26.6
4 3.49 0.511 25.7
5 2.57 0.624 28.4
6 1.27 0.506 30.0
7 2.60 0.369 25.4
8 2.89 0.398 25.2
9 1.94 0.596 29.2
10 1.34 0.186 25.3
# ℹ 4,990 more rowsPosterior conjunta de \((\beta_0, \beta_1)\)
posterior |>
ggplot(aes(beta0, beta1)) +
geom_point(alpha = 0.2, size = 0.7) +
geom_density_2d(color = "firebrick") +
geom_vline(xintercept = 0, lty = "dashed", color = "tomato") +
geom_hline(yintercept = 0, lty = "dashed", color = "tomato") +
labs(x = expression(beta[0]), y = expression(beta[1]))
Marginais e o LT-50
posterior |>
pivot_longer(everything(),
names_to = "parametro",
values_to = "valor") |>
group_by(parametro) |>
summarise(
media = mean(valor),
dp = sd(valor),
q025 = quantile(valor, 0.025),
q500 = quantile(valor, 0.5),
q975 = quantile(valor, 0.975)
)# A tibble: 3 × 6
parametro media dp q025 q500 q975
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 LT50 28.2 3.22 22.6 28.0 34.9
2 beta0 1.69 1.35 -0.490 1.54 4.74
3 beta1 0.375 0.172 0.114 0.353 0.768p0 <- posterior |>
ggplot(aes(beta0)) +
geom_density(fill = "steelblue", alpha = 0.5) +
geom_vline(xintercept = 0, lty = "dashed", col = "tomato") +
labs(x = expression(beta[0]), y = "Densidade")
p1 <- posterior |>
ggplot(aes(beta1)) +
geom_density(fill = "steelblue", alpha = 0.5) +
geom_vline(xintercept = 0, lty = "dashed", col = "tomato") +
labs(x = expression(beta[1]), y = "Densidade")
p2 <- posterior |>
ggplot(aes(LT50)) +
geom_density(fill = "darkorange", alpha = 0.5) +
labs(x = "LT-50 (cm)", y = "Densidade")
p0 + p1 + p2
A probabilidade posterior de \(\beta_1 > 0\) (comprimento de fato influencia positivamente a probabilidade de maturidade) é praticamente 1:
mean(posterior$beta1 > 0)[1] 1Comparando com o resumo posterior do LT-50:
tibble(
parametro = c("beta0_freq", "beta1_freq", "LT-50_freq",
"beta0_bayes", "beta1_bayes", "LT-50_bayes"),
estimativa = c(unname(beta_hat[1]),
unname(beta_hat[2]),
LT50_hat,
mean(posterior$beta0),
mean(posterior$beta1),
mean(posterior$LT50)),
EP = c(sqrt(V_beta[1, 1]),
sqrt(V_beta[2, 2]),
ep_LT50,
sd(posterior$beta0),
sd(posterior$beta1),
sd(posterior$LT50)),
inf95 = c(beta_hat[1] - 1.96 * sqrt(V_beta[1, 1]),
beta_hat[2] - 1.96 * sqrt(V_beta[2, 2]),
LT50_hat - 1.96 * ep_LT50,
quantile(posterior$beta0, 0.025),
quantile(posterior$beta1, 0.025),
quantile(posterior$LT50, 0.025)),
sup95 = c(beta_hat[1] + 1.96 * sqrt(V_beta[1, 1]),
beta_hat[2] + 1.96 * sqrt(V_beta[2, 2]),
LT50_hat + 1.96 * ep_LT50,
quantile(posterior$beta0, 0.975),
quantile(posterior$beta1, 0.975),
quantile(posterior$LT50, 0.975))
)# A tibble: 6 × 5
parametro estimativa EP inf95 sup95
<chr> <dbl> <dbl> <dbl> <dbl>
1 beta0_freq 1.16 1.14 -1.07 3.40
2 beta1_freq 0.267 0.146 -0.0188 0.553
3 LT-50_freq 28.1 3.11 22.1 34.2
4 beta0_bayes 1.69 1.35 -0.490 4.74
5 beta1_bayes 0.375 0.172 0.114 0.768
6 LT-50_bayes 28.2 3.22 22.6 34.9 ## Compara posterior com assintotica
LT50_grid <- seq(min(posterior$LT50), max(posterior$LT50), length.out = 200)
dens_freq <- tibble(
LT50 = LT50_grid,
densidade = dnorm(LT50_grid, mean = LT50_hat, sd = ep_LT50)
)
posterior |>
ggplot(aes(LT50)) +
geom_density(aes(color = "Bayesiana"),
fill = "darkorange", alpha = 0.5,
linewidth = 0.8) +
geom_line(data = dens_freq,
aes(x = LT50, y = densidade,
color = "Frequentista (assintótica)"),
linewidth = 0.8, linetype = "dashed") +
scale_color_manual(values = c("Bayesiana" = "darkorange",
"Frequentista (assintótica)" = "steelblue")) +
labs(x = "LT-50 (cm)", y = "Densidade", color = NULL) +
xlim(0, 60) +
theme(legend.position = "bottom")
Embora as estimativas pontuais sejam próximas, o intervalo frequentista é forçosamente simétrico em torno de \(\widehat{LT\text{-}50}\) e depende de uma aproximação assintótica, enquanto o intervalo de credibilidade bayesiano respeita a assimetria da posterior, o que é particularmente relevante aqui, dado o tamanho amostral pequeno (\(n_{total} = 17\)).
2.4 Tabela decisória e perdas esperadas
O LT-50 é discretizado nas seis classes de comprimento que definem os estados da natureza (\(\theta_1, \ldots, \theta_6\)) e calcula-se a frequência relativa em cada classe a partir da amostra posterior:
classes_lt <- cut(
posterior$LT50,
breaks = c(-Inf, 20, 25, 30, 35, 40, Inf),
right = FALSE
)
(prob_classes <- prop.table(table(classes_lt)))classes_lt
[-Inf,20) [20,25) [25,30) [30,35) [35,40) [40, Inf)
0.0058 0.1302 0.6112 0.2286 0.0200 0.0042 A matriz de perdas \(c_{i,j}\) associa a cada par (rede, classe) uma penalização que reflete o quanto aquela rede é inadequada para aquela faixa de tamanho. Os valores foram fixados ad hoc, penalizando mais fortemente a captura de imaturos:
perdas <- matrix(
c(1, 3, 4,
0, 2, 3,
2, 1, 2,
4, 0, 1,
6, 2, 0,
8, 4, 2),
nrow = 3,
dimnames = list(
decisao = c("A", "B", "C"),
classe = names(prob_classes)
)
)
perdas classe
decisao [-Inf,20) [20,25) [25,30) [30,35) [35,40) [40, Inf)
A 1 0 2 4 6 8
B 3 2 1 0 2 4
C 4 3 2 1 0 2mat_perdas_prob <- rbind(
perdas,
"p(theta|dados)" = round(as.numeric(prob_classes), 3)
)
knitr::kable(
mat_perdas_prob,
caption = "Tabela 1.2 -- Matriz de perdas $c_{i,j}$ associadas
a cada par (rede, classe de LT-50). A última linha contém as
probabilidades posteriores de cada classe.",
align = "c"
)| [-Inf,20) | [20,25) | [25,30) | [30,35) | [35,40) | [40, Inf) | |
|---|---|---|---|---|---|---|
| A | 1.000 | 0.00 | 2.000 | 4.000 | 6.00 | 8.000 |
| B | 3.000 | 2.00 | 1.000 | 0.000 | 2.00 | 4.000 |
| C | 4.000 | 3.00 | 2.000 | 1.000 | 0.00 | 2.000 |
| p(theta|dados) | 0.006 | 0.13 | 0.611 | 0.229 | 0.02 | 0.004 |
A perda esperada de cada decisão é o produto matricial das perdas pelas probabilidades posteriores:
perda_esp <- perdas %*% as.numeric(prob_classes)
res <- tibble(
decisao = c("A", "B", "C"),
perda_esperada = as.numeric(perda_esp)
)
res# A tibble: 3 × 2
decisao perda_esperada
<chr> <dbl>
1 A 2.30
2 B 0.946
3 C 1.87 A Decisão de Bayes é aquela que minimiza a perda esperada:
res |>
slice_min(perda_esperada, n = 1)# A tibble: 1 × 2
decisao perda_esperada
<chr> <dbl>
1 B 0.946A tabela decisória completa, contendo as perdas, as probabilidades posteriores dos estados da natureza e as perdas esperadas por decisão, sintetiza toda a análise:
## Adiciona coluna "Perda esperada" às linhas das decisoes
mat_dec <- cbind(perdas,
"Perda esperada" = round(as.numeric(perda_esp), 3))
## Adiciona linha com p(theta|dados); NA na ultima coluna
mat_res <- rbind(mat_dec,
c(round(as.numeric(prob_classes), 3), NA))
rownames(mat_res)[4] <- "p(theta|dados)"
knitr::kable(
mat_res,
caption = "Tabela decisória completa (similar à Tabela 6.1):
perdas $c_{i,j}$ no corpo, probabilidades posteriores dos
estados da natureza na última linha e perdas esperadas
$\\bar{c}(\\delta_i)$ na última coluna.",
align = "c"
)| [-Inf,20) | [20,25) | [25,30) | [30,35) | [35,40) | [40, Inf) | Perda esperada | |
|---|---|---|---|---|---|---|---|
| A | 1.000 | 0.00 | 2.000 | 4.000 | 6.00 | 8.000 | 2.296 |
| B | 3.000 | 2.00 | 1.000 | 0.000 | 2.00 | 4.000 | 0.946 |
| C | 4.000 | 3.00 | 2.000 | 1.000 | 0.00 | 2.000 | 1.873 |
| p(theta|dados) | 0.006 | 0.13 | 0.611 | 0.229 | 0.02 | 0.004 | NA |
Conclui-se que a rede B é a melhor escolha. É instrutivo notar que a rede A, mesmo tendo perda nula no estado \(\theta_2 = [20, 25)\), torna-se a pior alternativa quando se pondera pelas probabilidades posteriores – justamente porque o LT-50 está concentrado na faixa 25-30 cm, e a estrutura de perdas penaliza fortemente a captura de imaturos.
3 Parte II – Formalização
Componentes
Toda análise bayesiana de decisão organiza quatro elementos:
Decisões. Um conjunto exaustivo e exclusivo de alternativas \(\Delta = \{\delta_1, \ldots, \delta_m\}\). No exemplo, \(\Delta = \{A, B, C\}\).
Estados da natureza. Um conjunto exaustivo e exclusivo de eventos incertos \(\{\theta_1, \ldots, \theta_n\}\), ou seja, as hipóteses alternativas sobre o “verdadeiro” valor do parâmetro. No exemplo, as seis classes em que o LT-50 pode estar.
Função de perda. A cada par \((\delta_i, \theta_j)\) associa-se uma consequência \(c_{i,j} = L(\delta_i, \theta_j)\). Quando as consequências são desfavoráveis, fala-se em perdas; quando favoráveis, em utilidades.
Probabilidades posteriores. A descrição completa da incerteza sobre os \(\theta_j\) após observar os dados \(\mathbf{x}\) está contida na poosterior \(p(\theta_j \mid \mathbf{x})\).
A perda esperada
A perda esperada de uma decisão \(\delta_i\) é a média ponderada das perdas \(c_{i,j}\) pelas probabilidades posteriores dos estados:
\[ \bar{c}(\delta_i) \;=\; \sum_{j=1}^{n} c_{i,j}\, p(\theta_j \mid \mathbf{x}). \]
No caso contínuo, com \(\theta \in \Theta\) e função de perda \(L(\theta, \delta)\),
\[ \bar{c}(\delta) \;=\; \int_{\Theta} L(\theta, \delta)\, p(\theta \mid \mathbf{x})\, d\theta \;=\; \mathbb{E}_{p(\theta \mid \mathbf{x})}[L(\theta, \delta)]. \]
Definição da decisão de Bayes
NotaDefinição (Decisão de Bayes)
No conjunto \(\Delta = \{\delta_i;\, i = 1, 2, \ldots, m\}\) de todas as decisões possíveis, denomina-se decisão de Bayes \(\delta_B\) aquela que minimiza a perda esperada:
\[ \bar{c}(\delta_B) \;=\; \min_{\delta_i \in \Delta} \,\bar{c}(\delta_i). \]
Perda esperada como medida de risco
Note-se que \(\bar{c}(\delta_i)\) integra dois ingredientes distintos: as probabilidades dos estados da natureza e as magnitudes das consequências indesejáveis. Por isso, perdas esperadas são medidas naturais de risco – não basta um evento ser provável para ser arriscado; é preciso que suas consequências também sejam relevantes. A decisão de Bayes é, portanto, a alternativa que minimiza o risco.
Sensibilidade
Duas fontes podem alterar a decisão ótima:
- Mudanças nas probabilidades \(p(\theta_j \mid \mathbf{x})\), decorrentes de novos dados ou de prioris distintas;
- Mudanças na matriz de perdas \(c_{i,j}\), decorrentes de reavaliações das consequências (mercado, sustentabilidade, custos).
Uma vez estruturada a tabela decisória, qualquer alteração em qualquer um desses elementos é facilmente incorporada, garantindo coerência de longo prazo nas decisões.
4 Parte III – Aplicação à ANOVA: tubarão-martelo
4.1 O problema da CPUE
A pescaria de espinhel pelágico no sul do Brasil captura, entre outras espécies, o tubarão-martelo (Sphyrna spp.). Os dados de Kotas et al. (2005) registram capturas em três áreas oceânicas (9, 14 e 15) e o interesse é comparar a captura por unidade de esforço (CPUE), definida como kg capturados por 1000 anzóis. Como as variâncias da CPUE são fortemente heterogêneas entre áreas, adota-se como variável resposta o \(y = \log(\text{CPUE})\), em que a homocedasticidade é razoável.
O modelo de ANOVA com um fator fixo (área) é
\[ y_{g,i} = \mu_g + w_{g,i},\qquad w_{g,i} \sim N(0, \sigma^2),\quad g = 1, 2, 3,\ i = 1, \ldots, n_g, \]
com \(\mu_g\) sendo a média na escala log para a área \(g\).
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 do Exemplo 7.3 (Kotas et al., 2005)
## Tres areas de pesca (9, 14, 15); resposta: log(CPUE)
if(!file.exists("martelo.csv")) {
download.file(
"http://leg.ufpr.br/~fernandomayer/data/martelo.csv", "martelo.csv")
}
da <- read_csv("martelo.csv")
## 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))
## Calcula log(CPUE)
da$lcpue <- log(da$cpue)
## Estatisticas suficientes
fnsq <- function(x) sum((x - mean(x))^2)
(ygbar <- with(da, tapply(lcpue, areaf, mean)))
(sqg <- with(da, tapply(lcpue, area, fnsq)))
(ng_obs <- with(da, tapply(lcpue, area, length)))
(nG <- sum(ng_obs))
(G <- length(ng_obs))
(se2G <- sum(sqg) / (nG - G))
## Simulacao da posterior conjunta (m = 3000)
m <- 3000
(Vmu <- diag(1/ng_obs))
tauG <- 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))
}
## 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)## Probabilidade P(mu15 - mu9 > 0 | dados)
(prob15_9 <- mean(emuGpost$mu15 - emuGpost$mu9 > 0))
## Diferenca posterior
dif15_9 <- emuGpost$mu15 - emuGpost$mu9
mean(dif15_9 > 0)
mean(dif15_9)
icr(dif15_9)
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")4.2 Exercícios
4.2.1 Exercício 7.6 – Tabela decisória para escolha de área
Vimos que há duas áreas (9 e 15) com CPUE média semelhante, mas a área 15 é estocasticamente maior, \(P(\mu_{15} - \mu_{9} > 0 \mid \text{dados}) \approx 0{,}73\). À primeira vista, a área 15 seria a melhor opção. No entanto, a área 9 é mais próxima do porto e o custo operacional (combustível) é menor. Suponha que os critérios de avaliação foram somente o custo operacional e a diferença de CPUE. A tabela decisória resultante, com \(\theta = \mu_{15} - \mu_{9}\), é
| \(\theta > 0\) | \(\theta \leq 0\) | |
|---|---|---|
| \(\delta_9\) (área 9) | 4 | 0 |
| \(\delta_{15}\) (área 15) | 2 | 6 |
| \(p(\theta\mid\text{dados})\) | 0.73 | 0.27 |
A menor perda ocorreria na situação em que se usa a área 9 e ela de fato apresenta a maior CPUE média: o melhor cenário para o pescador. A pior é ir pescar longe (área 15) e a CPUE de lá ser, ao final, menor. As demais perdas representam situações intermediárias.
Calcule as perdas esperadas para cada decisão.
4.2.2 Resolução
Da amostra posterior obtemos \(P(\theta > 0)\) e \(P(\theta \leq 0)\):
theta <- dif15_9
(p_pos <- mean(theta > 0))[1] 0.72(p_neg <- 1 - p_pos)[1] 0.28## A matriz de perdas é
perdas_ex <- matrix(
c(4, 2,
0, 6),
nrow = 2,
dimnames = list(
decisao = c("delta_9", "delta_15"),
estado = c("theta > 0", "theta <= 0")
)
)
perdas_ex estado
decisao theta > 0 theta <= 0
delta_9 4 0
delta_15 2 6## A perda esperada de cada decisão é o produto das perdas pelas
## probabilidades posteriores:
prob_theta <- c(p_pos, p_neg)
perda_esp_ex <- perdas_ex %*% prob_theta
res_ex <- tibble(
decisao = c("delta_9 (area 9)", "delta_15 (area 15)"),
perda_esperada = as.numeric(perda_esp_ex)
)
res_ex# A tibble: 2 × 2
decisao perda_esperada
<chr> <dbl>
1 delta_9 (area 9) 2.88
2 delta_15 (area 15) 3.12## Tabela decisória completa para o exercício:
mat_ex <- cbind(perdas_ex,
"Perda esperada" = round(as.numeric(perda_esp_ex), 3))
mat_ex <- rbind(mat_ex,
c(round(prob_theta, 3), NA))
rownames(mat_ex)[3] <- "p(theta|dados)"
knitr::kable(
mat_ex,
caption = "Tabela decisória do Exercício 7.6:
escolha entre área 9 e área 15.",
align = "c"
)| theta > 0 | theta <= 0 | Perda esperada | |
|---|---|---|---|
| delta_9 | 4.00 | 0.00 | 2.88 |
| delta_15 | 2.00 | 6.00 | 3.12 |
| p(theta|dados) | 0.72 | 0.28 | NA |
## A decisão de Bayes:
res_ex |>
slice_min(perda_esperada, n = 1)# A tibble: 1 × 2
decisao perda_esperada
<chr> <dbl>
1 delta_9 (area 9) 2.88Embora a área 15 seja estocasticamente maior em CPUE, ao considerar o custo operacional adicional – representado pela perda \(c_{15,\theta\leq 0} = 6\) no cenário em que a aposta na área 15 não se confirma – a perda esperada da decisão \(\delta_{15}\) supera a de \(\delta_9\).
A decisão de Bayes recomenda a área 9, mais próxima do porto, com custo mais baixo, ainda que com CPUE média ligeiramente menor. O exemplo ilustra bem o princípio formalizado na Parte II: não basta um estado da natureza ser provável; é a combinação de probabilidade com magnitude da consequência que define a escolha ótima.
4.2.3 Exercício complementar – Conversão de perdas em utilidades
Em teoria da decisão, a decisão de Bayes é determinada apenas pela ordem relativa dos valores esperados, não por seus valores absolutos. Mais precisamente, se
\[ u(\theta, \delta) = a + b \cdot c(\theta, \delta), \]
com \(a, b\) constantes e \(b \neq 0\), então pela linearidade da esperança
\[ \mathbb{E}[u(\theta, \delta)] = a + b \cdot \mathbb{E}[c(\theta, \delta)]. \]
A ordenação entre decisões é preservada porque \(a\) é constante (soma-se igualmente a todas) e \(b\) é o mesmo escalar (multiplica todas igualmente). Logo:
- se \(b > 0\), maximizar \(\mathbb{E}[u]\) é equivalente a maximizar \(\mathbb{E}[c]\);
- se \(b < 0\), maximizar \(\mathbb{E}[u]\) é equivalente a minimizar \(\mathbb{E}[c]\).
Uma parametrização conveniente
Para converter a matriz de perdas \(c_{i,j}\) do Exercício 7.6 em uma matriz de utilidades, considere a transformação decrescente
\[ u_{i,j} = (M - c_{i,j}) \cdot k, \]
que corresponde a \(a = Mk\) e \(b = -k\), com:
\(M\): uma constante aditiva que define o “ponto zero” da utilidade – a perda contra a qual todos os outros cenários são comparados. Uma escolha conveniente é \(M = \max_{i,j} c_{i,j}\), pois assim o pior cenário recebe utilidade zero (lucro nulo) e todos os outros recebem utilidade positiva. É uma escolha interpretativa, não matemática: qualquer outro valor de \(M\) preservaria a decisão de Bayes.
\(k\): o fator de escala – converte as “unidades de perda” originais (que eram adimensionais, fixadas ad hoc na elaboração da tabela) em unidades monetárias. Uma escolha como \(k = 1000\) significa “uma unidade de perda equivale a R$ 1.000”. Também é arbitrário: a decisão de Bayes não muda com \(k\).
Sobre as unidades
A unidade resultante depende do que se interpreta como “uma unidade de perda” na tabela original. As perdas em \(c_{i,j}\) não têm unidade declarada, são índices ad hoc. Ao multiplicar por \(k = 1000\) R$, estamos implicitamente afirmando que uma unidade de perda equivale a R$ 1.000 de ganho/prejuízo líquido por viagem decisória – não R$ por kg nem R$ por 1000 anzóis, mas R$ por escolha entre as áreas 9 e 15 numa viagem completa, abrangendo combustível, tripulação, dias de trabalho e valor de mercado da captura. Sem ancoragem em dados reais, \(k\) é apenas uma escala nominal.
Pergunta
Com \(M = \max_{i,j} c_{i,j} = 6\) e \(k = 1000\), como você construiria a tabela de utilidades \(u_{i,j}\) a partir da matriz de perdas do Exercício 7.6 e calcularia a utilidade esperada de cada decisão? Qual seria a área de maior utilidade esperada? Verifique se a decisão de Bayes obtida via maximização da utilidade coincide com aquela obtida via minimização da perda.
4.2.4 Resolução
## Parametros da transformacao afim decrescente
M <- max(perdas_ex)
k <- 1000
## Matriz de utilidades: u_ij = (M - c_ij) * k
(util_ex <- (M - perdas_ex) * k) estado
decisao theta > 0 theta <= 0
delta_9 2000 6000
delta_15 4000 0## Utilidade esperada: produto matricial pelas probabilidades
(util_esp_ex <- util_ex %*% prob_theta)
decisao [,1]
delta_9 3120
delta_15 2880## Tabela completa
mat_util <- cbind(util_ex,
"Utilidade esperada" = round(as.numeric(util_esp_ex), 0))
mat_util <- rbind(mat_util,
c(round(prob_theta, 3), NA))
rownames(mat_util)[3] <- "p(theta|dados)"
knitr::kable(
mat_util,
caption = "Tabela decisória do Exercício 7.6 reformulada em
utilidades (em R$): ganhos $u_{i,j}$ no corpo, probabilidades
posteriores na última linha e utilidades esperadas na última
coluna.",
align = "c"
)| theta > 0 | theta <= 0 | Utilidade esperada | |
|---|---|---|---|
| delta_9 | 2000.00 | 6000.00 | 3120 |
| delta_15 | 4000.00 | 0.00 | 2880 |
| p(theta|dados) | 0.72 | 0.28 | NA |
## Decisao de Bayes: maximiza utilidade esperada
res_util <- tibble(
decisao = c("delta_9 (area 9)", "delta_15 (area 15)"),
utilidade_esperada = as.numeric(util_esp_ex)
)
res_util |>
slice_max(utilidade_esperada, n = 1)# A tibble: 1 × 2
decisao utilidade_esperada
<chr> <dbl>
1 delta_9 (area 9) 3120A área 9 é novamente a decisão de Bayes, agora com utilidade esperada de R$ 3120 contra R$ 2880 da área 15. A diferença entre as utilidades esperadas (R$ 0.24) é exatamente \(k\) vezes a diferença entre as perdas esperadas calculadas anteriormente (\(\bar{c}(\delta_{15}) - \bar{c}(\delta_9) = 3.12 - 2.88 \approx 0.24\)),
diff(as.numeric(perda_esp_ex))[1] 0.24diff(as.numeric(perda_esp_ex)) * k[1] 240confirmando numericamente que a ordenação entre alternativas é preservada, como previsto pela teoria. A escolha entre trabalhar com perdas ou utilidades é, portanto, apenas uma questão de interpretação: “menor risco” no primeiro caso, “maior lucro esperado” no segundo.