#********************************************************************************
# Description: To calculate the statistical quantities for the general A+B with
#              or without dose de-escalation discussed in the paper
#              "Statistical Properties of The Traditional Algorithm-based Designs 
#               for Phase I Cancer Clinical Trials" by Yong Lin and Weichung Jow Shih      
# Date: August 10, 1999  yl                                      
# Modified on: 10/99, 3/17/00 yl
# Reasons  of modification: 
#     10/99
#         Modify the program for A+B with and without de-escalation
#********************************************************************************
#
pmtd <-
function(ptox, design = c(3, 3, 1, 1, 1), deescalate = 1)
{
#
# ptox   = probability of toxicity at each dose level
# design = (A,B,C,D,E) the general A+B design with 
#           parameters C, D, and E (where 0<=C<=D<=E,D<=A,E<=A+B,A>0,B>0)
# deescalate = 0      without de-escalation
#            = 1      with de-escalation
#
	n <- length(ptox)
	A <- design[1]
	B <- design[2]
	CC <- design[3]
	DD <- design[4]
	EE <- design[5]
	if(!(0 <= CC & CC <= DD & DD <= EE & EE <= (A + B) & DD <= A & A > 0 & B > 0)) 
		{
		cat("The parameters should satisfy A>0,B>0,0<=C<=D<=E,D<=A,E<=A+B \n")
		invisible()
	}
	P0 <- pbinom(CC - 1, A, ptox)
	P1 <- pbinom(DD, A, ptox) - P0
	Q0 <- Q1 <- Q2 <- numeric(n)
	tmp1 <- tmp2 <- numeric(DD + 1)
	for(i in 1:n) {
		for(k in 0:DD) {
			tmp1[k + 1] <- dbinom(k, A, ptox[i]) * pbinom(EE - k, B, ptox[i
				])
			tmp2[k + 1] <- dbinom(k, A, ptox[i]) * (1 - pbinom(EE - k, B, 
				ptox[i]))
		}
		Q0[i] <- sum(tmp1[(CC + 1):(DD + 1)])
		Q1[i] <- sum(tmp1[1:CC])
		Q2[i] <- sum(tmp2[1:CC])
	}
	p.mtd <- numeric(n + 1)	# Probability of MTD
	p.go <- P0 + Q0
	p.stop <- 1 - p.go
	if(deescalate == 1) {
		p.dose <- Q0 + Q1
		pmtdnstp <- matrix(0, n, n)
		for(i in 1:(n - 1)) {
			for(k in (i + 1):n) {
				if(i == 1)
				  bfi <- 1
				else bfi <- prod(p.go[1:(i - 1)])
				if(k == i + 1)
				  btwn <- 1
				else btwn <- prod(Q2[(i + 1):(k - 1)])
				p.mtd[i] <- p.mtd[i] + bfi * p.dose[i] * btwn * p.stop[
				  k]
				pmtdnstp[k, i + 1] <- bfi * p.dose[i] * btwn * p.stop[k
				  ]
			}
		}
	}
	else if(deescalate == 0) {
		for(i in 1:(n - 1))
			p.mtd[i] <- prod(p.go[1:i]) * p.stop[i + 1]
	}
	p.mtd[n] <- prod(p.go[1:n])
	p.mtd[n + 1] <- p.stop[1]
	if(deescalate == 1) {
		pmtdnstp[1, 1] <- p.stop[1]
		for(k in 2:n) {
			p.mtd[n + 1] <- p.mtd[n + 1] + prod(Q2[1:(k - 1)]) * p.stop[k]
			pmtdnstp[k, 1] <- prod(Q2[1:(k - 1)]) * p.stop[k]
		}
	}
# find target toxicity level
	pmtd <- c(p.mtd[c( - n,  - (n + 1))], 0)
	theta <- sum(ptox * pmtd)/sum(pmtd)	
	# find estimated number of patients in each dose level
	exp.go <- (A * P0 + (A + B) * Q0)/(P0 + Q0)
	exp.stop <- (A * (1 - P0 - P1) + (A + B) * (P1 - Q0))/p.stop
	exp.num <- numeric(n)
	if(deescalate == 1) {
		exp.val <- array(0, rep(n, 3))
		for(j in 1:n) {
			for(k in j:n) {
				for(i in 1:k) {
				  if(k == j)
				    exp.val[k, i, j] <- exp.stop[j]
				  else if(k > j) {
				    if(i <= j + 1)
				      exp.val[k, i, j] <- A + B
				    else exp.val[k, i, j] <- exp.go[j]
				  }
				}
			}
		}
		for(j in 1:n) {
			exp.num[j] <- sum(exp.val[,  , j] * pmtdnstp) + exp.go[j] * 
				p.mtd[n]
		}
	}
	else if(deescalate == 0) {
		mtd2 <- c(p.mtd[n + 1], p.mtd[ - (n + 1)])
		for(j in 1:n)
			exp.num[j] <- exp.stop[j] * mtd2[j] + exp.go[j] * sum(mtd2[(j + 
				1):(n + 1)])
	}
	exp.tox <- ptox * exp.num
	pct.pat <- exp.num/sum(exp.num) * 100
	summ <- rbind(ptox, p.mtd[ - (n + 1)], exp.num, exp.tox, pct.pat)
	dimnames(summ) <- list(c("Prob. of Tox.: ", "Prob. of MTD: ", 
		"Exp. # of Pat.: ", "Exp. # of Tox.: ", "% of Pat.: "), 1:n)
	summ2 <- rbind(theta, p.mtd[n + 1], sum(exp.num), sum(exp.tox), sum(exp.tox)/
		sum(exp.num))
	dimnames(summ2) <- list(c("Target Toxicity Level: ", 
		"Prob. that MTD < dose 1: ", "Exp. All Pat.: ", "Exp. All Tox.: ", 
		"Exp. All Tox. Rate: "), "")
	names(design) <- c("A", "B", "C", "D", "E")
	return(list(Design = design, Deescalate = as.logical(deescalate), summary1 = 
		round(summ, 3), summary2 = round(summ2, 3)))
}