example R2OpenBUGS

Igor — Jun 4, 2015, 12:55 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


#### LOAD R2OPENBUGS LIBRARY:
library(R2OpenBUGS)

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

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


#### LOAD INITIAL VALUES
inits <- function()
{
  list(delta=0, taudelta=1)
}


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


#### ACTUALLY GENERATE THE CHAIN:
mcmc.simulation =
  bugs(model.data,
       inits,
       model.file="hemobugs.txt",
       parameters=parameters,
       n.chains=1,
       n.iter=20000,
       n.burnin=1000,
       n.thin=1,
       codaPkg=FALSE)


#### SHOW POSTERIOR STATISTICS
print(mcmc.simulation)

Inference for Bugs model at "hemobugs.txt", 
Current: 1 chains, each with 20000 iterations (first 1000 discarded)
Cumulative: n.sims = 19000 iterations saved
          mean  sd  2.5%   25%   50%   75% 97.5%
delta     -9.5 3.0 -15.4 -11.5  -9.5  -7.5  -3.5
taudelta   0.0 0.0   0.0   0.0   0.0   0.0   0.0
sigma     15.8 2.2  12.1  14.2  15.6  17.1  20.9
pH0        0.1 0.3   0.0   0.0   0.0   0.0   1.0
deviance 233.5 2.1 231.5 232.0 232.9 234.3 239.1

DIC info (using the rule, pD = Dbar-Dhat)
pD = 2.1 and DIC = 235.6
DIC is an estimate of expected predictive error (lower deviance is better).

#### PLOT POSTERIOR STATISTICS
plot(mcmc.simulation)

#### MCMC CHAIN IN ONE DATA FRAME
chain=mcmc.simulation$sims.list

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:
for(p_ in parameters)
  {
    plot(chain[[p_]][1:300], main=p_, type="l",
         ylab=NA, xlab=NA, col="red")
  }

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

#### PLOT THE HISTOGRAMS OF THE SAMPLED VALUES
for(p_ in parameters)
  {
    hist(chain[[p_]], main=p_,
         ylab=NA, xlab=NA,
         nclas=50, col="red")
  }