Nothing
R/Trial.Simulation.R
In phase1RMD: Repeated Measurement Design for Phase I Clinical Trial
###################################################################################
###########Simulate Data from the Matrix and Do Stage 1-3##########################
###################################################################################
###model: yij = beta0 + beta1xi + beta2tj + ui + epiij#############################
######### ei = alpha0 + alpha1xi + alpha2xi^2 + gamma0ui + epi2i###################
######### epi ~ N(0, sigma_e); ui ~ N(0, sigma_u); epi2 ~ N(0, sigma_f)############
###################################################################################
#library(rjags)
#library(R2WinBUGS)
#library(BiocParallel)
# SimRMDEff <- function(seed=2441, tox_matrix, strDose = 1, chSize = 3, trlSize = 36, sdose = 1:6, MaxCycle = 6,
# tox.target = 0.28, eff.structure = c(0.1, 0.2, 0.3, 0.4, 0.7, 0.9)){
# # SimRMDEff <- function(seed=2441, tox_matrix, StrDose = 1, chSize = 3, trlSize = 36, sdose = 1:6, MaxCycle = 6,
# # tox.target = 0.28, eff.structure = c(0.1, 0.2, 0.3, 0.4, 0.7, 0.9),
# # eff.sd = 0.2, p1 = 0.2, p2 = 0.2,
# # ps1 = 0.2,
# # proxy.thrd = 0.1, thrd1 = 0.28, thrd2 = 0.28,
# # wm = matrix(c(0, 0.5, 0.75, 1 , 1.5,
# # 0, 0.5, 0.75, 1 , 1.5,
# # 0, 0 , 0 , 0.5, 1 ),
# # byrow = T, ncol = 5),
# # toxmax = 2.5, toxtype = NULL, intercept.alpha = NULL,
# # coef.beta = NULL, cycle.gamma = NULL){
# #revision031917
# #dummy variable defination for winbugs parameters
# beta_other <- 0;
# beta_dose <- 0;
# N1 <- 0;
# inprod <- 0;
# u <- 0;
# N2 <- 0;
# Nsub <- 0;
# alpha <- 0;
# gamma0 <- 0;
# #revision031917
# n.iter <- 3000;
# eff.sd <- 0.2; p1 <- 0.2; p2 <- 0.2;
# ps1 <- 0.2; proxy.thrd <- 0.1; thrd1 <- 0.28; thrd2 <- 0.28;
# wm <- matrix(c(0, 0.5, 0.75, 1 , 1.5,
# 0, 0.5, 0.75, 1 , 1.5,
# 0, 0 , 0 , 0.5, 1 ),
# byrow = T, ncol = 5);
# toxmax <- 2.5; toxtype = NULL; intercept.alpha = NULL;
# coef.beta <- NULL; cycle.gamma <- NULL;
# set.seed(seed)
# trialSize <- trlSize;
# doses <- sdose;
# cycles <- 1:MaxCycle;
# listdo <- paste0("D", doses)
# if(length(eff.structure) != length(doses))
# stop("Check if you have specified the efficacy mean structure for the correct number of doses.")
# if(nrow(wm) != 3)
# stop("Make sure that you specified the clinical weight matrix for three types of toxicities.")
# MaxCycle <- length(cycles)
# Numdose <- length(doses)
# flag.tox.matrix.null <- 0
# if(is.null(tox_matrix)){
# flag.tox.matrix.null <- 1
# tox_matrix <- array(NA, dim = c(Numdose, MaxCycle, 3, 5))
# if(length(toxtype) != 3)
# stop("Right now we are only considering three toxicity types, but we will relax this constraint in the future work.")
# if(length(intercept.alpha) != 4){
# stop("Exactly four intercepts alpha are needed for grade 0--4 in simulation!")
# }
# if(min(diff(intercept.alpha, lag = 1, differences = 1)) < 0){
# stop("Intercepts alpha for simulation must be in a monotonic increasing order!")
# }
# if(length(coef.beta) != 3){
# stop("Exactly three betas need to be specified for three types of toxicities!")
# }
# }
# #############################################
# #########Determin Scenarios##################
# #############################################
# w1 <- wm[1, ]
# w2 <- wm[2, ]
# w3 <- wm[3, ]
# ####Array of the normalized scores########
# nTTP.array <- function(w1,w2,w3){
# nTTP <- array(NA,c(5,5,5))
# for (t1 in 1:5){
# for (t2 in 1:5){
# for (t3 in 1:5){
# nTTP[t1,t2,t3] <-
# sqrt(w1[t1]^2 + w2[t2]^2 + w3[t3]^2)/toxmax
# }
# }
# }
# return(nTTP)
# }
# nTTP.all <- nTTP.array(w1,w2,w3)
# ####compute pdlt given probability array####
# pDLT <- function(proba){
# return(sum(proba[c(4,5), , ])
# + sum(proba[ ,c(4,5), ])
# + sum(proba[ , ,5])
# - sum(proba[c(4,5),,])*sum(proba[,c(4,5),])
# - sum(proba[c(4,5),,])*sum(proba[,,5])
# - sum(proba[,c(4,5),])*sum(proba[,,5])
# + sum(proba[c(4,5),,])*sum(proba[,c(4,5),])*sum(proba[,,5]))
# }
# ####compute mean nTTP given probability array######
# mnTTP <-function(proba){
# return(sum(proba * nTTP.all))
# }
# sc.mat <- matrix(NA, MaxCycle, Numdose)
# pdlt <- matrix(NA, MaxCycle, Numdose)
# for(i in 1:MaxCycle){ ##loop through cycle
# for(j in 1:Numdose){ ##loop through dose
# if(flag.tox.matrix.null == 0){
# nTTP.prob <- tox_matrix[j, i, ,]
# }else{
# # initialize storing matrix
# cump <- matrix(NA, nrow = length(toxtype), ncol = 5) # cumulative probability
# celp <- matrix(NA, nrow = length(toxtype), ncol = 5) # cell probability
# for(k in 1:length(toxtype)){
# # proportional odds model
# logitcp <- intercept.alpha + coef.beta[k] * j + cycle.gamma * (i - 1)
# # cumulative probabilities
# cump[k, 1:4] <- exp(logitcp) / (1 + exp(logitcp))
# cump[k, 5] <- 1
# # cell probabilities
# celp[k, ] <- c(cump[k,1], diff(cump[k, ], lag = 1, differences = 1))
# }
# nTTP.prob <- celp
# tox_matrix[j, i, ,] <- celp
# }
# proba <- outer(outer(nTTP.prob[1, ], nTTP.prob[2, ]), nTTP.prob[3, ])
# sc.mat[i, j] <- round(mnTTP(proba), 4)
# pdlt[i, j] <- round(pDLT(proba), 4)
# }
# }
# mEFF <- eff.structure
# sc <- rbind(sc.mat, mEFF, pdlt[1, ])
# colnames(sc) <- listdo
# rownames(sc) <- c("mnTTP.1st", "mnTTP.2nd", "mnTTP.3rd",
# "mnTTP.4th", "mnTTP.5th", "mnTTP.6th",
# "mEFF", "pDLT")
# ####################################################################################
# k <- 1 ####Let us use k to denote the column of altMat####
# doseA <- strDose ####Let us use doseA to denote the current dose assigned
# cmPat <- chSize ####current patient size
# masterData <- list()
# rec_doseA <- NULL
# stage1.allow <- NULL
# altMat <- matrix(0, Numdose, k)
# altMat[doseA, k] <- chSize
# ###############################################################
# ########################Stage 1################################
# ###############################################################
# #print('Stage 1')
# while(cmPat <= trialSize/2){
# rec_doseA <- c(rec_doseA, doseA)
# PatData <- NULL
# n <- chSize
# K <- MaxCycle
# ##############################################################
# #############design matrix for toxicity#######################
# ##############################################################
# dose.pat <- rep(doseA, n)
# xx <- expand.grid(dose.pat, 1:K)
# X <- as.matrix(cbind(rep(1, n * K), xx))
# X <- cbind(X, ifelse(X[, 3] == 1, 0, 1))
# colnames(X) <- c("intercept", "dose", "cycle", "icycle")
# ##############################################################
# #######design matrix for efficacy#############################
# ##############################################################
# Xe <- cbind(1, rep(doseA, n), rep(doseA^2, n))
# colnames(Xe) <- c("intercept", "dose", "squared dose")
# #####simulate nTTP and Efficacy for the cohort########
# outcome <- apply(X, 1, function(a){
# tox.by.grade <- tox_matrix[a[2], a[3], ,]
# nttp.indices <- apply(tox.by.grade, 1, function(o){
# sample(1:5, 1, prob = o)
# })
# if(max(wm[cbind(1:nrow(wm), nttp.indices)]) >= 1){
# y.dlt <- 1
# }else{
# y.dlt <- 0
# }
# y.nTTP <- nTTP.all[nttp.indices[1],
# nttp.indices[2],
# nttp.indices[3]]
# return(c(y.dlt = y.dlt,
# y = y.nTTP))
# })
# y <- outcome[2, ]
# y.dlt <- outcome[1, ]
# beta.mean <- eff.structure[Xe[1, 2]]
# beta.sd <- eff.sd
# beta.shape1 <- ((1 - beta.mean)/beta.sd^2 - 1/beta.mean) * beta.mean^2
# beta.shape2 <- beta.shape1 * (1/beta.mean - 1)
# e <- matrix(rbeta(chSize, beta.shape1, beta.shape2), ncol = 1)
# #####construct master Dataset and extract from this dataset for each interim#########
# #####PatData stores the dataset used to estimate the model at each interim###########
# temp.mtdata <- data.frame(cbind(y.dlt, y, X, rep(cmPat/chSize, n * K), 1:n))
# temp.mtdata$subID <- paste0("cohort", temp.mtdata[, 7], "subject", temp.mtdata[, 8])
# masterData$toxicity <- data.frame(rbind(masterData$toxicity, temp.mtdata))
# temp.eff.mtdata <- data.frame(cbind(e, Xe, cmPat/chSize))
# masterData$efficacy <- data.frame(rbind(masterData$efficacy, temp.eff.mtdata))
# current_cohort <- cmPat/chSize
# PatData.index <- data.frame(cohort = 1:current_cohort, cycles = current_cohort:1)
# PatData.index <- within(PatData.index, cycles <- ifelse(cycles >= MaxCycle, MaxCycle, cycles))
# for(i in 1:nrow(PatData.index)){
# PatData <- rbind(PatData, masterData$toxicity[masterData$toxicity[, 7] == PatData.index[i, "cohort"] &
# masterData$toxicity[, 5] <= PatData.index[i, "cycles"],
# c(1, 2, 3, 4, 5, 6, 9)])
# }
# #####dropout.dlt is 0 means to keep the observation; 1 means to drop it due to dlt#######
# dlt.drop <- sapply(by(PatData, PatData$subID, function(a){a[, 1]}), function(o){
# if((length(o) > 1 & all(o[1: (length(o) - 1)] == 0)) | length(o) == 1){
# return(rep(0, length(o)))
# }else{
# o.temp <- rep(0, length(o))
# o.temp[(min(which(o == 1)) + 1) : length(o.temp)] <- 1
# return(o.temp)
# }
# })
# dropout.dlt <- NULL
# for(i in unique(PatData$subID)){
# dropout.dlt[PatData$subID == i] <- dlt.drop[[i]]
# }
# PatData <- PatData[dropout.dlt == 0, ]
# ###################################################################################
# ################Model Estimation given PatData#####################################
# ################Stage 1 only considers the toxicity model##########################
# ###################################################################################
# nTTP <- PatData[, 2]
# X_y <- as.matrix(PatData[, c("intercept", "dose", "icycle")])
# n.subj <- length(unique(PatData[, "subID"]))
# W_y <- matrix(0, nrow(PatData), n.subj)
# W_y[cbind(1:nrow(W_y), as.numeric(sapply(PatData$subID, function(a){
# which(a == unique(PatData$subID))
# })))] <- 1
# model.file <- function()
# {
# beta <- c(beta_other[1], beta_dose, beta_other[2])
# for(i in 1:N1){
# y[i] ~ dnorm(mu[i], tau_e)
# mu[i] <- inprod(X_y[i, ], beta) + inprod(W_y[i, ], u)
# }
# for(j in 1:N2){
# u[j] ~ dnorm(0, tau_u)
# }
# beta_other ~ dmnorm(p1_beta_other[], p2_beta_other[, ])
# beta_dose ~ dunif(p1_beta_dose, p2_beta_dose)
# tau_e ~ dunif(p1_tau_e, p2_tau_e)
# tau_u ~ dunif(p1_tau_u, p2_tau_u)
# }
# mydata <- list(N1 = length(nTTP), N2 = n.subj, y = nTTP, X_y = X_y,
# W_y = W_y, p1_beta_other = c(0, 0),
# p2_beta_other = diag(rep(0.001, 2)),
# p1_beta_dose = 0, p2_beta_dose = 1000,
# p1_tau_e = 0, p2_tau_e = 1000,
# p1_tau_u = 0, p2_tau_u = 1000)
# path.model <- file.path(tempdir(), "model.file.txt")
# R2WinBUGS::write.model(model.file, path.model)
# inits.list <- list(list(beta_other = c(0.1, 0.1),
# beta_dose = 0.1,
# tau_e = 0.1,
# tau_u = 0.1,
# .RNG.seed = sample(1:1e+06, size = 1),
# .RNG.name = "base::Wichmann-Hill"))
# jagsobj <- rjags::jags.model(path.model, data = mydata, n.chains = 1,
# quiet = TRUE, inits = inits.list)
# update(jagsobj, n.iter = n.iter, progress.bar = "none")
# post.samples <- rjags::jags.samples(jagsobj, c("beta_dose", "beta_other"),
# n.iter = n.iter,
# progress.bar = "none")
# if(cmPat == trialSize/2){
# ######################################################################
# #############define allowable doses###################################
# ######################################################################
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= p1 & prob2.doses >= p2)
# if(length(allow.doses) == 0){
# stage1.allow <- 0
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# stage1.allow <- allow.doses
# }else{
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1]))
# loss.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# abs(o[1] + o[2] * a - tox.target)
# }))
# })
# nxtdose <- doses[which.min(loss.doses)]
# if(as.numeric(nxtdose) > (doseA + 1)){
# doseA <- doseA + 1;
# }else{
# doseA <- as.numeric(nxtdose);
# }
# if (cmPat == chSize){
# dlt <- with(masterData$toxicity, y.dlt[cycle == 1 & V7 == cmPat/chSize])
# if (sum(dlt) == 0){
# doseA <- 2
# dose_flag <- 0
# }else{
# doseA <- 1
# dose_flag <- 1
# }
# }
# if ((cmPat >= chSize*2) & (dose_flag == 1)){
# dlt <- with(masterData$toxicity, y.dlt[cycle == 1 & V7 == cmPat/chSize])
# if (sum(dlt) == 0){
# doseA <- 2
# dose_flag <- 0
# }else{
# doseA <- 1
# dose_flag <- 1
# }
# }
# if(cmPat >= chSize*3){
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= ps1 & prob2.doses >= ps1)
# if(length(allow.doses) == 0){
# stage1.allow <- 0
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# }
# altMat[doseA, k] <- altMat[doseA, k] + chSize
# }
# cmPat <- cmPat + chSize # cmPat for the next cohort (loop)
# }
# #print('Stage 2')
# ###############################################################
# ########################Stage 2################################
# ###############################################################
# #######let us wait until the efficacy data for stage 1 is available######
# #########################################################################
# cmPat <- cmPat - chSize
# PatData <- list()
# PatData.index <- data.frame(cohort = current_cohort:1, cycles = 3:(3 + current_cohort - 1))
# PatData.index <- within(PatData.index, cycles <- ifelse(cycles >= MaxCycle, MaxCycle, cycles))
# for(i in 1:nrow(PatData.index)){
# PatData$toxicity <- rbind(PatData$toxicity, masterData$toxicity[masterData$toxicity[, 7] == PatData.index[i, "cohort"] &
# masterData$toxicity[, 5] <= PatData.index[i, "cycles"],
# c(1, 2, 3, 4, 5, 6, 9)])
# }
# #####change the ordering of the subjects--being consistent######
# temp.toxicity <- NULL
# for(i in 1:current_cohort){
# temp.toxicity <- rbind(temp.toxicity, PatData$toxicity[grepl(paste0("cohort",i,"subject"), PatData$toxicity$subID), ])
# }
# PatData$toxicity <- temp.toxicity
# rm(temp.toxicity)
# #####toxicity dropout########
# #####dropout.dlt is 0 means to keep the observation; 1 means to drop it due to dlt
# dlt.drop <- sapply(by(PatData$toxicity, PatData$toxicity$subID, function(a){a[, 1]}), function(o){
# if((length(o) > 1 & all(o[1: (length(o) - 1)] == 0)) | length(o) == 1){
# return(rep(0, length(o)))
# }else{
# o.temp <- rep(0, length(o))
# o.temp[(min(which(o == 1)) + 1) : length(o.temp)] <- 1
# return(o.temp)
# }
# })
# dropout.dlt <- NULL
# for(i in unique(PatData$toxicity$subID)){
# dropout.dlt[PatData$toxicity$subID == i] <- dlt.drop[[i]]
# }
# PatData$toxicity <- PatData$toxicity[dropout.dlt == 0, ]
# #####efficacy dropout########
# dropout.eff <- sapply(dlt.drop, function(a){
# if(sum(a) == 0){
# return(0)
# }else if(min(which(a == 1)) >= 4){
# return(0)
# }else{
# return(1)
# }
# })
# dropout.eff <- as.numeric(dropout.eff)
# PatData$efficacy <- masterData$efficacy[dropout.eff == 0, ]
# ############Model Estimation to randomize among allowable doses##########
# nTTP <- PatData$toxicity[, 2]
# X_y <- as.matrix(PatData$toxicity[, c("intercept", "dose", "icycle")])
# n.subj <- length(unique(PatData$toxicity[, "subID"]))
# W_y <- matrix(0, nrow(PatData$toxicity), n.subj)
# W_y[cbind(1:nrow(W_y), as.numeric(sapply(PatData$toxicity$subID, function(a){
# which(a == unique(PatData$toxicity$subID))
# })))] <- 1
# EFF <- PatData$efficacy[, 1]
# X_e <- as.matrix(PatData$efficacy[, c("intercept", "dose", "squared.dose")])
# keepeff.ind <- which(dropout.eff == 0)
# model.file <- function()
# {
# beta <- c(beta_other[1], beta_dose, beta_other[2])
# for(i in 1:N1){
# y[i] ~ dnorm(mu[i], tau_e)
# mu[i] <- inprod(X_y[i, ], beta) + inprod(W_y[i, ], u)
# }
# for(k in 1:Nsub){
# u[k] ~ dnorm(0, tau_u)
# }
# for(j in 1:N2){
# e[j] ~ dnorm(mu_e[j], tau_f)
# mu_e[j] <- inprod(X_e[j, ], alpha) + gamma0 * u[keepeff.ind[j]]
# }
# beta_other ~ dmnorm(p1_beta_other[], p2_beta_other[, ])
# beta_dose ~ dunif(p1_beta_dose, p2_beta_dose)
# alpha ~ dmnorm(p1_alpha[], p2_alpha[, ])
# gamma0 ~ dnorm(p1_gamma0, p2_gamma0)
# tau_e ~ dunif(p1_tau_e, p2_tau_e)
# tau_u ~ dunif(p1_tau_u, p2_tau_u)
# tau_f ~ dunif(p1_tau_f, p2_tau_f)
# }
# mydata <- list(N1 = length(nTTP), N2 = length(EFF), y = nTTP,
# Nsub = n.subj, X_y = X_y, e = EFF, keepeff.ind = keepeff.ind,
# W_y = W_y, X_e = X_e, p1_beta_other = c(0, 0),
# p2_beta_other = diag(rep(0.001, 2)),
# p1_beta_dose = 0, p2_beta_dose = 1000,
# p1_alpha = c(0, 0, 0), p2_alpha = diag(rep(0.001, 3)),
# p1_gamma0 = 0, p2_gamma0 = 0.001,
# p1_tau_e = 0, p2_tau_e = 1000,
# p1_tau_u = 0, p2_tau_u = 1000,
# p1_tau_f = 0, p2_tau_f = 1000)
# path.model <- file.path(tempdir(), "model.file.txt")
# R2WinBUGS::write.model(model.file, path.model)
# inits.list <- list(list(beta_other = c(0.1, 0.1),
# beta_dose = 0.1,
# alpha = c(0.1, 0.1, 0.1),
# gamma0 = 0.1,
# tau_e = 0.1,
# tau_u = 0.1,
# tau_f = 0.1,
# .RNG.seed = sample(1:1e+06, size = 1),
# .RNG.name = "base::Wichmann-Hill"))
# jagsobj <- rjags::jags.model(path.model, data = mydata, n.chains = 1,
# quiet = TRUE, inits = inits.list)
# update(jagsobj, n.iter = n.iter, progress.bar = "none")
# post.samples <- rjags::jags.samples(jagsobj, c("beta_dose", "beta_other", "alpha",
# "gamma0"),
# n.iter = n.iter,
# progress.bar = "none")
# ######################################################################
# #############update allowable doses###################################
# ######################################################################
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= p1 & prob2.doses >= p2)
# if(length(allow.doses) == 0){
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# sim.alphas <- as.matrix(rbind(post.samples$alpha[,,1]))
# RAND.EFF <- sapply(allow.doses, function(a){
# mean(apply(sim.alphas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * a^2)
# }))
# })
# RAND.EFF <- exp(RAND.EFF)/sum(exp(RAND.EFF))
# nxtdose <- sample(allow.doses, 1, prob = RAND.EFF)
# ####a condition for untried higher dose level that is randomized#### rec_doseA[length(rec_doseA)]
# if(nxtdose > max(rec_doseA) + 1 & !nxtdose %in% rec_doseA){
# nxtdose <- max(rec_doseA) + 1
# }
# #################################################################
# ########generate data for the new enrolled cohort in Stage 2#####
# #################################################################
# while(cmPat < trialSize - chSize){
# doseA <- nxtdose
# rec_doseA <- c(rec_doseA, doseA)
# altMat[doseA, k] <- altMat[doseA, k] + chSize
# cmPat <- cmPat + chSize
# PatData <- list()
# n <- chSize
# K <- MaxCycle
# #######################################################
# #################design matrix for toxicity############
# #######################################################
# dose.pat <- rep(doseA, n)
# xx <- expand.grid(dose.pat, 1:K)
# X <- as.matrix(cbind(rep(1, n*K), xx))
# X <- cbind(X, ifelse(X[, 3] == 1, 0, 1))
# colnames(X) <- c("intercept", "dose", "cycle", "icycle")
# #########################################################
# #######design matrix for efficacy########################
# #########################################################
# Xe <- cbind(1, rep(doseA, n), rep(doseA^2, n))
# colnames(Xe) <- c("intercept", "dose", "squared dose")
# #####simulate nTTP and Efficacy for the cohort########
# outcome <- apply(X, 1, function(a){
# tox.by.grade <- tox_matrix[a[2], a[3], ,]
# nttp.indices <- apply(tox.by.grade, 1, function(o){
# sample(1:5, 1, prob = o)
# })
# if(max(wm[cbind(1:nrow(wm), nttp.indices)]) >= 1){
# y.dlt <- 1
# }else{
# y.dlt <- 0
# }
# y.nTTP <- nTTP.all[nttp.indices[1],
# nttp.indices[2],
# nttp.indices[3]]
# return(c(y.dlt = y.dlt,
# y = y.nTTP))
# })
# y <- outcome[2, ]
# y.dlt <- outcome[1, ]
# beta.mean <- eff.structure[Xe[1, 2]]
# beta.sd <- eff.sd
# beta.shape1 <- ((1 - beta.mean)/beta.sd^2 - 1/beta.mean) * beta.mean^2
# beta.shape2 <- beta.shape1 * (1/beta.mean - 1)
# e <- matrix(rbeta(chSize, beta.shape1, beta.shape2), ncol = 1)
# #################################################################
# #################################################################
# #####add to the master Dataset and extract from this dataset for each interim#########
# #####PatData stores the dataset used to estimate the model at each interim############
# temp.mtdata <- data.frame(cbind(y.dlt, y, X, rep(cmPat/chSize, n*K), 1:n))
# temp.mtdata$subID <- paste0("cohort", temp.mtdata[, 7], "subject", temp.mtdata[, 8])
# masterData$toxicity <- data.frame(rbind(masterData$toxicity, temp.mtdata))
# temp.eff.mtdata <- data.frame(cbind(e, Xe, cmPat/chSize))
# masterData$efficacy <- data.frame(rbind(masterData$efficacy, temp.eff.mtdata))
# current_cohort <- cmPat/chSize
# PatData.index <- data.frame(cohort = 1:current_cohort, cycles = current_cohort:1)
# cycles.adj <- c(rep(2, trialSize/(2 * chSize)), rep(0, current_cohort - trialSize/(2 * chSize)))
# PatData.index <- within(PatData.index, {cycles <- cycles + cycles.adj
# cycles <- ifelse(cycles >= MaxCycle, MaxCycle, cycles)})
# for(i in 1:nrow(PatData.index)){
# PatData$toxicity <- rbind(PatData$toxicity, masterData$toxicity[masterData$toxicity[, 7] == PatData.index[i, "cohort"] &
# masterData$toxicity[, 5] <= PatData.index[i, "cycles"],
# c(1, 2, 3, 4, 5, 6, 9)])
# }
# #####toxicity dropout########
# #####dropout.dlt is 0 means to keep the observation; 1 means to drop it due to dlt
# dlt.drop <- sapply(by(PatData$toxicity, PatData$toxicity$subID, function(a){a[, 1]}), function(o){
# if((length(o) > 1 & all(o[1: (length(o) - 1)] == 0)) | length(o) == 1){
# return(rep(0, length(o)))
# }else{
# o.temp <- rep(0, length(o))
# o.temp[(min(which(o == 1)) + 1) : length(o.temp)] <- 1
# return(o.temp)
# }
# })
# dropout.dlt <- NULL
# for(i in unique(PatData$toxicity$subID)){
# dropout.dlt[PatData$toxicity$subID == i] <- dlt.drop[[i]]
# }
# PatData$toxicity <- PatData$toxicity[dropout.dlt == 0, ]
# cohort.eff.index <- with(PatData.index, which(cycles >= 3))
# #####efficacy dropout########
# dropout.eff <- sapply(dlt.drop, function(a){
# if(sum(a) == 0){
# return(0)
# }else if(min(which(a == 1)) >= 4){
# return(0)
# }else{
# return(1)
# }
# })
# dropout.eff <- as.numeric(dropout.eff)
# PatData$efficacy <- subset(masterData$efficacy, masterData$efficacy[, 5] %in% cohort.eff.index & dropout.eff == 0)
# #############################################################################################################################
# ######################Model Estimation to update toxicity information and randomization probs################################
# #############################################################################################################################
# nTTP <- PatData$toxicity[, 2]
# X_y <- as.matrix(PatData$toxicity[, c("intercept", "dose", "icycle")])
# n.subj <- length(unique(PatData$toxicity[, "subID"]))
# W_y <- matrix(0, nrow(PatData$toxicity), n.subj)
# W_y[cbind(1:nrow(W_y), as.numeric(sapply(PatData$toxicity$subID, function(a){
# which(a == unique(PatData$toxicity$subID))
# })))] <- 1
# EFF <- PatData$efficacy[, 1]
# X_e <- as.matrix(PatData$efficacy[, c("intercept", "dose", "squared.dose")])
# keepeff.ind <- which(dropout.eff == 0 & masterData$efficacy[, 5] %in% cohort.eff.index)
# model.file <- function()
# {
# beta <- c(beta_other[1], beta_dose, beta_other[2])
# for(k in 1:Nsub){
# u[k] ~ dnorm(0, tau_u)
# }
# for(i in 1:N1){
# y[i] ~ dnorm(mu[i], tau_e)
# mu[i] <- inprod(X_y[i, ], beta) + inprod(W_y[i, ], u)
# }
# for(j in 1:N2){
# e[j] ~ dnorm(mu_e[j], tau_f)
# mu_e[j] <- inprod(X_e[j, ], alpha) + gamma0 * u[keepeff.ind[j]]
# }
# beta_other ~ dmnorm(p1_beta_other[], p2_beta_other[, ])
# beta_dose ~ dunif(p1_beta_dose, p2_beta_dose)
# alpha ~ dmnorm(p1_alpha[], p2_alpha[, ])
# gamma0 ~ dnorm(p1_gamma0, p2_gamma0)
# tau_e ~ dunif(p1_tau_e, p2_tau_e)
# tau_u ~ dunif(p1_tau_u, p2_tau_u)
# tau_f ~ dunif(p1_tau_f, p2_tau_f)
# }
# mydata <- list(N1 = length(nTTP), N2 = length(EFF), Nsub = n.subj,
# y = nTTP, X_y = X_y, e = EFF, keepeff.ind = keepeff.ind,
# W_y = W_y, X_e = X_e, p1_beta_other = c(0, 0),
# p2_beta_other = diag(rep(0.001, 2)),
# p1_beta_dose = 0, p2_beta_dose = 1000,
# p1_alpha = c(0, 0, 0), p2_alpha = diag(rep(0.001, 3)),
# p1_gamma0 = 0, p2_gamma0 = 0.001,
# p1_tau_e = 0, p2_tau_e = 1000,
# p1_tau_u = 0, p2_tau_u = 1000,
# p1_tau_f = 0, p2_tau_f = 1000)
# path.model <- file.path(tempdir(), "model.file.txt")
# R2WinBUGS::write.model(model.file, path.model)
# inits.list <- list(list(beta_other = c(0.1, 0.1),
# beta_dose = 0.1,
# alpha = c(0.1, 0.1, 0.1),
# gamma0 = 0.1,
# tau_e = 0.1,
# tau_u = 0.1,
# tau_f = 0.1,
# .RNG.seed = sample(1:1e+06, size = 1),
# .RNG.name = "base::Wichmann-Hill"))
# jagsobj <- rjags::jags.model(path.model, data = mydata, n.chains = 1,
# quiet = TRUE, inits = inits.list)
# update(jagsobj, n.iter = n.iter, progress.bar = "none")
# post.samples <- rjags::jags.samples(jagsobj, c("beta_dose", "beta_other", "alpha",
# "gamma0"),
# n.iter = n.iter,
# progress.bar = "none")
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# sim.alphas <- as.matrix(rbind(post.samples$alpha[,,1]))
# ############redefine allowable doses#############
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= p1 & prob2.doses >= p2)
# if(length(allow.doses) == 0){
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# RAND.EFF <- sapply(allow.doses, function(a){
# mean(apply(sim.alphas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * a^2)
# }))
# })
# RAND.EFF <- exp(RAND.EFF)/sum(exp(RAND.EFF))
# nxtdose <- sample(allow.doses, 1, prob = RAND.EFF)
# ####a condition for untried higher dose level that is predicted to be efficacious####
# if(nxtdose > max(rec_doseA) + 1 & !nxtdose %in% rec_doseA){
# nxtdose <- max(rec_doseA) + 1
# }
# }
# #print('Stage 3');
# ############################################################################
# ###############################Stage 3######################################
# ############################################################################
# ###################when all the data are available##########################
# ############################################################################
# doseA <- nxtdose
# rec_doseA <- c(rec_doseA, doseA)
# altMat[doseA, k] <- altMat[doseA, k] + chSize
# cmPat <- cmPat + chSize
# PatData <- list()
# n <- chSize
# K <- MaxCycle
# ##################################################
# #######design matrix for toxicity#################
# ##################################################
# dose.pat <- rep(doseA, n)
# xx <- expand.grid(dose.pat, 1:K)
# X <- as.matrix(cbind(rep(1, n*K), xx))
# X <- cbind(X, ifelse(X[, 3] == 1, 0, 1))
# colnames(X) <- c("intercept", "dose", "cycle", "icycle")
# ##################################################
# #######design matrix for efficacy#################
# ##################################################
# Xe <- cbind(1, rep(doseA, n), rep(doseA^2, n))
# colnames(Xe) <- c("intercept", "dose", "squared dose")
# #####simulate nTTP and Efficacy for the cohort########
# outcome <- apply(X, 1, function(a){
# tox.by.grade <- tox_matrix[a[2], a[3], ,]
# nttp.indices <- apply(tox.by.grade, 1, function(o){
# sample(1:5, 1, prob = o)
# })
# if(max(wm[cbind(1:nrow(wm), nttp.indices)]) >= 1){
# y.dlt <- 1
# }else{
# y.dlt <- 0
# }
# y.nTTP <- nTTP.all[nttp.indices[1],
# nttp.indices[2],
# nttp.indices[3]]
# return(c(y.dlt = y.dlt,
# y = y.nTTP))
# })
# y <- outcome[2, ]
# y.dlt <- outcome[1, ]
# beta.mean <- eff.structure[Xe[1, 2]]
# beta.sd <- eff.sd
# beta.shape1 <- ((1 - beta.mean)/beta.sd^2 - 1/beta.mean) * beta.mean^2
# beta.shape2 <- beta.shape1 * (1/beta.mean - 1)
# e <- matrix(rbeta(chSize, beta.shape1, beta.shape2), ncol = 1)
# #################################################################
# #################################################################
# #####add to the master Dataset and extract from this dataset for final analysis###
# #####PatData stores the dataset used to estimate the model########################
# temp.mtdata <- data.frame(cbind(y.dlt, y, X, rep(cmPat/chSize, n*K), 1:n))
# temp.mtdata$subID <- paste0("cohort", temp.mtdata[, 7], "subject", temp.mtdata[, 8])
# masterData$toxicity <- data.frame(rbind(masterData$toxicity, temp.mtdata))
# temp.eff.mtdata <- data.frame(cbind(e, Xe, cmPat/chSize))
# masterData$efficacy <- data.frame(rbind(masterData$efficacy, temp.eff.mtdata))
# PatData$toxicity <- masterData$toxicity[, c(1, 2, 3, 4, 5, 6, 9)]
# #####toxicity dropout########
# #####dropout.dlt is 0 means to keep the observation; 1 means to drop it due to dlt
# dlt.drop <- sapply(by(PatData$toxicity, PatData$toxicity$subID, function(a){a[, 1]}), function(o){
# if((length(o) > 1 & all(o[1: (length(o) - 1)] == 0)) | length(o) == 1){
# return(rep(0, length(o)))
# }else{
# o.temp <- rep(0, length(o))
# o.temp[(min(which(o == 1)) + 1) : length(o.temp)] <- 1
# return(o.temp)
# }
# }, simplify = FALSE)
# dropout.dlt <- NULL
# for(i in unique(PatData$toxicity$subID)){
# dropout.dlt[PatData$toxicity$subID == i] <- dlt.drop[[i]]
# }
# PatData$toxicity <- PatData$toxicity[dropout.dlt == 0, ]
# #####efficacy dropout########
# dropout.eff <- sapply(dlt.drop, function(a){
# if(sum(a) == 0){
# return(0)
# }else if(min(which(a == 1)) >= 4){
# return(0)
# }else{
# return(1)
# }
# })
# dropout.eff <- as.numeric(dropout.eff)
# PatData$efficacy <- masterData$efficacy[dropout.eff == 0, ]
# #################################################################
# #################Model Estimation################################
# #################################################################
# nTTP <- PatData$toxicity[, 2]
# X_y <- as.matrix(PatData$toxicity[, c("intercept", "dose", "cycle")])
# n.subj <- length(unique(PatData$toxicity[, "subID"]))
# W_y <- matrix(0, nrow(PatData$toxicity), n.subj)
# W_y[cbind(1:nrow(W_y), as.numeric(sapply(PatData$toxicity$subID, function(a){
# which(a == unique(PatData$toxicity$subID))
# })))] <- 1
# EFF <- PatData$efficacy[, 1]
# X_e <- as.matrix(PatData$efficacy[, c("intercept", "dose", "squared.dose")])
# keepeff.ind <- which(dropout.eff == 0)
# model.file <- function()
# {
# beta <- c(beta_other[1], beta_dose, beta_other[2])
# for(k in 1:Nsub){
# u[k] ~ dnorm(0, tau_u)
# }
# for(i in 1:N1){
# y[i] ~ dnorm(mu[i], tau_e)
# mu[i] <- inprod(X_y[i, ], beta) + inprod(W_y[i, ], u)
# }
# for(j in 1:N2){
# e[j] ~ dnorm(mu_e[j], tau_f)
# mu_e[j] <- inprod(X_e[j, ], alpha) + gamma0 * u[keepeff.ind[j]]
# }
# beta_other ~ dmnorm(p1_beta_other[], p2_beta_other[, ])
# beta_dose ~ dunif(p1_beta_dose, p2_beta_dose)
# alpha ~ dmnorm(p1_alpha[], p2_alpha[, ])
# gamma0 ~ dnorm(p1_gamma0, p2_gamma0)
# tau_e ~ dunif(p1_tau_e, p2_tau_e)
# tau_u ~ dunif(p1_tau_u, p2_tau_u)
# tau_f ~ dunif(p1_tau_f, p2_tau_f)
# }
# mydata <- list(N1 = length(nTTP), N2 = length(EFF), Nsub = n.subj,
# y = nTTP, X_y = X_y, e = EFF, keepeff.ind = keepeff.ind,
# W_y = W_y, X_e = X_e, p1_beta_other = c(0, 0),
# p2_beta_other = diag(rep(0.001, 2)),
# p1_beta_dose = 0, p2_beta_dose = 1000,
# p1_alpha = c(0, 0, 0), p2_alpha = diag(rep(0.001, 3)),
# p1_gamma0 = 0, p2_gamma0 = 0.001,
# p1_tau_e = 0, p2_tau_e = 1000,
# p1_tau_u = 0, p2_tau_u = 1000,
# p1_tau_f = 0, p2_tau_f = 1000)
# path.model <- file.path(tempdir(), "model.file.txt")
# R2WinBUGS::write.model(model.file, path.model)
# inits.list <- list(list(beta_other = c(0.1, 0.1),
# beta_dose = 0.1,
# alpha = c(0.1, 0.1, 0.1),
# gamma0 = 0.1,
# tau_e = 0.1,
# tau_u = 0.1,
# tau_f = 0.1,
# .RNG.seed = sample(1:1e+06, size = 1),
# .RNG.name = "base::Wichmann-Hill"))
# jagsobj <- rjags::jags.model(path.model, data = mydata, n.chains = 1,
# quiet = TRUE, inits = inits.list)
# update(jagsobj, n.iter = n.iter, progress.bar = "none")
# post.samples <- rjags::jags.samples(jagsobj, c("beta_dose", "beta_other", "alpha",
# "gamma0"),
# n.iter = n.iter,
# progress.bar = "none")
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# sim.alphas <- as.matrix(rbind(post.samples$alpha[,,1]))
# ############redefine allowable doses#############
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(mean(sapply(2:K, function(m){
# as.numeric(o[1] + o[2] * a + o[3] * m)
# })) <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= p1 & prob2.doses >= p2)
# if(length(allow.doses) == 0){
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# effcy.doses <- sapply(allow.doses, function(a){
# mean(apply(sim.alphas, 2, function(o){
# o[1] + o[2] * a + o[3] * a^2}))
# })
# proxy.eff.doses <- allow.doses[which(sapply(effcy.doses, function(a){
# abs(a - max(effcy.doses)) <= proxy.thrd
# }))]
# ###recommend the lowest dose that is efficacious####
# recom.doses <- min(proxy.eff.doses)
# #print('output')
# return(results = list(sc = sc, opt.dose = recom.doses,
# altMat = altMat, masterData = masterData,
# stage1.allow = stage1.allow,
# effy.doses = proxy.eff.doses,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
Try the phase1RMD package in your browser
Any scripts or data that you put into this service are public.
phase1RMD documentation built on March 13, 2020, 1:40 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
###################################################################################
###########Simulate Data from the Matrix and Do Stage 1-3##########################
###################################################################################
###model: yij = beta0 + beta1xi + beta2tj + ui + epiij#############################
######### ei = alpha0 + alpha1xi + alpha2xi^2 + gamma0ui + epi2i###################
######### epi ~ N(0, sigma_e); ui ~ N(0, sigma_u); epi2 ~ N(0, sigma_f)############
###################################################################################
#library(rjags)
#library(R2WinBUGS)
#library(BiocParallel)
# SimRMDEff <- function(seed=2441, tox_matrix, strDose = 1, chSize = 3, trlSize = 36, sdose = 1:6, MaxCycle = 6,
# tox.target = 0.28, eff.structure = c(0.1, 0.2, 0.3, 0.4, 0.7, 0.9)){
# # SimRMDEff <- function(seed=2441, tox_matrix, StrDose = 1, chSize = 3, trlSize = 36, sdose = 1:6, MaxCycle = 6,
# # tox.target = 0.28, eff.structure = c(0.1, 0.2, 0.3, 0.4, 0.7, 0.9),
# # eff.sd = 0.2, p1 = 0.2, p2 = 0.2,
# # ps1 = 0.2,
# # proxy.thrd = 0.1, thrd1 = 0.28, thrd2 = 0.28,
# # wm = matrix(c(0, 0.5, 0.75, 1 , 1.5,
# # 0, 0.5, 0.75, 1 , 1.5,
# # 0, 0 , 0 , 0.5, 1 ),
# # byrow = T, ncol = 5),
# # toxmax = 2.5, toxtype = NULL, intercept.alpha = NULL,
# # coef.beta = NULL, cycle.gamma = NULL){
# #revision031917
# #dummy variable defination for winbugs parameters
# beta_other <- 0;
# beta_dose <- 0;
# N1 <- 0;
# inprod <- 0;
# u <- 0;
# N2 <- 0;
# Nsub <- 0;
# alpha <- 0;
# gamma0 <- 0;
# #revision031917
# n.iter <- 3000;
# eff.sd <- 0.2; p1 <- 0.2; p2 <- 0.2;
# ps1 <- 0.2; proxy.thrd <- 0.1; thrd1 <- 0.28; thrd2 <- 0.28;
# wm <- matrix(c(0, 0.5, 0.75, 1 , 1.5,
# 0, 0.5, 0.75, 1 , 1.5,
# 0, 0 , 0 , 0.5, 1 ),
# byrow = T, ncol = 5);
# toxmax <- 2.5; toxtype = NULL; intercept.alpha = NULL;
# coef.beta <- NULL; cycle.gamma <- NULL;
# set.seed(seed)
# trialSize <- trlSize;
# doses <- sdose;
# cycles <- 1:MaxCycle;
# listdo <- paste0("D", doses)
# if(length(eff.structure) != length(doses))
# stop("Check if you have specified the efficacy mean structure for the correct number of doses.")
# if(nrow(wm) != 3)
# stop("Make sure that you specified the clinical weight matrix for three types of toxicities.")
# MaxCycle <- length(cycles)
# Numdose <- length(doses)
# flag.tox.matrix.null <- 0
# if(is.null(tox_matrix)){
# flag.tox.matrix.null <- 1
# tox_matrix <- array(NA, dim = c(Numdose, MaxCycle, 3, 5))
# if(length(toxtype) != 3)
# stop("Right now we are only considering three toxicity types, but we will relax this constraint in the future work.")
# if(length(intercept.alpha) != 4){
# stop("Exactly four intercepts alpha are needed for grade 0--4 in simulation!")
# }
# if(min(diff(intercept.alpha, lag = 1, differences = 1)) < 0){
# stop("Intercepts alpha for simulation must be in a monotonic increasing order!")
# }
# if(length(coef.beta) != 3){
# stop("Exactly three betas need to be specified for three types of toxicities!")
# }
# }
# #############################################
# #########Determin Scenarios##################
# #############################################
# w1 <- wm[1, ]
# w2 <- wm[2, ]
# w3 <- wm[3, ]
# ####Array of the normalized scores########
# nTTP.array <- function(w1,w2,w3){
# nTTP <- array(NA,c(5,5,5))
# for (t1 in 1:5){
# for (t2 in 1:5){
# for (t3 in 1:5){
# nTTP[t1,t2,t3] <-
# sqrt(w1[t1]^2 + w2[t2]^2 + w3[t3]^2)/toxmax
# }
# }
# }
# return(nTTP)
# }
# nTTP.all <- nTTP.array(w1,w2,w3)
# ####compute pdlt given probability array####
# pDLT <- function(proba){
# return(sum(proba[c(4,5), , ])
# + sum(proba[ ,c(4,5), ])
# + sum(proba[ , ,5])
# - sum(proba[c(4,5),,])*sum(proba[,c(4,5),])
# - sum(proba[c(4,5),,])*sum(proba[,,5])
# - sum(proba[,c(4,5),])*sum(proba[,,5])
# + sum(proba[c(4,5),,])*sum(proba[,c(4,5),])*sum(proba[,,5]))
# }
# ####compute mean nTTP given probability array######
# mnTTP <-function(proba){
# return(sum(proba * nTTP.all))
# }
# sc.mat <- matrix(NA, MaxCycle, Numdose)
# pdlt <- matrix(NA, MaxCycle, Numdose)
# for(i in 1:MaxCycle){ ##loop through cycle
# for(j in 1:Numdose){ ##loop through dose
# if(flag.tox.matrix.null == 0){
# nTTP.prob <- tox_matrix[j, i, ,]
# }else{
# # initialize storing matrix
# cump <- matrix(NA, nrow = length(toxtype), ncol = 5) # cumulative probability
# celp <- matrix(NA, nrow = length(toxtype), ncol = 5) # cell probability
# for(k in 1:length(toxtype)){
# # proportional odds model
# logitcp <- intercept.alpha + coef.beta[k] * j + cycle.gamma * (i - 1)
# # cumulative probabilities
# cump[k, 1:4] <- exp(logitcp) / (1 + exp(logitcp))
# cump[k, 5] <- 1
# # cell probabilities
# celp[k, ] <- c(cump[k,1], diff(cump[k, ], lag = 1, differences = 1))
# }
# nTTP.prob <- celp
# tox_matrix[j, i, ,] <- celp
# }
# proba <- outer(outer(nTTP.prob[1, ], nTTP.prob[2, ]), nTTP.prob[3, ])
# sc.mat[i, j] <- round(mnTTP(proba), 4)
# pdlt[i, j] <- round(pDLT(proba), 4)
# }
# }
# mEFF <- eff.structure
# sc <- rbind(sc.mat, mEFF, pdlt[1, ])
# colnames(sc) <- listdo
# rownames(sc) <- c("mnTTP.1st", "mnTTP.2nd", "mnTTP.3rd",
# "mnTTP.4th", "mnTTP.5th", "mnTTP.6th",
# "mEFF", "pDLT")
# ####################################################################################
# k <- 1 ####Let us use k to denote the column of altMat####
# doseA <- strDose ####Let us use doseA to denote the current dose assigned
# cmPat <- chSize ####current patient size
# masterData <- list()
# rec_doseA <- NULL
# stage1.allow <- NULL
# altMat <- matrix(0, Numdose, k)
# altMat[doseA, k] <- chSize
# ###############################################################
# ########################Stage 1################################
# ###############################################################
# #print('Stage 1')
# while(cmPat <= trialSize/2){
# rec_doseA <- c(rec_doseA, doseA)
# PatData <- NULL
# n <- chSize
# K <- MaxCycle
# ##############################################################
# #############design matrix for toxicity#######################
# ##############################################################
# dose.pat <- rep(doseA, n)
# xx <- expand.grid(dose.pat, 1:K)
# X <- as.matrix(cbind(rep(1, n * K), xx))
# X <- cbind(X, ifelse(X[, 3] == 1, 0, 1))
# colnames(X) <- c("intercept", "dose", "cycle", "icycle")
# ##############################################################
# #######design matrix for efficacy#############################
# ##############################################################
# Xe <- cbind(1, rep(doseA, n), rep(doseA^2, n))
# colnames(Xe) <- c("intercept", "dose", "squared dose")
# #####simulate nTTP and Efficacy for the cohort########
# outcome <- apply(X, 1, function(a){
# tox.by.grade <- tox_matrix[a[2], a[3], ,]
# nttp.indices <- apply(tox.by.grade, 1, function(o){
# sample(1:5, 1, prob = o)
# })
# if(max(wm[cbind(1:nrow(wm), nttp.indices)]) >= 1){
# y.dlt <- 1
# }else{
# y.dlt <- 0
# }
# y.nTTP <- nTTP.all[nttp.indices[1],
# nttp.indices[2],
# nttp.indices[3]]
# return(c(y.dlt = y.dlt,
# y = y.nTTP))
# })
# y <- outcome[2, ]
# y.dlt <- outcome[1, ]
# beta.mean <- eff.structure[Xe[1, 2]]
# beta.sd <- eff.sd
# beta.shape1 <- ((1 - beta.mean)/beta.sd^2 - 1/beta.mean) * beta.mean^2
# beta.shape2 <- beta.shape1 * (1/beta.mean - 1)
# e <- matrix(rbeta(chSize, beta.shape1, beta.shape2), ncol = 1)
# #####construct master Dataset and extract from this dataset for each interim#########
# #####PatData stores the dataset used to estimate the model at each interim###########
# temp.mtdata <- data.frame(cbind(y.dlt, y, X, rep(cmPat/chSize, n * K), 1:n))
# temp.mtdata$subID <- paste0("cohort", temp.mtdata[, 7], "subject", temp.mtdata[, 8])
# masterData$toxicity <- data.frame(rbind(masterData$toxicity, temp.mtdata))
# temp.eff.mtdata <- data.frame(cbind(e, Xe, cmPat/chSize))
# masterData$efficacy <- data.frame(rbind(masterData$efficacy, temp.eff.mtdata))
# current_cohort <- cmPat/chSize
# PatData.index <- data.frame(cohort = 1:current_cohort, cycles = current_cohort:1)
# PatData.index <- within(PatData.index, cycles <- ifelse(cycles >= MaxCycle, MaxCycle, cycles))
# for(i in 1:nrow(PatData.index)){
# PatData <- rbind(PatData, masterData$toxicity[masterData$toxicity[, 7] == PatData.index[i, "cohort"] &
# masterData$toxicity[, 5] <= PatData.index[i, "cycles"],
# c(1, 2, 3, 4, 5, 6, 9)])
# }
# #####dropout.dlt is 0 means to keep the observation; 1 means to drop it due to dlt#######
# dlt.drop <- sapply(by(PatData, PatData$subID, function(a){a[, 1]}), function(o){
# if((length(o) > 1 & all(o[1: (length(o) - 1)] == 0)) | length(o) == 1){
# return(rep(0, length(o)))
# }else{
# o.temp <- rep(0, length(o))
# o.temp[(min(which(o == 1)) + 1) : length(o.temp)] <- 1
# return(o.temp)
# }
# })
# dropout.dlt <- NULL
# for(i in unique(PatData$subID)){
# dropout.dlt[PatData$subID == i] <- dlt.drop[[i]]
# }
# PatData <- PatData[dropout.dlt == 0, ]
# ###################################################################################
# ################Model Estimation given PatData#####################################
# ################Stage 1 only considers the toxicity model##########################
# ###################################################################################
# nTTP <- PatData[, 2]
# X_y <- as.matrix(PatData[, c("intercept", "dose", "icycle")])
# n.subj <- length(unique(PatData[, "subID"]))
# W_y <- matrix(0, nrow(PatData), n.subj)
# W_y[cbind(1:nrow(W_y), as.numeric(sapply(PatData$subID, function(a){
# which(a == unique(PatData$subID))
# })))] <- 1
# model.file <- function()
# {
# beta <- c(beta_other[1], beta_dose, beta_other[2])
# for(i in 1:N1){
# y[i] ~ dnorm(mu[i], tau_e)
# mu[i] <- inprod(X_y[i, ], beta) + inprod(W_y[i, ], u)
# }
# for(j in 1:N2){
# u[j] ~ dnorm(0, tau_u)
# }
# beta_other ~ dmnorm(p1_beta_other[], p2_beta_other[, ])
# beta_dose ~ dunif(p1_beta_dose, p2_beta_dose)
# tau_e ~ dunif(p1_tau_e, p2_tau_e)
# tau_u ~ dunif(p1_tau_u, p2_tau_u)
# }
# mydata <- list(N1 = length(nTTP), N2 = n.subj, y = nTTP, X_y = X_y,
# W_y = W_y, p1_beta_other = c(0, 0),
# p2_beta_other = diag(rep(0.001, 2)),
# p1_beta_dose = 0, p2_beta_dose = 1000,
# p1_tau_e = 0, p2_tau_e = 1000,
# p1_tau_u = 0, p2_tau_u = 1000)
# path.model <- file.path(tempdir(), "model.file.txt")
# R2WinBUGS::write.model(model.file, path.model)
# inits.list <- list(list(beta_other = c(0.1, 0.1),
# beta_dose = 0.1,
# tau_e = 0.1,
# tau_u = 0.1,
# .RNG.seed = sample(1:1e+06, size = 1),
# .RNG.name = "base::Wichmann-Hill"))
# jagsobj <- rjags::jags.model(path.model, data = mydata, n.chains = 1,
# quiet = TRUE, inits = inits.list)
# update(jagsobj, n.iter = n.iter, progress.bar = "none")
# post.samples <- rjags::jags.samples(jagsobj, c("beta_dose", "beta_other"),
# n.iter = n.iter,
# progress.bar = "none")
# if(cmPat == trialSize/2){
# ######################################################################
# #############define allowable doses###################################
# ######################################################################
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= p1 & prob2.doses >= p2)
# if(length(allow.doses) == 0){
# stage1.allow <- 0
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# stage1.allow <- allow.doses
# }else{
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1]))
# loss.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# abs(o[1] + o[2] * a - tox.target)
# }))
# })
# nxtdose <- doses[which.min(loss.doses)]
# if(as.numeric(nxtdose) > (doseA + 1)){
# doseA <- doseA + 1;
# }else{
# doseA <- as.numeric(nxtdose);
# }
# if (cmPat == chSize){
# dlt <- with(masterData$toxicity, y.dlt[cycle == 1 & V7 == cmPat/chSize])
# if (sum(dlt) == 0){
# doseA <- 2
# dose_flag <- 0
# }else{
# doseA <- 1
# dose_flag <- 1
# }
# }
# if ((cmPat >= chSize*2) & (dose_flag == 1)){
# dlt <- with(masterData$toxicity, y.dlt[cycle == 1 & V7 == cmPat/chSize])
# if (sum(dlt) == 0){
# doseA <- 2
# dose_flag <- 0
# }else{
# doseA <- 1
# dose_flag <- 1
# }
# }
# if(cmPat >= chSize*3){
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= ps1 & prob2.doses >= ps1)
# if(length(allow.doses) == 0){
# stage1.allow <- 0
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# }
# altMat[doseA, k] <- altMat[doseA, k] + chSize
# }
# cmPat <- cmPat + chSize # cmPat for the next cohort (loop)
# }
# #print('Stage 2')
# ###############################################################
# ########################Stage 2################################
# ###############################################################
# #######let us wait until the efficacy data for stage 1 is available######
# #########################################################################
# cmPat <- cmPat - chSize
# PatData <- list()
# PatData.index <- data.frame(cohort = current_cohort:1, cycles = 3:(3 + current_cohort - 1))
# PatData.index <- within(PatData.index, cycles <- ifelse(cycles >= MaxCycle, MaxCycle, cycles))
# for(i in 1:nrow(PatData.index)){
# PatData$toxicity <- rbind(PatData$toxicity, masterData$toxicity[masterData$toxicity[, 7] == PatData.index[i, "cohort"] &
# masterData$toxicity[, 5] <= PatData.index[i, "cycles"],
# c(1, 2, 3, 4, 5, 6, 9)])
# }
# #####change the ordering of the subjects--being consistent######
# temp.toxicity <- NULL
# for(i in 1:current_cohort){
# temp.toxicity <- rbind(temp.toxicity, PatData$toxicity[grepl(paste0("cohort",i,"subject"), PatData$toxicity$subID), ])
# }
# PatData$toxicity <- temp.toxicity
# rm(temp.toxicity)
# #####toxicity dropout########
# #####dropout.dlt is 0 means to keep the observation; 1 means to drop it due to dlt
# dlt.drop <- sapply(by(PatData$toxicity, PatData$toxicity$subID, function(a){a[, 1]}), function(o){
# if((length(o) > 1 & all(o[1: (length(o) - 1)] == 0)) | length(o) == 1){
# return(rep(0, length(o)))
# }else{
# o.temp <- rep(0, length(o))
# o.temp[(min(which(o == 1)) + 1) : length(o.temp)] <- 1
# return(o.temp)
# }
# })
# dropout.dlt <- NULL
# for(i in unique(PatData$toxicity$subID)){
# dropout.dlt[PatData$toxicity$subID == i] <- dlt.drop[[i]]
# }
# PatData$toxicity <- PatData$toxicity[dropout.dlt == 0, ]
# #####efficacy dropout########
# dropout.eff <- sapply(dlt.drop, function(a){
# if(sum(a) == 0){
# return(0)
# }else if(min(which(a == 1)) >= 4){
# return(0)
# }else{
# return(1)
# }
# })
# dropout.eff <- as.numeric(dropout.eff)
# PatData$efficacy <- masterData$efficacy[dropout.eff == 0, ]
# ############Model Estimation to randomize among allowable doses##########
# nTTP <- PatData$toxicity[, 2]
# X_y <- as.matrix(PatData$toxicity[, c("intercept", "dose", "icycle")])
# n.subj <- length(unique(PatData$toxicity[, "subID"]))
# W_y <- matrix(0, nrow(PatData$toxicity), n.subj)
# W_y[cbind(1:nrow(W_y), as.numeric(sapply(PatData$toxicity$subID, function(a){
# which(a == unique(PatData$toxicity$subID))
# })))] <- 1
# EFF <- PatData$efficacy[, 1]
# X_e <- as.matrix(PatData$efficacy[, c("intercept", "dose", "squared.dose")])
# keepeff.ind <- which(dropout.eff == 0)
# model.file <- function()
# {
# beta <- c(beta_other[1], beta_dose, beta_other[2])
# for(i in 1:N1){
# y[i] ~ dnorm(mu[i], tau_e)
# mu[i] <- inprod(X_y[i, ], beta) + inprod(W_y[i, ], u)
# }
# for(k in 1:Nsub){
# u[k] ~ dnorm(0, tau_u)
# }
# for(j in 1:N2){
# e[j] ~ dnorm(mu_e[j], tau_f)
# mu_e[j] <- inprod(X_e[j, ], alpha) + gamma0 * u[keepeff.ind[j]]
# }
# beta_other ~ dmnorm(p1_beta_other[], p2_beta_other[, ])
# beta_dose ~ dunif(p1_beta_dose, p2_beta_dose)
# alpha ~ dmnorm(p1_alpha[], p2_alpha[, ])
# gamma0 ~ dnorm(p1_gamma0, p2_gamma0)
# tau_e ~ dunif(p1_tau_e, p2_tau_e)
# tau_u ~ dunif(p1_tau_u, p2_tau_u)
# tau_f ~ dunif(p1_tau_f, p2_tau_f)
# }
# mydata <- list(N1 = length(nTTP), N2 = length(EFF), y = nTTP,
# Nsub = n.subj, X_y = X_y, e = EFF, keepeff.ind = keepeff.ind,
# W_y = W_y, X_e = X_e, p1_beta_other = c(0, 0),
# p2_beta_other = diag(rep(0.001, 2)),
# p1_beta_dose = 0, p2_beta_dose = 1000,
# p1_alpha = c(0, 0, 0), p2_alpha = diag(rep(0.001, 3)),
# p1_gamma0 = 0, p2_gamma0 = 0.001,
# p1_tau_e = 0, p2_tau_e = 1000,
# p1_tau_u = 0, p2_tau_u = 1000,
# p1_tau_f = 0, p2_tau_f = 1000)
# path.model <- file.path(tempdir(), "model.file.txt")
# R2WinBUGS::write.model(model.file, path.model)
# inits.list <- list(list(beta_other = c(0.1, 0.1),
# beta_dose = 0.1,
# alpha = c(0.1, 0.1, 0.1),
# gamma0 = 0.1,
# tau_e = 0.1,
# tau_u = 0.1,
# tau_f = 0.1,
# .RNG.seed = sample(1:1e+06, size = 1),
# .RNG.name = "base::Wichmann-Hill"))
# jagsobj <- rjags::jags.model(path.model, data = mydata, n.chains = 1,
# quiet = TRUE, inits = inits.list)
# update(jagsobj, n.iter = n.iter, progress.bar = "none")
# post.samples <- rjags::jags.samples(jagsobj, c("beta_dose", "beta_other", "alpha",
# "gamma0"),
# n.iter = n.iter,
# progress.bar = "none")
# ######################################################################
# #############update allowable doses###################################
# ######################################################################
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= p1 & prob2.doses >= p2)
# if(length(allow.doses) == 0){
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# sim.alphas <- as.matrix(rbind(post.samples$alpha[,,1]))
# RAND.EFF <- sapply(allow.doses, function(a){
# mean(apply(sim.alphas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * a^2)
# }))
# })
# RAND.EFF <- exp(RAND.EFF)/sum(exp(RAND.EFF))
# nxtdose <- sample(allow.doses, 1, prob = RAND.EFF)
# ####a condition for untried higher dose level that is randomized#### rec_doseA[length(rec_doseA)]
# if(nxtdose > max(rec_doseA) + 1 & !nxtdose %in% rec_doseA){
# nxtdose <- max(rec_doseA) + 1
# }
# #################################################################
# ########generate data for the new enrolled cohort in Stage 2#####
# #################################################################
# while(cmPat < trialSize - chSize){
# doseA <- nxtdose
# rec_doseA <- c(rec_doseA, doseA)
# altMat[doseA, k] <- altMat[doseA, k] + chSize
# cmPat <- cmPat + chSize
# PatData <- list()
# n <- chSize
# K <- MaxCycle
# #######################################################
# #################design matrix for toxicity############
# #######################################################
# dose.pat <- rep(doseA, n)
# xx <- expand.grid(dose.pat, 1:K)
# X <- as.matrix(cbind(rep(1, n*K), xx))
# X <- cbind(X, ifelse(X[, 3] == 1, 0, 1))
# colnames(X) <- c("intercept", "dose", "cycle", "icycle")
# #########################################################
# #######design matrix for efficacy########################
# #########################################################
# Xe <- cbind(1, rep(doseA, n), rep(doseA^2, n))
# colnames(Xe) <- c("intercept", "dose", "squared dose")
# #####simulate nTTP and Efficacy for the cohort########
# outcome <- apply(X, 1, function(a){
# tox.by.grade <- tox_matrix[a[2], a[3], ,]
# nttp.indices <- apply(tox.by.grade, 1, function(o){
# sample(1:5, 1, prob = o)
# })
# if(max(wm[cbind(1:nrow(wm), nttp.indices)]) >= 1){
# y.dlt <- 1
# }else{
# y.dlt <- 0
# }
# y.nTTP <- nTTP.all[nttp.indices[1],
# nttp.indices[2],
# nttp.indices[3]]
# return(c(y.dlt = y.dlt,
# y = y.nTTP))
# })
# y <- outcome[2, ]
# y.dlt <- outcome[1, ]
# beta.mean <- eff.structure[Xe[1, 2]]
# beta.sd <- eff.sd
# beta.shape1 <- ((1 - beta.mean)/beta.sd^2 - 1/beta.mean) * beta.mean^2
# beta.shape2 <- beta.shape1 * (1/beta.mean - 1)
# e <- matrix(rbeta(chSize, beta.shape1, beta.shape2), ncol = 1)
# #################################################################
# #################################################################
# #####add to the master Dataset and extract from this dataset for each interim#########
# #####PatData stores the dataset used to estimate the model at each interim############
# temp.mtdata <- data.frame(cbind(y.dlt, y, X, rep(cmPat/chSize, n*K), 1:n))
# temp.mtdata$subID <- paste0("cohort", temp.mtdata[, 7], "subject", temp.mtdata[, 8])
# masterData$toxicity <- data.frame(rbind(masterData$toxicity, temp.mtdata))
# temp.eff.mtdata <- data.frame(cbind(e, Xe, cmPat/chSize))
# masterData$efficacy <- data.frame(rbind(masterData$efficacy, temp.eff.mtdata))
# current_cohort <- cmPat/chSize
# PatData.index <- data.frame(cohort = 1:current_cohort, cycles = current_cohort:1)
# cycles.adj <- c(rep(2, trialSize/(2 * chSize)), rep(0, current_cohort - trialSize/(2 * chSize)))
# PatData.index <- within(PatData.index, {cycles <- cycles + cycles.adj
# cycles <- ifelse(cycles >= MaxCycle, MaxCycle, cycles)})
# for(i in 1:nrow(PatData.index)){
# PatData$toxicity <- rbind(PatData$toxicity, masterData$toxicity[masterData$toxicity[, 7] == PatData.index[i, "cohort"] &
# masterData$toxicity[, 5] <= PatData.index[i, "cycles"],
# c(1, 2, 3, 4, 5, 6, 9)])
# }
# #####toxicity dropout########
# #####dropout.dlt is 0 means to keep the observation; 1 means to drop it due to dlt
# dlt.drop <- sapply(by(PatData$toxicity, PatData$toxicity$subID, function(a){a[, 1]}), function(o){
# if((length(o) > 1 & all(o[1: (length(o) - 1)] == 0)) | length(o) == 1){
# return(rep(0, length(o)))
# }else{
# o.temp <- rep(0, length(o))
# o.temp[(min(which(o == 1)) + 1) : length(o.temp)] <- 1
# return(o.temp)
# }
# })
# dropout.dlt <- NULL
# for(i in unique(PatData$toxicity$subID)){
# dropout.dlt[PatData$toxicity$subID == i] <- dlt.drop[[i]]
# }
# PatData$toxicity <- PatData$toxicity[dropout.dlt == 0, ]
# cohort.eff.index <- with(PatData.index, which(cycles >= 3))
# #####efficacy dropout########
# dropout.eff <- sapply(dlt.drop, function(a){
# if(sum(a) == 0){
# return(0)
# }else if(min(which(a == 1)) >= 4){
# return(0)
# }else{
# return(1)
# }
# })
# dropout.eff <- as.numeric(dropout.eff)
# PatData$efficacy <- subset(masterData$efficacy, masterData$efficacy[, 5] %in% cohort.eff.index & dropout.eff == 0)
# #############################################################################################################################
# ######################Model Estimation to update toxicity information and randomization probs################################
# #############################################################################################################################
# nTTP <- PatData$toxicity[, 2]
# X_y <- as.matrix(PatData$toxicity[, c("intercept", "dose", "icycle")])
# n.subj <- length(unique(PatData$toxicity[, "subID"]))
# W_y <- matrix(0, nrow(PatData$toxicity), n.subj)
# W_y[cbind(1:nrow(W_y), as.numeric(sapply(PatData$toxicity$subID, function(a){
# which(a == unique(PatData$toxicity$subID))
# })))] <- 1
# EFF <- PatData$efficacy[, 1]
# X_e <- as.matrix(PatData$efficacy[, c("intercept", "dose", "squared.dose")])
# keepeff.ind <- which(dropout.eff == 0 & masterData$efficacy[, 5] %in% cohort.eff.index)
# model.file <- function()
# {
# beta <- c(beta_other[1], beta_dose, beta_other[2])
# for(k in 1:Nsub){
# u[k] ~ dnorm(0, tau_u)
# }
# for(i in 1:N1){
# y[i] ~ dnorm(mu[i], tau_e)
# mu[i] <- inprod(X_y[i, ], beta) + inprod(W_y[i, ], u)
# }
# for(j in 1:N2){
# e[j] ~ dnorm(mu_e[j], tau_f)
# mu_e[j] <- inprod(X_e[j, ], alpha) + gamma0 * u[keepeff.ind[j]]
# }
# beta_other ~ dmnorm(p1_beta_other[], p2_beta_other[, ])
# beta_dose ~ dunif(p1_beta_dose, p2_beta_dose)
# alpha ~ dmnorm(p1_alpha[], p2_alpha[, ])
# gamma0 ~ dnorm(p1_gamma0, p2_gamma0)
# tau_e ~ dunif(p1_tau_e, p2_tau_e)
# tau_u ~ dunif(p1_tau_u, p2_tau_u)
# tau_f ~ dunif(p1_tau_f, p2_tau_f)
# }
# mydata <- list(N1 = length(nTTP), N2 = length(EFF), Nsub = n.subj,
# y = nTTP, X_y = X_y, e = EFF, keepeff.ind = keepeff.ind,
# W_y = W_y, X_e = X_e, p1_beta_other = c(0, 0),
# p2_beta_other = diag(rep(0.001, 2)),
# p1_beta_dose = 0, p2_beta_dose = 1000,
# p1_alpha = c(0, 0, 0), p2_alpha = diag(rep(0.001, 3)),
# p1_gamma0 = 0, p2_gamma0 = 0.001,
# p1_tau_e = 0, p2_tau_e = 1000,
# p1_tau_u = 0, p2_tau_u = 1000,
# p1_tau_f = 0, p2_tau_f = 1000)
# path.model <- file.path(tempdir(), "model.file.txt")
# R2WinBUGS::write.model(model.file, path.model)
# inits.list <- list(list(beta_other = c(0.1, 0.1),
# beta_dose = 0.1,
# alpha = c(0.1, 0.1, 0.1),
# gamma0 = 0.1,
# tau_e = 0.1,
# tau_u = 0.1,
# tau_f = 0.1,
# .RNG.seed = sample(1:1e+06, size = 1),
# .RNG.name = "base::Wichmann-Hill"))
# jagsobj <- rjags::jags.model(path.model, data = mydata, n.chains = 1,
# quiet = TRUE, inits = inits.list)
# update(jagsobj, n.iter = n.iter, progress.bar = "none")
# post.samples <- rjags::jags.samples(jagsobj, c("beta_dose", "beta_other", "alpha",
# "gamma0"),
# n.iter = n.iter,
# progress.bar = "none")
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# sim.alphas <- as.matrix(rbind(post.samples$alpha[,,1]))
# ############redefine allowable doses#############
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= p1 & prob2.doses >= p2)
# if(length(allow.doses) == 0){
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# RAND.EFF <- sapply(allow.doses, function(a){
# mean(apply(sim.alphas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * a^2)
# }))
# })
# RAND.EFF <- exp(RAND.EFF)/sum(exp(RAND.EFF))
# nxtdose <- sample(allow.doses, 1, prob = RAND.EFF)
# ####a condition for untried higher dose level that is predicted to be efficacious####
# if(nxtdose > max(rec_doseA) + 1 & !nxtdose %in% rec_doseA){
# nxtdose <- max(rec_doseA) + 1
# }
# }
# #print('Stage 3');
# ############################################################################
# ###############################Stage 3######################################
# ############################################################################
# ###################when all the data are available##########################
# ############################################################################
# doseA <- nxtdose
# rec_doseA <- c(rec_doseA, doseA)
# altMat[doseA, k] <- altMat[doseA, k] + chSize
# cmPat <- cmPat + chSize
# PatData <- list()
# n <- chSize
# K <- MaxCycle
# ##################################################
# #######design matrix for toxicity#################
# ##################################################
# dose.pat <- rep(doseA, n)
# xx <- expand.grid(dose.pat, 1:K)
# X <- as.matrix(cbind(rep(1, n*K), xx))
# X <- cbind(X, ifelse(X[, 3] == 1, 0, 1))
# colnames(X) <- c("intercept", "dose", "cycle", "icycle")
# ##################################################
# #######design matrix for efficacy#################
# ##################################################
# Xe <- cbind(1, rep(doseA, n), rep(doseA^2, n))
# colnames(Xe) <- c("intercept", "dose", "squared dose")
# #####simulate nTTP and Efficacy for the cohort########
# outcome <- apply(X, 1, function(a){
# tox.by.grade <- tox_matrix[a[2], a[3], ,]
# nttp.indices <- apply(tox.by.grade, 1, function(o){
# sample(1:5, 1, prob = o)
# })
# if(max(wm[cbind(1:nrow(wm), nttp.indices)]) >= 1){
# y.dlt <- 1
# }else{
# y.dlt <- 0
# }
# y.nTTP <- nTTP.all[nttp.indices[1],
# nttp.indices[2],
# nttp.indices[3]]
# return(c(y.dlt = y.dlt,
# y = y.nTTP))
# })
# y <- outcome[2, ]
# y.dlt <- outcome[1, ]
# beta.mean <- eff.structure[Xe[1, 2]]
# beta.sd <- eff.sd
# beta.shape1 <- ((1 - beta.mean)/beta.sd^2 - 1/beta.mean) * beta.mean^2
# beta.shape2 <- beta.shape1 * (1/beta.mean - 1)
# e <- matrix(rbeta(chSize, beta.shape1, beta.shape2), ncol = 1)
# #################################################################
# #################################################################
# #####add to the master Dataset and extract from this dataset for final analysis###
# #####PatData stores the dataset used to estimate the model########################
# temp.mtdata <- data.frame(cbind(y.dlt, y, X, rep(cmPat/chSize, n*K), 1:n))
# temp.mtdata$subID <- paste0("cohort", temp.mtdata[, 7], "subject", temp.mtdata[, 8])
# masterData$toxicity <- data.frame(rbind(masterData$toxicity, temp.mtdata))
# temp.eff.mtdata <- data.frame(cbind(e, Xe, cmPat/chSize))
# masterData$efficacy <- data.frame(rbind(masterData$efficacy, temp.eff.mtdata))
# PatData$toxicity <- masterData$toxicity[, c(1, 2, 3, 4, 5, 6, 9)]
# #####toxicity dropout########
# #####dropout.dlt is 0 means to keep the observation; 1 means to drop it due to dlt
# dlt.drop <- sapply(by(PatData$toxicity, PatData$toxicity$subID, function(a){a[, 1]}), function(o){
# if((length(o) > 1 & all(o[1: (length(o) - 1)] == 0)) | length(o) == 1){
# return(rep(0, length(o)))
# }else{
# o.temp <- rep(0, length(o))
# o.temp[(min(which(o == 1)) + 1) : length(o.temp)] <- 1
# return(o.temp)
# }
# }, simplify = FALSE)
# dropout.dlt <- NULL
# for(i in unique(PatData$toxicity$subID)){
# dropout.dlt[PatData$toxicity$subID == i] <- dlt.drop[[i]]
# }
# PatData$toxicity <- PatData$toxicity[dropout.dlt == 0, ]
# #####efficacy dropout########
# dropout.eff <- sapply(dlt.drop, function(a){
# if(sum(a) == 0){
# return(0)
# }else if(min(which(a == 1)) >= 4){
# return(0)
# }else{
# return(1)
# }
# })
# dropout.eff <- as.numeric(dropout.eff)
# PatData$efficacy <- masterData$efficacy[dropout.eff == 0, ]
# #################################################################
# #################Model Estimation################################
# #################################################################
# nTTP <- PatData$toxicity[, 2]
# X_y <- as.matrix(PatData$toxicity[, c("intercept", "dose", "cycle")])
# n.subj <- length(unique(PatData$toxicity[, "subID"]))
# W_y <- matrix(0, nrow(PatData$toxicity), n.subj)
# W_y[cbind(1:nrow(W_y), as.numeric(sapply(PatData$toxicity$subID, function(a){
# which(a == unique(PatData$toxicity$subID))
# })))] <- 1
# EFF <- PatData$efficacy[, 1]
# X_e <- as.matrix(PatData$efficacy[, c("intercept", "dose", "squared.dose")])
# keepeff.ind <- which(dropout.eff == 0)
# model.file <- function()
# {
# beta <- c(beta_other[1], beta_dose, beta_other[2])
# for(k in 1:Nsub){
# u[k] ~ dnorm(0, tau_u)
# }
# for(i in 1:N1){
# y[i] ~ dnorm(mu[i], tau_e)
# mu[i] <- inprod(X_y[i, ], beta) + inprod(W_y[i, ], u)
# }
# for(j in 1:N2){
# e[j] ~ dnorm(mu_e[j], tau_f)
# mu_e[j] <- inprod(X_e[j, ], alpha) + gamma0 * u[keepeff.ind[j]]
# }
# beta_other ~ dmnorm(p1_beta_other[], p2_beta_other[, ])
# beta_dose ~ dunif(p1_beta_dose, p2_beta_dose)
# alpha ~ dmnorm(p1_alpha[], p2_alpha[, ])
# gamma0 ~ dnorm(p1_gamma0, p2_gamma0)
# tau_e ~ dunif(p1_tau_e, p2_tau_e)
# tau_u ~ dunif(p1_tau_u, p2_tau_u)
# tau_f ~ dunif(p1_tau_f, p2_tau_f)
# }
# mydata <- list(N1 = length(nTTP), N2 = length(EFF), Nsub = n.subj,
# y = nTTP, X_y = X_y, e = EFF, keepeff.ind = keepeff.ind,
# W_y = W_y, X_e = X_e, p1_beta_other = c(0, 0),
# p2_beta_other = diag(rep(0.001, 2)),
# p1_beta_dose = 0, p2_beta_dose = 1000,
# p1_alpha = c(0, 0, 0), p2_alpha = diag(rep(0.001, 3)),
# p1_gamma0 = 0, p2_gamma0 = 0.001,
# p1_tau_e = 0, p2_tau_e = 1000,
# p1_tau_u = 0, p2_tau_u = 1000,
# p1_tau_f = 0, p2_tau_f = 1000)
# path.model <- file.path(tempdir(), "model.file.txt")
# R2WinBUGS::write.model(model.file, path.model)
# inits.list <- list(list(beta_other = c(0.1, 0.1),
# beta_dose = 0.1,
# alpha = c(0.1, 0.1, 0.1),
# gamma0 = 0.1,
# tau_e = 0.1,
# tau_u = 0.1,
# tau_f = 0.1,
# .RNG.seed = sample(1:1e+06, size = 1),
# .RNG.name = "base::Wichmann-Hill"))
# jagsobj <- rjags::jags.model(path.model, data = mydata, n.chains = 1,
# quiet = TRUE, inits = inits.list)
# update(jagsobj, n.iter = n.iter, progress.bar = "none")
# post.samples <- rjags::jags.samples(jagsobj, c("beta_dose", "beta_other", "alpha",
# "gamma0"),
# n.iter = n.iter,
# progress.bar = "none")
# sim.betas <- as.matrix(rbind(post.samples$beta_other[1,,1], post.samples$beta_dose[,,1],
# post.samples$beta_other[2,,1]))
# sim.alphas <- as.matrix(rbind(post.samples$alpha[,,1]))
# ############redefine allowable doses#############
# ####condition 1####
# prob1.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(o[1] + o[2] * a + o[3] * 1 <= thrd1)
# }))
# })
# ####condition 2####
# prob2.doses <- sapply(doses, function(a){
# mean(apply(sim.betas, 2, function(o){
# as.numeric(mean(sapply(2:K, function(m){
# as.numeric(o[1] + o[2] * a + o[3] * m)
# })) <= thrd2)
# }))
# })
# allow.doses <- which(prob1.doses >= p1 & prob2.doses >= p2)
# if(length(allow.doses) == 0){
# return(results = list(sc = sc, opt.dose = 0,
# altMat = altMat,
# masterData = masterData,
# stage1.allow = stage1.allow,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
# effcy.doses <- sapply(allow.doses, function(a){
# mean(apply(sim.alphas, 2, function(o){
# o[1] + o[2] * a + o[3] * a^2}))
# })
# proxy.eff.doses <- allow.doses[which(sapply(effcy.doses, function(a){
# abs(a - max(effcy.doses)) <= proxy.thrd
# }))]
# ###recommend the lowest dose that is efficacious####
# recom.doses <- min(proxy.eff.doses)
# #print('output')
# return(results = list(sc = sc, opt.dose = recom.doses,
# altMat = altMat, masterData = masterData,
# stage1.allow = stage1.allow,
# effy.doses = proxy.eff.doses,
# dose_trace = rec_doseA,
# cmPat = cmPat))
# }
Try the phase1RMD package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.