R Code for Introductory Adaptive Design with R

# Chapter 2

## R Function: Dunnett Test for Multiple Arms with Common Contol

## us[j] = true response in the jth group

## sigma[j] = true standard deviation for the jth group

## Zalpha =critical point for rejection

## nArms = number of active groups

## ns[j] = sample size for the jth group

## powers[j] = Probability of rejection for the jth arm

## power = prob of at least one Hi si rejected.

powerDunnett=function(n0, u0, sigma0, nArms, cAlpha, nSims)

{

powers=rep(0,nArms); tValue=rep(0,nArms); yCtl=rep(0,nArms); yTest=rep(0,nArms)

power=0

for (i in 1: nSims) {

yCtl=rnorm(n0, mean=u0, sd=sigma0)

OverSig=0

for (j in 1: nArms) {

yTest=rnorm(ns[j],mean=us[j], sd=sigmas[j])

tValue=t.test(yTest,yCtl)\$statistic

if (tValue>=cAlpha) {

powers[j]=powers[j]+1/nSims

OverSig=1}

}

power=power+OverSig/nSims

}

return (c("power=", power, "Rejection Probs = ", powers))

}

## Determine critical value Zalpha for alpha (power) =0.025 ##

ns=c(288,288,288); us=c(0,0,0); sigmas=c(2,2,2)

powerDunnett(n0=288, u0=0, sigma0=2, nArms=3, cAlpha=2.361, nSims=100000)

## Determine Power ##

ns=c(288,288,288); us=c(0.3,0.4,0.5); sigmas=c(2,2,2)

powerDunnett(n0=288, u0=0, sigma0=2, nArms=3, cAlpha=2.361, nSims=10000)

Examples of Invoking the function:

## Determine critical value Zalpha for alpha (power) =0.025 ##

ns=c(288,288,288); us=c(0,0,0); sigmas=c(2,2,2)

powerDunnett(n0=288, u0=0, sigma0=2, nArms=3, cAlpha=2.361, nSims=100000)

## Determine Power ##

ns=c(288,288,288); us=c(0.3,0.4,0.5); sigmas=c(2,2,2)

powerDunnett(n0=288, u0=0, sigma0=2, nArms=3, cAlpha=2.361, nSims=10000)

# Chapter 3

## Chapter 3-powerTwoStageClassicGSD with small or larger N

## sigma[j] = true standard deviation for the jth group

## Zalpha =critical point for rejection

## nArms = number of active groups

## ns[j] = sample size for the jth group

## EEP = Early efficacy stopping Probability

## power = prob of at least one Hi si rejected.

powerTwoStageClassicGSD=function(u0,u1,sigma0,sigma1, n0Stg1, n1Stg1, n0Stg2, n1Stg2, alpha1,beta1, alpha2, nSims)

{

power=0; EEP=0; EFP=0

for (i in 1: nSims) {

y0Stg1=rnorm(n0Stg1, mean=u0, sd=sigma0)

y1Stg1=rnorm(n1Stg1, mean=u1, sd=sigma1)

p1=t.test(y1Stg1,y0Stg1, alternative ="greater")\$p.value

if(p1<=alpha1){EEP=EEP+1/nSims}

if(p1>=beta1) {EFP=EFP+1/nSims}

if(p1>alpha1 & p1<beta1) {

y0Stg2=rnorm(n0Stg2, mean=u0, sd=sigma0)

y1Stg2=rnorm(n1Stg2,mean=u1, sd=sigma1)

y0=c(y0Stg1, y0Stg2)

y1=c(y1Stg1, y1Stg2)

p2=t.test(y1,y0, alternative ="greater")\$p.value

if(p2<=alpha2) {power=power+1/nSims}

}

}

power=power+EEP

aveTotalN=n0Stg1+n1Stg1+(1-EEP-EFP)*(n0Stg2+n1Stg2)

return (c("Average Total N=", aveTotalN, "power=", power, "EEP= ", EEP, "EFP=", EFP))

}

Invoke the function:

## Checking classic design type-I error ##

powerTwoStageClassicGSD(u0=0,u1=0,sigma0=1,sigma1=1, n0Stg1=36, n1Stg1=36, n0Stg2=36, n1Stg2=36, alpha1=0.0,beta1=1, alpha2=0.025, nSims=10000)

# Checking Type-I error ##

powerTwoStageClassicGSD(u0=0,u1=0,sigma0=1,sigma1=1, n0Stg1=36, n1Stg1=36, n0Stg2=36, n1Stg2=36, alpha1=0.0024,beta1=1, alpha2=0.024, nSims=10000)

## Power of Classic design ##

powerTwoStageClassicGSD(u0=0,u1=0.5,sigma0=1,sigma1=1, n0Stg1=36, n1Stg1=36, n0Stg2=36, n1Stg2=36, alpha1=0.0,beta1=1, alpha2=0.025, nSims=10000)

## Power ##

powerTwoStageClassicGSD(u0=0,u1=0.5,sigma0=1,sigma1=1, n0Stg1=36, n1Stg1=36, n0Stg2=36, n1Stg2=36, alpha1=0.0024,beta1=1, alpha2=0.024, nSims=10000)

powerTwoStageClassicGSD(u0=0,u1=0.5,sigma0=1,sigma1=1, n0Stg1=58, n1Stg1=58, n0Stg2=36, n1Stg2=36, alpha1=0.025,beta1=1, alpha2=0.024, nSims=10000)

## Chapter 3-Two-Stage GSD of NormalEndpoint with MINP-MSP-MPP

## u0, u1 = means for two treatment groups

## sigma0, sigma1 = standard deviations for two treatment groups

## n0Stg1, n1Stg1, n0Stg2, n1Stg2 = sample sizes for the two groups at stages 1 and 2

## alpha1,beta1, alpha2 = efficacy and futility stopping boundaries

## method = adaptive design method, either MSP. MPP or MINP

## w1squared = weight squared for MINP only at stage 1

## nSims = the number of simulation runs

## ESP1, ESP2, FSP1 = efficacy and futility stopping probabilities

## power, aveTotalN = power and expected total sample size

TwoStageGSDwithNormalEndpoint=function(u0,u1,sigma0,sigma1, n0Stg1, n1Stg1, n0Stg2, n1Stg2, alpha1, beta1, alpha2, method, w1squared=0.5, nSims=100000)

{

ESP1=0; ESP2=0; FSP1=0; w1=sqrt(w1squared);

for (i in 1:nSims) {

y0Stg1=rnorm(1, u0, sigma0/sqrt(n0Stg1))

y1Stg1=rnorm(1, u1, sigma1/sqrt(n1Stg1))

z1=(y1Stg1-y0Stg1)/sqrt(sigma1**2/n1Stg1+sigma0**2/n0Stg1)

t1=1-pnorm(z1)

if(t1<=alpha1){ESP1=ESP1+1/nSims}

if(t1>=beta1) {FSP1=FSP1+1/nSims}

if(t1>alpha1 & t1<beta1) {

y0Stg2=rnorm(1, u0, sigma0/sqrt(n0Stg2))

y1Stg2=rnorm(1, u1, sigma1/sqrt(n1Stg2))

z2=(y1Stg2-y0Stg2)/sqrt(sigma1**2/n1Stg2+sigma0**2/n0Stg2)

if (method=="MINP"){t2=1-pnorm(w1*z1+sqrt(1-w1*w1)*z2)}

if (method=="MSP" ){t2=t1+pnorm(-z2)}

if (method=="MPP" ){t2=t1*pnorm(-z2)}

if(t2<=alpha2) {ESP2=ESP2+1/nSims}

}

}

power=ESP1+ESP2

aveTotalN=n0Stg1+n1Stg1+(1-ESP1-FSP1)*(n0Stg2+n1Stg2)

return (c("Average total sample size=", aveTotalN, "power=", power, "ESP1= ", ESP1, "FSP1=", FSP1))

}

Invoke the R function:

## Example 3.1

TwoStageGSDwithNormalEndpoint(u0=0.05,u1=0.12, sigma0=0.22, sigma1=0.22, n0Stg1=120, n1Stg1=120, n0Stg2=120, n1Stg2=120, alpha1=0.01, beta1=1, alpha2=0.18321, method="MSP")

TwoStageGSDwithNormalEndpoint(u0=0.05,u1=0.05, sigma0=0.22, sigma1=0.22, n0Stg1=120, n1Stg1=120, n0Stg2=120, n1Stg2=120, alpha1=0.01, beta1=1, alpha2=0.18321, method="MSP")

## Example 3.2

TwoStageGSDwithNormalEndpoint(u0=0.05,u1=0.12, sigma0=0.22, sigma1=0.22, n0Stg1=113, n1Stg1=113, n0Stg2=113, n1Stg2=113, alpha1=0.01, beta1=1, alpha2=0.0033, method="MPP")

TwoStageGSDwithNormalEndpoint(u0=0.05,u1=0.05, sigma0=0.22, sigma1=0.22, n0Stg1=113, n1Stg1=113, n0Stg2=113, n1Stg2=113, alpha1=0.01, beta1=1, alpha2=0.0033, method="MPP")

## Example 3.3

TwoStageGSDwithNormalEndpoint(u0=0.05,u1=0.12, sigma0=0.22, sigma1=0.22, n0Stg1=110, n1Stg1=110, n0Stg2=110, n1Stg2=110, alpha1=0.01, beta1=1, alpha2=0.0188, w1squared=0.5, method="MINP")

TwoStageGSDwithNormalEndpoint(u0=0.05,u1=0.05, sigma0=0.22, sigma1=0.22, n0Stg1=110, n1Stg1=110, n0Stg2=110, n1Stg2=110, alpha1=0.01, beta1=1, alpha2=0.0188, w1squared=0.5, method="MINP")

## Example 3.4

TwoStageGSDwithNormalEndpoint(u0=0.05, u1=0.12, sigma0=0.22, sigma1=0.22, n0Stg1=124, n1Stg1=124, n0Stg2=105, n1Stg2=105, alpha1=0.0026, beta1=0.2143, alpha2=0.2143, method="MSP")

TwoStageGSDwithNormalEndpoint(u0=0.05, u1=0.05, sigma0=0.22, sigma1=0.22, n0Stg1=124, n1Stg1=124, n0Stg2=105, n1Stg2=105, alpha1=0.0026, beta1=0.2143, alpha2=0.2143, method="MSP")

TwoStageGSDwithNormalEndpoint(u0=0.05, u1=0.12, sigma0=0.22, sigma1=0.22, n0Stg1=130, n1Stg1=130, n0Stg2=105, n1Stg2=105, alpha1=0.0026, beta1=0.2143, alpha2=0.0038, method="MPP")

TwoStageGSDwithNormalEndpoint(u0=0.05, u1=0.05, sigma0=0.22, sigma1=0.22, n0Stg1=130, n1Stg1=130, n0Stg2=105, n1Stg2=105, alpha1=0.0026, beta1=0.2143, alpha2=0.0038, method="MPP")

TwoStageGSDwithNormalEndpoint(u0=0.05, u1=0.12, sigma0=0.22, sigma1=0.22, n0Stg1=116, n1Stg1=116, n0Stg2=116, n1Stg2=116, alpha1=0.0026, beta1=0.2143, alpha2=0.024, w1squared=0.5, method="MINP")

TwoStageGSDwithNormalEndpoint(u0=0.05, u1=0.05, sigma0=0.22, sigma1=0.22, n0Stg1=116, n1Stg1=116, n0Stg2=116, n1Stg2=116, alpha1=0.0026, beta1=0.2143, alpha2=0.024, w1squared=0.5, method="MINP")

## Chapter 3-Two-Stage GSD of Various Endpoints with MINP-MSP-MPP

## Two-Stage GSD for Different Endpoints with MINP, MSP, and MPP

## u0, u1 = true means of the group

## NId = noninferiority margin

## sigma0, sigma1 = standard deviation deo group 0 and 1

## nStg1, nStg2 = sample sizes per group at stages 1 and 2

## alpha1, beta1, alpha2 = stopping boundaries for the two stages

## tStd, tAcr = study and accrual durations

## endpoint = normal, binary, or survival

## method = MINP, MSP, or MPP

## tStd, tAcr = time of study and time of acrural for survuval study

## u0, u1 = means, proportions, or hazard rates for the two groups

## w1quared = the weight w1-squared at stage 1 of MINP

## ESP1, ESP2, FSP1 = efficacy and futility stopping probability at stages 1 and 2

## power, aveN = power and average sample size per group

## number of simulation runs

