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

library(DirichletReg) 
library(compositions) 
library(colorspace)
library(MASS)
library(purrr)
library(coda)
library(patchwork)
library(tidyverse)   
library(DT)          
library(knitr)       
library(kableExtra)  
library(htmltools)   
library(plotly)      
library(scales)      
library(posterior)
library(bayesplot)

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

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

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(2004) # Para reprodutibilidade das 4 cadeias

# Algoritmo
for (i in 1:m) {
  # 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)
  
  # Passo de Metropolis
  aceitos <- 0
  for (t in 1:(nite - 1)) {
    beta_prop <- MASS::mvrnorm(1, mu = beta_cadeia[t, ], Sigma = Sigma_prop)
    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
  
  # Imprimir taxa de aceitação de cada cadeia
  cat("Cadeia", i, "- Taxa de aceitação:", round((aceitos / (nite - 1)) * 100, 2), "%\n")
}
## Cadeia 1 - Taxa de aceitação: 25.51 %
## Cadeia 2 - Taxa de aceitação: 25.33 %
## Cadeia 3 - Taxa de aceitação: 25.4 %
## Cadeia 4 - Taxa de aceitação: 25.51 %

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 Robustos de Convergência
Parâmetro R-Hat ESS Bulk ESS Tail
y1: Intercepto 1.000 4229 7908
y1: Slope v1 1.001 4054 7344
y1: Slope v2 1.001 4066 8274
y2: Intercepto 1.000 4667 7732
y2: Slope v1 1.001 4277 8320
y2: Slope v2 1.000 4460 8445
y3: Intercepto 1.001 4363 8421
y3: Slope v1 1.001 4412 7751
y3: Slope v2 1.001 4281 8511
y4: Intercepto 1.001 4567 8854
y4: Slope v1 1.001 4028 7508
y4: Slope v2 1.001 4291 7959
y5: Intercepto 1.001 4454 8652
y5: Slope v1 1.001 4516 9123
y5: Slope v2 1.001 4361 8068

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.442 0.442 0.049 0.345 0.535
y1: Slope v1 -0.8 -0.773 -0.773 0.024 -0.820 -0.726
y1: Slope v2 1.2 1.191 1.191 0.068 1.057 1.325
y2: Intercepto 0.0 -0.067 -0.065 0.048 -0.162 0.024
y2: Slope v1 0.6 0.570 0.570 0.024 0.523 0.616
y2: Slope v2 -0.9 -0.873 -0.873 0.068 -1.005 -0.741
y3: Intercepto -0.5 -0.491 -0.491 0.049 -0.586 -0.396
y3: Slope v1 0.3 0.329 0.329 0.024 0.281 0.376
y3: Slope v2 0.4 0.374 0.374 0.069 0.237 0.508
y4: Intercepto 0.2 0.123 0.124 0.048 0.026 0.215
y4: Slope v1 -0.1 -0.079 -0.079 0.024 -0.125 -0.032
y4: Slope v2 -0.7 -0.634 -0.634 0.069 -0.771 -0.497
y5: Intercepto -0.4 -0.376 -0.376 0.048 -0.471 -0.284
y5: Slope v1 0.5 0.512 0.512 0.023 0.467 0.558
y5: Slope v2 0.2 0.094 0.094 0.068 -0.038 0.226

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_total <- nrow(cadeia_final_df)
tamanho_alvo <- 2000

# Cria uma sequência de 2000 índices espaçados de forma equidistante
indices_thinning <- round(seq(1, S_total, length.out = tamanho_alvo))

# Extrai apenas as 2000 amostras selecionadas
cadeia_thin_df <- cadeia_final_df[indices_thinning, ]

# 1. Configurações Iniciais
S <- nrow(cadeia_thin_df)
n <- nrow(y)
K <- 3  # Número de covariáveis
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) 

# --- O SUPER LOOP MESTRE ---
for (s in 1:S) {
  
  # PASSO A: Reconstrução dos Parâmetros da iteração 's'
  vetor_betas_s <- as.numeric(cadeia_thin_df[s, nomes_param])
  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 (Os Y ajustados que você queria reaproveitar!)
  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 (Geometria do Simplex)
  # ----------------------------------------------------------------------------
  dist_n <- numeric(n)
  for(i in 1:n) {
    dist_n[i] <- norm(acomp(y[i, ]) - acomp(Y_hat_s[i, ]))
  }
  distancias_aitchison[s] <- mean(dist_n)
  
  # ----------------------------------------------------------------------------
  # MÉTRICA 3: Divergência Kullback-Leibler (Teoria da Informação)
  # ----------------------------------------------------------------------------
  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
  # ----------------------------------------------------------------------------
  Y_sim_s <- matrix(0, nrow = n, ncol = C)
  for (i in 1:n) {
    gamas_sim <- rgamma(C, shape = alpha_s[i, ], rate = 1)
    Y_sim_s[i, ] <- gamas_sim / sum(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
  }
}

# --- FIM DO LOOP: 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  # Matriz média dos valores esperados (Linha de tendência)
Y_sim_medio <- Y_sim_acumulado / S  # Matriz média das predições

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

par(mfrow = c(1, 2))

plot(acomp(Y_sim_unica), col = "black", pch = 16, cex = 0.6, 
     main = "Dados Simulados (Y_sim)")

plot(acomp(y), col = rgb(1, 0, 0, 0.6), 
     pch = 4, cex = 0.6, 
     main = "Dados Observados Real (Y)")