######################################################################### # # This R program calculates the least informative normal prior # for the Bayesian CRM, given a chosen working model skeleton # # *Lee SM, Cheung YK. (2011) Calibration of prior # variance in the Bayesian CRM. Stat Med;30: 2081-89. # ######################################################################### ###Install and load required packages library(dfcrm) ###Input d<-6 #number of dose levels theta<-0.25 #target DLT rate wmfun<-"empiric" #functional form of the working model delta<-0.05 #halfwidth of the indfference interval (i.e., spacing of the skeleton) skeleton <- round(getprior(halfwidth=delta, target=theta, nu=ifelse(d %% 2 == 0,d/2,(d+1)/2), nlevel=d, model=wmfun),2) foo <- crmsens(skeleton, theta, model=wmfun, detail=TRUE) sigma2<-seq(from = 0.1, to = 2, by = 0.001) kld<-rep(0,length(sigma2)) for(j in 1:length(sigma2)){ avar=sigma2[j] muj<-rep(0,length(skeleton )) for(i in 1:length(muj)){ if(i==1){ muj[i]<-pnorm(foo\$Hset[i,2]/avar)-pnorm(-10/avar) } else { muj[i]<-pnorm(foo\$Hset[i,2]/avar)-pnorm(foo\$Hset[(i-1),2]/avar) } } kld[j]<-sum(muj*log(length(skeleton )*muj)) } #plot(sigma2,kld) ###'sli' is the least infomative prior standard deviation given 'skeleton' sli<-sigma2[which.min(kld)] print(sli) ###'priorvar' is the least informative prior variance given 'skeleton' priorvar<-sli**2 print(round(priorvar,3)) ###'muj' is the prior mass on each dose level given 'skeleton' and 'priorvar' muj<-rep(0,length(skeleton)) asd=sli for(i in 1:length(skeleton)){ if(i==1){ muj[i]<-pnorm(foo\$Hset[i,2]/asd)-pnorm(-10/asd) } else { muj[i]<-pnorm(foo\$Hset[i,2]/asd)-pnorm(foo\$Hset[(i-1),2]/asd) } } print(round(muj,2))