# -----------------------------------------------------------------
# Filename: BLR-Dual-Post.odc
#
# Description: Bayesian Model for dual-agent phase I dose-escalation trials
#
# Model: 5 parameter model
#   odds12 = odds12.0*exp(eta x d1/d1* x d2/d2*)     
#   odds12.0 = odds1 + odds2 + odds1*odds2 (independence)
#   log(odds1) = logAlpha1 + Beta1*log(Dose1/DoseRef1)
#   log(odds2) = logAlpha2 + Beta2*log(Dose2/DoseRef2)
# 
# Priors:
#   BVN for (logAlpha1,logBeta1) and (logAlpha2,logBeta2)
#   Normal prior for eta
#
# Reference: 
#   Neuenschwander, Matano, Tang, Roychoudhury, Wandel and Bailey. 
#   A Bayesian Industry Approach to Phase I Combination Trials in Oncology. 
#   Statistical Methods in Drug Combination Studies, Boca Raton, FL:   
#   Chapman & Hall/CRC Press 2015. 
#
# Nodes to monitor: 
#    P12: probabiliy of toxicity, vector of length Ndoses1 x Ndoses2
#    pCat: category indicators, Ndoses1 x Ndoses2 x Ncat matrix
#    logAlphaBeta: model parameters, 
#    logAlphaBeta1[1]: logAlpha1,
#    logAlphaBeta1[2]: logBeta1
#    logAlphaBeta2[1]: logAlpha2,
#    logAlphaBeta2[2]: logBeta2
#    eta: interaction odds multiplier
# -----------------------------------------------------------------------

model{
	#  bivariate normal prior for (logAlpha1,logBeta1) [marginal model for Agent1]
	#  Prior1 = (mean.logAlpha1, mean.logBeta1, sd.logAlpha1, sd.logBeta1, corr1)

  #  Covariance and precision matrix for agent 1
  cova1[1,1] <- Prior1[3]*Prior1[3]
  cova1[2,2] <- Prior1[4]*Prior1[4]
  cova1[1,2] <- Prior1[3]*Prior1[4]*Prior1[5]
  cova1[2,1] <- cova1[1,2]
  prec1[1:2,1:2] <- inverse(cova1[,])
 
  #  Prior distribution for logAlpha1,logBeta1
  logAlphaBeta1[1:2] ~ dmnorm(Prior1[1:2],prec1[1:2,1:2])

  #  bivariate normal prior for (logAlpha2,logBeta2) [marginal model for Agent2]
	#  Prior2 = (mean.logAlpha2, mean.logBeta2, sd.logAlpha2, sd.logBeta2, corr2)
  cova2[1,1] <- Prior2[3]*Prior2[3]
  cova2[2,2] <- Prior2[4]*Prior2[4]
  cova2[1,2] <- Prior2[3]*Prior2[4]*Prior2[5]
  cova2[2,1] <- cova2[1,2]
  prec2[1:2,1:2] <- inverse(cova2[,])
  logAlphaBeta2[1:2] ~ dmnorm(Prior2[1:2],prec2[1:2,1:2])

  # Normal prior for interaction parameter (eta) 
     etaPriorPrec <- pow(etaPrior[2], -2)
   eta  ~ dnorm(etaPrior[1], etaPriorPrec)
   OddsFactor <- exp(eta)

 #  define Puctoffs1, with 1 as last entry
  for (j in 1:(Ncat-1)) { Pcutoffs1[j] <- Pcutoffs[j] }
  Pcutoffs1[Ncat] <- 1

  for (i in 1:Ndoses1){
 #  Marginal odds for Agent 1
    lin1.1[i] <- logAlphaBeta1[1] + exp(logAlphaBeta1[2])*log(Doses1[i]/DoseRef1)
    odds1[i] <- exp(lin1.1[i]) 
                       }
 for(i in 1:Ndoses2){      
 #  Marginal odds for Agent 2 
    lin2.1[i] <- logAlphaBeta2[1] + exp(logAlphaBeta2[2])*log(Doses2[i]/DoseRef2)
    odds2[i] <- exp(lin2.1[i])
                     } 

#  combination model and probability of each category for each combination
 for(i in 1:Ndoses1){
   for(j in 1:Ndoses2){
#  Odds without and with interaction 
      odds12.0[i,j] <- odds1[i] + odds2[j] + odds1[i] * odds2[j]
     odds12[i,j] <- odds12.0[i,j] * 
        exp(eta *(Doses1[i]/DoseRef1)*(Doses2[j]/DoseRef2))
     P12[i,j] <- odds12[i,j]/(1 + odds12[i,j])

  # for each dose, indicators for different toxicity categories
  # (means of these indicators correspond to category probabilities)   
    pCat[i,j,1] <- step(Pcutoffs1[1] - P12[i,j])
    for (k in 2:Ncat) {
      pCat[i,j,k] <- step(Pcutoffs1[k] - P12[i,j]) - sum(pCat[i,j,1:(k-1)])
                      }
                  }
               }

 #binomial likelihoods (data)
for (i in 1:Ncohorts){
   lin1.2[i] <- logAlphaBeta1[1] + exp(logAlphaBeta1[2])*log(DosesAdm1[i]/DoseRef1)
        lin1[i] <- lin1.2[i] + step(-20 - lin1.2[i]) * (-20 - lin1.2[i]) +  
                step(lin1.2[i] - 20) * (20 - lin1.2[i])
   odds11[i] <- exp(lin1[i]) 

   lin2.2[i] <- logAlphaBeta2[1] + exp(logAlphaBeta2[2])*log(DosesAdm2[i]/DoseRef2)
   lin2[i] <- lin2.2[i] + step(-20 - lin2.2[i]) * (-20 - lin2.2[i]) +  
                 step(lin2.2[i] - 20) * (20 - lin2.2[i]) 
   odds22[i] <- exp(lin2[i]) 
      
    odds121.0[i] <- odds11[i] + odds22[i] + odds11[i] * odds22[i]         
   odds121[i] <-  odds121.0[i] *  
                    exp(eta*(DosesAdm1[i]/DoseRef1)*(DosesAdm2[i]/DoseRef2))
      P12.1[i] <- odds121[i]/(1 + odds121[i])

      Ntox[i] ~ dbin(P12.1[i], Npat[i]) 
                   } 
}


