R/momApp_modalityTest.R

In CharlotteJana/pdmppoly: Moment Approximation for Polynomial PDMPs

#' Perform modality tests on approximated moments
#' 
#' Test if the moments generated by \code{\link{momApp}} 
#' come from a unimodal distribution with compact support.
#' The test used is given by function \code{\link[momcalc]{is.unimodal}}
#' from package \pkg{momcalc}. The lower and upper bounds of
#' the compact support may depend on a time value \code{t}.
#' 
#' @param momApp an object of class \code{\link{momApp}}.
#' @param vars a character vector specifiying the names of the variables
#' on which a modality test should be performed. The variables should be
#' part of the PDMP \code{momApp$model}.
#' @param lower numeric vector or matrix or data.frame specifying the lower
#'   bounds of the compact distribution that determines the law of the PDMP
#'   \code{momApp$model}. If \code{lower} is a vector, the i-th entry should
#'   give the lower bound of the i-th variable given in \code{vars}, independent
#'   of a time value. If \code{lower} is a matrix or data.frame, it should have
#'   a column named \code{time} containing all time values specified in
#'   \code{times(momApp$model)}. The other column names should be identical to
#'   parameter \code{vars} and contain the lower bounds for the corresponding
#'   variable and time value.
#' @param upper numeric vector or matrix or data.frame specifying the upper
#'   bounds of the compact distribution that determines the law of the PDMP
#'   \code{momApp$model}. If \code{upper} is a vector, the i-th entry should
#'   give the upper bound of the i-th variable given in \code{vars}, independent
#'   of a time value. If \code{upper} is a matrix or data.frame, it should have
#'   a column named \code{times} containing all time values specified in
#'   \code{times(momApp$model)}. The other column names should be identical to
#'   parameter \code{vars} and contain the upper bounds for the corresponding
#'   variable and time value.
#' @return A data.frame containing the test results.
#' @examples
#' data(genePolyBF)
#' support <- getSupport(genePolyBF)
#' modalityTest(momApp = momApp(genePolyBF), vars = "f",
#'              lower = support$lower, upper = support$upper)
#' @seealso \code{\link{plotModalities}} to plot the result.
#' @importFrom tidyr spread
#' @importFrom momcalc is.unimodal
#' @aliases modalitytest
#' @export
modalityTest <- function(momApp, lower, upper, 
                         vars = names(init(momApp$model))){
  
  if(is.vector(lower))
    stopifnot(length(lower) == length(vars))
  if(is.vector(upper))
    stopifnot(length(upper) == length(vars))
  if(is.matrix(lower) | is.data.frame(lower))
    stopifnot(colnames(lower) %in% c("time", names(init(momApp$model))))
  if(is.matrix(upper) | is.data.frame(upper))
    stopifnot(colnames(upper) %in% c("time", names(init(momApp$model))))
  
  modalities <- data.frame()
  
  for(calcMethod in unique(momApp$moments$method)){

    for(i in seq_along(vars)){
      
      # select moments
      m <- subset(momApp$moments, method == calcMethod & order <= 4, 
                  select = c("time", "order", vars[i]))
      
      # add lower to m
      if(is.matrix(lower) | is.data.frame(lower)){
        m2 <- lower[, c("time", vars[i])]
        m2$order <- "lower"
      }
      else{
        m2 <- expand.grid("time" = fromtoby(times(momApp$model)),
                         "order" = "lower",
                         "vars" = lower[i])
        colnames(m2)[3] <- vars[i]
      }
      m <- rbind(m, m2)
      
      # add upper to m
      if(is.matrix(upper) | is.data.frame(upper)){
        m2 <- upper[, c("time", vars[i])]
        m2$order <- "upper"
      }
      else{
        m2 <- expand.grid("time" = fromtoby(times(momApp$model)),
                          "order" = "upper",
                          "vars" = upper[i])
        colnames(m2)[3] <- vars[i]
      }
      m <- rbind(m, m2)
      
      # spread m
      m <- tidyr::spread(m, order, vars[i])
      m <- m[order(m$time),]

      # perform test
      m3 <- momcalc::is.unimodal(
        lower = m[, "lower"], upper = m[, "upper"],
        moments = m[ , !(names(m) %in% c("time", "lower", "upper"))]
      )
      
      modalityMethod <- data.frame(time = m$time,
                                   method = calcMethod,
                                   variable = vars[i],
                                   modality = factor(m3, levels = c( "4-b-unimodal",
                                                                     "not unimodal",
                                                                     "not existant",
                                                                     NA_character_)))
      modalities <- dplyr::bind_rows(modalities, modalityMethod)
    }
  }
  modalities$method <- as.factor(modalities$method)
  return(modalities)
}

#' Plot the result of \code{\link{modalityTest}}
#' 
#' Test if the moments generated by \code{\link{momApp}} 
#' come from a unimodal distribution with compact support
#' and plot the result. The test is done by function 
#' \code{\link{modalityTest}}. Already existing test results
#' can be passed via the argument \code{modalities}.
#' 
#' @param momApp an object of class \code{\link{momApp}}.
#' @param modalities The result of \code{modalityTest(momApp, ...)}. If
#'   \code{modalities = NULL}, function \code{\link{modalityTest}} will be
#'   called inside \code{plotModalities}.
#' @param ... additional arguments to \code{\link{modalityTest}}. They are
#' only used if \code{modalities = NULL}.
#' @param vars character vector with variable names. Only the test results for
#' these variables will appear in the plot. Vector \code{vars} should be equal
#' or a subset to \code{names(init(momApp$model))}.
#' @return an object of class \code{ggplot}
#' @aliases plotmodalities plotModality plotmodality
#' @importFrom ggplot2 ggplot geom_segment guides labs facet_wrap
#' @export
plotModalities <- function(momApp, modalities = NULL, 
                           vars = names(init(momApp$model)), ...){
  
  if(is.null(modalities))
    modalities <- modalityTest(momApp, vars, ...)
  
  modalities <- subset(modalities, variable %in% vars)
  
  timedist <- times(momApp$model)["by"]
  
  plot <- ggplot(data = modalities, aes(x = time, y = method, color = modality)) + 
          geom_segment(aes(xend = time + timedist, yend = method), size = 5, lineend = "butt") +
          guides(colour = ggplot2::guide_legend(override.aes = list(size=5))) +
          labs(
            title = descr(momApp$model),
            subtitle = format(momApp$model, slots = c("parms"), short = FALSE)) +
          facet_wrap(~ variable, ncol = 1)
  
  if(!is.null(momApp$seeds))
    plot <- plot + 
            labs(caption = paste("number of simulations:", length(momApp$seeds)))
  
  return(plot)
}