R/generateTest.R

In gMCP: Graph Based Multiple Comparison Procedures

#' generateTest
#' 
#' generates a test function for the multiple comparison procedure with
#' correlated test statistics defined by a graph
#' 
#' 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{cr} 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{cr[i,i] = 1} for diagonal elements,
#' \eqn{cr[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 cr[1,2]=0
#' indicates that the first and second test statistic are uncorrelated, whereas
#' cr[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,k) for
#' i!=j!=k are specified then cor(j,k) has to be specified as well) otherwise
#' the corresponding intersection null hypotheses tests are not uniquely
#' defined and an error is returned.
#' 
#' The parametric tests in (Bretz et al. (2011)) are defined such that the
#' tests of intersection null hypotheses always upscale the full alpha level
#' even if the sum of weights is strictly smaller than one. This has the
#' consequence that certain test procedures that do not test each intersection
#' null hypothesis at the full level alpha may not be implemented (e.g., a
#' single step Dunnett test). If \code{upscale} is set to \code{FALSE}
#' (default) the parametric tests are performed at a reduced level alpha of
#' sum(w) * alpha. If set to
#' \code{TRUE} the tests are performed as defined in Equation (3) of (Bretz et
#' al. (2011)).
#' 
#' @param g graph defined as a matrix, each element defines how much of the
#' local alpha reserved for the hypothesis corresponding to its row index is
#' passed on to the hypothesis corresponding to its column index
#' @param w vector of weights, defines how much of the overall alpha is
#' initially reserved for each elementary hypothesis
#' @param cr correlation matrix if p-values arise from one-sided tests with
#' multivariate normal distributed test statistics for which the correlation is
#' partially known. Unknown values can be set to NA. (See details for more
#' information)
#' @param al overall alpha level at which the family error is controlled
#' @return Returns a function that will take a vector of z-scores to which the
#' test will be applied. This function in turn will return a boolean vector
#' with elements false if the particular elementary hypothesis can not be
#' rejected and true otherwise.
#' @param upscale if \code{FALSE} (default) the parametric tests are performed 
#' at a reduced level alpha of sum(w) * alpha. (See details)
#' @author Florian Klinglmueller
#' @references Bretz F, Maurer W, Brannath W, Posch M; (2008) - A graphical
#' approach to sequentially rejective multiple testing procedures. - Stat Med -
#' 28/4, 586-604 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 - to appear
#' @keywords htest
#' @examples
#' 
#'  ## Define some graph as matrix
#'  g <- matrix(c(0,0,1,0,
#'                0,0,0,1,
#'                0,1,0,0,
#'                1,0,0,0), nrow = 4,byrow=TRUE)
#'  ## Choose weights
#'  w <- c(.5,.5,0,0)
#'  ## Some correlation (upper and lower first diagonal 1/2)
#'  c <- diag(4)
#'  c[1:2,3:4] <- NA
#'  c[3:4,1:2] <- NA
#'  c[1,2] <- 1/2
#'  c[2,1] <- 1/2
#'  c[3,4] <- 1/2
#'  c[4,3] <- 1/2
#' 
#'  ## Test function for further use:
#'  myTest <- generateTest(g,w,c,.05)
#'  myTest(c(3,2,1,2))
#' 
#' @export generateTest
#' 
generateTest <-
function(g,w,cr,al, upscale=FALSE){#, alternatives="less"){
  bounds <- generateBounds(g,w,cr,al,upscale=upscale)#, alternatives=alternatives)
  return(
         function(z){
           dm <- t(apply(bounds,1,function(b) {
             ## check whether z is larger than boundary
             d <- rep(NA,length(b))
             d[!is.na(b)] <- (b[!is.na(b)]<=z[!is.na(b)])
             return(d)
           }))
           d <- apply(dm,2,function(h) {
             ## closed testing
             #min(rowSums(dm[!is.na(h),],na.rm=TRUE)>0)>0
             all(apply(dm[!is.na(h),],1,any,na.rm=T))
           })
           d
         })
}