R/tree-simulation-functions.R

In pneuvial/sanssouci: Post Hoc Multiple Testing Inference

#' Generate an unstructured signal
#' 
#' @param m An integer value, the number of hypotheses
#' @param pi0 A numeric value in \eqn{[0,1]}, the proportion of true null signals
#' @param barmu A numeric value, the strength of the signal 
#' @return A numeric vector of length \code{m}, the signal values for each hypothesis
#' @export
#' @examples
#' m <- 100
#' pi0 <- 0.7
#' barmu <- 4
#' mu <- gen.mu.noleaves(m = m, pi0 = pi0, barmu = barmu) 
#' plot(mu)
gen.mu.noleaves <- function(m, pi0, barmu) {
    mu <- numeric(m)
    m0 <- floor(m * pi0)
    m1 <- m - m0
    mu[sample(seq(1, m), m1)] <- barmu
    mu
}

#' Generate a block of deterministic signal of Gaussian shape
#' 
#' @param barmu A numeric value, the strength of the signal 
#' @param len An integer value, the length of the block
#' @return A numeric vector of length \code{len}, the signal values in the block
gauss_bloc <- function(barmu, len) {
    # this choice of sigma2 is so that the values at the edge of the vector 'out' are = 1/200
    sigma2 <- len^2/(8 * log(200))
    out <- exp(-(seq_len(len) - floor(len/2))^2/(2 * sigma2))
    return(barmu * out/mean(out))
}

#' Generate the signal in the leaves of the tree
#' 
#' @param m An integer value, the number of hypotheses
#' @param K1 An integer value, the number of non-null leaves
#' @param d A numeric value in \eqn{[0,1]}, the proportion of non-null signals in each non-null leaf
#' @param grouped A boolean value, whether the non-null leaves are contiguous or not
#' @param setting A character value specifying the shape of the signal in each non-null leaf
#' @param barmu A numeric value, the strength of the signal 
#' @param leaf_list A list of leaves as generated by \code{dyadic.from.Nnn} functions
#' @return A numeric vector of length \code{m}, the signal values for each hypothesis
#' @details If \code{setting == "const"} the non-null signal is constant equal to \code{barmu}. If \code{setting == "gauss"} and \eqn{d=1} the signal in an active leaf has the shape of a Gaussian bell with mean \code{barmu}, see the \code{\link{gauss_bloc}} function which generates it. If \eqn{d<1} it has the same shape but randomly pruned. If \code{setting == "poisson"} the non-null signal is randomly drawn according to a Poisson distribution of mean \code{barmu} (slightly modified to not yield 0). If \code{setting == "rgauss"} the non-null signal is randomly drawn according to a Normal distribution of mean \code{barmu} and variance \code{barmu}.
#' @seealso \code{\link{gauss_bloc}}
#' @export
#' @importFrom stats rpois rnorm
#' @examples
#' m <- 160
#' s <- 10
#' K1 <- floor(m/(s * 4))
#' d <- 1
#' barmu <- 4
#' dd <- dyadic.from.window.size(m, s, method = 2)
#' leaf_list <- dd$leaf_list
#' muC <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = FALSE, 
#'                     setting = "const", barmu = barmu, leaf_list =leaf_list)
#' muC <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = TRUE, 
#'                     setting = "const", barmu = barmu, leaf_list =leaf_list)
#' muG <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = FALSE, 
#'                     setting = "gauss", barmu = barmu, leaf_list =leaf_list)
#' muP <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = FALSE, 
#'                     setting = "poisson", barmu = barmu, leaf_list =leaf_list)
#' mu <- cbind(constant = muC, Gaussian = muG, Poisson = muP)
#' matplot(mu, t = 's')
gen.mu.leaves <- function(m, K1, d, grouped, setting, barmu, leaf_list) {
    mu <- numeric(m)
    K <- length(leaf_list)
    if (K1 > K) 
        stop("K1>K,\nwe don't have so many leaves")
    active_leaves <- numeric(K1)
    if (grouped) {
        active_leaves <- seq(1, K1) # + sample(seq(0, K - K1), 1)
    } else {
        active_leaves <- sample(seq(1, K), K1)
    }
    for (i in active_leaves) {
        length_leaf <- length(leaf_list[[i]])
        m1loc <- floor(length_leaf * d)
        signal <- sample(length_leaf, m1loc)
        mu[leaf_list[[i]][signal]] <- switch(setting, const = {
            barmu
        }, rgauss = {
            rnorm(m1loc, mean = barmu, sd = sqrt(barmu))
        }, gauss = {
        	  gauss_bloc(barmu, m1loc)
            # gauss_bloc(barmu, length_leaf)[signal] # old code, broken because
        	  # that way the signal doesn't have mean barmu 
        	  # but lesser mean as soon as d<1
        }, poisson = {
            rpois(m1loc, 999 * barmu/1000) + barmu/1000
        })
    }
    return(mu)
}

#' Generate one-sided p-values associated to a given signal with equi-correlated noise
#' 
#' @param m An integer value, the number of hypotheses
#' @param mu A vector of \eqn{m} signal values
#' @param rho A numeric value in \eqn{[0,1]}, the level of equi-correlation between test statistics
#' @param alternative A character string specifying the alternative hypothesis.
#'   Must be one of "two.sided" (default), "greater" or "less".
#' @return A vector of \eqn{m} one-sided \eqn{p}-values
#' @export
#' @importFrom stats pnorm
#' @examples
#' m <- 100
#' s <- 10
#' K1 <- floor(m/(s * 4))
#' d <- 1
#' barmu <- 4
#' dd <- dyadic.from.window.size(m, s, method = 2)
#' leaf_list <- dd$leaf_list
#' mu <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = FALSE, 
#'                     setting = "const", barmu = barmu, leaf_list =leaf_list)
#' pvals <- gen.p.values(m = m, mu = mu, rho = 0)
#' plot(-log(pvals), t = 'b')
gen.p.values <- function(m, mu = rep(0, m), rho = 0, alternative = c("two.sided", "less", "greater")) {
    alternative <- match.arg(alternative)
    stopifnot(m == length(mu)) ## sanity check
    Z <- rnorm(m + 1, 0, 1)
    Y <- sqrt(1 - rho) * Z[seq_len(m)] + sqrt(rho) * Z[m + 1]
    if (alternative == "two.sided"){
      return(2 * pnorm(abs(Y + mu), lower.tail = FALSE))
    }
    else if (alternative == "less"){
    	return(pnorm(Y + mu, lower.tail = TRUE))
    }
    else if (alternative == "greater"){
    	return(pnorm(Y + mu, lower.tail = FALSE))
    }
}