*! Date : 1 Jul 2020 *! Version : 1.0 *! Author : Adrian Mander *! Email : mandera@cardiff.ac.uk *! The modified toxicity probability interval *! Ported from Yuan Ji's R code (7 Jan 2008) based on the mTPI method (Ji et al, 2009). /* 1Jul20 v1.00 The command is born */ /* START HELP FILE title[A command for the modified Toxicity Probability Interval design] desc[ The modified Toxicity Probability Interval design is an attempt to implement a method that is as simple as a 3+3 decision but based on statistical distributions rather than a pure rule-based approach. This method produces decision making tables for the fully trial so a statistician need not be involved in all the dose escalation decisions. The command {cmd:mtpi} performs simulations to find the operating characteristics alongside the decision matrix given a sample size and number of DLTs. ] opt[ep1() specifies the lower bound of the interval around the TTL that is safe.] opt[ep2() specifies the upper bound of the interval around the TTL that is safe.] opt[safe1() specifies the safety cutoff for declaring the lowest dose to be safe.] opt[safe2() specifies the safety cutoff for declating the dose to be safe.] opt[ttl() specifies the target toxicity level that defines the maximum tolerated dose.] opt[simn() specifes the number of simulations used] opt[sampsize() specifies the sample size of each trial in the simulations.] opt[csize() specifies the cohort size used in the trial. ] opt[startdose() specifies the starting dose. ] opt[priora() specifies the a parameter of the prior Beta distribution.] opt[priorb() specifies the a parameter of the prior Beta distribution.] opt[truep() specifies a vector of true probabilities of DLT for each dose level.] opt[decisiontable() specifies that the single cohort decision rules be displayed as a lookup table.] example[ The basic command will do a simulation {stata mtpi} The results show a table of each dose level with the true probabilities of a DLTs and the proportion of time those doses are recommended, experimented on and how many DLTs were observed on average. ] author[Prof Adrian Mander] institute[Cardiff University] email[mandera@cardiff.ac.uk] references[ Yuan Ji, Ping Liu, Yisheng Li and NB Bekele (2010) A modified toxicity probability interval method for dose-finding trials. Clin. Trials 7(6) 653-663. ] seealso[ {help crm} (if installed) {stata ssc install crm} (to install this command) ] END HELP FILE */ pr mtpi /* Allow use on earlier versions of stata that have not been fully tested */ local version = _caller() if `version' < 16.1 { di "{err}WARNING: Tested only for Stata version 16.1 and higher." di "{err}Your Stata version `version' is not officially supported." } else { version 16.1 } preserve syntax [varlist] [, ep1(real 0.05) ep2(real 0.05) safe1(real 0.8) safe2(real 0.8) ttl(real 0.3) /* */ simn(real 100) sampsize(real 30) csize(real 3) startdose(real 4) priora(real 1) priorb(real 1) /* */ truep(numlist ascending >=0 <=1) decisiontable ] /* create a string of true probabilities that Mata accepts */ if "`truep'"~= "" { local d 0 local troop "(" foreach value of numlist `truep' { local troop "`troop'`value'," local `d++' } local troop = substr("`troop'", 1, length("`troop'")-1) local troop "`troop')" } else { local troop "(.15, .25, .35, .45, .55, .65, .75, .85)" local d "8" } di di "{txt}modified Toxicity Probability Interval (mTPI) design" di "{dup 52:{c -}}" if "`decisiontable'"~="" { mata: singledose(`safe1', `safe2', `ep1', `ep2', `ttl', `d', `simn', `sampsize', `csize', `startdose', `priora' , `priorb') } mata: simulation(`safe1', `safe2', `ep1', `ep2', `ttl', `d', `simn', `sampsize', `csize', `startdose', `priora' , `priorb', `troop') restore end /**********************MATA STARTS ********************************************/ mata: /****************************************************************************** * PAVA is the pool adjacent violator algorithm to perform isotonic * transformation for the posterior means later ******************************************************************************/ real pava(x, wt) { n = rows(x) if (n <= 1) return(x) if ( missing(x)>0 | missing(wt)>0 ) { printf("Missing values in 'x' or 'wt' not allowed") exit(198) } lvlsets = (1::n) violator = 1 while(violator) { diffx = x[2::rows(x)] :- x[1::(rows(x)-1)] viol = (diffx :< 0) if (sum(viol)==0) violator=0 else { maxindex(viol, 1, i, w) /* select the first violator */ i = i[1] lvl1 = lvlsets[i] lvl2 = lvlsets[i + 1] ilvl = (lvlsets :== lvl1) + ( lvlsets :== lvl2) /*need full index..*/ maxindex(ilvl, 1, ii, ww) minlvl = min(ii) maxlvl = max(ii) x[minlvl::maxlvl] = J(maxlvl-minlvl+1,1,sum( x[minlvl::maxlvl] :* wt[minlvl::maxlvl]):/sum(wt[minlvl::maxlvl]) ) lvlsets[minlvl::maxlvl] = J(maxlvl-minlvl+1,1,minlvl) } } return(x) } // betavar computes variances of beta distributions real scalar betavar(a,b) { return( a*b/((a+b)^2*(a+b+1)) ) } /* Main simulation module */ void simulation(scalar xi1, scalar xi2, scalar eps1, scalar eps2, scalar pT, scalar D, scalar simN, scalar sampsize, scalar csize, scalar startdose, scalar a, scalar b, real rowvector p) { datan = J(simN,D, 0) datax = J(simN,D, 0) rez = J(simN,1,0) //############################## Start simulations ############### for(sim=1;sim<=simN;sim=sim+1) { x = J(D,1,0) n = J(D,1,0) pa=J(D,1,0) pb=J(D,1,0) q = J(3,1,0) d=startdose st=0 nodose=0 maxdose=1 toxdose=D+1 seldose=0 mindose=D while(st==0) { // ## st = 1 indicates the trial must be terminated maxdose = max( (maxdose, d) ) mindose = min( (mindose, d) ) // generate random toxicity response xx = sum( uniform(csize,1) :< p[d] ) x[d] = x[d]+xx n[d] = n[d]+csize // Update posterior beta distribution pa[d]=x[d]+a pb[d]=n[d]-x[d]+b if (dxi1) { tt=1 toxdose=d+1 } else tt=0 } // Compute the UPM for three intervals defined by the equivalence interval q[1] = (1-ibeta(pa[d], pb[d], eps2+pT) )/(1-eps2-pT) q[2] = (ibeta( pa[d], pb[d], eps2+pT) - ibeta( pa[d], pb[d],pT-eps1))/(eps2+eps1) q[3] = (ibeta( pa[d], pb[d], pT-eps1)/(pT-eps1))*(1-tt) // implement the dose-assignment rules based on the UPM if (d==1){ //## if the first dose is too toxic, the trial will be terminated if ((1-ibeta( 0.5+x[d], 0.5+n[d]-x[d], pT))>xi1) { /* Jeffreys prior */ st=1 nodose=1 } else { if ((q[2]>q[1])&(q[2]>q[3])) d=d /* Stay */ if ((q[3]>q[1])&(q[3]>q[2])) d=d+1 /* Escalate */ } } else { if (d==D) { //## if the last dose is much lower than the MTD, the trial will be terminated and no dose selected if (ibeta(0.5+x[d], 0.5+n[d]-x[d], pT)>xi2) { /* they used Jeffrey prior again*/ st=1 nodose=1 } else { if ((q[1]>q[2])&(q[1]>q[3])) d=d-1 /* de-escalate */ if ((q[2]>q[1])&(q[2]>q[3])) d=d /* Stay */ } } else { if ((d>1)&(dq[2])&(q[1]>q[3])) d=d-1 /* de-escalat */ if ((q[2]>q[1])&(q[2]>q[3])) d=d /* Stay */ if ((q[3]>q[1])&(q[3]>q[2])) d=d+1 /* Escalate */ } } } total=sum(n) if (total >= sampsize) st=1 } //### compute the posterior mean from the beta distribution if(nodose==0) { tdose = min( (maxdose, toxdose-1) ) pp = J(tdose, 1, -100) ppvar = J(tdose,1,0) for(i=1; i<= tdose; i++) { pp[i] = (x[i]+.005)/(n[i]+.01) ppvar[i] = betavar(x[i]+0.005, n[i]-x[i]+0.005) // here adding 0.005 is to use beta(0.005, 0.005) for estimating the MTD, which is different from the dose-finding. } pp=pava(pp, wt=1:/ppvar) // ##pp[maxdose:(toxdose-1)]<-pava(pp[maxdose:(toxdose-1)], wt=1/pp.var[maxdose:(toxdose-1)]) ## perform the isotonic transformation using PAVA with weights being posterior variances // ##for(i in maxdose:(toxdose-1)){ for(i=2;i<=tdose;i++) { pp[i] = pp[i] + i*1E-10 //## by adding an increasingly small number to tox prob at higher doses, it will break the ties and make the lower dose level the MTD if the ties were larger than pT or make the higher dose level the MTD if the ties are smaller than pT } minindex(abs(pp:-pT)[mindose::tdose],1, seldose, w) /* select from experimented doses!*/ //##seldose is the final MTD that is selected based on the order-transofromed posterior means rez[sim] = seldose+mindose-1; } for(i=1;i<=D;i++){ datan[sim,i] = n[i] datax[sim,i] = x[i] } } //##rez[is.na(rez)]<-0 aaa=J(D,1,0) //################## output results ################################ for(i=1;i<=D;i++){ aaa[i] = sum(rez:==i)/simN //### aaa is the propotion of be selected as the MTD } nodoses = 1-sum(aaa) printf("{txt}Simulation with the following settings \n") printf("{txt} Prior {res}Beta(a=%2.0f, b=%2.0f)\n", a,b) printf("{txt} Target toxicity level = {res}%f {txt}and interval {res}(eps1=%5.2f, eps2=%5.2f) \n", pT, eps1, eps2) printf("{txt} Safety cutoffs: xi1= {res}%5.2f {txt}and xi2={res}%5.2f \n", xi1, xi2) printf("{txt} Number of simulations {res}%f\n", simN) printf("{txt} Sample size {res}%f\n", sampsize) printf("{txt} Cohort size {res}%f\n", csize) printf("{txt} Starting dose {res}%f\n", startdose) printf("\n{txt}A table of the simulation results") printf("\n {c TLC}{dup 17:{c -}}{c TT}") for (i=1;in) { printf(" ") } else { // Posterior beta distribution pa=x+a pb=n-x+b // Compute the UPM for three intervals defined by the equivalence interval q[1] = (1-ibeta(pa, pb, eps2+pT) )/(1-eps2-pT) q[2] = (ibeta( pa, pb, eps2+pT) - ibeta( pa, pb,pT-eps1))/(eps2+eps1) q[3] = (ibeta( pa, pb, pT-eps1)/(pT-eps1)) if ((q[1]>q[2])&(q[1]>q[3])) printf("{res} D") /* de-escalat */ if ((q[2]>q[1])&(q[2]>q[3])) printf("{res} S") /* Stay */ if ((q[3]>q[1])&(q[3]>q[2])) printf("{res} E") /* Escalate */ /* unacceptable dose */ if ((1-ibeta( 0.5+x, 0.5+n-x, pT))>xi1) { /* Jeffreys prior */ printf("U") } else printf(" ") } } printf("\n{txt}") } printf("\n Key: E = Escalate D = De-escalate S = Stay U = Unsafe\n\n") } end /* end of Mata*/