Estudo Simulado do modelo de Regressão Dirichlet

Abordagem Bayesiana - Parte 3

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 a dados simulados. O cenário aqui considerado foi uma amostra com tamanho \(n = 750\) de vetores composicionais de dimensão 3. As covariáveis associadas, possuem, respectivamente, distribuição uniforme no conjunto \([0,10]\) e distribuição Bernoulli com parâmetro \(p = 1/2\), ou seja, geramos:\[Y_i = (Y_{i1}, Y_{i2}, Y_{i3}) \sim \text{Dirichlet}(\alpha_{i1}, \alpha_{i2}, \alpha_{i3})\]\[\alpha_{ij} = \exp(\beta_{0j} + \beta_{1j} X_{1i} + \beta_{2j} X_{2i})\]\[\beta \sim N(0, 100 I_d)\]

Para a simulação, os valores verdadeiros dos parâmetros foram organizados na matriz de coeficientes \(\boldsymbol{\beta}\), onde cada linha representa o efeito de uma covariável e cada coluna mapeia uma das componentes da resposta composicional:\[\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}\]Onde: Linha 1 (\(\beta_{0j}\)): Interceptos para as componentes \(Y_1, Y_2\) e \(Y_3\).\ Linha 2 (\(\beta_{1j}\)): Efeitos da covariável \(X_1\) sobre as componentes.\ Linha 3 (\(\beta_{2j}\)): Efeitos da covariável \(X_2\) sobre as componentes.

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(coda)

3 Geração dos dados

set.seed(2003)

# Tamanho da amostra
n <- 750

# 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
# Coluna 1 = y1 | Coluna 2 = y2 | Coluna 3 = y3
# Linha 1 = Interceptos | Linha 2 = Slopes de X1 | Linha 3 = Slopes de X2
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)

# Função de ligação
log_alpha <- X_mat %*% beta_real

# Cálculo dos alphas
alpha <- exp(log_alpha)

# Gerando as composições 
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_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")

4 Análise descritiva

5 Ajuste clássico do modelo do Maier

## Call:
## DirichReg(formula = AL ~ x1 + x2, data = dados_modelo)
## 
## Standardized Residuals:
##         Min       1Q   Median      3Q     Max
## v1  -6.0958  -0.6618   0.0797  0.8086  3.3222
## v2  -2.0036  -0.6772  -0.3793  0.4736  6.9896
## v3  -1.8228  -0.7178  -0.2997  0.5042  9.1979
## 
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 1: v1
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.63978    0.05609   11.41   <2e-16 ***
## x1          -0.86061    0.03132  -27.48   <2e-16 ***
## x2           1.21255    0.08391   14.45   <2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 2: v2
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.14520    0.05412   2.683  0.00729 ** 
## x1           0.46982    0.02942  15.972  < 2e-16 ***
## x2          -0.73160    0.07718  -9.479  < 2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 3: v3
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.43498    0.05179  -8.398  < 2e-16 ***
## x1           0.23914    0.02755   8.679  < 2e-16 ***
## x2           0.52891    0.07568   6.989 2.77e-12 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-likelihood: 2237 on 9 df (60 BFGS + 1 NR Iterations)
## AIC: -4457, BIC: -4415
## Number of Observations: 750
## Link: Log
## Parametrization: common

6 Rotina de Metropolis-Hastings

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

log_verossimilhanca <- function(beta_vec, Y, X) {
  N <- nrow(Y)
  C <- ncol(Y)
  K <- ncol(X)
  
  beta_mat <- matrix(beta_vec, nrow = K, ncol = C)
  alpha <- exp(X %*% beta_mat)
  
  termo1 <- sum(lgamma(rowSums(alpha)))
  termo2 <- sum(lgamma(alpha))
  termo3 <- sum((alpha - 1) * log(Y))
  
  ll <- termo1 - termo2 + termo3
  return(ll)
}

log_posteriori <- function(beta_vec, Y, X, sigma_prior = 100) {
  ll <- log_verossimilhanca(beta_vec, Y, X)
  log_prior <- sum(dnorm(beta_vec, mean = 0, sd = sigma_prior, log = TRUE))
  log_post <- ll + log_prior
  return(log_post)
}

6.2 Algoritmo

dados_modelo <- data.frame(x1 = x1, x2 = x2, 
                           y1 = y[, 1], y2 = y[, 2], y3 = y[, 3])
C <- 3 
K <- 3 
d <- C*K

Y_matriz <- as.matrix(dados_modelo[, c("y1", "y2", "y3")])
X_matriz <- cbind(1, dados_modelo$x1, dados_modelo$x2)

# Configurações para múltiplas cadeias
m <- 4         
nite <- 50000  

# Parâmetros de base do modelo clássico
beta_mle <- as.numeric(unlist(coef(modelo)))
Sigma_mle <- vcov(modelo)

# Fator de tuning
tuning <- 1.2
Sigma_prop <- tuning * (2.4^2 / d) * Sigma_mle

# Inicializando lista para armazenar o histórico de cada cadeia
cadeias <- list()
set.seed(1003) # Para reprodutibilidade das 4 cadeias

