Teoria da Credibilidade — Laboratório 2, Parte 1 (2026/1)

Previsão de Mortalidade em Fundo de Pensão via Bühlmann-Straub

Autor

Arthur Pontes Motta e Catarine Martins

Data de Publicação

26 de junho de 2026

Os dados referem-se a um fundo de pensão brasileiro do tipo Entidade Fechada de Previdência Complementar (EFPC), com registros de expostos e óbitos por idade para os anos de 2012, 2013 e 2014. As EFPCs são fundações sem fins lucrativos vinculadas a empresas patrocinadoras; os participantes acumulam reservas para aposentadoria complementar e o plano cobre ainda o risco de morte do segurado, cujo benefício é repassado aos beneficiários indicados (ABRAPP, 2022). A supervisão é exercida pela Superintendência Nacional de Previdência Complementar (PREVIC).

O objetivo desta parte é obter o número esperado de sinistros por classe de risco para o próximo ano com base no modelo de Bühlmann-Straub para a frequência de sinistros, seguindo o desenvolvimento de Bühlmann; Gisler (2005). O caminho percorrido é: (i) análise exploratória; (ii) perfis etários e faixas de risco; (iii) seleção e avaliação formal de tábuas de referência; (iv) definição dos grupos; (v) implementação do modelo e previsão.

Leitura dos dados

Ocultar código 
fundo <- read.csv("dadosfundopensao.csv", header = TRUE, sep = ";")
cat("Dimensão:", paste(dim(fundo), collapse = " x "), "\n")
Dimensão: 116 x 7
Ocultar código 
fundo |> head(8) |>
  gt() |>
  tab_header(title = "Primeiras linhas — dados do fundo de pensão") |>
  tab_options(table.font.size = 13)
Primeiras linhas — dados do fundo de pensão
idade EXP2012 OBT2012 EXP2013 OBT2013 EXP2014 OBT2014
0 0 0 0 0 0 0
1 1 0 0 0 0 0
2 3 0 1 0 0 0
3 3 0 3 0 1 0
4 0 0 3 0 4 0
5 4 0 0 0 3 0
6 2 0 5 0 0 0
7 2 0 2 0 6 0
Ocultar código 
fundo_long <- tibble(
  Idade    = rep(fundo$idade, 3),
  Expostos = c(fundo$EXP2012, fundo$EXP2013, fundo$EXP2014),
  Obitos   = c(fundo$OBT2012, fundo$OBT2013, fundo$OBT2014),
  Ano      = rep(c("2012","2013","2014"), each = nrow(fundo))
) |>
  mutate(
    tx     = if_else(Expostos > 0, Obitos / Expostos, NA_real_),
    log_tx = if_else(tx > 0, log(tx), NA_real_)
  )

sum_d2 <- fundo_long |>
  group_by(Idade) |>
  summarise(tot_obt = sum(Obitos), tot_exp = sum(Expostos), .groups = "drop")

# ── Métricas para uso em texto inline ────────────────────────────────────────
resumo_anual <- fundo_long |>
  group_by(Ano) |>
  summarise(
    Expostos         = sum(Expostos),
    Obitos           = sum(Obitos),
    taxa_bruta_mil   = round(1000 * sum(Obitos) / sum(Expostos), 1),
    .groups = "drop"
  )

tot_exp_tri  <- sum(sum_d2$tot_exp)
tot_obt_tri  <- sum(sum_d2$tot_obt)
n_idades     <- nrow(sum_d2)
taxa_2012    <- resumo_anual$taxa_bruta_mil[resumo_anual$Ano == "2012"]
taxa_2014    <- resumo_anual$taxa_bruta_mil[resumo_anual$Ano == "2014"]
cresc_taxa   <- round((taxa_2014 - taxa_2012) / taxa_2012 * 100, 1)
exp_2012     <- resumo_anual$Expostos[resumo_anual$Ano == "2012"]
exp_2014     <- resumo_anual$Expostos[resumo_anual$Ano == "2014"]
cresc_exp    <- round((exp_2014 - exp_2012) / exp_2012 * 100, 1)

resumo_anual |>
  select(Ano, Expostos, Obitos, `Taxa bruta (‰)` = taxa_bruta_mil) |>
  gt() |>
  tab_header(title = "Resumo anual da carteira") |>
  fmt_number(columns = c(Expostos, Obitos), decimals = 0) |>
  tab_options(table.font.size = 13)
Resumo anual da carteira
Ano Expostos Obitos Taxa bruta (‰)
2012 137,797 1,131 8.2
2013 138,919 1,192 8.6
2014 141,248 1,335 9.5

A base cobre 116 idades (0 a 115) ao longo de \(T = 3\) anos, totalizando 417.964 pessoas-ano de exposição e 3.658 óbitos. A taxa bruta passou de 8.2‰ em 2012 para 9.5‰ em 2014, crescimento de 15.9% com exposição praticamente estável (+2.5%), revelando envelhecimento progressivo da carteira. Vale registrar que a Portaria PREVIC nº 835/2020 (PREVIC, 2020) exige, para adequação da tábua de mortalidade geral, no mínimo cinco exercícios — horizonte não atingido com esses dados, limitação que deve ser declarada explicitamente e que reforça a escolha do modelo de credibilidade como alternativa robusta para períodos curtos.

Questão (i): Análise Exploratória

Caracterização da massa e métricas de maturidade

Ocultar código 
# Idade média ponderada pela exposição (triênio)
idade_media <- with(sum_d2, sum(Idade * tot_exp) / sum(tot_exp))

# Exposição e óbitos acima de 60 anos
pct_obt_60p <- sum(sum_d2$tot_obt[sum_d2$Idade >= 60]) / sum(sum_d2$tot_obt)
pct_exp_60p <- sum(sum_d2$tot_exp[sum_d2$Idade >= 60]) / sum(sum_d2$tot_exp)

# Índice de dependência (expostos 60+ / expostos 18-59)
dep_idx <- sum(sum_d2$tot_exp[sum_d2$Idade >= 60]) /
           sum(sum_d2$tot_exp[sum_d2$Idade >= 18 & sum_d2$Idade < 60])

# Variáveis para uso inline
idade_media_r  <- round(idade_media, 1)
dep_idx_r      <- round(dep_idx, 3)
pct_obt_60p_r  <- round(pct_obt_60p * 100, 1)
pct_exp_60p_r  <- round(pct_exp_60p * 100, 1)

tibble(
  Métrica = c(
    "Exposição total (pessoas-ano)",
    "Óbitos totais",
    "Idade média ponderada pela exposição",
    "Taxa bruta média 2012--2014 (‰)",
    "% óbitos acima de 60 anos",
    "% expostos acima de 60 anos",
    "Índice de dependência (60+)/(18--59)"
  ),
  Valor = c(
    sum(sum_d2$tot_exp),
    sum(sum_d2$tot_obt),
    idade_media_r,
    round(1000 * sum(sum_d2$tot_obt) / sum(sum_d2$tot_exp), 3),
    pct_obt_60p_r,
    pct_exp_60p_r,
    dep_idx_r
  )
) |>
  gt() |>
  tab_header(title = "Métricas de caracterização da carteira (triênio 2012--2014)") |>
  tab_options(table.font.size = 13)
Métricas de caracterização da carteira (triênio 2012--2014)
Métrica Valor
Exposição total (pessoas-ano) 417964.000
Óbitos totais 3658.000
Idade média ponderada pela exposição 45.200
Taxa bruta média 2012--2014 (‰) 8.752
% óbitos acima de 60 anos 89.700
% expostos acima de 60 anos 26.400
Índice de dependência (60+)/(18--59) 0.359

A idade média ponderada de \(\approx 45.2\) anos e o índice de dependência de \(\approx 0.359\) apontam para um plano em maturidade intermediária: a massa ainda é predominantemente ativa, mas o peso dos assistidos nas idades avançadas já responde por 89.7% dos óbitos. Planos com índice de dependência acima de \(0{,}5\) são considerados demograficamente maduros (Pitacco et al., 2009); este fundo está em trajetória ascendente, o que implica crescimento estrutural da sinistralidade mesmo sem piora intrínseca das taxas de mortalidade.

Pirâmide etária da exposição

Ocultar código 
intervalos <- c(0, 18, 40, 60, 80, 116)
labels_fx  <- c("[0, 18)", "[18, 40)", "[40, 60)", "[60, 80)", "80+")

sum_d2 <- sum_d2 |>
  mutate(faixa = cut(Idade, breaks = intervalos, right = FALSE, labels = labels_fx))

piramide_ind <- sum_d2 |>
  filter(Idade >= 18, Idade <= 100) |>
  mutate(
    faixa   = cut(Idade, breaks = intervalos, right = FALSE, labels = labels_fx),
    exp_neg = -tot_exp
  )

p_exp_ind <- ggplot(
  piramide_ind |>
    mutate(tooltip = sprintf("Idade: %d\nExpostos: %s",
                             Idade, formatC(tot_exp, format = "d", big.mark = ".", decimal.mark = ","))),
  aes(x = Idade, y = exp_neg, fill = faixa)
) +
  geom_col_interactive(aes(tooltip = tooltip, data_id = Idade),
                       width = 0.9, alpha = 0.85) +
  coord_flip() +
  scale_y_continuous(
    labels = function(x) scales::comma(abs(x), big.mark = "."),
    breaks = -seq(0, 20000, by = 5000)
  ) +
  scale_fill_scico_d(palette = "batlow", name = "Faixa") +
  labs(x = "Idade", y = "Expostos (pessoas-ano)") +
  theme_bw(base_size = 11) +
  theme(legend.position = "none")

p_obt_ind <- ggplot(
  piramide_ind |>
    mutate(tooltip = sprintf("Idade: %d\nÓbitos: %d", Idade, tot_obt)),
  aes(x = Idade, y = tot_obt, fill = faixa)
) +
  geom_col_interactive(aes(tooltip = tooltip, data_id = Idade),
                       width = 0.9, alpha = 0.85) +
  coord_flip() +
  scale_y_continuous(labels = scales::comma_format(big.mark = ".")) +
  scale_fill_scico_d(palette = "batlow", name = "Faixa") +
  labs(x = NULL, y = "Óbitos") +
  theme_bw(base_size = 11) +
  theme(legend.position = "bottom",
        axis.text.y  = element_blank(),
        axis.ticks.y = element_blank())

p_pir_static <- p_exp_ind + p_obt_ind +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

ggsave("index_files/figure-html/g-piramide-1.png",
       plot = p_pir_static, width = 10, height = 7, dpi = 150,
       create.dir = TRUE)

girafe(
  ggobj = p_pir_static,
  options = list(
    opts_hover(css = "opacity:1; stroke:white; stroke-width:0.5px;"),
    opts_hover_inv(css = "opacity:0.25;"),
    opts_tooltip(css = "background-color:white; color:black; border:1px solid #ccc; border-radius:4px; font-size:12px; padding:5px 8px;"),
    opts_toolbar(saveaspng = FALSE)
  )
)

Pirâmide etária: exposição (esq., valores absolutos) e óbitos (dir.) por idade individual — triênio 2012–2014. Cores por faixa de risco.

Ocultar código 
piramide <- sum_d2 |>
  group_by(faixa) |>
  summarise(
    Expostos = sum(tot_exp),
    Obitos   = sum(tot_obt),
    .groups  = "drop"
  ) |>
  mutate(pct_exp = Expostos / sum(Expostos) * 100,
         pct_obt = Obitos   / sum(Obitos)   * 100)

# Variáveis inline para o parágrafo
pct_exp_1840_4060 <- round(sum(piramide$pct_exp[piramide$faixa %in% c("[18, 40)", "[40, 60)")]), 1)
pct_obt_1840_4060 <- round(sum(piramide$pct_obt[piramide$faixa %in% c("[18, 40)", "[40, 60)")]), 1)
pct_exp_6080_80p  <- round(sum(piramide$pct_exp[piramide$faixa %in% c("[60, 80)", "80+")]), 1)
pct_obt_6080_80p  <- round(sum(piramide$pct_obt[piramide$faixa %in% c("[60, 80)", "80+")]), 1)

p_exp <- ggplot(piramide, aes(x = faixa, y = pct_exp, fill = faixa)) +
  geom_col(alpha = 0.85, width = 0.7) +
  geom_text(aes(label = sprintf("%.1f%%", pct_exp)), hjust = -0.1, size = 3.2) +
  coord_flip() +
  scale_y_continuous(limits = c(0, 52), labels = label_percent(scale = 1)) +
  scale_fill_scico_d(palette = "batlow") +
  labs(x = NULL, y = "% Expostos") +
  theme_bw(base_size = 11) + theme(legend.position = "none")

p_obt <- ggplot(piramide, aes(x = faixa, y = pct_obt, fill = faixa)) +
  geom_col(alpha = 0.85, width = 0.7) +
  geom_text(aes(label = sprintf("%.1f%%", pct_obt)), hjust = -0.1, size = 3.2) +
  coord_flip() +
  scale_y_continuous(limits = c(0, 55), labels = label_percent(scale = 1)) +
  scale_fill_scico_d(palette = "batlow") +
  labs(x = NULL, y = "% Óbitos") +
  theme_bw(base_size = 11) + theme(legend.position = "none",
                                    axis.text.y = element_blank(),
                                    axis.ticks.y = element_blank())

p_exp + p_obt +
  plot_annotation(title = "Distribuição relativa de expostos e óbitos por faixa etária")

Distribuição relativa de expostos e óbitos por faixa etária (triênio 2012–2014)

A pirâmide evidencia a assimetria estrutural do fundo: 73.4% dos expostos pertencem às faixas [18, 40) e [40, 60), mas apenas 10.3% dos óbitos. As faixas [60, 80) e 80+ reúnem 26.4% da exposição e 89.7% dos sinistros. Esse desacoplamento é o que justifica o modelo de credibilidade com pesos de volume distintos por faixa: a informação de cada grupo contribui de forma desigual para a estimativa do nível global \(\lambda_0\).

Distribuição etária da exposição e dos óbitos

Ocultar código 
p1 <- ggplot(fundo_long, aes(x = Idade, y = Expostos, colour = Ano)) +
  geom_line(linewidth = 0.9) +
  scale_x_continuous(breaks = seq(0, 110, by = 10)) +
  labs(x = "Idade", y = "Expostos", colour = NULL) +
  theme_bw(base_size = 11) + theme(legend.position = "bottom")

p2 <- ggplot(fundo_long, aes(x = Idade, y = Obitos, colour = Ano)) +
  geom_line(linewidth = 0.9) +
  scale_x_continuous(breaks = seq(0, 110, by = 10)) +
  labs(x = "Idade", y = "Óbitos", colour = NULL) +
  theme_bw(base_size = 11) + theme(legend.position = "bottom")

p1 + p2

Expostos e óbitos por idade e ano (2012–2014)
Ocultar código 
# Variáveis inline para o parágrafo
exp_19_trienio <- sum(fundo_long$Expostos[fundo_long$Idade == 19])
idade_max_obt  <- sum_d2$Idade[which.max(sum_d2$tot_obt)]
obt_max_trienio <- max(sum_d2$tot_obt)

A distribuição de expostos exibe dois picos: o principal entre 19 e 22 anos (entrada de jovens no plano, com idade 19 acumulando \(\approx 22\,000\) pessoas-ano no triênio) e um secundário entre 40 e 55 anos. Os óbitos concentram-se acima dos 60 anos, com máximo em torno de 84 anos (175 óbitos no triênio). A estabilidade temporal das curvas indica ausência de choques na composição da carteira no período.

Log da taxa de mortalidade e heatmap idade\(\times\)ano

Ocultar código 
ggplot(fundo_long |> filter(!is.na(log_tx)),
       aes(x = Idade, y = log_tx, colour = Ano)) +
  geom_line(linewidth = 0.9) +
  scale_x_continuous(breaks = seq(0, 110, by = 10)) +
  labs(x = "Idade", y = expression(log~hat(m)[x]), colour = NULL) +
  theme_bw(base_size = 12) + theme(legend.position = "bottom")

Log da taxa central de mortalidade por idade e ano (2012–2014)
Ocultar código 
heat_df <- fundo_long |>
  filter(!is.na(log_tx), Expostos >= 5) |>
  mutate(Ano = as.numeric(Ano))

p_heatmap <- ggplot(
  heat_df |> mutate(tooltip = sprintf("Idade: %d\nAno: %.0f\nlog(m\u0302\u2093) = %.2f",
                                      Idade, Ano, log_tx)),
  aes(x = Ano, y = Idade, fill = log_tx)
) +
  geom_tile_interactive(aes(tooltip = tooltip, data_id = paste(Idade, Ano)),
                        width = 0.95, height = 0.95) +
  scale_fill_scico(palette = "lajolla", direction = 1,
                   name = expression(log~hat(m)[x])) +
  scale_x_continuous(breaks = c(2012, 2013, 2014)) +
  scale_y_continuous(breaks = seq(0, 110, by = 10)) +
  labs(x = "Ano", y = "Idade") +
  theme_bw(base_size = 12)

ggsave("index_files/figure-html/g-heatmap-1.png",
       plot = p_heatmap, width = 8, height = 4.5, dpi = 150,
       create.dir = TRUE)

girafe(
  ggobj = p_heatmap,
  options = list(
    opts_hover(css = "stroke: white; stroke-width: 0.8px;"),
    opts_tooltip(css = "background-color:white; color:black; border:1px solid #ccc; border-radius:4px; font-size:12px; padding:5px 8px;"),
    opts_toolbar(saveaspng = FALSE)
  )
)

Heatmap de log(taxa central de mortalidade) por idade e ano — superfície de Lexis

