R/plots.R

In bhmevik/pls: Partial Least Squares and Principal Component Regression

### Plots for mvr objects.  Some of them also work for other
### objects, but that is not a priority.

###
### Plot method for mvr objects
###



#' @title Plot Method for MVR objects
#'
#' @description \code{plot.mvr} plots predictions, coefficients, scores, loadings, biplots,
#' correlation loadings or validation plots (RMSEP curves, etc.).
#'
#' @details The function is simply a wrapper for the underlying plot functions used to
#' make the selected plots.  See \code{\link{predplot.mvr}},
#' \code{\link{validationplot}}, \code{\link{coefplot}},
#' \code{\link{scoreplot}}, \code{\link{loadingplot}}, \code{\link{biplot.mvr}}
#' or \code{\link{corrplot}} for details.  Note that all arguments except
#' \code{x} and \code{plottype} must be named.
#'
#' @param x an object of class \code{mvr}.  The fitted model to plot.
#' @param plottype character.  What kind of plot to plot.
#' @param \dots further arguments, sent to the underlying plot functions.
#' @return \code{plot.mvr} returns whatever the underlying plot function
#' returns.
#' @author Ron Wehrens and Bjørn-Helge Mevik
#' @seealso \code{\link{mvr}}, \code{\link{predplot.mvr}},
#' \code{\link{validationplot}}, \code{\link{coefplot}},
#' \code{\link{scoreplot}}, \code{\link{loadingplot}},
#' \code{\link{biplot.mvr}}, \code{\link{corrplot}}
#' @keywords regression multivariate hplot
#' @examples
#'
#' data(yarn)
#' nir.pcr <- pcr(density ~ NIR, ncomp = 9, data = yarn, validation = "CV")
#' \dontrun{
#' plot(nir.pcr, ncomp = 5) # Plot of cross-validated predictions
#' plot(nir.pcr, "scores") # Score plot
#' plot(nir.pcr, "loadings", comps = 1:3) # The three first loadings
#' plot(nir.pcr, "coef", ncomp = 5) # Coefficients
#' plot(nir.pcr, "val") # RMSEP curves
#' plot(nir.pcr, "val", val.type = "MSEP", estimate = "CV") # CV MSEP
#' }
#'
#' @export
plot.mvr <- function(x, plottype = c("prediction", "validation",
                        "coefficients", "scores", "loadings", "biplot",
                        "correlation"),
                     ...)
{
    plottype <- match.arg(plottype)
    plotFunc <- switch(plottype,
                       prediction = predplot.mvr,
                       validation = validationplot,
                       coefficients = coefplot,
                       scores = scoreplot,
                       loadings = loadingplot,
                       biplot = biplot.mvr,
                       correlation = corrplot)
    plotFunc(x, ...)
}


###
### Scoreplot
###



#' @name scoreplot
#' @title Plots of Scores, Loadings and Correlation Loadings
#'
#' @description  Functions to make scatter plots of scores or correlation loadings, and
#' scatter or line plots of loadings.
#'
#' @details \code{plot.scores} is simply a wrapper calling \code{scoreplot}, passing all
#' arguments.  Similarly for \code{plot.loadings}.
#'
#' \code{scoreplot} is generic, currently with a default method that works for
#' matrices and any object for which \code{\link{scores}} returns a matrix.
#' The default \code{scoreplot} method makes one or more scatter plots of the
#' scores, depending on how many components are selected.  If one or two
#' components are selected, and \code{identify} is \code{TRUE}, the function
#' \code{\link{identify}} is used to interactively identify points.
#'
#' Also \code{loadingplot} is generic, with a default method that works for
#' matrices and any object where \code{\link{loadings}} returns a matrix.  If
#' \code{scatter} is \code{TRUE}, the default method works exactly like the
#' default \code{scoreplot} method.  Otherwise, it makes a lineplot of the
#' selected loading vectors, and if \code{identify} is \code{TRUE}, uses
#' \code{\link{identify}} to interactively identify points.  Also, if
#' \code{legendpos} is given, a legend is drawn at the position indicated.
#'
#' \code{corrplot} works exactly like the default \code{scoreplot} method,
#' except that at least two components must be selected.  The
#' \dQuote{correlation loadings}, i.e. the correlations between each variable
#' and the selected components (see References), are plotted as pairwise
#' scatter plots, with concentric circles of radii given by \code{radii}.  Each
#' point corresponds to a variable.  The squared distance between the point and
#' origin equals the fraction of the variance of the variable explained by the
#' components in the panel.  The default \code{radii} corresponds to 50\% and
#' 100\% explained variance.  By default, only the correlation loadings of the
#' \eqn{X} variables are plotted, but if \code{ploty} is \code{TRUE}, also the
#' \eqn{Y} correlation loadings are plotted.
#'
#' \code{scoreplot}, \code{loadingplot} and \code{corrplot} can also be called
#' through the plot method for \code{mvr} objects, by specifying
#' \code{plottype} as \code{"scores"}, \code{"loadings"} or
#' \code{"correlation"}, respectively.  See \code{\link{plot.mvr}}.
#'
#' The argument \code{labels} can be a vector of labels or one of
#' \code{"names"} and \code{"numbers"}.
#'
#' If a scatter plot is produced (i.e., \code{scoreplot}, \code{corrplot}, or
#' \code{loadingplot} with \code{scatter = TRUE}), the labels are used instead
#' of plot symbols for the points plotted.  If \code{labels} is \code{"names"}
#' or \code{"numbers"}, the row names or row numbers of the matrix (scores,
#' loadings or correlation loadings) are used.
#'
#' If a line plot is produced (i.e., \code{loadingplot}), the labels are used
#' as \eqn{x} axis labels.  If \code{labels} is \code{"names"} or
#' \code{"numbers"}, the variable names are used as labels, the difference
#' being that with \code{"numbers"}, the variable names are converted to
#' numbers, if possible.  Variable names of the forms \samp{"number"} or
#' \samp{"number text"} (where the space is optional), are handled.
#'
#' The argument \code{pretty.xlabels} is only used when \code{labels} is
#' specified for a line plot.  If \code{TRUE} (default), the code tries to use
#' a \sQuote{pretty} selection of labels.  If \code{labels} is
#' \code{"numbers"}, it also uses the numerical values of the labels for
#' horisontal spacing.  If one has excluded parts of the spectral region, one
#' might therefore want to use \code{pretty.xlabels = FALSE}.
#'
#' @aliases scoreplot scoreplot.default plot.scores loadingplot
#' loadingplot.default plot.loadings corrplot
#' @param object an object.  The fitted model.
#' @param comps integer vector.  The components to plot.
#' @param scatter logical.  Whether the loadings should be plotted as a scatter
#' instead of as lines.
#' @param labels optional.  Alternative plot labels or \eqn{x} axis labels.
#' See Details.
#' @param plotx locical.  Whether to plot the \eqn{X} correlation loadings.
#' Defaults to \code{TRUE}.
#' @param ploty locical.  Whether to plot the \eqn{Y} correlation loadings.
#' Defaults to \code{FALSE}.
#' @param radii numeric vector, giving the radii of the circles drawn in
#' \code{corrplot}.  The default radii represent 50\% and 100\% explained
#' variance of the \eqn{X} variables by the chosen components.
#' @param identify logical.  Whether to use \code{identify} to interactively
#' identify points.  See below.
#' @param type character.  What type of plot to make.  Defaults to \code{"p"}
#' (points) for scatter plots and \code{"l"} (lines) for line plots.  See
#' \code{\link{plot}} for a complete list of types (not all types are
#' possible/meaningful for all plots).
#' @param lty vector of line types (recycled as neccessary).  Line types can be
#' specified as integers or character strings (see \code{\link{par}} for the
#' details).
#' @param lwd vector of positive numbers (recycled as neccessary), giving the
#' width of the lines.
#' @param pch plot character.  A character string or a vector of single
#' characters or integers (recycled as neccessary).  See \code{\link{points}}
#' for all alternatives.
#' @param cex numeric vector of character expansion sizes (recycled as
#' neccessary) for the plotted symbols.
#' @param col character or integer vector of colors for plotted lines and
#' symbols (recycled as neccessary).  See \code{\link{par}} for the details.
#' @param legendpos Legend position.  Optional.  Ignored if \code{scatter} is
#' \code{TRUE}.  If present, a legend is drawn at the given position.  The
#' position can be specified symbolically (e.g., \code{legendpos =
#' "topright"}).  This requires >= 2.1.0.  Alternatively, the position can be
#' specified explicitly (\code{legendpos = t(c(x,y))}) or interactively
#' (\code{legendpos = \link{locator}()}).
#' @param xlab,ylab titles for \eqn{x} and \eqn{y} axes.  Typically character
#' strings, but can be expressions or lists.  See \code{\link{title}} for
#' details.
#' @param pretty.xlabels logical.  If \code{TRUE}, \code{loadingplot} tries to
#' plot the \eqn{x} labels more nicely.  See Details.
#' @param xlim optional vector of length two, with the \eqn{x} limits of the
#' plot.
#' @param x a \code{scores} or \code{loadings} object.  The scores or loadings
#' to plot.
#' @param \dots further arguments sent to the underlying plot function(s).
#' @return The functions return whatever the underlying plot function (or
#' \code{identify}) returns.
#' @note \code{\link{legend}} has many options.  If you want greater control
#' over the appearance of the legend, omit the \code{legendpos} argument and
#' call \code{legend} manually.
#'
#' Graphical parametres (such as \code{pch} and \code{cex}) can also be used
#' with \code{scoreplot} and \code{corrplot}.  They are not listed in the
#' argument list simply because they are not handled specifically in the
#' function (unlike in \code{loadingplot}), but passed directly to the
#' underlying plot functions by \code{\dots{}}.
#'
#' Tip: If the labels specified with \code{labels} are too long, they get
#' clipped at the border of the plot region.  This can be avoided by supplying
#' the graphical parameter \code{xpd = TRUE} in the plot call.
#'
#' The handling of \code{labels} and \code{pretty.xlabels} in \code{coefplot}
#' is experimental.
#' @author Ron Wehrens and Bjørn-Helge Mevik
#' @seealso \code{\link{mvr}}, \code{\link{plot.mvr}}, \code{\link{scores}},
#' \code{\link{loadings}}, \code{\link{identify}}, \code{\link{legend}}
#' @references Martens, H., Martens, M. (2000) Modified Jack-knife Estimation
#' of Parameter Uncertainty in Bilinear Modelling by Partial Least Squares
#' Regression (PLSR).  \emph{Food Quality and Preference}, \bold{11}(1--2),
#' 5--16.
#' @keywords regression multivariate hplot
#' @examples
#'
#' data(yarn)
#' mod <- plsr(density ~ NIR, ncomp = 10, data = yarn)
#' ## These three are equivalent:
#' \dontrun{
#' scoreplot(mod, comps = 1:5)
#' plot(scores(mod), comps = 1:5)
#' plot(mod, plottype = "scores", comps = 1:5)
#'
#' loadingplot(mod, comps = 1:5)
#' loadingplot(mod, comps = 1:5, legendpos = "topright") # With legend
#' loadingplot(mod, comps = 1:5, scatter = TRUE) # Plot as scatterplots
#'
#' corrplot(mod, comps = 1:2)
#' corrplot(mod, comps = 1:3)
#' }
#'
#' @export
scoreplot <- function(object, ...) UseMethod("scoreplot")

