plots.r

做主成分回归和偏最小二乘回归

### Plots for mvr objects.  Some of them also work for other
### objects, but that is not a priority.
###
### $Id: plots.R 148 2007-10-16 20:33:38Z bhm $

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

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
###

scoreplot <- function(object, ...) UseMethod("scoreplot")

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 (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:
plot.scores <- function(x, ...) scoreplot(x, ...)


###
### Loadingplot
###

loadingplot <- function(object, ...) UseMethod("loadingplot")

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 (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 (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 (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 (identify)
            identify(c(row(L)), c(L),
                     labels = paste(c(col(L)), rownames(L), sep = ": "))
    }                                   # if (scatter)
}

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


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

corrplot <- function(object, comps = 1:2, labels, radii = c(sqrt(1/2), 1),
                     identify = FALSE, type = "p", xlab, ylab, ...) {
    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]
        varlab <- colnames(cl)
    } else {
        S <- scores(object)[,comps, drop = FALSE]
        if (is.null(S))
            stop("`", deparse(substitute(object)), "' has no scores.")
        cl <- cor(model.matrix(object), S)
        varlab <- compnames(object, comps, explvar = TRUE)
    }
    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, ...)
        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, ...)
        if (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, ...)
    }
}


###
### prediction plot
###

## Generic:
predplot <- function(object, ...)
  UseMethod("predplot")

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

## Method for mvr objects:
predplot.mvr <- function(object, ncomp = object$ncomp, which, newdata,
                         nCols, nRows, xlab = "measured", ylab = "predicted",
                         main, ..., 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