O heatmap (superfície de Lexis discreta) apresenta o padrão típico de mortalidade humana: tons mais claros (taxas mais altas) nas idades acima de 70 anos e tons mais escuros (taxas baixas) nas idades jovens. A estrutura é consistente entre os três anos — não há coluna atípica que indique anomalia de registro em nenhum ano específico. Com apenas \(T = 3\) períodos, o heatmap serve para diagnóstico de consistência; a inferência de tendência temporal de melhoria de mortalidade requer horizontes mais longos (Rau et al., 2018).

A curva do log da taxa central exibe três segmentos: mínimo entre 5 e 15 anos (\(\log\hat{m}_x \approx -10\)), elevação suave entre 15 e 25 anos por causas externas, e crescimento aproximadamente linear acima dos 25 anos, compatível com a lei de Gompertz — \(m_x \propto e^{\beta x}\), com \(\beta \approx 0{,}09\) por ano. Nas idades acima de 90 anos as curvas anuais divergem, pois \(E_x \leq 5\) em muitas células: um único óbito adicional desloca \(\log\hat{m}_x\) em \(\log 2 \approx 0{,}7\) unidades, argumento direto para agregar essas idades num único grupo de risco.

Evolução dos óbitos anuais por faixa etária

Ocultar código 
if (!requireNamespace("gtExtras", quietly = TRUE))
  install.packages("gtExtras", repos = "https://cloud.r-project.org", quiet = TRUE)
library(gtExtras)

intervalos_q1 <- c(0, 18, 40, 60, 80, 116)
labels_fx_q1  <- c("[0, 18)", "[18, 40)", "[40, 60)", "[60, 80)", "80+")

box_df <- fundo_long |>
  mutate(faixa = cut(Idade, breaks = intervalos_q1, right = FALSE,
                     labels = labels_fx_q1)) |>
  filter(!is.na(faixa)) |>
  group_by(faixa, Ano) |>
  summarise(Obitos = sum(Obitos), .groups = "drop") |>
  mutate(faixa = factor(faixa, levels = labels_fx_q1))

spark_df <- box_df |>
  group_by(faixa) |>
  summarise(
    `2012`    = Obitos[Ano == "2012"],
    `2013`    = Obitos[Ano == "2013"],
    `2014`    = Obitos[Ano == "2014"],
    CV        = if_else(mean(Obitos) > 0,
                        round(sd(Obitos) / mean(Obitos) * 100, 1),
                        NA_real_),
    tendencia = list(Obitos),
    .groups   = "drop"
  )

spark_df |>
  gt() |>
  tab_header(
    title    = "Óbitos anuais por faixa etária (2012–2014)",
    subtitle = "A coluna Trajetória mostra a evolução 2012 → 2013 → 2014"
  ) |>
  cols_label(
    faixa     = "Faixa",
    `2012`    = "2012",
    `2013`    = "2013",
    `2014`    = "2014",
    CV        = "CV (%)",
    tendencia = "Trajetória"
  ) |>
  gt_plt_sparkline(
    tendencia,
    type       = "default",
    same_limit = FALSE,
    palette    = c("steelblue4", "steelblue4", "steelblue4", "steelblue4", "steelblue4")
  ) |>
  tab_style(
    style     = cell_fill(color = "white"),
    locations = cells_body(columns = tendencia)
  ) |>
  sub_missing(columns = CV, missing_text = "—") |>
  fmt_number(columns = c(`2012`, `2013`, `2014`), decimals = 0) |>
  tab_style(
    style     = cell_fill(color = "#FFF9C4"),
    locations = cells_body(rows = !is.na(CV) & CV > 15)
  ) |>
  tab_style(
    style     = list(cell_fill("#EBF5FB"), cell_text(weight = "bold")),
    locations = cells_body(rows = faixa == "80+")
  ) |>
  tab_style(
    style     = cell_text(color = "grey50", style = "italic"),
    locations = cells_body(rows = faixa == "[0, 18)")
  ) |>
  cols_align(align = "center", columns = c(`2012`, `2013`, `2014`, CV)) |>
  tab_options(
    table.font.size        = 13,
    table.width            = pct(80),
    column_labels.font.weight = "bold",
    heading.subtitle.font.size = 11
  )
Óbitos anuais por faixa etária (2012–2014)
A coluna Trajetória mostra a evolução 2012 → 2013 → 2014
Faixa 2012 2013 2014 CV (%) Trajetória
[0, 18) 0 0 0 — 0.00
[18, 40) 34 33 22 22.4 22.0
[40, 60) 97 92 99 3.8 99.0
[60, 80) 476 484 531 6.0 531.0
80+ 524 583 683 13.5 683.0

A tabela com sparklines comprime em uma linha por grupo toda a informação relevante: os três valores anuais, o CV e a trajetória visual. A faixa 80+ (destacada em azul) exibe crescimento monotônico de 524 para 583 para 683 óbitos, sinal do envelhecimento progressivo da carteira. A faixa \([18, 40)\) tem CV elevado (22.4%, destacado em amarelo) por oscilar entre 22 e 34 óbitos por ano. A faixa \([0, 18)\) registra zero óbitos nos três anos; sua taxa credibilizada virá inteiramente do nível coletivo \(\hat{\lambda}_0\).

Curva de sobrevivência acumulada

Ocultar código 
# lx construída a partir das mx do triênio acumulado
lx_df <- sum_d2 |>
  arrange(Idade) |>
  mutate(
    mx   = if_else(tot_exp > 0, tot_obt / tot_exp, 0),
    px   = exp(-mx),           # probabilidade de sobreviver 1 ano (força constante)
    lx   = cumprod(c(1, px[-n()])) * 100000   # radix = 100.000
  )

ggplot(lx_df |> filter(Idade <= 110, lx > 0), aes(x = Idade, y = lx)) +
  geom_line(linewidth = 1.0, colour = "steelblue4") +
  geom_area(alpha = 0.15, fill = "steelblue") +
  geom_vline(xintercept = c(18, 40, 60, 80),
             linetype = "dotted", colour = "grey40", linewidth = 0.7) +
  scale_x_continuous(breaks = seq(0, 110, by = 10)) +
  scale_y_continuous(labels = scales::comma_format(big.mark = ".")) +
  labs(x = "Idade", y = expression(italic(l)[x]~"(radix = 100.000)")) +
  theme_bw(base_size = 12)

Função de sobrevivência acumulada \(l_x\) a partir das taxas centrais observadas (triênio 2012–2014)

A curva \(l_x\) parte de 100.000 nascidos vivos e decresce de acordo com as taxas centrais observadas no triênio, sob a hipótese de força de mortalidade constante dentro de cada ano de idade (\(p_x = e^{-\hat{m}_x}\)). A queda mais acentuada ocorre acima dos 60 anos, com a curva atingindo \(\approx 64.100\) aos 80 anos — consistente com o índice de dependência de 0.359 calculado anteriormente. As linhas pontilhadas marcam as fronteiras das faixas etárias definidas na questão (ii).

Questão (ii): Perfil Etário e Faixas de Risco

Critérios de formação dos grupos

A definição dos \(m\) grupos para o modelo de Bühlmann-Straub obedece a três critérios simultâneos, seguindo Pitacco et al. (2009) e Bühlmann; Gisler (2005):

Homogeneidade interna. As taxas individuais dentro de cada faixa devem variar suavemente para que a taxa agregada seja representativa. A linearidade em log nas faixas adultas indica que intervalos de 20 anos mantêm homogeneidade razoável.

Exposição suficiente. O modelo pondera cada observação pelo volume \(v_{it} = E_{it}\) (pessoas-ano). Com \(T = 2\) anos de ajuste e \(m = 4\) grupos (excluindo \([0, 18)\) por ausência de sinistros), os parâmetros estruturais \(\hat{\lambda}_0\) e \(\hat{\tau}^2\) são estimados com \(m - 1 = 3\) graus de liberdade para a componente entre grupos. Grupos com \(v_{it} \approx 0\) recebem peso negligenciável automaticamente.

Relevância atuarial. Os grupos devem capturar os regimes de risco que impactam o passivo do plano: jovens ativos, núcleo produtivo, aposentados recentes e idosos avançados.

A segmentação inicial em \(\{18, 40, 60, 80\}\) parte do julgamento atuarial; a seção seguinte a valida objetivamente por detecção de pontos de mudança via PELT.

Faixas etárias e estatísticas descritivas

Ocultar código 
sum_faixas <- sum_d2 |>
  group_by(faixa) |>
  summarise(
    tot_obt = sum(tot_obt),
    tot_exp = sum(tot_exp),
    .groups = "drop"
  ) |>
  mutate(
    tx_media = tot_obt / tot_exp,
    pct_exp  = tot_exp / sum(tot_exp),
    pct_obt  = tot_obt / sum(tot_obt)
  )

# Coeficiente de variação dos óbitos anuais por faixa
cv_data <- tibble(
  faixa  = labels_fx,
  obt_12 = c(0, 34, 97,  476, 524),
  obt_13 = c(0, 33, 92,  484, 583),
  obt_14 = c(0, 22, 99,  531, 683)
) |>
  rowwise() |>
  mutate(
    media = mean(c(obt_12, obt_13, obt_14)),
    cv    = if_else(media > 0,
                    sd(c(obt_12, obt_13, obt_14)) / media, NA_real_)
  )

sum_faixas |>
  left_join(cv_data |> select(faixa, cv), by = "faixa") |>
  gt() |>
  tab_header(title = "Resumo por faixa etária (triênio 2012--2014)") |>
  cols_label(
    faixa    = "Faixa",
    tot_obt  = "Óbitos",
    tot_exp  = "Expostos",
    tx_media = md("$\\hat{m}$ média"),
    pct_exp  = "% Expostos",
    pct_obt  = "% Óbitos",
    cv       = "CV óbitos anuais"
  ) |>
  fmt_number(columns = c(tot_obt, tot_exp), decimals = 0) |>
  fmt_number(columns = tx_media, decimals = 6) |>
  fmt_percent(columns = c(pct_exp, pct_obt), decimals = 1) |>
  fmt_number(columns = cv, decimals = 3) |>
  tab_options(table.font.size = 13)
Resumo por faixa etária (triênio 2012--2014)
Faixa Óbitos Expostos \(\hat{m}\) média % Expostos % Óbitos CV óbitos anuais
[0, 18) 0 1,029 0.000000 0.2% 0.0% NA
[18, 40) 89 180,635 0.000493 43.2% 2.4% 0.224
[40, 60) 288 126,107 0.002284 30.2% 7.9% 0.038
[60, 80) 1,491 84,038 0.017742 20.1% 40.8% 0.060
80+ 1,790 26,155 0.068438 6.3% 48.9% 0.135
Ocultar código 
# Variáveis inline dos CVs
cv_4060 <- round(cv_data$cv[cv_data$faixa == "[40, 60)"] * 100, 1)
cv_1840 <- round(cv_data$cv[cv_data$faixa == "[18, 40)"] * 100, 1)
cv_80p  <- round(cv_data$cv[cv_data$faixa == "80+"]      * 100, 1)

O coeficiente de variação (CV) dos óbitos anuais quantifica a estabilidade temporal de cada grupo. A faixa \([40, 60)\) tem CV \(= 3.8\%\), a mais estável, enquanto \([18, 40)\) atinge 22.4%, reflexo de contagens baixas e raras em que um ou dois óbitos a mais ou a menos representam variação percentual expressiva. A faixa \(80+\) tem CV \(= 13.5\%\) por exposição reduzida em idades muito avançadas. No modelo de Bühlmann-Straub, esses grupos de alta variabilidade recebem menor peso \(v_{it}\) e, consequentemente, menor credibilidade \(\omega_i\) em favor do nível coletivo \(\lambda_0\).

Ocultar código 
bxp_df <- sum_d2 |>
  filter(tot_exp > 0, tot_obt > 0, Idade >= 18) |>
  mutate(
    mx    = tot_obt / tot_exp,
    faixa = cut(Idade, breaks = intervalos, right = FALSE, labels = labels_fx)
  ) |>
  filter(!is.na(faixa), faixa != "[0, 18)")

p_bxp <- ggplot(
  bxp_df |> mutate(tooltip = sprintf("Idade: %d\nm̂ₓ = %.5f", Idade, mx)),
  aes(x = faixa, y = mx, fill = faixa, colour = faixa)
) +
  geom_boxplot(alpha = 0.35, outlier.shape = NA, linewidth = 0.7) +
  geom_point_interactive(aes(tooltip = tooltip, data_id = Idade),
                         position = position_jitter(width = 0.15, seed = 42),
                         size = 1.4, alpha = 0.6) +
  scale_y_log10(
    breaks = 10^(-5:-1),
    labels = scales::trans_format("log10", scales::math_format(10^.x))
  ) +
  scale_fill_scico_d(palette   = "batlow", name = "Faixa") +
  scale_colour_scico_d(palette = "batlow", name = "Faixa") +
  labs(x = "Faixa etária",
       y = expression(hat(m)[x]~"(escala log)")) +
  theme_bw(base_size = 12) +
  theme(legend.position = "none")

ggsave("index_files/figure-html/g-boxplot-faixas-1.png",
       plot = p_bxp, width = 8, height = 4, dpi = 150,
       create.dir = TRUE)

girafe(
  ggobj = p_bxp,
  options = list(
    opts_hover(css = "opacity:1; r:4px;"),
    opts_hover_inv(css = "opacity:0.2;"),
    opts_tooltip(css = "background-color:white; color:black; border:1px solid #ccc; border-radius:4px; font-size:12px; padding:5px 8px;"),
    opts_toolbar(saveaspng = FALSE)
  )
)

Distribuição das taxas centrais de mortalidade \(\hat{m}_x\) por faixa etária (escala log). Cada ponto é uma idade individual no triênio agregado. A variabilidade intragrupo é um insumo direto de \(\hat{\sigma}^2\) no modelo de Bühlmann-Straub.

O boxplot revela três padrões relevantes para a formação dos grupos. Primeiro, a separação entre medianas é pronunciada — cada faixa ocupa uma região distinta da escala log, validando a homogeneidade interna das fronteiras adotadas. Segundo, a dispersão intragrupo decresce de \([18, 40)\) para \([40, 60)\) e volta a crescer em \(80+\): \([18, 40)\) tem alta dispersão por incluir o accident hump (idades 18–30) e o início do Gompertz (31–39); \([40, 60)\) é o grupo mais homogêneo, com mortalidade crescendo suavemente; \([60, 80)\) apresenta dispersão moderada; e \(80+\) exibe dispersão crescente por esparsidade nas idades extremas (\(D_x \in \{0, 1, 2\}\) acima de 90 anos). Essa dispersão intragrupo é exatamente o que o estimador \(\hat{\sigma}^2 = \hat{\lambda}_0\) captura no modelo de Bühlmann-Straub — grupos mais heterogêneos recebem menor \(\hat{\omega}_i\).

Curva de mortalidade e detecção de pontos de mudança (PELT Poisson)

A escolha de pontos de corte por julgamento atuarial introduz subjetividade. Uma alternativa estatisticamente fundamentada é formular a segmentação etária como um problema de particionamento ótimo de série ordenada: dada a série de pares \((D_x, E_x)\) para \(x = x_1 < \cdots < x_n\), encontrar os \(k\) pontos de corte \(\tau_0 = 0 < \tau_1 < \cdots < \tau_k < \tau_{k+1} = n\) que minimizam

\[ \mathcal{Q}(k, \boldsymbol{\tau}) = \sum_{j=1}^{k+1} \underbrace{\mathcal{C}\!\bigl((D_x, E_x)_{\tau_{j-1}+1:\tau_j}\bigr)}_{\text{custo intragrupo}} + \beta\,(k+1), \tag{PELT} \]

onde \(\beta > 0\) é uma penalidade que controla o número de grupos.

Função de custo Poisson. A abordagem anterior usava o modelo Normal sobre \(\log\hat{m}_x\), o que tem uma limitação estrutural: pelo método delta, \(\text{Var}[\log\hat{m}_x] \approx 1/D_x\), de modo que uma célula com \(D_x = 2\) tem variância \(\approx 0{,}5\) e uma com \(D_x = 200\) tem variância \(\approx 0{,}005\) — diferença de 100 vezes que o modelo Normal homoscedástico ignora, dando peso igual a células esparssas e densas. Adotamos em seu lugar o custo Poisson, que é a log-verossimilhança negativa do modelo de taxa constante dentro de cada segmento:

\[ \mathcal{C}([a, b]) = -2\left[ D_{a:b}\log\!\frac{D_{a:b}}{E_{a:b}} - D_{a:b} \right], \qquad D_{a:b} = \sum_{x=a}^{b} D_x, \quad E_{a:b} = \sum_{x=a}^{b} E_x. \]

Este custo tem três propriedades fundamentais para este problema: (i) pondera automaticamente cada idade pela sua informação, pois células com \(D_x\) pequeno contribuem com verossimilhança difusa e, portanto, menos influência; (ii) não requer a transformação \(\log\hat{m}_x\), eliminando problemas com \(D_x = 0\) e o viés de Jensen em células esparsas; (iii) é coerente com a hipótese \(N_x \sim \text{Poisson}(E_x m_x)\) do próprio modelo de Bühlmann-Straub.

A penalidade BIC é \(\beta = \log(n)\) (um parâmetro por segmento: a taxa \(m_j = D_{a:b}/E_{a:b}\)). O algoritmo PELT (Pruned Exact Linear Time, Killick; Fearnhead; Eckley (2012)) resolve o problema de forma exata em tempo \(O(n)\) amortizado via programação dinâmica com poda: um candidato \(s\) é descartado quando

\[ \mathcal{Q}(s) + \mathcal{C}(\mathbf{y}_{s+1:t}) + K \geq \mathcal{Q}(t) \]

