Estudo Simulado do modelo de Regressão Dirichlet

Abordagem Bayesiana - Parte 9

Pedro Frazão Dutra

22/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) aplicado a dados simulados com \(C = 5\) componentes composicionais. O cenário considerado consiste em uma amostra de tamanho \(n = 750\). 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}, Y_{i4}, Y_{i5}) \sim \text{Dirichlet}(\alpha_{i1}, \alpha_{i2}, \alpha_{i3}, \alpha_{i4}, \alpha_{i5}) \]

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}^{5} \mu_{ij} = 1 \] A expectância \(\boldsymbol{\mu}_i\) é mapeada no simplex através da função Softmax:

\[ \mu_{ij} = \frac{\exp(\eta_{ij})}{\sum_{k=1}^{5} \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}\)). O algoritmo MCMC estima os coeficientes relativos \(\boldsymbol{\beta}_j^* = \boldsymbol{\beta}_j - \boldsymbol{\beta}_{c_{\text{ref}}}\) para as \(C - 1 = 4\) componentes restantes, totalizando \(d = 3 \times 4 + 1 = 13\) parâmetros no vetor amostrado.Os valores verdadeiros utilizados na simulação foram \(\phi = 15\) (\(\log(\phi) \approx 2{,}708\)) e os coeficientes contidos na matriz \(\boldsymbol{\beta}_{3 \times 5}\):

\[ \boldsymbol{\beta} = \begin{bmatrix} \phantom{-}0{,}5 & \phantom{-}0{,}0 & -0{,}5 & \phantom{-}0{,}2 & -0{,}4 \\ -0{,}8 & \phantom{-}0{,}6 & \phantom{-}0{,}3 & -0{,}1 & \phantom{-}0{,}5 \\ \phantom{-}1{,}2 & -0{,}9 & \phantom{-}0{,}4 & -0{,}7 & \phantom{-}0{,}2 \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(2009)

# Tamanho da amostra e número de componentes
n <- 750
C <- 5

# 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 (3 linhas x 5 colunas)
beta_real <- matrix(c(
   0.5,  0.0, -0.5,  0.2, -0.4,  # Interceptos
  -0.8,  0.6,  0.3, -0.1,  0.5,  # Efeitos de X1 
   1.2, -0.9,  0.4, -0.7,  0.2   # Efeitos de X2 
), nrow = 3, byrow = TRUE)

# Mapeamento do preditor linear para a média mu no simplex (Softmax)
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 (5 componentes)
y <- matrix(0, nrow = n, ncol = C)

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

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                    # Intercepto, X1, X2
d <- K * (C - 1) + 1      # d = 3 * 4 + 1 = 13 parâmetros

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

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
)

params_mle <- as.numeric(unlist(coef(modelo_alt)))
Sigma_mle  <- vcov(modelo_alt)

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

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"),
  paste0(Nomes_Y[4], ": Intercepto"), paste0(Nomes_Y[4], ": Slope v1"), paste0(Nomes_Y[4], ": Slope v2"),
  paste0(Nomes_Y[5], ": Intercepto"), paste0(Nomes_Y[5], ": Slope v1"), paste0(Nomes_Y[5], ": Slope v2"),
  "log(phi)"
)

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

set.seed(1009) 

