# R/plot.spv.R In vdg: Variance Dispersion Graphs and Fraction of Design Space Plots

#' Plot VDGs or FDS plots
#'
#' Produce Variance Dispersion Graphs and/or Fraction of Design Space plots for
#' experimental designs. There are methods for the S3 classes \code{spv},
#' \code{spvlist}, \code{spvforlist} and \code{spvlistforlist} -- see
#'
#' @aliases plot.spv plot.spvlist plot.spvforlist plot.spvlistforlist
#' @param x an object of type \code{spv} for a single experimental design or an
#' object of type \code{spvlist} for multiple designs.
#' @param which either a numeric vector of integers or a character vector
#' indicating which plots to produce. The possible plots are:
#' \describe{
#' \item{\code{1} or \code{"fds"}}{A (variance ratio) FDS plot}
#' \item{\code{2} or \code{"vdgsim"}}{A VDG with only the simulated prediction variance points plotted}
#' \item{\code{3} or \code{"vdgquantile"}}{A VDG with only the quantile regression lines corresponding to \code{tau} shown}
#' \item{\code{4} or \code{"vdgboth"}}{A combination of \code{2} and \code{3}}
#' \item{\code{5} or \code{"boxplots"}}{Parallel boxplots of the prediction variance}
#' }
#' @param np scalar; the number of points to use for calculating the fraction of design space criterion.
#' @param alpha the alpha transparency coefficient for the plots
#' @param points.colour colour for points in scatterplot of SPV against the radius
#' @param points.size size for points in scatterplot of SPV against the radius
#' @param tau the tau parameter for \code{\link[quantreg]{rq}} (quantile
#' regression)
#' @param radii either a numeric vector containing the radii to use for
#' calculating the mean spherical SPV over the spherical design space, or an integer
#' (length one vector) giving the number of radii to use for calculationg
#' the mean spherical SPV. If missing, the mean spherical SPV is not used.
#' @param hexbin logical indicating whether hexagonal binning should be used to display
#' density instead of colour transparency
#' @param bins argument passed to \code{\link{stat_binhex}} to determine the
#' number of hexagons used for binning.
#' @param VRFDS logical indicating whether to construct a variance ratio FDS plot or not (only for class \code{spvlist}). The
#' first design is used as reference design in case of \code{VRFDS} is \code{TRUE}
#' @param df degrees-of-freedom parameter passed to \code{\link{bs}}
#' @param lines.size line size passed to \code{\link[ggplot2]{geom_line}}
#' @param origin numeric vector specifying the origin of the design space
#' @param method optional; passed to \code{\link[proxy]{dist}} to overwrite
#' defaults of "Euclidean" for spherical regions or "supremum" for cubiodal
#' regions
#' @param arrange Logical indicating whether to return a single graphical object arranging the
#' resulting plots in a single plot window via \code{\link{grid.arrange}}, or whether to return the
#' list of graphical objects containing the plots.
#' @return Returns a list of \code{\link{ggplot}} graphical objects (or grobs) with names corresponding
#' to the character version of \code{which}. These plot objects can be manipulated by changing plot
#' aesthetics and theme elements.
#' @keywords hplot
#' @author Pieter C. Schoonees
#' @references
#' Pieter C. Schoonees, Niel J. le Roux, Roelof L.J. Coetzer (2016). Flexible Graphical Assessment of
#' Experimental Designs in R: The vdg Package. \emph{Journal of Statistical Software}, 74(3), 1-22.
#' \doi{10.18637/jss.v074.i03}.
#' @method plot spv
#' @export
#' @import ggplot2
#' @import quantreg
#' @importFrom proxy dist
#' @importFrom splines bs
#' @importFrom gridExtra grid.arrange
#' @examples
#'
#' # Single design (class 'spv')
#' # Larger n should be used in actual cases
#' library(rsm)
#' bbd3 <- as.data.frame(bbd(3)[,3:5])
#' colnames(bbd3) <- paste0("x", 1:3)
#' quad.3f <- formula(~ x1*x2*x3 - x1:x2:x3 + I(x1^2) + I(x2^2) + I(x3^2))
#' set.seed(1234)
#' out <- spv(n = 1000, design = bbd3, type = "spherical", formula = quad.3f)
#' out
#' plot(out)
#'
#' # List of designs (class 'spvlist')
#' \dontrun{
#' data(SCDH5); data(SCDDL5)
#' des.list <- list(SCDH5 = SCDH5, SCDDL5 = SCDDL5)
#' quad.5f <- formula(~ x1 + x2 + x3 + x4 + x5 + x1:x2 + x1:x3 + x1:x4 + x1:x5
#'                     + x2:x3 + x2:x4 + x2:x5 + x3:x4 + x3:x5 + x4:x5
#'                    + I(x1^2) + I(x2^2) + I(x3^2) + I(x4^2) + I(x5^2))
#' out2 <- spv(n = 500, design = des.list, type = "spherical", formula = quad.5f)
#' out2
#' plot(out2)
#' }
#'
#' # List of formulae (class 'spvforlist')
#' \dontrun{
#' fact3 <- expand.grid(x1 = c(-1,1), x2 = c(-1, 1), x3 = c(-1,1))
#' lin.3f <- formula(~ x1 + x2 + x3)
#' int.3f <- formula(~ (x1+x2+x3)^2)
#' set.seed(4312)
#' out3 <- spv(n = 500, design = fact3, type = "cuboidal",
#'              formula = list(linear = lin.3f, interaction = int.3f))
#' out3
#' plot(out3)
#' }
#'
#' # List of formulae and designs (class 'spvlistforlist')
#' \dontrun{
#' fact3.n <- rbind(fact3, 0, 0, 0)
#' set.seed(4312)
#' out4 <- spv(n = 200, design = list(factorial = fact3, factorial.with.cntr = fact3.n),
#'              type = "cuboidal", formula = list(linear = lin.3f, interaction = int.3f))
#' out4
#' plot(out4)
#' }
plot.spv <- function (x, which = c("fds", "vdgsim", "vdgquantile", "vdgboth", "boxplots"),
np = 50, alpha = 7/sqrt(length(x$spv)), points.colour = "#39BEB1", points.size = 2, tau = c(0.05, 0.95), radii = 21, hexbin = FALSE, bins = 80, df = 10, lines.size = 1, origin = rep(0, ncol(x$sample)), method,
arrange = FALSE, ...) {
# Avoid global variable notes for R CMD check and ggplot2
Radius <- SPV <- Fraction <- Location <- NULL

# Handle which depending on whether it is numeric or character (gets transformed to numeric)
pnms <- c("fds", "vdgsim", "vdgquantile", "vdgboth", "boxplots")
show <- rep(FALSE, 5)
if (is.character(which)) {
which <- match.arg(which, several.ok = TRUE)
which <- sort(match(which, pnms))
}
if (!is.numeric(which)) stop("Argument 'which' is of incorrect type.")
show[which] <- TRUE
type <- x$type if (x$at && show[1L]){
show[1L] <- FALSE
which <- which[!(which %in% 1L)]
message("Plot 1 = 'fds' cannot be produced: 'at' is TRUE (inaccurate FDS plot)")
}

add.meanspv <- x$type == "spherical" && !is.null(radii) if (is.null(tau) & !add.meanspv){ if(any(3L:4L %in% which)) message("Plots 3 = 'vdgquantile' and/or 4 = 'vdgboth' cannot be produced: 'tau' is NULL and mean SPV not requested/possible") show[3L:4L] <- FALSE which <- which[!(which %in% 3L:4L)] } if (!x$at && show[5L]){
show[5L] <- FALSE
which <- which[!(which %in% 5L)]
message("Plot 5 = 'boxplots' cannot be produced: 'at' is FALSE")
}

pnms <- pnms[show]

if (missing(method)) method <- switch(type, spherical = "Euclidean", cuboidal = "supremum",
lhs = "supremum", mlhs = "supremum", slhs = "supremum",
rslhs = "supremum")
xvec <- proxy::dist(x$sample, matrix(origin, nrow = 1, ncol = ncol(x$sample)), method = method, ...)
method <- attr(xvec, "method")
xvec <- as.numeric(xvec)

mspv <- meanspv(formula = x$formula, radii = radii, FtF.inv = x$FtF.inv,
n = ifelse(x$unscaled, 1, x$ndes))
tmp3 <- as.data.frame(mspv)
tmp3$Location <- "Mean" } if (any(show[-1L])) tmp1 <- data.frame(Radius = xvec, SPV = x$spv)

