Estudo Simulado do modelo de Regressão Dirichlet

Abordagem Bayesiana - Parte 4

---

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 5. 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}, Y_{i4}, Y_{i5}) \sim \text{Dirichlet}(\alpha_{i1}, \alpha_{i2}, \alpha_{i3}, \alpha_{i4}, \alpha_{i5})\]

\[\alpha_{ij} = \exp(\beta_{0j} + \beta_{1j} X_{1i} + \beta_{2j} X_{2i}), \quad \text{para } j = 1, \dots, 5\]

\[\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_{04} & \beta_{05} \\ \beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} & \beta_{15} \\ \beta_{21} & \beta_{22} & \beta_{23} & \beta_{24} & \beta_{25} \end{bmatrix} = \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}\]

Onde: * Linha 1 (\(\beta_{0j}\)): Interceptos para as componentes \(Y_1, Y_2, Y_3, Y_4\) e \(Y_5\). * 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(2004)

# 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 (Agora com 5 colunas)
# Linha 1 = Interceptos | Linha 2 = Slopes de X1 | Linha 3 = Slopes de X2
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)

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

for(i in 1:n) {
  # Sorteia 5 gamas em vez de 3
  z <- rgamma(5, shape = alpha[i, ], rate = 1)
  
  # Fecha no simplex (soma 1)
  y[i, ] <- z / sum(z)
}

# Expandindo os nomes dos parâmetros para as 5 dimensões (15 parâmetros no total)
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",
                 "y4: Intercepto", "y4: Slope v1", "y4: Slope v2",
                 "y5: Intercepto", "y5: Slope v1", "y5: 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  -4.9783  -0.6715  -0.0209  0.6556  2.9596
## v2  -1.8580  -0.6470  -0.3865  0.2881  5.3855
## v3  -1.4260  -0.7303  -0.3296  0.3387  5.9248
## v4  -1.3778  -0.7587  -0.3060  0.6143  4.9006
## v5  -1.6584  -0.6920  -0.3752  0.4531  7.3423
## 
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 1: v1
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.54538    0.04727   11.54   <2e-16 ***
## x1          -0.82574    0.02357  -35.03   <2e-16 ***
## x2           1.22220    0.06555   18.64   <2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 2: v2
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.04019    0.04726    0.85    0.395    
## x1           0.44446    0.02433   18.27   <2e-16 ***
## x2          -0.71337    0.06789  -10.51   <2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 3: v3
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.39864    0.04890  -8.152 3.58e-16 ***
## x1           0.25123    0.02427  10.351  < 2e-16 ***
## x2           0.40859    0.06774   6.032 1.62e-09 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 4: v4
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.18935    0.04739   3.996 6.45e-05 ***
## x1          -0.12431    0.02336  -5.322 1.02e-07 ***
## x2          -0.57651    0.06751  -8.540  < 2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 5: v5
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.26693    0.04826  -5.531 3.19e-08 ***
## x1           0.41064    0.02391  17.172  < 2e-16 ***
## x2           0.14378    0.06727   2.137   0.0326 *  
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-likelihood: 5005 on 15 df (94 BFGS + 2 NR Iterations)
## AIC: -9980, BIC: -9911
## 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],
                           y4 = y[, 4], y5 = y[, 5])
C <- 5 
K <- 3 
d <- C*K

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

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

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

# Fator de tuning ajustado para maior estabilidade em 15D
tuning <- 1
Sigma_prop <- tuning * (2.4^2 / d) * Sigma_mle

# Inicializando lista para armazenar o histórico de cada cadeia
cadeias <- list()
set.seed(1004) 

# Algoritmo Otimizado
for (i in 1:m) {
  
  tempo_inicio <- Sys.time()
  
  # Inicialização ajustada para evitar regiões de densidade nula em alta dimensão
  beta_ini <- MASS::mvrnorm(1, mu = beta_mle, Sigma = 3 * 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)
  
  ruido_prop <- MASS::mvrnorm(nite - 1, mu = rep(0, d), Sigma = Sigma_prop)
  
  # Passo de Metropolis
  aceitos <- 0
  for (t in 1:(nite - 1)) {
    
    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, ]
    }
  }
  
  cadeias[[i]] <- beta_cadeia
  tempo_fim <- Sys.time()
  tempo_execucao <- round(difftime(tempo_fim, tempo_inicio, units = "secs"), 2)
  
  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: 25.27 % | Tempo de execução: 28.84 segundos
