Calculos para implementar o algoritimo na mão

Para a implementação do algoritmo, fez-se necessário o cálculo das distribuições condicionais completas a posteriori (DCCP), primeiramente veja os resultados obtidos para \(\tau\):

• DCCP de \(\tau\):

\[ \tau|y_1, ...,y_n,\beta_0, \beta_1 \sim Gama ( \frac{n}{2}+a,b+\frac{1}{2} \sum^n_{i=1}(y_i-\beta_0-\beta_1 x_i)^2 ) \]

Em seguida, veja o resultado obtido para \(\beta_0\), o coeficiente linear da reta, isto é, a altura em que a reta de regressão intercepta o eixo dos \(Y\)’s:

• DCCP de \(\beta_0\):

\[ \beta_0 | y_1,...,y_n , \tau,\beta_1 \sim N(\dfrac{(\tau\sum^n_{i=1}y_i - \tau\beta_1\sum^n_{i=1}x_i +\frac{m_0}{\sigma_0^2})}{ \tau n + \frac{1}{\sigma_0^2}}, (n\tau + \frac{1}{\sigma_0^2} )^{-1}) \]

Por fim, veja o resultado obtido para \(\beta_1\), é o coeficiente angular da reta, ou seja, é o a variação esperada na variável \(Y\) quando a variável explicativa é acrescida de 1 unidade:

• DCCP de \(\beta_1\):

\[ \beta_1 | y_1,...,y_n , \tau,\beta_0 \sim N(\frac{\tau\sum^n_{i=1}x_i y_i - \tau\beta_0\sum^n_{i=1}x_i + \frac{m_1}{\sigma_1^2}}{\tau \sum^n_{i=1}x_i^2 + \frac{1}{\sigma_1^2}}, ( \tau \sum^n_{i=1}x_i^2 + \frac{1}{\sigma_1^2} )^{-1}) \]

Agora que todas as distribuições condicionais completas estão calculadas o algorítimo já pode ser implementado, no R foi feito da seguinte maneira: (note que as linhas que foram comentadas executariam uma barra de carregamento, com ilustrado em seguida)

# pb <- txtProgressBar(min = 0, max = nsim, style = 3) # iniciando barra de processo
for (k in 2:nsim){
  
  #Cadeia tau
  cadeia.tau[k]   <-  rgamma(1, (n/2) + a, b + (sum((y - cadeia.b0[k-1] - (cadeia.b1[k-1]*x))^2)/2))
  
  # Cadeia B0
  c0              <-  (n*cadeia.tau[k]) + (1/sig0)
  m0              <-  (cadeia.tau[k]*sum(y) - (cadeia.tau[k]*cadeia.b1[k-1]*sum(x)) + (mu0/sig0))/c0
  cadeia.b0[k]    <-  rnorm(1, m0, 1/sqrt(c0))
  
  # Cadeia B1
  c1              <-   (sum(x^2)*cadeia.tau[k]) + (1/sig1)
  m1              <-   ((cadeia.tau[k]*sum(x*y)) - (cadeia.tau[k]*cadeia.b0[k]*sum(x)) + (mu1/sig1))/c1
  cadeia.b1[k]    <-   rnorm(1, m1, 1/sqrt(c1))
  
  # setTxtProgressBar(pb, k)
  
}# ;close(pb) #Encerrando barra de processo

Resultados da cadeia

A seguir definiremos a variável inds que indica os valores após a amostra de aquecimento (ou burn-in), a variável real que contém os valores reais utilizados para gerar a amostra para conferir se o modelo foi capaz de recuperá-los, os nomes dos parâmetros e os resultados das cadeias foram agregados em uma matriz:

# Juntando resultados:
inds    <- seq(burnin, nsim) # Definindo os indices
real    <- c(b0, b1, tau)
name    <- c(expression(beta[0]), expression(beta[1]), expression(tau))
results <- cbind(cadeia.b0, cadeia.b1, cadeia.tau) %>% as.data.frame() %>% .[inds, ] %T>% head

