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.
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:
Post a Comment