R/nmScatterPlot.R

#' Generates a set of scatter plots with features that are tailored for PK/PD data generated by NONMEM.
#' @name nmScatterPlot
#' @title NONMEM scatter plot
#' @param obj An object of class NMRun, NMProblem, or data.frame. The object from which data will be plotted.
#' @param xVars Single string with name of x-axis variable.
#' @param yVars Variables to plot on the y-axes, specified as a comma separated string or character vector
#' @param bVars "Trellis" variables on which to split data.
#' @param gVars "Grouping" variable -- used to group points by colour, for legends and so on.  
#' @param iVars Variable name of the "inherent grouping" (subject identifier) variable.  This affects plots whose type is "l", "i", "o" or "t".
#' Lines will be grouped by this variable in addition to any other grouping (e.g. in gVar), but colours will not vary in this grouping
#' @param addLegend Should legends be added?
#' @param addGrid should grids be added?
#' @param addLoess should a loess smoother line be added? 
#' @param titles Plot title 
#' @param logX Should the x-axis be logged?
#' @param logY similar to logX for y axis
#' @param idLines Should reference lines of slope 1 and y-intercept 0 be added?
#' @param abLines A list of reference lines.  Each element should be a length 1 or 2 numeric vector.  If length 1, a 
#' horizontal reference line will be added at that point, otherwise the intercept and slope of the line.  If parameter does not satisfy
#' these conditions, a warning will be emitted
#' @param xLab x-axis label.  By default, will use the names of the x-axis variable
#' @param yLab similar to the above, but for y-axis label
#' @param types Plot type to use.  Allowed types are "p" (for standard points), "l" (lines connected by the variable
#' specified in "iVars", "i" for points labelled by "iVars"), "o" (lines and points), and "t" for labels connected by lines grouped by "iVars"
#' @param overlaid Logical flag. If TRUE, for each fixed x, the y variables will be overlaid onto a single plot
#' @param layout A length 2 vector which is passed to the layout argument of xyplot
#' @param maxPanels Maximum number of panels that should appear on each page of a graph.
#' @param maxTLevels If a single numeric (or string), the maximum number of levels that a "by" variable can have before it is binned.
#'                      If a character vector or a vector of length greater than one, the explicit breakpoints.
#' @param xRotAngle Angle by which to rotate the x-axis tick marks
#' @param equalAxisScales single logical
#' @param problemNum Number of the problem (applicable to NMRun class only)
#' @param subProblems Number of the simulation subproblems to use (applicable to the NMSim* classes only)
#' @param uniqueX Number of unique values for the X axis.  If number of uniques is less than or equal to this, a boxplot is created instead
#' @param graphParams A full list of graphical parameter settings, such as that returned by \code{\link{getAllGraphParams}}.
#' Will override global settings. By default, these will be the global settings
#' @param xAxisScaleRelations X-axis scale relations when panels are displayed. One of \code{"same"}, \code{"free"} or \code{"sliced"}.
#' @param yAxisScaleRelations Y-axis scale relations when panels are displayed. One of \code{"same"}, \code{"free"} or \code{"sliced"}. 
#' @param \dots Additional variables passed to the xyplot function 
#' @return An object of class \code{multiTrellis} holding the plot
#' @examples
#' \dontrun{ 
#'  Theoph.df <- as.data.frame(Theoph)
#'  nmScatterPlot(Theoph.df, xVar = "Time", yVar = "conc", 
#'    iVar = "Subject", type = "l", title = "Theophiline", yLab = "conc")
#'  Indometh.df <- as.data.frame(Indometh)
#'  nmScatterPlot(Indometh.df, xVar = "time", yVar = "conc", bVar = "Subject")
#' }
#' @author Mango Solutions
#' @keywords hplot
#' @exportMethod nmScatterPlot

LOGAXISLIMIT <- 0.01

