library(tidyverse)
library(patchwork)
theme_set(theme_bw())
set.seed(42)Modelo Beta-Binomial Hierárquico
Simulação e Inferência Posterior
Especificação do modelo
O modelo Beta-Binomial hierárquico assume a seguinte estrutura:
\[ Y_i \mid \theta_i \sim \text{Binomial}(n_i, \theta_i), \quad i = 1, \ldots, N \]
\[ \theta_i \mid \omega, \kappa \sim \text{Beta}(\alpha, \beta) \]
\[ \omega \sim \text{Beta}(a, b) \]
onde \(\kappa\) é tratado como fixado. A reparametrização por moda e concentração relaciona \((\alpha, \beta)\) com \((\omega, \kappa)\):
\[ \omega = \frac{\alpha - 1}{\alpha + \beta - 2}, \qquad \kappa = \alpha + \beta \]
de modo que
\[ \alpha = (\kappa - 2)\,\omega + 1, \qquad \beta = (\kappa - 2)\,(1 - \omega) + 1. \]
Aqui, \(\omega\) representa a moda da distribuição populacional dos \(\theta_i\) e \(\kappa\) controla a concentração em torno dessa moda.
Simulação dos dados
Os dados são gerados seguindo a estrutura hierárquica do modelo. Primeiro, amostramos \(\omega\) da priori \(\text{Beta}(a, b)\). Em seguida, calculamos \(\alpha\) e \(\beta\) pela reparametrização e simulamos os parâmetros individuais \(\theta_i\) e, por fim, as observações \(Y_i\).
## Parâmetros fixados
N <- 8 # número de observações
kappa <- 20 # concentração (fixada)
a <- 2 # hiperparâmetro da prior de omega
b <- 2 # hiperparâmetro da prior de omega
(n <- sample(10:50, N, replace = TRUE)) # tamanhos amostrais n_i[1] 46 10 34 19 45 27 33 16## Passo 1: Simular w ~ Beta(a, b)
(w <- rbeta(1, shape1 = a, shape2 = b))[1] 0.3194048## Passo 2: Calcular alpha e beta via reparametrização
alpha_par <- (kappa - 2) * w + 1
beta_par <- (kappa - 2) * (1 - w) + 1
c(alpha = alpha_par, beta = beta_par) alpha beta
6.749286 13.250714 ## Passo 3: Simular theta_i ~ Beta(alpha, beta)
theta <- rbeta(N, shape1 = alpha_par, shape2 = beta_par)
## Passo 4: Simular Y_i ~ Binomial(n_i, theta_i)
y <- rbinom(N, size = n, prob = theta)
## Dataset resultante
da <- tibble(
i = 1:N,
n_i = n,
theta_i = theta,
y_i = y,
prop = y/n
)Os dados simulados são apresentados abaixo. A coluna prop contém a proporção observada \(\hat{p}_i = y_i / n_i\).
da# A tibble: 8 × 5
i n_i theta_i y_i prop
<int> <int> <dbl> <int> <dbl>
1 1 46 0.205 12 0.261
2 2 10 0.356 0 0
3 3 34 0.573 22 0.647
4 4 19 0.342 9 0.474
5 5 45 0.522 23 0.511
6 6 27 0.468 12 0.444
7 7 33 0.454 15 0.455
8 8 16 0.398 3 0.188Visualização da simulação
O painel esquerdo mostra a priori \(\omega \sim \text{Beta}(a, b)\) com o valor simulado (linha vermelha). O painel central mostra a densidade \(\text{Beta}(\alpha, \beta)\) da qual os \(\theta_i\) foram sorteados (pontos vermelhos). O painel direito compara os verdadeiros \(\theta_i\) com as proporções observadas \(\hat{p}_i\).
Inferência posterior — Modelo hierárquico
Nos modelos a seguir, \(\kappa = 20\) é tratado como fixado. O objetivo é inferir a moda populacional \(\omega\) e os parâmetros individuais \(\theta_i\).
Nos gráficos, o ponto/linha azul representa a estimativa posterior (média), a linha pontilhada azul indica a média posterior de \(\omega\), e o \(\times\) vermelho marca o valor verdadeiro simulado.
JAGS
No JAGS, \(\omega\) e cada \(\theta_i\) são amostrados explicitamente via MCMC. O parâmetro \(\kappa\) é passado como dado.
Resumo das estatísticas posteriores (JAGS):
Mean SD Naive SE Time-series SE
omega 0.3915665 0.05518700 0.0004506000 0.0009023320
theta[1] 0.3033539 0.05718258 0.0004668938 0.0006323449
theta[2] 0.2680786 0.08601996 0.0007023500 0.0011032047
theta[3] 0.5558269 0.07012064 0.0005725327 0.0007764277
theta[4] 0.4367178 0.08302786 0.0006779197 0.0009825717
theta[5] 0.4769087 0.06308557 0.0005150915 0.0007446719
theta[6] 0.4261183 0.07438395 0.0006073424 0.0008358378
theta[7] 0.4344975 0.06990150 0.0005707433 0.0007992912
theta[8] 0.3083573 0.08126922 0.0006635604 0.0009767386Estimativas posteriores dos \(\theta_i\) com intervalos de credibilidade de 95%:
# A tibble: 8 × 7
i n_i y_i theta_i mean_post lower upper
<int> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 1 46 12 0.205 0.303 0.196 0.421
2 2 10 0 0.356 0.268 0.119 0.453
3 3 34 22 0.573 0.556 0.417 0.691
4 4 19 9 0.342 0.437 0.279 0.603
5 5 45 23 0.522 0.477 0.356 0.599
6 6 27 12 0.468 0.426 0.285 0.576
7 7 33 15 0.454 0.434 0.300 0.573
8 8 16 3 0.398 0.308 0.163 0.482
Stan
O Stan estima \(\omega\) e \(\theta_i\) como parâmetros, com \(\alpha\) e \(\beta\) calculados no bloco transformed parameters.
Resumo dos parâmetros populacionais (Stan):
Inference for Stan model: anon_model.
3 chains, each with iter=4000; warmup=1000; thin=1;
post-warmup draws per chain=3000, total post-warmup draws=9000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
omega 0.39 0.00 0.06 0.28 0.35 0.39 0.43 0.50 9214 1
alpha_par 8.05 0.01 1.00 6.08 7.37 8.04 8.71 10.02 9214 1
beta_par 11.95 0.01 1.00 9.98 11.29 11.96 12.63 13.92 9214 1
Samples were drawn using NUTS(diag_e) at Wed Mar 25 21:14:30 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Resumo completo:
mean se_mean sd 2.5% 25%
omega 0.3914937 0.0005778848 0.05547175 0.2821015 0.3540462
theta[1] 0.3037390 0.0004807565 0.05859146 0.1928141 0.2633462
theta[2] 0.2683474 0.0007568800 0.08667727 0.1162772 0.2058847
theta[3] 0.5566788 0.0005575710 0.06887388 0.4193662 0.5105903
theta[4] 0.4368292 0.0006927684 0.08189060 0.2782478 0.3801587
theta[5] 0.4777244 0.0005176590 0.06331711 0.3539386 0.4344812
theta[6] 0.4264163 0.0006629314 0.07506989 0.2821794 0.3751180
theta[7] 0.4341967 0.0005591783 0.07032628 0.2981355 0.3866304
theta[8] 0.3067798 0.0007003285 0.08060128 0.1601622 0.2491753
alpha_par 8.0468859 0.0104019264 0.99849157 6.0778271 7.3728311
beta_par 11.9531141 0.0104019264 0.99849157 9.9754252 11.2935305
lp__ -158.1190617 0.0341722596 2.17876438 -163.2563427 -159.3970877
50% 75% 97.5% n_eff Rhat
omega 0.3910284 0.4281372 0.5013653 9214.277 1.0001456
theta[1] 0.3016535 0.3426389 0.4231404 14853.139 0.9998779
theta[2] 0.2622097 0.3240079 0.4552193 13114.641 1.0000111
theta[3] 0.5569401 0.6040399 0.6878695 15258.390 0.9997493
theta[4] 0.4363420 0.4921073 0.5979500 13973.072 1.0000399
theta[5] 0.4769579 0.5209598 0.5985619 14960.794 0.9998307
theta[6] 0.4255054 0.4760301 0.5759792 12823.141 0.9998393
theta[7] 0.4333382 0.4811060 0.5750683 15817.383 0.9999405
theta[8] 0.3037443 0.3586832 0.4738452 13245.864 0.9997781
alpha_par 8.0385110 8.7064695 10.0245748 9214.277 1.0001456
beta_par 11.9614890 12.6271689 13.9221729 9214.277 1.0001456
lp__ -157.7785814 -156.5002031 -154.9093665 4065.118 1.0000293Estimativas posteriores dos \(\theta_i\) (Stan):
# A tibble: 8 × 7
i n_i y_i theta_i mean_post lower upper
<int> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 1 46 12 0.205 0.304 0.193 0.423
2 2 10 0 0.356 0.268 0.116 0.455
3 3 34 22 0.573 0.557 0.419 0.688
4 4 19 9 0.342 0.437 0.278 0.598
5 5 45 23 0.522 0.478 0.354 0.599
6 6 27 12 0.468 0.426 0.282 0.576
7 7 33 15 0.454 0.434 0.298 0.575
8 8 16 3 0.398 0.307 0.160 0.474
Modelo não-hierárquico: \(\omega\) e \(\kappa\) fixados
Nesta seção, tanto \(\omega\) quanto \(\kappa\) são tratados como conhecidos. Os parâmetros \(\alpha\) e \(\beta\) ficam totalmente determinados:
\[ \alpha = (\kappa - 2)\,\omega + 1, \qquad \beta = (\kappa - 2)\,(1 - \omega) + 1. \]
A única inferência é sobre cada \(\theta_i\). A verossimilhança Binomial combinada com a priori \(\text{Beta}(\alpha, \beta)\) produz uma posterior conjugada exata:
\[ \theta_i \mid y_i \sim \text{Beta}(\alpha + y_i,\; \beta + n_i - y_i). \]
Usamos o valor verdadeiro simulado como \(\omega\) fixado; na prática, seria um valor escolhido pelo analista.
omega alpha beta
0.3194048 6.7492860 13.2507140 Analítico (exato)
A posterior conjugada permite calcular a média e os quantis diretamente, sem MCMC.
Estimativas posteriores dos \(\theta_i\) — método analítico (\(\omega\) fixo):
# A tibble: 8 × 7
i n_i y_i theta_i mean_post lower upper
<int> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 1 46 12 0.205 0.284 0.183 0.398
2 2 10 0 0.356 0.225 0.0971 0.387
3 3 34 22 0.573 0.532 0.400 0.663
4 4 19 9 0.342 0.404 0.257 0.560
5 5 45 23 0.522 0.458 0.339 0.579
6 6 27 12 0.468 0.399 0.265 0.541
7 7 33 15 0.454 0.410 0.283 0.544
8 8 16 3 0.398 0.271 0.141 0.425
JAGS (\(\omega\) fixado)
Os parâmetros \(\alpha\) e \(\beta\) são passados diretamente como dados. O modelo contém apenas os \(\theta_i\) como parâmetros.
Resumo das estatísticas posteriores (JAGS, \(\omega\) fixo):
Mean SD Naive SE Time-series SE
theta[1] 0.2833574 0.05481086 0.0004475288 0.0005624095
theta[2] 0.2258850 0.07391789 0.0006035370 0.0007903427
theta[3] 0.5332618 0.06730760 0.0005495642 0.0006891883
theta[4] 0.4037323 0.07782805 0.0006354634 0.0007929290
theta[5] 0.4586859 0.06116554 0.0004994146 0.0006203019
theta[6] 0.3987718 0.07020680 0.0005732362 0.0007138665
theta[7] 0.4103670 0.06717122 0.0005484507 0.0006960785
theta[8] 0.2698853 0.07245511 0.0005915935 0.0007598613Estimativas posteriores dos \(\theta_i\) (JAGS, \(\omega\) fixo):
# A tibble: 8 × 7
i n_i y_i theta_i mean_post lower upper
<int> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 1 46 12 0.205 0.283 0.180 0.395
2 2 10 0 0.356 0.226 0.0990 0.383
3 3 34 22 0.573 0.533 0.403 0.663
4 4 19 9 0.342 0.404 0.258 0.563
5 5 45 23 0.522 0.459 0.341 0.579
6 6 27 12 0.468 0.399 0.267 0.538
7 7 33 15 0.454 0.410 0.282 0.543
8 8 16 3 0.398 0.270 0.142 0.422
Stan (\(\omega\) fixado)
Nesta versão, \(\omega\) é passado como dado escalar. Os valores \(\alpha\) e \(\beta\) são calculados no bloco transformed data e não há parâmetros populacionais a estimar.
Resumo dos \(\theta_i\) (Stan, \(\omega\) fixo):
Inference for Stan model: anon_model.
3 chains, each with iter=4000; warmup=1000; thin=1;
post-warmup draws per chain=3000, total post-warmup draws=9000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta[1] 0.28 0 0.06 0.18 0.25 0.28 0.32 0.40 14129 1
theta[2] 0.22 0 0.08 0.10 0.17 0.22 0.27 0.39 14727 1
theta[3] 0.53 0 0.07 0.40 0.49 0.53 0.58 0.67 14628 1
theta[4] 0.40 0 0.08 0.26 0.35 0.40 0.46 0.56 15095 1
theta[5] 0.46 0 0.06 0.34 0.42 0.46 0.50 0.58 14224 1
theta[6] 0.40 0 0.07 0.27 0.35 0.40 0.45 0.54 15592 1
theta[7] 0.41 0 0.07 0.28 0.36 0.41 0.46 0.54 14970 1
theta[8] 0.27 0 0.07 0.14 0.22 0.27 0.32 0.42 15045 1
Samples were drawn using NUTS(diag_e) at Wed Mar 25 21:15:20 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1). mean se_mean sd 2.5% 25%
theta[1] 0.2840216 0.0004640256 0.05515605 0.18295001 0.2458039
theta[2] 0.2247145 0.0006223695 0.07552862 0.09682611 0.1692605
theta[3] 0.5319673 0.0005621011 0.06798385 0.39771643 0.4869593
theta[4] 0.4046748 0.0006392682 0.07854163 0.25640519 0.3495960
theta[5] 0.4572172 0.0005072923 0.06050172 0.33973721 0.4157865
theta[6] 0.3996578 0.0005734591 0.07160769 0.26692281 0.3489548
theta[7] 0.4100763 0.0005534089 0.06771014 0.28129451 0.3629778
theta[8] 0.2706741 0.0005873602 0.07204465 0.13976635 0.2197837
lp__ -256.4277223 0.0319167013 2.04621577 -261.32928336 -257.5783996
50% 75% 97.5% n_eff Rhat
theta[1] 0.2823004 0.3199977 0.3972347 14128.710 0.9999925
theta[2] 0.2180354 0.2731862 0.3877040 14727.415 0.9999048
theta[3] 0.5327747 0.5776296 0.6652741 14627.923 0.9997815
theta[4] 0.4032495 0.4584911 0.5616540 15095.017 0.9997556
theta[5] 0.4563727 0.4979507 0.5779065 14223.910 0.9997365
theta[6] 0.3978312 0.4490095 0.5408636 15592.450 1.0000576
theta[7] 0.4089764 0.4560760 0.5426223 14969.772 1.0001664
theta[8] 0.2679950 0.3174830 0.4224905 15045.077 0.9999032
lp__ -256.0803095 -254.9231425 -253.4843524 4110.237 0.9998810Estimativas posteriores dos \(\theta_i\) (Stan, \(\omega\) fixo):
# A tibble: 8 × 7
i n_i y_i theta_i mean_post lower upper
<int> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 1 46 12 0.205 0.284 0.183 0.397
2 2 10 0 0.356 0.225 0.0968 0.388
3 3 34 22 0.573 0.532 0.398 0.665
4 4 19 9 0.342 0.405 0.256 0.562
5 5 45 23 0.522 0.457 0.340 0.578
6 6 27 12 0.468 0.400 0.267 0.541
7 7 33 15 0.454 0.410 0.281 0.543
8 8 16 3 0.398 0.271 0.140 0.422
Comparação entre métodos (\(\omega\) fixo)
Os três métodos (Analítico, JAGS e Stan) produzem estimativas praticamente idênticas quando \(\omega\) é fixado, confirmando a consistência entre a solução conjugada exata e as aproximações MCMC.
Tabela consolidada e variâncias
A tabela abaixo reúne as estimativas de todas as abordagens. A proporção global \(\hat{\theta}_g = \sum y_i / \sum n_i\) serve como referência frequentista.
# A tibble: 8 × 8
n_i theta_i y_i prop theta_jags theta_anl_f theta_jags_f theta_g
<int> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 46 0.205 12 0.261 0.303 0.284 0.283 0.417
2 10 0.356 0 0 0.268 0.225 0.226 0.417
3 34 0.573 22 0.647 0.556 0.532 0.533 0.417
4 19 0.342 9 0.474 0.437 0.404 0.404 0.417
5 45 0.522 23 0.511 0.477 0.458 0.459 0.417
6 27 0.468 12 0.444 0.426 0.399 0.399 0.417
7 33 0.454 15 0.455 0.434 0.410 0.410 0.417
8 16 0.398 3 0.188 0.308 0.271 0.270 0.417Variância de cada coluna — observar o efeito de shrinkage (encolhimento) nas estimativas hierárquicas em relação às proporções observadas:
i n_i theta_i y_i prop
6.000000000 174.214285714 0.013400050 66.285714286 0.043169117
theta_jags theta_stan theta_anl_f theta_jags_f theta_stan_f
0.009792479 0.009869601 0.010836502 0.010906021 0.010832560
theta_g
0.000000000 Correlação entre estimativas
Modelo com \(\omega\) fixo
O gráfico de pares abaixo compara os verdadeiros \(\theta_i\), as proporções amostrais \(\hat{p}_i\) e as estimativas posteriores dos três métodos com \(\omega\) fixado.
Modelo fixo vs hierárquico
Comparação entre as estimativas do modelo com \(\omega\) fixo (JAGS) e do modelo hierárquico (JAGS), evidenciando o efeito de shrinkage do modelo hierárquico.
Comparação final: fixo vs hierárquico vs proporção observada
O gráfico abaixo sintetiza as três abordagens. No modelo com \(\omega\) fixo, cada \(\theta_i\) é estimado independentemente. No modelo hierárquico, as estimativas são “puxadas” em direção à média global (shrinkage), especialmente para as observações com tamanho amostral menor. A proporção observada \(\hat{p}_i\) não incorpora informação da priori.
