Modelo de Regressão Linear Múltipla - Ajustando MRLM sob a abordagem bayesiana

1 Simulando um conjunto de dados

Simule um conjunto de dados considerando n=1000 (tamanho amostral), \(\beta_0=2\), \(\beta_1=5\), \(\beta_2=3\), \(\sigma^2=1\) e que \(X_{i1} \stackrel{iid}{\sim} N(0,1)\), \(X_{i2} \stackrel{iid}{\sim} t_{(5)}\) e o seguinte modelo \[Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2}, \sigma^2).\]

# Configurações iniciais
n   = 1000  # Número de observações
# Simulando as covariáveis
x1  = rnorm(n, mean = 0, sd = sqrt(1))   # Covariável 1
x2  = rt(n, 5)   # Covariável 2
x   = cbind(1, x1=x1, x2=x2)

#definindo total de covariaveis
p  = ncol(x)-1

# Definindo os coeficientes reais
beta_0 = 2    # Intercepto
beta_1 = 5     # Coeficiente para x1
beta_2 = 3     # Coeficiente para x2
beta   = c(beta_0, beta_1, beta_2)

# Gerando a variável resposta com erro aleatório
sigma2  = 1
e       = rnorm(n, mean = 0, sd = sqrt(sigma2))  # Erro aleatório
y       = x %*% beta + e

# Criando um data frame com os dados
dados   = data.frame(y, x1=x1, x2=x2)

2 Analisando os dados

Analise o conjunto de dados simulado.

#analisando os dados gerados
par(mfrow=c(1,2))
hist(y, main="", freq=F, bty="n", ylab="densidade")
hist(e, main="", freq=F, bty="n", ylab="densidade")

par(mfrow=c(1,3))
plot(x1, y, bty="n")
plot(x2, y, bty="n")
plot(x1, x2, bty="n")

round(cor(dados),3)
##        y     x1     x2
## y  1.000  0.765  0.617
## x1 0.765  1.000 -0.008
## x2 0.617 -0.008  1.000

3 Ajuste usando rstanarm

Considere agora que os parâmetros são desconhecidos e independentes. Supondo que \(\beta_0 \sim N(0,1000)\), \(\beta_1 \sim N(0,1000)\), \(\beta_2 \sim N(0,1000)\) e \(\sigma^2 \sim Exp(1)\), ajuste o modelo usando o pacote rstanarm.

if(!require("rstanarm")) {install.packages("rstanarm")}
if(!require("parallel")) {install.packages("parallel")}
library(rstanarm)
library(parallel)

# Marcando o tempo inicial
tempo_inicial <- Sys.time()
# Fazendo o ajuste
ajusteBay_Stanarm = stan_glm(y ~ x1 + x2, 
                     data = dados,
                     family  = gaussian(), 
                     warmup = 500, chain = 2, iter = 2000, 
                     prior = normal(0, 1000, autoscale = FALSE),
                     prior_intercept = normal(0, 1000, autoscale = FALSE),
                     prior_aux = exponential(1) ,
                     refresh = 0  # Suprime a saída de progresso 
)

# Marcando o tempo final
tempo_final <- Sys.time()
# Calculando o tempo gasto
tempo_final - tempo_inicial
## Time difference of 4.890376 secs

3.1 Estimativas dos parâmetros

Obtenha as médias a posteriori, o desvio padrão e os quantis amostrais de 10%, 50% e 90%. Compare com os valores obtidos pelo summary do ajuste.

