logit <- function(p) return(log(p/(1-p))) expit <- function(x) exp(x)/(1+exp(x)) prob.calc <- function(alpha, beta, d) { prob.matrix <- NULL for (i in 1:length(d)) { numerator <- exp(alpha+beta*d[i]) denominator <- 1+numerator p <- numerator/denominator prob.matrix <- cbind(prob.matrix, apply(p, 2, mean, na.rm=T)) } return(prob.matrix) } quant.calc <- function(alpha, beta, d, gamma=0.05) { q.matrix <- NULL for (i in 1:length(d)) { numerator <- exp(alpha+beta*d[i]) denominator <- 1+numerator p <- numerator/denominator q.matrix <- cbind(q.matrix, apply(p, 2, quantile, probs=gamma, na.rm=T)) } return(q.matrix) } get.posterior <- function(Ndraw, mu.alpha, mu.beta, mu.delta, s.alpha, s.beta, y, d, A, B) { data <- list(mu.beta=mu.beta, mu.alpha=mu.alpha, mu.delta=mu.delta, s.alpha=1/s.alpha, s.beta=s.beta, d=d, A=A, B=B, y=y, n=length(y)) myjags <- jags.model(paste(dir, "/my model.txt", sep=""), data) mysamp <- as.matrix(coda.samples(myjags, c("alpha", "beta"), Ndraw)) mysamp[-(1:burnin),] } pick <- function(nsubj, target, Dlt.true, nA, nB, Ndraw, mu.alpha, mu.beta, mu.delta, s.alpha, s.beta) { N.Matrix <- Dlt.Matrix <- matrix(0, nB, nA) posB <- posA <- 1 N.Matrix[1,1] <- 1 y <- Dlt.Matrix[1,1] <- rbinom(1,1,Dlt.true[1,1]) dlt <- NULL B <- A <- 1 for (i in 1:nsubj) { #Record true DLT rate of each assignment dlt <- c(dlt, Dlt.true[posB,posA]) #Compute posterior estimates of DLT probabilties temp <- get.posterior(Ndraw, mu.alpha, mu.beta, mu.delta, s.alpha, s.beta, y, d, A, B) alpha <- temp[,-ncol(temp)] beta <- temp[,ncol(temp)] prob.matrix <- prob.calc(alpha, beta, d) q.matrix <- quant.calc(alpha, beta, d, target) #Compute prop of posterior distribution for DLT rate of combo (1,1) #and check if trial should be stopped num <- exp(alpha[,1]+beta*d[1]) den <- 1+num p11post <- num/den p11q <- mean(p11post>target) if (p11q>0.9) break # LCL <- binom.test(sum(Dlt.Matrix), sum(N.Matrix))$conf.int[1] # if (i>=3 && LCL>target) break #Determine current estimate of MTC in neighborhood of previous MTC nbrhdB <- posB + -1:1; nbrhdB <- nbrhdB[nbrhdB>0 & nbrhdB<=nB] nbrhdA <- posA + -1:1; nbrhdA <- nbrhdA[nbrhdA>0 & nbrhdA<=nA] dif <- NULL for (k in nbrhdB) for (l in nbrhdA) dif <- rbind(dif, c(l, k, abs(prob.matrix[k,l]-target))) posA <- dif[dif[,3]==min(dif[,3]),1] posB <- dif[dif[,3]==min(dif[,3]),2] N.Matrix[posB,posA] <- N.Matrix[posB,posA]+1 temp.y <- rbinom(1, 1, Dlt.true[posB,posA]) y <- c(y,temp.y) B <- c(B,posB) A <- c(A,posA) Dlt.Matrix[posB,posA] <- Dlt.Matrix[posB,posA]+temp.y } if(i==nsubj) { N.Matrix[posB,posA] <- N.Matrix[posB,posA]-1 dlt <- dlt[1:nsubj] } if (i < nsubj) posB <- posA <- NULL list(Dlt.estimate=prob.matrix, Dlt.true=Dlt.true, N.Matrix=N.Matrix, Dlt.Matrix=Dlt.Matrix, choose=c(posB,posA), dlt=dlt, y=y, A=A, B=B) }