R Code for Introductory Adaptive Design with R
## 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-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-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-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-ADwithImcompletePairs
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)
## Chapter 6-NIADwithPairedBinaryData
## 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-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[1]; 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-FourPlusOneAddArmDesign
##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[2]>xObs[3]){
SelectedArm=2
if (xObs[1]*sqrt(N1)/sigma>xObs[2]*sqrt(N1)/sigma-cR) {SelectedArm=1}
}
if (xObs[2]<=xObs[3]){
SelectedArm=3
if (xObs[4]*sqrt(N1)/sigma>xObs[3]*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-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[1]; ux[j]=u[2]
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-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[1]
for (i in 1:nArms){if (uObs[i]>=uMax){iMax=i} }
TSmax=(uObs[iMax]-uObs[1])*(Ns[1]+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-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 - 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 - 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-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)