# Extraindo as amostras sem diferenciar a cadeia
amostras_matriz = as.matrix(ajusteBay_Stanarm)  # linhas = iterações, colunas = parâmetros
X = model.matrix(~ x1 + x2, data = dados)  # inclui intercepto automaticamente
amostra_Stanarm_beta  = amostras_matriz[, colnames(X)]  # seleciona Intercept, x1, x2
amostra_Stanarm_sigma = amostras_matriz[, "sigma"]
#Calculando as medidas resumos
round(
  rbind(
  c(mean(amostra_Stanarm_beta[,1]), sd(amostra_Stanarm_beta[,1]), quantile(amostra_Stanarm_beta[,1], c(0.1, 0.5, 0.9))),
  c(mean(amostra_Stanarm_beta[,2]), sd(amostra_Stanarm_beta[,2]), quantile(amostra_Stanarm_beta[,2], c(0.1, 0.5, 0.9))),
  c(mean(amostra_Stanarm_beta[,3]), sd(amostra_Stanarm_beta[,3]), quantile(amostra_Stanarm_beta[,3], c(0.1, 0.5, 0.9))),
  c(mean(amostra_Stanarm_sigma), sd(amostra_Stanarm_sigma), quantile(amostra_Stanarm_sigma, c(0.1, 0.5, 0.9))))
,1)
##            10% 50% 90%
## [1,] 2.1 0 2.0 2.1 2.1
## [2,] 5.0 0 4.9 5.0 5.0
## [3,] 3.1 0 3.0 3.1 3.1
## [4,] 1.0 0 1.0 1.0 1.1
summary(ajusteBay_Stanarm)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      y ~ x1 + x2
##  algorithm:    sampling
##  sample:       3000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 1000
##  predictors:   3
## 
## Estimates:
##               mean   sd   10%   50%   90%
## (Intercept) 2.1    0.0  2.0   2.1   2.1  
## x1          5.0    0.0  4.9   5.0   5.0  
## x2          3.1    0.0  3.0   3.1   3.1  
## sigma       1.0    0.0  1.0   1.0   1.1  
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 2.2    0.0  2.2   2.2   2.3  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0   780 
## x1            0.0  1.0  3214 
## x2            0.0  1.0  2918 
## sigma         0.0  1.0   794 
## mean_PPD      0.0  1.0  1198 
## log-posterior 0.0  1.0   883 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

3.2 Média da distribuição preditiva a posteriori

Após o ajuste de um modelo bayesiano, utilizamos os parâmetros estimados para gerar possíveis valores da variável resposta. Sendo assim, para cada iteração \(l\) obtida, calculamos \[\hat{Y}_i^{(l)} = \hat{\beta}_0^{(l)} + \hat{\beta}_1^{(l)} X_{i1} + \hat{\beta}_2^{(l)} X_{i2} + \hat{e}_i^{(l)},\] sendo \(\hat{e}_i^{(l)} \stackrel{iid}{\sim} N\left(0 \;, \; (\hat{\sigma}^2)^{(l)}\right).\)

A média dessas previsões, \(\frac{\hat{Y}_i^{(1)}+\ldots + \hat{Y}_i^{(Nite)}}{Nite}\), sendo \(Nite\) o tamanho da amostra obtida para cada parâmetro, é chamada de média da distribuição preditiva a posteriori. Essa métrica serve como um diagnóstico rápido para verificar se o modelo consegue reproduzir, ao menos, a média dos dados observados. A ideia é simples:

  • Se a média da distribuição preditiva a posteriori for próxima da média amostral da variável resposta \(\bar{Y}\), isso sugere que o modelo consegue capturar o comportamento médio da variável resposta.
  • Se for muito diferente, isso indica um problema, como:
    • Erros de especificação do modelo.
    • Problemas com os dados (outliers, erro de medição).
    • Problemas computacionais (não convergência do MCMC, por exemplo).

Essa comparação não garante que o modelo seja bom, mas ajuda a detectar falhas grosseiras.

Sendo assim, calcule essa medida manualmente e compare com os valores obtidos pelo summary do ajuste.

# Calculando o valor medio previsto para cada iteração e observação
mu = amostra_Stanarm_beta %*% t(X)  # multiplicação de matriz

# Obtendo o y predito 
nite  = nrow(mu)
epred = matrix(rnorm(nite * n, 0, amostra_Stanarm_sigma), nite, n)
ypred = mu + epred