#' @rdname scoreplot
#' @export
scoreplot.default <- function(object, comps = 1:2, labels, identify = FALSE,
                              type = "p", xlab, ylab, ...)
{
    ## Check arguments
    nComps <- length(comps)
    if (nComps == 0) stop("At least one component must be selected.")
    ## Get the scores
    if (is.matrix(object)) {
        ## Assume this is already a score matrix
        S <- object[,comps, drop = FALSE]
    } else {
        ## Try to get the scores
        S <- scores(object)[,comps, drop = FALSE]
        if (is.null(S))
            stop("`", deparse(substitute(object)), "' has no scores.")
    }
    if (!missing(labels)) {
        ## Set up point labels
        if (length(labels) == 1) {
            labels <- switch(match.arg(labels, c("names", "numbers")),
                             names = rownames(S),
                             numbers = 1:nrow(S)
                             )
        }
        labels <- as.character(labels)
        type <- "n"
    }
    varlab <- compnames(object, comps, explvar = TRUE)
    if (nComps <= 2) {
        if (nComps == 1) {
            ## One component versus index
            if (missing(xlab)) xlab <- "observation"
            if (missing(ylab)) ylab <- varlab
        } else {
            ## Second component versus first
            if (missing(xlab)) xlab <- varlab[1]
            if (missing(ylab)) ylab <- varlab[2]
        }
        plot(S, xlab = xlab, ylab = ylab, type = type, ...)
        if (!missing(labels)) text(S, labels, ...)
        if (isTRUE(identify)) {
            if (!is.null(rownames(S))) {
                identify(S, labels = rownames(S))
            } else {
                identify(S)
            }
        }
    } else {
        ## Pairwise scatterplots of several components
        panel <- if (missing(labels))
            function(x, y, ...) points(x, y, type = type, ...) else
            function(x, y, ...) text(x, y, labels = labels, ...)
        pairs(S, labels = varlab, panel = panel, ...)
    }
}

## A plot method for scores:
#' @rdname scoreplot
#' @export
plot.scores <- function(x, ...) scoreplot(x, ...)


###
### Loadingplot
###

#' @rdname scoreplot
#' @export
loadingplot <- function(object, ...) UseMethod("loadingplot")

