Estudo Simulado do modelo de Regressão Dirichlet

Abordagem Bayesiana - Parte 6

Pedro Frazão Dutra

20/07/2026

1 Definição do Modelo

O presente documento apresenta os resultados da aplicação do modelo de Regressão Dirichlet sob a parametrização alternativa (média e precisão) a dados simulados. O cenário considerado consiste em uma amostra de tamanho \(n = 30\) vetores composicionais de dimensão \(C = 3\). As covariáveis associadas foram geradas a partir de uma distribuição Uniforme centralizada na média, \(X_{1i} \sim \text{Uniforme}(0, 5)\), e de uma distribuição Bernoulli, \(X_{2i} \sim \text{Bernoulli}(p = 0{,}5)\).

O modelo probabilístico é especificado por: \[Y_i = (Y_{i1}, Y_{i2}, Y_{i3}) \sim \text{Dirichlet}(\alpha_{i1}, \alpha_{i2}, \alpha_{i3})\]

Na parametrização alternativa, cada parâmetro de forma \(\alpha_{ij}\) é decomposto pelo produto entre a proporção esperada \(\mu_{ij}\) e o parâmetro de precisão constante \(\phi > 0\): \[\alpha_{ij} = \mu_{ij} \cdot \phi, \quad \text{onde } \sum_{j=1}^{3} \mu_{ij} = 1\]

A média esperada \(\boldsymbol{\mu}_i\) é mapeada no simplex através da função Softmax: \[\mu_{ij} = \frac{\exp(\eta_{ij})}{\sum_{k=1}^{3} \exp(\eta_{ik})}, \quad \text{com } \eta_{ij} = \beta_{0j} + \beta_{1j} X_{1i} + \beta_{2j} X_{2i}\]

Para garantir a identificabilidade dos parâmetros na estimação, fixa-se uma categoria como referência (\(c_{\text{ref}}\)) com vetor de coeficientes nulo (\(\boldsymbol{\beta}_{c_{\text{ref}}} = \mathbf{0}\)). Dessa forma, o algoritmo MCMC estima os coeficientes relativos \(\boldsymbol{\beta}_j^* = \boldsymbol{\beta}_j - \boldsymbol{\beta}_{c_{\text{ref}}}\) para as demais componentes (\(j \neq c_{\text{ref}}\)).

Os valores verdadeiros utilizados na simulação foram a precisão \(\phi = 15\) (\(\log(\phi) \approx 2{,}708\)) e os coeficientes originais dispostos na matriz \(\boldsymbol{\beta}\): \[\boldsymbol{\beta} = \begin{bmatrix} \beta_{01} & \beta_{02} & \beta_{03} \\ \beta_{11} & \beta_{12} & \beta_{13} \\ \beta_{21} & \beta_{22} & \beta_{23} \end{bmatrix} = \begin{bmatrix} \phantom{-}0{,}5 & \phantom{-}0{,}0 & -0{,}5 \\ -0{,}8 & \phantom{-}0{,}6 & \phantom{-}0{,}3 \\ \phantom{-}1{,}2 & -0{,}9 & \phantom{-}0{,}4 \end{bmatrix}\]

As distribuições a priori independentes adotadas foram: \[\boldsymbol{\beta}_j^* \sim \mathcal{N}_3(\mathbf{0}, 100 \mathbf{I}_3) \quad \text{para } j \neq c_{\text{ref}}, \qquad \log(\phi) \sim \mathcal{N}(0, 100)\]

2 Pacotes utilizados

# Análise composicional
library(DirichletReg) 
library(compositions) 


# Produção de gráficos e tabelas
library(tidyverse)
library(ggtern)
library(colorspace)
library(purrr)
library(patchwork)
library(DT)          
library(knitr)       
library(kableExtra)  
library(htmltools)   
library(scales)

# Inferência Bayesiana
library(posterior)
library(MASS)
library(coda)

3 Geração dos dados

set.seed(2006)

# Tamanho da amostra
n <- 30

# Simulando e centralizando a covariável X1 (Uniforme)
x1_bruto <- runif(n, 0, 5)
x1 <- x1_bruto - mean(x1_bruto)

# Simulando a covariável X2 (Bernoulli com p = 0.5)
x2 <- rbinom(n, 1, 0.5)

# Matriz de design X (Intercepto, x1 e x2)
X_mat <- cbind(1, x1, x2)

