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. 



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

Comments

Popular posts from this blog

Alder/korrekt århundrede udfra cpr nummer

De fleste, der arbejder med registre eller databaser, står ofte med problemstillingen, at alder er uoplyst, medens cpr-nummer er kendt. Hvordan regner man den ud? Følgende regel er gældende: Hvis syvende ciffer er 0, 1, 2 eller 3 er man født i det 20. århunderede (1900-tallet) Ligeledes, hvis syvende ciffer er 4 eller 9, og årstallet (femte og sjette ciffer) er større end eller lig 37. Endelig er man født i det 19. århundrede (1800-tallet) hvis syvende ciffer er 5, 6, 7 eller 8 og årstallet er større end eller lig 58. Nedenfor finder du eksempel i SAS kode: En lille makro, der udover fødselsdato også udregner køn samt den præcise alder givet datovariabel. Kilde: Opbygning af CPR nummeret, cpr.dk proc format library=work; value gender 0="Female" 1="Male" ; run; %macro agefromCPR(cpr,datevar=inddto,birthvar=birth,agevar=age); dy_temp=input(substrn(&cpr,1,2),2.); mt_temp=input(substrn(&cpr,3,2),2.); yr_temp=input(substrn(&cpr,5,2),

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

Comorbidity indexes in SQL

Generating Elixhauser comorbidity index from Danish National Health Register as relational database. ( ICD 10 Coding  in SAS) A lookup-table based version of Charlson comorbidity index I made in SQL. A similar approach can be applied to Elixhauser. SELECT V_CPR, MAX(EI1)+MAX(EI2)+MAX(EI3)+MAX(EI4)+MAX(EI5)+ MAX(EI6)+MAX(EI7)+MAX(EI8)+MAX(EI9)+MAX(EI10)+ MAX(EI11)+MAX(EI12)+MAX(EI13)+MAX(EI14)+MAX(EI15)+ MAX(EI16)+MAX(EI17)+MAX(EI18)+MAX(EI19)+MAX(EI20)+ MAX(EI21)+MAX(EI22)+MAX(EI23)+MAX(EI24)+MAX(EI25)+ MAX(EI26)+MAX(EI27)+MAX(EI28)+MAX(EI29)+MAX(EI30)+MAX(EI31) AS Elixhauser FROM (SELECT V_CPR, -- Congestive Heart Failure CASE WHEN DIAG LIKE 'DI099%' OR DIAG LIKE 'DI110%' OR DIAG LIKE 'DI130%' OR DIAG LIKE 'DI132%' OR DIAG LIKE 'DI255%' OR DIAG LIKE 'DI420%' OR DIAG LIKE 'DI425%' OR DIAG LIKE 'DI426%' OR DIAG LIKE 'DI427%' OR DIAG LIKE 'DI428%' OR DIAG LIKE 'DI429%' OR D