# Calculando a media final
ypred_medio = rowMeans(ypred) # média de cada iteração (linha)
c(round(mean(ypred_medio),1),
round(sd(ypred_medio),1),
round(quantile(ypred_medio,c(0.1, 0.5, 0.9)),1))
##         10% 50% 90% 
## 2.2 0.0 2.2 2.2 2.3
summary(ajusteBay_Stanarm)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      y ~ x1 + x2
##  algorithm:    sampling
##  sample:       3000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 1000
##  predictors:   3
## 
## Estimates:
##               mean   sd   10%   50%   90%
## (Intercept) 2.1    0.0  2.0   2.1   2.1  
## x1          5.0    0.0  4.9   5.0   5.0  
## x2          3.1    0.0  3.0   3.1   3.1  
## sigma       1.0    0.0  1.0   1.0   1.1  
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 2.2    0.0  2.2   2.2   2.3  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0   780 
## x1            0.0  1.0  3214 
## x2            0.0  1.0  2918 
## sigma         0.0  1.0   794 
## mean_PPD      0.0  1.0  1198 
## log-posterior 0.0  1.0   883 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
cat("Média observada: ", round(mean(y), 2), "\n")
## Média observada:  2.25
cat("Média da PPD:    ", round(mean(ypred_medio), 2), "\n")
## Média da PPD:     2.25

3.3 Diagnóstico da Cadeia MCMC

Após o ajuste do modelo, é fundamental avaliar a qualidade das amostras geradas pelo algoritmo MCMC. A tabela de diagnósticos MCMC inclui três medidas principais: o erro padrão de Monte Carlo (mcse), o fator de redução de escala potencial (\hat{R}) e o tamanho efetivo da amostra (n_{\text{eff}}). A seguir, é descrito cada uma delas com sua interpretação e formulação matemática.

3.4 n_eff - Número Efetivo de Amostras

O número efetivo de amostras ajusta o número total de iterações levando em conta a autocorrelação da cadeia.

A fórmula é:

\[ n_{\text{eff}} = \frac{m \times n_{ite}}{1 + 2\sum_{t=1}^\infty \rho_t} \]

Onde:

  • \(m\) é o número de cadeias,

  • \(n_{ite}\) é o número de iterações por cadeia,

  • \(\rho_t\) é a autocorrelação da cadeia no lag \(t\).

Na prática, o somatório é truncado quando \(\rho_t\) começa a oscilar ou ficar próximo de zero.

Interpretação: Quanto maior o \(n_{\text{eff}}\), mais “informação útil” temos da posteriori.

3.5 MCSE - Erro Padrão Monte Carlo

O Erro Padrão Monte Carlo (MCSE) mede a incerteza introduzida ao estimar uma estatística (como a média) com um número finito de amostras autocorrelacionadas da posteriori.

A fórmula clássica para a média é:

\[ \text{MCSE}(\bar{\theta}) \approx \sqrt{\frac{\hat{\sigma}^2}{n_{\text{eff}}}} \]

Onde:

  • \(\hat{\sigma}^2\) é a variância estimada da amostra a posteriori,

  • \(n_{\text{eff}}\) é o número efetivo de amostras independentes.

Interpretação: Quanto menor o MCSE, maior a precisão da estimativa.

3.6 Rhat - Estatística de Convergência de Gelman-Rubin

O Rhat (ou \(\hat{R}\)) compara a variabilidade entre cadeias com a variabilidade dentro das cadeias. Ele é calculado como:

\[ \hat{R} = \sqrt{ \frac{\hat{V}}{W} } \]

Onde:

  • \(W\) é a variância média dentro das cadeias,

  • \(\hat{V}\) é uma estimativa da variância marginal do parâmetro com base na combinação das cadeias:

\[ W = \frac{1}{m} \sum_{j=1}^m s_j^2, \quad B = \frac{n_{ite}}{m - 1} \sum_{j=1}^m (\bar{\theta}_j - \bar{\theta})^2 \]

\[ \hat{V} = \frac{n_{ite} - 1}{n_{ite}} W + \frac{1}{n_{ite}} B \]

Onde:

  • \(\bar{\theta}_j\) é a média da cadeia \(j\),

  • \(\bar{\theta}\) é a média geral das médias.

Regra prática: Valores de \(\hat{R} < 1,01\) indicam boa convergência.

3.7 Conclusão

Essas medidas ajudam a garantir que as estimativas obtidas via MCMC sejam confiáveis. Devemos buscar:

  • MCSE pequeno em relação à variabilidade a posteriori,
  • Rhat próximo de 1 (idealmente < 1.01),
  • n_eff alto, indicando boa mistura e independência.