para algum \(t > s\), garantindo que \(s\) não pode integrar nenhuma solução ótima futura. A qualidade da segmentação é avaliada pela razão de taxas entre segmentos adjacentes: \(\hat{m}_{j+1}/\hat{m}_j\), onde \(\hat{m}_j = D_{a:b}/E_{a:b}\) é a taxa agregada do segmento \(j\). Razões \(\gg 1\) confirmam separação genuína de regimes de mortalidade; razões próximas de 1 indicam subdivisão sem ganho informacional para o modelo de credibilidade.

Ocultar código 
if (!requireNamespace("fastcpd", quietly = TRUE))
  install.packages("fastcpd", repos = "https://cloud.r-project.org", quiet = TRUE)
library(fastcpd)

# Dados: pares (D_x, E_x) para idades com exposição positiva >= 18
dados_cpt <- sum_d2 |>
  filter(tot_exp > 0, Idade >= 18) |>
  arrange(Idade) |>
  mutate(
    log_mx = if_else(tot_obt > 0, log(tot_obt / tot_exp), NA_real_)
  )

# ── Custo Poisson analítico ───────────────────────────────────────────────────
# C([a,b]) = -2 * [D_{a:b} * log(D_{a:b}/E_{a:b}) - D_{a:b}]
# Retorna NA se D_{a:b} = 0 (segmento sem óbitos → custo infinito → nunca ótimo)
poisson_cost <- function(data) {
  # data é uma submatriz com colunas [D_x, E_x]
  D <- sum(data[, 1])
  E <- sum(data[, 2])
  if (D == 0 || E == 0) return(NA_real_)
  -2 * (D * log(D / E) - D)
}

# Matriz de entrada: [óbitos, exposição]
mat_pelt <- dados_cpt |> select(tot_obt, tot_exp) |> as.matrix()
n_obs    <- nrow(mat_pelt)

# PELT com custo Poisson customizado, penalidade BIC (log(n) por segmento)
resultado_pelt <- fastcpd(
  formula   = ~ . - 1,
  data      = as.data.frame(mat_pelt),
  cost      = poisson_cost,
  beta      = log(n_obs),       # penalidade BIC
  segment_count = 10,           # máximo de segmentos candidatos
  min_segment_length = 8,       # mínimo 8 idades por segmento
  r.progress = FALSE
)

idx_corte    <- resultado_pelt@cp_set
idades_cpt   <- dados_cpt$Idade
idades_corte <- idades_cpt[idx_corte]

# Segmento PELT e faixa adotada para cada idade
dados_cpt <- dados_cpt |>
  mutate(
    seg_pelt = cut(Idade,
      breaks = c(min(Idade) - 1, idades_corte, max(Idade)),
      labels = FALSE),
    faixa_adotada = cut(Idade,
      breaks = c(18, 40, 60, 80, 116), right = FALSE,
      labels = c("[18, 40)", "[40, 60)", "[60, 80)", "80+"))
  )

# ── Qualidade: razão de deviance (redução de custo relativa) ─────────────────
# Nota: o R² Poisson (1 - custo_seg/custo_nulo) não é análogo ao R² Normal
# porque os custos têm escalas distintas. Usamos a razão de taxas entre
# segmentos adjacentes como métrica de separação.
stats_pelt <- dados_cpt |>
  filter(!is.na(seg_pelt)) |>
  group_by(seg_pelt) |>
  summarise(
    faixa    = sprintf("[%d, %d]", min(Idade), max(Idade)),
    n_idades = n(),
    D        = sum(tot_obt),
    E        = round(sum(tot_exp)),
    mx_grp   = round(sum(tot_obt) / sum(tot_exp), 5),
    .groups  = "drop"
  )

# Razão de taxas entre segmentos adjacentes (medida de separação)
razao_adj <- round(stats_pelt$mx_grp[-1] / stats_pelt$mx_grp[-nrow(stats_pelt)], 2)

# ── Output ────────────────────────────────────────────────────────────────────
cat("Pontos de mudança PELT-Poisson:", paste(idades_corte, collapse = ", "), "\n")
cat(sprintf("N grupos PELT: %d\n\n", length(idades_corte) + 1))

cat("=== Estatísticas por segmento PELT-Poisson ===\n")
print(stats_pelt)
cat("\nRazão de taxas entre segmentos adjacentes:",
    paste(razao_adj, collapse = " | "), "\n")

cat("\n=== Estatísticas por faixa adotada {18,40,60,80} ===\n")
print(dados_cpt |>
  filter(!is.na(faixa_adotada)) |>
  group_by(faixa_adotada) |>
  summarise(
    n_idades = n(),
    D        = sum(tot_obt),
    E        = round(sum(tot_exp)),
    mx_grp   = round(sum(tot_obt) / sum(tot_exp), 5),
    .groups  = "drop"
  ))
Pontos de mudança PELT-Poisson: 41, 58, 78, 91, 97 
N grupos PELT: 6

=== Estatísticas por segmento PELT-Poisson ===
# A tibble: 6 × 6
  seg_pelt faixa     n_idades     D      E mx_grp
     <int> <chr>        <int> <int>  <dbl>  <dbl>
1        1 [18, 41]        24    95 189987 0.0005
2        2 [42, 58]        17   257 111928 0.0023
3        3 [59, 78]        20  1398  86082 0.0162
4        4 [79, 91]        13  1650  26401 0.0625
5        5 [92, 97]         6   226   1938 0.117 
6        6 [98, 114]       17    32    599 0.0534

Razão de taxas entre segmentos adjacentes: 4.6 | 7.06 | 3.85 | 1.87 | 0.46 

=== Estatísticas por faixa adotada {18,40,60,80} ===
# A tibble: 4 × 5
  faixa_adotada n_idades     D      E  mx_grp
  <fct>            <int> <int>  <dbl>   <dbl>
1 [18, 40)            22    89 180635 0.00049
2 [40, 60)            20   288 126107 0.00228
3 [60, 80)            20  1491  84038 0.0177 
4 80+                 35  1790  26155 0.0684
Ocultar código 
p_faixas_pelt <- ggplot(
  dados_cpt |>
    filter(!is.na(log_mx)) |>
    mutate(
      tx_tot  = tot_obt / tot_exp,
      tooltip = sprintf(
        "Faixa: %s\nIdade: %d\nm\u0302\u2093 = %.5f\nSegmento PELT: %d",
        faixa_adotada, Idade, tx_tot, seg_pelt)
    ),
  aes(x = Idade, y = tot_obt / tot_exp, colour = faixa_adotada)
) +
  geom_point_interactive(aes(tooltip = tooltip, data_id = faixa_adotada),
                         size = 1.8) +
  geom_vline(xintercept = idades_corte,
             linetype = "solid", colour = "grey30", linewidth = 0.85) +
  scale_colour_scico_d(palette = "batlow", name = "Faixa") +
  scale_y_log10(
    breaks = 10^(-5:-1),
    labels = scales::trans_format("log10", scales::math_format(10^.x))
  ) +
  scale_x_continuous(breaks = seq(20, 100, by = 10)) +
  labs(x = "Idade", y = expression(hat(m)[x]~"(escala log)")) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

ggsave("index_files/figure-html/g-faixas-pelt-1.png",
       plot = p_faixas_pelt, width = 9, height = 4.2, dpi = 150,
       create.dir = TRUE)

girafe(
  ggobj = p_faixas_pelt,
  options = list(
    opts_hover(css = "r:5px;"),
    opts_hover_inv(css = "opacity:0.2;"),
    opts_tooltip(css = "background-color:white; color:black; border:1px solid #ccc; border-radius:4px; font-size:12px; padding:5px 8px;"),
    opts_toolbar(saveaspng = FALSE)
  )
)

Curva \(\hat{m}_x\) por faixa adotada (escala log) com pontos de mudança PELT-Poisson (linhas verticais cinzas).

O PELT-Poisson identificou 6 segmentos com pontos de mudança em \(x \in \{41, 58, 78, 91, 97\}\). A métrica relevante não é um \(R^2\) global — cujo valor é determinado pela escala do custo Poisson e não tem interpretação direta — mas as razões de taxa entre segmentos adjacentes, que medem a separação real entre grupos: \(4.6\times, 7.06\times, 3.85\times, 1.87\times, 0.46 \times\) (da menor para a maior taxa). Razões \(\gg 1\) confirmam quebras genuínas de regime; razões próximas de 1 indicam subdivisão sem ganho informacional.

O PELT-Poisson opera diretamente nos pares \((D_x, E_x)\) sem transformação logarítmica, eliminando a influência desigual de células esparsas que afetava a versão Normal. Os cortes em 41 e 58 capturam a transição do accident hump para o regime Gompertz e a maior quebra de regime da série (razão 7×, a mais pronunciada). Os cortes em 78 e 91 isolam a zona de alta mortalidade e o início da zona esparsa acima de 90 anos — a inversão de taxa em \([98, 114]\) (\(\hat{m} = 0.053 < \hat{m}_{[92,97]} = 0.117\)) com apenas \(D = 32\) óbitos é artefato de esparsidade, não desaceleração biológica.

Sobre a fronteira em 58 vs. 60. O PELT identifica 58 como a maior quebra individual da série. Uma segmentação alternativa \(\{18, 40, 58, 80\}\) poderia ser considerada, mas apresenta um custo metodológico relevante com \(T = 2\) anos: o grupo \([58, 80)\) incorporaria as idades 58–59, que têm taxa de mortalidade \(\approx 0{,}007\) — substancialmente abaixo da taxa do grupo \([60, 80)\) atual (\(\hat{m} = 0{,}0177\)). Isso reduziria a taxa agregada de \([58, 80)\) em aproximadamente 20% em relação a \([60, 80)\), alterando \(\hat{\lambda}_0\) e os prêmios de credibilidade sem ganho proporcional em identificação de \(\hat{\tau}^2\) — que com \(T = 2\) já opera com apenas \(m - 1 = 3\) graus de liberdade. Com horizonte \(T \geq 5\) anos, a fronteira em 58 seria preferível; com \(T = 2\), a diferença de duas idades não justifica a instabilidade adicional.

Adotamos portanto \(m = 4\) grupos com fronteiras \(\{18, 40, 60, 80\}\). Os grupos terminais do PELT ([92,97] com \(D = 226\) e [98,114] com \(D = 32\)) teriam volumes \(v_{it}\) insuficientes para credibilidades distintas — estão corretamente absorvidos no grupo \(80+\) como zona terminal de baixa credibilidade (\(\hat{\omega}_{80+} \approx 1\) por \(v_{80+} \gg \hat{\kappa}\), mas com alta incerteza nas estimativas individuais acima de 90 anos). O PELT-Poisson cumpre seu papel de validação objetiva: os cortes em 41, 58 e 78 confirmam que as três principais quebras de regime desta carteira ocorrem nas regiões \(\{40, 60, 80\}\) adotadas, com margem de \(\pm 2\) anos.

Questão (iii): Seleção de Tábuas de Referência

Contexto regulatório e das tábuas

A escolha da tábua de referência é relevante para o modelo de Bühlmann-Straub porque ela define o nível coletivo \(\hat{\lambda}_0\) ao qual as estimativas individuais de cada grupo são encolhidas. Uma tábua sistematicamente inadequada introduz viés em todos os prêmios de credibilidade (Bühlmann; Gisler, 2005).

No Brasil, a referência regulatória para EFPCs é estabelecida pela Resolução CNPC nº 30/2018 e pela Portaria PREVIC nº 835/2020 (PREVIC, 2020), que determinam a realização de pelo menos dois testes estatísticos de aderência. Avaliamos:

  • AT-2000 masculina (Society of Actuaries, 2000): tábua americana de rentistas, adotada historicamente no Brasil por ausência de alternativa nacional. A principal limitação é o viés de tradução demográfica: a mortalidade americana de rentistas de 2000 difere substancialmente da mortalidade brasileira atual. Adicionalmente, trata-se de uma tábua de seguros de vida individuais, população distinta de participantes de fundos de pensão.
  • BR-EMSsb-m v.2021 (FenaPrevi, 2021): tábua construída com dados de apólices do mercado segurador brasileiro, publicada pela FenaPrevi em 2021. Representa a experiência nacional mais recente e é a referência regulatória atual.
  • Pub-2010 Amount-Weighted Retiree masculina (Society of Actuaries, 2014): tábua desenvolvida pela SOA especificamente para participantes de fundos de pensão públicos americanos, com base em experiência de 2010. É metodologicamente mais próxima do perfil de uma EFPC do que a AT-2000 — ambas são tábuas americanas, mas a Pub-2010 foi calibrada para aposentados de planos de benefício definido, população com perfil demográfico similar ao desta carteira. A comparação permite avaliar em que medida a adoção histórica da AT-2000 no Brasil foi adequada mesmo no contexto internacional.
Ocultar código 
at2000 <- read.csv("at2000.iba.csv",    header = TRUE, dec = ",")
brems  <- read.csv("brems2021.iba.csv", header = TRUE, dec = ",")

# ── Pub-2010 Amount-Weighted Retiree (masculina) via mortSOA ─────────────────
if (!requireNamespace("mortSOA", quietly = TRUE))
  install.packages("mortSOA", repos = "https://cloud.r-project.org", quiet = TRUE)
library(mortSOA)

# ID 3400: Pub-2010 Amount-Weighted General Retiree — Male (SOA, 2014)
pub2010_raw <- read_mort_soa(3400)
pub2010_m   <- pub2010_raw[[1]] |>        # tibble com age e qx
  rename(idade = age) |>
  filter(idade <= 115) |>
  arrange(idade)

# Completar idades 0-115 com qx = NA para idades não cobertas pela tábua
# (Pub-2010 começa em 18; idades abaixo recebem qx = NA)
pub2010_full <- tibble(idade = 0:115) |>
  left_join(pub2010_m, by = "idade")

# Os dados fornecem a taxa central mx = Dx/Ex. As tábuas de referência estão
# em qx. Convertemos: qx_obs = 1 - exp(-mx), para que a comparação e a razão
# A/E sejam conceitualmente coerentes (observado e esperado na mesma escala).
frame_comp <- sum_d2 |>
  select(Idade, tot_obt, tot_exp) |>
  mutate(
    mx_obs   = if_else(tot_exp > 0, tot_obt / tot_exp, NA_real_),
    qx_obs   = if_else(!is.na(mx_obs), 1 - exp(-mx_obs), NA_real_),
    at2000_m = at2000$at2000.m,
    brems_m  = brems$brems2021sb.m[1:116],
    pub2010_m = pub2010_full$qx,
    obt_at   = tot_exp * at2000_m,
    obt_br   = tot_exp * brems_m,
    obt_pub  = tot_exp * pub2010_m
  )

Comparação visual das tábuas

Ocultar código 
ggplot(frame_comp |> filter(!is.na(qx_obs), qx_obs > 0, Idade <= 110), aes(x = Idade)) +
  geom_point(aes(y = qx_obs,    colour = "Observado"),          size = 1.2, alpha = 0.55) +
  geom_line(aes(y = at2000_m,   colour = "AT-2000-m"),          linewidth = 0.95) +
  geom_line(aes(y = brems_m,    colour = "BR-EMSsb-m v.2021"),  linewidth = 0.95) +
  geom_line(aes(y = pub2010_m,  colour = "Pub-2010 Retiree-m"), linewidth = 0.95,
            linetype = "dashed", na.rm = TRUE) +
  scale_colour_manual(
    values = c(
      "Observado"           = "grey40",
      "AT-2000-m"           = "#E74C3C",
      "BR-EMSsb-m v.2021"  = "#2980B9",
      "Pub-2010 Retiree-m" = "#27AE60"
    ),
    name = NULL
  ) +
  scale_y_log10(
    breaks = 10^(-5:-1),
    labels = scales::trans_format("log10", scales::math_format(10^.x))
  ) +
  scale_x_continuous(breaks = seq(0, 110, by = 10)) +
  labs(x = "Idade", y = expression(q[x]~"(escala log)")) +
  theme_bw(base_size = 12) + theme(legend.position = "bottom")

Comparação das tábuas com a experiência observada (escala log; observado convertido para \(\hat{q}_x = 1 - e^{-\hat{m}_x}\))

Razão A/E por faixa etária

Os dados do fundo fornecem diretamente a taxa central de mortalidade \(\hat{m}_x = D_x / E_x\), enquanto as tábuas AT-2000 e BR-EMSsb-m estão calibradas em probabilidade de morte \(q_x\). Para que a comparação e a razão A/E sejam conceitualmente coerentes (observado e esperado na mesma escala), convertemos a taxa central observada por meio da relação de força constante:

\[ \hat{q}_x = 1 - e^{-\hat{m}_x}. \]

A razão A/E por faixa agrega os esperados ponderados pela exposição:

\[ \text{A/E}_j = \frac{\sum_{x \in j} D_x}{\sum_{x \in j} E_x \cdot q_x^{\text{tábua}}} \]

Ocultar código 
frame_ae <- frame_comp |>
  mutate(faixa = cut(Idade, breaks = intervalos, right = FALSE, labels = labels_fx)) |>
  group_by(faixa) |>
  summarise(
    obt_obs = sum(tot_obt),
    esp_at  = sum(obt_at,  na.rm = TRUE),
    esp_br  = sum(obt_br,  na.rm = TRUE),
    esp_pub = sum(obt_pub, na.rm = TRUE),
    .groups = "drop"
  ) |>
  mutate(
    ae_at  = obt_obs / esp_at,
    ae_br  = obt_obs / esp_br,
    ae_pub = obt_obs / esp_pub
  )

