Multiplicity Adjustment in Plot method

[davidaarmstrong]

I watched a presentation today where someone used your package and it seemed really neat. I asked about whether there was any sort of multiple testing adjustment in the plot and she suggested that wasn't an option even though she wished it was. After looking at your plotting function for the mwa objects, it seems like a trivial option to add. Below is an option.

#' Alternative plot method for MWA
#' 
#' @param x Results from `matchedwake`
#' @param zlim Manually sets the range of the color map of the contour plot, required format is c(MINIMUM,MAXIMUM). Default = NA, i.e. the range is automatically set from the MINIMUM and MAXIMUM values of the estimates.
#' @param plotNAs Boolean indicating whether or not to visualize NA estimates as “no effect” (i.e. 0). Default = TRUE.
#' @param adjust p-value adjustment method, see `?p.adjust` for additional details. 
#' @param ... further arguments passed to or from other methods.
plot_mwa <- function (x, zlim = NA, plotNAs = TRUE, 
                      adjust = c("none", "holm", "hochberg", "hommel", 
                                 "bonferroni", "BH", "BY", 
                                 "fdr"), 
                      ...) 
{
  adj <- match.arg(adjust)
  density <- 3
  lty <- 2
  lwd <- 2
  pdata <- x$estimates
  t_unit <- x$parameters$t_unit
  alpha1 <- x$parameters$alpha1
  alpha2 <- x$parameters$alpha2
  yAxis <- sort(unique(pdata$t_window))
  xAxis <- sort(unique(pdata$spat_window))
  pswcplot <- matrix(nrow = length(xAxis) + 2, ncol = length(yAxis) + 
                       2)
  eswcplot <- matrix(nrow = length(xAxis) + 2, ncol = length(yAxis) + 
                       2)
  pdata$pvalue <- p.adjust(pdata$pvalue, method=adj)
  for (y in 1:length(yAxis)) {
    for (x in 1:length(xAxis)) {
      eswcplot[x + 1, y + 1] <- as.numeric(pdata$estimate[pdata$t_window == 
                                                            yAxis[y] & pdata$spat_window == xAxis[x]])
      pswcplot[x + 1, y + 1] <- as.numeric(pdata$pvalue[pdata$t_window == 
                                                          yAxis[y] & pdata$spat_window == xAxis[x]])
    }
  }
  if (plotNAs) {
    eswcplot[is.na(eswcplot)] <- 0
  }
  for (y in 1:(length(yAxis) + 2)) {
    eswcplot[1, y] <- eswcplot[2, y]
    eswcplot[length(xAxis) + 2, y] <- eswcplot[length(xAxis) + 
                                                 1, y]
  }
  for (x in 1:(length(xAxis) + 2)) {
    eswcplot[x, 1] <- eswcplot[x, 2]
    eswcplot[x, length(yAxis) + 2] <- eswcplot[x, length(yAxis) + 
                                                 1]
  }
  if (length(xAxis) == 1) {
    xvals <- c(0, 0.5, 1)
    xlims <- c(0, 1)
    xat <- 0.5
    xhalfRectSize <- 1
    xstepsize <- 1
  }
  else {
    xvals <- seq(-1/(length(xAxis) - 1), 1 + 1/(length(xAxis) - 
                                                  1), length.out = length(xAxis) + 2)
    xlims <- c(-1/(length(xAxis) - 1)/2, 1 + 1/(length(xAxis) - 
                                                  1)/2)
    xat <- seq(0, 1, by = 1/(length(xAxis) - 1))
    xhalfRectSize <- (1/(length(xAxis) - 1))/2
    xstepsize <- 1/(length(xAxis) - 1)
  }
  if (length(yAxis) == 1) {
    yvals <- c(0, 0.5, 1)
    ylims <- c(0, 1)
    yat <- 0.5
    yhalfRectSize <- 1
    ystepsize <- 1
  }
  else {
    yvals <- seq(-1/(length(yAxis) - 1), 1 + 1/(length(yAxis) - 
                                                  1), length.out = length(yAxis) + 2)
    ylims <- c(-1/(length(yAxis) - 1)/2, 1 + 1/(length(yAxis) - 
                                                  1)/2)
    yat <- seq(0, 1, by = 1/(length(yAxis) - 1))
    yhalfRectSize <- (1/(length(yAxis) - 1))/2
    ystepsize <- 1/(length(yAxis) - 1)
  }
  if (is.na(zlim[1])) {
    zlims <- range(eswcplot[!is.na(eswcplot)], finite = TRUE)
  }
  else {
    zlims <- zlim
  }
  filled.contour(x = xvals, y = yvals, eswcplot, xlab = "Spatial window in kilometers", 
                 ylab = paste("Temporal window in ", t_unit, sep = ""), 
                 nlevels = 20, color.palette = gray.colors, xlim = xlims, 
                 ylim = ylims, zlim = zlims, plot.axes = {
                   axis(1, labels = xAxis, at = xat)
                   axis(2, at = yat, labels = yAxis)
                   for (y in 1:length(yAxis)) {
                     for (x in 1:length(xAxis)) {
                       if (is.na(pswcplot[x + 1, y + 1]) | is.nan(pswcplot[x + 
                                                                           1, y + 1]) | pswcplot[x + 1, y + 1] > alpha2) {
                         rect(c(((x - 1) * xstepsize) - xhalfRectSize, 
                                ((x - 1) * xstepsize) - xhalfRectSize), 
                              c(((y - 1) * ystepsize) - yhalfRectSize, 
                                ((y - 1) * ystepsize) - yhalfRectSize), 
                              c(((x - 1) * xstepsize) + xhalfRectSize, 
                                ((x - 1) * xstepsize) + xhalfRectSize), 
                              c(((y - 1) * ystepsize) + yhalfRectSize, 
                                ((y - 1) * ystepsize) + yhalfRectSize), 
                              density = density, lwd = lwd)
                       }
                       else if (is.na(pswcplot[x + 1, y + 1]) | is.nan(pswcplot[x + 
                                                                                1, y + 1]) | pswcplot[x + 1, y + 1] > alpha1 & 
                                pswcplot[x + 1, y + 1] <= alpha2) {
                         rect(c(((x - 1) * xstepsize) - xhalfRectSize, 
                                ((x - 1) * xstepsize) - xhalfRectSize), 
                              c(((y - 1) * ystepsize) - yhalfRectSize, 
                                ((y - 1) * ystepsize) - yhalfRectSize), 
                              c(((x - 1) * xstepsize) + xhalfRectSize, 
                                ((x - 1) * xstepsize) + xhalfRectSize), 
                              c(((y - 1) * ystepsize) + yhalfRectSize, 
                                ((y - 1) * ystepsize) + yhalfRectSize), 
                              density = density, lty = lty, lwd = lwd)
                       }
                     }
                   }
                 })
}