Usar essas métricas junto com inspeções gráficas (traceplots, densidades) é essencial em qualquer análise Bayesiana.

# analisando o ajuste
summary(ajusteBay_Stanarm)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      y ~ x1 + x2
##  algorithm:    sampling
##  sample:       3000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 1000
##  predictors:   3
## 
## Estimates:
##               mean   sd   10%   50%   90%
## (Intercept) 2.1    0.0  2.0   2.1   2.1  
## x1          5.0    0.0  4.9   5.0   5.0  
## x2          3.1    0.0  3.0   3.1   3.1  
## sigma       1.0    0.0  1.0   1.0   1.1  
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 2.2    0.0  2.2   2.2   2.3  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0   780 
## x1            0.0  1.0  3214 
## x2            0.0  1.0  2918 
## sigma         0.0  1.0   794 
## mean_PPD      0.0  1.0  1198 
## log-posterior 0.0  1.0   883 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
#prior_summary(ajusteBay_Stanarm)

# Extraindo amostras
amostras_array = as.array(ajusteBay_Stanarm)

# Acessando as cadeias 
amostra_Stanarm_beta0_c1 = amostras_array[, 1, "(Intercept)"]  # Primeira cadeia de beta0
amostra_Stanarm_beta0_c2 = amostras_array[, 2, "(Intercept)"]  # Segunda cadeia de beta0
amostra_Stanarm_beta1_c1 = amostras_array[, 1, "x1"]  # Primeira cadeia de beta1
amostra_Stanarm_beta1_c2 = amostras_array[, 2, "x1"]  # Segunda cadeia de beta1
amostra_Stanarm_beta2_c1 = amostras_array[, 1, "x2"]  # Primeira cadeia de beta1
amostra_Stanarm_beta2_c2 = amostras_array[, 2, "x2"]  # Segunda cadeia de beta1
amostra_Stanarm_sigma2_c1 = amostras_array[, 1, "sigma"]^2  # Primeira cadeia de sigma2
amostra_Stanarm_sigma2_c2 = amostras_array[, 2, "sigma"]^2  # Segunda cadeia de sigma2


param <- "x1"  # escolha um parâmetro
param <- "x2"  # escolha um parâmetro
param <- "(Intercept)"  # escolha um parâmetro
samples <- amostras_array[, , param]  # matriz n x m (iteração x cadeia)
n <- nrow(samples)
m <- ncol(samples)
medias_cadeia <- colMeans(samples)
media_global <- mean(samples)
# B: variância entre as cadeias
B <- n / (m - 1) * sum((medias_cadeia - media_global)^2)
# W: variância dentro das cadeias
variancias_dentro <- apply(samples, 2, var)
W <- mean(variancias_dentro)
# Variância estimada
var_hat <- ((n - 1) / n) * W + (1 / n) * B
# Rhat
Rhat <- sqrt(var_hat / W)
Rhat
## [1] 1.000136
# Combina todas as cadeias em um único vetor
samples_all <- c(samples)
# Autocorrelação
acf_vals <- acf(samples_all, lag.max = 1000, plot = FALSE)$acf[-1]  # Remove lag 0
# Soma truncada das autocorrelações (até onde elas permanecem positivas)
g <- 0
for (k in seq(1, length(acf_vals), by = 2)) {
  if (k+1 > length(acf_vals)) break
  soma <- acf_vals[k] + acf_vals[k + 1]
  if (soma < 0) break
  g <- g + soma
}
# n_eff
n_eff <- m * n / (1 + 2 * g)
n_eff
## [1] 778.4701
var_total <- var(samples_all)
mcse <- sqrt(var_total / n_eff)
mcse
## [1] 0.0012065
summary(ajusteBay_Stanarm)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      y ~ x1 + x2
##  algorithm:    sampling
##  sample:       3000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 1000
##  predictors:   3
## 
## Estimates:
##               mean   sd   10%   50%   90%
## (Intercept) 2.1    0.0  2.0   2.1   2.1  
## x1          5.0    0.0  4.9   5.0   5.0  
## x2          3.1    0.0  3.0   3.1   3.1  
## sigma       1.0    0.0  1.0   1.0   1.1  
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 2.2    0.0  2.2   2.2   2.3  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0   780 
## x1            0.0  1.0  3214 
## x2            0.0  1.0  2918 
## sigma         0.0  1.0   794 
## mean_PPD      0.0  1.0  1198 
## log-posterior 0.0  1.0   883 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
round(c(mcse, Rhat, n_eff),1)
## [1]   0.0   1.0 778.5
par(mfrow=c(3,4))
plot(amostra_Stanarm_beta0_c1, xlab="iteração", ylab=expression(beta[0]), type="l",bty="n")
lines(amostra_Stanarm_beta0_c2, col=2)
#
plot(amostra_Stanarm_beta1_c1, xlab="iteração", ylab=expression(beta[1]), type="l",bty="n")
lines(amostra_Stanarm_beta1_c2, col=2)
#
plot(amostra_Stanarm_beta2_c1, xlab="iteração", ylab=expression(beta[2]), type="l",bty="n")
lines(amostra_Stanarm_beta2_c2, col=2)
#
plot(amostra_Stanarm_sigma2_c1, xlab="iteração", ylab=expression(sigma^2), type="l",bty="n")
lines(amostra_Stanarm_sigma2_c2, col=2)

