source("/Users/pm/s/first.s") dir <- "/Users/pm/book/ctrial/pm/3/" setwd(dir) ## dose range n <- 16 # n doses, including placebo, dose=1 doses <- 1:n ## truth, using parameters in f(.) below scenario0 <- c(0,6,1.5) # curve 1, flop scenario1 <- c(2,4,1.5) # moderate effect scenario4 <- c(4,6,1.5) # curve 4, with incerase 4 points over placebo scenario5 <- c(4,10,1) f <- function(d,b) { ## true mean response: ## b[1]=scale, b[2]=centering of logistic, b[3]=logistic slope return(b[1]*1/(1+exp(-(d-b[2])*b[3])) ) } ## model p <- 2 G <- rbind(c(1,1),c(0,1)) Z <- matrix(c(1,0),nrow=1,ncol=2) sig <- 1.0 tau <- 0.1 ## prior with fake data d0 <- doses y0 <- f(d0,c(2,8,1.5)) n0 <- rep(0.1,n) # weight of prior data (1.0 = one data point) sim.trial <- function(n.mx=100,b=scenario4,plt=F) {# simulate one realization of the entire trial ## 0. initialize N <- 0 M <- 100 Y <- NULL if (plt) # start plot plot(doses,f(doses,b),type="l",lty=3,xlab="DOSES", ylab="", bty="l",ylim=c(-1,6)) trajED <- NULL trajm <- NULL pi <- rep(1,n)/n # start with uniform allocaction prob's trajpi <- pi k <- 3 # batch size repeat{ ## 1. next batch of patients di <- sample(doses, k, prob=pi, replace=T) m <- f(di,b) yi <- rnorm(k, m=m, s=sqrt(sig)) Y <- rbind(Y,cbind(di,yi)) if (plt) points(Y[,1],Y[,2],pch=19,cex=0.5) N <- N+k ## 2. update posterior inference fb <- ffbs(n,p,Y,G,Z,sig,tau,y0=y0, n0=n0, M=M) if (plt) lines(doses,fb$m,lty=2,lwd=1) ED <- apply(fb$X[1,,],2,ed95) EDs <- modex(ED) trajED <- c(trajED, EDs) trajm <- rbind(trajm, fb$m) ## 3. check for stopping criteria if (N>n.mx){ stp <- "max"; break } p3 <- 1-pnorm(3,mean=fb$m,s=fb$s) # prob of success if (N>10){ if (max(p3) < 0.05){ # stop for futility cat("\n Stopping for futility.\n"); stp <- "fut"; break } if (p3[EDs] > 0.95){ # stop for victory cat("\n Stopping for victory.\n"); stp <- "vic"; break } } # allocation probs: dose j proportional to p(ED95=j | Y) pi <- sqrt( hist(ED,breaks <- c(0.5,doses +.5), plot=F)$dens ) pi <- pi/sum(pi) if (pi[1] < 0.1){ ## enforce min 10% at placebo pi[-1] <- pi[-1]/sum(pi[-1])*0.9 pi[1] <- 0.1 } trajpi <- rbind(trajpi,pi) } if (plt) lines(doses,fb$m,lty=1,lwd=3) return(list(Y=Y,b=b,ED=trajED,mn=trajm,alloc=trajpi,stp=stp)) } ############## some auxiliary functions ################# modex <- function(x) { # returns mode (lowest if multiple solutions) tbl <- table(x) i <- which.max(tbl) return(as.numeric(names(tbl)[i])) } ed95 <- function(x) { # computes ED95 among a list of mean responses fstar <- min(x)+0.95*(max(x)-min(x)) k <- which.min(abs(x-fstar)) return(k) } plt.st <- function(st,plt=1) { # plots summary of *one* simulated trial attach(st) if (plt==1){ ## plot truth, estimated curves & data yl <- c(-1,6) plot(doses,f(doses,b),type="l",lty=3,xlab="DOSES", ylab="", bty="l",ylim=yl) matlines(doses,t(mn),col=1,lty=2) mm <- nrow(mn) lines(doses,mn[mm,],lty=1,lwd=3) points(Y[,1],Y[,2],pch=19,cex=.5) abline(v=ED[mm],col=2) } if (plt==2) ## plot allocations plot(Y[,1],xlab="PATIENT",ylab="DOSE",bty="l",pch=19,ylim=c(1,16)) detach(st) } ############### FFBS ##################################### ffbs <- function(n,p,Y,G,Z,sig,tau,m0=rep(0,p),C0=diag(p), y0=rep(0,n), n0=rep(0.25,n), M=0,plt=F,lw=1) {## implements basic FFBS ## Y = (N x 2) matrix of dose (1st col) & response (2nd column) ## n = number doses ## x = (n x p) matrix of (latent) states ## OBS EQ: y = Z x + eps eps ~ N(0,sig) ## EVOL EQ: x[t+1] = G x[t] + eta, eta ~ N(0,tau*I) ## (m0,C0): prior moments of x[0] ## (y0[j],n0[j]): prior by n0[j] fake data point points with mean y0[j] ## ## use plt and lw for debugging ## plausibility check of paramters if (nrow(G) != p){ cat("\n *** Error: dim(G) != p. p=",p,", dim(G)=",nrow(G),"\n") return(-1) } if (!is.null(Y)) if (max(Y[,1]) > n){ cat("\n *** Error: max dose in Y > n; mx=",max(Y[,1]),", n=",n,"\n") return(-1) } if ( (ncol(Z) != p) | (nrow(Z) !=1)){ cat("\n *** Error: dim of Z do not match p:",ncol(Z),nrow(Z),"\n") return(-1) } ## arrange the data by dose: record mean yj and n obs nj per dose yj <- rep(0,n) nj <- rep(0,n) for(i in 1:n){ idx <- which(Y[,1]==i) nji <- length(idx) if (nji==0) yji <- 0 else yji <- mean(Y[idx,2]) nj[i] <- nji+n0[i] yj[i] <- (nji*yji+n0[i]*y0[i])/(nj[i]) } forw <- ff(yj,nj,G,Z,sig,tau,C0,m0,plt=F) backw <- bs(forw$C,forw$m,G,tau,M,plt,lw=lw) return(list(m=backw$mn[1,],s=sqrt( backw$Cn[1,1,] ),X=backw$X)) } ff <- function(yj,nj,G,Z,sig,tau,C0,m0,plt=F,lw=1) { ## forward filtering ## let Y[t]=(y1..yt) ## plausibility check of paramters ## use plt and lw for debugging n <- length(yj) p <- nrow(G) doses <- 1:n if ( (length(m0) != p) | (nrow(C0) !=p)) { cat("\n *** Error: length(m0)=",length(m0)," != p=",p,"\n") return(-1) } if (nrow(C0) != p){ cat("\n *** Error: nrow(C0)=",nrow(C0)," != p=",p,"\n") return(-1) } m <- matrix(0,ncol=n,nrow=p) # E(x[t] | Y[t]), C <- array(0,dim=c(p,p,n)) # V(x[t] | Y[t]) Hi <- t(Z)%*%Z/sig # likelihood precision Om <- tau*diag(p) ## FF: Forward filtering ## R = G C G' + tau*I' Clag <- C0 mlag <- m0 for(i in 1:n){ R <- G %*% Clag %*% t(G) + Om # R = GCG' + Omega r <- G %*% mlag # r = Gm if (nj[i]==0){ C[,,i] <- R m[ ,i] <- r } else { Ri <- solve(R) Ci <- Ri + nj[i]*Hi C[,,i] <-solve(Ci) m[ ,i] <- C[,,i] %*% (Ri%*%r + nj[i]/sig* t(Z) %*% yj[i]) } Clag <- C[,,i] # prepare for next iteration mlag <- m[ ,i] } if (plt) lines(doses,m[1,],col=3,lty=2,lwd=lw) return(list(m=m,C=C)) } bs <- function(C,m,G,tau,M=0,plt=F,lw=1) { ## BS: backward smoothing ## use plt and lw for debugging p <- nrow(m) n <- ncol(m) doses <- 1:n mn <- matrix(0,nrow=p,ncol=n) # E(x[t] | Y[n]) Cn <- array(0,dim=c(p,p,n)) # V(x[t] | Y[n]) if (M>0) X <- array(0,dim=c(p,n,M)) else X <- 0 # dummy for use in final return statement ## p(xn | Yn) = N(m[n],C[n]) ## p(xi | Yn, x[i+1..n]) = c* ## p(xi | Yi, x[i+1..n]) * p(y[i+1..n] | Yi,xi,x[i+1..n]) = ## p(xi | Yi) * p(x[i+1] | x[i]) = ## N(m[i],C[i])*N(x[i+1]; G x[i], Om) ## => 1/Cn[i] = 1/C[i] + G'G/tau ## mn[i] = Cn[i]*( 1/C[i]*m[i] + G' x[i+1] ) mn[,n] <- m[,n] Cn[,,n] <- C[,,n] if (M>0) # generate states L <- t( chol(Cn[,,n]) ) for(j in 1:M) X[,n,j] <- L %*% rnorm(2) + mn[,n] for(i in (n-1):1){ Ci <- solve(C[,,i]) Cni <- Ci + t(G)%*%G/tau Cn[,,i] <- solve(Cni) mn[ ,i] <- Cn[,,i] %*% (Ci %*% m[,i]+1/tau*t(G) %*% mn[,i+1]) if (M>0){ # generate states L <- t( chol(Cn[,,i]) ) for(j in 1:M){ mu <- Cn[,,i] %*% (Ci %*% m[,i]+1/tau*t(G) %*% X[,i+1,j]) X[,i,j] <- L %*% rnorm(2) + mu }# for j }#M>0 }# for i if (plt) lines(doses,mn[1,],col=2,lty=2,lwd=lw) return(list(mn=mn,Cn=Cn,X=X)) } ############### run examples ############################# examples <- function() { ## execute these lines to run the example st <- sim.trial(b=scenario4,plt=T) ps("35-4m") plt.st(st,1) devoff() ps("35-4d") plt.st(st,2) devoff() st <- sim.trial(b=scenario1,plt=T) ps("35-1m") plt.st(st,1) devoff() ps("35-1d") plt.st(st,2) devoff() st <- sim.trial(b=scenario0,plt=T) ps("35-0m") plt.st(st,1) devoff() ps("35-0d") plt.st(st,2) devoff() }