g1 <- hist_den(results[,1],name = name[1], p = real[1])
g2 <- hist_den(results[,2],name = name[2], p = real[2])
g3 <- hist_den(results[,3],name = name[3], p = real[3])
grid.arrange(g1,g2,g3,ncol=1)

# Cadeia
cadeia(results, name, real)

# ACF
FAC(results)

coef <- coeficientes(results, real = real) %>% as.data.frame()

print(coef)

##      Média Desv. Pad.  IC inf  IC sup Real
## b0  1.0244     0.0230  0.9792  1.0697  1.0
## b1  0.4933     0.0241  0.4464  0.5409  0.5
## tau 1.9001     0.0850  1.7371  2.0695  2.0

Comparando com o modelo linear clássico

Agora que os resultados sob o paradigma bayesiano já foram conferidos será ajustado um modelo de regressão linear simples pelo método dos mínimos quadrados através da função lm() sob o paradigma clássico para comparar com os resultados de um modelo de regressão linear simples sob o paradigma bayesiano utilizando os resultados calculados.

# Reta do modelo classico
plot(x, y)
modelo.classico <- lm(y ~ 1 + x)
a.classico      <- modelo.classico$coefficients[1]
b.classico      <- modelo.classico$coefficients[2]
abline(a        <- a.classico, b = b.classico, col = "blue")

O modelo estimado para estes dados sob o paradigma da inferência clássica foi o seguinte: \(\hat{y} = 1.0245 x + 0,4933\), o que mostra que as estimativas de \(\beta_0\) e \(\beta_1\) foram muito parecidas com as estimativas sob o paradigma da inferência bayesiana.

# Reta do modelo bayesiano
plot(x, y)
a.bayes  <-  mean(results[, 1])
b.bayes  <-  mean(results[, 2])
abline(a = a.bayes, b = b.bayes, col = "red")

A figura apresenta o gráfico de dispersão entre as variáveis da amostra simulada e as retas dos ajustes de ambos os modelos:

library(stringr)
library(ggplot2)
library(ggExtra)

# Texto da imagem
text.classico <- str_c("Modelo Classico: ","y = ",round(a.classico,4)," x + ",round(b.classico,4))
text.bayes    <- str_c("Modelo Bayesiano: ","y = ",round(a.bayes,4)," x + ",round(b.bayes,4))

# Gerando o e ambos:
cbind(y, x) %>%
  as.data.frame %>%
    ggplot(aes(y = y, x = x)) +
    geom_point() +
    geom_smooth(method = "lm", se = F, col = "red") +
    theme_classic() +
    geom_abline(slope = b.bayes,
    intercept = a.bayes,
    col = "blue") +
    labs(title = "",
    x = "Covariável",
    y = "Reposta")

Agora que os resultados no algoritmo já foram conferidos e avaliados de maneira satisfatória utilizando os dados simulados, é a vez de fazer o ajuste para dados reais.

Histograma e densidade

A figura abaixo exibe os histogramas com as densidades de três cadeias obtidas ao se iniciar o amostrador em pontos diferentes de todos os parâmetros \(\theta\) mas dessa vez sem a linha vermelha que indicava o valor do parâmetro real pois agora ele é desconhecido.

g1 <- hist_den(results[, 1], name = name[1])
g2 <- hist_den(results[, 2], name = name[2])
g3 <- hist_den(results[, 3], name = name[3])
grid.arrange(g1, g2, g3, ncol = 1)

Nota-se que ambas as cadeias convergiram uma mesma distribuição e que as últimas três cadeias apresentaram valores próximos.

Cadeias

A figura abaixo apresenta os traços das cadeias dos parâmetros amostrados. Note que há indícios de convergência.

cadeia(results,name)

Autocorrelação

A Figura abaixo apresenta os gráficos de autocorrelação dos parâmetros amostrados.

FAC(results)

É possível notar que apenas nas primeiras defasagens das cadeias das estimativas para os parâmetros \(\beta_0\) e \(\beta_1\) se apresentaram de forma autocorrelacionada e que a partir dessa defasagem o gráfico de autocorrelação se apresentou de forma desejável.

Estimativas