hist(amostra_Stanarm_beta0_c1, ylab="densidade", xlab=expression(beta[0]), freq=F,main="")
abline(v=beta_0, lwd=2, col="blue")
abline(v=mean(amostra_Stanarm_beta0_c1), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_Stanarm_beta0_c1, 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_Stanarm_beta0_c1, 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_Stanarm_beta1_c1, ylab="densidade", xlab=expression(beta[1]), freq=F,main="")
abline(v=beta_1, lwd=2, col="blue")
abline(v=mean(amostra_Stanarm_beta1_c1), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_Stanarm_beta1_c1, 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_Stanarm_beta1_c1, 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_Stanarm_beta2_c1, ylab="densidade", xlab=expression(beta[2]), freq=F,main="")
abline(v=beta_2, lwd=2, col="blue")
abline(v=mean(amostra_Stanarm_beta2_c1), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_Stanarm_beta2_c1, 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_Stanarm_beta2_c1, 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_Stanarm_sigma2_c1, ylab="densidade", xlab=expression(sigma^2), freq=F,main="")
abline(v=sigma2, lwd=2, col="blue")
abline(v=mean(amostra_Stanarm_sigma2_c1), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_Stanarm_sigma2_c1, 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_Stanarm_sigma2_c1, 0.975), col="red", lty=3, lwd=2)

acf(amostra_Stanarm_beta0_c1, xlab="defasagem", main="")
acf(amostra_Stanarm_beta1_c1, xlab="defasagem", main="")
acf(amostra_Stanarm_beta2_c1, xlab="defasagem", main="")
acf(amostra_Stanarm_sigma2_c1, xlab="defasagem", main="")

4 Ajuste usando brms

Considere agora que os parâmetros são desconhecidos e independentes. Supondo que \(\beta_0 \sim N(0,1000)\), \(\beta_1 \sim N(0,1000)\), \(\beta_2 \sim N(0,1000)\) e \(\sigma^2 \sim GI(2,1)\), ajuste o modelo usando o pacote brms.

# Instalar o pacote brms se necessário
if (!requireNamespace("brms", quietly = TRUE)) {
  install.packages("brms")
}
library(brms)