frame_ae |>
  gt() |>
  tab_header(title = "Razão A/E por faixa etária (triênio 2012--2014)") |>
  cols_label(
    faixa   = "Faixa",    obt_obs = "Óbitos obs.",
    esp_at  = "Esp. AT-2000",       esp_br  = "Esp. BR-EMS",
    esp_pub = "Esp. Pub-2010",
    ae_at   = "A/E AT-2000",        ae_br   = "A/E BR-EMS",
    ae_pub  = "A/E Pub-2010"
  ) |>
  fmt_number(columns = c(obt_obs, esp_at, esp_br, esp_pub), decimals = 1) |>
  fmt_number(columns = c(ae_at, ae_br, ae_pub), decimals = 3) |>
  # Verde na menor distância de 1 entre as três tábuas
  tab_style(
    style     = cell_fill(color = "#D5F5E3"),
    locations = cells_body(
      columns = ae_br,
      rows    = abs(ae_br - 1) < abs(ae_at - 1) & abs(ae_br - 1) <= abs(ae_pub - 1)
    )
  ) |>
  tab_style(
    style     = cell_fill(color = "#D5F5E3"),
    locations = cells_body(
      columns = ae_pub,
      rows    = !is.na(ae_pub) &
                abs(ae_pub - 1) < abs(ae_at - 1) & abs(ae_pub - 1) < abs(ae_br - 1)
    )
  ) |>
  tab_style(
    style     = cell_fill(color = "#D5F5E3"),
    locations = cells_body(
      columns = ae_at,
      rows    = abs(ae_at - 1) < abs(ae_br - 1) &
                (is.na(ae_pub) | abs(ae_at - 1) < abs(ae_pub - 1))
    )
  ) |>
  tab_options(table.font.size = 13)
Razão A/E por faixa etária (triênio 2012--2014)
Faixa Óbitos obs. Esp. AT-2000 Esp. BR-EMS Esp. Pub-2010 A/E AT-2000 A/E BR-EMS A/E Pub-2010
[0, 18) 0.0 0.5 0.5 0.0 0.000 0.000 NaN
[18, 40) 89.0 119.9 150.9 0.0 0.742 0.590 Inf
[40, 60) 288.0 441.5 362.6 284.5 0.652 0.794 1.012
[60, 80) 1,491.0 1,638.7 1,445.1 1,389.3 0.910 1.032 1.073
80+ 1,790.0 2,361.7 2,160.1 2,589.6 0.758 0.829 0.691
Ocultar código 
ae_long <- frame_ae |>
  filter(faixa != "[0, 18)") |>
  select(faixa, ae_at, ae_br, ae_pub) |>
  pivot_longer(c(ae_at, ae_br, ae_pub), names_to = "tabua", values_to = "ae") |>
  filter(!is.na(ae)) |>
  mutate(
    tabua   = recode(tabua,
                     ae_at  = "AT-2000-m",
                     ae_br  = "BR-EMSsb-m v.2021",
                     ae_pub = "Pub-2010 Retiree-m"),
    tooltip = sprintf("Faixa: %s\nTábua: %s\nA/E = %.3f", faixa, tabua, ae)
  )

p_ae <- ggplot(ae_long, aes(x = faixa, y = ae, colour = tabua, group = tabua)) +
  geom_line_interactive(linewidth = 0.9) +
  geom_point_interactive(aes(tooltip = tooltip, data_id = paste(faixa, tabua)),
                         size = 3.5) +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "grey40") +
  scale_colour_manual(
    values = c(
      "AT-2000-m"           = "#E74C3C",
      "BR-EMSsb-m v.2021"  = "#2980B9",
      "Pub-2010 Retiree-m" = "#27AE60"
    ),
    name = NULL
  ) +
  labs(x = "Faixa etária", y = "A/E") +
  theme_bw(base_size = 12) + theme(legend.position = "bottom")

ggsave("index_files/figure-html/g-ae-1.png",
       plot = p_ae, width = 8, height = 3.6, dpi = 150,
       create.dir = TRUE)

girafe(
  ggobj = p_ae,
  options = list(
    opts_hover(css = "r:6px;"),
    opts_hover_inv(css = "opacity:0.25;"),
    opts_tooltip(css = "background-color:white; color:black; border:1px solid #ccc; border-radius:4px; font-size:12px; padding:5px 8px;"),
    opts_toolbar(saveaspng = FALSE)
  )
)

Razão A/E por faixa etária

A AT-2000 superestima a mortalidade em todas as faixas (A/E \(< 1\)). A BR-EMSsb-m v.2021 é mais aderente: apresenta A/E próximo de 1 nas faixas intermediárias, subestima ligeiramente em [60, 80) (A/E \(= 1.032\)) e superestima em \(80+\) (A/E \(= 0.829\)). A Pub-2010 Retiree apresenta A/E próximo de 1.073 em \([60, 80)\) e 0.691 em \(80+\) — valores mais próximos de 1 do que a AT-2000, confirmando que uma tábua de aposentados americanos de fundos de pensão é mais comparável a esta carteira do que a AT-2000 de seguros de vida individuais. Ambas as tábuas americanas, contudo, refletem a mortalidade de uma população estrangeira e são inferiores à BR-EMS em aderência global.

Testes estatísticos de aderência

A comparação visual e a razão A/E fornecem diagnóstico qualitativo; para a exigência regulatória da Portaria PREVIC nº 835/2020 (PREVIC, 2020), são necessários pelo menos dois testes formais. Utilizamos o pacote mortalityAdherence de Melo; Graziadei; Targino (2026), que implementa 11 testes em quatro famílias:

  • Não-paramétricos (KS, Qui-quadrado): testam se a distribuição empírica dos óbitos por idade é compatível com a esperada pela tábua, sem hipótese paramétrica para o processo de contagem.
  • Poisson Tipo I (Wald e LRT): testam \(H_0: \text{A/E global} = 1\), com \(D_x \sim \text{Poisson}(E_x \cdot q_x)\) sob a tábua.
  • Poisson Tipo II (Wald e LRT): testam ausência de viés etário, ou seja, se a razão observado/esperado é constante ao longo das idades.
  • Binomial Negativo (Tipo I e II): versões robustas à superdispersão. Rejeição nos testes NB confirma que o desajuste é real e não artefato da hipótese Poisson.
  • Intervalo bayesiano (BayesCI): verifica se a A/E observada está dentro do intervalo credível a \(95\%\) sob priori de Jeffreys para \(\theta = q_{\text{obs}}/q_{\text{tábua}}\).
Ocultar código 
fund_adh <- sum_d2 |> select(age = Idade, exposed = tot_exp, deaths = tot_obt)
tab_at   <- data.frame(age = at2000$idade,       qx = at2000$at2000.m)
tab_br   <- data.frame(age = brems$idade[1:116], qx = brems$brems2021sb.m[1:116])
tab_pub  <- data.frame(age = pub2010_full$idade,  qx = pub2010_full$qx) |>
  filter(!is.na(qx))

res_at  <- testAdherence(fund_adh, table = tab_at,  alpha = 0.05)
res_br  <- testAdherence(fund_adh, table = tab_br,  alpha = 0.05)
res_pub <- testAdherence(fund_adh, table = tab_pub, alpha = 0.05)
Ocultar código 
res_at$summary |>
  select(test, description, family, statistic, p_value, reject_h0) |>
  gt() |>
  tab_header(
    title    = "Testes de aderência --- AT-2000 masculina",
    subtitle = md(sprintf("A/E = %.4f &nbsp;|&nbsp; Obs: %d &nbsp;|&nbsp; Esp: %.1f",
                          res_at$ae_ratio, res_at$observed_deaths, res_at$expected_deaths))
  ) |>
  cols_label(test="Teste", description="Descrição", family="Família",
             statistic="Estatística", p_value=md("$p$-valor"), reject_h0=md("Rejeita $H_0$?")) |>
  fmt_number(columns = c(statistic, p_value), decimals = 4) |>
  tab_style(style=cell_fill("#FDECEA"), locations=cells_body(rows=reject_h0=="Yes")) |>
  tab_style(style=cell_fill("#E8F5E9"), locations=cells_body(rows=reject_h0=="No"))  |>
  tab_style(style=cell_text(weight="bold"), locations=cells_body(columns=reject_h0)) |>
  tab_options(table.font.size = 13)
Testes de aderência --- AT-2000 masculina
A/E = 0.8018  |  Obs: 3658  |  Esp: 4562.4
Teste Descrição Família Estatística \(p\)-valor Rejeita \(H_0\)?
KS Kolmogorov-Smirnov Non-parametric 0.0508 NA Yes
ChiSquare Chi-Square Goodness-of-Fit Non-parametric 469.3741 0.0000 Yes
WaldI_P Wald Type I (Poisson) Poisson −13.3619 0.0000 Yes
WaldII_P Wald Type II — age slope (Poisson) Poisson NA 0.0000 Yes
LRTI_P LRT Type I (Poisson) Poisson 192.4492 0.0000 Yes
LRTII_P LRT Type II — age slope (Poisson) Poisson 199.4474 0.0000 Yes
BayesCI Bayesian Credibility Interval Bayesian NA NA Yes
WaldI_NB Wald Type I (Negative Binomial) Negative Binomial −8.1180 0.0000 Yes
WaldII_NB Wald Type II — age slope (Neg. Binomial) Negative Binomial NA 0.0000 Yes
LRTI_NB LRT Type I (Negative Binomial) Negative Binomial 44.2823 0.0000 Yes
LRTII_NB LRT Type II — age slope (Neg. Binomial) Negative Binomial 52.6410 0.0000 Yes
Ocultar código 
res_br$summary |>
  select(test, description, family, statistic, p_value, reject_h0) |>
  gt() |>
  tab_header(
    title    = "Testes de aderência --- BR-EMSsb-m v.2021",
    subtitle = md(sprintf("A/E = %.4f &nbsp;|&nbsp; Obs: %d &nbsp;|&nbsp; Esp: %.1f",
                          res_br$ae_ratio, res_br$observed_deaths, res_br$expected_deaths))
  ) |>
  cols_label(test="Teste", description="Descrição", family="Família",
             statistic="Estatística", p_value=md("$p$-valor"), reject_h0=md("Rejeita $H_0$?")) |>
  fmt_number(columns = c(statistic, p_value), decimals = 4) |>
  tab_style(style=cell_fill("#FDECEA"), locations=cells_body(rows=reject_h0=="Yes")) |>
  tab_style(style=cell_fill("#E8F5E9"), locations=cells_body(rows=reject_h0=="No"))  |>
  tab_style(style=cell_text(weight="bold"), locations=cells_body(columns=reject_h0)) |>
  tab_options(table.font.size = 13)
Testes de aderência --- BR-EMSsb-m v.2021
A/E = 0.8880  |  Obs: 3658  |  Esp: 4119.1
Teste Descrição Família Estatística \(p\)-valor Rejeita \(H_0\)?
KS Kolmogorov-Smirnov Non-parametric 0.0726 NA Yes
ChiSquare Chi-Square Goodness-of-Fit Non-parametric 408.5956 0.0000 Yes
WaldI_P Wald Type I (Poisson) Poisson −7.1810 0.0000 Yes
WaldII_P Wald Type II — age slope (Poisson) Poisson NA 0.0000 Yes
LRTI_P LRT Type I (Poisson) Poisson 53.6692 0.0000 Yes
LRTII_P LRT Type II — age slope (Poisson) Poisson 60.4761 0.0000 Yes
BayesCI Bayesian Credibility Interval Bayesian NA NA Yes
WaldI_NB Wald Type I (Negative Binomial) Negative Binomial −5.2910 0.0000 Yes
WaldII_NB Wald Type II — age slope (Neg. Binomial) Negative Binomial NA 0.0000 Yes
LRTI_NB LRT Type I (Negative Binomial) Negative Binomial 23.2139 0.0000 Yes
LRTII_NB LRT Type II — age slope (Neg. Binomial) Negative Binomial 28.1877 0.0000 Yes
Ocultar código 
res_pub$summary |>
  select(test, description, family, statistic, p_value, reject_h0) |>
  gt() |>
  tab_header(
    title    = "Testes de aderência --- Pub-2010 Amount-Weighted Retiree masculina",
    subtitle = md(sprintf("A/E = %.4f &nbsp;|&nbsp; Obs: %d &nbsp;|&nbsp; Esp: %.1f",
                          res_pub$ae_ratio, res_pub$observed_deaths, res_pub$expected_deaths))
  ) |>
  cols_label(test="Teste", description="Descrição", family="Família",
             statistic="Estatística", p_value=md("$p$-valor"), reject_h0=md("Rejeita $H_0$?")) |>
  fmt_number(columns = c(statistic, p_value), decimals = 4) |>
  tab_style(style=cell_fill("#FDECEA"), locations=cells_body(rows=reject_h0=="Yes")) |>
  tab_style(style=cell_fill("#E8F5E9"), locations=cells_body(rows=reject_h0=="No"))  |>
  tab_style(style=cell_text(weight="bold"), locations=cells_body(columns=reject_h0)) |>
  tab_options(table.font.size = 13)
Testes de aderência --- Pub-2010 Amount-Weighted Retiree masculina
A/E = 0.8207  |  Obs: 3499  |  Esp: 4263.3
Teste Descrição Família Estatística \(p\)-valor Rejeita \(H_0\)?
KS Kolmogorov-Smirnov Non-parametric 0.1110 NA Yes
ChiSquare Chi-Square Goodness-of-Fit Non-parametric 421.2508 0.0000 Yes
WaldI_P Wald Type I (Poisson) Poisson −11.6866 0.0000 Yes
WaldII_P Wald Type II — age slope (Poisson) Poisson NA 0.0000 Yes
LRTI_P LRT Type I (Poisson) Poisson 146.0334 0.0000 Yes
LRTII_P LRT Type II — age slope (Poisson) Poisson 285.7796 0.0000 Yes
BayesCI Bayesian Credibility Interval Bayesian NA NA Yes
WaldI_NB Wald Type I (Negative Binomial) Negative Binomial −4.5149 0.0000 Yes
WaldII_NB Wald Type II — age slope (Neg. Binomial) Negative Binomial NA 0.0000 Yes
LRTI_NB LRT Type I (Negative Binomial) Negative Binomial 17.0362 0.0000 Yes
LRTII_NB LRT Type II — age slope (Neg. Binomial) Negative Binomial 57.5949 0.0000 Yes

Os campos NA na coluna de estatística correspondem a testes cujo resultado é expresso por um \(p\)-valor direto sem estatística de teste escalar (Wald Tipo II, que usa estatística matricial) ou que não produzem \(p\)-valor analítico (KS assintótico, BayesCI baseado em intervalo). Isso não indica falha computacional; é uma característica do design dos testes Melo; Graziadei; Targino (2026).

As três tábuas rejeitam \(H_0\) nos 11 testes. A rejeição universal é esperada quando a tábua não foi calibrada especificamente para a população analisada (Melo; Graziadei; Targino, 2026): com 3658 óbitos, qualquer desvio sistemático de A/E gera rejeição praticamente certa para qualquer \(\alpha\) usual. O que importa é a magnitude das estatísticas. A BR-EMS produz as menores estatísticas de teste (Qui-quadrado: \(469 \to 421 \to 409\); LRT Poisson: \(192 \to 146 \to 54\) — AT-2000, Pub-2010 e BR-EMS respectivamente), confirmando que a BR-EMS é a tábua mais aderente. A Pub-2010 apresenta desajuste intermediário entre a AT-2000 e a BR-EMS: sendo uma tábua de aposentados de fundos de pensão americanos, captura melhor o perfil etário desta carteira do que a AT-2000 de seguros de vida individuais, mas ainda assim reflete uma população estrangeira com mortalidade estruturalmente diferente. Os testes NB também rejeitam para as três tábuas, descartando superdispersão como explicação alternativa.

A BR-EMSsb-m v.2021 é adotada como tábua de referência para o modelo de Bühlmann-Straub por ser a tábua com menor grau de desajuste, construída com experiência do mercado brasileiro e atualmente exigida pela regulação PREVIC.

Questão (iv): Definição dos Perfis de Risco

Grupos adotados

Ocultar código 
labels_fx_m4  <- c("[18, 40)", "[40, 60)", "[60, 80)", "80+")

tibble(
  Grupo = 1:4,
  Faixa = labels_fx_m4,
  Descrição = c(
    "Jovens trabalhadores; mortalidade baixa; CV = 22,4%",
    "Trabalhadores ativos; mortalidade crescente; CV = 3,8%",
    "Aposentados; 40,8% dos óbitos; CV = 6,0%",
    "Idosos avançados; mortalidade alta e instável; CV = 13,5%"
  )
) |>
  left_join(sum_faixas |> select(faixa, tot_exp, tot_obt, tx_media),
            by = c("Faixa" = "faixa")) |>
  gt() |>
  tab_header(
    title    = md("Perfis de risco para o modelo de Bühlmann-Straub ($m = 4$)"),
    subtitle = "A faixa [0, 18) é excluída do modelo por ausência de sinistros em 2012--2013"
  ) |>
  cols_label(
    Grupo="Grupo", Faixa="Faixa", Descrição="Caracterização",
    tot_exp="Expostos totais", tot_obt="Óbitos totais",
    tx_media=md("$\\hat{m}$ média")
  ) |>
  fmt_number(columns = c(tot_exp, tot_obt), decimals = 0) |>
  fmt_number(columns = tx_media, decimals = 6) |>
  tab_options(table.font.size = 13)
Perfis de risco para o modelo de Bühlmann-Straub (\(m = 4\))
A faixa [0, 18) é excluída do modelo por ausência de sinistros em 2012--2013
Grupo Faixa Caracterização Expostos totais Óbitos totais \(\hat{m}\) média
1 [18, 40) Jovens trabalhadores; mortalidade baixa; CV = 22,4% 180,635 89 0.000493
2 [40, 60) Trabalhadores ativos; mortalidade crescente; CV = 3,8% 126,107 288 0.002284
3 [60, 80) Aposentados; 40,8% dos óbitos; CV = 6,0% 84,038 1,491 0.017742
4 80+ Idosos avançados; mortalidade alta e instável; CV = 13,5% 26,155 1,790 0.068438

Questão (v): Modelo de Bühlmann-Straub — Frequência de Sinistros

