Thursday, January 26, 2012

Hello, Bugs world!

Here is my first Bayesian code modified from the sample code for the Chapter 3 of Bayesian Adaptive Methods for Clinical Trials. The R code depends on BRugs that interfaces with Openbugs (manual is here).
It appears that with the availability of Bugs, the knowledge of actual implementation of MCMC is less important to a user than the experiences of convergence diagnosis. I found a good slide deck here.

setwd("C:/paper/clinicaltrials/bayesian/CRM")    
library(BRugs)

#generate input data;
std.dose <- c(-1.47, -1.1, -0.69, -0.42, 0 , 0.42);
dose.index <- c(3, 4,4, 5);
y <- c(0, 0, 0, 1);

curr.n <- length(y);


mydata <- list(N1 = curr.n, s = y[1:curr.n, drop=FALSE],
n = rep(1,times=curr.n, drop=FALSE) ,
d = std.dose[dose.index, drop=FALSE], pm = 1);

bugsData(mydata, file="wei.txt")  

 modelCheck("HTexp.txt")
 modelData("wei.txt")
   modelCompile()

   modelInits("initdat.txt") # initial value of a is set to 0.5

 
   samplesSet(c("a"))
   modelUpdate(1000);

# convergence check;
   samplesAutoC('a', chain=1, mfrow=c(1,1))
   samplesAutoC('a', chain=2, mfrow=c(1,1))
   samplesAutoC('a', chain=3, mfrow=c(1,1))
   
   r1 <- samplesBgr('a', mfrow=c(1,1)) #R is around 1.01 after 500 iterations
   
   samplesDensity('a', beg=500, mfrow=c(1,1))
   
   samplesHistory('a', mfrow=c(1,1))
   samplesHistory('a', mfrow=c(1,1), thin=10) #need thin to see individual lines
   
   samplesStats('a', beg=500) #MC error is 3% of SD
   
   
   modelUpdate(10000)    
   a.p<- samplesSample("a")
   
   for (i in 1:length(std.dose)){
    cat("when dose =",std.dose[i], "expected DLT rate is", signif(mean( (exp(std.dose[i])/(exp(std.dose[i])+ exp(std.dose[i]*-1)))^a.p),3),"\n");
   }

No comments: