R/gMCP.R

In gMCP: Graph Based Multiple Comparison Procedures

#' Graph based Multiple Comparison Procedures
#' 
#' Performs a graph based multiple test procedure for a given graph and
#' unadjusted p-values.
#' 
#' For the Bonferroni procedure the p-values can arise from any statistical
#' test, but if you improve the test by specifying a correlation matrix, the
#' following assumptions apply:
#' 
#' It is assumed that under the global null hypothesis
#' \eqn{(\Phi^{-1}(1-p_1),...,\Phi^{-1}(1-p_m))} follow a multivariate normal
#' distribution with correlation matrix \code{correlation} where
#' \eqn{\Phi^{-1}} denotes the inverse of the standard normal distribution
#' function.
#' 
#' For example, this is the case if \eqn{p_1,..., p_m} are the raw p-values
#' from one-sided z-tests for each of the elementary hypotheses where the
#' correlation between z-test statistics is generated by an overlap in the
#' observations (e.g. comparison with a common control, group-sequential
#' analyses etc.). An application of the transformation \eqn{\Phi^{-1}(1-p_i)}
#' to raw p-values from a two-sided test will not in general lead to a
#' multivariate normal distribution. Partial knowledge of the correlation
#' matrix is supported. The correlation matrix has to be passed as a numeric
#' matrix with elements of the form: \eqn{correlation[i,i] = 1} for diagonal
#' elements, \eqn{correlation[i,j] = \rho_{ij}}, where \eqn{\rho_{ij}} is the
#' known value of the correlation between \eqn{\Phi^{-1}(1-p_i)} and
#' \eqn{\Phi^{-1}(1-p_j)} or \code{NA} if the corresponding correlation is
#' unknown. For example correlation[1,2]=0 indicates that the first and second
#' test statistic are uncorrelated, whereas correlation[2,3] = NA means that
#' the true correlation between statistics two and three is unknown and may
#' take values between -1 and 1. The correlation has to be specified for
#' complete blocks (ie.: if cor(i,j), and cor(i,j') for i!=j!=j' are specified
#' then cor(j,j') has to be specified as well) otherwise the corresponding
#' intersection null hypotheses tests are not uniquely defined and an error is
#' returned.
#' 
#' For further details see the given references.
#' 
#' @param graph A graph of class \code{\link{graphMCP}}.
#' @param pvalues A numeric vector specifying the p-values for the graph based
#' MCP. Note the assumptions in the details section for the parametric tests, 
#' when a correlation is specified.
#' @param test Should be either \code{"Bonferroni"}, \code{"Simes"} or \code{"parametric"}.
#' If not specified by default the Bonferroni-based test procedure is used if no
#' correlation is specified or the algorithm from Bretz et al. 2011 if a
#' correlation is specified. If \code{test} is set to \code{"Simes"} the weighted
#' Simes test will be performed for each subset of hypotheses.
#' @param correlation Optional correlation matrix.  If the weighted Simes test
#' is performed, it is checked whether type I error rate can be ensured and a
#' warning is given if this is not the case.  For parametric tests the p-values
#' must arise from one-sided tests with multivariate normal distributed test
#' statistics for which the correlation is (partially) known. In that case a
#' weighted parametric closed test is performed (also see
#' \code{\link{generatePvals}}). Unknown values can be set to NA. (See details
#' for more information)
#' @param alpha A numeric specifying the maximal allowed type one error rate.
#' @param approxEps A boolean specifying whether epsilon values should be
#' substituted with the value given in the parameter \code{eps}.
#' @param eps A numeric scalar specifying a value for epsilon edges.
#' @param ...  Test specific arguments can be given here.
#' @param upscale Logical. If \code{upscale=FALSE} then for each intersection 
#' of hypotheses (i.e. each subgraph) a weighted test is performed at the 
#' possibly reduced level alpha of sum(w)*alpha, 
#' where sum(w) is the sum of all node weights in this subset.
#' If \code{upscale=TRUE} all weights are upscaled, so that sum(w)=1.
#' 
#' For backward comptibility the default value is TRUE if a the parameter \code{test}
#' is missing, but parameter \code{correlation} is specified or if \code{test=="Bretz2011"}.
#' @param useC Logical scalar. If \code{TRUE} neither adjusted p-values nor
#' intermediate graphs are returned, but the calculation is sped up by using
#' code written in C. THIS CODE IS NOT FOR PRODUCTIVE USE YET!  If approxEps is
#' \code{FALSE} and the graph contains epsilon edges, a warning is thrown and
#' \code{useC} will be ignored.
#' @param verbose Logical scalar. If \code{TRUE} verbose output is generated.
#' @param keepWeights Logical scalar. If \code{FALSE} the weight of a node
#' without outgoing edges is set to 0 if it is removed.  Otherwise it keeps its
#' weight.
#' @param adjPValues Logical scalar. If \code{FALSE} no adjusted p-values will
#' be calculated.  Especially for the weighted Simes test this will result in
#' significantly less calculations in most cases.
#' @return An object of class \code{gMCPResult}, more specifically a list with
#' elements
#' \describe{
#' \item{\code{graphs}}{list of graphs}
#' \item{\code{pvalues}}{p-values}
#' \item{\code{rejected}}{logical whether hyptheses could be rejected}
#' \item{\code{adjPValues}}{adjusted p-values}
#' }
#' @author Kornelius Rohmeyer \email{rohmeyer@@small-projects.de}
#' @seealso \code{\link{graphMCP}} \code{\link[multcomp:contrMat]{graphNEL}}
#' @references Frank Bretz, Willi Maurer, Werner Brannath, Martin Posch: A
#' graphical approach to sequentially rejective multiple test procedures.
#' Statistics in Medicine 2009 vol. 28 issue 4 page 586-604.
#' \url{https://www.meduniwien.ac.at/fwf_adaptive/papers/bretz_2009_22.pdf}
#' 
#' Bretz F., Posch M., Glimm E., Klinglmueller F., Maurer W., Rohmeyer K.
#' (2011): Graphical approaches for multiple endpoint problems using weighted
#' Bonferroni, Simes or parametric tests. Biometrical Journal 53 (6), pages 894-913, Wiley.
#' \doi{10.1002/bimj.201000239}
#' 
#' Strassburger K., Bretz F.: Compatible simultaneous lower confidence bounds
#' for the Holm procedure and other Bonferroni based closed tests. Statistics
#' in Medicine 2008; 27:4914-4927.
#' 
#' Hommel G., Bretz F., Maurer W.: Powerful short-cuts for multiple testing
#' procedures with special reference to gatekeeping strategies. Statistics in
#' Medicine 2007; 26:4063-4073.
#' 
#' Guilbaud O.: Simultaneous confidence regions corresponding to Holm's
#' stepdown procedure and other closed-testing procedures. Biometrical Journal
#' 2008; 50:678-692.
#' @keywords htest graphs
#' @examples
#' 
#' 
#' g <- BonferroniHolm(5)
#' gMCP(g, pvalues=c(0.01, 0.02, 0.04, 0.04, 0.7))
#' # Simple Bonferroni with empty graph:
#' g2 <- matrix2graph(matrix(0, nrow=5, ncol=5))
#' gMCP(g2, pvalues=c(0.01, 0.02, 0.04, 0.04, 0.7))
#' # With 'upscale=TRUE' equal to BonferroniHolm:
#' gMCP(g2, pvalues=c(0.01, 0.02, 0.04, 0.04, 0.7), upscale=TRUE)
#' 
#' # Entangled graphs:
#' g3 <- Entangled2Maurer2012()
#' gMCP(g3, pvalues=c(0.01, 0.02, 0.04, 0.04, 0.7), correlation=diag(5))
#' 
#' @export gMCP