for (i in 1:m) {
  tempo_inicio <- Sys.time()
  
  params_ini <- MASS::mvrnorm(1, mu = params_mle, Sigma = 10 * 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: 22.91 % | Tempo: 34.42 s
## Cadeia 2 - Ref: Y1 | Aceitação: 22.85 % | Tempo: 33.54 s
## Cadeia 3 - Ref: Y1 | Aceitação: 22.8 % | Tempo: 33.8 s
## Cadeia 4 - Ref: Y1 | Aceitação: 22.89 % | Tempo: 33.34 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	3585	6350
Y2: Slope v1	1.001	4005	7247
Y2: Slope v2	1.002	3756	6853
Y3: Intercepto	1.001	3739	6811
Y3: Slope v1	1.002	3882	7209
Y3: Slope v2	1.001	3671	6319
Y4: Intercepto	1.001	4007	6706
Y4: Slope v1	1.001	3755	6904
Y4: Slope v2	1.001	3866	6988
Y5: Intercepto	1.001	3762	7336
Y5: Slope v1	1.002	3720	7001
Y5: Slope v2	1.001	3982	6925
log(phi)	1.001	4174	7646

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%
Parâmetro	Valor Real Relativo	Média Post.	Mediana	Desvio Padrão	Quantil 2.5%	Quantil 97.5%
Y2: Intercepto	-0.500	-0.492	-0.492	0.040	-0.571	-0.413
Y2: Slope v1	1.400	1.391	1.391	0.022	1.348	1.434
Y2: Slope v2	-2.100	-2.128	-2.128	0.059	-2.245	-2.012
Y3: Intercepto	-1.000	-0.977	-0.977	0.044	-1.064	-0.890
Y3: Slope v1	1.100	1.114	1.114	0.023	1.070	1.159
Y3: Slope v2	-0.800	-0.803	-0.803	0.057	-0.914	-0.688
Y4: Intercepto	-0.300	-0.302	-0.302	0.037	-0.375	-0.230
Y4: Slope v1	0.700	0.726	0.726	0.022	0.683	0.768
Y4: Slope v2	-1.900	-1.894	-1.894	0.058	-2.007	-1.780
Y5: Intercepto	-0.900	-0.929	-0.929	0.044	-1.016	-0.843
Y5: Slope v1	1.300	1.327	1.327	0.022	1.283	1.371
Y5: Slope v2	-1.000	-0.951	-0.950	0.057	-1.061	-0.838
log(phi)	2.708	2.731	2.731	0.027	2.679	2.784

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
# ==================================================================
S <- nrow(cadeia_thin_df)
n <- nrow(Y_matriz)

K <- 3  # Número de preditores (Intercepto, X1, X2)
C <- 5  # Número de componentes (Y1 a Y5)
d <- K * (C - 1) + 1  # Total de parâmetros amostrados (d = 13)

# Vetores para acumular métricas
distancias_aitchison <- numeric(S)
rmse_amostras        <- numeric(S)
kl_amostras          <- numeric(S)

# Matrizes acumuladoras
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_mcmc])

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

# ==============================================================================
# SUPER LOOP MESTRE
# ==============================================================================
for (s in 1:S) {
  set.seed(semente_da_simulacao + s)
  
  vetor_params_s <- matriz_betas_loop[s, ]
  
  # Betas das 4 categorias não-referência (Matriz 3x4)
  beta_est_s <- matrix(vetor_params_s[1:(K * (C - 1))], nrow = K, ncol = C - 1)
  
  # Matriz completa com 0 na coluna 1 (Matriz 3x5)
  beta_full_s <- cbind(0, beta_est_s)
  
  log_phi_s <- vetor_params_s[d]
  phi_s     <- exp(log_phi_s)
  
  # 1. RMSE nos 13 parâmetros estimados
  rmse_amostras[s] <- sqrt(mean((valores_reais_mcmc - vetor_params_s)^2))
  
  # 2. Média esperada no simplex (Softmax)
  eta_s     <- X_matriz %*% beta_full_s
  exp_eta_s <- exp(eta_s)
  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
  
  # 3. Distância de Aitchison
  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)))
  
  # 4. Divergência KL
  kl_individual  <- rowSums(Y_matriz * log((Y_matriz + 1e-10) / (Y_hat_s + 1e-10)))
  kl_amostras[s] <- mean(kl_individual)
  
  # 5. Amostragem Preditiva via Gamma
  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)

distancia_media_aitchison <- mean(distancias_aitchison)
rmse_global               <- mean(rmse_amostras)
kl_media                  <- mean(kl_amostras)

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

## Tempo de execução do Super Loop: 1.15 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

# Comparação em formato longo para os 5 componentes
df_sim_long <- as.data.frame(Y_sim_unica) %>%
  setNames(Nomes_Y) %>%
  mutate(Tipo = "Simulado") %>%
  pivot_longer(cols = -Tipo, names_to = "Componente", values_to = "Proporcao")

df_obs_long <- as.data.frame(Y_matriz) %>%
  setNames(Nomes_Y) %>%
  mutate(Tipo = "Observado") %>%
  pivot_longer(cols = -Tipo, names_to = "Componente", values_to = "Proporcao")

df_ppc <- rbind(df_sim_long, df_obs_long)

ggplot(df_ppc, aes(x = Componente, y = Proporcao, fill = Tipo)) +
  geom_boxplot(alpha = 0.7, outlier.size = 1) +
  scale_fill_manual(values = c("Observado" = "firebrick", "Simulado" = "steelblue")) +
  theme_minimal(base_size = 14) +
  labs(
    x = "Componentes",
    y = "Proporção",
    fill = "Origem dos Dados"
  ) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    legend.position = "bottom"
  )