Dose response model and biased coin up-and-down design estimating EV90. A simulation study with a boundary condition and addition of auxilliary data points. PAVA estimators and 2 parameter isotonic regression.
Experimental Spreadsheet from study with reporting sheet and random numbers sheet, latter should be blinded during the course of the experiment. The study had a minimum of 25 patients.
#Upper boundary of 30 (mL), steps of 2 (mL)
#Estimat for number of patients
ANTAL<-NULL
MU_1<-NULL
MU_2<-NULL
MU_3<-NULL
EV90<-NULL
library(Iso)
library(drc)
iter<-1
while(iter<=10000){
#MEV
MEV<-0.9
#Start dose
V<-20
#From previous study start dose is probably around MEV92
P<-0.92-(30-V)*0.014
#Defining vectors for simulation
succes<-0
coins<-0
kt<-1
#Increase in dose if failure
dosedeltaNeg<-2
dosedeltaPos<-2
while(V[kt]<30){
coins[kt]<-runif(1)
if(coins[kt]>P[kt]){
#If we observe a failure then increase dosis
succes[kt]<-0
V[kt+1]<-V[kt]+dosedeltaNeg
#We estimate effect a dose increase based on empirical data
P[kt+1]<-P[kt]+0.014
}
else {
#If we do not observe a failure then
#there is a (1-MEV)/MEV probability of a dose decrease
succes[kt]<-1
V[kt+1]<-V[kt]-dosedeltaPos*(runif(1)<=((1-MEV)/MEV))
#We estimate effect a dose decrease based on empirical data
P[kt+1]<-P[kt]+(V[kt+1]-V[kt])*0.014
}
kt<-kt+1
}
V<-c(V[1:(length(V)-1)],rep(30,25))
ssum<-sum(succes)
succes<-c(succes,rep(1,23),0,0)
datf<-data.frame(succes=succes,V=V)
#par(mfrow=c(2,1))
#plot(V,type="l",main=cat("Number of successes:", ssum))
tbl<-aggregate(x=datf$succes,by=list(datf$V),FUN=mean)
tbln<-aggregate(x=datf$succes,by=list(datf$V),FUN=length)
tbls<-aggregate(x=datf$succes,by=list(datf$V),FUN=sum)
tbl$cnt<-tbln$x
tbl$sum<-tbls$x
#Construction of PAVA estimates
#Volume
tbl$Group.1
#Isotone regression estimate
isoreg<-pava(y=tbl$x,w=tbl$cnt,decreasing=FALSE)
x_mu1<-max((1:dim(tbl)[1])[isoreg<=MEV])
mu1<-tbl$Group.1[x_mu1]
x_mu2<-which.min(abs(isoreg-MEV))
mu2<-tbl$Group.1[x_mu2]
mu3<-((MEV-isoreg[x_mu1])/(isoreg[x_mu1+1]-isoreg[x_mu1]))*(tbl$Group.1[x_mu1+1]-mu1)+mu1
MU_1[iter]<-mu1
MU_2[iter]<-mu2
MU_3[iter]<-mu3
tbl$id<-rep(1,length(tbl$cnt))
fourpar <- drm(sum/cnt ~ Group.1, id,weight=cnt,data = tbl, fct = LL2.2(names = c("Slope", "ED90")), type = "binomial")
try(EV90[iter]<-exp(ED(fourpar, c(90), interval = "delta",display="FALSE")[1]),silent="TRUE")
#plot(isoreg,xaxt="n",main=paste("PAVA, mu1:",round(mu1,2),"mu2",round(mu2,2),"mu3",round(mu3,2)),sub=paste("Number of successes:", ssum),xlab="Doses", ylab="EV")
#isoregplot<-pava(y=tbl$x,w=tbl$cnt,decreasing=FALSE,stepfun=TRUE)
#plot(isoregplot,add=TRUE,verticals=FALSE,pch=20,col.points="red",xaxt="n")
#axis(1, labels=as.character(tbl$Group.1), at=1:length(tbl$Group.1))
ANTAL[iter]<-length(P)
iter<-iter+1
}
summary(ANTAL)
quantile(ANTAL,c(0.9,0.95,0.99),na.rm=TRUE)
summary(MU_3)
quantile(MU_3,c(0.025,0.975),na.rm=TRUE)
#Example plot
par(mfrow=c(2,1))
ssum<-sum(succes)
plot(V,type="l",main=paste("Number of successes:",ssum))
plot(isoreg,xaxt="n",main=paste("PAVA, mu1:",round(mu1,2),"mu2:",round(mu2,2),"mu3:",round(mu3,2)),sub=paste("Number of successes:", ssum),xlab="Doses", ylab="EV")
isoregplot<-pava(y=tbl$x,w=tbl$cnt,decreasing=FALSE,stepfun=TRUE)
plot(isoregplot,add=TRUE,verticals=FALSE,pch=20,col.points="red",xaxt="n")
axis(1, labels=as.character(tbl$Group.1), at=1:length(tbl$Group.1))