plots.r

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

        ## 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 (nRows * nCols < nPlots && dev.interactive()) par(ask = TRUE)
    } 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:
predplotXy <- function(x, y, line = FALSE, main = "Prediction plot",
                       xlab = "measured response",
                       ylab = "predicted response", line.col = par("col"),
                       line.lty = NULL, line.lwd = NULL, ...)
{
    plot(y ~ x, main = main, xlab = xlab, ylab = ylab, ...)
    if (line) abline(0, 1, col = line.col, lty = line.lty, lwd = line.lwd)
    invisible(cbind(measured = x, predicted = as.vector(y)))
}


###
### Coefficient plot
###

coefplot <- function(object, ncomp = object$ncomp, comps, intercept = FALSE,
                     separate = FALSE, nCols, nRows, labels,
                     type = "l", lty = 1:nLines, lwd = NULL,
                     pch = 1:nLines, cex = NULL, col = 1:nLines, legendpos,
                     xlab = "variable", ylab = "regression coefficient",
                     main, pretty.xlabels = TRUE, xlim, ...)
{
    ## This simplifies code below:
    if (missing(comps)) comps <- NULL

    ## 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 (nRows * nCols < nPlots && dev.interactive()) par(ask = TRUE)
    } else {
        nCols <- nRows <- 1
    }
    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]]

    ## 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) {
                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, ...)
            } else {
                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, ...)
                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)
###

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:
plot.mvrVal <- function(x, nCols, nRows, type = "l", lty = 1:nEst, lwd = NULL,
                        pch = 1:nEst, cex = NULL, col = 1:nEst, legendpos,
                        xlab = "number of components", ylab = x$type, main,
                        ...)
{
    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 (nRows * nCols < nResp && dev.interactive()) par(ask = TRUE)
    } 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
###

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) && is.logical(xlabs) && !xlabs) # i.e. xlabs = FALSE
        mc$xlabs <- rep("o", nrow(objects))
    if (!missing(ylabs) && is.logical(ylabs) && !ylabs) # i.e. ylabs = FALSE
        mc$ylabs <- rep("o", nrow(vars))
    mc[[1]] <- as.name("biplot")
    ## Evaluate the call:
    eval(mc, parent.frame())
}