# Matriz de Betas Reais 
beta_real <- matrix(c(
  0.5,   0.0, -0.5,  # Interceptos
  -0.8,   0.6,  0.3,  # Efeitos de X1 
  1.2,  -0.9,  0.4   # Efeitos de X2 
), nrow = 3, byrow = TRUE)

# Mapeamento do preditor linear para a média mu no simplex 
eta <- X_mat %*% beta_real
exp_eta <- exp(eta)
mu <- exp_eta / rowSums(exp_eta)  

# Definição da precisão real
phi_real <- 15
log_phi_real <- log(phi_real)

# Reconstrução dos alphas para geração no suporte Dirichlet
alpha <- mu * phi_real

# Gerando as composições no simplex
y <- matrix(0, nrow = n, ncol = 3)

for(i in 1:n) {
  z <- rgamma(3, shape = alpha[i, ], rate = 1)
  y[i, ] <- z / sum(z)
}

# Nomes para identificação dos parâmetros no MCMC / Tabela
nomes_param <- c("y1: Intercepto", "y1: Slope v1", "y1: Slope v2",
                 "y2: Intercepto", "y2: Slope v1", "y2: Slope v2",
                 "y3: Intercepto", "y3: Slope v1", "y3: Slope v2",
                 "log(phi)")

4 Análise descritiva

5 Seleção da Categoria de Referência

## -> Categoria selecionada como referência (Maior alfa médio): Y1

6 Rotina de Metropolis-Hastings

6.1 Funções de Log-Verossimilhança e Log-Posteriori

# ==============================================================================
# 1. LOG-VEROSSIMILHANÇA (Considerando y1 como categoria de referência)
# ==============================================================================
log_lik <- function(params, X, Y) {
  N <- nrow(Y)
  K <- ncol(X) # K = 3
  C <- ncol(Y) # C = 3
  
  # Extrai os betas para as categorias 2 e 3 (matriz K x (C-1))
  beta_estimados <- matrix(params[1:(K * (C - 1))], nrow = K, ncol = C - 1)
  
  # Inclui a coluna de zeros para a categoria 1 (Referência)
  beta_full <- cbind(0, beta_estimados)
  
  log_phi <- params[K * (C - 1) + 1]
  phi <- exp(log_phi)
  
  # Preditor linear e Softmax
  eta <- X %*% beta_full
  exp_eta <- exp(eta)
  mu <- exp_eta / rowSums(exp_eta)
  
  alpha <- mu * phi
  
  termo1 <- N * lgamma(phi)
  termo2 <- - sum(lgamma(alpha))
  termo3 <- sum((alpha - 1) * log(Y))
  
  return(termo1 + termo2 + termo3)
}

# ==============================================================================
# 2. LOG-PRIORIS
# ==============================================================================
log_prior <- function(params, K, C, prior_specs) {
  beta_estimados <- matrix(params[1:(K * (C - 1))], nrow = K, ncol = C - 1)
  log_phi <- params[K * (C - 1) + 1]
  
  log_prior_beta <- 0
  for (c in 1:(C - 1)) {
    diff_beta <- beta_estimados[, c] - prior_specs$mu_beta[, c]
    log_prior_beta <- log_prior_beta - 0.5 * as.numeric(t(diff_beta) %*% prior_specs$Sigma_beta_inv[[c]] %*% diff_beta)
  }
  
  diff_log_phi <- log_phi - prior_specs$mu_log_phi
  log_prior_log_phi <- - 0.5 * (diff_log_phi^2) / prior_specs$var_log_phi
  
  return(log_prior_beta + log_prior_log_phi)
}

# ==============================================================================
# 3. LOG-POSTERIORI
# ==============================================================================
log_post <- function(params, X, Y, K, C, prior_specs) {
  ll <- log_lik(params, X, Y)
  lp <- log_prior(params, K, C, prior_specs)
  return(ll + lp)
}

6.2 Algoritmo

K <- 3  # Dimensão do preditor linear (Intercepto, X1, X2)
d <- K * (C - 1) + 1  # d = 7 (6 betas estipulados + 1 log_phi)

X_matriz <- cbind(1, dados_modelo$x1, dados_modelo$x2)

# Prioris para as 2 categorias estimadas
var_beta <- 100
mu_beta <- matrix(0, nrow = K, ncol = C - 1)
Sigma_beta_inv <- lapply(1:(C - 1), function(c) diag(1 / var_beta, K))

prior_specs <- list(
  mu_beta = mu_beta,
  Sigma_beta_inv = Sigma_beta_inv,
  mu_log_phi = 0,
  var_log_phi = 100
)