gMCP <- function(graph, pvalues, test, correlation, alpha=0.05, 
		approxEps=TRUE, eps=10^(-3), ..., upscale=ifelse(missing(test)&&!missing(correlation)||!missing(test)&&test=="Bretz2011", TRUE, FALSE),
    useC=FALSE, verbose=FALSE, keepWeights=FALSE, adjPValues=TRUE) {
#		, alternatives="less") {	
  # Temporary translation of old syntax:
  if (!missing(test)) {
    if (test=="Bretz2011") {
      test <- "parametric"
      if (!upscale) warning("For 'Bretz2011' upscale is 'TRUE'.")
      upscale <- TRUE
    } else if (test=="simple-parametric") {
      test <- "parametric"
      if (upscale) warning("For 'simple-parametric' upscale is 'FALSE'.")
      upscale <- FALSE
    }
  }
  # Check for entangled graphs and yet unsupported tests:
  if (!missing(test) && test=="Simes" && "entangledMCP" %in% class(graph)) {
    #TODO
    stop("Simes test is not yet supported for entangled graphs in this version.")
  }
  
  if (approxEps) {
    graph <- substituteEps(graph, eps=eps)
  }
  if ("graphMCP" %in% class(graph) && !is.numeric(graph@m)) {
    stop("Graph seems to contain variables - please use function replaceVariables.")
    #graph <- parse2numeric(graph) # TODO ask for variables
  }
  
  if ((missing(test) && missing(correlation) || !missing(test) && test=="Bonferroni") && "entangledMCP" %in% class(graph)) {		
		out <- graphTest(pvalues=pvalues, weights=getWeights(graph), alpha=alpha*graph@weights, G=getMatrices(graph))
		result <- new("gMCPResult", graphs=list(graph), alpha=alpha, pvalues=pvalues, rejected=(out==1), adjPValues=numeric(0))
    attr(result, "call") <- call2char(match.call())
		return(result)
	}
	output <- ""
	callFromGUI <- !is.null(list(...)[["callFromGUI"]])
  #TODO - Disable checkOptimal as long as it takes > 1 min for some example graphs.
	#if (verbose) {
	#	output <- paste(output, checkOptimal(graph, verbose=FALSE), sep="\n")
	#}
	sequence <- list(graph)
	if (length(pvalues)!=length(getNodes(graph))) {
		stop("Length of pvalues must equal number of nodes.")
	}
	if (is.null(names(pvalues))) {
		names(pvalues) <- getNodes(graph)
	}
	if ((missing(test) && missing(correlation)) || !missing(test) && test == "Bonferroni") {
		if (!missing(correlation)) stop("Bonferroni test can not take correlation into account. Please specify test procedure.")
		# Bonferroni-based test procedure		
		if (useC && !upscale) {
			w <- getWeights(graph)			
			result <- fastgMCP(m=graph2matrix(graph), w=w, p=pvalues, a=alpha, keepWeights=keepWeights)
			row.names(result$m) <- getNodes(graph)
			lGraph <- matrix2graph(result$m)
			lGraph <- setWeights(lGraph, result$w)			
			lGraph <- setRejected(lGraph, getNodes(lGraph), result$rejected)
			sequence <- c(sequence, lGraph)
			result <- new("gMCPResult", graphs=sequence, alpha=alpha, pvalues=pvalues, rejected=getRejected(lGraph), adjPValues=numeric(0))
			attr(result, "call") <- call2char(match.call())
			return(result)
		} else {
		  if (useC) warning("Parameter UseC==TRUE is not supported for upscale==TRUE.")
			while(!is.null(node <- getRejectableNode(graph, alpha, pvalues))) {
				# if (verbose) cat(paste("Node \"",node,"\" can be rejected.\n",sep=""))
				graph <- rejectNode(graph, node, verbose, keepWeights=keepWeights, upscale=upscale)
				sequence <- c(sequence, graph)
			}	
			adjPValues <- adjPValues(sequence[[1]], pvalues, upscale=upscale, verbose=verbose)@adjPValues
			result <- new("gMCPResult", graphs=sequence, alpha=alpha, pvalues=pvalues, rejected=getRejected(graph), adjPValues=adjPValues)
			attr(result, "call") <- call2char(match.call())
			return(result)
		}
	} else if ((missing(test) && !missing(correlation)) || !missing(test) && test == "parametric") {				
		if (missing(correlation) || !is.matrix(correlation)) {
			stop("Procedure for correlated tests, expects a correlation matrix as parameter \"correlation\".")
		} else {
      check <- checkCorrelation(correlation, TRUE)
      if (!isTRUE(check)) {
        stop(paste("Parameter 'correlation' is no correlation matrix:", check, sep="\n"))
      }
#			if (is.character(correlation)) {
#				samplesize <- list(...)[["samplesize"]]
#				if (is.null(samplesize)) samplesize <- getBalancedDesign(correlation, length(pvalues))				
#				x <- contrMat(samplesize, type = correlation) # balanced design up to now and only Dunnett will work with n+1
#				var <- x %*% diag(length(samplesize)) %*% t(x)
#				correlation <- diag(1/sqrt(diag(var)))%*%var%*%diag(1/sqrt(diag(var)))
#			}                       			
			w <- getWeights(graph)
			if( all(w==0) ) {
			  adjP <- rep(1,length(w))
			  rejected <- rep(FALSE,length(w))
			  names(rejected) <- getNodes(graph)
			  names(adjP) <- getNodes(graph)
			} else {
			  #myTest <- generateTest(Gm, w, correlation, alpha)
			  #zScores <- -qnorm(pvalues)
			  #rejected <- myTest(zScores)
			  adjP <- generatePvals(g=graph, cr=correlation, p=pvalues, upscale=upscale) #, alternatives=alternatives)
			  rejected <- adjP <= alpha
			  names(adjP) <- getNodes(graph)
			  names(rejected) <- getNodes(graph)
			}
			result <- new("gMCPResult", graphs=sequence, alpha=alpha, pvalues=pvalues, rejected=rejected, adjPValues=adjP)
			attr(result, "call") <- call2char(match.call())
			return(result)
		}
	} else if (!missing(test) && test=="Simes") {		
		m <- graph2matrix(graph)
		if (!adjPValues) {
			if (all(pvalues>alpha)) {
				result <- new("gMCPResult", graphs=sequence, alpha=alpha, pvalues=pvalues, rejected=getRejected(graph))
				if (verbose) {
					output <- paste(output, "All p-values above alpha.", sep="\n")
					if (!callFromGUI) cat(output,"\n")
					attr(result, "output") <- output
				}
				attr(result, "call") <- call2char(match.call())
				return(result)
			}
			if (all(pvalues<=alpha)) {			
				rejected <- rep(TRUE, dim(m)[1])
				names(rejected) <- getNodes(graph)
				result <- new("gMCPResult", graphs=sequence, alpha=alpha, pvalues=pvalues, rejected=rejected)
				if (verbose) {
					output <- paste(output, "All p-values below alpha.", sep="\n")
					if (!callFromGUI) cat(output,"\n")
					attr(result, "output") <- output
				}
				attr(result, "call") <- call2char(match.call())
				return(result)
			}
			while(!is.null(node <- getRejectableNode(graph, alpha, pvalues))) {
				output <- paste(output, paste("Hypothesis \"",node,"\" can be rejected by the weighted Bonferroni based test and therefore by weighted Simes test.\n",sep=""), sep="\n")
				graph <- rejectNode(graph, node, upscale=upscale, verbose=verbose)
				sequence <- c(sequence, graph)
			}
		}
		n <- sum(!getRejected(graph))
		if (n<3 && !adjPValues) {
			result <- new("gMCPResult", graphs=sequence, alpha=alpha, pvalues=pvalues, rejected=getRejected(graph))
			if (verbose) {
				if (n==2) output <- paste(output, "Only two hypotheses remaining.", sep="\n")
				if (n==1) output <- paste(output, "Only one hypothesis remaining.", sep="\n")
				if (n==0) output <- paste(output, "Everything already rejected.", sep="\n")
				if (!callFromGUI) cat(output,"\n")
				attr(result, "output") <- output
			}
			attr(result, "call") <- call2char(match.call())
			return(result)
		} else {			
			graph2 <- subgraph(graph, !getRejected(graph))
			output <- paste(output, "Remaining hypotheses (new numbering):", paste(1:length(getNodes(graph2)),": ", getNodes(graph2), sep="",collapse="\n"), sep="\n")
			pvalues2 <- pvalues[!getRejected(graph)]
			allSubsets <- permutations(length(getNodes(graph2)))[-1,]
			result <- cbind(allSubsets, 0, Inf)
			weights <- generateWeights(graph2@m, getWeights(graph2))[,(n+(1:n))]
			if (verbose) explanation <- rep("not rejected", dim(allSubsets)[1])
			for (i in 1:dim(allSubsets)[1]) {
				subset <- allSubsets[i,]
				if(!all(subset==0)) {
					J <- which(subset!=0)	
					if (verbose) explanation[i] <- paste("Subset {",paste(J,collapse=","),"}: ",explanation[i], sep="")			
					mJ <- Inf					
					for (j in J) {
						Jj <- subset!=0 & (pvalues2 <= pvalues2[j]) # & (1:n)!=j
						if (adjPValues) {
							mJt <- pvalues2[j]/sum(weights[i, Jj])	
							if (is.na(mJt)) { # this happens only if pvalues2[j] is 0
								mJt <- 0
							}
							#cat("pvalues2:\n", pvalues2, "\nmJt: ", mJt, "\nmJ: ", mJ, "\nj: ", j, "\nJ: ", J, "\nsubset: ", subset, "\n")
							if (mJt<mJ) {
								mJ <- mJt
							}
						}
						#cat("j: ",j, ", Jj: ",Jj,"\n")
						#cat("p_",j,"=",pvalues2[j],"<=a*(w_",paste(which(Jj),collapse ="+w_"),")=",alpha,"*(",paste(weights[i, Jj],collapse ="+"),")=",sum(weights[i, Jj]),"\n")
						if (pvalues2[j]<=alpha*sum(weights[i, Jj])) {
							result[i, n+1] <- 1
							if (verbose) {
								explanation[i] <- paste("Subset {",paste(J,collapse=","),"}",ifelse(adjPValues,paste(" (padj: ",mJ,")",sep=""),""),": p_",j,"=", pvalues2[j],"<=a*(w_",paste(which(Jj),collapse ="+w_"),")\n     =",alpha,"*(",paste(weights[i, Jj],collapse ="+"),")=",alpha*sum(weights[i, Jj]),sep="")
							}
						}
					}	
					result[i, n+2] <- mJ
				} 
			}
			adjPValuesV <- rep(NA, n)
			for (i in 1:n) {
				if (all(result[result[,i]==1,n+1]==1)) {
					graph <- rejectNode(graph, getNodes(graph2)[i], upscale=upscale)
				}
				adjPValuesV[i] <- max(result[result[,i]==1,n+2])
			}
			result <- new("gMCPResult", graphs=sequence, alpha=alpha, pvalues=pvalues, rejected=getRejected(graph), adjPValues=adjPValuesV)
			if (verbose) {
				output <- paste(output, paste(explanation, collapse="\n"), sep="\n")
				if (!callFromGUI) cat(output,"\n")
				attr(result, "output") <- output
			}
			attr(result, "call") <- call2char(match.call())
			return(result)
		}
	} else {
		stop(paste("Specified test \"",test,"\" is not known.",sep=""))
	}
}

call2char <- function(call) {
  paste(capture.output(print(call)), collapse="\n")
}

createGMCPCall <- function(graph, pvalues, test, correlation, alpha=0.05, 
		approxEps=TRUE, eps=10^(-3), ..., upscale=FALSE, useC=FALSE, 
		verbose=FALSE, keepWeights=TRUE, adjPValues=TRUE) {	
	command <- dputGraph(graph, "graph")
	command <- paste(command, "pvalues <- ",dput2(unname(pvalues)),"\n", sep="")
	if (!missing(correlation)) {
		command <- paste(command, dputMatrix(correlation, name="cr", indent=12),"\n", sep="")
	}
	command <- paste(command, "gMCP(graph, pvalues", sep="")
	if (!missing(test)) {
		command <- paste(command, ", test=\"",test,"\"", sep="")
	}
	if (upscale) {
	  command <- paste(command, ", upscale=TRUE", sep="")
	}	
	if (!missing(correlation)) {
		command <- paste(command, ", correlation=cr", sep="")
	}
	command <- paste(command, ", alpha=",dput2(alpha), sep="")
	command <- paste(command, ")\n", sep="")
	return(command)
}

#TODO: Set rejected.
dputGraph <- function(g, name="graph") {
	# Entangled graphs and recursive calls:
	if ("entangledMCP" %in% class(g)) {
		s <- c()
		i <- 1
		for (graph in g@subgraphs) {
			s <- paste(s, dputGraph(graph, paste("subgraph",i,sep="")), "\n", sep="")
			i <- i + 1
		}
		s <- paste(s, "weights <- ",dput2(unname(g@weights)),"\n", sep="")
		s <- paste(s, name, " <- new(\"entangledMCP\", subgraphs=list(", paste(paste("subgraph",1:length(g@subgraphs),sep=""), collapse=", "), "), weights=weights)\n", sep="")
		return(s)
	}
	# Normal graphs:
	s <- dputMatrix(g@m, name="m", indent=11, rowNames=TRUE)
	s <- paste(s, "weights <- ",dput2(unname(g@weights)),"\n", sep="")
	s <- paste(s, name, " <- new(\"graphMCP\", m=m, weights=weights)\n", sep="")
	return(s)
}