#' @rdname scoreplot
#' @export
loadingplot.default <- function(object, comps = 1:2, scatter = FALSE, labels,
                                identify = FALSE, type, lty, lwd = NULL, pch,
                                cex = NULL, col, legendpos, xlab, ylab,
                                pretty.xlabels = TRUE, xlim, ...)
{
    ## Check arguments
    nComps <- length(comps)
    if (nComps == 0) stop("At least one component must be selected.")
    if (!missing(type) &&
        (length(type) != 1 || is.na(nchar(type, "c")) || nchar(type, "c") != 1))
        stop("Invalid plot type.")
    ## Get the loadings
    if (is.matrix(object)) {
        ## Assume this is already a loading matrix
        L <- object[,comps, drop = FALSE]
    } else {
        ## Try to get the loadings:
        L <- loadings(object)[,comps, drop = FALSE]
        if (is.null(L))
            stop("`", deparse(substitute(object)), "' has no loadings.")
    }
    varlab <- compnames(object, comps, explvar = TRUE)
    if (isTRUE(scatter)) {
        ## Scatter plots
        if (missing(type)) type <- "p"
        if (!missing(labels)) {
            ## Set up point/tick mark labels
            if (length(labels) == 1) {
                labels <- switch(match.arg(labels, c("names", "numbers")),
                                 names = {
                                     if (is.null(rnames <- rownames(L))) {
                                         stop("The loadings have no row names.")
                                     } else {
                                         rnames
                                     }},
                                 numbers = 1:nrow(L)
                                 )
            }
            labels <- as.character(labels)
            type <- "n"
        }
        if (missing(lty)) lty <- NULL
        if (missing(pch)) pch <- NULL
        if (missing(col)) col <- par("col") # `NULL' means `no colour'
        if (nComps <= 2) {
            if (nComps == 1) {
                ## One component versus index
                if (missing(xlab)) xlab <- "variable"
                if (missing(ylab)) ylab <- varlab
            } else {
                ## Second component versus first
                if (missing(xlab)) xlab <- varlab[1]
                if (missing(ylab)) ylab <- varlab[2]
            }
            plot(L, xlab = xlab, ylab = ylab, type = type, lty = lty,
                 lwd = lwd, pch = pch, cex = cex, col = col, ...)
            if (!missing(labels)) text(L, labels, cex = cex, col = col, ...)
            if (isTRUE(identify))
                identify(L, labels = paste(1:nrow(L), rownames(L), sep = ": "))
        } else {
            ## Pairwise scatterplots of several components
            panel <- if (missing(labels)) {
                function(x, y, ...)
                    points(x, y, type = type, lty = lty, lwd = lwd,
                           pch = pch, col = col, ...)
            } else {
                function(x, y, ...)
                    text(x, y, labels = labels, col = col, ...)
            }
            pairs(L, labels = varlab, panel = panel, cex = cex, ...)
        }
    } else {                            # if (isTRUE(scatter))
        ## Line plots
        if (missing(type)) type <- "l"
        if (missing(lty))  lty  <- 1:nComps
        if (missing(pch))  pch  <- 1:nComps
        if (missing(col))  col  <- 1:nComps
        if (missing(xlab)) xlab <- "variable"
        if (missing(ylab)) ylab <- "loading value"
        xnum <- 1:nrow(L)
        if (missing(labels)) {
            xaxt <- par("xaxt")
        } else {
            xaxt <- "n"
            if (length(labels) == 1) {
                xnam <- rownames(L)
                switch(match.arg(labels, c("names", "numbers")),
                       names = {        # Simply use the names as is
                           labels <- xnam
                       },
                       numbers = {      # Try to use them as numbers
                           if (length(grep("^[-0-9.]+[^0-9]*$", xnam)) ==
                               length(xnam)) {
                               ## Labels are on "num+text" format
                               labels <- sub("[^0-9]*$", "", xnam)
                               if (isTRUE(pretty.xlabels)) {
                                   xnum <- as.numeric(labels)
                                   xaxt <- par("xaxt")
                               }
                           } else {
                               stop("Could not convert variable names to numbers.")
                           }
                       }
                       )
            } else {
                labels <- as.character(labels)
            }
        }
        if (missing(xlim)) xlim <- xnum[c(1, length(xnum))] # Needed for reverted scales
        matplot(xnum, L, xlab = xlab, ylab = ylab, type = type,
                lty = lty, lwd = lwd, pch = pch, cex = cex, col = col,
                xaxt = xaxt, xlim = xlim, ...)
        if (!missing(labels) && xaxt == "n") {
            if (isTRUE(pretty.xlabels)) {
                ticks <- axTicks(1)
                ticks <- ticks[ticks >= 1 & ticks <= length(labels)]
            } else {
                ticks <- 1:length(labels)
            }
            axis(1, ticks, labels[ticks], ...)
        }
        if (!missing(legendpos)) {
            ## Are we plotting lines?
            dolines <- type %in% c("l", "b", "c", "o", "s", "S", "h")
            ## Are we plotting points?
            dopoints <- type %in% c("p", "b", "o")
            if (length(lty) > nComps) lty <- lty[1:nComps]
            do.call("legend", c(list(legendpos, varlab, col = col),
                                if (dolines) list(lty = lty, lwd = lwd),
                                if (dopoints) list(pch = pch, pt.cex = cex,
                                                   pt.lwd = lwd)))
        }
        if (isTRUE(identify))
            identify(c(row(L)), c(L),
                     labels = paste(c(col(L)), rownames(L), sep = ": "))
    }                                   # if (isTRUE(scatter))
}

## A plot method for loadings (loadings, loading.weights or Yloadings):
#' @rdname scoreplot
#' @export
plot.loadings <- function(x, ...) loadingplot(x, ...)


###
### Correlation loadings plot
###

#' @rdname scoreplot
#' @export
corrplot <- function(object, comps = 1:2, labels, plotx = TRUE, ploty = FALSE,
                     radii = c(sqrt(1/2), 1), identify = FALSE,
                     type = "p", xlab, ylab, col, ...)
{
    nComps <- length(comps)
    if (nComps < 2) stop("At least two components must be selected.")
    if (is.matrix(object)) {
        ## Assume this is already a correlation matrix
        cl <- object[,comps, drop = FALSE]
        numX <- nrow(cl)
        numY <- 0
        varlab <- colnames(cl)
    } else {
        S <- scores(object)[,comps, drop = FALSE]
        if (is.null(S))
            stop("`", deparse(substitute(object)), "' has no scores.")

        if (isTRUE(plotx)) {
            clX <- cor(model.matrix(object), S)
            numX <- nrow(clX)
        } else {
            clX <- NULL
            numX <- 0
        }
        if (isTRUE(ploty)) {
            clY <- cor(model.response(model.frame(object)), S)
            numY <- nrow(clY)
            if (numY == 1) {
                ## Add response name for single response models
                rownames(clY) <- all.vars(formula(object))[1]
            }
        } else {
            clY <- NULL
            numY <- 0
        }
        cl <- rbind(clX, clY)
        varlab <- compnames(object, comps, explvar = TRUE)
    }
    if (missing(col)) {
        ## Set up default colours:
        col <- c(rep(1, numX), rep(2, numY))
    }
    if (!missing(labels)) {
        ## Set up point labels
        if (length(labels) == 1) {
            labels <- switch(match.arg(labels, c("names", "numbers")),
                             names = rownames(cl),
                             numbers = 1:nrow(cl)
                             )
        }
        labels <- as.character(labels)
        type <- "n"
    }
    ## Build the expression to add circles:
    if (length(radii)) {
        addcircles <- substitute(symbols(cent, cent, circles = radii,
                                         inches = FALSE, add = TRUE),
                                 list(cent = rep(0, length(radii))))
    } else {
        addcircles <- expression()
    }
    if (nComps == 2) {
        ## Second component versus first
        if (missing(xlab)) xlab <- varlab[1]
        if (missing(ylab)) ylab <- varlab[2]
        plot(cl, xlim = c(-1,1), ylim = c(-1,1), asp = 1,
             xlab = xlab, ylab = ylab, type = type, col = col, ...)
        eval(addcircles)
        segments(x0 = c(-1, 0), y0 = c(0, -1), x1 = c(1, 0), y1 = c(0, 1))
        if (!missing(labels)) text(cl, labels, col = col, ...)
        if (isTRUE(identify)) {
            if (!is.null(rownames(cl))) {
                identify(cl, labels = rownames(cl))
            } else {
                identify(cl)
            }
        }
    } else {
        ## Pairwise scatterplots of several components
        pointsOrText <- if (missing(labels)) {
            function(x, y, ...) points(x, y, type = type, ...)
        } else {
            function(x, y, ...) text(x, y, labels = labels, ...)
        }
        panel <- function(x, y, ...) {
            ## Ignore the leading `ghost points':
            pointsOrText(x[-(1:2)], y[-(1:2)], ...)
            eval(addcircles)
            segments(x0 = c(-1, 0), y0 = c(0, -1), x1 = c(1, 0),
                     y1 = c(0, 1))
        }
        ## Call `pairs' with two leading `ghost points', to get
        ## correct xlim and ylim:
        pairs(rbind(-1, 1, cl), labels = varlab, panel = panel, asp = 1,
              col = col, ...)
    }
}