TwoStageGSDwithVariousEndpoints=function(u0,u1,sigma0=1, sigma1=1, tStd=24, tAcr=9, nStg1, nStg2, alpha1, beta1=0.5, alpha2, method="MINP", endpoint="normal", NId=0, w1squared=0.5, nSims=100000)

{

w1=sqrt(w1squared); w2=sqrt(1-w1squared)

if (endpoint=="binary") {sigma0=sqrt(u0*(1-u0)); sigma1=sqrt(u1*(1-u1))}

if (endpoint=="survival") {

Expu0d=exp(-u0*tStd); Expu1d=exp(-u1*tStd);

sigma0=u0*(1+Expu0d*(1-exp(u0*tAcr))/(tAcr*u0))**(-0.5);

sigma1=u1*(1+Expu1d*(1-exp(u1*tAcr))/(tAcr*u1))**(-0.5);

}

aveN=0; ESP1=0; ESP2=0; FSP1=0

for (i in 1:nSims) {

y0Stg1=rnorm(1, u0, sigma0/sqrt(nStg1))

y1Stg1=rnorm(1, u1, sigma1/sqrt(nStg1))

z1=(y1Stg1-y0Stg1+NId)/sqrt(sigma1**2/nStg1+sigma0**2/nStg1)

t1=1-pnorm(z1)

if(t1<=alpha1){ESP1=ESP1+1/nSims}

if(t1>=beta1) {FSP1=FSP1+1/nSims}

aveN=aveN+nStg1/nSims

if(t1>alpha1 & t1<beta1) {

y0Stg2=rnorm(1, u0, sigma0/sqrt(nStg2))

y1Stg2=rnorm(1, u1, sigma1/sqrt(nStg2))

z2=(y1Stg2-y0Stg2+NId)/sqrt(sigma1**2/nStg2+sigma0**2/nStg2)

if (method=="MINP"){t2=1-pnorm(w1*z1+w2*z2)}

if (method=="MSP" ){t2=t1+pnorm(-z2)}

if (method=="MPP" ){t2=t1*pnorm(-z2)}

if(t2<=alpha2) {ESP2=ESP2+1/nSims}

aveN=aveN+nStg2/nSims

}

}

power=ESP1+ESP2

return (c("Average n per group=", aveN, "power=", power, "ESP1= ", ESP1, "FSP1=", FSP1))

}

Invoke the function:

## Example 3.5

TwoStageGSDwithVariousEndpoints(u0=0.12,u1=0.14, nStg1=2580, nStg2=2580, alpha1=0.0026, beta1=0.5, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="binary")

TwoStageGSDwithVariousEndpoints(u0=0.14,u1=0.14, nStg1=2580, nStg2=2580, alpha1=0.0026, beta1=0.5, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="binary")

## Example 3.6

TwoStageGSDwithVariousEndpoints(u0=0.06601, u1=0.08664, nStg1=190, nStg2=190, alpha1=0.0026, beta1=1, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="survival")

TwoStageGSDwithVariousEndpoints(u0=0.08664, u1=0.08664, nStg1=190, nStg2=190, alpha1=0.0026, beta1=1, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="survival")

## Example 3.7

TwoStageGSDwithVariousEndpoints(u0=0.06601, u1=0.08664, nStg1=123, nStg2=123, alpha1=0.0026, beta1=1, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="survival", NId=0.005)

TwoStageGSDwithVariousEndpoints(u0=0.08664+0.005, u1=0.08664, nStg1=123, nStg2=123, alpha1=0.0026, beta1=1, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="survival", NId=0.005)

## Chapter 3-Two-Stage GSD for Different Endpoints with MINP-MSP-MPP for Superiority and NI Trial

## Two-Stage SSR for Different Endpoints with MINP, MSP, and MPP

## sigma0, sigma1 = standard deviation deo group 0 and 1

## nStg1 = sample size per group at stage 1

## alpha1, beta1, alpha2 = stopping boundaries for the two stages

## tStd, tAcr = study and accrual durations

## endpoint = normal, binary, or survival

## power = prob of at least one Hi si rejected.

## method = MINP, MSP, or MPP

## tStd, tAcr = time of study and time of acrural for survuval study

## u0, u1 = means, proportions, or hazard rates for the two groups

## w1quared = the weight w1-squared at stage 1 of MINP

## NId = noninferiority margin, NId=0 for superiority design

## ESP1, ESP2, FSP1 = efficacy and futility stopping probability at stages 1 and 2

## power, aveN = power and average sample size per group

## number of simulation runs

TwoStageGSDwithVariousEndpoints=function(u0,u1,sigma0=1, sigma1=1, tStd=24, tAcr=9, nStg1, nStg2, alpha1, beta1, alpha2, method, endpoint="normal", w1squared=0.5, NId=0, nSims=100000)

{

w1=sqrt(w1squared); w2=sqrt(1-w1squared)

if (endpoint=="binary") {sigma0=sqrt(u0*(1-u0)); sigma1=sqrt(u1*(1-u1))}

if (endpoint=="survival") {

Expu0d=exp(-u0*tStd); Expu1d=exp(-u1*tStd);

sigma0=u0*(1+Expu0d*(1-exp(u0*tAcr))/(tAcr*u0))**(-0.5);

sigma1=u1*(1+Expu1d*(1-exp(u1*tAcr))/(tAcr*u1))**(-0.5);

}

aveN=0; ESP1=0; ESP2=0; FSP1=0 ; u1=u1+NId

for (i in 1:nSims) {

y0Stg1=rnorm(1, u0, sigma0/sqrt(nStg1))

y1Stg1=rnorm(1, u1, sigma1/sqrt(nStg1))

z1=(y1Stg1-y0Stg1)/sqrt(sigma1**2/nStg1+sigma0**2/nStg1)

t1=1-pnorm(z1)

if(t1<=alpha1){ESP1=ESP1+1/nSims}

if(t1>=beta1) {FSP1=FSP1+1/nSims}

aveN=aveN+nStg1/nSims

if(t1>alpha1 & t1<beta1) {

y0Stg2=rnorm(1, u0, sigma0/sqrt(nStg2))

y1Stg2=rnorm(1, u1, sigma1/sqrt(nStg2))

z2=(y1Stg2-y0Stg2)/sqrt(sigma1**2/nStg2+sigma0**2/nStg2)

if (method=="MINP"){t2=1-pnorm(w1*z1+w2*z2)}

if (method=="MSP" ){t2=t1+pnorm(-z2)}

if (method=="MPP" ){t2=t1*pnorm(-z2)}

if(t2<=alpha2) {ESP2=ESP2+1/nSims}

aveN=aveN+nStg2/nSims

}

}

power=ESP1+ESP2

return (c("Average N per group=", aveN, "power=", power, "ESP1= ", ESP1, "FSP1=", FSP1))

}

Invoke the function:

## Example 3.5

TwoStageGSDwithVariousEndpoints(u0=0.12,u1=0.14, nStg1=2580, nStg2=2580, alpha1=0.0026, beta1=0.5, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="binary")

TwoStageGSDwithVariousEndpoints(u0=0.14,u1=0.14, nStg1=2580, nStg2=2580, alpha1=0.0026, beta1=0.5, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="binary")

## Example 3.6

TwoStageGSDwithVariousEndpoints(u0=0.06601, u1=0.08664, nStg1=190, nStg2=190, alpha1=0.0026, beta1=1, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="survival")

TwoStageGSDwithVariousEndpoints(u0=0.08664, u1=0.08664, nStg1=190, nStg2=190, alpha1=0.0026, beta1=1, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="survival")

## Example 3.7

TwoStageGSDwithVariousEndpoints(u0=0.06601, u1=0.08664, nStg1=123, nStg2=123, alpha1=0.0026, beta1=1, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="survival", NId=0.005)

TwoStageGSDwithVariousEndpoints(u0=0.08664+0.005, u1=0.08664, nStg1=123, nStg2=123, alpha1=0.0026, beta1=1, alpha2=0.024, method="MINP", w1squared=0.5, endpoint="survival", NId=0.005)

# Chapter 4

## Chapter 4-K-Stage Group-Sequential DesignwithMINP for NormalEndpoint

## ux, uy = true treatment means for groups x and y

## sigmax, sigmay = true standard deviation for groups x and y

## alphas = efficacy stopping boundaries

## betas = efficacy stopping boundaries

## Nx[i],Ny[i] = subsample size at stage i for groups x and y

## ESP = efficacy stopping Probability

## FSP = efficacy stopping Probability

## nStgs = number of stages of the trial

## w[k,i] = weights, calculated based on Info time

## power, aveN = power and the average total sample size

TwoArmKStgGSDwithMINPNormalEndpoint = function (nSims=100000, nStgs=2, ux=0, uy=1, NId=0, sigmax=1,sigmay=1) {

AveTotalN=0; FSP=rep(0, nStgs); ESP=rep(0, nStgs); sumWs=rep(0, nStgs)

w=matrix(rep(0,nStgs), nrow=nStgs, ncol=nStgs, byrow = TRUE)

for (k in 1:nStgs) {

for (i in 1:k) {sumWs[k]=sumWs[k]+(Nx[i]+Ny[i])**2}

sumWs[k]=sqrt(sumWs[k])

}

for (k in 1:nStgs){for (i in 1:k){w[k,i]=(Nx[i]+Ny[i])/sumWs[k]}}

for (iSim in 1:nSims) {

TS=0; TSc=rep(0, nStgs)

for (i in 1:nStgs) {

uxObs=rnorm(1, ux, sigmax/sqrt(Nx[i]))

uyObs=rnorm(1, uy, sigmay/sqrt(Ny[i]))

TS0=(uyObs-uxObs+NId)/sqrt(sigmax**2/Nx[i]+sigmay**2/Ny[i])

for (k in 1:nStgs) {TSc[k]=TSc[k]+w[k,i]*TS0}

TS=1-pnorm(TSc[i])

AveTotalN=AveTotalN+(Nx[i]+Ny[i])/nSims

if (TS<=alphas[i]){ESP[i]=ESP[i]+1/nSims; break}

if (TS>betas[i]) { FSP[i]=FSP[i]+1/nSims; break}

} # End of i

} # End of iSim

power=0; for (i in 1:nStgs) {power=power+ESP[i]}

return (c("power=", power, "Average total N =", AveTotalN, "ESP=", ESP, "FSP=", FSP))

}

Invoke the function:

# Example 4.4: Three-Stage Design

Nx = c(92,92,92); Ny = c(92,92,92); alphas = c(0.0025, 0.00575, 0.022); betas = c(0.5, 0.5, 0.022)

TwoArmKStgGSDwithMINPNormalEndpoint(nSims=100000, nStgs=3, ux=0, uy=0, NId=0, sigmax=0.22, sigmay=0.22)

TwoArmKStgGSDwithMINPNormalEndpoint(nSims=100000, nStgs=3, ux=0.05, uy=0.12, NId=0, sigmax=0.22, sigmay=0.22)

TwoArmKStgGSDwithMINPNormalEndpoint(nSims=100000, nStgs=3, ux=0.05, uy=0.10, NId=0, sigmax=0.22, sigmay=0.22)

## Chapter 4-K-Stage GSD for Different Endpoints with MINP-MSP-MPP for Superiority and NI Trial

## ux, uy = means, proportions, or hazard rates for the two groups

## sigmax, sigmay = true standard deviation for groups x and y

## alphas = efficacy stopping boundaries

## betas = efficacy stopping boundaries

## Nx[i],Ny[i] = subsample size at stage i for groups x and y

## ESP = efficacy stopping Probability

## FSP = efficacy stopping Probability

## nStgs = number of stages of the trial

## w[k,i] = weights, calculated based on Info time

## power, aveN = power and the average total sample size

## endpoint = normal, binary, or survival

## method = MINP, MSP, or MPP

## tStd, tAcr = time of study and time of acrural for survuval study

KStgGSDwithVariousEndpointsForMspMppMinp=function (ux=0.2, uy=0.3, sigmax=1, sigmay=1, tStd=24, tAcr=9, nStgs=3, method="MINP", endpoint="normal", NId=0, nSims=10000)

