Software - Exemplo 1

Recurso ao R2OpenBUGS

Depois de abrir a sessão no R (ou R-Studio), deve instalar-se o software R2OpenBUGS (Lunn et al., 2009) usando a instrução

install.packages("R2OpenBUGS", dependencies=TRUE,
repos="http://cran.us.r-project.org")

Para simplificar notação vai-se sempre admitir que todos os ficheiros necessários para a execução dos programas estão guardados no diretório em que foi aberta a sessão de R. Se não, defina a sua diretoria de trabalho usando e.g.

setwd("C:/.../EBA/Exercicios")

Os passos necessários para executar o modelo no R2OpenBUGS são os seguintes:

B1. Escrever o código do modelo e guardá-lo num ficheiro com a extensão .txt. Neste caso foi guardado com o nome Example_Line_Bugs1.txt).

model
{
for(i in 1:n){
   Y[i] ~ dnorm(mu[i],tau)
   mu[i] <- beta[1] + beta[2] * x[i]
 }

for(j in 1:2){
   beta[j] ~ dnorm(0,0.0001)
 }
tau ~ dgamma(0.001,0.001)
sigma <- 1 / sqrt(tau)
}

O primeiro ciclo de for corresponde à formulação do modelo probabilístico indicado no ponto (A1). No BUGS o segundo parâmetro (\(tau\)) da distribuição é a precisão (inverso da variância). Neste ciclo está também incluída a definição do preditor linear \(mu\) que aparecem no ponto (A2). O número de observações é \(N=5\).
O segundo ciclo de for corresponde à formulação da distribuição a priori para os parâmetros de regressão. Seguidamente está formulada a distribuição a priori do último parâmetro do modelo. Todas as distribuições a priori formuladas refletem a natureza vaga da informação a priori existente. Para poder monitorizar o desvio padrão ele também é definido, por fim, em função da respetiva precisão. De notar que a ordem destas instruções é perfeitamente arbitrária.
Para se verificar o ficheiro que define o modelo a ser usado pelo OpenBUGS pode escrever-se a instrução

file.show("Example_Line_Bugs1.txt")

B2. Fazer a leitura dos dados vide ficheiro Example_Line_Data1.txt).
É aconselhável que as covariáveis que aparecem na definição do preditor linear sejam centradas na correspondente média para efeitos de aceleração do processo de convergência das cadeias. Se o ficheiro dos dados não estiver construído desse modo, esse procedimento pode ser depois feito no R.

Example <- read.table("Example_Line_Data1.txt",header=T)
dim(Example)

## [1] 5 2

head(round(Example,3))

##     x Y
## 1 1.1 1
## 2 1.2 3
## 3 1.3 3
## 4 1.4 3
## 5 1.5 5

B3. Definir os vetores da matriz dos dados a usar no modelo. Os dados têm de ser fornecidos na forma de vetor, matriz ou lista.
Criação de objetos separados para cada variável. Criação de uma lista com os dados que serão fornecidos ao OpenBUGS.

Y <- Example$Y
x <- Example$x
n <- nrow(Example)

Example.data1 <-  list("Y","x","n")

B4. Definir os parâmetros que se pretende monitorizar.

Example.params1 <- c("beta","tau","sigma")

Se se quisesse monitorizar o mu, ele também deveria aparecer na lista acima.

B5. Definir os valores iniciais dos parâmetros do modelo. Neste caso há que definir valores iniciais para beta e tau. Eles têm de ser definidos através de uma lista, devendo aparecer uma lista de valores iniciais para cada cadeia.

Example.Inits1 <- list(tau=1,beta=c(0,0))

Se se tiver mais do que uma cadeia deve-se criar uma lista do mesmo género para cada cadeia, com valores iniciais diferentes, colocando em objetos com nomes, por exemplo, Inits1, Inits2,… e dar a instrução

Example.Inits1 <- list(Inits1,Inits2,...)

sendo portanto uma lista de listas.

B6. A execução do software OpenBUGS pode ser agora feita através do R utilizando a função bugs(), depois de ter garantido que o pacote OpenBUGS se encontra carregado na sessão do R.

library(R2OpenBUGS)

## Warning: package 'R2OpenBUGS' was built under R version 3.5.3

Example_BUGS.fit1 <- bugs(data=Example.data1, inits=list(Example.Inits1), parameters.to.save=Example.params1, "Example_Line_Bugs1.txt", n.chains=1, n.iter=20000, n.burnin=10000, debug=FALSE, save.history=FALSE, DIC=TRUE)

É conveniente averiguar quais as componentes que constituem a função bugs() através da instrução ?bugs().

B7. Para obter um sumário da distribuição a posteriori marginal dos parâmetros que foram declarados no vetor guardado no objeto Example.params1, escreve-se

Example_BUGS.fit1$summary

##                mean        sd        2.5%     25%      50%      75%
## beta[1]  -7.5296287 3.7661993 -15.1500000 -9.5940 -7.54600 -5.50575
## beta[2]   8.0969876 2.8776831   2.1198500  6.5470  8.09900  9.67300
## tau       1.9527595 1.5460008   0.1778975  0.8291  1.56900  2.64600
## sigma     0.9475351 0.5613027   0.4099975  0.6147  0.79825  1.09800
## deviance 12.4606655 3.2071924   8.7419250 10.0900 11.62000 13.95000
##               97.5%
## beta[1]   0.2520825
## beta[2]  13.8802500
## tau       5.9490750
## sigma     2.3710250
## deviance 20.6900000