###
### prediction plot
###

## Generic:


#' @name predplot
#' @title Prediction Plots
#'
#' @description Functions to plot predicted values against measured values for a fitted
#' model.
#'
#' @details \code{predplot} is a generic function for plotting predicted versus measured
#' response values, with default and \code{mvr} methods currently implemented.
#' The default method is very simple, and doesn't handle multiple responses or
#' new data.
#'
#' The \code{mvr} method, handles multiple responses, model sizes and types of
#' predictions by making one plot for each combination.  It can also be called
#' through the plot method for \code{mvr}, by specifying \code{plottype =
#' "prediction"} (the default).
#'
#' The argument \code{main} can be used to specify the main title of the plot.
#' It is handled in a non-standard way.  If there is only on (sub) plot,
#' \code{main} will be used as the main title of the plot.  If there is
#' \emph{more} than one (sub) plot, however, the presence of \code{main} will
#' produce a corresponding \sQuote{global} title on the page.  Any graphical
#' parametres, e.g., \code{cex.main}, supplied to \code{coefplot} will only
#' affect the \sQuote{ordinary} plot titles, not the \sQuote{global} one.  Its
#' appearance can be changed by setting the parameters with \code{\link{par}},
#' which will affect \emph{both} titles (with the exception of \code{font.main}
#' and \code{cex.main}, which will only affect the \sQuote{global} title when
#' there is more than one plot).  (To have different settings for the two
#' titles, one can override the \code{par} settings with arguments to
#' \code{predplot}.)
#'
#' \code{predplotXy} is an internal function and is not meant for interactive
#' use.  It is called by the \code{predplot} methods, and its arguments, e.g,
#' \code{line}, can be given in the \code{predplot} call.
#'
#' @aliases predplot predplot.default predplot.mvr predplotXy
#' @param object a fitted model.
#' @param ncomp integer vector.  The model sizes (numbers of components) to use
#' for prediction.
#' @param which character vector.  Which types of predictions to plot.  Should
#' be a subset of \code{c("train", "validation", "test")}.  If not specified,
#' \code{plot.mvr} selects test set predictions if \code{newdata} is supplied,
#' otherwise cross-validated predictions if the model has been cross-validated,
#' otherwise fitted values from the calibration data.
#' @param newdata data frame.  New data to predict.
#' @param nCols,nRows integer.  The number of coloumns and rows the plots will
#' be laid out in.  If not specified, \code{plot.mvr} tries to be intelligent.
#' @param xlab,ylab titles for \eqn{x} and \eqn{y} axes.  Typically character
#' strings, but can be expressions or lists.  See \code{\link{title}} for
#' details.
#' @param labels optional.  Alternative plot labels to use.  Either a vector of
#' labels, or \code{"names"} or \code{"numbers"} to use the row names or row
#' numbers of the data as labels.
#' @param type character.  What type of plot to make.  Defaults to \code{"p"}
#' (points).  See \code{\link{plot}} for a complete list of types.  The
#' argument is ignored if \code{labels} is specified.
#' @param main optional main title for the plot.  See Details.
#' @param ask logical.  Whether to ask the user before each page of a plot.
#' @param font.main font to use for main titles.  See \code{\link{par}} for
#' details.  Also see Details below.
#' @param cex.main numeric.  The magnification to be used for main titles
#' relative to the current size.  Also see Details below.
#' @param x numeric vector.  The observed response values.
#' @param y numeric vector.  The predicted response values.
#' @param line logical.  Whether a target line should be drawn.
#' @param line.col,line.lty,line.lwd character or numeric.  The \code{col},
#' \code{lty} and \code{lwd} parametres for the target line.  See
#' \code{\link{par}} for details.
#' @param \dots further arguments sent to underlying plot functions.
#' @return The functions invisibly return a matrix with the (last) plotted
#' data.
#' @note The \code{font.main} and \code{cex.main} must be (completely) named.
#' This is to avoid that any argument \code{cex} or \code{font} matches them.
#'
#' Tip: If the labels specified with \code{labels} are too long, they get
#' clipped at the border of the plot region.  This can be avoided by supplying
#' the graphical parameter \code{xpd = TRUE} in the plot call.
#' @author Ron Wehrens and Bjørn-Helge Mevik
#' @seealso \code{\link{mvr}}, \code{\link{plot.mvr}}
#' @keywords regression multivariate hplot
#' @examples
#'
#' data(yarn)
#' mod <- plsr(density ~ NIR, ncomp = 10, data = yarn[yarn$train,], validation = "CV")
#' \dontrun{
#' predplot(mod, ncomp = 1:6)
#' plot(mod, ncomp = 1:6) # Equivalent to the previous
#' ## Both cross-validated and test set predictions:
#' predplot(mod, ncomp = 4:6, which = c("validation", "test"),
#'          newdata = yarn[!yarn$train,])
#' }
#'
#' data(oliveoil)
#' mod.sens <- plsr(sensory ~ chemical, ncomp = 4, data = oliveoil)
#' \dontrun{plot(mod.sens, ncomp = 2:4) # Several responses gives several plots}
#'
#' @export
predplot <- function(object, ...)
  UseMethod("predplot")

## Default method:
#' @rdname predplot
#' @export
predplot.default <- function(object, ...) {
    measured <- model.response(model.frame(object))
    predicted <- predict(object)
    predplotXy(measured, predicted, ...)
}