{

w=matrix(rep(0,nStgs), nrow=nStgs, ncol=nStgs, byrow = TRUE)

sumWs=rep(0, nStgs)

for (k in 1:nStgs) {

for (i in 1:k) {sumWs[k]=sumWs[k]+(Nx[i]+Ny[i])**2}

sumWs[k]=sqrt(sumWs[k])

}

## Determine the equivalent std dev.

if (endpoint=="binary") {

sigmax=sqrt(ux*(1-ux));

sigmay=sqrt(uy*(1-uy))+0.000000001} # to avoid numerical overflow

if (endpoint=="survival") {

Expuxd=exp(-ux*tStd); Expuyd=exp(-uy*tStd);

sigmax=ux*(1+Expuxd*(1-exp(ux*tAcr))/(tAcr*ux))**(-0.5);

sigmay=uy*(1+Expuyd*(1-exp(uy*tAcr))/(tAcr*uy))**(-0.5);

}

for (k in 1:nStgs){for (i in 1:k){w[k,i]=(Nx[i]+Ny[i])/sumWs[k]}}

AveTotalN=0; FSP=rep(0, nStgs); ESP=rep(0, nStgs); uy=uy+NId

for (iSim in 1:nSims) {

TS=0; TSc=rep(0, nStgs); if (method=="MPP") {TSc=rep(1, nStgs)}

for (i in 1:nStgs) {

uxObs=rnorm(1, ux, sigmax/sqrt(Nx[i]))

uyObs=rnorm(1, uy, sigmay/sqrt(Ny[i]))

TS0=(uyObs-uxObs+NId)/sqrt(sigmax**2/Nx[i]+sigmay**2/Ny[i])

for (k in 1:nStgs) {

if (method=="MINP") {TSc[k]=TSc[k]+w[k,i]*TS0}

if (method=="MSP") {TSc[k]=TSc[k]+pnorm(-TS0)}

if (method=="MPP") {TSc[k]=TSc[k]*pnorm(-TS0)}

}

if (method=="MINP") {TSc[i]=pnorm(-TSc[i])} ## Convert to p-scale

AveTotalN=AveTotalN+(Nx[i]+Ny[i])/nSims

if (TSc[i]<=alphas[i]){ESP[i]=ESP[i]+1/nSims; break}

if (TSc[i]>betas[i]) { FSP[i]=FSP[i]+1/nSims; break}

} # End of i

} # End of iSim

power=0; for (i in 1:nStgs) {power=power+ESP[i]}

return (c("power=", power, "Average total N =", AveTotalN, "ESP=", ESP, "FSP=", FSP))

}

Invoke the function:

# Example 4.4: Three-Stage Design Normal endpoint

Nx = c(92,92,92); Ny = c(92,92,92); alphas = c(0.0025, 0.00575, 0.022); betas = c(0.5, 0.5, 0.022)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0, uy=0, sigmax=0.22, sigmay=0.22, nStgs=3, method="MINP", endpoint="normal", NId=0, nSims=10000)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0.05, uy=0.12, sigmax=0.22, sigmay=0.22, nStgs=3, method="MINP", endpoint="normal", NId=0, nSims=10000)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0.05, uy=0.1, sigmax=0.22, sigmay=0.22, nStgs=3, method="MINP", endpoint="normal", NId=0, nSims=10000)

Nx = c(92,92,92); Ny = c(92,92,92); alphas = c(0.0025, 0.1025, 0.49257); betas = c(0.49257, 0.49257, 0.49257)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0, uy=0, sigmax=0.22, sigmay=0.22, nStgs=3, method="MSP", endpoint="normal", NId=0, nSims=10000)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0.05, uy=0.12, sigmax=0.22, sigmay=0.22, nStgs=3, method="MSP", endpoint="normal", NId=0, nSims=10000)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0.05, uy=0.10, sigmax=0.22, sigmay=0.22, nStgs=3, method="MSP", endpoint="normal", NId=0, nSims=10000)

# Example 4.5: Three-Stage Design binary endpoint

Nx = c(1860,1860,1860); Ny = c(1860,1860,1860); alphas = c(0.005, 0.00653, 0.0198); betas = c(0.45, 0.45, 0.0198)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0, uy=0, nStgs=3, method="MINP", endpoint="binary", NId=0, nSims=10000)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0.12, uy=0.14, nStgs=3, method="MINP", endpoint="binary", NId=0, nSims=10000)

Nx = c(1860,1860,1860); Ny = c(1860,1860,1860); alphas = c(0.0025, 0.00083452, 0.00071363); betas = c(1,1,0.00071363)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0, uy=0, nStgs=3, method="MPP", endpoint="binary", NId=0, nSims=10000)

KStgGSDwithVariousEndpointsForMspMppMinp(ux=0.12, uy=0.14, nStgs=3, method="MPP", endpoint="binary", NId=0, nSims=10000)

# Chapter 5

## Chapter 5-Two-Stage MaxInfo Design for Normal Endpoint

## Sample Size Reestimation with the Mixed Method

## SSR Design with Small or Larger sample Size

## sigma[j] = true standard deviation for the jth group

## alpha1, alpha2, beta1 = stopping boundaries on p-sacle

## u0, u1 = true means for the two groups

## delta0 = estimated mean difference u1-u0

## nStg1, nStg2 = sample sizes for stage 1 and 2

## Nmin, Nmax = minimum and maximum sample sizes allowed for SSR

## ESP1, ESP2, FSP = efficacy and futility stopping Probabilities

## aveN = average sample size per group

TwoStageMaxInfoSSRDesign=function(u0,u1,sigma0,sigma1, delta0, nStg1, Nmin=72, Nmax=200, alpha=0.025, cPower=0.9, nSims=10000)

{

ESP1=0; ESP2=0; FSP1=0; aveN=0; nFinal=Nmin

for (i in 1: nSims) {

y0Stg1=rnorm(nStg1, u0, sigma0)

y1Stg1=rnorm(nStg1, u1, sigma1)

p1=t.test(y1Stg1,y0Stg1, alternative ="greater")\$p.value

aveN=aveN+nStg1/nSims

## New n2 by mixed method

yStg1=c(y0Stg1, y1Stg1)

sigmaBlind=sd(yStg1)

if (sigmaBlind>0) {

nFinal=2*((sigmaBlind/delta0)**2-1/4)*(qnorm(1-alpha)+qnorm(cPower))**2-nStg1

## Force Nmin<=nFinal<=Nmax

nFinal=round(min(max(nFinal,Nmin),Nmax))

}

nStg2=nFinal-nStg1; aveN=aveN+nStg2/nSims

y0Stg2=rnorm(nStg2, u0, sigma0)

y1Stg2=rnorm(nStg2, u1, sigma1)

y0=c(y0Stg1, y0Stg2); y1=c(y1Stg1, y1Stg2)

p2=t.test(y1,y0, alternative ="greater")\$p.value

if(p2<=alpha) {ESP2=ESP2+1/nSims}

}

power=ESP1+ESP2

return (c(sigmaBlind, nFinal, "Power=", power, "Average N=", aveN, "ESP1= ", ESP1, "FSP1=", FSP1))

}

Invoke the function:

## Checking classic design type-I error ##

TwoStageMaxInfoSSRDesign(u0=0,u1=0,sigma0=3,sigma1=3, delta0=0.3, nStg1=1050, Nmin=2100, Nmax=4000, alpha=0.025, cPower=0.9)

## Power of MaxInfo SSR deisgn

TwoStageMaxInfoSSRDesign(u0=0,u1=.25,sigma0=3,sigma1=3, delta0=0.3, nStg1=1050, Nmin=2100, Nmax=4000, alpha=0.025, cPower=0.9)

TwoStageMaxInfoSSRDesign(u0=0,u1=.30,sigma0=3,sigma1=3, delta0=0.3, nStg1=1050, Nmin=2100, Nmax=4000, alpha=0.025, cPower=0.9)

TwoStageMaxInfoSSRDesign(u0=0,u1=.30,sigma0=4,sigma1=4, delta0=0.3, nStg1=1050, Nmin=2100, Nmax=4000, alpha=0.025, cPower=0.9)

TwoStageMaxInfoSSRDesign(u0=0,u1=.35,sigma0=3,sigma1=3, delta0=0.3, nStg1=1050, Nmin=2100, Nmax=4000, alpha=0.025, cPower=0.9)

## Chapter 5-Two-Stage Unblind SSR for Normal Endpoints with MINP-MSP-MPP

## u0 and u1 = parameter means for the two groups

## sigma = common standard deviation

## alpha1, alpha2, beta1 = stopping boundaries on p-scale

## nStg1, nStg2 = sample size for stages 1 and 2

## nMax = maximum sample size per gorup

## w1= stage-1 weight in MINP only

## ESP = efficacy stopping Probability

## FSP = futility stopping Probability

## power = prob of rejecting Ha.

## NId = noninferiority margin

## cPower = the targeted conditional power

## n2 = new sample size per group for stage 2

## aveN = average sample size per group

## method = MINP, MSP, or MPP

powerTwoStageUnblindSSRnormalEndpoint=function(u0,u1,sigma, nStg1, nStg2, alpha1,beta1=0.205, alpha2,NId=0, cPower, method, nMax=100, w1=0.7071, nSims=100000)

{

power=0; aveN=0; ESP=0; FSP=0; w2=sqrt(1-w1**2); n2=nStg2

for (i in 1: nSims) {

y0Stg1=rnorm(1, mean=u0, sd=sigma/sqrt(nStg1))

y1Stg1=rnorm(1, mean=u1, sd=sigma/sqrt(nStg1))

z1=(y1Stg1-y0Stg1)/sigma*sqrt(nStg1/2)

t1=1-pnorm(z1)

aveN=aveN+nStg1/nSims

if(t1<=alpha1){ESP=ESP+1/nSims}

if(t1>=beta1) {FSP=FSP+1/nSims}

if(t1>alpha1 & t1<beta1) {

# n2 by conditional power

if (method=="MSP") {BFun = qnorm(1-max(0.00001,alpha2-t1))}

if (method=="MPP") {BFun = qnorm(1- min(0.99999,alpha2/t1))}

if (method=="MINP"){BFun = (qnorm(1-alpha2)- w1*qnorm(1-t1))/w2}

eSize=(y1Stg1-y0Stg1+NId)/sigma

if (eSize >0) {n2=min(nMax-nStg1, max(2*((BFun-qnorm(1-cPower))/eSize)**2, nStg2))}

y0Stg2=rnorm(1, mean=u0, sd=sigma/sqrt(n2))

y1Stg2=rnorm(1, mean=u1, sd=sigma/sqrt(n2))

z2=(y1Stg2-y0Stg2)/sigma*sqrt(n2/2)

if (method=="MSP" ){t2=t1+pnorm(-z2)}

if (method=="MPP" ){t2=t1*pnorm(-z2)}

if (method=="MINP"){t2=1-pnorm(w1*z1+w2*z2)}

if(t2<=alpha2) {power=power+1/nSims}

aveN=aveN+n2/nSims

}

} # End of iSim

power=power+ESP

return (c("Average N=", aveN, "power=", power, "ESP= ", ESP, "FSP=", FSP))

}

Invoke the function:

# Example 5.1

# Checking Type-I error with MINP ##

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0,sigma=1, nStg1=80, nStg2=80, alpha1=0.0026, beta1=0.205, alpha2=0.024, cPower=0.9, method="MINP", nMax=320, w1=0.7071)

## Power with MINP under Ha

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.33,sigma=1, nStg1=80, nStg2=80, alpha1=0.0026, beta1=0.205, alpha2=0.024, cPower=0.9, method="MINP", nMax=320, w1=0.7071)

## Power with MINP Under Hs

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.28,sigma=1, nStg1=80, nStg2=80, alpha1=0.0026, beta1=0.205, alpha2=0.024, cPower=0.9, method="MINP", nMax=320, w1=0.7071)

# Checking Type-I error with MSP ##

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0,sigma=1, nStg1=80, nStg2=80, alpha1=0.0026, beta1=0.205, alpha2=0.21463, cPower=0.9, method="MSP", nMax=320)

## Power with MSP under Ha

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.33,sigma=1, nStg1=80, nStg2=80, alpha1=0.0026, beta1=0.205, alpha2=0.21463, cPower=0.9, method="MSP", nMax=320)

## Power with MSP under Hs

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.28,sigma=1, nStg1=80, nStg2=80, alpha1=0.0026, beta1=0.205, alpha2=0.21463, cPower=0.9, method="MSP", nMax=320)

## Checking Type-I error with MPP##

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0,sigma=1, nStg1=80, nStg2=80, alpha1=0.0025, beta1=0.205, alpha2=0.0038, cPower=0.9, method="MPP",nMax=320)

## Power with MPP##

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.33,sigma=1,nStg1=80, nStg2=80, alpha1=0.0025, beta1=0.205, alpha2=0.0038, cPower=0.9, method="MPP",nMax=320)

## Power with MPP##

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.28,sigma=1,nStg1=80, nStg2=80, alpha1=0.0025, beta1=0.205, alpha2=0.0038, cPower=0.9, method="MPP",nMax=320)