# Extração dos estimadores MLE para proposta do MCMC (Tamanho d = 7)
params_mle <- as.numeric(unlist(coef(modelo_alt)))
Sigma_mle  <- vcov(modelo_alt)

tuning <- 1.5
Sigma_prop <- tuning * (2.4^2 / d) * Sigma_mle

# Nomeação automática dos parâmetros baseada na nova referência
colnames_params <- c(
  paste0(Nomes_Y[2], ": Intercepto"), paste0(Nomes_Y[2], ": Slope v1"), paste0(Nomes_Y[2], ": Slope v2"),
  paste0(Nomes_Y[3], ": Intercepto"), paste0(Nomes_Y[3], ": Slope v1"), paste0(Nomes_Y[3], ": Slope v2"),
  "log(phi)"
)

m <- 4         
nite <- 60000 
cadeias <- list()

set.seed(1006) 

for (i in 1:m) {
  tempo_inicio <- Sys.time()
  
  params_ini <- MASS::mvrnorm(1, mu = params_mle, Sigma = 2 * Sigma_mle)
  
  params_cadeia <- matrix(0, nrow = nite, ncol = d)
  colnames(params_cadeia) <- colnames_params
  params_cadeia[1, ] <- params_ini
  
  log_post_atual <- log_post(
    params = params_ini, 
    X = X_matriz, 
    Y = Y_matriz, 
    K = K, 
    C = C, 
    prior_specs = prior_specs
  )
  
  ruido_prop <- MASS::mvrnorm(nite - 1, mu = rep(0, d), Sigma = Sigma_prop)
  aceitos <- 0
  
  for (t in 1:(nite - 1)) {
    params_prop <- params_cadeia[t, ] + ruido_prop[t, ]
    
    log_post_prop <- log_post(
      params = params_prop, 
      X = X_matriz, 
      Y = Y_matriz, 
      K = K, 
      C = C, 
      prior_specs = prior_specs
    )
    
    if (is.na(log_post_prop) || is.infinite(log_post_prop)) {
      log_alfa <- -Inf 
    } else {
      log_alfa <- log_post_prop - log_post_atual
    }
    
    if (!is.na(log_alfa) && log(runif(1)) <= log_alfa) {
      params_cadeia[t + 1, ] <- params_prop
      log_post_atual <- log_post_prop
      aceitos <- aceitos + 1
    } else {
      params_cadeia[t + 1, ] <- params_cadeia[t, ]
    }
  }
  
  cadeias[[i]] <- params_cadeia
  tempo_fim <- Sys.time()
  tempo_execucao <- round(difftime(tempo_fim, tempo_inicio, units = "secs"), 2)
  
  cat("Cadeia", i, 
      "- Ref:", Nomes_Y[1],
      "| Aceitação:", round((aceitos / (nite - 1)) * 100, 2), "%",
      "| Tempo:", tempo_execucao, "s\n")
}

## Cadeia 1 - Ref: Y1 | Aceitação: 24.29 % | Tempo: 6.09 s
## Cadeia 2 - Ref: Y1 | Aceitação: 24.36 % | Tempo: 6.16 s
## Cadeia 3 - Ref: Y1 | Aceitação: 24.06 % | Tempo: 6 s
## Cadeia 4 - Ref: Y1 | Aceitação: 23.7 % | Tempo: 6.14 s

6.3 Diagnósticos de convergência

6.3.1 Cadeias piloto

Nesta etapa, simulamos 4 cadeias para cada um dos parâmetos, todas com inicialização superdispersa. Com isso, podemos checar se existe multimodalidade em cada uma das distribuições posteriori. Além disso, para avaliar a convergência de cada grupo de cadeias para a uma mesma distribuição comum, calculamos a estatística potencial de redução de escala, \(\hat{R}\).

6.3.2 Funções de Autocorrelação

Plotamos as funções de Autocorrelação das cadeias para verificar a eficiência com que o algoritmo explora o espaço paramétrico da distribuições posteriori. Além disso, é importante descobrir o menor lag k tal que as ACF’s apresentem todos os seus valores dentro do intervalo de confiança centrado em zero. Com isso, se for preciso, poderemos efetuar um espaçamento de tamanho \(k\) em todas as cadeia, de modo a restarem somente amostras estatísticamente independentes para realizarmos estimações. Como existem muitos parâmetros, plotamos apenas as ACF’s da primeira cadeia. Isto é adequado, uma vez que as quatro cadeias apresentaram um comportamento empíricamente semelhante.

