1.5.3 Bayesian Estimation using WinBUGS or OPENBUGS in R

As we have seen in the previous section with the two files that generates the BayesianModeling can implement the Bayesian Estimation using WinBUGS or OpenBugs programs. Also Bayesian estimation can be implemented by using R2WinBUGS or R2OpenBUGS (Sturtz, Ligges and Gelman, 2005), p ackages for Running WinBUGS and OpeBUGS from R respectively or BRugs (Thomas, et al, 2006) a collection of R functions that allow users to analyze graphical models using MCMC techniques.

Here we will need the file Beatles.txt and the skewlogitModel.txt syntax of the model generated in the BayesianModeling.

We copy all the syntax before “Inits” and save it in a file, by example as modelbr.txt. Then the file modelbr.txt would remain as

model {

for(i in 1:n) {

m[i] <- beta[1] + beta[2]*x[i]

muz[i]<- m[i]+delta*sigma*V[i]

zs[i]~dlogis(muz[i],1)I(lo[y[i]+1],up[y[i]+1])

V[i] ~ dnorm(0,1)I(0,)

}

for (j in 1:k) {

beta[j] ~ dnorm(0.0,1.0E-3)

}

delta ~ dunif(-1,1)

lambda<-delta/sqrt(1-pow(delta,2))

sigma<-1/sqrt(1-pow(delta,2))

lo[1]<- -50; lo[2]<- 0;

up[1]<- 0; up[2]<-50;

}

To implement the Bayesian estimation in R we follow the following steps to use the library R2WinBUGS.

1. In R, load the library R2WinBUGS, installed previously, with the following commando:

library(R2WinBUGS)

1. Read the data (the file beetles.txt for this Example is in F:\MILUS\beetles.txt)

datos <- read.table("F:/MILUS/beetles.txt", header=TRUE, sep="", na.strings="NA", dec=".",strip.white=TRUE)

2. Create a list that contain the data and the information of the size of the data set typing:

n=nrow(datos)

k=ncol(datos)

data<-c(as.list(datos),list(n=n,k=k))

3. Create a program that generates the initials values typing:

inits<-function(){list(beta=rep(0,k),delta=0.5)}

4. Finally with the command bugs implements the Bayesian estimation. Here will explain in brief the syntax of the command Bugs

parameters.to.save = is a vector with the names of the parameters of the model which simulations we wish to save.

model.file = is the direction where finds the file of the model

n.chains = is the number of chains to be generated.

n.iter = is the number of the total iterations of each chain.

n.burnin = is the number of iterations to be discharged as burn-in.

program = is the program that will be used to implement the Bayesian estimation

With the following command the Bayesian estimation is implemented and the simulations are saved in the object out.

out<-bugs(data,inits,parameters.to.save=c("beta","lambda"),

model.file="F:/MILUS/modelbr.txt", n.chains=1, n.iter=44000,

n.burnin=4000,program="WinBUGS")

6. If we type out in the line of commands of the R a summary of the simulation is obtained

> out

Inference for Bugs model at "G:/MILUS/modelbr.txt", fit using WinBUGS,

1 chains, each with 44000 iterations (first 4000 discarded), n.thin = 40

n.sims = 1000 iterations saved

mean sd 2.5% 25% 50% 75% 97.5%

beta[1] -62.2 3.1 -67.8 -64.4 -62.3 -60.1 -55.9

beta[2] 35.0 1.6 31.4 33.9 35.2 36.2 38.0

lambda 0.2 0.8 -1.2 -0.3 0.1 0.7 1.7

deviance 1029.3 15.7 986.7 1024.7 1035.0 1040.0 1044.0

DIC info (using the rule, pD = Dbar-Dhat)

pD = -9.3 and DIC = 1020.0

DIC is an estimate of expected predictive error (lower deviance is better).

4. Finally, to greater details about the command bugs you can obtain help typing in the line of the commands the following

?bugs

Note. You can specify Bug.directory. The directory that contains the WinBUGS executable. If the global option `R2WinBUGS.bugs.directory` is not `NULL`, it will be used as the default. Also you can specify the program to use, either `winbugs`/`WinBUGS` or `openbugs`/`OpenBUGS`, the latter makes use of function `openbugs` and requires the CRAN package BRugs.