# Table 5.4

# Checking Type-I error with MINP ##

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0,sigma=3, nStg1=800, nStg2=800, alpha1=0, beta1=0.5, alpha2=0.025, cPower=0.9, method="MINP", nMax=4000, w1=0.7071)

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.25,sigma=3, nStg1=800, nStg2=800, alpha1=0, beta1=0.5, alpha2=0.025, cPower=0.9, method="MINP", nMax=3000, w1=0.7071)

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.30,sigma=3, nStg1=800, nStg2=800, alpha1=0, beta1=0.5, alpha2=0.025, cPower=0.9, method="MINP", nMax=3000, w1=0.7071)

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.30,sigma=4, nStg1=800, nStg2=800, alpha1=0, beta1=0.5, alpha2=0.025, cPower=0.9, method="MINP", nMax=3000, w1=0.7071)

powerTwoStageUnblindSSRnormalEndpoint(u0=0,u1=0.35,sigma=3, nStg1=800, nStg2=800, alpha1=0, beta1=0.5, alpha2=0.025, cPower=0.9, method="MINP", nMax=3000, w1=0.7071)

## Chapter 5-Two-Stage Unblind SSR for Different Endpoints with MINP-MSP-MPP

## Two-Stage SSR for Different Endpoints with MINP, MSP, and MPP

## sigma0, sigma1 = standard deviation deo group 0 and 1

## nStg1 = sample size per group at stage 1

## alpha1, beta1, alpha2 = stopping boundaries for the two stages

## tStd, tAcr = study and accrual durations

## endpoint = normal, binary, or survival

## power = prob of at least one Hi si rejected.

## method = MINP, MSP, or MPP

## tStd, tAcr = time of study and time of acrural for survuval study

## u0, u1 = means, proportions, or hazard rates for the two groups

## w1quared = the weight w1-squared at stage 1 of MINP

## Nmin, Nmax = minimum and maximum sample size per group allowed in SSR

## cPower = the targeted conditional power in SSR

## ESP1, ESP2, FSP1 = efficacy and futility stopping probability at stages 1 and 2

## power, aveN = power and average sample size per group

## number of simulation runs

TwoStageSSRwithVariousEndpoints=function(u0,u1,sigma0=1, sigma1=1, tStd=24, tAcr=9, nStg1, alpha1, beta1=0.5, alpha2, method="MINP", endpoint="normal", NId=0, w1squared=0.5, Nmin, Nmax, cPower=0.9, nSims=100000)

{

w1=sqrt(w1squared); w2=sqrt(1-w1squared); n2=Nmin-nStg1

if (endpoint=="binary") {sigma0=sqrt(u0*(1-u0)); sigma1=sqrt(u1*(1-u1))}

if (endpoint=="survival") {

Expu0d=exp(-u0*tStd); Expu1d=exp(-u1*tStd);

sigma0=u0*(1+Expu0d*(1-exp(u0*tAcr))/(tAcr*u0))**(-0.5)

sigma1=u1*(1+Expu1d*(1-exp(u1*tAcr))/(tAcr*u1))**(-0.5)

}

sigma=sqrt(sigma0**2+sigma1**2)

aveN=0; ESP1=0; ESP2=0; FSP1=0

for (i in 1:nSims) {

y0Stg1=rnorm(1, u0, sigma0/sqrt(nStg1))

y1Stg1=rnorm(1, u1, sigma1/sqrt(nStg1))

z1=(y1Stg1-y0Stg1+NId)/sqrt(sigma1**2/nStg1+sigma0**2/nStg1)

t1=1-pnorm(z1)

if(t1<=alpha1){ESP1=ESP1+1/nSims}

if(t1>=beta1) {FSP1=FSP1+1/nSims}

aveN=aveN+nStg1/nSims

if(t1>alpha1 & t1<beta1) {

# new n2 by conditional power

if (method=="MSP") {BFun=qnorm(1-max(0.000001,alpha2-t1))}

if (method=="MPP") {BFun=qnorm(1- min(0.999999,alpha2/t1))}

if (method=="MINP"){BFun=(qnorm(1-alpha2)- w1*qnorm(1-t1))/w2}

eSize=(y1Stg1-y0Stg1+NId)/sigma

n2=round(min(Nmax-nStg1, max(2*((BFun-qnorm(1-cPower))/eSize)^2, Nmin-nStg1)))

y0Stg2=rnorm(1, u0, sigma0/sqrt(n2))

y1Stg2=rnorm(1, u1, sigma1/sqrt(n2))

z2=(y1Stg2-y0Stg2+NId)/sqrt(sigma1**2/n2+sigma0**2/n2)

if (method=="MINP"){t2=1-pnorm(w1*z1+w2*z2)}

if (method=="MSP" ){t2=t1+pnorm(-z2)}

if (method=="MPP" ){t2=t1*pnorm(-z2)}

if(t2<=alpha2) {ESP2=ESP2+1/nSims}

aveN=aveN+n2/nSims

}

}

power=ESP1+ESP2

return (c("Average n per group=", aveN, "power=", power, "ESP1= ", ESP1, "FSP1=", FSP1))

}

Invoke the function:

## Example 5.2

## Classic design

TwoStageSSRwithVariousEndpoints(u0=0.22,u1=0.22, nStg1=145, alpha1=0, beta=1, alpha2=0.025, method="MINP", endpoint="binary", w1squared=0.5, Nmin=290, Nmax=290, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.10,u1=0.22, nStg1=145, alpha1=0, beta=1, alpha2=0.025, method="MINP", endpoint="binary", w1squared=0.5, Nmin=290, Nmax=290, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.11,u1=0.22, nStg1=145, alpha1=0, beta=1, alpha2=0.025, method="MINP", endpoint="binary", w1squared=0.5, Nmin=290, Nmax=290, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.14,u1=0.22, nStg1=145, alpha1=0, beta=1, alpha2=0.025, method="MINP", endpoint="binary", w1squared=0.5, Nmin=290, Nmax=290, nSims=100000)

## GSD with OBF boundary

TwoStageSSRwithVariousEndpoints(u0=0.22,u1=0.22, nStg1=150, alpha1=0.0026, beta=0.5, alpha2=0.024, method="MINP", endpoint="binary", w1squared=0.5, Nmin=300, Nmax=300, cPower=0.9, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.10,u1=0.22, nStg1=150, alpha1=0.0026, beta=0.5, alpha2=0.024, method="MINP", endpoint="binary", w1squared=0.5, Nmin=300, Nmax=300, cPower=0.9, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.11,u1=0.22, nStg1=150, alpha1=0.0026, beta=0.5, alpha2=0.024, method="MINP", endpoint="binary", w1squared=0.5, Nmin=300, Nmax=300, cPower=0.9, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.14,u1=0.22, nStg1=150, alpha1=0.0026, beta=0.5, alpha2=0.024, method="MINP", endpoint="binary", w1squared=0.5, Nmin=300, Nmax=300, cPower=0.9, nSims=100000)

## SSR Design with OBF boundary

TwoStageSSRwithVariousEndpoints(u0=0.22,u1=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, method="MINP", endpoint="binary", w1squared=0.5, Nmin=180, Nmax=410, cPower=0.9, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.10,u1=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, method="MINP", endpoint="binary", w1squared=0.5, Nmin=180, Nmax=410, cPower=0.9, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.11,u1=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, method="MINP", endpoint="binary", w1squared=0.5, Nmin=180, Nmax=410, cPower=0.9, nSims=100000)

TwoStageSSRwithVariousEndpoints(u0=0.14,u1=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, method="MINP", endpoint="binary", w1squared=0.5, Nmin=180, Nmax=410, cPower=0.9, nSims=100000)

## Chapter 5-Two-Stage Mixed SSR Design with Normal Endpoint

## Sample Size Reestimation with the Mixed Method

## SSR Design with Small or Larger sample Size

## sigma[j] = true standard deviation for the jth group

## alpha1, alpha2, beta1 = stopping boundaries on p-sacle

## u0, u1 = true means for the two groups

## delta0 = estimated mean difference u1-u0

## nStg1, nStg2 = sample sizes for stage 1 and 2

## Nmin, Nmax = minimum and maximum sample sizes allowed for SSR

## ESP1, ESP2, FSP = efficacy and futility stopping Probabilities

## aveN = average sample size per group

TwoStageMixedSSRDesignWithNormalEndpoint=function(u0,u1,sigma0,sigma1, delta0, nStg1, Nmin=72, Nmax=200, alpha1,beta1, alpha2, cPower=0.9, nSims=10000)

{

ESP1=0; ESP2=0; FSP1=0; aveN=0; nFinal=Nmin

for (i in 1: nSims) {

y0Stg1=rnorm(nStg1, u0, sigma0)

y1Stg1=rnorm(nStg1, u1, sigma1)

p1=t.test(y1Stg1,y0Stg1, alternative ="greater")\$p.value

aveN=aveN+nStg1/nSims

if(p1<=alpha1){ESP1=ESP1+1/nSims}

if(p1>=beta1) {FSP1=FSP1+1/nSims}

if(p1>alpha1 & p1<beta1) {

## New n2 by mixed method

yStg1=c(y0Stg1, y1Stg1)

sigmaBlind=sd(yStg1)

if (sigmaBlind>0) {

eSize=delta0/sigmaBlind

BFun=sqrt(2)*qnorm(1-alpha2)-qnorm(1-p1);

nFinal=nStg1+2*((BFun-qnorm(1-cPower))/eSize)**2

## Force Nmin<=nFinal<=Nmax

nFinal=round(min(max(nFinal,Nmin),Nmax))

}

nStg2=nFinal-nStg1; aveN=aveN+nStg2/nSims

y0Stg2=rnorm(nStg2, u0, sigma0)

y1Stg2=rnorm(nStg2, u1, sigma1)

y0=c(y0Stg1, y0Stg2); y1=c(y1Stg1, y1Stg2)

p2=t.test(y1,y0, alternative ="greater")\$p.value

if(p2<=alpha2) {ESP2=ESP2+1/nSims}

}

}

power=ESP1+ESP2

return (c("Power=", power, "Average N=", aveN, "ESP1= ", ESP1, "FSP1=", FSP1))

}

Invoke the function:

## Table 5.4

## Checking Mixed SRR design type-I error ##

TwoStageMixedSSRDesignWithNormalEndpoint(u0=0,u1=0,sigma0=3,sigma1=3, delta0=0.3, nStg1=900, Nmin=1800, Nmax=4000, alpha1=0,beta1=0.5, alpha2=0.025, cPower=0.9, nSims=100000)

## Power of Classic design ##

TwoStageMixedSSRDesignWithNormalEndpoint(u0=0,u1=0.25,sigma0=3,sigma1=3, delta0=0.3, nStg1=900, Nmin=1800, Nmax=4000, alpha1=0,beta1=0.5, alpha2=0.025, cPower=0.9, nSims=10000)

TwoStageMixedSSRDesignWithNormalEndpoint(u0=0,u1=0.3,sigma0=3,sigma1=3, delta0=0.3, nStg1=900, Nmin=1800, Nmax=4000, alpha1=0,beta1=0.5, alpha2=0.025, cPower=0.9, nSims=10000)

TwoStageMixedSSRDesignWithNormalEndpoint(u0=0,u1=0.3,sigma0=4,sigma1=4, delta0=0.3, nStg1=900, Nmin=1800, Nmax=4000, alpha1=0,beta1=0.5, alpha2=0.025, cPower=0.9, nSims=10000)

TwoStageMixedSSRDesignWithNormalEndpoint(u0=0,u1=0.35,sigma0=3,sigma1=3, delta0=0.3, nStg1=900, Nmin=1800, Nmax=4000, alpha1=0,beta1=0.5, alpha2=0.025, cPower=0.9, nSims=10000)

## Chapter 5-Two-Stage Mixed and unblind SSR Design with Binary Endpoint

## Two-Stage SSR for Binary Endpoint with mixed and unblinded MINP

## nStg1 = sample size per group at stage 1

## alpha1, beta1, alpha2 = stopping boundaries for the two stages

## power = prob of at least one Hi si rejected.

## px, py = proportions for the two groups

## Nmin, Nmax = minimum and maximum sample size per group allowed in SSR

## cPower = the targeted conditional power in SSR

## ESP1, ESP2, FSP1 = efficacy and futility stopping probability at stages 1 and 2

## power, aveN = power and average sample size per group

## number of simulation runs

