*! Date : 19 May 2018 *! Version : 1.00 *! Author : Adrian Mander *! Email : adrian.mander@mrc-bsu.cam.ac.uk *! The fully flexible design /* START HELP FILE title[ A command to find the optimal flexible two stage single-arm binary outcome design using the discrete conditional error function. ] desc[ This command creates a two-stage fully flexible design with a binary outcome variable, often response or not, that uses a conditional error function. The design usually is interested in rejecting the null hypothesis that response is {bf:H0: p=p0}. The methods are published by S. Englert and M. Keiser (2012) Improving the flexibility and efficiency of phase 2 designs for oncology trials, Biometrics. A discrete conditional error function, D(i) is defined for the number of responders, i, in the first phase of the trial. If the second stage p-value is < D(i) then the null hypothesis is rejected. The function D() must be non-decreasing with i but this command allows no other restrictions on what D() can take. However, there is no point picking values other than the discrete set of second stage p-values. However, there are a huge number of possible first stage and second stage sample sizes but this command will search over a range of values and find the optimal one. It can take a huge amount of time if you pick large sample sizes but the search uses a branch and bound algorithm, so is efficient. Currently the optimal one is one that has the smallest expected sample size at p0. ] opt[p1() specifies the alternative probability that the study is powered for.] opt[p0() specifies the null probability that is used in the null hypothesis.] opt[alpha() specifies the significance level.] opt[beta() specifies the type II eror level.] opt[n1min() specifies the minimum first stage sample size.] opt[n1max() specifies the maximum first stage sample size.] opt[ n2min() specifies the minimum second stage sample size.] opt[n1n1max() specifies the maximum total sample size to search over.] example[ Assuming that p0=0.1, p1=0.5,we can run the command {stata cedesign, p0(0.1) p1(0.5) n1min(5) n1n2max(30)} ] author[Dr Adrian Mander] institute[MRC Biostatistics Unit, University of Cambridge] email[adrian.mander@mrc-bsu.cam.ac.uk] seealso[ {help simon2stage} (if installed) install by clicking {stata ssc install simon2stage} {help crm} (if installed) install by clicking {stata ssc install crm} ] END HELP FILE */ /* 11May18 v1.00 The command is born */ pr cedesign version 15.1 preserve syntax [, p1(real 0.3) p0(real 0.1) alpha(real 0.05) beta(real 0.1) n1min(integer 2) /// n1max(integer 22) n2min(integer 2) n1n2max(integer 50)] mata: find_design(`p0', `p1', `alpha', `beta', `n1min', `n1max', `n2min', `n1n2max') end mata: void find_design(real p0, real p1, real alpha, real beta, n1min, n1max, n2min, n1n2max) { printf("{txt}\nFinding the flexible two stage design using conditional error functions\n") printf("{dup 71:{c -}}\n") printf("{txt}p0 = {res}%f\n", p0) printf("{txt}p1 = {res}%f\n", p1) printf("{txt}alpha = {res}%f\n", alpha) printf("{txt}beta = {res}%f\n", beta) printf("{txt}\nStarting the branch and bound algorithm \n") FD=. FOC=. Fn2=. Fn1=. FRDpv2 = . /* loop over all designs conditional on n1 */ minEN = . /* set min Expected sample size as missing*/ printf("Starting loop of n1\n") for (n1=n1min;n1<=n1max; n1++) { if (n1>minEN) continue printf("{res}.%f",n1) displayflush() for (n2=n2min;n2+n1<=n1n2max;n2++) { /* some initialisations*/ pv2 = binomialtail(n2,(0..n2),p0) /* vector of P(X2>=l) there are n2+1 pvalues */ pmf1 = binomialp(n1,(0..n1),p0) pmf2 = binomialp(n2,(0..n2),p0) pmf1h1 = binomialp(n1,(0..n1),p1) pmf2h1 = binomialp(n2,(0..n2),p1) RDpv2 = (0,binomialtail(n2,(n2..0),p0)) /* D refers to all possible p-values and 0 in reverse order */ /****************************************************************************** * D is a conditional error function based on the number of responders in * stage 1 i.e. is dimension n1+1 but can take any value of the 2nd stage p-values * Start at the lowest D=(1,1,...,1) note RDpv2(D) gives you the CE function ******************************************************************************/ D=J(1,n1+1,1) while ( (D~=J(1,n1+1,n2+2)) & (D~=J(1,n1+1,n2+3))) { /*if (!checkorder(D)) { /* not sure this is a problem anymore */ printf("{err}order problem\n") D = createD(D) /* might need a check*/ } */ for(j=1;j<=n1;j++){ OC = partial_oc(j,D, pmf1, RDpv2, pmf2h1, pmf1h1, n1, n2, pv2) if (OC[1,1]>alpha) { /* failed on alpha control */ D =getnextD(D,j,n2+1) j=n1+1 continue } else if (OC[1,2]>beta) { /* failed on beta control */ D = getnextD(D,j,n2+1) j=n1+1 continue } else if (OC[1,3]>minEN) { /* failed on beating another design on expected sample size */ D = getnextD(D,j,n2+1) j=n1+1 continue } else if (j==n1) { /* Design has been found and they are stored */ FD=D FOC=OC FRDpv2 = RDpv2 Fn2 = n2 Fn1 = n1 minEN = OC[1,3] D=nearnextD(D, n2) } } } } /*end of n2 loop */ } /* end of n1 loop */ printf("\n{txt}The best conditional error function\n") FRDpv2[FD]' printf("The operating characteristics (alpha,beta, E(N))\n") FOC printf("Sample sizes n1={res}%f {txt}n2={res}%f\n", Fn1, Fn2) }/* end of find_design()*/ /*********************************************************************** * Figure out the next design after having failed to find a design at element j ***********************************************************************/ real vector getnextD(D,j,n) { len = cols(D) if(D[1,j]~=n+1) { /* check jth element isn't max n2+2 */ if (j>1) nextD = D[1..j-1], J(1,len+1-j,D[1,j]+1) else nextD = J(1,len,D[1,j]+1) } else { if (D[1,1]==n+1) return(D:+1) if ((j==2) & (D[1,2]==n+1) ) nextD = J(1,len, D[1,1]+1) if (j>2) nextD = D[1..j-2], J(1,len+2-j,D[1,j-1]+1) } return(nextD) } /****************************************************************** * Figure out the next D by adding one to the next non-max element * and making sure it is ordered *******************************************************************/ real vector nearnextD(D,n) /* n here is the sample size of 2nd stage */ { for(i=cols(D);i>1;i--) { if(D[1,i]~=n+2) { return( D[1..i-1], J(1,cols(D)-i+1,D[1,i]+1) ) } } return( J(1,cols(D),n+2) ) } /************************************************* * Check vector non-decreasing and returns 1 or 0 *************************************************/ real checkorder(real vector chk) { for(i=1;ichk[i+1]) return(0) } return(1) } /************************************************* * Check vector non-decreasing and returns 1 or 0 *************************************************/ real createD(real vector cretD) { retD = cretD for(i=1;icretD[i+1]) { retD[1,i+1..cols(retD)]=J(1,cols(retD)-i,retD[1,i]) continue } } return(retD) } /********************************************************************** * This code calculates the alpha and beta for a given D with it only * partially created * here it is the element-th value that is truncated * D is defined for first 0 to element responses **********************************************************************/ real vector partial_oc(element, subD, pmf1, Dpv2, pmf2h1, pmf1h1, n1, n2, pv2) { /* For alpha just fill in the missing values as the last value */ if (element==n2) { /* This does direct calculation on current D()*/ alpha = sum(Dpv2[subD] :* pmf1) /* Expected sample size */ EN = n1+ sum( (1:-((Dpv2[subD] :==0) :+ (Dpv2[subD] :==1 ))) :* pmf1) * n2 /* Power P_H1(P2<= D(p1)) is pv2 <= D(p1) */ power=0 for (i=1;i<=cols(subD);i++) { power = power + sum(pmf2h1:*(pv2:<=Dpv2[subD[1,i]]) ) * pmf1h1[i] } } else { /* This is the bound part of the algorithm */ newD=subD[1..element],J(1,n1+1-element,subD[1,element]) newD2=subD[1..element],J(1,n1+1-element,n2+2) /* add max possible i.e. n2+2 to min beta */ newD3=subD[1..element],J(1,n1+1-element,1) /* add in more terminations */ alpha = sum(Dpv2[newD] :* pmf1) /* Expected sample size */ EN = n1+ sum( (1:-((Dpv2[newD3] :==0) :+ (Dpv2[newD3] :==1 ))) :* pmf1) * n2 /* Power P_H1(P2<= D(p1)) is pv2 <= D(p1) */ power=0 for (i=1;i<=cols(newD2);i++) { power = power + sum(pmf2h1:*(pv2:<=Dpv2[newD2[1,i]]) ) * pmf1h1[i] } } return(alpha,1-power, EN) } /* end of partial_oc() */ end /*end of mata*/