Formulação do modelo

Seguindo Bühlmann; Gisler (2005) e Klugman; Panjer; Willmot (2012), define-se para cada grupo de risco \(i = 1, \ldots, m\) e período \(t = 1, 2\):

\[ F_{it} = \frac{N_{it}}{v_{it}}, \qquad E[N_{it} \mid \theta_i] = \text{Var}[N_{it} \mid \theta_i] = v_{it}\,\lambda(\theta_i), \quad \text{com } \lambda(\theta_i) = \lambda_0\,\theta_i, \]

onde \(N_{it}\) é o número de óbitos do grupo \(i\) no ano \(t\) e \(v_{it} = E_{it}\) é a exposição (pessoas-ano). A hipótese \(E[N_{it} \mid \theta_i] = \text{Var}[N_{it} \mid \theta_i]\) decorre diretamente da distribuição Poisson condicional, consequência natural da dinâmica de contagem do diagrama de Lexis (Rau et al., 2018).

Este modelo é um caso particular do modelo de Bühlmann-Straub com \(X_{it} = F_{it}\) e \(\sigma^2(\theta_i) = \lambda(\theta_i)\). Os parâmetros estruturais reduzem-se a apenas dois: o nível global \(\lambda_0\) e a variância entre grupos \(\tau^2\), com \(\sigma^2 = \lambda_0\) (Bühlmann; Gisler, 2005, cap. 4).

A faixa \([0, 18)\) é excluída da análise de credibilidade pois não apresenta sinistros em 2012–2013, tornando a frequência observada identicamente nula e inviabilizando a estimação de uma taxa própria. O modelo opera com \(m = 4\) grupos: \([18, 40)\), \([40, 60)\), \([60, 80)\) e \(80+\).

O modelo é ajustado com os dados de 2012 e 2013 (\(T = 2\)) e a previsão para 2014 é confrontada com os óbitos efetivamente observados nesse ano, servindo como validação fora da amostra.

Estimação iterativa de \(\lambda_0\) e \(\tau^2\)

Como os estimadores de \(\lambda_0\) e \(\tau^2\) são mutuamente dependentes, a estimação é feita iterativamente. Definindo as quantidades auxiliares (Bühlmann; Gisler, 2005):

\[ c_1 = \frac{m-1}{m}\left\{\sum_i \frac{v_i}{v}\left(1 - \frac{v_i}{v}\right)\right\}^{-1}, \qquad c_2 = \frac{m}{m-1}\sum_i \frac{v_i}{v}\left(\bar{F}_i^v - \bar{\bar{F}}^v\right)^2, \]

onde \(v_i = \sum_{t=1}^{2} v_{it}\) e \(v = \sum_i v_i\) (total de pessoas-ano em 2012–2013 para os \(m = 4\) grupos); a frequência volume ponderada do grupo \(i\) é \(\bar{F}_i^v = \sum_t (v_{it}/v_i)\,F_{it}\); e \(\bar{\bar{F}}^v = \sum_i (v_i/v)\,\bar{F}_i^v\) é a média global.

O algoritmo inicia com \(\hat{\lambda}_0^{(0)} = \bar{\bar{F}}^v\) e itera:

\[ \begin{aligned} \hat{\kappa}^{(k)} &= \frac{\hat{\lambda}_0^{(k)}}{\hat{\tau}^{2(k)}}, \qquad \hat{\omega}_i^{(k)} = \frac{v_i}{v_i + \hat{\kappa}^{(k)}},\\[6pt] \hat{\lambda}_0^{(k+1)} &= \sum_i \frac{\hat{\omega}_i^{(k)}}{\hat{\omega}^{(k)}}\,\bar{F}_i^v, \qquad \hat{\tau}^{2(k+1)} = \max\!\left\{0,\; c_1\!\left(c_2 - \frac{m\hat{\lambda}_0^{(k)}}{v}\right)\right\}. \end{aligned} \]

Ocultar código 
intervalos_num <- c(0, 18, 40, 60, 80, 116)

# Dados por grupo e ano — TODOS os 3 anos (necessário para previsão e validação)
bs_data <- fundo_long |>
  mutate(faixa = cut(Idade, breaks = intervalos_num, right = FALSE,
                     labels = labels_fx)) |>
  filter(!is.na(faixa)) |>
  group_by(faixa, Ano) |>
  summarise(vit = sum(Expostos), Nit = sum(Obitos), .groups = "drop") |>
  mutate(Fit = if_else(vit > 0, Nit / vit, 0))

# ── AJUSTE: apenas 2012-2013 e m=4 grupos (exclui [0,18) — zero sinistros) ──
bs_ajuste <- bs_data |>
  filter(Ano %in% c("2012", "2013"), faixa != "[0, 18)")

vi_df <- bs_ajuste |>
  group_by(faixa) |>
  summarise(
    vi     = sum(vit),
    Ni     = sum(Nit),
    Fbar_v = sum(vit * Fit) / sum(vit),
    .groups = "drop"
  )
v_total <- sum(vi_df$vi)
m       <- nrow(vi_df)

# Fbarbar (média global volume ponderada, 2012-2013)
Fbarbar <- sum(vi_df$vi * vi_df$Fbar_v) / v_total

# c1, c2
c1_denom <- sum((vi_df$vi / v_total) * (1 - vi_df$vi / v_total))
c1 <- ((m - 1) / m) / c1_denom
c2 <- (m / (m - 1)) * sum((vi_df$vi / v_total) * (vi_df$Fbar_v - Fbarbar)^2)

cat(sprintf("Ajuste: 2012-2013 | v = %d | m = %d | F̄̄ = %.6f\n",
            v_total, m, Fbarbar))
Ajuste: 2012-2013 | v = 276007 | m = 4 | F̄̄ = 0.008416
Ocultar código 
cat(sprintf("c1 = %.6f | c2 = %.8f\n", c1, c2))
c1 = 1.111937 | c2 = 0.00033687
Ocultar código 
# Algoritmo iterativo
lam0 <- Fbarbar
tau2 <- max(0, c1 * (c2 - m * lam0 / v_total))
iter_log <- tibble(k = 0, lam0 = lam0, tau2 = tau2)

for (k in seq_len(50)) {
  if (tau2 <= 0) break
  kappa    <- lam0 / tau2
  omega_i  <- vi_df$vi / (vi_df$vi + kappa)
  lam0_new <- sum(omega_i / sum(omega_i) * vi_df$Fbar_v)
  tau2_new <- max(0, c1 * (c2 - m * lam0_new / v_total))
  iter_log <- bind_rows(iter_log, tibble(k = k, lam0 = lam0_new, tau2 = tau2_new))
  if (abs(lam0_new - lam0) < 1e-10 && abs(tau2_new - tau2) < 1e-12) {
    lam0 <- lam0_new; tau2 <- tau2_new; break
  }
  lam0 <- lam0_new; tau2 <- tau2_new
}

cat(sprintf("\nConvergência em %d iterações\n", nrow(iter_log) - 1))

Convergência em 4 iterações
Ocultar código 
cat(sprintf("lambda0_hat = %.8f\n", lam0))
lambda0_hat = 0.02131241
Ocultar código 
cat(sprintf("tau2_hat    = %.8f\n", tau2))
tau2_hat    = 0.00037424
Ocultar código 
cat(sprintf("kappa_hat   = %.4f\n",  lam0 / tau2))
kappa_hat   = 56.9488
Ocultar código 
iter_log |>
  head(8) |>
  gt() |>
  tab_header(title = "Trajetória do algoritmo iterativo (ajuste 2012--2013)") |>
  cols_label(k    = "Iteração",
             lam0 = md("$\\hat{\\lambda}_0^{(k)}$"),
             tau2 = md("$\\hat{\\tau}^{2(k)}$")) |>
  fmt_number(columns = c(lam0, tau2), decimals = 8) |>
  tab_options(table.font.size = 13)
Trajetória do algoritmo iterativo (ajuste 2012--2013)
Iteração \(\hat{\lambda}_0^{(k)}\) \(\hat{\tau}^{2(k)}\)
0 0.00841645 0.00037445
1 0.02133048 0.00037424
2 0.02131239 0.00037424
3 0.02131241 0.00037424
4 0.02131241 0.00037424

Parâmetros, credibilidade e validação (2014)

Ocultar código 
kappa_hat <- lam0 / tau2

resultados <- vi_df |>
  mutate(
    omega_i  = vi / (vi + kappa_hat),
    lam_cred = omega_i * Fbar_v + (1 - omega_i) * lam0
  )

# Dados de 2014 para validação
val_2014 <- bs_data |>
  filter(Ano == "2014", faixa != "[0, 18)") |>
  select(faixa, v_2014 = vit, N_obs_2014 = Nit)

# Óbitos esperados pelas três tábuas em 2014
tabua_2014 <- fundo_long |>
  filter(Ano == "2014", Idade >= 18) |>
  left_join(
    data.frame(Idade = brems$idade[1:116],  qx_br  = brems$brems2021sb.m[1:116]),
    by = "Idade"
  ) |>
  left_join(
    data.frame(Idade = at2000$idade,        qx_at  = at2000$at2000.m),
    by = "Idade"
  ) |>
  left_join(
    data.frame(Idade = pub2010_full$idade,  qx_pub = pub2010_full$qx),
    by = "Idade"
  ) |>
  mutate(
    faixa  = cut(Idade, breaks = intervalos_num, right = FALSE, labels = labels_fx),
    D_br   = Expostos * qx_br,
    D_at   = Expostos * qx_at,
    D_pub  = Expostos * replace_na(qx_pub, 0)
  ) |>
  filter(!is.na(faixa)) |>
  group_by(faixa) |>
  summarise(
    D_br  = sum(D_br,  na.rm = TRUE),
    D_at  = sum(D_at,  na.rm = TRUE),
    D_pub = sum(D_pub, na.rm = TRUE),
    .groups = "drop"
  ) |>
  rename(D_tab = D_br)   # manter compatibilidade: D_tab = BR-EMS

# Tabela unificada: credibilidade + previsão + validação + três tábuas
previsao <- resultados |>
  left_join(val_2014,   by = "faixa") |>
  left_join(tabua_2014, by = "faixa") |>
  mutate(
    N_prev    = v_2014 * lam_cred,
    erro_bs   = if_else(N_obs_2014 > 0, (N_prev - N_obs_2014) / N_obs_2014 * 100, NA_real_),
    erro_br   = if_else(N_obs_2014 > 0, (D_tab  - N_obs_2014) / N_obs_2014 * 100, NA_real_),
    erro_at   = if_else(N_obs_2014 > 0, (D_at   - N_obs_2014) / N_obs_2014 * 100, NA_real_),
    erro_pub  = if_else(N_obs_2014 > 0 & D_pub > 0,
                        (D_pub - N_obs_2014) / N_obs_2014 * 100, NA_real_)
  )

# Helper: cor verde na célula com menor erro absoluto entre múltiplas colunas
previsao |>
  select(faixa, vi, Fbar_v, omega_i, lam_cred,
         v_2014, N_prev, N_obs_2014,
         D_tab, D_at, D_pub,
         erro_bs, erro_br, erro_at, erro_pub) |>
  bind_rows(
    tibble(
      faixa      = "Total",
      vi         = sum(previsao$vi),
      Fbar_v     = NA_real_,
      omega_i    = NA_real_,
      lam_cred   = NA_real_,
      v_2014     = sum(previsao$v_2014),
      N_prev     = sum(previsao$N_prev),
      N_obs_2014 = sum(previsao$N_obs_2014),
      D_tab      = sum(previsao$D_tab),
      D_at       = sum(previsao$D_at),
      D_pub      = sum(previsao$D_pub),
      erro_bs    = (sum(previsao$N_prev)  - sum(previsao$N_obs_2014)) /
                   sum(previsao$N_obs_2014) * 100,
      erro_br    = (sum(previsao$D_tab)   - sum(previsao$N_obs_2014)) /
                   sum(previsao$N_obs_2014) * 100,
      erro_at    = (sum(previsao$D_at)    - sum(previsao$N_obs_2014)) /
                   sum(previsao$N_obs_2014) * 100,
      erro_pub   = (sum(previsao$D_pub)   - sum(previsao$N_obs_2014)) /
                   sum(previsao$N_obs_2014) * 100
    )
  ) |>
  gt() |>
  tab_header(
    title    = "Bühlmann-Straub: estimação e validação (2014)",
    subtitle = md(sprintf(
      "$\\hat{\\lambda}_0 = %.5f$ &nbsp;|&nbsp; $\\hat{\\tau}^2 = %.6f$ &nbsp;|&nbsp; $\\hat{\\kappa} = %.2f$",
      lam0, tau2, kappa_hat
    ))
  ) |>
  cols_label(
    faixa      = "Grupo",
    vi         = md("$v_i$"),
    Fbar_v     = md("$\\bar{F}_i^v$"),
    omega_i    = md("$\\hat{\\omega}_i$"),
    lam_cred   = md("$\\hat{\\lambda}_i^H$"),
    v_2014     = md("$v_{i,2014}$"),
    N_prev     = md("$\\hat{N}_{i,2014}$ (BS)"),
    N_obs_2014 = md("$N_i^{\\text{obs}}$"),
    D_tab      = md("$\\hat{D}^{\\text{BR-EMS}}$"),
    D_at       = md("$\\hat{D}^{\\text{AT-2000}}$"),
    D_pub      = md("$\\hat{D}^{\\text{Pub-2010}}$"),
    erro_bs    = "Erro BS (%)",
    erro_br    = "Erro BR-EMS (%)",
    erro_at    = "Erro AT-2000 (%)",
    erro_pub   = "Erro Pub-2010 (%)"
  ) |>
  tab_spanner(label = "Ajuste 2012--2013",
              columns = c(vi, Fbar_v, omega_i, lam_cred)) |>
  tab_spanner(label = "Validação 2014",
              columns = c(v_2014, N_prev, N_obs_2014)) |>
  tab_spanner(label = "Tábuas de referência",
              columns = c(D_tab, D_at, D_pub)) |>
  tab_spanner(label = "Erros de previsão (%)",
              columns = c(erro_bs, erro_br, erro_at, erro_pub)) |>
  fmt_number(columns = vi,                    decimals = 0) |>
  fmt_number(columns = c(Fbar_v, lam_cred),   decimals = 6) |>
  fmt_number(columns = omega_i,               decimals = 4) |>
  fmt_number(columns = v_2014,                decimals = 0) |>
  fmt_number(columns = c(N_prev, N_obs_2014, D_tab, D_at, D_pub), decimals = 1) |>
  fmt_number(columns = c(erro_bs, erro_br, erro_at, erro_pub),     decimals = 1) |>
  sub_missing(missing_text = "—") |>
  tab_style(
    style     = list(cell_fill("#EBF5FB"), cell_text(weight = "bold")),
    locations = cells_body(rows = faixa == "Total")
  ) |>
  # Verde: BS vs. todas as tábuas
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_bs,
      rows    = !is.na(erro_bs) &
                abs(erro_bs) < abs(replace_na(erro_br, Inf)) &
                abs(erro_bs) < abs(replace_na(erro_at, Inf)) &
                abs(erro_bs) < abs(replace_na(erro_pub, Inf))
    )
  ) |>
  # Verde: melhor tábua por faixa (lógica inline — evita problema de tamanho de vetor)
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_br,
      rows    = !is.na(erro_br) & faixa != "Total" &
                abs(erro_br) <= abs(replace_na(erro_at, Inf)) &
                abs(erro_br) <= abs(replace_na(erro_pub, Inf))
    )
  ) |>
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_at,
      rows    = !is.na(erro_at) & faixa != "Total" &
                abs(erro_at) < abs(replace_na(erro_br, Inf)) &
                abs(erro_at) <= abs(replace_na(erro_pub, Inf))
    )
  ) |>
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_pub,
      rows    = !is.na(erro_pub) & faixa != "Total" &
                abs(erro_pub) < abs(replace_na(erro_br, Inf)) &
                abs(erro_pub) < abs(replace_na(erro_at, Inf))
    )
  ) |>
  tab_options(table.font.size = 11)
Bühlmann-Straub: estimação e validação (2014)
\(\hat{\lambda}_0 = 0.02131\)  |  \(\hat{\tau}^2 = 0.000374\)  |  \(\hat{\kappa} = 56.95\)
Grupo
Ajuste 2012--2013
Validação 2014
Tábuas de referência
Erros de previsão (%)
\(v_i\) \(\bar{F}_i^v\) \(\hat{\omega}_i\) \(\hat{\lambda}_i^H\) \(v_{i,2014}\) \(\hat{N}_{i,2014}\) (BS) \(N_i^{\text{obs}}\) \(\hat{D}^{\text{BR-EMS}}\) \(\hat{D}^{\text{AT-2000}}\) \(\hat{D}^{\text{Pub-2010}}\) Erro BS (%) Erro BR-EMS (%) Erro AT-2000 (%) Erro Pub-2010 (%)
[18, 40) 120,440 0.000556 0.9995 0.000566 60,195 34.1 22.0 50.4 40.1 0.0 54.9 129.0 82.3 —
[40, 60) 83,446 0.002265 0.9993 0.002278 42,661 97.2 99.0 125.5 152.8 100.7 −1.8 26.8 54.3 1.7
[60, 80) 55,125 0.017415 0.9990 0.017419 28,913 503.6 531.0 500.2 567.3 481.5 −5.2 −5.8 6.8 −9.3
80+ 16,996 0.065133 0.9967 0.064987 9,159 595.2 683.0 783.0 851.5 938.8 −12.9 14.6 24.7 37.5
Total 276,007 — — — 140,928 1,230.1 1,335.0 1,459.1 1,611.6 1,521.0 −7.9 9.3 20.7 13.9

