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

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

# Tamanho da amostra
n <- 30

# 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  -3.8008  -0.4976  -0.0242  0.9058  1.6252
## v2  -1.7917  -0.6980  -0.4717  0.4601  1.9766
## v3  -1.2532  -0.7638  -0.2860  0.5221  4.6809
## 
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 1: v1
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.4815     0.2900   1.660 0.096880 .  
## x1           -0.4712     0.1591  -2.961 0.003065 ** 
## x2            1.7349     0.4746   3.655 0.000257 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 2: v2
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.2113     0.2563   0.825    0.410    
## x1            0.6283     0.1401   4.483 7.37e-06 ***
## x2           -0.4618     0.4095  -1.128    0.259    
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 3: v3
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.2315     0.2446  -0.946   0.3440  
## x1            0.4408     0.1836   2.400   0.0164 *
## x2            0.5725     0.3993   1.434   0.1516  
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-likelihood: 84.83 on 9 df (43 BFGS + 1 NR Iterations)
## AIC: -151.7, BIC: -139
## Number of Observations: 30
## 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
Sigma_prop <- tuning * (2.4^2 / d) * Sigma_mle

# Inicializando lista para armazenar o histórico de cada cadeia
cadeias <- list()
set.seed(2001) # 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: 26.59 %
## Cadeia 2 - Taxa de aceitação: 26.26 %
## Cadeia 3 - Taxa de aceitação: 26.48 %
## Cadeia 4 - Taxa de aceitação: 26.71 %

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.001 3166 4491
y1: Slope v1 1.002 2795 4903
y1: Slope v2 1.001 3109 5055
y2: Intercepto 1.000 3100 5026
y2: Slope v1 1.001 3324 6235
y2: Slope v2 1.001 3163 4923
y3: Intercepto 1.001 3278 5588
y3: Slope v1 1.002 3440 5242
y3: Slope v2 1.002 3424 5848

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.404 0.421 0.299 -0.231 0.953
y1: Slope v1 -0.8 -0.476 -0.478 0.166 -0.793 -0.141
y1: Slope v2 1.2 1.652 1.655 0.489 0.678 2.587
y2: Intercepto 0.0 0.135 0.148 0.269 -0.427 0.628
y2: Slope v1 0.6 0.620 0.618 0.141 0.349 0.904
y2: Slope v2 -0.9 -0.515 -0.513 0.431 -1.375 0.313
y3: Intercepto -0.5 -0.312 -0.303 0.251 -0.833 0.155
y3: Slope v1 0.3 0.430 0.429 0.187 0.061 0.800
y3: Slope v2 0.4 0.508 0.512 0.408 -0.306 1.292

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

set.seed(100)

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