# ---------------------------------------------------------
# Data (Application 2: Cohort 5 Example from Book Chapter)
# ---------------------------------------------------------

#Analysis 1

list(
  Prior1     = c(-3.146, 0.388,1.214, 0.886,-0.63),
  Prior2     = c(-1.961, 0.690,0.541, 0.650, 0.04),
  etaPrior   = c(0,1.121),

# actual data: number of cohorts, administered dose of the respective
# cohort, number of toxicities, number of patients
# for data scenarios, change input in this part only
  Ncohorts = 3,    
  DosesAdm1 = c(    3,   3,   6),
  DosesAdm2 = c(  400, 800, 400),
  Npat = c(         3,   3,   3),
  Ntox = c(         0,   1,   1),
# doses of interest (i.e. for which inference will be shown) and reference dose
   Ndoses1 = 4,
   Ndoses2 = 4,
   Doses1 = c(0, 3, 4.5, 6),
   Doses2 = c(0, 400, 600, 800),
   DoseRef1 = 3,
   DoseRef2 = 960,
# interval categories
    Ncat = 3,
    Pcutoffs = c(0.16,0.35)
)

#Analysis 2

# Data list
list(
  Prior1     = c(-3.146, 0.388,1.214, 0.886,-0.63),
  Prior2     = c(-1.961, 0.690,0.541, 0.650, 0.04),
  etaPrior   = c(0,1.121),

# actual data: number of cohorts, administered dose of the respective
# cohort, number of toxicities, number of patients
# for data scenarios, change input in this part only
  Ncohorts = 4,    
  DosesAdm1 = c(    3,   3,   6, 4.5),
  DosesAdm2 = c(  400, 800, 400, 600),
  Npat = c(         3,   6,  10,  10),
  Ntox = c(         0,   2,   3,   2),
# doses of interest (i.e. for which inference will be shown) and reference dose
   Ndoses1 = 4,
   Ndoses2 = 4,
   Doses1 = c(0, 3, 4.5, 6),
   Doses2 = c(0, 400, 600, 800),
   DoseRef1 = 3,
   DoseRef2 = 960,
# interval categories
    Ncat = 3,
    Pcutoffs = c(0.16,0.35)
)

# ---------------------------------------------------------
# Initial Values
# ---------------------------------------------------------
# 4 chains with the initial values given below
# 25'000 iterations each chain
#  5'000 burn-in followed by 20'000 updates
# thin = 1
# seed set to 314159
list( logAlphaBeta1 = c(-1.13, -2.12),logAlphaBeta2 = c(-1.13, -2.12),eta=0)
list( logAlphaBeta1 = c(-0.49, -2.30),logAlphaBeta2 = c(-0.49, -2.30),eta=0)
list( logAlphaBeta1 = c(-4.55, -1.84),logAlphaBeta2 = c(-4.55, -1.84),eta=0)
list( logAlphaBeta1 = c(-1.34, -2.29),logAlphaBeta2 = c(-1.34, -2.29),eta=0)