dputMatrix <- function(m, name, indent=6, rowNames=FALSE) {
  if (isTRUE(indent)) {
    indent <- 6
    if (!missing(name)) {
      indent <- 10+nchar(name)
    }
  }
	s <- "rbind("
	if (!missing(name)) s <- paste(name,"<- rbind(") 
	for (i in 1:(dim(m)[1])) {
		nameLater <- FALSE
		if (any(make.names(row.names(m))!=row.names(m))) {
			rowNames <- FALSE
			nameLater <- TRUE
		}
		rName <- ifelse(rowNames, paste(row.names(m)[i],"=",sep=""), "")
		s <- paste(s, 
			ifelse(i==1,"",paste(rep(" ",indent),collapse="")),
			rName,
			dput2(unname(m[i,])),
			ifelse(i==dim(m)[1],")\n",",\n"),
			sep="")	    
	}
	if (nameLater) {
		if (missing(name)) {
			warning("Can't set row names if no name for matrix is given.")
			return(s)
		}
		s <- paste(s, 
				"row.names(",name,") <- ", dput2(row.names(m)), "\n", sep="")
	}
	return(s)
}

dput2 <- function(x) {
	paste(capture.output(dput(x)), collapse=" ")
}

# This function calculates the number of uncorrelated test statistics given the correlation structure and the number of p-values.
getBalancedDesign <- function (correlation, numberOfPValues) {
	if (correlation == "Dunnett") {
		return (rep(10, numberOfPValues+1))
	}
	stop(paste("The string \"",correlation,"\" does not specify a supported correlation.", sep=""))
}

adjPValues <- function(graph, pvalues, upscale=FALSE, verbose=FALSE) {
	if (length(pvalues)!=length(getNodes(graph))) {
		stop("Length of pvalues must equal number of nodes.")
	}
	if (is.null(names(pvalues))) {
		names(pvalues) <- getNodes(graph)
	}
	# TODO for graphs with sum(weights)<1 (do we want to allow this) this will give wrong results
	if (sum(getWeights(graph))>0) setWeights(graph, getWeights(graph)/sum(getWeights(graph)), getNodes(graph))
	adjPValues <- rep(0, length(getNodes(graph)))
	names(adjPValues) <- getNodes(graph)	
	J <- getNodes(graph)
	names(J) <- getNodes(graph)
	sequence <- list(graph)
	pmax <- 0
	while(length(J) >= 1) {
		fraction <- pvalues[J]/getWeights(graph)[J]
		fraction[pvalues[J]==0] <- 0
		j <- which.min(fraction)
		node <- J[j]
		adjPValues[node] <- max(min(ifelse(pvalues[node]==0,0,pvalues[node]/getWeights(graph)[node]), 1), pmax)
		pmax <- adjPValues[node]
		# if (verbose) cat(paste("We will update the graph with node \"",node,"\".\n",sep=""))
		graph <- rejectNode(graph, node, upscale=upscale, verbose=verbose)
		J <- J[J!=node]
		sequence <- c(sequence, graph)
	}	
	return(new("gMCPResult", graphs=sequence, pvalues=pvalues, adjPValues=adjPValues))
}

