Estudo Simulado do modelo de Regressão Dirichlet

Abordagem Bayesiana - Parte 10

---

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 = 7\) 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)\).

\[ Y_i = (Y_{i1}, Y_{i2}, Y_{i3}, Y_{i4}, Y_{i5}, Y_{i6}, Y_{i7}) \sim \text{Dirichlet}(\alpha_{i1}, \alpha_{i2}, \alpha_{i3}, \alpha_{i4}, \alpha_{i5}, \alpha_{i6}, \alpha_{i7}) \]

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}^{7} \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}^{7} \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 = 6\) componentes restantes, totalizando \(d = 3 \times 6 + 1 = 19\) 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 7}\):

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

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

# 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 7 colunas)
beta_real <- matrix(c(
   0.5,  0.0, -0.5,  0.2, -0.4,  0.1, -0.2,  # Interceptos
  -0.8,  0.6,  0.3, -0.1,  0.5, -0.4,  0.2,  # Efeitos de X1 
   1.2, -0.9,  0.4, -0.7,  0.2,  0.8, -0.3   # 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 (7 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 * 6 + 1 = 19 parâmetros

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

# Prioris para as C - 1 = 6 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
)

# Estimadores MLE para proposta do MCMC (Dimensão d = 19)
params_mle <- as.numeric(unlist(coef(modelo_alt)))
Sigma_mle  <- vcov(modelo_alt)

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

# Nomeação dinâmica das 6 componentes não-referência
colnames_params <- c(
  as.vector(sapply(Nomes_Y[2:C], function(nm) paste0(nm, ": ", c("Intercepto", "Slope v1", "Slope v2")))),
  "log(phi)"
)

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

set.seed(1010) 

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: 23.46 % | Tempo: 66.75 s
## Cadeia 2 - Ref: Y1 | Aceitação: 23.09 % | Tempo: 66.62 s
## Cadeia 3 - Ref: Y1 | Aceitação: 22.84 % | Tempo: 61.5 s
## Cadeia 4 - Ref: Y1 | Aceitação: 22.98 % | Tempo: 58.84 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	3361	6597
Y2: Slope v1	1.001	3617	7086
Y2: Slope v2	1.002	3433	6507
Y3: Intercepto	1.001	3442	6521
Y3: Slope v1	1.001	3287	6864
Y3: Slope v2	1.002	3175	6061
Y4: Intercepto	1.001	3093	6171
Y4: Slope v1	1.003	3387	5975
Y4: Slope v2	1.001	3081	5949
Y5: Intercepto	1.001	2906	6166
Y5: Slope v1	1.001	3614	6576
Y5: Slope v2	1.001	3138	5498
Y6: Intercepto	1.001	3395	6749
Y6: Slope v1	1.002	3486	6397
Y6: Slope v2	1.001	3460	5798
Y7: Intercepto	1.002	3345	6251
Y7: Slope v1	1.001	3329	6320
Y7: Slope v2	1.001	3415	6683
log(phi)	1.001	3659	6914

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.489	-0.489	0.043	-0.572	-0.405
Y2: Slope v1	1.400	1.420	1.419	0.023	1.374	1.465
Y2: Slope v2	-2.100	-2.126	-2.126	0.062	-2.249	-2.006
Y3: Intercepto	-1.000	-1.011	-1.012	0.048	-1.105	-0.917
Y3: Slope v1	1.100	1.116	1.116	0.023	1.071	1.163
Y3: Slope v2	-0.800	-0.722	-0.721	0.062	-0.847	-0.600
Y4: Intercepto	-0.300	-0.243	-0.243	0.041	-0.325	-0.164
Y4: Slope v1	0.700	0.718	0.718	0.022	0.674	0.762
Y4: Slope v2	-1.900	-1.963	-1.963	0.064	-2.089	-1.835
Y5: Intercepto	-0.900	-0.853	-0.853	0.047	-0.945	-0.763
Y5: Slope v1	1.300	1.298	1.298	0.023	1.252	1.342
Y5: Slope v2	-1.000	-0.996	-0.996	0.061	-1.115	-0.874
Y6: Intercepto	-0.400	-0.373	-0.372	0.041	-0.455	-0.293
Y6: Slope v1	0.400	0.395	0.395	0.020	0.355	0.434
Y6: Slope v2	-0.400	-0.380	-0.381	0.052	-0.482	-0.278
Y7: Intercepto	-0.700	-0.740	-0.740	0.047	-0.831	-0.649
Y7: Slope v1	1.000	1.007	1.007	0.024	0.960	1.054
Y7: Slope v2	-1.500	-1.461	-1.461	0.064	-1.587	-1.336
log(phi)	2.708	2.723	2.723	0.021	2.680	2.764

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 <- 7  # Número de componentes (Y1 a Y7)
d <- K * (C - 1) + 1  # Total de parâmetros amostrados (d = 19)

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

# ==============================================================================
# SUPER LOOP MESTRE
# ==============================================================================
for (s in 1:S) {
  set.seed(semente_da_simulacao + s)
  
  vetor_params_s <- matriz_betas_loop[s, ]
  
  # Betas das 6 categorias não-referência (Matriz 3x6)
  beta_est_s <- matrix(vetor_params_s[1:(K * (C - 1))], nrow = K, ncol = C - 1)
  
  # Matriz completa com 0 na coluna 1 (Matriz 3x7)
  beta_full_s <- cbind(0, beta_est_s)
  
  log_phi_s <- vetor_params_s[d]
  phi_s     <- exp(log_phi_s)
  
  # 1. RMSE nos 19 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.51 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"
  )