6.3.3 R-hat & ESS

Por fim, calculamos a estatística de redução de escala potencial, \(\hat{R}\) e o tamanho efetivo das amostras para estimação de quantidades próximas do centro de massa da distribuição posteriori e para estimação de quantidades próximas das caudas. Para cada parâmetro, juntamos as cadeias aquecidas e realizamos ambos os cálculos.

Diagnósticos de Convergência
Parâmetro	R-Hat	ESS Bulk	ESS Tail
Y2: Intercepto	1.001	4213	7446
Y2: Slope v1	1.001	5490	8782
Y2: Slope v2	1.001	4726	7911
Y3: Intercepto	1.001	4890	7501
Y3: Slope v1	1.001	4747	7637
Y3: Slope v2	1.001	4696	7525
log(phi)	1.001	5565	8688

7 Inferência Bayesiana

7.1 Distribuições marginais

Assumido com segurança o bom desempenho do algoritmo, plotamos os histogramas das distribuições marginais e calculamos medidas resumos usuais, como média, mediana e intervalos de 95% de credibilidade.

##Estimativas pontuais e intervalares

Estimativas Pontuais e Intervalos de 95% de Credibilidade
Parâmetro	Valor Real Relativo	Média Post.	Mediana	Desvio Padrão	Quantil 2.5%	Quantil 97.5%
Y2: Intercepto	-0.500	-0.404	-0.405	0.216	-0.829	0.023
Y2: Slope v1	1.400	1.534	1.535	0.103	1.329	1.733
Y2: Slope v2	-2.100	-2.244	-2.243	0.285	-2.803	-1.686
Y3: Intercepto	-1.000	-0.631	-0.629	0.215	-1.058	-0.214
Y3: Slope v1	1.100	0.974	0.973	0.104	0.767	1.182
Y3: Slope v2	-0.800	-1.019	-1.022	0.280	-1.562	-0.464
log(phi)	2.708	2.715	2.721	0.198	2.308	3.080

8 Avaliando a qualidade do ajuste do modelo

8.1 Métricas de desempenho

Nesta seção, constam os códigos referente a amostragem da distribuição das 3 métricas de desempenho do modelo utilizadas: Distância de Aitchson, Erro quadrático Médio Padrão e Divergência de Kullback-Leibler.

# ==============================================================================
# CONFIGURAÇÃO DE DIMENSÕES E PARÂMETROS DO MCMC
# ==============================================================================
S <- nrow(cadeia_thin_df)
n <- nrow(Y_matriz)

K <- 3  # Número de preditores no modelo (Intercepto, X1, X2)
C <- 3  # Número de componentes composicionais (Y1, Y2, Y3)
d <- K * (C - 1) + 1  # Total de parâmetros amostrados (d = 7)

# Vetores para guardar as distribuições das métricas
distancias_aitchison <- numeric(S)
rmse_amostras        <- numeric(S)
kl_amostras          <- numeric(S)

# Matrizes para acumular a análise preditiva
Y_hat_acumulado <- matrix(0, nrow = n, ncol = C)
Y_sim_acumulado <- matrix(0, nrow = n, ncol = C)
Y_sim_unica     <- matrix(0, nrow = n, ncol = C) 

nomes_param_mcmc <- colnames(cadeias[[1]])
matriz_betas_loop <- as.matrix(cadeia_thin_df[, nomes_param_mcmc])

# ---------------------------------------------------------------------------
# PRÉ-COMPUTAÇÃO: Transformação CLR dos dados observados Y_matriz
# ---------------------------------------------------------------------------
log_y_estavel <- log(Y_matriz + 1e-10)
clr_y <- log_y_estavel - rowMeans(log_y_estavel)

tempo_inicio <- Sys.time()
semente_da_simulacao <- 6666