# Algoritmo Otimizado
for (i in 1:m) {
  
  tempo_inicio <- Sys.time()
  
  # Inicialização Superdispersa 
  beta_ini <- MASS::mvrnorm(1, mu = beta_mle, Sigma = 10 * Sigma_mle)
  
  beta_cadeia <- matrix(0, nrow = nite, ncol = d)
  beta_cadeia[1, ] <- beta_ini
  
  log_post_atual <- log_posteriori(beta_ini, Y = Y_matriz, X = X_matriz)
  
  # ------------------------------------------------------------------------
  # O PULO DO GATO: Pre-gera todo o ruído da cadeia de uma só vez
  # Fazemos a decomposição da matriz uma única vez em vez de 50.000 vezes!
  # ------------------------------------------------------------------------
  ruido_prop <- MASS::mvrnorm(nite - 1, mu = rep(0, d), Sigma = Sigma_prop)
  
  # Passo de Metropolis
  aceitos <- 0
  for (t in 1:(nite - 1)) {
    
    # Em vez de chamar MASS::mvrnorm, fazemos apenas uma soma vetorial ultra rápida
    beta_prop <- beta_cadeia[t, ] + ruido_prop[t, ]
    
    log_post_prop <- log_posteriori(beta_prop, Y = Y_matriz, X = X_matriz)
    
    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) {
      beta_cadeia[t+1, ] <- beta_prop
      log_post_atual <- log_post_prop
      aceitos <- aceitos + 1
    } else {
      beta_cadeia[t+1, ] <- beta_cadeia[t, ]
    }
  }
  
  # Armazenamos a matriz da cadeia i na lista
  cadeias[[i]] <- beta_cadeia
  
  tempo_fim <- Sys.time()
  tempo_execucao <- round(difftime(tempo_fim, tempo_inicio, units = "secs"), 2)
  
  # Imprimir resultados
  cat("Cadeia", i, 
      "- Taxa de aceitação:", round((aceitos / (nite - 1)) * 100, 2), "%",
      "| Tempo de execução:", tempo_execucao, "segundos\n")
}

## Cadeia 1 - Taxa de aceitação: 22.4 % | Tempo de execução: 12.42 segundos
## Cadeia 2 - Taxa de aceitação: 22.69 % | Tempo de execução: 12.49 segundos
## Cadeia 3 - Taxa de aceitação: 22.95 % | Tempo de execução: 12.35 segundos
## Cadeia 4 - Taxa de aceitação: 22.56 % | Tempo de execução: 12.39 segundos

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
y1: Intercepto	1.002	2988	5045
y1: Slope v1	1.002	3516	6136
y1: Slope v2	1.001	3052	5705
y2: Intercepto	1.002	3203	5112
y2: Slope v1	1.002	3772	6739
y2: Slope v2	1.001	3210	5990
y3: Intercepto	1.001	3634	5852
y3: Slope v1	1.002	3737	6173
y3: Slope v2	1.003	3688	5658

6.4 Inferência Bayesiana

6.4.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 Intervalos de 95% de Credibilidade
Parâmetro	Valor Real	Média Post.	Mediana	Desvio Padrão	Quantil 2.5%	Quantil 97.5%
y1: Intercepto	0.5	0.545	0.546	0.057	0.432	0.654
y1: Slope v1	-0.8	-0.787	-0.787	0.032	-0.850	-0.725
y1: Slope v2	1.2	1.119	1.119	0.086	0.952	1.289
y2: Intercepto	0.0	0.086	0.087	0.054	-0.020	0.189
y2: Slope v1	0.6	0.616	0.616	0.027	0.563	0.670
y2: Slope v2	-0.9	-1.017	-1.018	0.076	-1.162	-0.868
y3: Intercepto	-0.5	-0.508	-0.508	0.052	-0.612	-0.408
y3: Slope v1	0.3	0.323	0.323	0.026	0.270	0.375
y3: Slope v2	0.4	0.447	0.446	0.075	0.302	0.597

7 Avaliando a qualidade do ajuste do modelo

7.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.

S <- nrow(cadeia_thin_df)
n <- nrow(y)
K <- 3  # Número de covariáveis
C <- 3  # Número de componentes

# 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) 

matriz_betas_loop <- as.matrix(cadeia_thin_df[, nomes_param])

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

# CRONÔMETRO: Dispara o relógio antes de entrar no loop
tempo_inicio <- Sys.time()

semente_da_simulacao <- 3333

# --- O SUPER LOOP MESTRE 
for (s in 1:S) {
  set.seed(semente_da_simulacao + s)
  
  # PASSO A: Reconstrução dos Parâmetros
  vetor_betas_s <- matriz_betas_loop[s, ]
  beta_s        <- matrix(vetor_betas_s, nrow = K, ncol = C, byrow = FALSE)
  
  # MÉTRICA 1: RMSE (Espaço dos Parâmetros)
  rmse_amostras[s] <- sqrt(mean((beta_real - beta_s)^2))
  
  # PASSO B: Projeção no Simplex
  alpha_s <- exp(X_mat %*% beta_s)
  Y_hat_s <- alpha_s / rowSums(alpha_s)
  
  # Acumula para calcular a tendência média no final
  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 * log((y + 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)
  
  # Acumula a simulação e guarda a última para o gráfico de dispersão
  Y_sim_acumulado <- Y_sim_acumulado + Y_sim_s
  if (s == S) {
    Y_sim_unica <- Y_sim_s
  }
}

# CRONÔMETRO: Para o relógio imediatamente após o fim do loop
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  

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

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

7.1.1 Distância de Aitchson

7.1.2 Erro quadrático médio padrão

7.2 Divergência de Kullback-Leibler

7.3 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)