############################################################################
#  Accelerated Titration Design Model for Phase I Clinical Trials          #
############################################################################
# The program fits the model from "Accelerated Titration Designs for PhaseI#
# Clinical Trials in Oncology" by  Simon R., Freidlin B., Rubinstein L.,   #
# Arbuck S., Collins G., Christian M.                                      #
#  Journal of the National Cancer Institute, 1997, V89, No. 15, p 1138-1147#
#                                                                          #
# Program was written by Boris Freidlin,  bfreidlin@emmes.com (301) 2998655#
# in collaboration with  Richard Simon,    rich@brb.nci.nih.gov            #
#                        Larry Rubinstein  rubinsteinl@ctep.nci.nih.gov    #
############################################################################
# The program is written in S-PLUS  (updated on 04/24/98)                  #
# Before running the program, load the functions                           #
#                          by typing source('...../ph1atd').               #
# To run the program,                                                      #
# type ph1atd('inputfile','doselist') at the S-PLUS prompt,                #
# where the "inputfile" is an ASCII file which contains the patient data   #
# in the following format, one record(line) per course:                    #
# Column 1:        patient number                                          #
# Column 2:        course number                                           #
# Column 3:        dose                                                    #
# Column 4:        toxicity                                                #
#      records should be ordered by patient and course number              #
# and "doselist" is  an ASCII file which contains all the doses planned    #
# to be used in the trial, listed in a single column in ascending order    #
############################################################################
# The output contains:                                                     #
# 1) max likelihood estimates of the model parameters                      #
# 2) confidence intervals for parameters with non-zero MLE                 #
#    both MLE and CI are sent to the screen and file "ph1atd.out"          #
# 3) graph of probabilities of grade2+,3+ and 4+ toxicities (for the first #
#    course) averaged over the population of patients as in Figure 5       #
#    in the publication.(the graph is stored in file "ph1atd1.wmf")        #                    #
# 4) graphs of probabilities of grade2+,3+ and 4+ toxicities (for the      #
#    first course) for the patients with beta equal to MLE-SD, MLE         #
#    MLE+SD as in Figure 6 in the publication. (the graph is stored        #
#    in file "ph1atd2.wmf")                                                #
#    Note: if the MLE of sigma beta is 0 these graphs are not produced.    #
############################################################################
pco_function(bi,vp,d,dc,tox)
{

##  the following matrix contains log(d+alpha*D)+Bi:
##  rows:     correspond to the treatment courses
##  columns:  correspond to ss randomly generated Bi for the MC
mc1_t(matrix(bi,length(bi),length(d)))+
               matrix(log(d+vp[4]*dc),length(d),length(bi))

## first level of integration: integration over epsilon from Ki to K(i+1)
mc1[,]_pnorm(ifelse(tox<2,vp[1],ifelse(tox==2,sum(vp[c(1:2)]),
           ifelse(tox==3,sum(vp[c(1:3)]),Inf))),mean=mc1,sd=vp[6])-
             pnorm(ifelse(tox<2,-Inf,ifelse(tox==2,vp[1],ifelse(tox==3,
                     sum(vp[c(1:2)]),sum(vp[c(1:3)])))),mean=mc1,sd=vp[6])
mc1
}