#' Rejects a node/hypothesis and updates the graph accordingly.
#' 
#' Rejects a node/hypothesis and updates the graph accordingly.
#' 
#' For details see the given references.
#' 
#' @param graph A graph of class \code{\link{graphMCP}} or \code{\link{entangledMCP}}.
#' @param node A character string specifying the node to reject.
#' @param upscale Logical. If \code{upscale=TRUE} then the weights of all non-rejected
#' nodes are scaled so that the sum is equal to 1. This forces \code{keepWeights=FALSE}
#' to reduce confusion, since otherwise the sum of weights could become bigger than 1.
#' @param verbose Logical scalar.  If \code{TRUE} verbose output is generated
#' during sequentially rejection steps.
#' @param keepWeights Logical scalar. If \code{FALSE} the weight of a node
#' without outgoing edges is set to 0 if it is removed.  Otherwise it keeps its
#' weight.
#' @return An updated graph of class \code{\link{graphMCP}} or \code{\link{entangledMCP}}.
#' @author Kornelius Rohmeyer \email{rohmeyer@@small-projects.de}
#' @seealso \code{\link{graphMCP}}
#' @references Frank Bretz, Willi Maurer, Werner Brannath, Martin Posch: A
#' graphical approach to sequentially rejective multiple test procedures.
#' Statistics in Medicine 2009 vol. 28 issue 4 page 586-604.
#' \url{https://www.meduniwien.ac.at/fwf_adaptive/papers/bretz_2009_22.pdf}
#' @keywords htest graphs
#' @examples
#' 
#' g <- BonferroniHolm(5)
#' 
#' rejectNode(g, "H1")
#' 
#' @export rejectNode
#' 
rejectNode <- function(graph, node, upscale=FALSE, verbose=FALSE, keepWeights=FALSE) {
  if (keepWeights && upscale) {
    warning("Parameter keepWeights is set to FALSE for upscale=TRUE.")
    keepWeights <- FALSE
  }
	# Entangled graphs
	if ("entangledMCP" %in% class(graph)) {
		for(i in 1:length(graph@subgraphs)) {
			graph@subgraphs[[i]] <- rejectNode(graph@subgraphs[[i]], node, upscale=upscale, verbose=verbose, keepWeights=keepWeights)
		}
		return(graph)
	}
	# Normal graphs
	weights <- graph@weights
	graph@weights <- weights+weights[node]*graph@m[node,]
	m <- graph@m	
	for (i in getNodes(graph)) {
		if (m[i, node]*m[node, i]<1) {
			graph@m[i,] <- (m[i,]+m[i,node]*m[node,])/(1-m[i,node]*m[node,i]) 
		} else {
			graph@m[i,] <- 0
		}
	}
	diag(graph@m) <- 0
	graph@m[node,] <- 0
	graph@m[, node] <- 0
	if (!all(m[node,]==0) || !keepWeights) graph@weights[node] <- 0	
	if (upscale && !isTRUE(all.equal(graph@weights, rep(0, length(graph@weights)), check.attributes = FALSE))) {
	  #cat("Scaling weights: ", graph@weights, "->", graph@weights/sum(graph@weights))
	  graph@weights <- graph@weights/sum(graph@weights)
	}
	graph@nodeAttr$rejected[node] <- TRUE	
	return(graph)
}

getRejectableNode <- function(graph, alpha, pvalues) {
	x <- getWeights(graph)*alpha/pvalues
	x[pvalues==0] <- 1
	x[graph@nodeAttr$rejected] <- NaN
	i <- which.max(x)
	if (length(i)==0) return(NULL)
	if (x[i]>1 | all.equal(unname(x[i]),1)[1]==TRUE) {return(getNodes(graph)[i])}
	return(NULL)	 
}