## Method for mvr objects:
#' @rdname predplot
#' @export
predplot.mvr <- function(object, ncomp = object$ncomp, which, newdata,
                         nCols, nRows, xlab = "measured", ylab = "predicted",
                         main,
                         ask = nRows * nCols < nPlots && dev.interactive(),
                         ..., font.main, cex.main)
{
    ## Select type(s) of prediction
    if (missing(which)) {
        ## Pick the `best' alternative.
        if (!missing(newdata)) {
            which <- "test"
        } else {
            if (!is.null(object$validation)) {
                which <- "validation"
            } else {
                which <- "train"
            }
        }
    } else {
        ## Check the supplied `which'
        allTypes <- c("train", "validation", "test")
        which <- allTypes[pmatch(which, allTypes)]
        if (length(which) == 0 || any(is.na(which)))
            stop("`which' should be a subset of ",
                 paste(allTypes, collapse = ", "))
    }

    ## Help variables
    nEst  <- length(which)
    nSize <- length(ncomp)
    nResp <- dim(object$fitted.values)[2]

    ## Set plot parametres as needed:
    dims <- c(nEst, nSize, nResp)
    dims <- dims[dims > 1]
    nPlots <- prod(dims)
    if (nPlots > 1) {
        ## Set up default font.main and cex.main for individual titles:
        if (missing(font.main)) font.main <- 1
        if (missing(cex.main)) cex.main <- 1.1
        ## Show the *labs in the margin:
        mXlab <- xlab
        mYlab <- ylab
        xlab <- ylab <- ""
        if(missing(nCols)) nCols <- min(c(3, dims[1]))
        if(missing(nRows))
            nRows <- min(c(3, ceiling(prod(dims[1:2], na.rm = TRUE) / nCols)))
        opar <- par(no.readonly = TRUE)
        on.exit(par(opar))
        par(mfrow = c(nRows, nCols),
            oma = c(1, 1, if(missing(main)) 0 else 2, 0) + 0.1,
            mar = c(3,3,3,1) + 0.1)
        if (isTRUE(ask)) {
            oask <- devAskNewPage(TRUE)
            on.exit(devAskNewPage(oask))
        }
    } else {
        ## Set up default font.main and cex.main for the main title:
        if (missing(font.main)) font.main <- par("font.main")
        if (missing(cex.main)) cex.main <- par("cex.main")
        nCols <- nRows <- 1
    }

    ## Set up measured and predicted for all model sizes, responses and
    ## estimates:
    if ("train" %in% which) {
        train.measured <- as.matrix(model.response(model.frame(object)))
        train.predicted <- object$fitted.values[,,ncomp, drop = FALSE]
    }
    if ("validation" %in% which) {
        if (is.null(object$validation)) stop("`object' has no `validation' component.")
        if(!exists("train.measured"))
            train.measured <- as.matrix(model.response(model.frame(object)))
        validation.predicted <- object$validation$pred[,,ncomp, drop = FALSE]
    }
    if ("test" %in% which) {
        if (missing(newdata)) stop("Missing `newdata'.")
        test.measured <- as.matrix(model.response(model.frame(formula(object),
                                                              data = newdata)))
        test.predicted <- predict(object, ncomp = ncomp, newdata = newdata)
    }

    ## Do the plots
    plotNo <- 0
    for (resp in 1:nResp) {
        for (size in 1:nSize) {
            for (est in 1:nEst) {
                plotNo <- plotNo + 1
                if (nPlots == 1 && !missing(main)) {
                    lmain <- main
                } else {
                    lmain <- sprintf("%s, %d comps, %s",
                                     respnames(object)[resp],
                                     ncomp[size], which[est])
                }
                sub <- which[est]
                switch(which[est],
                       train = {
                           measured <- train.measured[,resp]
                           predicted <- train.predicted[,resp,size]
                       },
                       validation = {
                           measured <- train.measured[,resp]
                           predicted <- validation.predicted[,resp,size]
                       },
                       test = {
                           measured <- test.measured[,resp]
                           predicted <- test.predicted[,resp,size]
                       }
                       )
                xy <- predplotXy(measured, predicted, main = lmain,
                                 font.main = font.main, cex.main = cex.main,
                                 xlab = xlab, ylab = ylab, ...)
                if (nPlots > 1 &&
                    (plotNo %% (nCols * nRows) == 0 || plotNo == nPlots)) {
                    ## Last plot on a page; add outer margin text and title:
                    mtext(mXlab, side = 1, outer = TRUE)
                    mtext(mYlab, side = 2, outer = TRUE)
                    if (!missing(main)) title(main = main, outer = TRUE)
                }
            }
        }
    }
    invisible(xy)
}

## The workhorse function:
#' @rdname predplot
#' @export
predplotXy <- function(x, y, line = FALSE, labels, type = "p",
                       main = "Prediction plot", xlab = "measured response",
                       ylab = "predicted response", line.col = par("col"),
                       line.lty = NULL, line.lwd = NULL, ...)
{
    if (!missing(labels)) {
        ## Set up point labels
        if (length(labels) == 1) {
            labels <- switch(match.arg(labels, c("names", "numbers")),
                             names = names(y),
                             numbers = as.character(1:length(y))
                             )
        }
        ## Override plot type:
        type <- "n"
    }
    plot(y ~ x, type = type, main = main, xlab = xlab, ylab = ylab, ...)
    if (!missing(labels)) text(x, y, labels, ...)
    if (isTRUE(line)) abline(0, 1, col = line.col, lty = line.lty, lwd = line.lwd)
    invisible(cbind(measured = x, predicted = as.vector(y)))
}


###
### Coefficient plot
###



