R/toolsSvenja.R

In dkaschek/dMod: Dynamic Modeling and Parameter Estimation in ODE Models

#' Plot an array of trajectories along the profile of a parameter
#' 
#' @param par Character of parameter name for which the array should be generated.
#' @param profs Lists of profiles as being returned by \link{profile}. 
#' @param prd Named list of matrices or data.frames, usually the output of a prediction function
#' as generated by \link{Xs}.
#' @param times Numeric vector of time points for the model prediction.
#' @param direction Character "up" or "down" indicating the direction the value should be traced along the profile starting at the bestfit value.
#' @param covtable Optional covariate table or condition.grid necessary if subsetting is required.
#' @param ... Further arguments for subsetting the plot.
#' @param nsimus Number of trajectories/ simulation to be calculated.
#' 
#' @return A plot object of class \code{ggplot}.
#' @author Svenja Kemmer, \email{svenja.kemmer@@fdm.uni-freiburg.de}
#' @examples
#' \dontrun{
#'  plotArray("myparameter", myprofiles, g*x*p, seq(0, 250, 1), 
#'     "up", condition.grid, name == "ProteinA" & condition == "c1") 
#' }
#' @export
#' @import data.table
plotArray <- function (par, profs, prd, times, direction = c("up", "down"), covtable, ..., nsimus = 4) {
  
  # select subframe from profiles
  mysub <- profs %>% as.data.table() %>% .[whichPar == par, ]
  mysub[, ID := 1:nrow(mysub)]
  
  # get ID of bestfit (constraint is 0 for bestfit)
  bestID <- mysub[constraint == 0.00]$ID
  if(direction == "up") mysubF <- mysub[ID >= bestID]  
  if(direction == "down") mysubF <- mysub[ID <= bestID]
  
  # select rows according to simulation number
  partable <- mysubF[seq(1, nrow(mysubF), (round(nrow(mysubF)/nsimus)))]
  
  # remove non_parameter names
  no_pars <- c("value", "constraint", "stepsize", "gamma", "whichPar", "data", "condition_obj", "AIC", "BIC", "prior", "ID", "chisquare")
  partable %>% .[, (no_pars) := NULL]
  
  # make predictions
  predictionDT <- predict_array(prd, times, pars = partable, whichpar = par)
  out_plot <- copy(predictionDT)
  
  # use covtable for subsetting of the plot
  if(!is.null(covtable)) {
    if(!"condition" %in% names(covtable)){
      covtable <- as.data.table(covtable, keep.rownames = "condition")
    } else covtable <- as.data.table(covtable)
    out_plot <- merge(out_plot, covtable, by = "condition")
    out_plot <- out_plot[...]
  }
  
  # plot
  P <- ggplot(out_plot , aes(x = time, y = value, group = ParValue, color = ParValue)) +
    facet_grid(name~condition, scales = "free_y") +
    geom_line(size = 1) + 
    theme_dMod(base_size = 18) + scale_color_viridis_c() +
    theme(legend.position = "top", legend.key.size = unit(0.6,"cm")) + 
    theme(axis.line = element_line(colour = "black"), 
          panel.grid.major = element_line(colour = "grey97"), 
          panel.grid.minor = element_line(colour = "grey97"), 
          panel.background = element_blank()) +
    xlab("time") +
    ylab(paste0("value"))
  
  return(P)
}

predict_array <- function (prd, times, pars = partable, whichpar = par, keep_names = NULL, FLAGverbose = FALSE, FLAGverbose2 = FALSE, FLAGbrowser = FALSE, ...) {
  if (FLAGverbose2) cat("Simulating", "\n")
  out <- lapply(1:nrow(pars), function(i) {
    if (FLAGverbose) cat("Parameter set", i, "\n")
    if (FLAGbrowser) browser()
    mypar <- pars[i,] %>% as.numeric()
    parval <- round(pars[i,][[whichpar]], digits = 2)
    names(mypar) <- names(pars)
    mypar <- as.parvec(mypar)
    prediction <- try(prd(times, mypar, deriv = FALSE, ...))
    if (inherits(prediction, "try-error")) {
      warning("parameter set ", i, " failed\n")
      return(NULL)
    }
    prediction <- purrr::imap(prediction, function(.x,.y){
      .x <- data.table(.x)
      if (!is.null(keep_names))
        .x[, (setdiff(names(.x), c(keep_names, "time"))) := NULL]
      .x[, `:=`(condition = .y, ParValue = parval)]
      .x
    })
    melt(rbindlist(prediction), variable.name = "name", value.name = "value", id.vars = c("time", "condition", "ParValue"))
  })
  if (FLAGverbose2) cat("postprocessing", "\n")
  out <- rbindlist(out[!is.null(out)])
  out
}