nmScatterPlot <- function(obj, xVars, yVars, bVars = NULL, gVars = NULL, iVars = "ID", 
		addLegend = FALSE, addGrid = TRUE, addLoess = FALSE, titles = "", 
		logX = NULL, logY = NULL, 
		idLines = FALSE, abLines = NULL, xLab = NULL, yLab = NULL, 
		types = "p", overlaid = FALSE, layout = NULL, maxPanels = NULL, maxTLevels = Inf, 
		xRotAngle = 0, equalAxisScales = FALSE,
		problemNum = 1, subProblems = 1, uniqueX = 7, graphParams = getAllGraphParams(),
		xAxisScaleRelations = c("same", "free", "sliced"), yAxisScaleRelations = c("same", "free", "sliced"), ...)
{
	RNMGraphicsStop("Not implemented for this class yet!")
}

setGeneric("nmScatterPlot")

nmScatterPlot.NMRun <- function(obj, xVars, yVars, bVars = NULL, gVars = NULL, iVars = "ID", 
		addLegend = FALSE, addGrid = TRUE, addLoess = FALSE, titles = "", 
		logX = FALSE, logY = FALSE,
		idLines = FALSE,  abLines = NULL, xLab = NULL, yLab = NULL, 
		types = "p", overlaid = FALSE, layout = NULL, maxPanels = NULL, maxTLevels = Inf, 
		xRotAngle = 0, equalAxisScales = FALSE, 
		problemNum = 1, subProblems = 1, uniqueX = 7, graphParams = getAllGraphParams(), 
		xAxisScaleRelations = c("same", "free", "sliced"), yAxisScaleRelations = c("same", "free", "sliced"), ...)
{
	prob <- getProblem(obj, problemNum)
	x <- as.list(match.call())
	x$obj <- prob
	do.call(nmScatterPlot, x[-1])
}

setMethod("nmScatterPlot", signature(obj = "NMRun"), nmScatterPlot.NMRun)


nmScatterPlot.NMProblem <- function(obj, xVars, yVars, bVars = NULL, gVars = NULL, iVars = "ID", 
		addLegend = FALSE, addGrid = TRUE, addLoess = FALSE, titles = "", 
		logX = FALSE, logY = FALSE,
		idLines = FALSE,  abLines = NULL, xLab = NULL, yLab = NULL, 
		types = "p", overlaid = FALSE, layout = NULL, maxPanels = NULL, maxTLevels = Inf, 
		xRotAngle = 0, equalAxisScales = FALSE, 
		problemNum = 1, subProblems = 1, uniqueX = 7, graphParams = getAllGraphParams(), 
		xAxisScaleRelations = c("same", "free", "sliced"), yAxisScaleRelations = c("same", "free", "sliced"), ...)
{
	
	dataSet <- applyGraphSubset(nmData(obj, subProblemNum = subProblems), graphSubset(obj))
	x <- as.list(match.call())
	x$obj <- dataSet
	
	do.call(nmScatterPlot, x[-1])
	
}