## Cadeia 2 - Taxa de aceitação: 25.46 % | Tempo de execução: 28.94 segundos
## Cadeia 3 - Taxa de aceitação: 25.31 % | Tempo de execução: 28.59 segundos
## Cadeia 4 - Taxa de aceitação: 25.41 % | Tempo de execução: 28.47 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.001	3326	5886
y1: Slope v1	1.002	3500	6744
y1: Slope v2	1.001	3367	6500
y2: Intercepto	1.001	3629	6822
y2: Slope v1	1.001	3876	7846
y2: Slope v2	1.001	3605	6643
y3: Intercepto	1.001	3604	6624
y3: Slope v1	1.002	3760	7001
y3: Slope v2	1.001	3673	6917
y4: Intercepto	1.002	3686	6789
y4: Slope v1	1.003	3229	7362
y4: Slope v2	1.002	3549	6888
y5: Intercepto	1.002	3300	6552
y5: Slope v1	1.001	3608	6498
y5: Slope v2	1.001	3553	6224

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.441	0.441	0.048	0.345	0.533
y1: Slope v1	-0.8	-0.774	-0.774	0.024	-0.822	-0.727
y1: Slope v2	1.2	1.192	1.192	0.068	1.060	1.325
y2: Intercepto	0.0	-0.067	-0.066	0.047	-0.160	0.026
y2: Slope v1	0.6	0.570	0.570	0.023	0.524	0.615
y2: Slope v2	-0.9	-0.873	-0.872	0.068	-1.005	-0.739
y3: Intercepto	-0.5	-0.491	-0.490	0.049	-0.588	-0.395
y3: Slope v1	0.3	0.328	0.328	0.024	0.281	0.376
y3: Slope v2	0.4	0.373	0.373	0.068	0.242	0.507
y4: Intercepto	0.2	0.124	0.125	0.047	0.029	0.215
y4: Slope v1	-0.1	-0.079	-0.079	0.023	-0.124	-0.033
y4: Slope v2	-0.7	-0.635	-0.635	0.068	-0.769	-0.500
y5: Intercepto	-0.4	-0.376	-0.376	0.048	-0.471	-0.281
y5: Slope v1	0.5	0.512	0.512	0.023	0.466	0.557
y5: Slope v2	0.2	0.093	0.094	0.067	-0.040	0.224

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 parâmetros
C <- 5  # 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 <- 4444

# --- O SUPER LOOP MESTRE (VERSÃO VETORIZADA - SEM LOOPS INTERNOS)
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: 1.09 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. Garante a nomeação correta para as 7 componentes atuais
colnames(y) <- paste0("Y", 1:5)
colnames(Y_sim_unica) <- paste0("Y", 1:5)

# 2. Estrutura os dados observados e simulados
df_obs <- as.data.frame(y) %>% mutate(Tipo = "Observado")
df_sim <- as.data.frame(Y_sim_unica) %>% mutate(Tipo = "Simulado")

# 3. Unifica e transforma para o formato longo (long format)
df_ppc_long <- rbind(df_obs, df_sim) %>%
  pivot_longer(
    cols = starts_with("Y"),
    names_to = "Componente",
    values_to = "Proporcao"
  )

# 4. Gráfico de Distribuição Facetado (Violin + Boxplot)
grafico_ppc_dist <- ggplot(df_ppc_long, aes(x = Tipo, y = Proporcao, fill = Tipo)) +
  # O violin mostra o formato da densidade dos dados
  geom_violin(alpha = 0.5, color = NA, position = position_dodge(0.8)) +
  # O boxplot traz os quartis e mediana de forma direta
  geom_boxplot(width = 0.22, color = "grey20", outlier.alpha = 0.1, 
               position = position_dodge(0.8), lwd = 0.4) +
  
  # Faceta para cada uma das 7 componentes
  facet_wrap(~ Componente, ncol = 4, scales = "free_y") +
  
  # Cores com alto contraste e fáceis de ler
  scale_fill_manual(values = c("Observado" = "#d95f02", "Simulado" = "#7570b3")) +
  
  theme_minimal(base_size = 11) +
  labs(
    title = "Checagem Preditiva Posterior (PPC)",
    subtitle = "Comparação das distribuições das proporções observadas vs. simuladas",
    x = NULL,
    y = "Proporção (Simplex)",
    fill = "Origem dos Dados"
  ) +
  theme(
    plot.title = element_text(face = "bold", size = 13, hjust = 0.5),
    plot.subtitle = element_text(size = 10, color = "grey40", hjust = 0.5),
    strip.background = element_rect(fill = "gray95", color = NA),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

print(grafico_ppc_dist)