PlotPaths <- function(profs=myprofiles, ..., whichPar, sort = FALSE, relative = TRUE, scales = "fixed", multi = TRUE, n_pars = 5) {
  
  if ("parframe" %in% class(profs)) {
    arglist <- list(profs)
  } else {
    arglist <- as.list(profs)
  }
  
  if (is.null(names(arglist))) {
    profnames <- 1:length(arglist)
  } else {
    profnames <- names(arglist)
  }
  
  
  data <- do.call(rbind, lapply(1:length(arglist), function(i) {
    # choose a proflist
    proflist <- as.data.frame(arglist[[i]])
    parameters <- attr(arglist[[i]], "parameters")
    
    if (is.data.frame(proflist)) {
      whichPars <- unique(proflist$whichPar)
      proflist <- lapply(whichPars, function(n) {
        with(proflist, proflist[whichPar == n, ])
      })
      names(proflist) <- whichPars
    }
    
    if (is.null(whichPar)) whichPar <- names(proflist)
    if (is.numeric(whichPar)) whichPar <- names(proflist)[whichPar]
    
    subdata <- do.call(rbind, lapply(whichPar, function(n) {
      # matirx
      paths <- as.matrix(proflist[[n]][, parameters])
      values <- proflist[[n]][, "value"]
      origin <- which.min(abs(proflist[[n]][, "constraint"]))
      if (relative) 
        for(j in 1:ncol(paths)) paths[, j] <- paths[, j] - paths[origin, j]
      
      combinations <- expand.grid.alt(whichPar, colnames(paths))
      if (sort) combinations <- apply(combinations, 1, sort) else combinations <- apply(combinations, 1, identity)
      combinations <- submatrix(combinations, cols = -which(combinations[1,] == combinations[2,]))
      combinations <- submatrix(combinations, cols = !duplicated(paste(combinations[1,], combinations[2,])))
      
      
      
      
      path.data <- do.call(rbind, lapply(1:dim(combinations)[2], function(j) {
        data.frame(chisquare = values, 
                   name = n,
                   proflist = profnames[i],
                   combination = paste(combinations[,j], collapse = " - \n"),
                   x = paths[, combinations[1,j]],
                   y = paths[, combinations[2,j]])
      }))
      
      if(multi) path.data <- path.data %>% as.data.table %>% .[, partner := tstrsplit(as.character(combination), "\n", fixed=TRUE, keep = 2)]
      
      
      return(path.data)
      
    }))
    
    return(subdata)
    
  }))
  
  data$proflist <- as.factor(data$proflist)
  
  if (relative){
    axis.labels <- c(expression(paste(Delta, "parameter 1")), expression(paste(Delta, "parameter 2")))  
  } else {
    axis.labels <- c("parameter 1", "parameter 2")
  }
  
  data <- droplevels(subset(data, ...))
  
  suppressMessages(
    p <- ggplot(data, aes(x = x, y = y, group = interaction(name, proflist), color = name, lty = proflist)) + 
      facet_wrap(~combination, scales = scales) + 
      geom_path() + #geom_point(aes=aes(size=1), alpha=1/3) +
      xlab(axis.labels[1]) + ylab(axis.labels[2]) +
      scale_linetype_discrete(name = "profile\nlist") +
      scale_color_manual(name = "profiled\nparameter", values = dMod_colors)
  )
  if(multi){
    
    # determine strength of change
    data[, max.dev := max(c(abs(max(y)), abs(min(y)))), by = "partner"]
    setorder(data, name, -max.dev)
    # max.devis <- unique(data$max.dev)[1:n_pars]
    data[!(max.dev %in% unique(max.dev)[1:n_pars]), partner := "others"]
    
    data$combination <- as.factor(data$combination)
    data$partner <- factor(data$partner, levels = unique(data$partner))
    
    suppressMessages(
      p <- ggplot(data, aes(x = x, y = y, color = partner)) + 
        geom_path() + #geom_point(aes=aes(size=1), alpha=1/3) +
        xlab(paste0("log(", whichPar, ")")) + ylab("relative change of\n other paramters") +
        scale_linetype_discrete(name = "profile\nlist") +
        scale_color_manual(values = c(dMod_colors[2:(n_pars+1)], rep("gray", 100))) + theme_dMod() +
        theme(legend.position="bottom",
              legend.title = element_blank(),
              legend.box.background = element_rect(colour = "black"),
              legend.key.size = unit(0.4, "cm"))
    )
  }
  
  attr(p, "data") <- data
  return(p)
  
}