TwoStageMixedSSRwithBinaryEndpoint=function(px, py, nStg1, alpha1, beta1=0.5, alpha2, px0=0.11, py0=0.22, NId=0, SSR="mixed", Nmin, Nmax, cPower=0.9, nSims=100000)

{

w1=sqrt(nStg1/Nmin); w2=sqrt(1-w1*w1); n2=Nmin-nStg1

aveN=0; ESP1=0; ESP2=0; FSP1=0

for (i in 1:nSims) {

pxStg1=rbinom(1, nStg1, px)/nStg1

pyStg1=rbinom(1, nStg1, py)/nStg1

sigma=sqrt(pxStg1*(1-pxStg1)+pyStg1*(1-pyStg1))

sigmaBlind=sqrt(2*(pxStg1+pyStg1)/2*(1-(pxStg1+pyStg1)/2))

z1=(pyStg1-pxStg1+NId)/sqrt(sigma**2/nStg1)

t1=1-pnorm(z1)

if(t1<=alpha1){ESP1=ESP1+1/nSims}

if(t1>=beta1) {FSP1=FSP1+1/nSims}

aveN=aveN+nStg1/nSims

if(t1>alpha1 & t1<beta1) {

# new n2 by conditional power

BFun=(qnorm(1-alpha2)- w1*qnorm(1-t1))/w2

if (SSR=="unblind"){eSize=(pyStg1-pxStg1+NId)/sigma}

if (SSR=="mixed"){eSize=(py0-px0+NId)/sigmaBlind}

n2=round(min(Nmax-nStg1, max(2*((BFun-qnorm(1-cPower))/eSize)**2, Nmin-nStg1)))

pxStg2=rbinom(1, n2, px)/n2

pyStg2=rbinom(1, n2, py)/n2

sigma=sqrt(pxStg2*(1-pxStg2)+pyStg2*(1-pyStg2))

z2=(pyStg2-pxStg2+NId)/sqrt(sigma**2/n2)

t2=1-pnorm(w1*z1+w2*z2)

if(t2<=alpha2) {ESP2=ESP2+1/nSims}

aveN=aveN+n2/nSims

}

}

power=ESP1+ESP2

return (c("Average n per group=", aveN, "power=", power, "ESP1= ", ESP1, "FSP1=", FSP1))

}

Invoke function:

## Example 5.2

## Classic design

TwoStageMixedSSRwithBinaryEndpoint(px=0.22,py=0.22, nStg1=145, alpha1=0, beta=1, alpha2=0.025, SSR="unblind", Nmin=290, Nmax=290, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.10,py=0.22, nStg1=145, alpha1=0, beta=1, alpha2=0.025, SSR="unblind", Nmin=290, Nmax=290, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.11,py=0.22, nStg1=145, alpha1=0, beta=1, alpha2=0.025, SSR="unblind", Nmin=290, Nmax=290, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.14,py=0.22, nStg1=145, alpha1=0, beta=1, alpha2=0.025, SSR="unblind", Nmin=290, Nmax=290, nSims=100000)

## GSD with OBF boundary

TwoStageMixedSSRwithBinaryEndpoint(px=0.22,py=0.22, nStg1=150, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=300, Nmax=300, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.10,py=0.22, nStg1=150, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=300, Nmax=300, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.11,py=0.22, nStg1=150, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=300, Nmax=300, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.14,py=0.22, nStg1=150, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=300, Nmax=300, cPower=0.9, nSims=100000)

## Unblind SSR Design with OBF boundary

TwoStageMixedSSRwithBinaryEndpoint(px=0.22,py=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=200, Nmax=410, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.10,py=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=200, Nmax=410, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.11,py=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=200, Nmax=410, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.14,py=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=200, Nmax=410, cPower=0.9, nSims=100000)

## Mixed SSR Design with OBF boundary

TwoStageMixedSSRwithBinaryEndpoint(px=0.22,py=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="mixed", Nmin=200, Nmax=410, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.10,py=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="mixed", Nmin=200, Nmax=410, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.11,py=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="mixed", Nmin=200, Nmax=410, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.14,py=0.22, nStg1=100, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="mixed", Nmin=200, Nmax=410, cPower=0.9, nSims=100000)

## GSD with OBF boundary

TwoStageMixedSSRwithBinaryEndpoint(px=0.22,py=0.22, nStg1=200, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=400, Nmax=400, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.10,py=0.22, nStg1=200, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=400, Nmax=400, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.11,py=0.22, nStg1=200, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=400, Nmax=400, cPower=0.9, nSims=100000)

TwoStageMixedSSRwithBinaryEndpoint(px=0.14,py=0.22, nStg1=200, alpha1=0.0026, beta=0.5, alpha2=0.024, SSR="unblind", Nmin=400, Nmax=400, cPower=0.9, nSims=100000)

# Chapter 6

ADwithIncompletePairs=function(alpha=0.025, alpha1=0.00153, alpha2=0.02454, beta1=1, sigma=0.433, NewSigma=0.489, N1=100, N2=100, Nmin=200, Nmax=250, u=0.1, cPower=0.9, nSims=100000)

{

w1=sqrt(N1/(N1+N2)); w2=sqrt(1-w1**2); FSP1=0; ESP1=0; ESP2=0; aveN=0;

for (iSim in 1:nSims) {

uObs=rnorm(1, u, sigma/sqrt(N1))

z1=uObs*sqrt(N1)/sigma

p1=1-pnorm(z1)

aveN=aveN+N1/nSims

if (p1>beta1){FSP1=FSP1+1/nSims}

if (p1<=alpha1){ESP1=ESP1+1/nSims}

if (p1>alpha1 & p1<=beta1){

nFinal=2*((NewSigma/u)**2-1/4)*(qnorm(1-alpha)+qnorm(cPower))**2

nFinal=round(min(max(nFinal,Nmin),Nmax))

NewN2=nFinal-N1

z2=rnorm(1,u*sqrt(NewN2)/NewSigma,1)

TS2=w1*z1+w2*z2; p2=1-pnorm(TS2);

if (p2<=alpha2) {ESP2=ESP2+1/nSims}

aveN=aveN+NewN2/nSims

}

} ## End of iSim

Power = ESP1+ESP2

return (c("power=", Power, "Average N=", aveN, "ESP1=", ESP1, "FSP1=", FSP1))

}

Invoke the function

## Classical under Ha Power

ADwithIncompletePairs(alpha=0.025, alpha1=.025, alpha2=0, beta1=0.5, sigma=0.433, NewSigma=0.433, N1=200, N2=0, Nmax = 200, u=0.1, cPower=0.9,)

ADwithIncompletePairs(alpha=0.025, alpha1=.025, alpha2=0, beta1=0.5, sigma=0.489, NewSigma=0.489, N1=200, N2=0, Nmax = 200, u=0.1, cPower=0.9,)

## Two-stage Adaptive Design with Blind SSR

ADwithIncompletePairs(alpha=0.025, alpha1=0.0026, alpha2=0.024, beta1=0.5, sigma=0.433, NewSigma=0.489, N1=100, N2=100, Nmin=200, Nmax = 250, u=0.1, cPower=0.9)

## Chapter 6-GSDwithCoprimaryEndpoints

## Two-group, two-endpoint with two stage design

## muij= standardized means of j endpoint in group i

## rho= the true correlation between the two endpoints

## c1 and c2 = stopping boundaries

## N= maximum sample size per group

## tau = Info time for the interim analysis

library(mvtnorm)

GSDwithCoprimaryEndpoints=function(mu11,mu12, mu21, mu22,rho,tau,c1,c2,N,nSims=10000)

{

ESP1=0; ESP2=0

for (i in 1:nSims)

{

n=round(N*tau)

varcov=matrix(c(1,rho,rho,1),2,2)

trtStg1=rmvnorm(n,mean=c(mu11,mu12), sigma=varcov)

ctStg1 =rmvnorm(n,mean=c(mu21,mu22),sigma=varcov)

t11=t.test(trtStg1[,1],ctStg1[,1])\$statistic

t12=t.test(trtStg1[,2],ctStg1[,2])\$statistic

trtStg2=rmvnorm(N-n, mean=c(mu11,mu12), sigma=varcov)

ctStg2 =rmvnorm(N-n, mean=c(mu21,mu22), sigma=varcov)

trt1=c(trtStg1[,1],trtStg2[,1]); trt2=c(trtStg1[,2],trtStg2[,2])

ct1 =c(ctStg1[,1] ,ctStg2[,1]); ct2 =c(ctStg1[,2] ,ctStg2[,2])

t21=t.test(trt1,ct1)\$statistic; t22=t.test(trt2,ct2)\$statistic

## if ((t11>c1 & t12>c1) | (t21>c2 & t22>c2)) {power=power+1/nSims}

if (t11>c1 & t12>c1){ESP1=ESP1+1/nSims}

else {if (t21>c2 & t22>c2){ESP2=ESP2+1/nSims}}

power=ESP1+ESP2; aveN=n+(1-ESP1)*(N-n)

}

return (c("Power=", power, "Average N=", aveN))

}

Invoke the function:

GSDwithCoprimaryEndpoints(mu11=0.2,mu12=0.25, mu21=0.005, mu22=0.015, rho=0.25, tau=0.5, c1=2.80, c2=1.98, N=584, nSims=10000)

GSDwithCoprimaryEndpoints(mu11=0.5,mu12=0.5, mu21=0, mu22=0, rho=0, tau=0.5, c1=2.80, c2=1.98, N=84, nSims=10000)

GSDwithCoprimaryEndpoints(mu11=0.5,mu12=0.5, mu21=0, mu22=0, rho=1, tau=0.5, c1=2.80, c2=1.98, N=84, nSims=10000)

## Adaptive Noninferiority Design with Paired Binary Data */

## Ho: p10-p01-delNI <= 0 */

## p10 and p01 are the % of disconcordant pairs */

## Sample size: nPairs = nPairs1 + nPairs2 from stage 1 and 2 */

## avenPairs = the expected sample size (pairs)

## alpha = one-sided significance level

## ESP1, ESP2 = Rejection probability at stages 1 and 2

## FSP1 = acceptance probability at Stage 1

## alpha1, alpha2=, beta1 = Stopping boundaries on p-scale

## nPairsMax = max number of pairs, TargetcPow=targeted conditional power

## w1 = the first stage weight, w2=sqrt(1-w1**2)

## delNI = the noninferiority margin

McNemarADwithSSR=function(alpha1=0.0026, alpha2=0.024, beta1=1, p10=0.125, p01=0.125, delNI=0.1, nPairs1=154, nPairs2=154, nPairsMax=600, TargetcPow=0.90, w1=0.707, nSims=10000)

{

w2=sqrt(1-w1**2); ESP1=0; ESP2=0; FSP1=0; avenPairs=0; nPairs20=nPairs2

for (iSim in 1:nSims) {

avenPairs=avenPairs+nPairs1/nSims

n10Stg1=rbinom(1, nPairs1,p10); n01Stg1=rbinom(1,nPairs1,p01)

p10obsStg1=n10Stg1/nPairs1; p01obsStg1=n01Stg1/nPairs1

epsStg1=p10obsStg1-p01obsStg1-delNI

b=(2+p01obsStg1-p10obsStg1)*delNI-p01obsStg1-p10obsStg1

c=-p01obsStg1*delNI*(1-delNI)

p01Wave=(-b+sqrt(b*b-8*c))/4

sigma2Stg1=2*p01Wave+delNI-delNI**2

z1=0; if (sigma2Stg1 != 0) {z1=epsStg1*sqrt(nPairs1/sigma2Stg1)}

pValue1=1-pnorm(z1); t1=pValue1

if (t1<=alpha1) {ESP1=ESP1+1/nSims}

if (t1>beta1) {FSP1=FSP1+1/nSims}

if (alpha1<t1 & t1<=beta1) {

## Sample size reestimation based on conditional power

Bval=(qnorm(1-alpha2)-w1*z1)/w2

nPairs2=nPairs20

if (epsStg1>0) {

nPairs2=sigma2Stg1/epsStg1**2*(Bval-qnorm(1-TargetcPow))**2

nPairs2=min(nPairsMax-nPairs1,nPairs2)

nPairs2=round(max(nPairs2, nPairs20))

} ## End of SSR

n10Stg2=rbinom(1,nPairs2,p10); n01Stg2=rbinom(1,nPairs2,p01)

p10obsStg2=n10Stg2/nPairs2; p01obsStg2=n01Stg2/nPairs2

epsStg2=p10obsStg2-p01obsStg2-delNI

b=(2+p01obsStg2-p10obsStg2)*delNI-p01obsStg2-p10obsStg2

c=-p01obsStg2*delNI*(1-delNI);

p01Wave=(-b+sqrt(b*b-8*c))/4;

sigma2Stg2=2*p01Wave+delNI-delNI**2;

z2=0; if (sigma2Stg2 !=0) {z2=epsStg2*sqrt(nPairs2/sigma2Stg2)}

pValue2=1-pnorm(z2); t2=1-pnorm(w1*z1+w2*z2)

if (t2<=alpha2) ESP2=ESP2+1/nSims

avenPairs=avenPairs+nPairs2/nSims;

} # End of Stage 2

} # End of isim

Power=ESP1+ESP2

return (c("Power=", Power, "Average N=", avenPairs, "ESP1=", ESP1, "ESP2=", ESP2, "FSP1=", FSP1))

}