#' @title Plot Regression Coefficients of PLSR and PCR models
#'
#' @description  Function to plot the regression coefficients of an \code{mvr} object.
#'
#' @details \code{coefplot} handles multiple responses by making one plot for each
#' response.  If \code{separate} is \code{TRUE}, separate plots are made for
#' each combination of model size and response.  The plots are laid out in a
#' rectangular fashion.
#'
#' If \code{legendpos} is given, a legend is drawn at the given position
#' (unless \code{separate} is \code{TRUE}).
#'
#' The argument \code{labels} can be a vector of labels or one of
#' \code{"names"} and \code{"numbers"}.  The labels are used as \eqn{x} axis
#' labels.  If \code{labels} is \code{"names"} or \code{"numbers"}, the
#' variable names are used as labels, the difference being that with
#' \code{"numbers"}, the variable names are converted to numbers, if possible.
#' Variable names of the forms \samp{"number"} or \samp{"number text"} (where
#' the space is optional), are handled.
#'
#' The argument \code{main} can be used to specify the main title of the plot.
#' It is handled in a non-standard way.  If there is only on (sub) plot,
#' \code{main} will be used as the main title of the plot.  If there is
#' \emph{more} than one (sub) plot, however, the presence of \code{main} will
#' produce a corresponding \sQuote{global} title on the page.  Any graphical
#' parametres, e.g., \code{cex.main}, supplied to \code{coefplot} will only
#' affect the \sQuote{ordinary} plot titles, not the \sQuote{global} one.  Its
#' appearance can be changed by setting the parameters with \code{\link{par}},
#' which will affect \emph{both} titles.  (To have different settings for the
#' two titles, one can override the \code{par} settings with arguments to
#' \code{coefplot}.)
#'
#' The argument \code{pretty.xlabels} is only used when \code{labels} is
#' specified.  If \code{TRUE} (default), the code tries to use a
#' \sQuote{pretty} selection of labels.  If \code{labels} is \code{"numbers"},
#' it also uses the numerical values of the labels for horisontal spacing.  If
#' one has excluded parts of the spectral region, one might therefore want to
#' use \code{pretty.xlabels = FALSE}.
#'
#' When \code{separate} is \code{TRUE}, the arguments \code{lty}, \code{col},
#' and \code{pch} default to their \code{par()} setting.  Otherwise, the
#' default for all of them is \code{1:nLines}, where \code{nLines} is the
#' number of model sizes specified, i.e., the length of \code{ncomp} or
#' \code{comps}.
#'
#' The function can also be called through the \code{mvr} plot method by
#' specifying \code{plottype = "coefficients"}.
#'
#' @param object an \code{mvr} object.  The fitted model.
#' @param ncomp,comps vector of positive integers.  The components to plot.
#' See \code{\link{coef.mvr}} for details.
#' @param separate logical.  If \code{TRUE}, coefficients for different model
#' sizes are blotted in separate plots.
#' @param se.whiskers logical.  If \code{TRUE}, whiskers at plus/minus 1
#' estimated standard error are added to the plot.  This is only available if
#' the model was cross-validated with \code{jackknife = TRUE}.  Also, in the
#' current implementation, \code{intercept} must be \code{FALSE}, and
#' \code{separate} must be \code{TRUE} if \code{length(ncomp) > 1}.
#' @param intercept logical.  Whether coefficients for the intercept should be
#' plotted.  Ignored if \code{comps} is specified.  Defaults to \code{FALSE}.
#' See \code{\link{coef.mvr}} for details.
#' @param nCols,nRows integer.  The number of coloumns and rows the plots will
#' be laid out in.  If not specified, \code{coefplot} tries to be intelligent.
#' @param labels optional.  Alternative \eqn{x} axis labels.  See Details.
#' @param type character.  What type of plot to make.  Defaults to \code{"l"}
#' (lines).  Alternative types include \code{"p"} (points) and \code{"b"}
#' (both).  See \code{\link{plot}} for a complete list of types.
#' @param lty vector of line types (recycled as neccessary).  Line types can be
#' specified as integers or character strings (see \code{\link{par}} for the
#' details).
#' @param lwd vector of positive numbers (recycled as neccessary), giving the
#' width of the lines.
#' @param pch plot character.  A character string or a vector of single
#' characters or integers (recycled as neccessary).  See \code{\link{points}}
#' for all alternatives.
#' @param cex numeric vector of character expansion sizes (recycled as
#' neccessary) for the plotted symbols.
#' @param col character or integer vector of colors for plotted lines and
#' symbols (recycled as neccessary).  See \code{\link{par}} for the details.
#' @param legendpos Legend position.  Optional.  Ignored if \code{separate} is
#' \code{TRUE}.  If present, a legend is drawn at the given position.  The
#' position can be specified symbolically (e.g., \code{legendpos =
#' "topright"}).  This requires >= 2.1.0.  Alternatively, the position can be
#' specified explicitly (\code{legendpos = t(c(x,y))}) or interactively
#' (\code{legendpos = \link{locator}()}).  This only works well for plots of
#' single-response models.
#' @param xlab,ylab titles for \eqn{x} and \eqn{y} axes.  Typically character
#' strings, but can be expressions (e.g., \code{expression(R^2)} or lists.  See
#' \code{\link{title}} for details.
#' @param main optional main title for the plot.  See Details.
#' @param pretty.xlabels logical.  If \code{TRUE}, \code{coefplot} tries to
#' plot the \eqn{x} labels more nicely.  See Details.
#' @param xlim,ylim optional vector of length two, with the \eqn{x} or \eqn{y}
#' limits of the plot.
#' @param ask logical.  Whether to ask the user before each page of a plot.
#' @param \dots Further arguments sent to the underlying plot functions.
#' @note \code{\link{legend}} has many options.  If you want greater control
#' over the appearance of the legend, omit the \code{legendpos} argument and
#' call \code{legend} manually.
#'
#' The handling of \code{labels} and \code{pretty.xlabels} is experimental.
#' @author Ron Wehrens and Bjørn-Helge Mevik
#' @seealso \code{\link{mvr}}, \code{\link{plot.mvr}}, \code{\link{coef.mvr}},
#' \code{\link{plot}}, \code{\link{legend}}
#' @keywords regression multivariate hplot
#' @examples
#'
#' data(yarn)
#' mod.nir <- plsr(density ~ NIR, ncomp = 8, data = yarn)
#' \dontrun{
#' coefplot(mod.nir, ncomp = 1:6)
#' plot(mod.nir, plottype = "coefficients", ncomp = 1:6) # Equivalent to the previous
#' ## Plot with legend:
#' coefplot(mod.nir, ncom = 1:6, legendpos = "bottomright")
#' }
#'
#' data(oliveoil)
#' mod.sens <- plsr(sensory ~ chemical, ncomp = 4, data = oliveoil)
#' \dontrun{coefplot(mod.sens, ncomp = 2:4, separate = TRUE)}
#'
#' @export
coefplot <- function(object, ncomp = object$ncomp, comps, intercept = FALSE,
                     separate = FALSE, se.whiskers = FALSE,
                     nCols, nRows, labels,
                     type = "l", lty, lwd = NULL,
                     pch, cex = NULL, col, legendpos,
                     xlab = "variable", ylab = "regression coefficient",
                     main, pretty.xlabels = TRUE, xlim, ylim,
                     ask = nRows * nCols < nPlots && dev.interactive(), ...)
{
    ## This simplifies code below:
    if (missing(comps)) comps <- NULL
    separate <- isTRUE(separate)
    se.whiskers <- isTRUE(se.whiskers)

    ## Help variables
    nLines <- if (is.null(comps)) length(ncomp) else length(comps)
    nSize <- if (separate) nLines else 1
    nResp <- dim(object$fitted.values)[2]

    ## Set plot parametres as needed:
    dims <- c(nSize, nResp)
    dims <- dims[dims > 1]
    nPlots <- prod(dims)
    if (nPlots > 1) {
        ## Show the *labs in the margin:
        mXlab <- xlab
        mYlab <- ylab
        xlab <- ylab <- ""
        if (missing(nCols)) nCols <- min(c(3, dims[1]))
        if (missing(nRows))
            nRows <- min(c(3, ceiling(prod(dims[1:2], na.rm = TRUE) / nCols)))
        opar <- par(no.readonly = TRUE)
        on.exit(par(opar))
        par(mfrow = c(nRows, nCols),
            oma = c(1, 1, if(missing(main)) 0 else 2, 0) + 0.1,
            mar = c(3,3,3,1) + 0.1)
        if (isTRUE(ask)) {
            oask <- devAskNewPage(TRUE)
            on.exit(devAskNewPage(oask))
        }
    } else {
        nCols <- nRows <- 1
    }
    if (missing(col)) col <- if (separate) par("col") else 1:nLines
    if (missing(pch)) pch <- if (separate) par("pch") else 1:nLines
    if (missing(lty)) lty <- if (separate) par("lty") else 1:nLines
    if (length(lty) > nLines) lty <- lty[1:nLines] # otherwise legend chokes
    if (length(type) != 1 || is.na(nchar(type, "c")) || nchar(type, "c") != 1)
        stop("Invalid plot type.")
    ## Are we plotting lines?
    dolines <- type %in% c("l", "b", "c", "o", "s", "S", "h")
    ## Are we plotting points?
    dopoints <- type %in% c("p", "b", "o")

    ## Get the coefficients:
    coefs <- coef(object, ncomp = ncomp, comps = comps, intercept = intercept)
    complabs <- dimnames(coefs)[[3]]

    ## Optionally, get the standard errors:
    if (se.whiskers) {
        if (isTRUE(intercept))
            stop(sQuote("se.whiskers"), " not supported when ",
                 sQuote("intercept"), " is TRUE")
        if (!is.null(comps))
            stop(sQuote("se.whiskers"), " not supported when ",
                 sQuote("comps"), " is specified")
        if (dim(coefs)[3] > 1 && !separate)
            stop(sQuote("se.whiskers"), " not supported when ",
                 sQuote("separate"), " is FALSE and length(ncomp) > 1")
        SEs <- sqrt(var.jack(object, ncomp = ncomp))
        npred <- dim(SEs)[1]
        if (!hasArg("ylim")) {
            ## Calculate new ylims:
            miny <- apply(coefs - SEs, 2:3, min)
            maxy <- apply(coefs + SEs, 2:3, max)
        }
    }

    ## Set up the x labels:
    xnum <- 1:dim(coefs)[1]
    if (missing(labels)) {
        xaxt <- par("xaxt")
    } else {
        xaxt <- "n"
        if (length(labels) == 1) {
            xnam <- prednames(object, intercept = intercept)
            switch(match.arg(labels, c("names", "numbers")),
                   names = {            # Simply use the names as is
                       labels <- xnam
                   },
                   numbers = {          # Try to use them as numbers
                       if (length(grep("^[-0-9.]+[^0-9]*$", xnam)) ==
                           length(xnam)) {
                           ## Labels are on "num+text" format
                           labels <- sub("[^0-9]*$", "", xnam)
                           if (isTRUE(pretty.xlabels)) {
                               xnum <- as.numeric(labels)
                               xaxt <- par("xaxt")
                           }
                       } else {
                           stop("Could not convert variable names to numbers.")
                       }
                   }
                   )
        } else {
            labels <- as.character(labels)
        }
    }
    if (missing(xlim)) xlim <- xnum[c(1, length(xnum))] # Needed for reverted scales
    ## Do the plots
    plotNo <- 0
    for (resp in 1:nResp) {
        respname <- respnames(object)[resp]
        for (size in 1:nSize) {
            plotNo <- plotNo + 1

            if (nPlots == 1 && !missing(main)) {
                lmain <- main
            } else if (separate) {
                lmain <- paste(respname, complabs[size], sep = ", ")
            } else {
                lmain <- respname
            }
            if (separate) {
                if (missing(ylim)) {
                    if (se.whiskers) {
                        ylims <- c(miny[resp,size], maxy[resp,size])
                    } else {
                        ylims <- range(coefs[,resp,size])
                    }
                } else {
                    ylims <- ylim
                }
                plot(xnum, coefs[,resp,size],
                     main = lmain, xlab = xlab, ylab = ylab, type = type,
                     lty = lty, lwd = lwd, pch = pch, cex = cex,
                     col = col, xaxt = xaxt, xlim = xlim, ylim = ylims, ...)
                if (se.whiskers) {
                    arrows(1:npred, (coefs - SEs)[,resp,size],
                           1:npred, (coefs + SEs)[,resp,size], length = 0.05,
                           angle = 90, code = 3, col = 2)
                }
            } else {
                if (missing(ylim)) {
                    if (se.whiskers) {
                        ylims <- c(miny[resp,], maxy[resp,])
                    } else {
                        ylims <- range(coefs[,resp,])
                    }
                } else {
                    ylims <- ylim
                }
                matplot(xnum, coefs[,resp,], main = lmain, xlab = xlab,
                        ylab = ylab, type = type, lty = lty, lwd = lwd,
                        pch = pch, cex = cex, col = col, xaxt = xaxt,
                        xlim = xlim, ylim = ylims, ...)
                if (se.whiskers) {
                    arrows(1:npred, (coefs - SEs)[,resp,],
                           1:npred, (coefs + SEs)[,resp,], length = 0.05,
                           angle = 90, code = 3, col = 2)
                }
                if (!missing(legendpos)) {
                    do.call("legend", c(list(legendpos, complabs, col = col),
                                        if(dolines) list(lty = lty, lwd = lwd),
                                        if(dopoints) list(pch = pch,
                                                          pt.cex = cex,
                                                          pt.lwd = lwd)))
                }
            }
            if (!missing(labels) && xaxt == "n") {
                if (isTRUE(pretty.xlabels)) {
                    ticks <- axTicks(1)
                    ticks <- ticks[ticks >= 1 & ticks <= length(labels)]
                } else {
                    ticks <- 1:length(labels)
                }
                axis(1, ticks, labels[ticks], ...)
            }
            abline(h = 0, col = "gray")

            if (nPlots > 1 &&
                (plotNo %% (nCols * nRows) == 0 || plotNo == nPlots)) {
                ## Last plot on a page; add outer margin text and title:
                mtext(mXlab, side = 1, outer = TRUE)
                mtext(mYlab, side = 2, outer = TRUE)
                if (!missing(main)) title(main, outer = TRUE)
            }
        }
    }
}