#' Profile likelihood: plot all parameter paths belonging to one profile in one plot
#' 
#' @param profs Lists of profiles as being returned by \link{profile}. 
#' @param whichpars Character vector of parameter names for which the profile paths should be generated.
#' @param npars Numeric vector of number of colored and named parameter paths.
#' 
#' @return A plot object of class \code{ggplot} for length(whichpars) = 1 and otherwise an object of class \code{cowplot}.
#' @author Svenja Kemmer, \email{svenja.kemmer@@fdm.uni-freiburg.de}
#' @examples
#' \dontrun{
#'  plotPathsMulti(myprofiles, c("mypar1", "mypar2"), npars = 5) 
#' }
#' @export
#' @import data.table
plotPathsMulti <- function(profs, whichpars, npars = 5) {
  if(length(whichpars) == 1){
    p <- PlotPaths(profs=profs, whichPar = whichpars, n_pars = npars)
    return(p)
  } else {
    PlotList <- NULL
    for(i in 1:length(whichpars)){
      par <- whichpars[i]
      p <- PlotPaths(profs=profs, whichPar = par, n_pars = npars)
      PlotList[[i]] <- p
    }
    pl <- cowplot::plot_grid(plotlist = PlotList)
    return(pl)
  }
}



expand.grid.alt <- function(seq1, seq2) {
  cbind(Var1=rep.int(seq1, length(seq2)), Var2=rep(seq2, each=length(seq1)))
}

#' Profile likelihood: plot profiles along with their parameter paths
#' 
#' @param profs Lists of profiles as being returned by \link{profile}. 
#' @param whichpars Character vector of parameter names for which the profile paths should be generated.
#' @param npars Numeric vector of colored and named parameter paths.
#' 
#' @return A plot object of class \code{ggplot}.
#' @author Svenja Kemmer, \email{svenja.kemmer@@fdm.uni-freiburg.de}
#' @examples
#' \dontrun{
#'  plotProfilesAndPaths(myprofiles, c("mypar1", "mypar2"), npars = 5) 
#' }
#' @export
#' @import data.table
plotProfilesAndPaths <- function(profs, whichpars, npars = 5){
  if(length(whichpars)<2) {
    profs <- profs[profs$whichPar %in% whichpars]
    pl1 <- plotProfile(profs, mode == "data")
    pl2 <- plotPathsMulti(profs, whichpars, npars)
    pl <- cowplot::plot_grid(pl1,cowplot::plot_grid(NULL,pl2, nrow = 1, rel_widths = c(0.2,1)),nrow = 2, rel_heights = c(1,0.7))
  } else{
    plotList <- NULL
    for(z in 1:length(whichpars)){
      prof_sub <- profs[profs$whichPar == whichpars[z]]
      pl1 <- plotProfile(prof_sub, mode == "data") + theme(legend.position = "none")
      pl2 <- plotPathsMulti(prof_sub, whichpars[z], npars)
      pl <- cowplot::plot_grid(pl1,cowplot::plot_grid(NULL,pl2, nrow = 1, rel_widths = c(0.2,1)),nrow = 2, rel_heights = c(1,0.7))
      plotList[[z]] <- pl
    }
    plot <- cowplot::plot_grid(plotlist = plotList, nrow = 1)
    print(plot)
  }
}