if (show[1L]){
maxmin <- range(x$spv) pts <- 0:np/np tmp2 <- data.frame(Fraction = pts, SPV = quantile(x$spv, probs = pts, type = 1))
}

if (show[3L] || show[4L]){
if (!is.null(tau)){
pts <- seq(from = min(tmp1$Radius), to = max(tmp1$Radius), length = np)
fits <- lapply(tau, function(x, data) quantreg::rq(SPV ~ bs(Radius, df = df), tau = x,
data = data), data = tmp1)
newdf <- data.frame(Radius = rep(pts, length(tau)), SPV = as.numeric(
sapply(fits, predict, newdata = data.frame(Radius = pts))),
Location = rep(paste("tau =", tau), each = np))
if (exists("tmp3", inherits = FALSE)) tmp3 <- rbind(tmp3, newdf)
else tmp3 <- newdf
}
}

if (exists("tmp3", inherits = FALSE)) tmp3$Location <- as.factor(tmp3$Location)

if (show[1L]){
plot1 <- ggplot(tmp2, aes(x = Fraction, y = SPV)) + ggtitle("Fraction of Design Space Plot") +
xlab("Fraction of Design Space") +
geom_line(size = lines.size, colour = points.colour) +
theme(plot.title = element_text(vjust = 1))
}