###
### Validation plot (MSEP/RMSEP/R2)
###



#' @name validationplot
#' @title Validation Plots
#'
#' @description Functions to plot validation statistics, such as RMSEP or \eqn{R^2}, as a
#' function of the number of components.
#'
#' @details \code{validationplot} calls the proper validation function (currently
#' \code{\link{MSEP}}, \code{\link{RMSEP}} or \code{\link{R2}}) and plots the
#' results with \code{plot.mvrVal}.  \code{validationplot} can be called
#' through the \code{mvr} plot method, by specifying \code{plottype =
#' "validation"}.
#'
#' \code{plot.mvrVal} creates one plot for each response variable in the model,
#' laid out in a rectangle.  It uses \code{\link{matplot}} for performing the
#' actual plotting.  If \code{legendpos} is given, a legend is drawn at the
#' given position.
#'
#' The argument \code{main} can be used to specify the main title of the plot.
#' It is handled in a non-standard way.  If there is only on (sub) plot,
#' \code{main} will be used as the main title of the plot.  If there is
#' \emph{more} than one (sub) plot, however, the presence of \code{main} will
#' produce a corresponding \sQuote{global} title on the page.  Any graphical
#' parametres, e.g., \code{cex.main}, supplied to \code{coefplot} will only
#' affect the \sQuote{ordinary} plot titles, not the \sQuote{global} one.  Its
#' appearance can be changed by setting the parameters with \code{\link{par}},
#' which will affect \emph{both} titles.  (To have different settings for the
#' two titles, one can override the \code{par} settings with arguments to the
#' plot function.)
#'
#' @aliases validationplot plot.mvrVal
#' @param object an \code{mvr} object.
#' @param val.type character.  What type of validation statistic to plot.
#' @param estimate character.  Which estimates of the statistic to calculate.
#' See \code{\link{RMSEP}}.
#' @param newdata data frame.  Optional new data used to calculate statistic.
#' @param ncomp,comps integer vector.  The model sizes to compute the statistic
#' for.  See \code{\link{RMSEP}}.
#' @param intercept logical.  Whether estimates for a model with zero
#' components should be calculated as well.
#' @param x an \code{mvrVal} object.  Usually the result of a
#' \code{\link{RMSEP}}, \code{\link{MSEP}} or \code{\link{R2}} call.
#' @param nCols,nRows integers.  The number of coloumns and rows the plots will
#' be laid out in.  If not specified, \code{plot.mvrVal} tries to be
#' intelligent.
#' @param type character.  What type of plots to create.  Defaults to
#' \code{"l"} (lines).  Alternative types include \code{"p"} (points) and
#' \code{"b"} (both).  See \code{\link{plot}} for a complete list of types.
#' @param lty vector of line types (recycled as neccessary).  Line types can be
#' specified as integers or character strings (see \code{\link{par}} for the
#' details).
#' @param lwd vector of positive numbers (recycled as neccessary), giving the
#' width of the lines.
#' @param pch plot character.  A character string or a vector of single
#' characters or integers (recycled as neccessary).  See \code{\link{points}}
#' for all alternatives.
#' @param cex numeric vector of character expansion sizes (recycled as
#' neccessary) for the plotted symbols.
#' @param col character or integer vector of colors for plotted lines and
#' symbols (recycled as neccessary).  See \code{\link{par}} for the details.
#' @param legendpos Legend position.  Optional.  If present, a legend is drawn
#' at the given position.  The position can be specified symbolically (e.g.,
#' \code{legendpos = "topright"}).  This requires >= 2.1.0.  Alternatively, the
#' position can be specified explicitly (\code{legendpos = t(c(x,y))}) or
#' interactively (\code{legendpos = \link{locator}()}).  This only works well
#' for plots of single-response models.
#' @param xlab,ylab titles for \eqn{x} and \eqn{y} axes.  Typically character
#' strings, but can be expressions (e.g., \code{expression(R^2)} or lists.  See
#' \code{\link{title}} for details.
#' @param main optional main title for the plot.  See Details.
#' @param ask logical.  Whether to ask the user before each page of a plot.
#' @param \dots Further arguments sent to underlying plot functions.
#' @note \code{\link{legend}} has many options.  If you want greater control
#' over the appearance of the legend, omit the \code{legendpos} argument and
#' call \code{legend} manually.
#' @author Ron Wehrens and Bjørn-Helge Mevik
#' @seealso \code{\link{mvr}}, \code{\link{plot.mvr}}, \code{\link{RMSEP}},
#' \code{\link{MSEP}}, \code{\link{R2}}, \code{\link{matplot}},
#' \code{\link{legend}}
#' @keywords regression multivariate hplot
#' @examples
#'
#' data(oliveoil)
#' mod <- plsr(sensory ~ chemical, data = oliveoil, validation = "LOO")
#' \dontrun{
#' ## These three are equivalent:
#' validationplot(mod, estimate = "all")
#' plot(mod, "validation", estimate = "all")
#' plot(RMSEP(mod, estimate = "all"))
#' ## Plot R2:
#' plot(mod, "validation", val.type = "R2")
#' ## Plot R2, with a legend:
#' plot(mod, "validation", val.type = "MSEP", legendpos = "top") # R >= 2.1.0
#' }
#'
#' @export
validationplot <- function(object, val.type = c("RMSEP", "MSEP", "R2"),
                           estimate, newdata, ncomp, comps, intercept, ...)
{
    cl <- match.call(expand.dots = FALSE)
    cl[[1]] <- as.name(match.arg(val.type))
    cl$val.type <- NULL
    x <- eval(cl, parent.frame())
    plot(x, ...)
}

