# This file contains revised versions of these functions:
#   seqMonitor
#   summary.seqDesign
#
# To install these patches, load this file into S-PLUS after
# the SeqTrial module is loaded, using the `source' function from
# the command line or by loading this file as a script file.
#
# The patched versions of the functions will be temporarily loaded
# in the current working directory, then moved to the location
# of the SeqTrial module.
#
# Alternately, to just load the functions, without moving them to the
# location of the SeqTrial module, load only the parts of this file
# between "Begin loading here" and "End loading here".

.SeqTrial.where <- find("seqDesign")
if(length(.SeqTrial.where) == 0)
  stop("The SeqTrial module must be loaded before installing this patch")

if(length(.SeqTrial.where) > 1){
  warning("The seqDesign function exists in multiple locations.  The patch will be loaded into the first location")
  .SeqTrial.where <- .SeqTrial.where[1]
}

.SeqTrial.objectsInPatch <- c(".SeqTrial.patch",
			      "seqMonitor",
			      "summary.seqDesign")


#--------------------------------------------------------------------
# Begin loading here if you do not want objects moved automatically.

.SeqTrial.patch <- "SeqTrial patch version of 3/29/2000"


"seqMonitor" <-
function(x, response, treatment = rep(1, length(response)), future.analyses, 
	constraint.scale = x$monitor$constraint.scale, min.increment = 0.1, 
	maxiter = 20)
{
	if(!is.seqDesign(x))
		stop("x must be a seqMonitor or seqDesign object")
	call <- match.call()
	if(is.seqMonitor(x) && x$monitor$action == "Stop")
		warning("Recommendation at previous analysis was to stop")
	xDes <- x$design
	if(is.null(xDes)) {
		xDes <- x
		x <- c(x, list(design = xDes))
	}
	xMod <- x$model
	prob.model <- xMod$prob.model
	arms <- xMod$arms
	log.transform <- xMod$log.transform
	power <- x$hypothesis$power
	theta0 <- xDes$hypothesis$theta.null
	theta.alt <- xDes$hypothesis$theta.alt
	upper <- theta.alt > theta0
	xPar <- x$parameters
	xStat <- x$statistics
	xMon <- x$monitor	#
# check treatment
	if((arms == 0 && length(unique(treatment)) < 2) || (arms == 1 && 
		length(unique(treatment)) != 1) || (arms == 2 && !all(sort(
		unique(treatment)) == c(0, 1)))) stop(
			"treatment is not of correct form for number of arms")
		#
# check constraint.scale
	if(is.null(constraint.scale))
		constraint.scale <- ifelse1(xDes$parameters$design.family == 
			"E", "E", "X")
	if(!is.element(constraint.scale, c("E", "X", "S", "Z", "P")))
		stop("Unrecognized value for constraint.scale")
	if(constraint.scale != "E" && xDes$parameters$design.family == "E")
		stop("Only error spending monitoring allowed with error spending design family"
			)
	if(prob.model == "N" || prob.model == "L") {
# analyze data for normal and lognormal models
		if(prob.model == "L") {
			response <- log(response)
			if(!log.transform)
				theta0 <- log(theta0)
		}
		if(arms == 1) {
			test <- t.test(response, mu = theta0, alternative = 
				"greater")
			X <- test$estimate
			Z <- test$statistic
			P <- test$p.value
			variance <- var(response)
			ratio <- 1
		}
		else if(arms == 2) {
			if(sum(treatment == 1) < 2 || sum(treatment == 0) < 2
				)
				stop(
				  "One or more arms has too few observations to estimate variance"
				  )
			test <- t.test(response[treatment == 1], response[
				treatment == 0], mu = theta0, alternative = 
				"greater", var.equal = F)
			X <- test$estimate[1] - test$estimate[2]
			Z <- test$statistic
			P <- test$p.value
			variance <- c(var(response[treatment == 1]), var(
				response[treatment == 0]))
			ratio <- rev(table(treatment))
		}
		else if(arms == 0) {
			test <- summary(lm(response ~ treatment))
			X <- test$coeff[2, 1]
			Z <- (X - theta0)/test$coeff[2, 2]
			P <- 1 - pt(Z, df = test$df[2])
			if(missing(variance))
				variance <- test$sigma^2
			ratio <- treatment
		}
		else stop("Invalid probability model for monitoring")
		if(arms > 0 && length(response) < 20)
			warning("variance may be inaccurately estimated from this small sample"
				)
		if(prob.model == "L" && !log.transform) {
			X <- exp(X)
			theta0 <- exp(theta0)
		}
	}
	else if(prob.model == "B") {
# analyze data for binomial proportions
		if(arms == 1) {
			X <- mean(response)
			variance <- X * (1 - X)
			Z <- (X - theta0) * sqrt(length(response)/variance)
			P <- 1 - pnorm(Z)
			ratio <- 1
		}
		else if(arms == 2) {
			p1 <- mean(response[treatment == 1])
			p0 <- mean(response[treatment == 0])
			variance <- c(p1 * (1 - p1), p0 * (1 - p0))
			ratio <- rev(table(treatment))
			X <- p1 - p0
			Z <- (X - theta0)/sqrt(sum(variance/ratio))
			P <- 1 - pnorm(Z)
		}
		else stop("Invalid probability model for monitoring")
	}
	else if(prob.model == "O") {
# analyze data for binomial odds
		if(!log.transform) theta0 <- log(theta0)
		if(arms == 1) {
			p1 <- mean(response)
			variance <- 1/(p1 * (1 - p1))
			test <- summary(glm(response ~ -1 + rep(1, length(
				response)), family = "binomial"))
			X <- test$coeff[1, 1]
			Z <- (test$coeff[1, 1] - theta0)/test$coeff[1, 2]
			P <- 1 - pnorm(Z)
			ratio <- 1
		}
		else if(arms == 2) {
			p1 <- mean(response[treatment == 1])
			p0 <- mean(response[treatment == 0])
			variance <- 1/c(p1 * (1 - p1), p0 * (1 - p0))
			ratio <- rev(table(treatment))
			test <- summary(glm(response ~ treatment, family = 
				"binomial"))
			X <- test$coeff[2, 1]
			Z <- (test$coeff[2, 1] - theta0)/test$coeff[2, 2]
			P <- 1 - pnorm(Z)
		}
		else if(arms == 0) {
			ratio <- treatment
			test <- summary(glm(response ~ treatment, family = 
				"binomial"))
			variance <- test$coeff[2, 2]^2 * length(response)
			X <- test$coeff[2, 1]
			Z <- (test$coeff[2, 1] - theta0)/test$coeff[2, 2]
			P <- 1 - pnorm(Z)
		}
		else stop("Invalid probability model for monitoring")
		if(!log.transform) {
			X <- exp(X)
			theta0 <- exp(theta0)
		}
	}
	else if(prob.model == "P") {
# analyze data for Poisson
		if(!log.transform) theta0 <- log(theta0)
		if(arms == 1) {
			m1 <- mean(response)
			variance <- 1/m1
			test <- summary(glm(response ~ -1 + rep(1, length(
				response)), family = "poisson"))
			X <- test$coeff[1, 1]
			Z <- (test$coeff[1, 1] - theta0)/test$coeff[1, 2]
			P <- 1 - pnorm(Z)
			ratio <- 1
		}
		else if(arms == 2) {
			m1 <- mean(response[treatment == 1])
			m0 <- mean(response[treatment == 0])
			variance <- 1/c(m1, m0)
			ratio <- rev(table(treatment))
			test <- summary(glm(response ~ treatment, family = 
				"poisson"))
			X <- test$coeff[2, 1]
			Z <- (test$coeff[2, 1] - theta0)/test$coeff[2, 2]
			P <- 1 - pnorm(Z)
		}
		else if(arms == 0) {
			ratio <- treatment
			test <- summary(glm(response ~ treatment, family = 
				"poisson"))
			variance <- test$coeff[2, 2]^2 * length(response)
			X <- test$coeff[2, 1]
			Z <- (test$coeff[2, 1] - theta0)/test$coeff[2, 2]
			P <- 1 - pnorm(Z)
		}
		else stop("Invalid probability model for monitoring")
		if(!log.transform) {
			X <- exp(X)
			theta0 <- exp(theta0)
		}
	}
	else if(prob.model == "H") {
# analyze data for proportional hazards
		if(!is.Surv(response)) stop(
				"Sequential proportional hazards model response must be Surv object"
				)
		if(!is.element(attr(response, "type"), c("right", "counting"))
			)
			stop("response must be of type right or counting")
		if(arms == 2) {
			if(!log.transform)
				theta0 <- log(theta0)
			test <- coxph(response ~ treatment, init = theta0)
			X <- test$coeff[1]
			Z <- sqrt(test$score) * ifelse1(X > 0, 1, -1)
			P <- 1 - pnorm(Z)
			ratio <- rev(table(treatment))
			variance <- 0.25
			if(!log.transform) {
				X <- exp(X)
				theta0 <- exp(theta0)
			}
		}
		else stop("Invalid probability model for monitoring")
	}
	else stop("Invalid probability model for monitoring")
	sigma.sqr <- ifelse1(arms != 2 || prob.model == "H", variance, sum(
		variance/ratio) * ratio[2])	#
# Preliminary processing of design constraints
	nbr.analyses <- xPar$nbr.analyses
	minConstraint <- xDes$specification$minimum.constraint
	maxConstraint <- xDes$specification$maximum.constraint
	exactConstraint <- xDes$specification$exact.constraint
	haveExact <- !is.null(exactConstraint)
	haveMin <- !is.null(minConstraint)
	haveMax <- !is.null(maxConstraint)
	if(haveExact) {
# Fill in values in min and max, to avoid NA's when interpolating
		if(haveMin) minConstraint[is.na(minConstraint)] <- 
				exactConstraint[is.na(minConstraint)]
		if(haveMax)
			maxConstraint[is.na(maxConstraint)] <- exactConstraint[
				is.na(maxConstraint)]
	}
	interpolateConstraint <- function(lambda, y)
	{
# Return y[lambda], except interpolate if lambda is non-integer
		lambda.i <- floor(lambda)
		lambda.r <- lambda - lambda.i
		z <- y[lambda.i]
		w <- lambda.r > 0
		z[w] <- z[w] + lambda.r[w] * (y[lambda.i[w] + 1] - y[lambda.i[
			w]])
		z
	}
	minimum.constraint <- NULL
	maximum.constraint <- NULL	# overwrite later if haveMin or haveMax
#
# current analysis time
	Nj <- ifelse1(prob.model == "H", sum(response[, "status"]), length(
		response))	#
# append current statistics to previous analysis results
	N <- c(xStat$N, Nj)
	X <- c(xStat$X, X)
	Z <- c(xStat$Z, Z)
	P <- c(xStat$P, P)
	variance <- c(xStat$variance, list(variance))
	ratio <- c(xStat$ratio, list(ratio))
	response <- c(xStat$response, list(response))
	treatment <- c(xStat$treatment, list(treatment))
	j <- length(N)	#
# check future.analyses, delete times too close after present
	if(missing(future.analyses))
		future.analyses <- ifelse1(is.seqMonitor(x), x$monitor$
			future.analyses, xPar$sample.size/max(xPar$sample.size
			))
	if(length(future.analyses)) {
		if(max(future.analyses, na.rm = T) < 1)
			future.analyses <- c(future.analyses, 1)
		increments <- diff(future.analyses[!is.na(future.analyses)])
		if(length(increments) && any(increments <= 0))
			stop("future.analyses must be increasing values")
		N.max <- max(future.analyses, na.rm = T)
		N.min <- ifelse1(N.max > 1, Nj + min.increment * N.max, Nj/max(
			xPar$sample.size) + min.increment)
		bad <- !is.na(future.analyses) & (future.analyses < N.min)
		if(any(bad)) {
			if(!missing(future.analyses)) {
# Stop or give a warning about deleted observations
				if(any(future.analyses > 1, na.rm = T) && bad[
				  length(bad)] && future.analyses[length(bad)] !=
				  Nj) stop(
				    "Specified maximum sample size must fall after current time + min.increment%"
				    )
				warning(paste(sum(bad), 
				  "specified future analysis time(s) were within min.increment of the current time or earlier, and are deleted"
				  ))
			}
			future.analyses <- future.analyses[!bad]
		}
	}
	UseFutureAnalyses <- (length(future.analyses) > 0)
	if(!UseFutureAnalyses) {
		constrainNJ <- T
		NJ <- Nj
	}
	else {
#
# define constrainNJ: if the last element of future.analyses is NA or
# if all elements of future.analyses are less than or equal to 1, 
# then the maximal sample size is to be adjusted to maintain
# the desired power to detect the alternative. Otherwise, 
# the maximal sample size is to be held constant at the value
# of the last element in future.analyses.
# NJ is a trial value for the maximal sample size
# If NJ is unconstrained, we take the first guess as the
# value estimated at the last analysis (or from the
# design if this is the first analysis) modified to
# allow for errors in estimating sigma.sqr
		NJ <- future.analyses[length(future.analyses)]
		constrainNJ <- !(is.na(NJ) || NJ == 1)
		if(constrainNJ)
			NJ <- NJ
		else {
			NJ.min <- ceiling(Nj/(1 - min.increment))
			NJ <- max(NJ.min, ceiling((max(xPar$sample.size) * 
				sigma.sqr)/xMod$sigma.sqr))
		}
	}
	interpolate <- function(x, y, xout)
	{
# local function to interpolate error spending function, etc.
		wna <- which.na(y)
		if(length(wna)) {
			x <- x[ - wna]
			y <- y[ - wna]
		}
		approx(x, y, xout)$y
	}
	replaceNAs <- function(x)
	{
# function to interpolate missing values in sample.size
		wna <- which.na(x)
		if(!length(wna))
			return(x)
		J <- length(x)
		approx((1:J)[ - wna], x[ - wna], 1:J)$y
	}
#
# iterate until monitoring design is found
	lowNJ <- NULL
	highNJ <- NULL
	for(iteration in 1:(maxiter + 1)) {
		if(iteration > maxiter) {
			warning("Failed to converge in seqMonitor")
			if(!is.null(highNJ))
				newMon <- previousHigh
			break
		}
		if(UseFutureAnalyses) {
# compute sample.size; delete elements if any increments too small
			N2 <- ifelse1(any(future.analyses <= 1, na.rm = T), 
				future.analyses * NJ, future.analyses)
			if(is.na(N2[length(N2)]))
				N2[length(N2)] <- NJ
			abs.increment <- NJ * min.increment
			J <- length(N2)
			violatesIncrement <- (diff(replaceNAs(c(Nj, N2))) < 
				abs.increment)
			while(J > 1 && any(violatesIncrement)) {
# delete first time that violates the abs.increment
				N2 <- N2[ - min(J - 1, ((1:J)[
				  violatesIncrement])[1])]
				J <- length(N2)
				violatesIncrement <- (diff(replaceNAs(c(Nj, N2
				  ))) < abs.increment)
			}
#
# Replace missing values, 
# Round up to nearest integers, except possibly if constrainNJ
			sample.size <- replaceNAs(c(N, N2))
			if(!constrainNJ || all(N2 == round(N2), na.rm = T))
				sample.size <- ceiling(sample.size - 1e-005)
		}
		else sample.size <- N
		J <- length(sample.size)	#
# fill in constraints from earlier analyses
		exact <- seqBoundary(as.numeric(rep(NA, J * 4)), 
			constraint.scale)
		if(j > 1) {
			if(constraint.scale == "X")
				exact[1:(j - 1),  ] <- xMon$Xbnd
			else if(constraint.scale == "Z")
				exact[1:(j - 1),  ] <- xMon$Zbnd
			else if(constraint.scale == "P")
				exact[1:(j - 1),  ] <- 1 - pnorm(xMon$Zbnd)
			else if(constraint.scale == "E")
				exact[1:(j - 1),  ] <- xMon$Ebnd
		}
#
# Fill in design constraints
		if(nbr.analyses == 1) {
			lambda <- 1
		}
		else {
			lambda <- approx(xPar$sample.size/max(xPar$sample.size
				), 1:nbr.analyses, xout = sample.size/NJ, rule
				 = 2)$y
			if(any(close <- abs(lambda - round(lambda)) < 0.01))
				lambda[close] <- round(lambda[close])
		}
		if(haveExact)
			for(i in 1:4)
				exact[j:J, i] <- interpolateConstraint(lambda, 
				  exactConstraint[, i])
		if(haveMin) {
			minimum.constraint <- NA * exact
			for(i in 1:4)
				minimum.constraint[j:J, i] <- 
				  interpolateConstraint(lambda, minConstraint[,
				  i])
		}
		if(haveMax) {
			maximum.constraint <- NA * exact
			for(i in 1:4)
				maximum.constraint[j:J, i] <- 
				  interpolateConstraint(lambda, maxConstraint[,
				  i])
		}
		if(constraint.scale == "E") {
# handle error spending constraint.scale
			PIcum <- sample.size/NJ
			error.spend <- xDes$error.spend
			dPIcum <- c(0, xDes$parameters$sample.size)
			dPIcum <- dPIcum/max(dPIcum)
			for(col in 1:4) {
				if(xDes$parameters$design.family == "E") {
				  exact[j:J, col] <- ifelse1(xDes$parameters$P[
				    col] != 0, PIcum[j:J]^( - xDes$parameters$
				    P[col]), 1 - (1 - PIcum[j:J])^xDes$
				    parameters$R[col])
				}
				else exact[j:J, col] <- interpolate(dPIcum, c(
				    0, error.spend[, col]), PIcum[j:J])
				for(row in 2:J)
				  exact[row, col] <- max(exact[row, col], 
				    exact[row - 1, col])
			}
			exact[J, !is.na(exact[J,  ])] <- 1
			exact[!is.na(exact) & exact < 0] <- 0
			newMon <- update(x$design, design.family = "E", 
				exact.constraint = exact, P = , R = , A = , 
				sample.size = sample.size, nbr.analyses = 
				length(sample.size), alt.hypothesis = , power
				 = power, ratio = ratio[[j]], variance = 
				variance[[j]], minimum.constraint = 
				minimum.constraint, maximum.constraint = 
				maximum.constraint, byValue = T)
		}
		else {
# handle all other constraint.scales
			newMon <- update(x$design, exact.constraint = exact, 
				sample.size = sample.size, nbr.analyses = 
				length(sample.size), alt.hypothesis = , power
				 = power, ratio = ratio[[j]], variance = 
				variance[[j]], minimum.constraint = 
				minimum.constraint, maximum.constraint = 
				maximum.constraint, byValue = T)
		}
		if(constrainNJ) break	#
# assign current results to previousHigh or previousLow
		largeEnough <- ifelse1(upper, newMon$hypothesis$theta.alt <= 
			theta.alt, newMon$hypothesis$theta.alt >= theta.alt)
		if(largeEnough) {
			if(!UseFutureAnalyses) {
# power is adequate with no future analyses
				break
			}
			if(NJ == NJ.min) {
# NJ is is already as low as possible given min.increment;
# see if power is adequate with no future analyses
				UseFutureAnalyses <- F
			}
			highNJ <- NJ
			previousHigh <- newMon
		}
		else {
			if(!UseFutureAnalyses) {
# Power adequate with NJ at minimum, and need at least one 
# future analysis; stop, use results from NJ at minimum
				newMon <- previousHigh
				break
			}
			lowNJ <- NJ
			previousLow <- newMon
		}
		if(!is.null(lowNJ) && !is.null(highNJ) && highNJ - lowNJ <= 1
			) {
			newMon <- previousHigh
			break
		}
#
# set next NJ
		if(is.null(lowNJ))
			NJ <- max(NJ.min, min(NJ - 1, floor(0.98 * NJ)))
		else if(is.null(highNJ))
			NJ <- ceiling(NJ * (1.01 + iteration/64))
		else NJ <- min(highNJ - 1, max(lowNJ + 1, round(interpolate(c(
				previousLow$hypothesis$theta.alt, previousHigh$
				hypothesis$theta.alt), c(lowNJ, highNJ), 
				theta.alt))))
	}
# end of for(iteration) loop.
	Xbnd <- ifelse1(attr(newMon$boundary, "scale") == "X", newMon$boundary,
		changeSeqScale(newMon, "X"))
	Zbnd <- ifelse1(attr(newMon$boundary, "scale") == "Z", newMon$boundary,
		changeSeqScale(newMon, "Z"))
	Ebnd <- newMon$error.spend
	newMon <- c(newMon, list(seqDesignCall = as.call(c(as.name("seqDesign"
		), newMon$specification)), design = xDes, statistics = list(N
		 = N, variance = variance, X = X, Z = Z, P = P, response = 
		response, treatment = treatment, ratio = ratio), monitor = 
		list(constraint.scale = constraint.scale, future.analyses = 
		future.analyses, Xbnd = rbind(xMon$Xbnd, Xbnd[j,  ]), Zbnd = 
		rbind(xMon$Zbnd, Zbnd[j,  ]), Ebnd = rbind(xMon$Ebnd, Ebnd[j,  
		]))))
	newMon$call <- call
	if(!constrainNJ && length(newMon$parameters$sample.size) == length(N)
		)
		cat("\n*** No future analyses are needed to attain desired power. ***\n"
			)
	else if(!missing(future.analyses) && length(newMon$parameters$
		sample.size) < length(c(N, future.analyses)))
		warning("One or more analysis times were deleted because of the min.increment requirement"
			)
	BNDj <- newMon$monitor$Xbnd[j,  ]
	Eps <- xPar$epsilon
	action <- c("Stop", "Continue", "Stop", "Continue", "Stop")
	decision <- c("Lower Alternative Hypothesis", "", "Null Hypothesis", 
		"", "Upper Alternative Hypothesis")
	if(Eps[1] == 0) {
		action[3] <- action[2]
		decision <- decision[c(3, 2, 2, 2, 5)]
	}
	else if(Eps[2] == 0) {
		action[3] <- action[2]
		decision <- decision[c(1, 2, 2, 2, 3)]
	}
	else if(abs(1 - sum(Eps)) <= .Machine$single.eps) {
		action[3] <- action[2]
		decision <- decision[c(1, 2, 2, 2, 3)]
	}
	indx <- sum(X[j] > BNDj) + 1
	decision <- decision[indx]
	action <- action[indx]
	newMon$monitor <- c(newMon$monitor, list(action = action, decision = 
		decision))
	class(newMon) <- c("seqMonitor", "seqDesign")
	newMon
}
# Changes:
# 2/28/2000 Case with alternative hypothesis < null hypothesis 
#           will now converge.