\(\hat{\kappa} = 56.9\) pequeno indica heterogeneidade entre grupos grande relativa ao ruído de processo (\(\hat{\tau}^2/\hat{\sigma}^2 = 1/\hat{\kappa} \approx 1.8\%\)): as taxas variam muito mais entre faixas do que dentro de cada faixa ao longo dos anos. Como todos os \(v_i \gg \hat{\kappa}\), os fatores \(\hat{\omega}_i \approx 1\) — a experiência própria de cada faixa domina. O BS prevê 1230 sinistros para 2014 (erro global de -7.9%). A célula verde nas colunas de erro marca o melhor método por faixa. Nas colunas das tábuas, a célula verde indica qual tábua foi mais aderente em cada faixa — permitindo identificar um frankenstein atuarial: a combinação por faixa da tábua com menor erro absoluto, prática comum no mercado para otimizar o nível coletivo \(\hat{\lambda}_0\) do BS. O erro de 54.9% em \([18, 40)\) é consequência de uma queda atípica nos óbitos de 2014 (22 vs. média histórica de 33.5), que nenhum método conseguiria prever com apenas dois anos de dados.

Questão (vi): Inferência Bayesiana Completa para o Modelo Bühlmann-Straub

Motivação e formulação

O modelo de Bühlmann-Straub da questão (v) realiza inferência parcialmente especificada: estima \(\hat{\lambda}_0\) e \(\hat{\tau}^2\) por momentos, sem atribuir distribuições probabilísticas a \(\lambda_0\) e \(\theta_i\). A questão (vi) faz a inferência completamente bayesiana, especificando distribuições para todas as quantidades aleatórias do modelo, seguindo os slides da Semana 12.

O modelo hierárquico Poisson-Gamma é:

\[ N_{it} \mid \lambda_{it} \sim \text{Poisson}(\lambda_{it}), \qquad \lambda_{it} = v_{it}\,\theta_i\,\lambda_0, \]

\[ \theta_i \sim \text{Gamma}(\theta_0,\, \theta_0), \qquad \lambda_0 \sim \text{Gamma}(a, b), \qquad \theta_0 \sim \text{Gamma}(a, b). \]

O parâmetro \(\theta_i\) é o fator multiplicativo de risco do grupo \(i\) em relação à taxa coletiva \(\lambda_0\): note que \(E[\theta_i] = 1\) e \(\text{Var}[\theta_i] = \theta_0^{-1}\). A escolha \(\theta_i \sim \text{Gamma}(\theta_0, \theta_0)\) garante conjugação: a condicional completa de cada \(\theta_i\) é também Gamma, o que permitiria um amostrador de Gibbs exato. Aqui usamos o Stan (NUTS/HMC) para maior flexibilidade.

O ganho em relação ao BS clássico é a propagação completa da incerteza: em vez de substituir \(\hat{\lambda}_0\) e \(\hat{\tau}^2 = \theta_0/\lambda_0\) por estimativas pontuais, o modelo integra sobre toda a distribuição posterior, produzindo intervalos credíveis para os fatores \(\omega_i\) e para o número esperado de sinistros.

O estimador de credibilidade resulta diretamente dos transformed parameters do Stan:

\[ \hat{\omega}_i = \frac{v_i}{v_i + \theta_0/\lambda_0}, \qquad \hat{\lambda}_i^H = \hat{\omega}_i\,\bar{F}_i^v + (1 - \hat{\omega}_i)\,\lambda_0. \]

Preparação dos dados

Ocultar código 
# Matrizes v (exposição) e n (óbitos) para 2012-2013, m=4 grupos
# Ordem dos grupos: [18-40), [40-60), [60-80), 80+
labels_m4 <- c("[18, 40)", "[40, 60)", "[60, 80)", "80+")

bs_12_13 <- bs_data |>
  filter(Ano %in% c("2012", "2013"), faixa != "[0, 18)") |>
  mutate(faixa = droplevels(faixa)) |>
  arrange(Ano, faixa)

# Matrizes Ti × m (linhas = anos, colunas = grupos)
v_mat <- bs_12_13 |>
  select(Ano, faixa, vit) |>
  pivot_wider(names_from = faixa, values_from = vit) |>
  select(-Ano) |>
  as.matrix()

n_mat <- bs_12_13 |>
  select(Ano, faixa, Nit) |>
  pivot_wider(names_from = faixa, values_from = Nit) |>
  select(-Ano) |>
  as.matrix()

fi_mat <- n_mat / v_mat  # frequência por ano e grupo

# Frequência volume-ponderada por grupo (Find no enunciado)
Find_vec <- colSums(n_mat) / colSums(v_mat)

m_stan  <- ncol(v_mat)
Ti_stan <- nrow(v_mat)

# Dados de 2014 para validação
v_2014_vec <- bs_data |>
  filter(Ano == "2014", faixa != "[0, 18)") |>
  mutate(faixa = droplevels(faixa)) |>
  arrange(faixa) |>
  pull(vit)

n_2014_vec <- bs_data |>
  filter(Ano == "2014", faixa != "[0, 18)") |>
  mutate(faixa = droplevels(faixa)) |>
  arrange(faixa) |>
  pull(Nit)

Ajuste com Stan

O código Stan segue exatamente o modelo da Semana 12: verossimilhança Poisson com \(\lambda_{it} = v_{it}\theta_i\lambda_0\), e prioris vagas \(\text{Gamma}(0{,}001; 0{,}001)\) para \(\lambda_0\) e \(\theta_0\). O fator de credibilidade e o estimador de credibilidade são calculados no bloco transformed parameters.

Ocultar código 
library(rstan)
options(mc.cores       = parallel::detectCores())
options(rstan.progress = FALSE)
options(rstan.verbose  = FALSE)
rstan_options(auto_write = TRUE)

stan_code_bs <- "
data {
  int<lower=1> m;                        // número de grupos
  int<lower=1> Ti;                       // anos de ajuste
  vector<lower=0>[m] Find;              // frequência volume-ponderada por grupo
  array[Ti, m] int<lower=0> n;          // óbitos por ano e grupo
  matrix<lower=0>[Ti, m] v;             // exposição por ano e grupo
  matrix<lower=0>[Ti, m] fi;            // taxas observadas por ano e grupo
}
parameters {
  vector<lower=0>[m] theta;  // fator de risco por grupo (E[theta_i]=1)
  real<lower=0> theta0;      // parâmetro de forma da priori de theta
  real<lower=0> lambda0;     // taxa coletiva global
}
transformed parameters {
  vector<lower=0>[m] w;       // fator de credibilidade
  vector<lower=0>[m] mucred;  // estimador de credibilidade
  matrix<lower=0>[Ti, m] lambda;
  for (i in 1:m) {
    for (t in 1:Ti) {
      lambda[t, i] = v[t, i] * theta[i] * lambda0;
    }
  }
  for (i in 1:m) {
    w[i]      = sum(v[, i]) / (sum(v[, i]) + (theta0 / lambda0));
    mucred[i] = w[i] * Find[i] + (1 - w[i]) * lambda0;
  }
}
model {
  // Prioris vagas
  for (i in 1:m)
    theta[i] ~ gamma(theta0, theta0);
  lambda0 ~ gamma(0.001, 0.001);
  theta0  ~ gamma(0.001, 0.001);
  // Verossimilhança
  for (i in 1:m)
    for (t in 1:Ti)
      n[t, i] ~ poisson(lambda[t, i]);
}
"

stan_dados <- list(
  m     = m_stan,
  Ti    = Ti_stan,
  Find  = Find_vec,
  n     = n_mat,
  v     = v_mat,
  fi    = fi_mat
)

rds_path <- "fit_stan_efpc.rds"

# Verificar compatibilidade: theta tem dimensão m
rds_valido <- FALSE
if (file.exists(rds_path)) {
  fit_salvo <- tryCatch(readRDS(rds_path), error = function(e) NULL)
  if (!is.null(fit_salvo)) {
    dims_theta <- tryCatch(
      dim(rstan::extract(fit_salvo, pars = "theta")$theta),
      error = function(e) NULL
    )
    rds_valido <- !is.null(dims_theta) && dims_theta[2] == stan_dados$m
  }
  if (!rds_valido)
    message("Stan: rds incompatível — reajustando com m=", stan_dados$m)
}

if (rds_valido) {
  fit_stan <- fit_salvo
  message("Stan: carregando fit salvo (m=", stan_dados$m, ")")
} else {
  fit_stan <- stan(
    model_code = stan_code_bs,
    data       = stan_dados,
    warmup     = 10000,
    iter       = 30000,
    thin       = 30,
    chains     = 4,
    refresh    = 0
  )
  saveRDS(fit_stan, file = rds_path)
  message("Stan: fit concluído e salvo")
}
Ocultar código 
print(fit_stan,
      pars   = c("lambda0", "theta0", "theta"),
      digits = 4,
      probs  = c(0.025, 0.5, 0.975))
Inference for Stan model: anon_model.
4 chains, each with iter=30000; warmup=10000; thin=30; 
post-warmup draws per chain=667, total post-warmup draws=2668.

           mean se_mean      sd   2.5%    50%   97.5% n_eff   Rhat
lambda0  1.0890  0.3790 19.6483 0.0063 0.0267  0.9414  2687 0.9998
theta0   0.4332  0.0053  0.2761 0.0857 0.3715  1.1340  2672 0.9990
theta[1] 0.0273  0.0005  0.0286 0.0006 0.0206  0.0930  2725 0.9996
theta[2] 0.1102  0.0021  0.1112 0.0024 0.0843  0.3630  2767 0.9996
theta[3] 0.8469  0.0164  0.8656 0.0188 0.6487  2.7517  2777 0.9994
theta[4] 3.1633  0.0611  3.2135 0.0689 2.4243 10.3534  2762 0.9995

Samples were drawn using NUTS(diag_e) at Fri Jun 26 14:40:16 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).

Diagnóstico de convergência

Ocultar código 
s_lab   <- rstan::extract(fit_stan)
M_draws <- length(s_lab$lambda0)

theta_s  <- s_lab$theta
w_s      <- s_lab$w
mucred_s <- s_lab$mucred

colnames(theta_s)  <- labels_m4
colnames(w_s)      <- labels_m4
colnames(mucred_s) <- labels_m4

posterior_theta <- as.array(fit_stan,
  pars = c("theta[1]", "theta[2]", "theta[3]", "theta[4]",
           "lambda0", "theta0"))

mcmc_trace(posterior_theta) +
  theme_bw(base_size = 11) +
  theme(strip.background = element_rect(fill = "white", colour = "grey80"))

Trajetórias das cadeias MCMC para \(\theta_i\), \(\lambda_0\) e \(\theta_0\)

O \(\hat{R}\) máximo entre todos os parâmetros é 1 e o ESS mínimo é 2672, indicando convergência satisfatória das 4 cadeias. Os trace plots mostram picos esporádicos de alta amplitude — especialmente em \(\lambda_0\) e \(\theta[4]\) — que são consequência da fraca identificabilidade do produto \(\theta_i \lambda_0\): o modelo não consegue distinguir entre \(\lambda_0\) grande com \(\theta_i\) pequeno e vice-versa. Esse é um comportamento esperado no modelo Poisson-Gamma com prioris vagas e poucos grupos (\(m = 4\), \(T = 2\)). O ESS elevado e \(\hat{R} \approx 1\) confirmam que as cadeias exploram esse espaço de forma adequada e que as marginais posteriores são estáveis.

A mediana posterior de \(\lambda_0 = 0.02669\) (IC 95%: \([0.00626;\; 0.94137]\)) é a taxa coletiva de mortalidade — note que o IC é amplo por causa da não-identificabilidade. A interpretação relevante para credibilidade é a de \(\mu_i^{\text{cred}} = \omega_i \bar{F}_i^v + (1-\omega_i)\lambda_0\), que é identificável porque o produto \(\theta_i \lambda_0\) é o que determina os sinistros previstos. A mediana de \(\theta_0 = 0.37\) controla a dispersão dos fatores: quanto maior \(\theta_0\), menor \(\text{Var}[\theta_i] = 1/\theta_0\) e mais os grupos convergem para o nível coletivo.

Fatores de credibilidade e estimadores posteriores

Ocultar código 
w_df <- as.data.frame(w_s) |>
  pivot_longer(everything(), names_to = "Grupo", values_to = "omega") |>
  mutate(Grupo = factor(Grupo, levels = labels_m4))

ggplot(w_df, aes(x = omega, fill = Grupo, colour = Grupo)) +
  geom_density(alpha = 0.35, linewidth = 0.8) +
  geom_vline(xintercept = 1, linetype = "dashed",
             colour = "grey30", linewidth = 0.7) +
  scale_fill_scico_d(palette = "batlow", name = "Faixa") +
  scale_colour_scico_d(palette = "batlow", name = "Faixa",
                       guide = guide_legend(
                         override.aes = list(
                           fill   = scico::scico(4, palette = "batlow"),
                           shape  = 22,
                           size   = 4,
                           colour = "grey40",
                           alpha  = 0.8
                         ))) +
  labs(x = expression(hat(omega)[i]), y = "Densidade") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

Distribuição posterior de \(\omega_i\) por grupo. Linha tracejada: \(\omega_i = 1\) (experiência própria domina integralmente).
Ocultar código 
mucred_summary <- tibble(
  faixa   = labels_m4,
  mediana = apply(mucred_s, 2, median),
  lo95    = apply(mucred_s, 2, quantile, 0.025),
  hi95    = apply(mucred_s, 2, quantile, 0.975),
  w_med   = apply(w_s,      2, median)
) |>
  left_join(resultados |> select(faixa, lam_cred), by = "faixa")

mucred_summary |>
  gt() |>
  tab_header(
    title    = "Estimadores de credibilidade: BS vs. Bayesiano completo",
    subtitle = md("$\\mu_i^{\\text{cred}} = \\omega_i\\bar{F}_i^v + (1-\\omega_i)\\lambda_0$")
  ) |>
  cols_label(
    faixa    = "Grupo",
    mediana  = md("$\\tilde{\\mu}_i^{\\text{cred}}$ (Bayes, mediana)"),
    lo95     = "IC 95% inf.",
    hi95     = "IC 95% sup.",
    w_med    = md("$\\tilde{\\omega}_i$ (mediana)"),
    lam_cred = md("$\\hat{\\lambda}_i^H$ (BS)")
  ) |>
  fmt_number(columns = c(mediana, lo95, hi95, lam_cred), decimals = 6) |>
  fmt_number(columns = w_med, decimals = 4) |>
  tab_options(table.font.size = 13)
Estimadores de credibilidade: BS vs. Bayesiano completo
\(\mu_i^{\text{cred}} = \omega_i\bar{F}_i^v + (1-\omega_i)\lambda_0\)
Grupo \(\tilde{\mu}_i^{\text{cred}}\) (Bayes, mediana) IC 95% inf. IC 95% sup. \(\tilde{\omega}_i\) (mediana) \(\hat{\lambda}_i^H\) (BS)
[18, 40) 0.000559 0.000557 0.000565 0.9999 0.000566
[40, 60) 0.002269 0.002266 0.002277 0.9998 0.002278
[60, 80) 0.017417 0.017404 0.017425 0.9997 0.017419
80+ 0.065101 0.064894 0.065148 0.9991 0.064987

Os fatores \(\tilde{\omega}_i\) são próximos de 1 para todos os grupos — as exposições \(v_i\) são muito maiores que \(\hat{\kappa} = \theta_0/\lambda_0\), de modo que a experiência própria de cada faixa domina. Os estimadores bayesianos coincidem com os do BS clássico nas medianas; o diferencial está nos intervalos credíveis que propagam a incerteza sobre \(\lambda_0\) e \(\theta_0\).

Verificação preditiva posterior (PPC)

Antes de fazer previsões fora da amostra, verificamos se o modelo é capaz de reproduzir os dados de ajuste (2012–2013). A PPC simula 2.000 réplicas de \(N_{it} \mid \theta_i, \lambda_0 \sim \text{Poisson}(v_{it}\,\tilde{\theta}_i\,\lambda_0)\) e compara com os valores observados.

Ocultar código 
set.seed(42)
n_rep   <- 2000
idx_rep <- sample(seq_len(M_draws), n_rep)

grupos_anos <- expand.grid(grupo = labels_m4, ano = c("2012", "2013"),
                           stringsAsFactors = FALSE) |>
  mutate(
    i     = match(grupo, labels_m4),
    t     = match(ano, c("2012", "2013")),
    v     = mapply(function(ii, tt) v_mat[tt, ii], i, t),
    N_obs = mapply(function(ii, tt) n_mat[tt, ii], i, t)
  )

pred_draws <- matrix(NA_integer_, nrow = n_rep, ncol = nrow(grupos_anos))
for (r in seq_len(n_rep)) {
  k <- idx_rep[r]
  mu <- grupos_anos$v * s_lab$theta[k, grupos_anos$i] * s_lab$lambda0[k]
  pred_draws[r, ] <- rpois(nrow(grupos_anos), mu)
}

ppc_summary <- grupos_anos |>
  mutate(
    pred_med  = apply(pred_draws, 2, median),
    pred_lo95 = apply(pred_draws, 2, quantile, 0.025),
    pred_hi95 = apply(pred_draws, 2, quantile, 0.975),
    grupo     = factor(grupo, levels = labels_m4)
  )

ggplot(ppc_summary, aes(x = ano, colour = grupo)) +
  geom_linerange(aes(ymin = pred_lo95, ymax = pred_hi95),
                 linewidth = 4, alpha = 0.25,
                 position = position_dodge(width = 0.4)) +
  geom_point(aes(y = pred_med), size = 2.5,
             position = position_dodge(width = 0.4)) +
  geom_point(aes(y = N_obs), shape = 4, size = 3, stroke = 1.5,
             colour = "tomato",
             position = position_dodge(width = 0.4)) +
  scale_colour_scico_d(palette = "batlow", name = "Faixa") +
  guides(colour = guide_legend(override.aes = list(shape = 16, size = 3))) +
  labs(x = "Ano", y = "Óbitos",
       caption = "Ponto: mediana preditiva  •  Barra: IC 95%  •  ×: observado") +
  facet_wrap(~grupo, scales = "free_y", nrow = 1) +
  theme_bw(base_size = 12) +
  theme(legend.position = "none",
        strip.background = element_rect(fill = "white", colour = "grey80"))

