#########################################################################
#									
# 	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))