R/plotTimeSeries.R

In BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics

#' Plots a time series, with the option to include confidence and prediction band
#' @author Florian Hartig
#' @param x optional values for x axis (time)
#' @param observed observed values
#' @param predicted predicted values
#' @param confidenceBand matrix with confidenceBand
#' @param predictionBand matrix with predictionBand
#' @param xlab a title for the x axis
#' @param ylab a title for the y axis
#' @param ... further arguments passed to \code{\link[graphics]{plot}}
#' @details Values for confidence and prediction bands can be generated with \code{\link{getPredictiveIntervals}}. For a more elaborate version of this plot, see \code{\link{plotTimeSeriesResults}}
#' @seealso \code{\link{marginalPlot}}, \code{\link{tracePlot}}, \code{\link{correlationPlot}}
#' @example /inst/examples/plotTimeSeriesHelp.R
#' @export
plotTimeSeries <- function(observed = NULL, predicted = NULL, x = NULL, confidenceBand = NULL, predictionBand = NULL, xlab = "Time", ylab = "Observed / predicted values", ...){
  
  ylim = range(observed, predicted, confidenceBand, predictionBand,na.rm=TRUE)
  
  if (is.null(x)){
    if(!is.null(observed)) x = 1:length(observed)
    else if(!is.null(predicted)) x = 1:length(predicted)
    else stop("either observed or predicted must be supplied")
  }
  
  len = length(x)
  
  plot(x, y = rep(0, len), ylim = ylim, type = "n", xlab = xlab, ylab = ylab, ...)
  
  if(!is.null(predictionBand)) polygon(c(x,rev(x)),c(predictionBand[1,],predictionBand[2,len:1]),col="moccasin",border=NA)
  if(!is.null(confidenceBand)) polygon(c(x,rev(x)),c(confidenceBand[1,],confidenceBand[2,len:1]),col="#99333380",border=NA)    
    
  if(!is.null(predicted)) lines(x, predicted, col = "red")
  if(!is.null(observed)) points(x, observed, col = "black", pch = 3, cex = 0.6)
}



#' Plots residuals of a time series
#' @author Florian Hartig
#' @param residuals x
#' @param x optional values for x axis (time)
#' @param main title of the plot
#' @export
plotTimeSeriesResiduals <- function(residuals, x = NULL, main = "residuals"){
  
  ylim = range(residuals)
  
  if (is.null(x)){
    x  = 1:length(residuals)
  }
  barplot(residuals)
}