Invoke the function:

## Classical 1-Stage Design : Type-I error rate p10-p01-delNI=0

McNemarADwithSSR(alpha1=0.025, alpha2=0, beta1=0, p10=0.1, p01=0.175, delNI=-0.075, nPairs1=322, nPairs2=0);

## Power for design with PF (rho=2), beta1=0.5 and SSR (Nmax>N0)

McNemarADwithSSR(alpha1=0.00625, alpha2=0.02173, beta1=0.5, p10=0.1, p01=0.1, delNI=-0.075, nPairs1=161, nPairs2=161, nPairsMax=500);

## Chapter 6-Two Stage Farrington-Manning NI Trial with SSR

## u0 and u1 = parameter means for the two groups

## sigma = common standard deviation

## alpha1, alpha2, beta1 = stopping boundaries on p-scale

## nStg1, nStg2 = sample size for stages 1 and 2

## nMax = maximum sample size per gorup

## w1= stage-1 weight in MINP only

## ESP = efficacy stopping Probability

## FSP = futility stopping Probability

## power = prob of rejecting Ha.

## NId = NI margin, NId=0.1 means retain 90% effect size

## cPower = the targeted conditional power

## n2 = new sample size per group for stage 2

## aveN = average sample size per group

## method = MINP, MSP, or MPP

powerTwoStageFarringtonManningNIwithSSRnormalEndpoint=function(u0,u1,sigma, nStg1, nStg2, alpha1,beta1=0.205, alpha2,NId=0.1, cPower, method, nMax=100, w1=0.7071, nSims=100000)

{

power=0; aveN=0; ESP=0; FSP=0; w2=sqrt(1-w1**2); n2=nStg2

for (i in 1: nSims) {

y0Stg1=rnorm(1, mean=u0, sd=sigma/sqrt(nStg1))

y1Stg1=rnorm(1, mean=u1, sd=sigma/sqrt(nStg1))

z1=(y1Stg1-(1-NId)*y0Stg1)/sigma*sqrt(nStg1/(2-NId))

t1=1-pnorm(z1)

aveN=aveN+nStg1/nSims

if(t1<=alpha1){ESP=ESP+1/nSims}

if(t1>=beta1) {FSP=FSP+1/nSims}

if(t1>alpha1 & t1<beta1) {

# n2 by conditional power

BFun = (qnorm(1-alpha2)- w1*qnorm(1-t1))/w2

eSize=(y1Stg1-y0Stg1+NId)/sigma

if (eSize >0) {n2=min(nMax-nStg1, max(2*((BFun-qnorm(1-cPower))/eSize)**2, nStg2))}

y0Stg2=rnorm(1, mean=u0, sd=sigma/sqrt(n2))

y1Stg2=rnorm(1, mean=u1, sd=sigma/sqrt(n2))

z2=(y1Stg2-(1-NId)*y0Stg2)/sigma*sqrt(n2/(2-NId))

t2=1-pnorm(w1*z1+w2*z2)

if(t2<=alpha2) {power=power+1/nSims}

aveN=aveN+n2/nSims

}

} # End of iSim

power=power+ESP

return (c("Average N=", aveN, "power=", power, "ESP= ", ESP, "FSP=", FSP))

}

Invoke the function:

# Example 6.2: Adaptive Design with Farrington-Manning NI Margin

powerTwoStageFarringtonManningNIwithSSRnormalEndpoint(u0=0.20,u1=0.20,sigma=0.08, nStg1=150, nStg2=150, alpha1=0.0026, beta1=0.205, alpha2=0.024, NId=0.1, cPower=0.9, method="MINP", nMax=320, w1=0.7071)

## Chapter 6-Two-Stage Equivalaence GSD of NormalEndpoint MINP

## du = mean difference between the two treatment groups

## sigma = common standard deviations

## nStg1, nStg2 = sample sizes per group at stages 1 and 2

## alpha1,beta1, alpha2 = efficacy and futility stopping boundaries

## w1squared = weight squared for MINP only at stage 1

## nSims = the number of simulation runs

## ESP1, ESP2, FSP1 = efficacy and futility stopping probabilities

## power, aveN = power and expected sample size per group

TwoStageEquivGSDwithNormalEndpoint=function(du,sigma, Eqd, nStg1, nStg2, alpha1, beta1, alpha2, w1squared=0.5, nSims=100000)

{

ESP1=0; ESP2=0; FSP1=0; w1=sqrt(w1squared); w2=sqrt(1-w1*w1)

for (i in 1:nSims) {

dyStg1=rnorm(1, du, sigma/sqrt(nStg1/2))

z1=(Eqd-abs(dyStg1))/sigma*sqrt(nStg1/2)

t1=1-pnorm(z1)

if(t1<=alpha1){ESP1=ESP1+1/nSims}

if(t1>=beta1) {FSP1=FSP1+1/nSims}

if(t1>alpha1 & t1<beta1) {

dyStg2=rnorm(1, du, sigma/sqrt(nStg2/2))

z2=(Eqd-abs(dyStg2))/sigma*sqrt(nStg2/2)

t2=1-pnorm(w1*z1+w2*z2)

if(t2<=alpha2) {ESP2=ESP2+1/nSims}

}

}

power=ESP1+ESP2

aveN=nStg1+(1-ESP1-FSP1)*nStg2

return (c("Average N per group=", aveN, "power=", power, "ESP1= ", ESP1, "FSP1=", FSP1))

}

Invoke the function:

## Example 6.1: Type I erro rcontrol with OBrien-Fleming Boundary is convervative since the pdf of the test statistics has heavy tails.

TwoStageEquivGSDwithNormalEndpoint(du=0.2, sigma=0.72, Eqd=0.2, nStg1=200, nStg2=200, alpha1=0.0052, beta1=1, alpha2=0.048, w1squared=0.5, nSims=1000000)

## Example 6.1: Power GSD for Equivalance Trial

TwoStageEquivGSDwithNormalEndpoint(du=0, sigma=0.72, Eqd=0.2, nStg1=200, nStg2=200, alpha1=0.0052, beta1=1, alpha2=0.048, w1squared=0.5, nSims=100000)

## Example 6.1: Type I erro rcontrol

TwoStageEquivGSDwithNormalEndpoint(du=0.2, sigma=0.72, Eqd=0.2, nStg1=200, nStg2=200, alpha1=0.0085, beta1=1, alpha2=0.048, w1squared=0.5, nSims=1000000)

## Example 6.1: Power GSD for Equivalance Trial

TwoStageEquivGSDwithNormalEndpoint(du=0, sigma=0.72, Eqd=0.2, nStg1=200, nStg2=200, alpha1=0.0085, beta1=1, alpha2=0.048, w1squared=0.5, nSims=100000)

# Chapter 7

Chapter 7-Classic K+1 Drop-Arm Design

## The Classical ACTIVE K+1 Arm dropping-Loser Design

## Large sample size assumption

## H0: all means are equal. Ha: at least one mean > mean0.

## For Dunnett's test, nStage2=0.

WinnerDesign=function(nSims=1000,NumOfArms=5, mu0=0, sigma=1, Z_alpha=2.407, nStg1=100, nStg2=100){

xObs=rep(0,NumOfArms); power=0

for (iSim in 1:nSims) {

for (i in 1:NumOfArms){xObs[i]=rnorm(1,mu[i],sigma/sqrt(nStg1))}

MaxRsp =xObs; SelectedArm=1

for (i in 1:NumOfArms) {

if (xObs[i]> MaxRsp) {SelectedArm=i; MaxRsp=xObs[i]}

}

x2=rnorm(1,mu[SelectedArm],sigma/max(1,sqrt(nStg2)))

FinalxAve=(MaxRsp*nStg1+x2*nStg2)/(nStg1+nStg2)

x0Ave=rnorm(1,mu0,sigma/sqrt(nStg1+nStg2))

TestZ=(FinalxAve-x0Ave)*sqrt((nStg1+nStg2)/2)/sigma

if (TestZ >= Z_alpha){power=power+1/nSims}

}

TotalN=(NumOfArms+1)*nStg1+2*nStg2

return (c("Power=", power, "Total N=", TotalN))

}

Invoke the function:

## Determine Critical Value Z_alpha for 4+1 Arms Winner Design

mu=c(0,0,0,0)

WinnerDesign(nSims=1000000,NumOfArms=4, mu0=0, sigma=1, Z_alpha=2.352, nStg1=60, nStg2=60)

## Determine Power for 4+1 Arms Winner Design

mu=c(0.12, 0.13, 0.14, 0.15)

WinnerDesign(nSims=100000,NumOfArms=4, mu0=0.06, sigma=0.18, Z_alpha=2.352, nStg1=60, nStg2=60)

## Determine Critical Value Z_alpha for 4+1 Arms Winner Design

mu=c(0,0,0,0)

WinnerDesign(nSims=1000000,NumOfArms=4, mu0=0, sigma=1, Z_alpha=2.3, nStg1=43, nStg2=86)

## Determine Power for 4+1 Arms Winner Design

mu=c(0.12, 0.13, 0.14, 0.15)

WinnerDesign(nSims=100000,NumOfArms=4, mu0=0.06, sigma=0.18, Z_alpha=2.3, nStg1=43, nStg2=86)

## Determine Critical Value Z_alpha for Dunnett Test

mu=c(0,0,0,0)

WinnerDesign(nSims=1000000,NumOfArms=4, mu0=0, sigma=1, Z_alpha=2.441, nStg1=102, nStg2=0);

## Determine Power for Dunnett Test

mu=c(0.12, 0.13, 0.14, 0.15);

WinnerDesign(nSims=100000,NumOfArms=4, mu0=0.06, sigma=0.18, Z_alpha=2.441, nStg1=102, nStg2=0)

# Chapter 8

##SelProb[i] = probability of selecting arm i as the best arm

FourPlus1AddArmDesign=function(nSims=10000, N1=100, N2=100, c_alpha = 2.267,cR = 0.55652, mu0=0, sigma=1)

{

Power=0; SelProb=rep(0,4); xObs=rep(0,4)

for (iSim in 1:nSims){

for (i in 1:4){xObs[i]=rnorm(1,mu[i],sigma/sqrt(N1))}

if (xObs>xObs){

SelectedArm=2

if (xObs*sqrt(N1)/sigma>xObs*sqrt(N1)/sigma-cR) {SelectedArm=1}

}

if (xObs<=xObs){

SelectedArm=3

if (xObs*sqrt(N1)/sigma>xObs*sqrt(N1)/sigma-cR) {SelectedArm=4}

}

MaxRsp = xObs[SelectedArm]

x2=rnorm(1,mu[SelectedArm], sigma/sqrt(N2))

FinalxAve = (MaxRsp*N1+x2*N2)/(N1+N2)

x0Ave = rnorm(1,mu0,sigma/sqrt(N1+N2))

TestZ=(FinalxAve-x0Ave)*sqrt((N1+N2)/2)/sigma

if (TestZ >= c_alpha) {Power=Power+1/nSims}

for (i in 1:4){if (SelectedArm==i){SelProb[i]=SelProb[i]+1/nSims}}

} #End of iSim

TotalN=4*N1+2*N2

return (c("power=", Power, "Sample Size=", TotalN, "Selection Prob=", SelProb))

}

Invoke the function:

mu=c(0,0,0,0)

FourPlus1AddArmDesign(nSims=100000, N1=100, N2=100, c_alpha = 2.267,cR = 0.55652, mu0=0, sigma=1)

mu=c(0.4,0.58,0.7,0.45)

FourPlus1AddArmDesign(nSims=10000, N1=116, N2=116, c_alpha = 2.267,cR = 0.55652, mu0=0.35, sigma=0.9)

mu=c(0,0,0,0)

FourPlus1AddArmDesign(nSims=100000, N1=30, N2=70, c_alpha = 2.267,cR = 0.55652, mu0=0, sigma=1)