nmScatterPlot.data.frame <- function(obj, xVars, yVars, bVars = NULL, gVars = NULL, iVars = "ID", 
		addLegend = FALSE, addGrid = TRUE, addLoess = FALSE, titles = "", 
		logX = FALSE, logY = FALSE,
		idLines = FALSE,  abLines = NULL, xLab = NULL, yLab = NULL, 
		types = "p", overlaid = FALSE, layout = NULL, maxPanels = NULL, maxTLevels = Inf, 
		xRotAngle = 0, equalAxisScales = FALSE, 
		problemNum = 1, subProblems = 1, uniqueX = 7, graphParams = getAllGraphParams(), 
		xAxisScaleRelations = c("same", "free", "sliced"), yAxisScaleRelations = c("same", "free", "sliced"), ...)
{
	
	if(length(maxPanels) > 0) layout <- NULL
	# ensure that maxPanels is numeric, even if empty
	else maxPanels <- numeric(0)
	# extract xVars - currently only one is allowed
	## agott (15/11/13)
	## include removeEmpty option to prevent empty string errors
	xVars <- CSLtoVector(xVars, removeEmpty = TRUE)
	RNMGraphicsStopifnot(length(xVars) == 1, msg = "Multiple x variables are not allowed at the moment\n")
	
	## agott (15/11/13)
	## include removeEmpty option to prevent empty string errors
	yVars <- CSLtoVector(yVars, removeEmpty = TRUE)
	RNMGraphicsStopifnot(length(yVars) >= 1, msg = "At least one y variable must be selected\n")
	
	
	# extract the desired subset
	dataSet <- applyGraphSubset(obj, graphSubset(obj))
	
	# check iVars
	iVars <- if(!is.null(iVars)) { CSLtoVector(iVars) }
	
	iVars <- iVars[iVars %in% colnames(dataSet)]
	
	if (length(iVars) == 0) { iVars <- NULL }
	
	# Check to see whether there are less than "uniqueX" values on the X axis.  If there are, switch instead to a boxplot
	RNMGraphicsStopifnot(length(uniqueX) == 1 && uniqueX %in% c(0:100, Inf), msg = "'uniqueX' input must be an integer in range [1-100], or Inf\n")
	xUni <- length(unique(dataSet[[xVars]]))
	if (xUni <= uniqueX) {
		# Need to switch to use the boxplot function instead
		return(nmBoxPlot(obj = dataSet, contVar = yVars, factVar = xVars, 
						bVars = bVars, iVar = iVars, titles = titles, xLabs = xLab, yLabs = yLab, 
						xRotAngle = xRotAngle, layout = layout, maxPanels = maxPanels, 
						maxTLevels = maxTLevels, yAxisScaleRelations = yAxisScaleRelations, xAxisScaleRelations = xAxisScaleRelations, 
						problemNum = problemNum, 
						subProblems = subProblems, ...)
		)
	}
	
	# if the y-axis variables should be overlaid (rather than displayed side-by-side on lattice panels),
	# go straight to the .overlaidScatter function and return the result
	if(overlaid && length(yVars) > 1)
	{
		return(.overlaidScatter(obj =dataSet, xVars = xVars, yVars = yVars, bVars = bVars, 
						gVars = gVars, iVars = iVars, 
						addLegend = addLegend, addGrid = addGrid,
						addLoess = addLoess, titles = titles, logX = logX, 
						logY = logY, idLines = idLines, xLab = xLab, abLines = abLines,
						yLab = yLab, types = types, equalAxisScales = equalAxisScales,
						maxPanels = maxPanels, layout = layout, maxTLevels = maxTLevels, xRotAngle = xRotAngle,
						graphParams = graphParams, xAxisScaleRelations = match.arg(xAxisScaleRelations),
						yAxisScaleRelations = match.arg(yAxisScaleRelations), ...))
	}
	
	# take all combinations of x variables against y variables
	# NOTE: this logic is outdated and is related to the fact that lattice's extended formula functionality
	# was not being used to plot multiple x/y variables, but is now.  numCombos should always be 1 at the
	# moment
	# It is benign at the moment but confusing.
	# TODO: remove this
	
	varCombos <- varComboMatrix(xVars, yVars)
	numCombos <- nrow(varCombos)
	
	gVars <- if(!is.null(gVars)) CSLtoVector(gVars) else "NULL"
	# At the moment, only one gVar handled
	if (!all(length(gVars) == 1)) { warning("gVars must be length 1; using", gVars[1], "only") } 
	gVars <- gVars[1]
	titles <- rep(titles,numCombos)
	
	# assign and y labels, taking care to handle the case where they are missing
	# TODO: There is a quirk here in that xLab and yLab _cannot_ be a comma seperated list since empty strings
	# are currently deleted by CSLtoVector.  A workaround is needed for this, and documentation should
	# make clear thet CSLs are not allowed.  In any case, this is irrelevant since numCombos can only be 1
	if(!is.null(xLab))
		xLab <- rep(xLab, length.out = numCombos)
	else
		xLab <- varCombos[,1]
	if(!is.null(yLab))
		yLab <- rep(yLab, length.out = numCombos)
	else
		yLab <- varCombos[,2]
	
	plotFormulas <- apply(cbind(varCombos[,2], varCombos[,1] ), 1, paste, collapse = "~")
	# if there are "by variables", adjust the plotting formulas passed to xyplot
	if(!is.null(bVars))
	{
		bVars <- CSLtoVector(bVars)
		# bin the by-variables as necessary
		temp <- processTrellis(dataSet, bVars, maxLevels = maxTLevels, exemptColumns = iVars)
		bVars <- temp$columns
		# coerce each "by" variable to a factor
		dataSet <- coerceToFactors(temp$data, bVars)
		# post-process the plot formulas by adding the trellis variables
		plotFormulas <- sapply(seq_len(numCombos), function(i) paste(plotFormulas[i], paste(bVars, collapse = "+"), sep = "|"))
	}
	plotList <- vector(mode = "list", length = numCombos)
	
	par.settings <- mapTopar.settings(graphParams)
	stripfn <- getStripFun()
	
	if(addLegend & gVars != "NULL")
	{
		plotKey <- scatterPlotKey(getVarLabel(gVars), dataSet[[gVars]], type = types, sortLevels = TRUE,
				graphParams = graphParams)
		# list(title = getVarDescription(gVars)$Label, rows = 10, cex=.7, space="right")
	}
	else plotKey <- NULL
	
	# get idLabels.  Note that these will have to be repeated if there is more than one yvar
	pastes <- function(...) { paste(..., sep = "/") }
	
	idLabels <- if(is.null(iVars)) { NULL } else { rep(do.call("pastes", dataSet[iVars]), times = length(yVars)) }
	
	# initialize x axis label rotation to xRotAngle		
	scales <- list(x = list(rot = xRotAngle), y = list())
	
	# log axes if necessary
	# TODO: handle errors in case xVars or yVars is a factor
	# note that logging is only done against the first y variable if there is more than one
	if(logX) scales$x <- c(scales$x, list(log = "e", at = pretty(dataSet[[xVars[1]]])))
	if(logY) scales$y <- c(scales$y, list(log = "e", at = pretty(dataSet[[yVars[1]]])))
	
	# Here we check for integer data on the x-axis and y-axis.  If so, use integer tick marks.  
	# We remove missing data first just in case prior to performing this.  Note that the y-axis
	# is only handled if there is a single y variable, and this code is only executed if equal axis-scales
	# is set to false.
	
	if(!equalAxisScales)
	{
		
		if(!is.factor(dataSet[[xVars[1]]]) && all(na.omit(dataSet[[xVars[1]]] == round(dataSet[[xVars[1]]]))))
		{
			if(length(unique(dataSet[[xVars[1]]])) > 1)
				scales$x$at <- unique(round(pretty(dataSet[[xVars[1]]])))
			else 
				scales$x$at <- unique(dataSet[[xVars[1]]])
		}
		if(length(yVars) == 1 &&  !is.factor(dataSet[[yVars]]) &&	
				all(na.omit(dataSet[[yVars]] == round(dataSet[[yVars]])))) 
		{
			if(length(unique(dataSet[[yVars]])) > 1)
				scales$y$at <- unique(round(pretty(dataSet[[yVars]])))
			else 
				scales$y$at <- unique(dataSet[[yVars]])	
			
		}
		
	}
	# if we are forcing the axis scales to be identical, we must pad out the range of the data since otherwise
	# clipping occurs.  However, care must be taken if logged axes were requested	
	else
	{
		
		if(logX | logY) 
		{
			
			# Here we use a minimum to avoid 0 or negative axes with logging
			
			scales$limits <- padLimits(range(unlist(dataSet[,c(xVars[1],yVars)]), na.rm=TRUE), min = LOGAXISLIMIT)
			RNMGraphicsWarning("Matching scales and logging requested - axis limits will be bounded from below by 0.01 which may caused data to be clipped\n")			
		}
		else
			scales$limits <- padLimits(range(unlist(dataSet[,c(xVars[1],yVars)]), na.rm=TRUE))
		
	}
	
	# only apply the yAxisScaleRelations if there is reason to 
	
	if(length(yVars) > 1 || length(bVars) > 0) 
	{
		scales$y$relation <- match.arg(yAxisScaleRelations)      
	}
	
	# only apply the xAxisScaleRelations if there is reason to
	
	if(length(bVars) > 0)
	{        
		scales$x$relation <- match.arg(xAxisScaleRelations)      
	}
	
	featuresToAdd <- c("grid" = addGrid, "loess" = addLoess, "idLine" = idLines)
	plotList[[1]] <- 
			with(graphParams, 
					xyplot(as.formula(plotFormulas[[1]]), groups = eval(parse(text = gVars)),  
							data = dataSet, panel = panel.nmScatterPlot, outer = TRUE, stack = FALSE, 
							featuresToAdd = featuresToAdd, key = plotKey, main = titles[[1]], idLabels = idLabels, 
							xlab = xLab[1],	ylab = yLab[1],	type = types[1], scales =scales,
							par.settings = par.settings, strip = stripfn, layout = layout, 
							multiYVar = length(yVars) > 1, logX = logX, logY = logY, 
							graphParams = graphParams, abLines = abLines, ...)
			) # end with
	gridDims <- stdGridDims(numCombos,3 )
	result <- multiTrellis(plotList, gridDims, maxPanels = maxPanels)
	result
}