summary.seqDesign <- function(object,  theta=NULL, fixed=F,
			      ordering="m")
{
  if(is.null(theta)){
    # Create default theta, 5 values.
    theta.null <- seqExtract(object, "theta.null")
    theta.alt <- seqExtract(object, "theta.alt")
    ThetaisExp <- ifelse1(is.element(object$model$prob.model, c("N","B")),
			  0, !object$model$log.transform)
    # use geometric progression if ThetaisExp
    # For 2-sided case put theta.null in the middle.
    end1 <- ifelse1(is.element(object$hypothesis$test.type,
			       c("less", "greater")),
		    theta.null,
		    ifelse1(ThetaisExp,
			    theta.null^2 / theta.alt,
			    2 * theta.null - theta.alt))
    theta <- sort(ifelse1(ThetaisExp,
			  exp(seq(log(end1), log(theta.alt), length=5)),
			  seq(end1, theta.alt, length=5)))
  }
  if(fixed){
    fixed <- seqFixDesign(object)
    fixed.oc <- seqOperatingChar(fixed, theta=theta)
  }
  else {
    fixed <- NULL
    fixed.oc <- NULL
  }

  obj <- list(operating.char = seqOperatingChar(object, theta=theta),
	      fixed.operating.char = fixed.oc,
	      inference = seqInference(object, ordering=ordering),
	      design = object,
	      fixed = fixed,
	      theta = theta)
  class(obj) <- "summary.seqDesign"
  obj
}
# Changes:
# 3/29/2000 Default theta: 
#           make null hypothesis the middle value for two-sided design,
#           use geometric progression when appropriate.



# End loading here if you do not want objects moved automatically.
#--------------------------------------------------------------------

if(.SeqTrial.where != search()[1]){
  # Save objects to the location of the SeqTrial module
  for(i in .SeqTrial.objectsInPatch)
    assign(i, get(i), where=.SeqTrial.where)
  remove(.SeqTrial.objectsInPatch)
}
rm(.SeqTrial.objectsInPatch, .SeqTrial.where)