## A plot method for mvrVal objects:
#' @rdname validationplot
#' @export
plot.mvrVal <- function(x, nCols, nRows, type = "l", lty = 1:nEst,
                        lwd = par("lwd"), pch = 1:nEst, cex = 1, col = 1:nEst,
                        legendpos, xlab = "number of components",
                        ylab = x$type, main,
                        ask = nRows * nCols < nResp && dev.interactive(), ...)
{
    if (!is.null(x$call$cumulative) && eval(x$call$cumulative) == FALSE)
        stop("`cumulative = FALSE' not supported.")
    ## Set plot parametres as needed:
    nResp <- dim(x$val)[2]              # Number of response variables
    if (nResp > 1) {
        ## Show the *labs in the margin:
        mXlab <- xlab
        mYlab <- ylab
        xlab <- ylab <- ""
        if(missing(nCols)) nCols <- min(c(3, nResp))
        if(missing(nRows)) nRows <- min(c(3, ceiling(nResp / nCols)))
        opar <- par(no.readonly = TRUE)
        on.exit(par(opar))
        par(mfrow = c(nRows, nCols),
            oma = c(1, 1, if(missing(main)) 0 else 2, 0) + 0.1,
            mar = c(3,3,3,1) + 0.1)
        if (isTRUE(ask)) {
            oask <- devAskNewPage(TRUE)
            on.exit(devAskNewPage(oask))
        }
    } else {
        nCols <- nRows <- 1
    }
    ynames <- dimnames(x$val)[[2]]      # Names of response variables
    estnames <- dimnames(x$val)[[1]]    # Names of estimators
    nEst <- length(estnames)
    if (length(lty) > nEst) lty <- lty[1:nEst] # otherwise legend chokes
    if (length(type) != 1 || is.na(nchar(type, "c")) || nchar(type, "c") != 1)
        stop("Invalid plot type.")
    ## Are we plotting lines?
    dolines <- type %in% c("l", "b", "c", "o", "s", "S", "h")
    ## Are we plotting points?
    dopoints <- type %in% c("p", "b", "o")

    for (resp in 1:nResp) {
        if (nResp == 1 && !missing(main)) {
            lmain <- main
        } else {
            lmain <- ynames[resp]
        }
        y <- x$val[,resp,]
        if (is.matrix(y)) y <- t(y)
        if (isTRUE(all.equal(x$comps, min(x$comps):max(x$comps)))) {
            matplot(x$comps, y, xlab = xlab, ylab = ylab, main = lmain,
                    type = type, lty = lty, lwd = lwd, pch = pch, cex = cex,
                    col = col, ...)
        } else {
            ## Handle irregular x$comps:
            matplot(y, xlab = xlab, ylab = ylab, main = lmain,
                    xaxt = "n", type = type, lty = lty, lwd = lwd,
                    pch = pch, cex = cex, col = col, ...)
            axis(1, seq(along = x$comps), x$comps)
        }
        if (!missing(legendpos)) {
            do.call("legend", c(list(legendpos, estnames, col = col),
                                if (dolines) list(lty = lty, lwd = lwd),
                                if (dopoints) list(pch = pch, pt.cex = cex,
                                                   pt.lwd = lwd)))
        }
        if (nResp > 1 && (resp %% (nCols * nRows) == 0 || resp == nResp)) {
            ## Last plot on a page; add outer margin text and title:
            mtext(mXlab, side = 1, outer = TRUE)
            mtext(mYlab, side = 2, outer = TRUE)
            if (!missing(main)) title(main, outer = TRUE)
        }
    }
}


###
### biplot
###



#' @name biplot.mvr
#' @title Biplots of PLSR and PCR Models.
#'
#' @description Biplot method for \code{mvr} objects.
#'
#' @details \code{biplot.mvr} can also be called through the \code{mvr} plot method by
#' specifying \code{plottype = "biplot"}.
#'
#' @param x an \code{mvr} object.
#' @param comps integer vector of length two.  The components to plot.
#' @param which character.  Which matrices to plot.  One of \code{"x"} (X
#' scores and loadings), \code{"y"} (Y scores and loadings), \code{"scores"} (X
#' and Y scores) and \code{"loadings"} (X and Y loadings).
#' @param var.axes logical.  If \code{TRUE}, the second set of points have
#' arrows representing them.
#' @param xlabs either a character vector of labels for the first set of
#' points, or \code{FALSE} for no labels.  If missing, the row names of the
#' first matrix is used as labels.
#' @param ylabs either a character vector of labels for the second set of
#' points, or \code{FALSE} for no labels.  If missing, the row names of the
#' second matrix is used as labels.
#' @param main character.  Title of plot.  If missing, a title is constructed
#' by \code{biplot.mvr}.
#' @param \dots Further arguments passed on to \code{biplot.default}.
#' @author Ron Wehrens and Bjørn-Helge Mevik
#' @seealso \code{\link{mvr}}, \code{\link{plot.mvr}},
#' \code{\link{biplot.default}}
#' @keywords regression multivariate hplot
#' @examples
#'
#' data(oliveoil)
#' mod <- plsr(sensory ~ chemical, data = oliveoil)
#' \dontrun{
#' ## These are equivalent
#' biplot(mod)
#' plot(mod, plottype = "biplot")
#'
#' ## The four combinations of x and y points:
#' par(mfrow = c(2,2))
#' biplot(mod, which = "x") # Default
#' biplot(mod, which = "y")
#' biplot(mod, which = "scores")
#' biplot(mod, which = "loadings")
#' }
#'
#' @export
biplot.mvr <- function(x, comps = 1:2,
                       which = c("x", "y", "scores", "loadings"),
                       var.axes = FALSE, xlabs, ylabs, main, ...)
{
    if (length(comps) != 2) stop("Exactly 2 components must be selected.")
    which <- match.arg(which)
    switch(which,
           x = {
               objects <- x$scores
               vars <- x$loadings
               title <- "X scores and X loadings"
           },
           y = {
               objects <- x$Yscores
               vars <- x$Yloadings
               title <- "Y scores and Y loadings"
           },
           scores = {
               objects <- x$scores
               vars <- x$Yscores
               title <- "X scores and Y scores"
           },
           loadings = {
               objects <- x$loadings
               vars <- x$Yloadings
               title <- "X loadings and Y loadings"
           }
           )
    if (is.null(objects) || is.null(vars))
        stop("'x' lacks the required scores/loadings.")
    ## Build a call to `biplot'
    mc <- match.call()
    mc$comps <- mc$which <- NULL
    mc$x <- objects[,comps, drop = FALSE]
    mc$y <- vars[,comps, drop = FALSE]
    if (missing(main)) mc$main <- title
    if (missing(var.axes)) mc$var.axes = FALSE
    if (!missing(xlabs) && isFALSE(xlabs))
        mc$xlabs <- rep("o", nrow(objects))
    if (!missing(ylabs) && isFALSE(ylabs))
        mc$ylabs <- rep("o", nrow(vars))
    mc[[1]] <- as.name("biplot")
    ## Evaluate the call:
    eval(mc, parent.frame())
}