if (show[2L]) {
plot2 <- ggplot(tmp1, aes(x = Radius, y = SPV)) + ggtitle("Variance Dispersion Graph") +
geom_point(alpha = alpha, colour = points.colour, size = points.size) +
theme(plot.title = element_text(vjust = 1)) +
xlab(paste0("Distance to Origin (", method,")"))
if (hexbin) {
plot2 <- plot2 + geom_hex(bins = bins) +
scale_fill_gradientn(colours = rev(topo.colors(5)[-(4:5)]), name = "Frequency", na.value = NA)
}
}

if (show[3L]) {
plot3 <- ggplot(tmp1, aes(x = Radius, y = SPV)) + ggtitle("Variance Dispersion Graph") +
theme(plot.title = element_text(vjust = 1), legend.text.align = 0.5)
if(hexbin){
plot3 <- plot3 + geom_hex(bins = bins) +
scale_fill_gradientn(colours = rev(topo.colors(5)[-(4:5)]), name = "Frequency", na.value = NA)
}
plot3 <- plot3 + geom_line(mapping = aes(x = Radius, y = SPV, linetype = Location), data = tmp3,
size = lines.size) + xlab(paste0("Distance to Origin (", method,")"))
}

if (show[4L]) {
plot4 <- ggplot(tmp1, aes(x = Radius, y = SPV)) + ggtitle("Variance Dispersion Graph") +
geom_point(alpha = alpha, colour = points.colour, size = points.size) +
theme(plot.title = element_text(vjust = 1), legend.text.align = 0.5)
if (hexbin) plot4 <- plot4 + geom_hex(bins = bins) +
scale_fill_gradientn(colours = rev(topo.colors(5)[-(4:5)]), name = "Frequency", na.value = NA)
plot4 <- plot4 + geom_line(mapping = aes(x = Radius, y = SPV, linetype = Location), data = tmp3,
size = lines.size) + xlab(paste0("Distance to Origin (", method,")"))
}

if (show[5L]) {
plot5 <- ggplot(tmp1, aes(x = Radius, y = SPV)) + ggtitle("Boxplots") +
xlab(paste0("Distance to Origin (", method,")")) +
theme(plot.title = element_text(vjust = 1), legend.text.align = 0.5)
}

if (length(which)) {
out <- mget(paste0("plot", which))
names(out) <- pnms
if(x\$unscaled) out <- lapply(out, '+', ylab("Unscaled Prediction Variance (UPV)"))
else out <- lapply(out, '+', ylab("Scaled Prediction Variance (SPV)"))
if (arrange) do.call(gridExtra::grid.arrange, out)
else return(out)
}
}


## Try the vdg package in your browser

Any scripts or data that you put into this service are public.

vdg documentation built on July 8, 2022, 1:08 a.m.