ph1atd_function(d0in,dosein)
{
set.seed(7)
# ss is the # of Monte-Carlo replications for the second level integration
ss_1000
d00_matrix(scan(d0in),ncol=4,byrow=T)
d00_d00[,2:4]
dosev_matrix(scan(dosein),ncol=1,byrow=T)
flik_function(vp,dvf=dddm3,vr1=vr)
{

fl_1
vr1_vr1*vp[5]

for (pn in (1:dvf[dim(dvf)[1],4]))
{

dv_dvf[(dvf[,4]==pn),(1:3)]
## one patient (all courses) contribution to the likelihood
fl_fl+
  log(if (length(dv)>4) mean(apply(pco(vr1,vp,dv[,1],dv[,2],dv[,3]),2,prod))
                   else mean(pco(vr1,vp,dv[1],dv[2],dv[3])))
                    }
         -fl
}

# following modifications of the above function use log scale
# for sb and se
# this is needed for the CI calculations
flikd_function(vp,dvf=dddm3,vr1=vr)
{
fl_1
# log scale
vr1_vr1*exp(vp[5])
vp[6]_exp(vp[6])
for (pn in (1:dvf[dim(dvf)[1],4]))
{
dv_dvf[(dvf[,4]==pn),(1:3)]
## one patient (all courses) contribution to the likelihood
fl_fl+
  log(if (length(dv)>4) mean(apply(pco(vr1,vp,dv[,1],dv[,2],dv[,3]),2,prod))
                   else mean(pco(vr1,vp,dv[1],dv[2],dv[3])))
}
-fl
}

## this part takes the raw data (d00) in the form
##  row1  course
##  row2  dose level
##  row3  toxicity
## and  converts to file ready for fitting (ddd):
##  row1 dose
##  row2 cum dose
##  toxicity
##  patient number
pc_0
ddd_matrix(0,dim(d00)[1],4)
d00[,3]_ifelse(d00[,3]==0,1,d00[,3])
ddd[,3]_d00[,3]
ddd[,1]_d00[,2]
pn_0

for (i in (1:dim(d00)[1]))
{
if(d00[i,1]==1) ddd[i,2]_0
if(d00[i,1]!=1) ddd[i,2]_ddd[i-1,2]+ddd[i-1,1]
if(d00[i,1]==1) pn_pn+1
ddd[i,4]_pn
}
vpi_matrix(c(0,0,0,0,.5,.5),1,6)
## initial guess for k1: mean of log dose levels with grade2 tox at course1
vpi[1]_mean(log(d00[(d00[,1]==1 & d00[,3]==2),2]))
## initial guesses for k2 and k3 correspond to 200% increments (log3)
ttd_sort(d00[,2])
vpi[c(2:3)]_log(3)
vpl_c(log(range(ttd)[1]),0,0,0,0,0)
vpu_log(max(range(ttd)[2],vpi[2]))*c(1,1,1,1,1,1)
print(date())
vr_matrix(0,ss,1)
vr_rnorm(ss,mean=0,sd=1)
mle_nlminb(vpi,flik,lower=vpl,upper=vpu,dvf=ddd,vr1=vr,
  control=nlminb.control(rel.tol=.00001,
         x.tol=.0001,step.min=.00001))$parameters


## confidence intervals part
mled_matrix(NA,6,3)
mled[1:3,1]_cumsum(mle[1:3])
mled[4:6,1]_mle[4:6]
dimnames(mled)_list(c(' K1   ',' K2   ',' K3   ', 'Alpha  ','Sigma Beta  ',
    'Sigma Epsilon '), c(' MLE  ', ' Lower Bound', ' Upper bound'))
mler_round(mle,3)
mle[5:6]_log(mle[5:6])

sde_matrix(-1,1,6)

del_mle[1]*diag(.00001,6,6)
# CIs are calculated only for non zero parameters (alpha, sigma beta)
del_del[,mler!=0]
npr_dim(del)[2]
# this loop  estimates of the matrix of second derivatives of the
# log likelihood until it is nonsingular
repeat {
dd2_matrix(0,npr,npr)

for (i2 in (1:npr)) {
for (j2 in (i2:npr)) {

dd2[i2,j2]_(flikd((mle+del[,i2]+del[,j2]),dvf=ddd,vr1=vr)-
               flikd((mle+del[,i2]),dvf=ddd,vr1=vr)-
                    flikd((mle+del[,j2]),dvf=ddd,vr1=vr)+
                         flikd((mle),dvf=ddd,vr1=vr))/(del[1,1]^2)
dd2[j2,i2]_dd2[i2,j2]
}}

# get out if dd2 is not singular reduce del otherwise
dd2i_apply(dd2,1,sum)
if (prod(dd2i)!=0 | max(diag(del))>.01) break
diag(del)[seq(along=dd2i)[dd2i==0]]_10*diag(del)[seq(along=dd2i)[dd2i==0]]
}
if (prod(dd2i)!=0 )
{
fish_solve(dd2)
sde[1]_sqrt(fish[1,1])
sde[2]_sqrt(matrix(1,1,2)%*%fish[1:2,1:2]%*%matrix(1,2,1))
sde[3]_sqrt(matrix(1,1,3)%*%fish[1:3,1:3]%*%matrix(1,3,1))
if (mler[4]<=.1 & mler[5]==0 ) sde[6]_sqrt(fish[4,4])
   else {if (mler[4]<=.1 ) { sde[5]_sqrt(fish[4,4])
                          sde[6]_sqrt(fish[5,5])}
         else {if (mler[5]==0 ) { sde[4]_sqrt(fish[4,4])
                          sde[6]_sqrt(fish[5,5])}
                   else { sde[4]_sqrt(fish[4,4])
                          sde[5]_sqrt(fish[5,5])
                          sde[6]_sqrt(fish[6,6]) }
                       } }



mled[1:3,2]_cumsum(mle[1:3])-1.96*sde[1:3]
mled[1:3,3]_cumsum(mle[1:3])+1.96*sde[1:3]
if (sde[4]>0) {mled[4,2]_(mle[4]-1.96*sde[4])
               mled[4,3]_(mle[4]+1.96*sde[4]) }
if (sde[5]>0) {mled[5,2]_exp(mle[5]-1.96*sde[5])
               mled[5,3]_exp(mle[5]+1.96*sde[5]) }
if (sde[6]>0) {mled[6,2]_exp(mle[6]-1.96*sde[6])
               mled[6,3]_exp(mle[6]+1.96*sde[6]) }
}
else print("WARNING: Unable to calculate confidence intervals")
print("Accelerated Titration Design Model for Phase I Clinical Trials")
print("        MLE and 95%CI  for the model parameters     ")
print(round(mled,3))
# put the output into the file
sink('ph1atd.out')
print("Accelerated Titration Design Model for Phase I Clinical Trials")
print("        MLE and 95%CI  for the model parameters     ")
print(round(mled,3))
sink()
#dose  scale
#dls_sort(d00[d00[,1]==1,2])
#nd_length(dls)
#dls1_matrix(0,1,nd)
#dls1[1:(nd-1)]_dls[2:nd]
#dls_dls[dls!=dls1]
nd_length(dosev)
dls_dosev

prbt_matrix(0,nd,3)
doselevel_c(1:nd)
dv_sqrt(sum(mled[5:6,1]^2))
prbt[1:(nd),1]_(round((1-pnorm((mled[1,1]-log(dls)-0)/dv)),4))
prbt[1:(nd),2]_(round((1-pnorm((mled[2,1]-log(dls)-0)/dv)),4))
prbt[1:(nd),3]_(round((1-pnorm((mled[3,1]-log(dls)-0)/dv)),4))

win.printer(format="placeable metafile",file='ph1atd1.wmf')

matplot(doselevel,prbt,type="l",lty=c(2,9,5),lwd=1,col=1,ylim=c(0,1),
xlab='Dose Level',ylab="Probability",lab=c(25,10,9),axes=F)
title('Probabilities of Grade 2+, Grade 3+, Grade 4+ Toxicities at the MLE')
axis(2)
axis(1,at=doselevel)
box()
legend(1,1,legend=c("Grade 2+","Grade 3+","Grade 4+"),lty=c(2,9,5))
dev.off()


# individual plots a done only if MLE(sigma beta)>0

if (mled[5,1]>0) {
win.printer(format="placeable metafile",file='ph1atd2.wmf',
                     orientation='portrait')
par(mfrow=c(3,1))
prbt[1:(nd),1]_(round((1-pnorm((mled[1,1]-log(dls)+mled[5,1])/mled[6,1])),4))
prbt[1:(nd),2]_(round((1-pnorm((mled[2,1]-log(dls)+mled[5,1])/mled[6,1])),4))
prbt[1:(nd),3]_(round((1-pnorm((mled[3,1]-log(dls)+mled[5,1])/mled[6,1])),4))
matplot(doselevel,prbt,type="l",lty=c(2,9,5),lwd=1,col=1,ylim=c(0,1),
xlab='Dose Level',ylab='Probability',lab=c(25,10,9),las=1,axes=F)
legend(1,1,legend=c("Grade2+","Grade3+","Grade4+"),lty=c(2,9,5))
title('Probabilities of the Toxicities at the MLE(beta)-SD')
axis(2,at=c(0,.5,1))
 axis(1,at=doselevel)
box()

prbt[1:(nd),1]_(round((1-pnorm((mled[1,1]-log(dls))/mled[6,1])),4))
prbt[1:(nd),2]_(round((1-pnorm((mled[2,1]-log(dls))/mled[6,1])),4))
prbt[1:(nd),3]_(round((1-pnorm((mled[3,1]-log(dls))/mled[6,1])),4))
matplot(doselevel,prbt,type="l",lty=c(2,9,5),lwd=1,col=1,ylim=c(0,1),
xlab='Dose Level',ylab='Probability',lab=c(25,10,9),las=1,axes=F)
legend(1,1,legend=c("Grade 2+","Grade 3+","Grade 4+"),lty=c(2,9,5))
title('Probabilities of the Toxicities at the MLE(beta)')
axis(2,at=c(0,.5,1))
 axis(1,at=doselevel)
box()



prbt[1:(nd),1]_(round((1-pnorm((mled[1,1]-log(dls)-mled[5,1])/mled[6,1])),4))
prbt[1:(nd),2]_(round((1-pnorm((mled[2,1]-log(dls)-mled[5,1])/mled[6,1])),4))
prbt[1:(nd),3]_(round((1-pnorm((mled[3,1]-log(dls)-mled[5,1])/mled[6,1])),4))
matplot(doselevel,prbt,type="l",lty=c(2,9,5),lwd=1,col=1,ylim=c(0,1),
xlab='Dose Level',ylab='Probability',lab=c(25,10,9),las=1,axes=F)
legend(1.5,1,legend=c("Grade 2+","Grade 3+","Grade 4+"),lty=c(2,9,5))
title('Probabilities of the Toxicities at the MLE(beta)+SD')
axis(2,at=c(0,.5,1))
axis(1,at=doselevel)
box()
dev.off()
    }

}