O que aqui se obtém são medidas sumárias das distribuições a posteriori marginais dos parâmetros do modelo, como seja a média, desvio padrão e quantis de probabilidades 0.025, 0.25, 0.5, 0.75 e 0.975. Por exemplo, ao olhar para estes elementos pode dizer-se que uma estimativa de \(\beta_2\) (declive da equação de regressão) é 8.097 e um intervalo de credibilidade de 95% é (2.12,13.88). Para obter intervalos HPD tem de se recorrer ao CODA ou BOA. Ir-se-á ver na subsecção concernente a tais pacotes como estes intervalos podem ser obtidos.

B8. Diversos outros elementos estão presentes em Example_BUGS.fit1. Usando a instrução

names(Example_BUGS.fit1)

##  [1] "n.chains"        "n.iter"          "n.burnin"       
##  [4] "n.thin"          "n.keep"          "n.sims"         
##  [7] "sims.array"      "sims.list"       "sims.matrix"    
## [10] "summary"         "mean"            "sd"             
## [13] "median"          "root.short"      "long.short"     
## [16] "dimension.short" "indexes.short"   "last.values"    
## [19] "isDIC"           "DICbyR"          "pD"             
## [22] "DIC"             "model.file"

obtendo-se a lista dos elementos desse objeto. Por exemplo, obtém-se informação sobre o DIC escrevendo

Example_BUGS.fit1$DIC

## [1] 15.95

Example_BUGS.fit1$pD

## [1] 3.491

B9. Por fim é imprescindível fazer um estudo de convergência do processo com recurso a métodos apropriados. Quer o pacote CODA, quer o pacote BOA do R podem ser usados para o efeito.

library(coda)

cadeia1 <- as.mcmc(Example_BUGS.fit1$sims.matrix[,c(1:2,4)])
head(cadeia1,9)    # 10 primeiros elementos da cadeia 

## Markov Chain Monte Carlo (MCMC) output:
## Start = 1 
## End = 10 
## Thinning interval = 1 
##       beta[1] beta[2]  sigma
##  [1,]  -9.878  10.020 0.7215
##  [2,] -10.410  10.790 0.5270
##  [3,]  -5.480   6.556 0.6396
##  [4,]  -8.665   9.179 0.6282
##  [5,]  -7.567   8.158 0.7123
##  [6,]  -4.989   6.369 0.9848
##  [7,]  -8.210   9.372 2.2260
##  [8,]  -7.308   7.769 0.7219
##  [9,]  -8.683   9.148 1.5180
## [10,]  -6.630   7.260 0.5698

# Gráfico do traço
par(mfrow=c(3,1))
traceplot(cadeia1,col="blue",sub="cadeia 1")

dev.off()

## null device 
##           1

# Testes de diagnóstico
geweke.diag(cadeia1,frac1=0.1,frac2=0.5)

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
## beta[1] beta[2]   sigma 
##  0.4016 -0.4819 -1.0657

heidel.diag(cadeia1)

##                                       
##         Stationarity start     p-value
##         test         iteration        
## beta[1] passed       1         0.988  
## beta[2] passed       1         0.989  
## sigma   passed       1         0.781  
##                                   
##         Halfwidth Mean   Halfwidth
##         test                      
## beta[1] passed    -7.530 0.0738   
## beta[2] passed     8.097 0.0564   
## sigma   passed     0.948 0.0107

raftery.diag(cadeia1)

## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
##                                                
##          Burn-in  Total Lower bound  Dependence
##          (M)      (N)   (Nmin)       factor (I)
##  beta[1] 2        3693  3746         0.986     
##  beta[2] 2        3620  3746         0.966     
##  sigma   2        3802  3746         1.010

# Intervalos HPD 
HPDinterval(cadeia1,prob=0.95)

##           lower   upper
## beta[1] -14.870  0.4144
## beta[2]   2.105 13.8500
## sigma     0.329  1.9340
## attr(,"Probability")
## [1] 0.95

B10. Uma análise completa do modelo irá envolver a seleção do melhor modelo e um estudo de adequabilidade. Todo esse estudo pode ser feito através do R. Para tal basta ter acesso aos valores simulados dos parâmetros do modelo. Esses valores encontram-se na forma de lista, array ou matriz. Esta última, por exemplo, obtém-se através da instrução

SimulatedData1 <- Example_BUGS.fit1$sims.matrix
dim(SimulatedData1)

## [1] 10000     5

head(SimulatedData1)  # mostra as seis primeiras linhas

##      beta[1] beta[2]   tau  sigma deviance
## [1,]  -9.878  10.020 1.921 0.7215   10.000
## [2,] -10.410  10.790 3.600 0.5270   18.340
## [3,]  -5.480   6.556 2.444 0.6396    9.163
## [4,]  -8.665   9.179 2.534 0.6282    9.858
## [5,]  -7.567   8.158 1.971 0.7123    8.970
## [6,]  -4.989   6.369 1.031 0.9848   11.400

De notar que só aparecem os valores simulados dos parâmetros que se declararam no objeto Example.params1. Para fazer um estudo dos resíduos, por exemplo, é importante ter monitorizado o mu.

Referências:

Lunn, D., Spiegelhalter, D., Thomas, A., Best, N. (2009). The BUGS project: Evolution, critique, and future directions, Statistics in Medicine, 28, 3049-3067. OpenBugs webpage

Plummer, M., Best, N., Cowles, K., Vines, K. (2006). CODA: Convergence Diagnosis and Output Analysis for MCMC. R News, 6, 7-11. CODA webpage

Paulino, C.D., Amaral Turkman, M.A., Murteira, B., Silva, G.L. (2018). Estatística Bayesiana, 2ª edição. Fundação Calouste Gulbenkian, Lisboa. Book webpage