# ==============================================================================
# SUPER LOOP MESTRE
# ==============================================================================
for (s in 1:S) {
  set.seed(semente_da_simulacao + s)
  
  # PASSO A: Reconstrução dos Parâmetros da Amostra s
  vetor_params_s <- matriz_betas_loop[s, ]
  
  # 1. Betas estimados para as categorias não-referência (dimensão 3x2)
  beta_est_s <- matrix(vetor_params_s[1:(K * (C - 1))], nrow = K, ncol = C - 1)
  
  # 2. Reconstrução da matriz completa adicionando a referência (0) na coluna 1 (dimensão 3x3)
  beta_full_s <- cbind(0, beta_est_s)
  
  # 3. Extração da precisão phi
  log_phi_s <- vetor_params_s[d]
  phi_s     <- exp(log_phi_s)
  
  # MÉTRICA 1: RMSE no espaço dos 7 parâmetros estimados
  rmse_amostras[s] <- sqrt(mean((valores_reais_mcmc - vetor_params_s)^2))
  
  # PASSO B: Projeção no Simplex (Parametrização Alternativa)
  eta_s     <- X_matriz %*% beta_full_s
  exp_eta_s <- exp(eta_s)
  
  # Média esperada no simplex (Softmax)
  mu_s    <- exp_eta_s / rowSums(exp_eta_s) 
  alpha_s <- mu_s * phi_s
  
  Y_hat_s <- mu_s
  Y_hat_acumulado <- Y_hat_acumulado + Y_hat_s
  
  # -------------------------------------------------------------------------
  # MÉTRICA 2: Distância de Aitchison Vetorizada
  # -------------------------------------------------------------------------
  log_Y_hat_estavel <- log(Y_hat_s + 1e-10)
  clr_Y_hat <- log_Y_hat_estavel - rowMeans(log_Y_hat_estavel)
  distancias_aitchison[s] <- mean(sqrt(rowSums((clr_y - clr_Y_hat)^2)))
  
  # -------------------------------------------------------------------------
  # MÉTRICA 3: Divergência Kullback-Leibler
  # -------------------------------------------------------------------------
  kl_individual  <- rowSums(Y_matriz * log((Y_matriz + 1e-10) / (Y_hat_s + 1e-10)))
  kl_amostras[s] <- mean(kl_individual)
  
  # -------------------------------------------------------------------------
  # ANÁLISE PREDITIVA: Amostragem Dirichlet via Gamma Vetorizada
  # -------------------------------------------------------------------------
  gamas_sim <- matrix(rgamma(n * C, shape = alpha_s, rate = 1), nrow = n, ncol = C)
  Y_sim_s   <- gamas_sim / rowSums(gamas_sim)
  
  Y_sim_acumulado <- Y_sim_acumulado + Y_sim_s
  if (s == S) {
    Y_sim_unica <- Y_sim_s
  }
}

tempo_fim <- Sys.time()
tempo_loop <- round(difftime(tempo_fim, tempo_inicio, units = "secs"), 2)

# EXTRAÇÃO DAS MÉDIAS GLOBAIS
distancia_media_aitchison <- mean(distancias_aitchison)
rmse_global               <- mean(rmse_amostras)
kl_media                  <- mean(kl_amostras)

Y_hat_medio <- Y_hat_acumulado / S  
Y_sim_medio <- Y_sim_acumulado / S  

cat("Tempo de execução do Super Loop:", tempo_loop, "segundos\n")

## Tempo de execução do Super Loop: 0.18 segundos

8.1.1 Distância de Aitchson

8.1.2 Erro quadrático médio padrão

8.1.3 Divergência de Kullback-Leibler

8.2 Análise preditiva a posteriori

Aqui, simulados dados de acordo com os parâmetros gerados pelo modelo e comparamos com os dados observados

# 1. UNIFICAÇÃO DOS DADOS
df_sim <- data.frame(
  var1 = Y_sim_unica[, 1], 
  var2 = Y_sim_unica[, 2], 
  var3 = Y_sim_unica[, 3], 
  Tipo = "Simulado"
)

df_obs <- data.frame(
  var1 = y[, 1],           
  var2 = y[, 2],           
  var3 = y[, 3],           
  Tipo = "Observado"
)

df_combinado <- rbind(df_sim, df_obs)

# 2. GRÁFICO TERNÁRIO
plot_sobreposto <- ggtern(data = df_combinado, aes(x = var1, y = var2, z = var3, color = Tipo, shape = Tipo)) +
  geom_point(alpha = 1, size = 1.5) +
  
  scale_color_manual(values = c("Observado" = "firebrick", "Simulado" = "steelblue")) +
  scale_shape_manual(values = c("Observado" = 4, "Simulado" = 16)) + 
  
  theme_rgbw() + 
  labs(
    title = "Análise Preditiva Posterior: Dados Observados vs. Simulados",
    x = "Y1", y = "Y2", z = "Y3",
    color = "Origem", shape = "Origem"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
    legend.position = "bottom",          
    legend.title = element_text(face = "bold"),
    panel.grid.tern = element_blank()
  )

# Exibe o gráfico limpo
print(plot_sobreposto)