Above, the first line of the function adj <- match.arg(adjust) and the 15th line pdata$pvalue <- p.adjust(pdata$pvalue, method=adj) are new. Everything else is the same within the function. The function definition also includes the adjust argument now with the relevant options. Here is an example using the function based on the example data from the mwa package. As written above, it allows the methods available in stats::p.adjust().

library(mwa)
data(mwa_data)

# Specify required parameters:
# - 2 to 10 days in steps of 2
t_window <- c(2,10,2)
# - 2 to 10 kilometers in steps of 2
spat_window <- c(2,10,2)
# - column and entries that indicate treatment events 
treatment <- c("type","treatment")
# - column and entries that indicate control events 
control  <- c("type","control")
# - column and entries that indicate dependent events 
dependent <- c("type","dependent")
# - columns to match on
matchColumns <- c("match1","match2")

# Specify optional parameters:
# - use weighted regression (default estimation method is "lm")
weighted <- TRUE
# - temporal units
t_unit <- "days" 
# - match on counts of previous treatment and control events
TCM <- TRUE

results <- matchedwake(mwa_data, t_window, spat_window, treatment, control, dependent,
                       matchColumns, weighted = weighted, t_unit = t_unit, TCM = TCM, 
                       alpha1 = .001, alpha2=.005)
plot_mwa(results, adjust="none")

plot_mwa(results, adjust="bonferroni")

^{Created on 2024-11-18 with reprex v2.1.0}

[kdonnay]

Hi Dave,

Thanks a lot for opening the issue and for providing such a comprehensive example/input. We had it on our list for a long time to implement multiple testing corrections but then our attention shifted elsewhere.

I'd be happy to implement your suggestion in the next mwa release I was planning to do in the coming months and, of course, credit you with the addition. I'm sure this would be something that others also find very useful!

It's always a little bit of hustle to prep packages for CRAN but I think I'll have time before the end of the year to turn this and a few other smaller updates around and make a new version available here and on CRAN.

Would that work for you? Many thanks again for looking into this and for already working out the alternative implementation, Sebastian and I really appreciate it!

Best, Karsten