# Chapter 9

## Chapter 9-BiomarkerEnrichmentDesign

## mup, mun = treatment effects for biomarker positive and negative populations

## sigma = common standard deviation for treatment effect

## Np1, Np2, Nn1, Nn2 = sample size for biomarlker positive and negative at stages 1 and 2

## c_alpha = critical value for rejection

## powerP, powerO, power = powers for biomarker positive, negative, and either of them

BiomarkerEnrichmentDesign=function(nSims=10000,mup=0.5, mun=0, sigma=1, Np1=100, Np2=100, Nn1=100, Nn2=100, c_alpha = 2.2)

{

Np=Np1+Np2; Nn=Nn1+Nn2; TotalN=Np+Nn; N1=Np1+Nn1; N2=Np2+Nn2

power=0; powerO=0; powerP=0; aveN=0

for (iSim in 1:nSims){

xp1=rnorm(1,mup,sigma/sqrt(Np1))

xn1=rnorm(1,mun,sigma/sqrt(Nn1))

x1=(xp1*Np1+xn1*Nn1)/N1

tp1=xp1*sqrt(Np1)/sigma

t1=x1*sqrt(N1)/sigma

xp2=rnorm(1,mup,sigma/sqrt(Np2))

xn2=rnorm(1,mun,sigma/sqrt(Nn2))

x2=(xp2*Np2+xn2*Nn2)/N2

tp2=xp2*sqrt(Np2)/sigma

t2=x2*sqrt(Np)/sigma

tp=sqrt(Np1/Np)*tp1+sqrt(Np2/Np)*tp2

t=sqrt(Np1/N1)*t1+sqrt(Np2/N2)*t2

aveN=aveN+(Np1+Np2+Nn1+Nn2)/nSims

if (tp1>t1) aveN=aveN-Nn2/nSims

if (tp1>t1 & tp>c_alpha) {powerP=powerP+1/nSims}

if (tp1<=t1 & t>c_alpha) {powerO=powerO+1/nSims}

if ((tp1>t1 & tp>c_alpha)|(tp1<=t1 & t>c_alpha)) {power=power+1/nSims}

} #End of iSim

return (c("Average N=", aveN, "power=", power, "powerP=", powerP, "powerO=", powerO))

}

Invoke the function

# # Simulation of Global H0 Design 1 (Nn1=Nn2=200)

BiomarkerEnrichmentDesign(nSims=1000000,mup=0, mun=0, sigma=1.2, Np1=200, Np2=200, Nn1=200, Nn2=200, c_alpha = 2.13)

# Simulation of Ha: Design 1 (Nn1=Nn2=200)

BiomarkerEnrichmentDesign(nSims=100000,mup=0.2, mun=0.05, sigma=1.2, Np1=200, Np2=200, Nn1=200, Nn2=200, c_alpha = 2.13)

# Simulation of Global H0 Design 2 (Nn1=Nn2=105)

BiomarkerEnrichmentDesign(nSims=1000000,mup=0, mun=0, sigma=1.2, Np1=210, Np2=210, Nn1=105, Nn2=105, c_alpha = 2.39)

# Simulation of H: Design 2 (Nn1=Nn2=105)

BiomarkerEnrichmentDesign(nSims=100000,mup=0.2, mun=0.05, sigma=1.2, Np1=210, Np2=210, Nn1=105, Nn2=105, c_alpha = 2.39)

## Chapter 9-Biomarker-Informed Design with Hierarchical Model

## R-Function 9.2: Biomarker-Informed Design with Hierarchical Model

## u0y[j], u0x[j], rhou, suy sux = parameters of the parent model in group j

## rho, sy, sx = parameters in the lower level model

## Zalpha =critical point for rejection

## nArms = number of active groups

## N1 and N = sample size per group at interim and final analyses

## probWinners = Probability of selecting the arm as the winner

library(mvtnorm)

powerBInfo=function(uCtl, u0y, u0x, rhou, suy, sux ,rho, sy, sx, Zalpha, N1, N, nArms, nSims)

{

uy=rep(0,nArms); ux=rep(0,nArms); probWinners=rep(0,nArms); power = 0

varcov0=matrix(c(suy^2,rhou*suy*sux,rhou*suy*sux, sux^2),2,2)

varcov=matrix(c(sy^2, rho*sy*sx, rho*sx*sy, sx^2),2,2)

for (i in 1: nSims) {

winnerMarker= -Inf

for (j in 1: nArms) {

u=rmvnorm(1,mean=c(u0y[j],u0x[j]), sigma=varcov0)

uy[j]=u; ux[j]=u

dataStg1=rmvnorm(N1, mean=c(uy[j], ux[j]), sigma=varcov)

meanxMarker=mean(dataStg1[,2])

if (meanxMarker>winnerMarker)

{winner=j; winnerMarker=meanxMarker; winnerY=dataStg1[,1]}

} ## End of j ##

for (j in 1:nArms) {if (winner==j) {probWinners[j]=probWinners[j]+1/nSims}}

yStg1=winnerY

yStg2=rnorm(N-N1, mean=uy[winner], sd=sy)

yTrt=c(yStg1+yStg2)

yCtl=rnorm(N, mean=uCtl, sd=sy)

tValue=t.test(yTrt,yCtl)\$statistic

if (tValue>=Zalpha) {power=power+1/nSims}

} ## End of i ##

return (c(power, probWinners))

}

Invoke the function

## Determine critical value Zalpha for alpha (power) =0.025 ##

u0y=c(0,0,0); u0x=c(0,0,0)

powerBInfo(uCtl=0, u0y, u0x, rhou=1, suy=0, sux=0 ,rho=1, sy=4, sx=4, Zalpha=2.772, N1=100, N=300, nArms=3, nSims=100000)

## Power simulation ##

u0y=c(1,0.5,0.2)

u0x=c(2,1,0.5)

powerBInfo(uCtl=0, u0y, u0x, rhou=0.2, suy=0.2, sux=0.2 ,rho=0.2, sy=4, sx=4, Zalpha=2.772, N1=100, N=300, nArms=3, nSims=5000)

# Chapter 10

##Chapter 10-Randomized Play-the-Winner Design with Binary endpoint

## Randomized Play-the-Winner Design with Binary endpoint

## a0, b0 = the initial numbers of balls for the two treatments A and B

## a1, b1 = balls added to the urn if a response is observed in arm A or arm B.

## RR1, RR2 = the response rates in groups 1 and 2

## nSbjs = total number of subjects (two groups combined)

## nMin (>0) = the minimum sample-size per group

## nAnlys = number of analyses

## interim analyses are designed for randomization adjustment

## aveP1 and aveP2 = the average response rates in groups 1 and 2

## Power = probability of the test statistic > Zc

RPWBinary = function(nSims=1000, Zc=1.96, nSbjs=200, nAnlys=3,

RR1=0.2, RR2=0.3, a0=1, b0=1, a1=1, b1=1, nMin=1) {

set.seed(21823)

Power=0; aveP1=0; aveP2=0; aveN1=0; aveN2=0

for (isim in 1:nSims) {

nResp1=0; nResp2=0; N1=0; N2=0

nMax=nSbjs-nMin

a=a0; b=b0; r0=a/(a+b)

for (iSbj in 1:nSbjs) {

nIA=round(nSbjs/nAnlys)

if (iSbj/nIA==round(iSbj/nIA)) {r0=a/(a+b)}

if ((rbinom(1,1,r0)==1 & N1<nMax) | N2>=nMax) {

N1=N1+1

if (rbinom(1,1,RR1)==1) {nResp1=nResp1+1; a=a+a1}

}

else

{

N2=N2+1

if (rbinom(1,1,RR2)==1) { nResp2=nResp2+1; b=b+b1 }

}

}

aveN1=aveN1+N1/nSims; aveN2=aveN2+N2/nSims

p1=nResp1/N1; p2=nResp2/N2

aveP1=aveP1+p1/nSims; aveP2=aveP2+p2/nSims

stdErr2=p1*(1-p1)/N1+p2*(1-p2)/N2

TS = (p2-p1)/sqrt(stdErr2)

if (TS>Zc) {Power=Power+1/nSims}

}

return (cbind(nSbjs, aveN1, aveN2, aveP1, aveP2, Zc, Power))

}

RPWBinary(nSims=10000, Zc=2.01, nSbjs=200, nAnlys=1, RR1=0.4, RR2=0.4, a0=1, b0=1,a1=0, b1=0, nMin=1)

RPWBinary(nSims=1000, Zc=2.01, nSbjs=200, nAnlys=1, RR1=0.3, RR2=0.5, a0=1, b0=1,a1=0, b1=0, nMin=1)

RPWBinary(nSims=10000, Zc=2.035, nSbjs=200, nAnlys=5, RR1=0.4, RR2=0.4, a0=2, b0=2,a1=1, b1=1, nMin=1)

RPWBinary(nSims=1000, Zc=2.035, nSbjs=200, nAnlys=5, RR1=0.3, RR2=0.5, a0=2, b0=2,a1=1, b1=1, nMin=1)

RPWBinary(nSims=10000, Zc=2.010, nSbjs=200, nAnlys=200, RR1=0.4, RR2=0.4, a0=1, b0=1,a1=0, b1=0, nMin=1)

RPWBinary(nSims=1000, Zc=2.010, nSbjs=200, nAnlys=200, RR1=0.3, RR2=0.5, a0=1, b0=1,a1=0, b1=0, nMin=1)

## Chapter 10-Randomized Play-the-Winner Design with Normal Endpoint

RARNorm=function(nSims=1000, nPts=100, nArms=5, b=1, m=1, CrtMax=1.96){

AveN=rep(0,nArms); rP=rep(0,nArms); AveU=rep(0,nArms); PowerMax=0

for (isim in 1:nSims){

uObs=rep(0,nArms); cuObs=rep(0,nArms); Ns=rep(0,nArms); Crp=rep(0,nArms)

for (iSubject in 1:nPts) {

for (i in 1:nArms){rP[i]=a0[i]+b*uObs[i]**m}

Suma=0; for (i in 1:nArms){Suma=Suma+rP[i]}

for (i in 1:nArms){rP[i]=rP[i]/Suma}

CrP=rep(0,nArms)

for (iArm in 1:nArms){

for (i in 1:iArm){CrP[iArm]=CrP[iArm]+rP[i]}

}

rn=runif(1,0,1); cArm=1

for (iArm in 2:nArms){

if (CrP[iArm-1]<rn & rn<CrP[iArm]){cArm=iArm}

}

Ns[cArm]= Ns[cArm]+1

## For Normal response

u=rnorm(1, us[cArm],s[cArm])

cuObs[cArm]=cuObs[cArm]+u

for (i in 1:nArms){uObs[i]=cuObs[i]/max(Ns[i],1)}

} ## End of iSubject

se2=0

## Assume sigma known for simplicity

for (i in 1:nArms){se2=se2+s[i]**2/max(Ns[i],1)*2/nArms}

uMax=uObs

for (i in 1:nArms){if (uObs[i]>=uMax){iMax=i} }

TSmax=(uObs[iMax]-uObs)*(Ns+Ns[iMax])/2/(nPts/nArms)/se2**0.5

if (TSmax>CrtMax){PowerMax=PowerMax+1/nSims}

for (i in 1:nArms) {

AveU[i]=AveU[i]+uObs[i]/nSims

AveN[i]=AveN[i]+Ns[i]/nSims

}

} # End of iSim

return (c("nPts=", nPts, "Power=", PowerMax, "AveN=", AveN, "AveU=", AveU))

}

Invike the function

a0=rep(1,5); us=rep(0.06,5); s=rep(.18,5)

RARNorm(nSims=10000, nPts=375, nArms=5, b=1, m=1, CrtMax=2.01)

a0=rep(1,5); us=c(0.06, 0.12, 0.13, 0.14, 0.15); s=rep(.18,5)

RARNorm(nSims=10000, nPts=375, nArms=5, b=1, m=1, CrtMax=2.01)

# Chapter 11

## Chapter 11-CRM     ## Trial simulation trial using CRM.

## b = model parameter in (11.1)

## aMin and aMax = the upper and lower limits of parameter a.

CRM = function(nSims=100, nPts=30, nLevels=10, b=100, aMin=0.1,      aMax=0.3, MTRate=0.3, nIntPts=100, ToxModel="Skeleton", nPtsForStop=6)

