R/plausibleAlts.R

In PoolBal: Balancing Central and Marginal Rejection of Pooled p-Values

 ##' @title Identify a region of plausible alternative hypotheses
##' in the proportion, strength of non-null evidence space
##' @description This function provides a convenient way to interact
##' with simulations performed over a grid of possible alternatives
##' spanning the proportion (eta) and strength (KL divergence) of
##' evidence against the null hypothesis under beta alternatives.
##' @details The simulation this function summarized used a range of
##' eta, w, and KL divergence values to generate thousands of
##' potential alternative distributions. The power of each
##' chi-squared pooled p-value for 161 kappa values ranging from
##' exp(-8) to exp(8) selected uniformly on the log scale was then
##' computed for each alternative using 10,000 simulated examples.
##' Every choice of kappa was compared to the maximum power across
##' all kappas for each setting using a binomial test of differences.
##' This same simulation was repeated twice: once for w values
##' selected uniformly from 0 to 1 and another where selection was
##' uniform on the log scale. The internal data summarizes the results
##' by reporting the count of instances in w (or logw) where a given
##' kappa value was most powerful for a given eta and KL divergence.
##'
##' Though the simulation data is not exported to users and so cannot
##' be accessed directly, this function allows a user to query the
##' data with a range of kappa values (corresponding to those where
##' a given sample seems most powerful) and returns the count of
##' cases in w where a kappa in the corresponding kappa range was most
##' powerful given the eta, KL-divergence combination with beta
##' alternatives. The simulations only spanned kappa values from
##' exp(-8) to exp(8), so providing values outside this range will
##' give very inaccurate results.
##' @param logKappaRange pair of numeric values
##' @param logW logical, should the log scale simulation be used?
##' @return An 81 by 81 matrix giving summarized counts of cases.
##' @examples
##' altFrequencyMat(c(-1, 1), logW = FALSE)
##' altFrequencyMat(c(-1, 1), logW = TRUE)
##' @author Chris Salahub
altFrequencyMat <- function(logKappaRange, logW = FALSE) {
    if (logW) {
        maskMat <- kappaPowersMaskedLogW
    } else {
        maskMat <- kappaPowersMasked
    }
    kapSeq <- as.numeric(dimnames(maskMat)$logk)
    kapInd <- kapSeq <= logKappaRange[2] & kapSeq >= logKappaRange[1]
    if (logKappaRange[2] > max(kapSeq) |
        logKappaRange[1] < min(kapSeq)) {
        warning(paste("kappaRange extends beyond the simulation range",
                      "(exp(-8) to exp(8)), map reported is only",
                      "accurate for values within this range"))
    }
    if (!any(kapInd)) { # if interval is between simulation steps
        kapInd <- logical(length(kapSeq))
        kapInd[which.min(abs(kapSeq - mean(logKappaRange)))] <- TRUE
    }
    ## aggregate
    apply(maskMat[,,kapInd], c(1,2), sum)
}

##' @title Heatmap with marginal histograms
##' @description Display a matrix using a heatmap with marginal
##' histograms.
##' @details This function accepts a matrix of values and plots the
##' matrix with saturation/hue determined by a provided palette
##' argument generated by colorRampPalette, for example. Marginal
##' histograms summarizing the relative frequencies along both
##' dimensions are also plotted to give a complete sense of the
##' individual distributions alongside their joint distribution.
##' This was designed to summarize the alternative distribution
##' space summarized by altFrequencyMat, and the defaults reflect
##' this.
##' @param mat numeric matrix to be plotted
##' @param main title
##' @param xlab x axis label
##' @param ylab y axis label
##' @param pal palette for heatmap
##' @param histFill colour to fill histogram bars
##' @param ... additional arguments to image
##' @return Plot the data using a heatmap and marginal histograms and
##' return nothing.
##' @examples
##' marHistHeatMap(altFrequencyMat(c(0, 2)))
##' @author Chris Salahub
marHistHeatMap <- function(mat, main = "", ylab = expression(eta),
                           xlab = "lnD(a,w)", pal = NULL,
                           histFill = adjustcolor("firebrick", 0.5),
                           ...) {
    if (is.null(pal)) {
        pal <- colorRampPalette(c("white",
                                  "firebrick"))(20)
    }
    image(mat, xaxt = "n", yaxt = "n", col = pal, ylab = "", xlab = "",
          main = main, ...)
    mtext(main, side = 3, line = 0, cex = 0.8) # main
    mtext(ylab, side = 2, line = 1, cex = 0.8) # ylab
    mtext(xlab, side = 1, line = 1, padj = 0, cex = 0.8) # xlab
    ## add ticks
    mtext(side = 1, at = seq(0, 1, by = 0.25), text = "|", line = 0,
          cex = 0.5, padj = -2)
    mtext(text = seq(-5, 5, by = 2.5), at = seq(0, 1, by = 0.25),
          side = 1, cex = 0.8)
    mtext(side = 2, at = seq(0, 1, by = 0.2), text = "|", line = 0,
          cex = 0.5, padj = 1)
    mtext(text = seq(0, 1, by = 0.2), at = seq(0, 1, by = 0.2),
          side = 2, cex = 0.8)
    bds <- par()$usr
    rowDist <- rowSums(mat, na.rm = TRUE)
    colDist <- colSums(mat, na.rm = TRUE) # marginal distributions
    vboxBds <- c(bds[2], bds[2] + 0.1, bds[3:4])
    hboxBds <- c(bds[1:2], bds[4], bds[4] + 0.1)
    ## density boxes
    rect(vboxBds[1], vboxBds[3], vboxBds[2], vboxBds[4], xpd = NA)
    rect(hboxBds[1], hboxBds[3], hboxBds[2], hboxBds[4], xpd = NA)
    ## add marginal histograms
    vseq <- seq(vboxBds[3], vboxBds[4],
                length.out = length(colDist) + 1)
    rect(vboxBds[1], vseq[1:length(colDist)],
         vboxBds[1] + 0.9*diff(vboxBds[1:2])*(colDist/max(colDist)),
         vseq[2:(length(colDist) + 1)], xpd = NA, col = histFill)
    hseq <- seq(hboxBds[1], hboxBds[2],
                length.out = length(rowDist) + 1)
    rect(hseq[1:length(rowDist)], hboxBds[3],
         hseq[2:(length(rowDist) + 1)], xpd = NA,
         hboxBds[3] + 0.9*diff(hboxBds[3:4])*(rowDist/max(rowDist)),
         col = histFill)
}