Como todas as características da cadeia gerada foram avaliadas de maneira satisfatória agora será possível conferir o ajuste dos parâmetros de maneira mais segura pois já foi constatada a convergência da cadeia

Comparando com o modelo linear clássico

Agora que os resultados sob o paradigma bayesiano já foram conferidos novamente será ajustado um modelo de regressão linear simples pelo método dos mínimos quadrados sob o paradigma clássico para comparar com os resultados do um modelo de regressão linear simples sob o paradigma bayesiano utilizando os resultados calculados na seção.

# Reta do modelo classico 
plot(x, y)
modelo.classico <- lm(y ~ 1 + x)
a.classico      <- modelo.classico$coefficients[1]
b.classico      <- modelo.classico$coefficients[2]
abline(a        <- a.classico, b = b.classico, col = "blue")

# Reta do modelo bayesiano
plot(x, y)
a.bayes <- mean(results[, 1])
b.bayes <- mean(results[, 2])
abline(a = a.bayes, b = b.bayes, col = "red")

A Tabela abaixo apresenta o resumo a posteriori dos parâmetros estimados da cadeia e note que esta tabela não conta com a coluna dos valores reais como no exemplo anterior e sim as estimativas sob o paradigma clássico.

coef <-
  coeficientes(results,real = classico) %>% as.data.frame()

print(coef)

##      Média Desv. Pad.  IC inf  IC sup   Real
## b0  8.2374     0.8481  6.5848  9.9239 8.2839
## b1  0.1663     0.0170  0.1326  0.1997 0.1656
## tau 0.1083     0.0214  0.0699  0.1542 0.1025

O modelo estimado sob este paradigma pode ser escrito da seguinte maneira: \(\hat{y} = 8,2839 x + 0,1656\), ou seja, os valores de \(\beta_0\) e de \(\beta_1\) novamente foram muito próximos dos parâmetros obtidos ao estimar sob o paradigma clássico.

Comparando de forma visual

A Figura ilustra o gráfico de dispersão dos dados citados acima, com a intenção de exibir quanto uma variável é afetada por outra, onde no eixo vertical representa a velocidade do carro e no eixo horizontal a distância tomada para parar.

Além do comportamento das variáveis, neste gráfico é exibido também os resultados obtidos do ajuste ao se utilizar o método de mínimos quadrados (representada pela linha em vermelho) para estimar os parâmetros e o ajuste do modelo ao se utilizar o método apresentado acima em (representada pela linha azul).

# Texto da imagem
text.classico <- str_c("Modelo Classico: ","y = ",round(a.classico,4)," x + ",round(b.classico,4))
text.bayes    <- str_c("Modelo Bayesiano: ","y = ",round(a.bayes,4)," x + ",round(b.bayes,4))

#Gerando o scatter.plot
cbind(y, x) %>%
  as.data.frame %>%
  ggplot(aes(y = y, x = x)) +
  geom_point() +
  geom_smooth(method = "lm", se = F, col = "red") +
  theme_classic() +
  geom_abline(slope = b.bayes,
              intercept = a.bayes,
              col = "blue") +
  labs(title = "Relação entre a Distância e a Velocidade com \nreta do modelo linear clássico vs bayesiano",
       x = "Distância",
       y = "Velocidade")

É possível notar que os coeficientes calculados foram muito parecidos, mesmo apresentando pequenas diferenças decimais no valor dos coeficientes ainda é possível notar que as retas estão basicamente sobrepostas, ou seja, os valores estimados em ambas as abordagens foram praticamente os mesmos.

Apesar dos valores dos ajustes terem apresentado basicamente os mesmo resultados, a maneira de se conferir a qualidade do ajuste é diferente em ambas as abordagens. Enquanto sob o paradigma clássico o ajuste do modelo pode ser checado ao avaliar os pre-supostos quanto à distribuição dos resíduos, como recomenda @GaussClarice, ao utilizar um método de MCMC faz-se necessário conferir também outros aspectos como por exemplo se houve convergência da cadeias além do comportamento das autocorrelações, vide @migon. Para outras técnicas de diagnóstico de modelos sob o paradigma clássico, veja este post sobre pacotes do R para avaliar o ajuste de modelos.