#' Creates a time series plot typical for an MCMC / SMC fit
#' @author Florian Hartig
#' @param sampler Either a) a matrix b) an MCMC object (list or not), or c) an SMC object
#' @param model function that calculates model predictions for a given parameter vector
#' @param observed observed values as vector
#' @param error function with signature f(mean, par) that generates observations with error (error = stochasticity according to what is assumed in the likelihood) from mean model predictions. Par is a vector from the matrix with the parameter samples (full length). f needs to know which of these parameters are parameters of the error function. See example in \code{\link{VSEM}}
#' @param start numeric start value for the plot (see \code{\link{getSample}})
#' @param plotResiduals logical determining whether residuals should be plotted
#' @param prior if a prior sampler is implemented, setting this parameter to TRUE will draw model parameters from the prior instead of the posterior distribution
#' @param ... further arguments passed to \code{\link[graphics]{plot}}
#' @seealso \code{\link{getPredictiveIntervals}}
#' @example /inst/examples/plotTimeSeriesHelp.R
#' @export
plotTimeSeriesResults <- function(sampler, model, observed, error = NULL, plotResiduals = TRUE, start = 1, prior = FALSE, ...){
  oldPar = par(no.readonly = TRUE)
  
  if (plotResiduals == TRUE && is.null(error)) {
    warning("Can not plot residuals without an error function.")
  }
  
  if(!is.vector(observed)) stop("wrong type given to observed. Must be a vector")
  
  if (plotResiduals == TRUE && !is.null(error)) {
    layout(matrix(c(1, 1, 1, 2, 3, 4), 2, 3, byrow = TRUE))
    par(mar = c(3, 3, 3, 3), oma = c(2, 2, 2, 2))
  }
  
  # ... can we pass on to both getSample and plot? 
  
  if(prior == FALSE){
    if(inherits(sampler,"bayesianOutput")) parMatrix = getSample(sampler, start = start)
    else if (inherits(sampler, "matrix")) parMatrix = sampler
    else if ("mcmc.list" %in% class(sampler) || "mcmc" %in% class(sampler)) parMatrix = getSample(sampler, start = start)
    else stop("wrong type given to variable sampler")    
  }else if (prior == TRUE){
    if(inherits(sampler,"bayesianOutput")) {
      if(inherits(sampler, "mcmcSamplerList")) parMatrix = sampler[[1]]$setup$prior$sampler(1000)
      else parMatrix <- sampler$setup$prior$sampler(1000)
    } else {
      stop("prior==TRUE is only available for sampler of type bayesianOutput") 
    }
  }else stop("BayesianTools::plotTimeSeriesResults - wrong argument to prior")

  numSamples = min(1000, nrow(parMatrix))
  
  pred <- getPredictiveIntervals(parMatrix = parMatrix,
                                 model = model,
                                 numSamples = numSamples,
                                 quantiles = c(0.025, 0.5, 0.975),
                                 error = error)
  
  if(!is.null(error)) plotTimeSeries(observed = observed,
                                     predicted = pred$posteriorPredictivePredictionInterval[2,],
                                     confidenceBand = pred$posteriorPredictiveCredibleInterval[c(1,3),],
                                     predictionBand = pred$posteriorPredictivePredictionInterval[c(1,3),],
                                     ...)
  else plotTimeSeries(observed = observed, predicted = pred$posteriorPredictiveSimulations,
                      confidenceBand = pred$posteriorPredictiveSimulations[c(1,3),],
                      ...)
  
  if (plotResiduals && !is.null(error)) {
    dh = getDharmaResiduals(model = model,
                            parMatrix = parMatrix,
                            numSamples = numSamples,
                            observed = observed,
                            error = error,
                            plot = FALSE)
    
    
    # qq-plot
    gap::qqunif(dh$scaledResiduals, pch=2, bty="n", logscale = F, col = "black", cex = 0.6, main = "QQ plot residuals", cex.main = 1)
    
    # residuals vs fitted
    noNa <- which(!is.na(dh$scaledResiduals))
    
    DHARMa::plotResiduals(dh$fittedPredictedResponse[noNa], dh$scaledResiduals[noNa], main = "Residual vs. predicted\n quantile lines should be\n horizontal lines at 0.25, 0.5, 0.75", cex.main = 1, xlab = "Predicted value", ylab = "Standardized residual")
    
    # residuals vs time
    t <- 1:length(dh$fittedPredictedResponse[noNa])
    DHARMa::plotResiduals(t, dh$scaledResiduals[noNa], xlab = "Time", ylab = "Standardized residual", main = "Residual vs. time\n quantile lines should be\n horizontal lines at 0.25, 0.5, 0.75", cex.main = 1)
    
    message("DHARMa::plotTimeSeriesResults called with posterior predictive (residual) diagnostics. Type vignette(\"DHARMa\", package=\"DHARMa\") for a guide on how to interpret these plots")

  }
  
  
  par(oldPar)
}

#' Creates a DHARMa object
#' @author Tankred Ott
#' @param model function that calculates model predictions for a given parameter vector
#' @param parMatrix a parameter matrix from which the simulations will be generated
#' @param numSamples the number of samples
#' @param observed a vector of observed values
#' @param error function with signature f(mean, par) that generates error expectations from mean model predictions. Par is a vector from the matrix with the parameter samples (full length). f needs to know which of these parameters are parameters of the error function
#' @param plot logical, determining whether the simulated residuals should be plotted
# #' @export
getDharmaResiduals <- function(model, parMatrix, numSamples, observed, error, plot = TRUE){

  predDistr <- getPredictiveDistribution(parMatrix = parMatrix,
                                         model = model,
                                         numSamples = numSamples)
  # apply error to predictions
  for (i in 1:nrow(predDistr)){
    predDistr[i,] = error(mean = predDistr[i,], par = parMatrix[i,]) 
  }
  
  fittedPars = apply(parMatrix, 2, median)
  fittedPredictedResponse = model(fittedPars)
  
  dh = DHARMa::createDHARMa(simulatedResponse = t(predDistr),
                            observedResponse = observed,
                            fittedPredictedResponse = fittedPredictedResponse)
  if (plot == TRUE) {
    DHARMa::plotResiduals(dh)
  }

  return(dh)
}