{

g0 = c(1); PercentMTD= c(1)      DLTs=0; AveMTD=0; VarMTD=0; AveTotalN=0

dx=(aMax-aMin)/nIntPts     for (k in 1:nIntPts) {g0[k]=g[k]}

PercentMTD=rep(0,nLevels)

for (iSim in 1:nSims) {

for (k in 1:nIntPts) {g[k]=g0[k]}

nPtsAt=rep(0,nLevels); nRsps=rep(0,nLevels)

iLevel=1; PreLevel=1     for (iPtient in 1:nPts) {

TotalN=iPtient

TotalN=iPtient

iLevel=min(iLevel, PreLevel+1); #Avoid dose-jump

iLevel=min(iLevel, nLevels)

PreLevel=iLevel;

Rate=RRo[iLevel]

nPtsAt[iLevel]=nPtsAt[iLevel]+1

r= rbinom(1, 1, Rate)

nRsps[iLevel]=nRsps[iLevel]+r

# Posterior distribution of a

c=0

for (k in 1:nIntPts) {

ak=aMin+k*dx

if (ToxModel=="Logist") { Rate=1/(1+b*exp(-ak*doses[iLevel]))}

if (ToxModel=="Skeleton") { Rate=skeletonP[iLevel]^exp(ak)}

if (r>0) {L=Rate}

if (r <= 0) {L=1-Rate}     g[k]=L*g[k]; c=c+g[k]*dx

}

for (k in 1:nIntPts) {g[k]=g[k]/c}     # Predict response rate and current MTD

MTD=iLevel; MinDR=1; RR=rep(0, nLevels)

for (i in 1:nLevels) {

for (k in 1:nIntPts) {     ak=aMin+k*dx

if (ToxModel=="Logist") { RR[i]= RR[i]+1/(1+b*exp(-ak*doses[i]))*g[k]*dx}     if (ToxModel=="Skeleton") { RR[i]= RR[i]+skeletonP[i]^exp(ak)*g[k]*dx}

}

DR=abs(MTRate-RR[i])

if (DR<MinDR){MinDR=DR; iLevel=i;MTD = i}

}

if (nPtsAt[iLevel] >= nPtsForStop) {

PercentMTD[iLevel] =PercentMTD[iLevel]+1/nSims

break()

}

}      for (i in 1:nLevels) {DLTs=DLTs+nRsps[i]/nSims}

AveMTD=AveMTD+MTD/nSims

VarMTD=VarMTD+MTD^2/nSims

AveTotalN=AveTotalN+TotalN/nSims

}

SdMTD=sqrt(VarMTD-AveMTD^2)

return (cbind(AveTotalN, nLevels, AveMTD, SdMTD, DLTs))

}

#Logistic Model

g= c(1); doses = c(1)

RRo = c(0.01,0.02,0.03,0.05,0.12,0.17,0.22,0.4)

for (k in 1:100) {g[k]=1}; # Flat prior     for (i in 1:8) {doses[i]=i}

CRM(nSims=500, nPts=30, nLevels=8, b=150, aMin=0, aMax=3,

MTRate=0.17, ToxModel="Logist", nPtsForStop=6)

#Skeleton Model

g = c(1); doses = c(1)

RRo = c(0.01,0.02,0.03,0.05,0.12,0.17,0.22,0.4)

skeletonP=c(0.01,0.02,0.04,0.08,0.16,0.32,0.4,0.5)

for (k in 1:100) {g[k]=1}; # Flat prior

CRM(nSims=500, nPts=20, nLevels=8, aMin=-0.1, aMax=2,

MTRate=0.17, ToxModel="Skeleton", nPtsForStop=6)

# Chapter 12

## Chapter 12 - Conditional Power

## EP = endpoint, "normal" or "binary"

## Model = "MSP", "MPP", or "MINP"

## alpha2 = 2nd stage stopping boundary on p-scale

## nStg2 = sample size at 2nd stage

## ux, uy = observed responses in the two groups

## p1 = interim p-value

## w1 = weight in MINP only

ConPower=function(EP="normal", Model="MINP", alpha2=0.0226, ux=0.2, uy=0.4, sigma=1, nStg2=100, p1=0.8, w1=0.707)

{

u=(ux+uy)/2

if (EP=="binary"){sigma=(u*(1-u))**0.5}

w2=sqrt(1-w1*w1)

eSize=(uy-ux)/sigma

if (Model=="MSP"){BFun=qnorm(1-max(0.0000001,alpha2-p1))}

if (Model=="MPP"){BFun=qnorm(1-alpha2/p1)}

if (Model=="MINP"){BFun=(qnorm(1-alpha2)- w1*qnorm(1-p1))/w2}

cPower=1-pnorm(BFun-eSize*sqrt(nStg2/2))

return (cPower)

}

Invoke the function:     ConPower(EP="binary", Model="MSP", alpha2=0.2050, ux=0.2, uy=0.4, nStg2=100, p1=0.1)     ConPower(EP="binary", Model="MPP", alpha2=0.0043, ux=0.2, uy=0.4, nStg2=100, p1=0.1)

ConPower(EP="binary", Model="MINP", alpha2=0.0226, ux=0.2, uy=0.4, nStg2=100, p1=0.1, w1=0.707)

ConPower(EP="normal", Model="MSP", alpha2=0.2050, ux=0.2, uy=0.4, sigma=1, nStg2=200, p1=0.1)

ConPower(EP="normal", Model="MPP", alpha2=0.0043, ux=0.2, uy=0.4, sigma=1, nStg2=200, p1=0.1)

ConPower(EP="normal", Model="MINP", alpha2=0.0226, ux=0.2, uy=0.4, sigma=1, nStg2=200, p1=0.1, w1=0.707)

## Chapter 12 - New sample size at stage 2

## New Sample Size for stage 2

## trtDiff = observed treatment difference

## sigma = observed standard deviation

## p1 = pvalue at the interim analysis

## NId = noninferiority margin

## cPower = targeted conditional power

## nStg2, n2New = initial and the new sample sizes for stage 2

NewN2byConditionalPowerForTwoStageSSR=function(trtDiff=0.5, sigma=1, NId, p1, alpha2, nStg2, cPower, method="MINP", w1)     {

w2=sqrt(1-w1*w1)

if (method=="MSP") {BFun = qnorm(1-max(0.000001,alpha2-p1))}

if (method=="MPP") {BFun = qnorm(1- min(0.999999,alpha2/p1))}

if (method=="MINP"){BFun = (qnorm(1-alpha2)- w1*qnorm(1-p1))/w2}

eSize=(trtDiff+NId)/sigma     n2New=max(2*((BFun-qnorm(1-cPower))/eSize)^2, nStg2)

return (c("New sample at stage 2=", round(n2New)))

}

NewN2byConditionalPowerForTwoStageSSR(trtDiff=0.2, sigma=1.2, NId=0, p1=0.1, alpha2=0.0226, nStg2=50, cPower=0.9, method="MINP", w1=0.707)

# Chapter 13

## Chapter 13 - CRM Monitoring

## Trial monitoring using CRM.

## b = model parameter in (11.1)

## aMin and aMax = the upper and lower limits of parameter a.

CRM = function(nSims=100, nPts=30, nLevels=10, b=100, aMin=0.1,

aMax=0.3, MTRate=0.3, nIntPts=100, ToxModel="Skeleton", nPtsForStop=6)      {

g0 = c(1); dx=(aMax-aMin)/nIntPts

for (k in 1:nIntPts) {g0[k]=g[k]}

## for (iSim in 1:nSims) {

for (k in 1:nIntPts) {g[k]=g0[k]}

nPtsAt=rep(0, nLevels); nRsps=rep(0, nLevels);

iLevel=1; PreLevel=1

for (iPtient in 1:nPts) {

TotalN=iPtient

TotalN=iPtient

iLevel=min(iLevel, PreLevel+1); #Avoid dose-jump

iLevel=min(iLevel, nLevels)

PreLevel=iLevel

Rate=RRo[iLevel]

nPtsAt[iLevel]=nPtsAt[iLevel]+1

## r= rbinom(1, 1, Rate) ## random Response

r=ObsResponses[iPtient] ## Observed response

nRsps[iLevel]=nRsps[iLevel]+r

## Posterior distribution of a

c=0

for (k in 1:nIntPts) {

ak=aMin+k*dx

if (ToxModel=="Logist") { Rate=1/(1+b*exp(-ak*doses[iLevel]))}

if (ToxModel=="Skeleton") { Rate=skeletonP[iLevel]^exp(ak)}

if (r>0) {L=Rate}

if (r <= 0) {L=1-Rate}

g[k]=L*g[k]; c=c+g[k]*dx

}

for (k in 1:nIntPts) {g[k]=g[k]/c}

# Predict response rate and current MTD

MTD=iLevel; MinDR=1

RR = rep(0,nLevels)

for (i in 1:nLevels) {

for (k in 1:nIntPts) {

ak=aMin+k*dx

if (ToxModel=="Logist") { RR[i]= RR[i]+1/(1+b*exp(-ak*doses[i]))*g[k]*dx}

if (ToxModel=="Skeleton") { RR[i]= RR[i]+skeletonP[i]^exp(ak)*g[k]*dx}

}

DR=abs(MTRate-RR[i])

if (DR<MinDR){MinDR=DR; iLevel=i;MTD = i}

}

if (nPtsAt[iLevel] >= nPtsForStop) {

break()

}

}

## } ## End of iSim

return (c("Dose Level=", PreLevel, "Predicted MTD=", MTD))

}

Invoke the function:

#Logistic Model

g= c(1); doses = c(1)

RRo = c(0.01,0.02,0.03,0.05,0.12,0.17,0.22,0.4)

ObsResponses=c(0,0,0,0,0,1,0,0,1,0,0)

for (k in 1:100) {g[k]=1}; # Flat prior

for (i in 1:8) {doses[i]=i}

CRM(nSims=500, nPts=11, nLevels=8, b=150, aMin=0, aMax=3,

MTRate=0.17, ToxModel="Logist", nPtsForStop=6)

#Skeleton Model

g = c(1); doses = c(1)

RRo = c(0.01,0.02,0.03,0.05,0.12,0.17,0.22,0.4)

skeletonP=c(0.01,0.02,0.04,0.08,0.16,0.32,0.4,0.5)

ObsResponses=c(0,0,0,0,0,1,0,0,1,0,0)

for (k in 1:100) {g[k]=1}; # Flat prior

CRM(nSims=500, nPts=11, nLevels=8, aMin=-0.1, aMax=2,

MTRate=0.17, ToxModel="Skeleton", nPtsForStop=6)

# Chapter 14

## Chapter 14-GSD pValue by Simulation

## GSD P-value by Simulation based stagewise ordering

## u0, u1 = means for two treatment groups

## sigma0, sigma1 = standard deviations for two treatment groups

## n0Stg1, n1Stg1, n0Stg2, n1Stg2 = sample sizes for the two groups at stages 1 and 2

## alpha1 = efficacy stopping boundary at stage 1

## w1squared = weight squared for MINP only at stage 1

## nSims = the number of simulation runs

##pValue, pValue2 = p-value and conditional p-value at stage 2

Pvalue TwoStageGSDwithNormalEndpoint=function(u0,u1,sigma0,sigma1, n0Stg1, n1Stg1, n0Stg2, n1Stg2, alpha1, w1squared=0.5, t2Obs=1.96, nSims=1000000)

{

pValue2=0; M=0; w1=sqrt(w1squared);

for (i in 1:nSims) {

y0Stg1=rnorm(1, u0, sigma0/sqrt(n0Stg1))

y1Stg1=rnorm(1, u1, sigma1/sqrt(n1Stg1))

z1=(y1Stg1-y0Stg1)/sqrt(sigma1**2/n1Stg1+sigma0**2/n0Stg1)     t1=1-pnorm(z1)

if(t1>alpha1) {

M=M+1

y0Stg2=rnorm(1, u0, sigma0/sqrt(n0Stg2))

y1Stg2=rnorm(1, u1, sigma1/sqrt(n1Stg2))

z2=(y1Stg2-y0Stg2)/sqrt(sigma1**2/n1Stg2+sigma0**2/n0Stg2)

t2=w1*z1+sqrt(1-w1*w1)*z2

if(t2>=t2Obs) {pValue2=pValue2+1}     }

}

pValue=alpha1+pValue2/nSims     pValue2=pValue2/M

return (c("pValue=", pValue, "pValue2=", pValue2))

}

Invoke the function:

## Example 14.1

PvalueTwoStageGSDwithNormalEndpoint(u0=0.05,u1=0.05, sigma0=0.22, sigma1=0.22, n0Stg1=120, n1Stg1=120, n0Stg2=120, n1Stg2=120, alpha1=0.01, t2Obs=1.96)