BRugs Example

Igor — Jun 4, 2015, 12:53 PM

##!/usr/bin/R
## -*- coding: utf-8 -*-

rm(list=ls()) # clear out the memory from previous data
par(mfrow = c(1,1)) # one plot per page
options(digits=7) # number of digits to print on screen

#### SET WORKING DIRECTORY
## script.dir <- dirname(sys.frame(1)$ofile)
## setwd(script.dir)

#### LOAD BRUGS LIBRARY:
library(BRugs)

Welcome to BRugs connected to OpenBUGS version 3.2.3

#### =================================================

#### LOAD MODEL FILE INTO BRUGS AND CHECK ITS SYNTAX:
modelCheck("hemobugs.txt")

model is syntactically correct

#### LOAD THE DATA FOR BRUGS MODEL FROM A DAT FILE:
in.data  <- read.table("hemodialysis.dat", header = FALSE)
prepre  <- in.data[[3]]
postpre <- in.data[[4]]

#### SPECIFY THE DATA FOR BRUGS MODEL, AS A LIST:
datalist = list(N=length(postpre),
  prepre=prepre, postpre=postpre)

#### GET THE DATA INTO BRUGS:
modelData(bugsData(datalist))

data loaded

#### COMPILE THE MODEL AND INITIALIZE THE CHAINS:
modelCompile()

model compiled

#### LOAD OR GENERATE INITIAL VALUES
initlist = list(delta=0, taudelta=1)
modelInits(bugsData(initlist))

Initializing chain 1: 
model is initialized

modelGenInits()

model is already initialized

#### =================================================
#### RUN THE MCMC SIMULATIONS:

#### BURN IN, TO REACH THE TERMINAL DISTRIBUTION:
N.burn <-1000
modelUpdate(N.burn)

1000 updates took 0 s

#### SPECIFY, WHICH PARAMETERS TO TRACE:
parameters <- c("delta","taudelta","sigma","pH0")
samplesSet(parameters)

monitor set for variable 'delta'
monitor set for variable 'taudelta'
monitor set for variable 'sigma'
monitor set for variable 'pH0'

#### SET LENGTH OF THE MAIN MARKOV CHAIN:
chainlength = 20000


#### ACTUALLY GENERATE THE CHAIN:
modelUpdate(chainlength)

20000 updates took 0 s

#### =================================================

#### SHOW POSTERIOR STATISTICS
sample.statistics=samplesStats("*")
print(sample.statistics)

              mean       sd  MC_error  val2.5pc    median val97.5pc start
delta    -9.471000 3.041000 2.226e-02 -15.46000 -9.464000 -3.532000  1001
pH0       0.068300 0.252300 1.754e-03   0.00000  0.000000  1.000000  1001
sigma    15.810000 2.246000 1.487e-02  12.12000 15.570000 20.950000  1001
taudelta  0.004234 0.001158 7.384e-06   0.00228  0.004126  0.006809  1001
         sample
delta     20000
pH0       20000
sigma     20000
taudelta  20000

#### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
chain <- data.frame(delta = samplesSample("delta"),
                    taudelta = samplesSample("taudelta"),
                    sigma = samplesSample("sigma"),
                    pH0 = samplesSample("pH0"))


par(mfrow = c(2,2)) # plot 4 plots per pa
par(mar=c(3,3,4,1)); # set marging of the plots


#### PLOT THE MCMC CHAINS:
## samplesHistory("*", beg=N.burn, end=N.burn+300,
##                mfrow = c(2,2), ask=NULL)

for(p_ in parameters)
  {
    plot(chain[[p_]][1:300], main=p_, type="l",
         ylab=NA, xlab=NA, col="red")
  }

#### PLOT AUTOCORRELATIONS:
## samplesAutoC("*", 1,  mfrow = c(2,2), ask=NULL)

for(p_ in parameters)
  {
    acf(chain[[p_]], main=p_,lwd=4,col="red")
  }

#### PLOT THE HISTOGRAMS OF THE SAMPLED VALUES
## samplesDensity("*", 1,  mfrow = c(2,2), ask=NULL)

for(p_ in parameters)
  {
    hist(chain[[p_]], main=p_,
         ylab=NA, xlab=NA,
         nclas=50, col="red")
  }