Skip to main content

Up-and-down dose response study design

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. 



 #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))  

Comments

Popular posts from this blog

HackRF on Windows 8

This technical note is based on an extract from thread. I have made several changes and added recommendations. I have experienced lot of latency using GnuRadio and HackRF on Pentoo Linux, so I wanted to try out GnuRadio on Windows.



HackRF One is a transceiver, so besides SDR capabilities, it can also transmit signals, inkluding sweeping a given range, uniform and Gaussian signals. Pentoo Linux provides the most direct access to HackRF and toolboxes. Install Pentoo Linux on a separate drive, then you can use osmocom_siggen from a terminal to transmit signals such as near-field GSM bursts, which will only be detectable within a meter.









Installation of MGWin and cmake: Download and install the following packages:
- MinGW Setup (Go to the Installer directory and download setup file)
- CMake (I am using CMake 3.2.2 and I installed it in C:\CMake, this path is important in the commands we must send in the MinGW shell)
Download and extract the packages respectively in the path C:\MinGW\msys\…

Example: Beeswarm plot in R

library(foreign)

data <- read.dta("C:/Users/hellmund/Documents/MyStataDataFile.dta")

names(data)

install.packages('beeswarm')

library(beeswarm)

levels(data$group)

png(file="C:/Users/hellmund/Documents/il6.png", bg="transparent")

beeswarm(data$il6~data$group,data=data, method=c("swarm"),pch=16,pwcol=data$Gender,xlab='',ylab='il6',ylim=c(0,20))

legend('topright',legend=levels(data$Gender),title='Gender',pch=16,col=2:1)

boxplot(data$il6~data$group, data=data, add = T, names = c("","",""), col="#0000ff22")

dev.off()

Real world split-plot designs

Google Earth picture from a blog on statistics. A real world example near Christchurch (NZ) of a split-plot design. Today things have completely changed on location as the forest has grown considerably. Google Earth coordinate link.