setMethod("nmScatterPlot", signature(obj = "data.frame"), nmScatterPlot.data.frame)

setMethod("nmScatterPlot", signature(obj = "NMProblem"), nmScatterPlot.NMProblem)


# nmScatterPlot panel function - this is used only when overlaid = FALSE. 
# @name panel.nmScatterPlot
# @title nmScatterPlot panel function (1 of 2)
# @param x (usual)
# @param y (usual)
# @param subscripts (usual) 
# @param featuresToAdd named logical vector of features to add (grid, loess, idLine) 
# @param idLabels Identifier labels 
# @param type one of "p", "i", "l", "o", "t"
# @param groups (usual)
# @param multiYVar currently unused
# @param ... additional parameters to panel.xyplot / panel.superpose
# @param graphParams Full list of RNMGraphics graphical settings
# @return none
# @keywords panel

panel.nmScatterPlot <- function(x, y, subscripts = seq_along(x), featuresToAdd =  c("grid" = FALSE, "loess" = FALSE, "idLine" = FALSE), 
		idLabels = NULL, type = c("p", "i", "l", "o","t"), groups = NULL, multiYVar = FALSE, 
		graphParams, abLines = NULL, ...)
{
	type <- match.arg(type)
	
	# call panel.xyplot to setup first
	panel.xyplot(x, y, type = "n",...)
	if(featuresToAdd["grid"])
	{
		gridOpts <- graphParams$"grid"
		panel.grid(h = -1, v = -1, col = gridOpts$col, alpha = gridOpts$alpha, lty = gridOpts$lty, 
				lwd = gridOpts$lwd)
	}
	
	if(featuresToAdd["idLine"])
		panel.abline(a = 0, b = 1)
	
	# grab "superpose.line" styles
	superpose.line.col <- graphParams$superpose.line$col
	
	
	if(type == "p")
	{
		# points only
		panel.xyplot(x = x, y = y, subscripts = subscripts, type = type, groups,
				...)
		
	}
	
	else if(type == "l")
	{		
		groupInfo <- subjectGrouping(idLabels, groups, graphParams$superpose.line)
		if(is.null(groups))
		{
			plot.line <- graphParams$plot.line
			panel.superpose(x, y, groups = idLabels, type = type, subscripts = subscripts, 
					col.line = plot.line$col, col.lwd = plot.line$lwd, lty = plot.line$lty, ...)
		}
		else
		{
			panel.superpose(x, y, groups = groupInfo$grouping, type = type, subscripts = subscripts, 
					col.line = groupInfo$elements$col, lty = groupInfo$elements$lty, lwd = groupInfo$elements$lwd, ...)
		}
		
	}
	# lines + points connecting subjects
	else if(type == "o")
	{
		RNMGraphicsStopifnot(!is.null(idLabels))
		groupInfo <- subjectGrouping(idLabels, groups, graphParams$superpose.line)
		
		if(is.null(groups))
		{
			plot.line <- graphParams$plot.line
			panel.xyplot(x, y, type = "p", subscripts = subscripts, ...)
			panel.superpose(x, y, groups = idLabels, 
					type = "l", subscripts = subscripts, col = plot.line$col, 
					lty = plot.line$lty, lwd = plot.line$lwd, ...)
		}
		else
		{
			panel.superpose(x, y, groups = groups, type = "p", subscripts = subscripts, ...)
			panel.superpose(x, y, groups = groupInfo$grouping, type = "l", subscripts = subscripts, 
					col.line = groupInfo$elements$col, lty = groupInfo$elements$lty, lwd = groupInfo$elements$lwd, ...)
		}
	}
	# subject identifiers
	else if(type == "i")
	{
		# TODO: allow use of "gVar"
		RNMGraphicsStopifnot(!is.null(idLabels))
		if(!is.null(groups))
		{
			textopt <- graphParams$"superpose.text"
			groupInfo <- subjectGrouping(idLabels, groups, textopt, expandElements = TRUE)
			ltext(x, y, idLabels[subscripts], col = groupInfo$elements$col[subscripts] , cex = groupInfo$elements$cex[subscripts] , ...)
		}
		else
		{
			textopt <- graphParams$"plot.text"
			ltext(x, y, idLabels[subscripts], col = textopt$col , cex = textopt$cex , ...)
		}
	}
	# lines connecting identifiers
	else if(type == "t")
	{
		RNMGraphicsStopifnot(!is.null(idLabels))
		
		if(!is.null(groups)) 
		{
			textopt <- graphParams$"superpose.text"
			groupInfo <- subjectGrouping(idLabels, groups, graphParams$superpose.line)
			groupInfo2 <- subjectGrouping(idLabels, groups, textopt, expandElements = TRUE)
			ltext(x, y, idLabels[subscripts], col = groupInfo2$elements$col[subscripts] ,
					groupInfo2$elements$cex[subscripts], ...)		
			panel.superpose(x, y, groups = groupInfo$grouping, type = "l", 
					subscripts = subscripts, col.line = groupInfo$elements$col,
					lty = groupInfo$elements$lty , lwd = groupInfo$elements$lwd,  ...)
		}
		else
		{
			plot.line <- graphParams$"plot.line"
			textopt <- graphParams$"plot.text"
			ltext(x, y, idLabels[subscripts], col = textopt$col , cex = textopt$cex , ...)
			groups <- idLabels
			panel.superpose(x = x, y = y, subscripts = subscripts, type = "l", 
					groups = idLabels, col.line = plot.line$col, lty = plot.line$lty, lwd = plot.line$lwd, ...) 
		}
	}
	
	if(featuresToAdd["loess"])
	{
		loessOpts <- graphParams$"loess.line"
		# implement a try-catch just in case loess curve fails to compute correctly
		tryLoess <- try(panel.loess(x,y, col = loessOpts$col, lwd = loessOpts$lwd), silent = TRUE)
		if(inherits(tryLoess, "try-error"))
			RNMGraphicsWarning("Failed to calculate loess curve, omitting from this panel\n")
	}
	
	if(!is.null(abLines))
	{
		# check that abLines is a list of vectors of length 1 or 2
		if(!(is(abLines, "list") & all(sapply(abLines, length) %in% c(1,2))))
		{
			RNMGraphicsWarning("Reference line parameter (abLines) is invalid")
			return()
		}
		
		for(x in abLines)
		{
			tryCatch(panel.abline(x), 
					error = function(e) RNMGraphicsWarning(paste("Failed to plot a reference line, error message is:", e)))
		}
		
	}
}
MangoTheCat/RNMGraphics documentation built on May 8, 2019, 3:51 p.m.