# Definindo total de cadeias
ncadeias = 2
# Marcando o tempo inicial
tempo_inicial <- Sys.time()
# Realizando o ajuste
ajusteBay_brms = brm(
  formula = y ~ x1+x2,                    # Fórmula do modelo
  data    = dados,
  family  = gaussian(),                 # Família de distribuição para y
  prior   = c(set_prior("normal(0, 1000)", class = "b"),  #para coeficientes
            #set_prior("normal(0, 1000)", class = "b", coef = "x2"), # Prior para beta2
            set_prior("normal(0, 1000)", class = "Intercept"),  #para o intercepto
            #set_prior("exponential(1)", class = "sigma")), #Para o desvio padrao
            set_prior("inv_gamma(2, 1)", class = "sigma")), #Para o desvio padrao
  chains  = ncadeias,                          # Número de cadeias MCMC
  iter    = 2000,                         # Número de iterações
  warmup  = 500,                       # Número de iterações de aquecimento
  refresh = 0  # Suprime a saída de progresso
)
# Marcando o tempo final
tempo_final <- Sys.time()
# Calculando o tempo gasto
tempo_final - tempo_inicial
## Time difference of 3.239807 mins
# Resumo do modelo ajustado
summary(ajusteBay_brms)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ x1 + x2 
##    Data: dados (Number of observations: 1000) 
##   Draws: 2 chains, each with iter = 2000; warmup = 500; thin = 1;
##          total post-warmup draws = 3000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     2.06      0.03     2.00     2.13 1.00     3430     2300
## x1            4.96      0.03     4.90     5.03 1.00     3459     2417
## x2            3.05      0.02     3.00     3.10 1.00     4280     2508
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.04      0.02     1.00     1.09 1.00     3353     2540
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Extraindo as amostras dos parâmetros
amostras = as.data.frame(as.matrix(ajusteBay_brms))

# Salvando as amostras em matrizes separadas
amostra_brms_sigma2 = matrix(amostras$sigma^2, ncol = ncadeias, byrow = FALSE)
amostra_brms_beta0 = matrix(amostras$`b_Intercept`, ncol = ncadeias, byrow = FALSE)
amostra_brms_beta1 = matrix(amostras$`b_x1`, ncol = ncadeias, byrow = FALSE)
amostra_brms_beta2 = matrix(amostras$`b_x2`, ncol = ncadeias, byrow = FALSE)

# Renomeando as colunas para indicar a cadeia
colnames(amostra_brms_sigma2) <- paste0("Cadeia_", 1:ncadeias)
colnames(amostra_brms_beta0) <- paste0("Cadeia_", 1:ncadeias)
colnames(amostra_brms_beta1) <- paste0("Cadeia_", 1:ncadeias)
colnames(amostra_brms_beta2) <- paste0("Cadeia_", 1:ncadeias)

par(mfrow=c(3,4))
plot(amostra_brms_beta0[,1], xlab="iteração", ylab=expression(beta[0]), type="l",bty="n")
lines(amostra_brms_beta0[,2], col=2)
plot(amostra_brms_beta1[,1], xlab="iteração", ylab=expression(beta[1]), type="l",bty="n")
lines(amostra_brms_beta1[,2], col=2)
plot(amostra_brms_beta2[,1], xlab="iteração", ylab=expression(beta[2]), type="l",bty="n")
lines(amostra_brms_beta2[,2], col=2)
plot(amostra_brms_sigma2[,1], xlab="iteração", ylab=expression(sigma^2), type="l",bty="n")
lines(amostra_brms_sigma2[,2], col=2)
#
hist(amostra_brms_beta0[,1], ylab="densidade", xlab=expression(beta[0]), freq=F,main="")
abline(v=beta_0, lwd=2, col="blue")
abline(v=mean(amostra_brms_beta0[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_brms_beta0[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_brms_beta0[,1], 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_brms_beta1[,1], ylab="densidade", xlab=expression(beta[1]), freq=F,main="")
abline(v=beta_1, lwd=2, col="blue")
abline(v=mean(amostra_brms_beta1[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_brms_beta1[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_brms_beta1[,1], 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_brms_beta2[,1], ylab="densidade", xlab=expression(beta[2]), freq=F,main="")
abline(v=beta_2, lwd=2, col="blue")
abline(v=mean(amostra_brms_beta2[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_brms_beta2[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_brms_beta2[,1], 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_brms_sigma2[,1], ylab="densidade", xlab=expression(sigma^2), freq=F,main="")
abline(v=sigma2, lwd=2, col="blue")
abline(v=mean(amostra_brms_sigma2[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_brms_sigma2[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_brms_sigma2[,1], 0.975), col="red", lty=3, lwd=2)

acf(amostra_brms_beta0[,1], main="", xlab="defasagem")
acf(amostra_brms_beta1[,1], main="", xlab="defasagem")
acf(amostra_brms_beta2[,1], main="", xlab="defasagem")
acf(amostra_brms_sigma2[,1], main="", xlab="defasagem")