Verificação preditiva posterior: IC 95% preditivo (barra) e mediana (ponto) vs. observado (×) por grupo e ano de ajuste.

Em todos os grupos e nos dois anos de ajuste, o valor observado (×) cai dentro do IC 95% preditivo — o modelo reproduz adequadamente a estrutura de contagem dos dados de ajuste. Note que a escala Y é livre por painel: [18,40) opera em dezenas de óbitos enquanto 80+ opera em centenas, mas o ajuste relativo é semelhante em todos os grupos.

Fatores de risco individuais (\(\theta_i\))

Ocultar código 
theta_df <- as.data.frame(s_lab$theta) |>
  setNames(labels_m4) |>
  pivot_longer(everything(), names_to = "Grupo", values_to = "theta") |>
  mutate(Grupo = factor(Grupo, levels = labels_m4))

theta_stats <- theta_df |>
  group_by(Grupo) |>
  summarise(med = median(theta), lo = quantile(theta, 0.025),
            hi  = quantile(theta, 0.975), .groups = "drop")

p_theta <- ggplot(theta_df, aes(x = theta, fill = Grupo, colour = Grupo)) +
  geom_density(alpha = 0.35, linewidth = 0.8) +
  geom_vline(xintercept = 1, linetype = "dashed",
             colour = "grey30", linewidth = 0.7) +
  geom_text(data = theta_stats,
            aes(x = med, y = Inf,
                label = sprintf("med = %.2f", med)),
            vjust = 1.5, hjust = 0.5, size = 3.0, inherit.aes = FALSE,
            colour = "grey20") +
  scale_fill_scico_d(palette   = "batlow", name = "Faixa") +
  scale_colour_scico_d(palette = "batlow", name = "Faixa") +
  labs(x = expression(theta[i]), y = "Densidade") +
  facet_wrap(~Grupo, scales = "free", nrow = 2) +
  theme_bw(base_size = 11) +
  theme(legend.position = "none",
        strip.background = element_rect(fill = "white", colour = "grey80"),
        axis.text.x      = element_text(size = 8),
        panel.spacing    = unit(0.8, "lines"))

p_theta

Distribuições posteriores de \(\theta_i\) por grupo (escala livre). Linha tracejada: \(\theta_i = 1\) (nível coletivo \(\lambda_0\)).

Os fatores \(\theta_i\) quantificam o risco relativo de cada grupo em relação ao nível coletivo \(\lambda_0\): \(\theta_i < 1\) indica mortalidade abaixo da média, \(\theta_i > 1\) acima. As medianas posteriores são \(\tilde{\theta}_{[18,40)} = 0.02\), \(\tilde{\theta}_{[40,60)} = 0.08\), \(\tilde{\theta}_{[60,80)} = 0.65\) e \(\tilde{\theta}_{80+} = 2.42\). A razão entre os extremos é \(\tilde{\theta}_{80+} / \tilde{\theta}_{[18,40)} \approx 117.5\), coerente com a razão direta de taxas brutas (\(\hat{m}_{80+}/\hat{m}_{[18,40)} \approx 139\)) ponderada pelos pesos de volume — a diferença entre as duas razões reflete o encolhimento bayesiano em direção a \(\lambda_0\). As distribuições posteriores são assimétricas à direita, com caudas pesadas especialmente em \([60,80)\) e \(80+\), consequência da não-identificabilidade parcial do produto \(\theta_i\lambda_0\) com prioris vagas.

Análise de sensibilidade às prioris

Para verificar a robustez das estimativas de credibilidade à especificação das prioris, comparamos dois cenários para o hiperparâmetro \(a = b\) da priori de \(\lambda_0 \sim \text{Gamma}(a, b)\). A priori \(\text{Gamma}(a, b)\) tem \(E[\lambda_0] = a/b\) e \(\text{Var}[\lambda_0] = a/b^2\); com \(a = b\) tem-se \(E[\lambda_0] = 1\) e \(\text{DP}[\lambda_0] = 1/\sqrt{a}\).

O cenário vago (\(a = b = 0{,}001\)) é a priori adotada no modelo principal: moda em 0, média em 1 e DP em \(\approx 31{,}6\) — essencialmente não-informativa, colocando massa em toda a semirreta positiva. O cenário moderado (\(a = b = 1\)) corresponde a \(\text{Gamma}(1, 1) = \text{Exponencial}(1)\), com média 1 e DP 1 — ainda vago para uma taxa de mortalidade da ordem de \(0{,}02\), mas com mais concentração próxima à origem. Uma priori genuinamente informativa sobre \(\lambda_0\) exigiria conhecimento externo independente dos dados (por exemplo, experiência histórica de outras carteiras), o que não está disponível; além disso, centrar a priori no \(\hat{\lambda}_0 = 0{,}021\) estimado pelo BS introduziria circularidade, usando os dados duas vezes. Por essa razão, limitamos a análise a prioris vagas.

Ocultar código 
cenarios <- list(
  "Vaga ($a = b = 0{,}001$)"    = 0.001,
  "Moderada ($a = b = 1$)"      = 1.000
)

sens_results <- list()

for (nome in names(cenarios)) {
  ab <- cenarios[[nome]]

  stan_code_sens <- sprintf("
data {
  int<lower=1> m;
  int<lower=1> Ti;
  vector<lower=0>[m] Find;
  array[Ti, m] int<lower=0> n;
  matrix<lower=0>[Ti, m] v;
  matrix<lower=0>[Ti, m] fi;
}
parameters {
  vector<lower=0>[m] theta;
  real<lower=0> theta0;
  real<lower=0> lambda0;
}
transformed parameters {
  vector<lower=0>[m] w;
  vector<lower=0>[m] mucred;
  matrix<lower=0>[Ti, m] lambda;
  for (i in 1:m)
    for (t in 1:Ti)
      lambda[t, i] = v[t, i] * theta[i] * lambda0;
  for (i in 1:m) {
    w[i]      = sum(v[, i]) / (sum(v[, i]) + (theta0 / lambda0));
    mucred[i] = w[i] * Find[i] + (1 - w[i]) * lambda0;
  }
}
model {
  for (i in 1:m) theta[i] ~ gamma(theta0, theta0);
  lambda0 ~ gamma(%.3f, %.3f);
  theta0  ~ gamma(0.001, 0.001);
  for (i in 1:m)
    for (t in 1:Ti)
      n[t, i] ~ poisson(lambda[t, i]);
}
", ab, ab)

  rds_sens <- sprintf("fit_stan_sens_%s.rds",
                     gsub("[^a-zA-Z0-9]", "_", nome))
  if (file.exists(rds_sens)) {
    fit_s <- readRDS(rds_sens)
    message("Sensibilidade: carregando ", rds_sens)
  } else {
    fit_s <- stan(model_code = stan_code_sens, data = stan_dados,
                  warmup = 5000, iter = 15000, thin = 15,
                  chains = 2, refresh = 0, verbose = FALSE)
    saveRDS(fit_s, file = rds_sens)
    message("Sensibilidade: salvo ", rds_sens)
  }
  ex <- rstan::extract(fit_s)
  sens_results[[nome]] <- tibble(
    cenario = nome,
    faixa   = labels_m4,
    mucred_med = apply(ex$mucred, 2, median),
    mucred_lo  = apply(ex$mucred, 2, quantile, 0.025),
    mucred_hi  = apply(ex$mucred, 2, quantile, 0.975)
  )
}

sens_df <- bind_rows(sens_results) |>
  mutate(cenario = factor(cenario, levels = names(cenarios)))
Ocultar código 
nomes_cenarios <- names(cenarios)
# Renomear colunas para nomes simples antes do gt, depois aplicar md() no cols_label
col1 <- nomes_cenarios[1]
col2 <- nomes_cenarios[2]

sens_df |>
  select(cenario, faixa, mucred_med) |>
  pivot_wider(names_from = cenario, values_from = mucred_med) |>
  left_join(resultados |> select(faixa, lam_cred), by = "faixa") |>
  rename(vaga = all_of(col1), moderada = all_of(col2)) |>
  gt() |>
  tab_header(
    title    = md("Análise de sensibilidade: $\\tilde{\\mu}_i^{\\text{cred}}$ sob duas prioris para $\\lambda_0$"),
    subtitle = md("$\\lambda_0 \\sim \\text{Gamma}(a, b)$ com $a = b$ variando")
  ) |>
  cols_label(
    faixa    = "Grupo",
    vaga     = md("Vaga ($a = b = 0{,}001$)"),
    moderada = md("Moderada ($a = b = 1$)"),
    lam_cred = md("BS ($\\hat{\\lambda}_i^H$)")
  ) |>
  fmt_number(columns = where(is.numeric), decimals = 6) |>
  tab_spanner(label   = md("Bayesiano — mediana posterior de $\\mu_i^{\\text{cred}}$"),
              columns = c(vaga, moderada)) |>
  tab_options(table.font.size = 13)
Análise de sensibilidade: \(\tilde{\mu}_i^{\text{cred}}\) sob duas prioris para \(\lambda_0\)
\(\lambda_0 \sim \text{Gamma}(a, b)\) com \(a = b\) variando
Grupo
Bayesiano — mediana posterior de \(\mu_i^{\text{cred}}\)
BS (\(\hat{\lambda}_i^H\))
Vaga (\(a = b = 0{,}001\)) Moderada (\(a = b = 1\))
[18, 40) 0.000559 0.000559 0.000566
[40, 60) 0.002269 0.002268 0.002278
[60, 80) 0.017417 0.017418 0.017419
80+ 0.065103 0.065137 0.064987

A variação máxima das estimativas de \(\tilde{\mu}_i^{\text{cred}}\) entre os três cenários é de 0.12% — confirmando que os prêmios de credibilidade são robustos à especificação da priori de \(\lambda_0\). Isso é esperado: como \(\tilde{\omega}_i \approx 1\) para todos os grupos, \(\tilde{\mu}_i^{\text{cred}} \approx \bar{F}_i^v\) independentemente de \(\lambda_0\), e a priori sobre \(\lambda_0\) tem influência negligenciável sobre a quantidade de interesse.

Previsão bayesiana para 2014 e comparação final

O número esperado de sinistros para 2014 é obtido com os \(\theta_i\) posteriores e a exposição observada:

\[ \hat{N}_{i,2014}^{\text{Bayes}} = v_{i,2014} \cdot \tilde{\theta}_i \cdot \lambda_0. \]

Ocultar código 
prev_bayes |>
  select(faixa, N_BS, N_bayes, N_lo95, N_hi95, N_br, N_at, N_pub, N_obs) |>
  bind_rows(
    tibble(
      faixa   = "Total",
      N_BS    = sum(prev_bayes$N_BS),
      N_bayes = sum(prev_bayes$N_bayes),
      N_lo95  = sum(prev_bayes$N_lo95),
      N_hi95  = sum(prev_bayes$N_hi95),
      N_br    = sum(prev_bayes$N_br),
      N_at    = sum(prev_bayes$N_at),
      N_pub   = sum(prev_bayes$N_pub),
      N_obs   = sum(prev_bayes$N_obs)
    )
  ) |>
  gt() |>
  tab_header(
    title    = "Previsão para 2014: credibilidade e tábuas de referência",
    subtitle = md("Ajuste com 2012--2013; $v_i^* = E_{i,2014}$")
  ) |>
  cols_label(
    faixa   = "Grupo",
    N_BS    = md("$\\hat{N}_i^*$ (BS)"),
    N_bayes = md("$\\hat{N}_i^{\\text{Bayes}}$ (med.)"),
    N_lo95  = "IC 95% inf.",
    N_hi95  = "IC 95% sup.",
    N_br    = "BR-EMS",
    N_at    = "AT-2000",
    N_pub   = "Pub-2010",
    N_obs   = md("$N_i^{\\text{obs}}$ (2014)")
  ) |>
  fmt_number(columns = c(N_BS, N_bayes, N_lo95, N_hi95, N_br, N_at, N_pub, N_obs),
             decimals = 1) |>
  sub_missing(missing_text = "—") |>
  tab_style(
    style     = list(cell_fill("#EBF5FB"), cell_text(weight = "bold")),
    locations = cells_body(rows = faixa == "Total")
  ) |>
  tab_spanner(label = "Bühlmann-Straub",           columns = N_BS) |>
  tab_spanner(label = "Bayesiano completo (Stan)",  columns = c(N_bayes, N_lo95, N_hi95)) |>
  tab_spanner(label = "Tábuas de referência",       columns = c(N_br, N_at, N_pub)) |>
  tab_spanner(label = "Observado",                  columns = N_obs) |>
  tab_options(table.font.size = 13)
Previsão para 2014: credibilidade e tábuas de referência
Ajuste com 2012–2013; \(v_i^* = E_{i,2014}\)
Grupo
Bühlmann-Straub
Bayesiano completo (Stan)
Tábuas de referência
Observado
\(\hat{N}_i^*\) (BS) \(\hat{N}_i^{\text{Bayes}}\) (med.) IC 95% inf. IC 95% sup. BR-EMS AT-2000 Pub-2010 \(N_i^{\text{obs}}\) (2014)
[18, 40) 34.1 33.6 26.3 42.3 50.4 40.1 0.0 22.0
[40, 60) 97.2 96.5 84.1 111.2 125.5 152.8 100.7 99.0
[60, 80) 503.6 503.2 472.5 536.4 500.2 567.3 481.5 531.0
80+ 595.2 596.0 561.3 632.6 783.0 851.5 938.8 683.0
Total 1,230.1 1,229.2 1,144.3 1,322.4 1,459.1 1,611.6 1,521.0 1,335.0

Distribuição preditiva do total de sinistros (2014)

A distribuição preditiva posterior é obtida por simulação de Monte Carlo: para cada draw \((\tilde{\theta}_i, \tilde{\lambda}_0)\) da posterior, simula-se \(N_{i,2014}^{\text{rep}} \sim \text{Poisson}(v_{i,2014}\,\tilde{\theta}_i\,\tilde{\lambda}_0)\) e soma-se \(N_{\text{total}}^{\text{rep}} = \sum_i N_{i,2014}^{\text{rep}}\). O resultado é uma distribuição completa sobre o número total de óbitos em 2014, que propaga simultaneamente a incerteza sobre a taxa coletiva \(\lambda_0\), os fatores individuais \(\theta_i\) e o processo Poisson de contagem — algo que o BS clássico não oferece, pois usa estimativas pontuais de \(\hat{\lambda}_0\) e \(\hat{\tau}^2\).

Ocultar código 
ggplot(tibble(total = total_draws_bayes), aes(x = total)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 60, fill = "steelblue", alpha = 0.7, colour = NA) +
  geom_density(colour = "steelblue4", linewidth = 1.0) +
  geom_vline(xintercept = p_med_b, colour = "steelblue4",
             linetype = "solid",  linewidth = 1.0) +
  geom_vline(xintercept = p_lo_b, colour = "steelblue4",
             linetype = "dashed", linewidth = 0.8) +
  geom_vline(xintercept = p_hi_b, colour = "steelblue4",
             linetype = "dashed", linewidth = 0.8) +
  geom_vline(xintercept = obs_2014, colour = "#E74C3C",
             linetype = "dotted", linewidth = 0.9) +
  annotate("text", x = p_med_b - 5, y = 0.020,
           label = sprintf("Mediana\n%.0f", p_med_b),
           hjust = 1, size = 3.2, colour = "steelblue4") +
  annotate("text", x = obs_2014 + 5, y = 0.018,
           label = sprintf("Observado\n%d", obs_2014),
           hjust = 0, size = 3.2, colour = "#E74C3C") +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.12))) +
  labs(x = "Total de sinistros previstos para 2014", y = "Densidade") +
  theme_bw(base_size = 12)

Distribuição preditiva posterior do total de sinistros em 2014. Linha sólida azul: mediana preditiva. Linhas tracejadas: IC 95%. Linha pontilhada vermelha: observado (1335 óbitos).

A mediana da distribuição preditiva é \(\approx 1230\) sinistros, com IC 95% de \([1180;\; 1282]\). O valor observado de 1335 situa-se fora do IC 95% preditivo, com erro de \(-7.9\%\) em relação à mediana. Isso não indica falha do modelo — indica que a aceleração de 2014 (crescimento de 8,6‰ para 9,5‰) não era previsível a partir de apenas dois anos de dados com trajetória estável. Um IC 95% que exclui o observado em \(\approx 1\) em cada 20 validações é esperado por construção; o que seria preocupante seria o observado estar sistematicamente acima do IC superior em múltiplos anos consecutivos.

O P95 da distribuição (\(\approx 1274\) sinistros) representa a reserva de solvência mínima para cobrir os sinistros de 2014 com probabilidade \(95\%\), sendo 3.6% acima da estimativa central. Em termos atuariais, esse carregamento corresponde ao prêmio de risco acima do prêmio puro — a incerteza epistêmica sobre \((\lambda_0, \theta_i)\) e a variância de processo Poisson contribuem em proporções distintas para essa margem, sendo a segunda dominante quando \(\tilde{\omega}_i \approx 1\).

Tabela comparativa final: BS, Bayesiano completo e Tábuas

Ocultar código 
comp_plot <- prev_bayes |>
  mutate(grupo_num = seq_along(faixa))

