R/plot_per_mod.R

#' @title Plotting function for periodic curves model
#' 
#' @description 
#' Plots graph of periodic curves with confidence intervals.
#' Data should be included in the model.
#' 
#' @param Model The built model 
#' @param XvarName Name of the x variable in the dataset (column name)
#' @param Ylab Label on vertical (y) axis
#' @param Xlab Label on horizontal (x) axis
#' @param Xlim Limits of x axis
#' @param Ylim Limits of y axis
#' @param Knots Locations of knots of the splines
#' @param Title Title of the plot
#' @param Vlines Where to plot vertical lines
#' @param Hlines Where to plot horizontal lines
#' @param Cex.lab Character expansion (aka "size of font") for the labels
#' @param Cex.main Character expansion for main text
#' @param Cex.axis Character expansion for the axis text
#' @param Axes Plot axes
#' @param Add Add to existing plot
#' @param Col Color of the plotted lines
#' @param PlotCI Plot confidence intervals
#' @param Xmin The min X of data to be predicted 
#' @param Xmax The max X of data to be predicted 
#' @param nPoints Number of points to use on the x axis
#' 
#' @import graphics stats
#' @importFrom graphics abline
#' @importFrom stats predict predict.glm na.omit
#' 
#' @export
plot_per_mod <- function(Model,
                         XvarName,
                         Ylab="Response",
                         Xlab="Covariate",
                         Ylim=NULL,
                         Xlim=NULL,
                         Xmin=NULL,
                         Xmax=NULL,
                         Knots=NULL,
                         Title=NULL,
                         Vlines=NULL,
                         Hlines=NULL,
                         Cex.lab=NULL,
                         Cex.main=NULL,
                         Cex.axis=NULL, 
                         Axes=TRUE,
                         Add=FALSE,
                         Col="black",
                         PlotCI=TRUE,
                         nPoints=100) {
  
  # NOTE: depends on how the model was built - if NOT by data=some_data, but by direct reference to
  # some_data$variable ~ someothervariable
  # then logic should be different

  # if model does not contain data, stop execution
  stopifnot(!is.null(Model$data))
  
  # load the X variable from the data in the model
  Xvar <- na.omit(Model$data[[XvarName]])
  Xvar <- Xvar[order(Xvar)]
  Xvar <- unique(Xvar)
  
  # determine where to plot vertical lines
  Intervals <- c(0.0001, 0.001, .01, .1, 1, 10, 100, 1000, 10000, 100000)
  my.range <- diff(range(Xvar, na.rm=TRUE))
  if(is.null(Vlines)){ # if vertical line are not user defined
    By.v <- Intervals[which(my.range/Intervals < 10)[1]]
    Vlines <- seq(-100000, 10000, by=By.v) # TODO: do this more generaly?
  }
  
  # determine where to plot horizontal lines
  Prediction <- predict(Model, type="response", se.fit = TRUE) # determine possible y
  my.range.y <- diff(range(Prediction, na.rm = TRUE))
  if(is.null(Hlines)){ # if horizontal lines are not defined
    By.h<-Intervals[which(my.range.y/Intervals < 10)[1]]
    Hlines <- seq(-100000, 10000, by=By.h)
  }
  
  # make new X values to make predictions smoother
  if (is.null(Xmin)) {Xmin <- min(Model$data[[XvarName]], na.rm = TRUE)}
  if (is.null(Xmax)) {Xmax <- max(Model$data[[XvarName]], na.rm = TRUE)}
  Xvar <- seq(Xmin, Xmax, length.out = nPoints ) # make a sequence of Xses to plot smoothly
  NewData <- data.frame((Xvar))
  colnames(NewData) <- c(XvarName)
  Prediction <- predict.glm(Model, type = "response", se.fit = TRUE, newdata = NewData) 

  # plot main curve
  if(Add==FALSE) {
    matplot( rbind(Xvar, Xvar, Xvar), # make a new plot
             rbind(Prediction$fit, Prediction$fit-1.96*Prediction$se, Prediction$fit+1.96*Prediction$se),
             pch=rep(1, length(Prediction$fit)*3),
             type="n",
             xlab=Xlab, ylab=Ylab, xlim=Xlim, ylim=Ylim,
             main=Title, cex.lab=Cex.lab,
             axes=Axes, cex.main=Cex.main, cex.axis=Cex.axis, col=Col) 
    } 
  # plot the lines in any case
  lines( Xvar[order(Xvar)], Prediction$fit[order(Xvar)], col=Col) # just add line to existing plot
      
  # plot confidence intervals
  if (PlotCI==TRUE) {lines( Xvar[order(Xvar)], (Prediction$fit+1.96*Prediction$se)[order(Xvar)], lty=2, col=Col)}
  if (PlotCI==TRUE) {lines( Xvar[order(Xvar)], (Prediction$fit-1.96*Prediction$se)[order(Xvar)], lty=2, col=Col)}
  
  # plot horizontal and vertical lines
  if(!is.null(Hlines) | !is.null(Vlines)) abline(h=Hlines, v= Vlines, lty=3, col="light grey") else grid()
  if(!is.null(Knots)) axis(1, at=Knots, line=2, cex.axis=.65)
  
}
crtahlin/peRiodic documentation built on May 14, 2019, 12:05 p.m.