# Dados no formato longo para as linhas
linhas_df <- bind_rows(
  comp_plot |> transmute(grupo_num, faixa, valor = N_BS,    metodo = "BS"),
  comp_plot |> transmute(grupo_num, faixa, valor = N_bayes, metodo = "Bayes (mediana)"),
  comp_plot |> transmute(grupo_num, faixa, valor = N_br,    metodo = "BR-EMS v.2021"),
  comp_plot |> transmute(grupo_num, faixa, valor = N_at,    metodo = "AT-2000"),
  comp_plot |> transmute(grupo_num, faixa, valor = N_pub,   metodo = "Pub-2010"),
  comp_plot |> transmute(grupo_num, faixa, valor = N_obs,   metodo = "Observado")
) |>
  mutate(metodo = factor(metodo, levels = c(
    "Observado", "BS", "Bayes (mediana)", "BR-EMS v.2021", "AT-2000", "Pub-2010"
  )))

cores_comp <- c(
  "Observado"      = "black",
  "BS"             = "#2980B9",
  "Bayes (mediana)"= "#1ABC9C",
  "BR-EMS v.2021"  = "#27AE60",
  "AT-2000"        = "#E74C3C",
  "Pub-2010"       = "#E67E22"
)

# Filtrar NA (Pub-2010 sem cobertura em [18,40))
linhas_df_clean <- linhas_df |> filter(!is.na(valor), valor > 0)

ggplot(linhas_df_clean,
       aes(x = faixa, y = valor, colour = metodo, group = metodo)) +
  geom_ribbon(
    data = comp_plot |> filter(!is.na(N_lo95)),
    aes(x = faixa, ymin = N_lo95, ymax = N_hi95, group = 1),
    fill = "#1ABC9C", alpha = 0.15, inherit.aes = FALSE
  ) +
  geom_line(aes(linetype = metodo),   linewidth = 0.9) +
  geom_point(aes(shape = metodo),     size = 3.0) +
  scale_colour_manual(values = cores_comp, name = NULL) +
  scale_linetype_manual(
    values = c("Observado" = "dashed", "BS" = "solid", "Bayes (mediana)" = "solid",
               "BR-EMS v.2021" = "solid", "AT-2000" = "solid", "Pub-2010" = "solid"),
    name = NULL
  ) +
  scale_shape_manual(
    values = c("Observado" = 16, "BS" = 15, "Bayes (mediana)" = 17,
               "BR-EMS v.2021" = 18, "AT-2000" = 4, "Pub-2010" = 8),
    name = NULL
  ) +
  scale_y_continuous(
    trans  = "log10",
    labels = function(x) formatC(x, format = "fg", big.mark = ".", decimal.mark = ",")
  ) +
  labs(x = "Faixa etária",
       y = "Óbitos previstos (escala log)") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom",
        legend.key.width  = unit(1.2, "cm"),
        legend.box        = "horizontal") +
  guides(colour   = guide_legend(nrow = 2),
         linetype = guide_legend(nrow = 2),
         shape    = guide_legend(nrow = 2))

Previsão de óbitos para 2014 por grupo de risco: modelos de credibilidade vs. tábuas de referência. Linha tracejada preta: observado. Banda azul: IC 95% bayesiano.
Ocultar código 
prev_bayes |>
  bind_rows(
    tibble(
      faixa      = "Total",
      N_BS       = sum(prev_bayes$N_BS),
      N_bayes    = sum(prev_bayes$N_bayes),
      N_lo95     = sum(prev_bayes$N_lo95),
      N_hi95     = sum(prev_bayes$N_hi95),
      N_br       = sum(prev_bayes$N_br),
      N_at       = sum(prev_bayes$N_at),
      N_pub      = sum(prev_bayes$N_pub),
      N_obs      = sum(prev_bayes$N_obs),
      erro_bayes = (sum(prev_bayes$N_bayes) - sum(prev_bayes$N_obs)) /
                   sum(prev_bayes$N_obs) * 100,
      erro_BS    = (sum(prev_bayes$N_BS)    - sum(prev_bayes$N_obs)) /
                   sum(prev_bayes$N_obs) * 100,
      erro_br    = (sum(prev_bayes$N_br)    - sum(prev_bayes$N_obs)) /
                   sum(prev_bayes$N_obs) * 100,
      erro_at    = (sum(prev_bayes$N_at)    - sum(prev_bayes$N_obs)) /
                   sum(prev_bayes$N_obs) * 100,
      erro_pub   = (sum(prev_bayes$N_pub)   - sum(prev_bayes$N_obs)) /
                   sum(prev_bayes$N_obs) * 100
    )
  ) |>
  gt() |>
  tab_header(
    title    = "Comparação final: erros de previsão por método (validação 2014)",
    subtitle = md(sprintf("Observado total: **%d** óbitos | Ajuste: 2012--2013",
                          obs_2014))
  ) |>
  cols_label(
    faixa      = "Grupo",
    N_BS       = "BS",
    N_bayes    = "Bayes (med.)",
    N_lo95     = "IC 95% inf.",
    N_hi95     = "IC 95% sup.",
    N_br       = "BR-EMS",
    N_at       = "AT-2000",
    N_pub      = "Pub-2010",
    N_obs      = md("$N_i^{\\text{obs}}$"),
    erro_BS    = "BS (%)",
    erro_bayes = "Bayes (%)",
    erro_br    = "BR-EMS (%)",
    erro_at    = "AT-2000 (%)",
    erro_pub   = "Pub-2010 (%)"
  ) |>
  tab_spanner(label = "Credibilidade",
              columns = c(N_BS, N_bayes, N_lo95, N_hi95)) |>
  tab_spanner(label = "Tábuas de referência",
              columns = c(N_br, N_at, N_pub)) |>
  tab_spanner(label = "Erros (%)",
              columns = c(erro_BS, erro_bayes, erro_br, erro_at, erro_pub)) |>
  fmt_number(columns = c(N_BS, N_bayes, N_lo95, N_hi95, N_br, N_at, N_pub, N_obs),
             decimals = 1) |>
  fmt_number(columns = c(erro_BS, erro_bayes, erro_br, erro_at, erro_pub),
             decimals = 1) |>
  sub_missing(columns = c(erro_BS, erro_bayes, erro_br, erro_at, erro_pub),
              missing_text = "—") |>
  tab_style(
    style     = list(cell_fill("#EBF5FB"), cell_text(weight = "bold")),
    locations = cells_body(rows = faixa == "Total")
  ) |>
  # Verde nos métodos de credibilidade com menor erro absoluto vs. tábuas
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_BS,
      rows    = !is.na(erro_BS) & !is.na(erro_bayes) &
                abs(erro_BS) <= abs(replace_na(erro_br, Inf)) &
                abs(erro_BS) <= abs(replace_na(erro_at, Inf)) &
                abs(erro_BS) <= abs(replace_na(erro_pub, Inf)) &
                abs(erro_BS) <= abs(erro_bayes)
    )
  ) |>
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_bayes,
      rows    = !is.na(erro_bayes) &
                abs(erro_bayes) <= abs(replace_na(erro_br, Inf)) &
                abs(erro_bayes) <= abs(replace_na(erro_at, Inf)) &
                abs(erro_bayes) <= abs(replace_na(erro_pub, Inf)) &
                abs(erro_bayes) <= abs(replace_na(erro_BS, Inf))
    )
  ) |>
  # Verde na melhor tábua por faixa
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_br,
      rows    = !is.na(erro_br) & faixa != "Total" &
                abs(erro_br) <= abs(replace_na(erro_at, Inf)) &
                abs(erro_br) <= abs(replace_na(erro_pub, Inf))
    )
  ) |>
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_at,
      rows    = !is.na(erro_at) & faixa != "Total" &
                abs(erro_at) < abs(replace_na(erro_br, Inf)) &
                abs(erro_at) <= abs(replace_na(erro_pub, Inf))
    )
  ) |>
  tab_style(
    style     = cell_fill("#D5F5E3"),
    locations = cells_body(
      columns = erro_pub,
      rows    = !is.na(erro_pub) & faixa != "Total" &
                abs(erro_pub) < abs(replace_na(erro_br, Inf)) &
                abs(erro_pub) < abs(replace_na(erro_at, Inf))
    )
  ) |>
  tab_options(table.font.size = 11)
Comparação final: erros de previsão por método (validação 2014)
Observado total: 1335 óbitos | Ajuste: 2012–2013
Grupo \(N_i^{\text{obs}}\)
Credibilidade
Tábuas de referência
Erros (%)
BS Bayes (med.) IC 95% inf. IC 95% sup. BR-EMS AT-2000 Pub-2010 BS (%) Bayes (%) BR-EMS (%) AT-2000 (%) Pub-2010 (%)
[18, 40) 22.0 34.1 33.6 26.3 42.3 50.4 40.1 0.0 54.9 52.5 129.0 82.3 —
[40, 60) 99.0 97.2 96.5 84.1 111.2 125.5 152.8 100.7 −1.8 −2.6 26.8 54.3 1.7
[60, 80) 531.0 503.6 503.2 472.5 536.4 500.2 567.3 481.5 −5.2 −5.2 −5.8 6.8 −9.3
80+ 683.0 595.2 596.0 561.3 632.6 783.0 851.5 938.8 −12.9 −12.7 14.6 24.7 37.5
Total 1,335.0 1,230.1 1,229.2 1,144.3 1,322.4 1,459.1 1,611.6 1,521.0 −7.9 −7.9 9.3 20.7 13.9

A célula verde nas colunas de credibilidade (BS e Bayes) marca quando esse método supera todas as tábuas naquela faixa. Nas colunas de tábuas, a célula verde indica a tábua com menor erro absoluto por faixa.

Análise por faixa:

\([18, 40)\): todos os métodos erram muito (BS: \(-54{,}9\%\), tábuas: \(+82\%\) a \(+129\%\)) porque 2014 teve apenas 22 óbitos — queda de 34 para 22 em relação a 2012. Nenhum modelo com \(T = 2\) anos preveria esse evento; o erro é ruído puro de processo Poisson em baixa contagem.

\([40, 60)\): a Pub-2010 é a tábua com menor erro (\(+1{,}7\%\)), melhor até que o BS (\(-1{,}8\%\)) nessa faixa. Isso é coerente: a Pub-2010 foi calibrada para aposentados de fundos de pensão, população com perfil etário próximo ao de trabalhadores ativos de 40-60 anos em EFPCs brasileiras. A AT-2000 superestima em \(+54\%\) — reflexo do viés americano de seguros de vida individuais.

\([60, 80)\): BS e Bayes dominam (\(-5{,}2\%\)), com BR-EMS próxima (\(-5{,}8\%\)). Essa é a faixa com maior massa de óbitos (531 em 2014) e mais informação — o modelo de credibilidade captura bem a experiência histórica de 2012-2013. AT-2000 e Pub-2010 superestimam (\(+6{,}8\%\) e \(-9{,}3\%\) respectivamente), confirmando que tábuas americanas subestimam a velocidade de convergência da mortalidade brasileira na transição epidemiológica.

\(80+\): todas as tábuas superestimam substancialmente (BR-EMS: \(+14{,}6\%\), AT-2000: \(+24{,}7\%\), Pub-2010: \(+37{,}5\%\)). O BS e o Bayes subestimam em \(-12{,}9\%\) — o crescimento de 524 para 683 óbitos em dois anos é aceleração estrutural por envelhecimento que os dados de 2012-2013 não capturam completamente.

O frankenstein atuarial ótimo — combinando a tábua com menor erro por faixa — usaria: BS em \([18, 40)\), Pub-2010 em \([40, 60)\), BS em \([60, 80)\), BS em \(80+\). Erro total resultante: \(\approx |{-54{,}9}| \cdot (22/1335) + |{+1{,}7}| \cdot (99/1335) + |{-5{,}2}| \cdot (531/1335) + |{-12{,}9}| \cdot (683/1335) \approx 8{,}2\%\) — ligeiramente inferior ao BS puro (\(-7{,}9\%\) em valor absoluto), pois a Pub-2010 melhora \([40, 60)\) sem custo nas demais faixas.

Esse exercício ilustra o valor do estudo de experiência própria: com \(T \geq 5\) anos, a credibilidade do BS aumentaria substancialmente e a comparação com as tábuas seria mais informativa. Com \(T = 2\), o BS já supera as tábuas de referência em três das quatro faixas comparáveis — resultado notável dado o horizonte curto.

Informações da Sessão

Informações da Sessão 
sessionInfo()
R version 4.5.3 (2026-03-11 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United Kingdom.utf8 
[2] LC_CTYPE=English_United Kingdom.utf8   
[3] LC_MONETARY=English_United Kingdom.utf8
[4] LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.utf8    

time zone: America/Sao_Paulo
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] mortSOA_0.1.0            fastcpd_1.0.0            gtExtras_0.6.2          
 [4] ggiraph_0.9.6            scico_1.5.0              bayesplot_1.15.0        
 [7] rstan_2.36.0.9000        StanHeaders_2.36.0.9000  mortalityAdherence_0.1.0
[10] kableExtra_1.4.0         patchwork_1.3.2          gt_1.3.0                
[13] scales_1.4.0             lubridate_1.9.5          forcats_1.0.1           
[16] stringr_1.6.0            dplyr_1.2.1              purrr_1.2.2             
[19] readr_2.2.0              tidyr_1.3.2              tibble_3.3.1            
[22] ggplot2_4.0.3            tidyverse_2.0.0         

loaded via a namespace (and not attached):
 [1] httr2_1.2.2             gridExtra_2.3           inline_0.3.21          
 [4] rematch2_2.1.2          rlang_1.2.0             magrittr_2.0.5         
 [7] otel_0.2.0              matrixStats_1.5.0       compiler_4.5.3         
[10] loo_2.9.0               reshape2_1.4.5          systemfonts_1.3.2      
[13] vctrs_0.7.3             pkgconfig_2.0.3         shape_1.4.6.1          
[16] fastmap_1.2.0           backports_1.5.1         fontawesome_0.5.3      
[19] labeling_0.4.3          utf8_1.2.6              rmarkdown_2.31         
[22] markdown_2.0            tzdb_0.5.0              ragg_1.5.2             
[25] xfun_0.58               glmnet_5.0              litedown_0.9           
[28] jsonlite_2.0.0          parallel_4.5.3          R6_2.6.1               
[31] stringi_1.8.7           RColorBrewer_1.1-3      car_3.1-5              
[34] Rcpp_1.1.1-1.1          iterators_1.0.14        knitr_1.51             
[37] base64enc_0.1-6         splines_4.5.3           Matrix_1.7-5           
[40] timechange_0.4.0        tidyselect_1.2.1        rstudioapi_0.18.0      
[43] dichromat_2.0-0.1       abind_1.4-8             yaml_2.3.12            
[46] codetools_0.2-20        curl_7.1.0              pkgbuild_1.4.8         
[49] plyr_1.8.9              lattice_0.22-9          withr_3.0.2            
[52] S7_0.2.2                posterior_1.7.0         evaluate_1.0.5         
[55] survival_3.8-6          RcppParallel_5.1.11-2   xml2_1.5.2             
[58] pillar_1.11.1           tensorA_0.36.2.1        carData_3.0-6          
[61] checkmate_2.3.4         foreach_1.5.2           stats4_4.5.3           
[64] distributional_0.7.0    generics_0.1.4          paletteer_1.7.0        
[67] hms_1.1.4               commonmark_2.0.0        glue_1.8.1             
[70] gdtools_0.5.1           tools_4.5.3             fs_2.1.0               
[73] grid_4.5.3              QuickJSR_1.10.0         sfsmisc_1.1-24         
[76] Formula_1.2-5           cli_3.6.6               rappdirs_0.3.4         
[79] textshaping_1.0.5       fontBitstreamVera_0.1.1 viridisLite_0.4.3      
[82] svglite_2.2.2           V8_8.2.0                gtable_0.3.6           
[85] sass_0.4.10             digest_0.6.39           fontquiver_0.2.1       
[88] htmlwidgets_1.6.4       farver_2.1.2            htmltools_0.5.9        
[91] lifecycle_1.0.5         fontLiberation_0.1.0    MASS_7.3-65

Referências

ABRAPP. Consolidado Estatístico — Associação Brasileira das Entidades Fechadas de Previdência Complementar. https://www.abrapp.org.br, 2022.
BÜHLMANN, Hans; GISLER, Alois. A Course in Credibility Theory and Its Applications. Berlin: Springer, 2005.
FENAPREVI. Tábua BR-EMSsb-m v.2021 — Experiência do Mercado Segurador Brasileiro. São PauloFederação Nacional de Previdência Privada e Vida, 2021.
KILLICK, Rebecca; FEARNHEAD, Paul; ECKLEY, Idris A. Optimal Detection of Changepoints With a Linear Computational Cost. Journal of the American Statistical Association, v. 107, n. 500, p. 1590–1598, 2012.
KLUGMAN, Stuart A.; PANJER, Harry H.; WILLMOT, Gordon E. Loss Models: From Data to Decisions. 4. ed. Hoboken, NJ: John Wiley & Sons, 2012.
MELO, Eduardo Fraga L. de; GRAZIADEI, Helton; TARGINO, Rodrigo. mortalityAdherence: Mortality Table Adherence Tests. R package. https://github.com/eduardoflm/mortalityAdherence, 2026.
PITACCO, Ermanno et al. Modelling Longevity Dynamics for Pensions and Annuity Business. Oxford: Oxford University Press, 2009.
PREVIC. Portaria PREVIC no 835, de 3 de setembro de 2020. Superintendência Nacional de Previdência Complementar, 2020.
RAU, Roland et al. Visualizing Mortality Dynamics in the Lexis Diagram. Cham: Springer, 2018.
SOCIETY OF ACTUARIES. The 2000 Mortality Table. Schaumburg, ILSociety of Actuaries, 2000.
SOCIETY OF ACTUARIES. Pub-2010 Public Pension Mortality Tables. Society of Actuaries. https://mort.soa.org, 2014.