R/plot.vgam.R

In VGAM: Vector Generalized Linear and Additive Models

# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.














plotvgam <-
plot.vgam <-
  function(x, newdata = NULL,
           y = NULL, residuals = NULL,
           rugplot = TRUE,
           se = FALSE, scale = 0,
           raw = TRUE, offset.arg = 0,
           deriv.arg = 0, overlay = FALSE,
           type.residuals = c("deviance",
                              "working",
                              "pearson",
                              "response"),
           plot.arg = TRUE, which.term = NULL,
           which.cf = NULL,
           control = plotvgam.control(...),
           varxij = 1, ...) {

  missing.control <- missing(control)

  na.act <- x@na.action
  x@na.action <- list()

  if (!is.Numeric(varxij, integer.valued = TRUE,
                  length.arg = 1, positive = TRUE))
      stop("bad input for the 'varxij' argument")
  if (any(slotNames(x) == "control")) {
    x@control$varxij <- varxij
  }


  missing.type.residuals <- missing(type.residuals)
  if (mode(type.residuals) != "character" &&
      mode(type.residuals) != "name")
  type.residuals <-
    as.character(substitute(type.residuals))
  if (!missing.type.residuals)
  type.residuals <-
      match.arg(type.residuals,
                c("deviance", "working",
                  "pearson", "response"))[1]


  if (!is.Numeric(deriv.arg,
                  integer.valued = TRUE,
                  length.arg = 1) ||
      deriv.arg < 0)
    stop("bad input for the 'deriv' argument")

  if (se && deriv.arg > 0) {
    warning("standard errors not available ",
            "with derivatives. ",
            "Setting 'se = FALSE'")
    se <- FALSE
  }

      if (deriv.arg > 0 &&
          any(slotNames(x) == "ospsslot")) {
    stop("derivatives are not available for ",
         "objects ",
         " with 'sm.os()' or 'sm.ps()' terms")
  }

  preplot.object <- x@preplot
  if (!length(preplot.object)) {
      preplot.object <- preplotvgam(x,
                          newdata = newdata,
                          raw = raw,
                          deriv.arg = deriv.arg,
                          se = se,
                          varxij = varxij)
  }

  x@preplot <- preplot.object


    if (!is.null(residuals) &&
        length(residuals) == 1) {
    if (residuals) {
      if (missing.type.residuals) {
        for (rtype in type.residuals)
          if (!is.null(residuals <-
              resid(x, type = rtype))) break
      } else {
       residuals = resid(x, type = type.residuals)
        if (!length(residuals))
          warning("residuals are NULL. Ignoring",
                  " 'residuals = TRUE'")
      }
    } else {
      residuals <- NULL
    }
  }

  if (!missing.control) {
    control <-
      c(plotvgam.control( .include.dots = FALSE,
                         ...),
        control, plotvgam.control(...))
  }

  x@post$plotvgam.control <- control

  if (plot.arg)
    plotpreplotvgam(preplot.object,
                    residuals = residuals,
                    rugplot = rugplot,
                    scale = scale, se = se,
                    offset.arg = offset.arg,
                    deriv.arg = deriv.arg,
                    overlay = overlay,
                    which.term = which.term,
                    which.cf = which.cf,
                    control = control)

  x@na.action <- na.act  # Restore its orig. value
  invisible(x)
}  # plotvgam and plot.vgam






ylim.scale <- function(ylim, scale = 0) {
  if (length(ylim) != 2 ||
      ylim[2] < ylim[1])
    stop("error in 'ylim'")
  try <- ylim[2] - ylim[1]
  if (try > scale) ylim else
    c(ylim[1] + ylim[2] - scale,
      ylim[1] + ylim[2] + scale) / 2
}






getallresponses <- function(xij) {
  if (!is.list(xij))
    return("")

  allterms <- lapply(xij, terms)
  allres <- NULL
  for (ii in seq_along(xij))
    allres <- c(allres,
         as.character(attr(allterms[[ii]], "variables"))[2])
  allres
}






 headpreplotvgam <-
  function(object, newdata = NULL,
           terms = attr((object@terms)$terms,
                        "term.labels"),
           raw = TRUE, deriv.arg = deriv.arg,
           se = FALSE,
           varxij = 1) {
  Terms <- terms(object)  # 20030811; object@terms$terms
  aa <- attributes(Terms)
  all.terms <- labels(Terms)
  xvars <- parse(text = all.terms)




  names(xvars) <- all.terms
  terms <- sapply(terms, match.arg, all.terms)

  Interactions <- aa$order > 1
  if (any(Interactions)) {
    stop("cannot handle interactions")
  }

  xvars <- xvars[terms]
  xnames <- as.list(terms)
  names(xnames) <- terms
  modes <- sapply(xvars, mode)
  for (term in terms[modes != "name"]) {
    evars <- all.names(xvars[term], functions = FALSE,
                       unique = TRUE)
    if (!length(evars))
      next
    xnames[[term]] <- evars
    evars <- parse(text=evars)
    if (length(evars) == 1) {
      evars <- evars[[1]]
    } else
    if (length(evars) > 1 &&
       length(intersect(getallresponses(object@control$xij),
                        names(xnames)))) {




      evars <- evars[[varxij]]
    } else {
      evars <- c(as.name("list"), evars)
      mode(evars) <- "call"
    }
    xvars[[term]] <- evars
  }


  xvars <- c(as.name("list"), xvars)
  mode(xvars) <- "call"
  if (length(newdata)) {
    xvars <- eval(xvars, newdata)
  } else {
    Call <- object@call
    if (!is.null(Call$subset) | !is.null(Call$na.action) |
        !is.null(options("na.action")[[1]])) {
      Rownames <- names(fitted(object))
      if (!(Rl <- length(Rownames)))
        Rownames <- dimnames(fitted(object))[[1]]

      if (length(object@x) && !(Rl <- length(Rownames)))
        Rownames <- (dimnames(object@x))[[1]]
      if (length(object@y) && !(Rl <- length(Rownames)))
        Rownames <- (dimnames(object@y))[[1]]

      if (!(Rl <- length(Rownames)))
        stop("need to have names for fitted.values ",
             "when call has a 'subset' or 'na.action' argument")

        form <- paste("~", unlist(xnames),
                      collapse = "+")
        Mcall <- c(as.name("model.frame"),
                   list(formula =
                 terms(as.formula(form)),
                 subset = Rownames,
                 na.action = function(x) x))
      mode(Mcall) <- "call"
      Mcall$data <- Call$data
      xvars <- eval(xvars, eval(Mcall))
    } else {
      ecall <- substitute(eval(expression(xvars)))
      ecall$local <- Call$data
      xvars <- eval(ecall)
    }
  }
  list(xnames = xnames, xvars = xvars)
}  # headpreplotvgam






preplotvgam <-
  function(object, newdata = NULL,
           terms = attr((object@terms)$terms,
                        "term.labels"),
           raw = TRUE, deriv.arg = deriv.arg,
           se = FALSE,
           varxij = 1) {

  result1 <- headpreplotvgam(object,
                             newdata = newdata,
                             terms = terms,
                             raw = raw,
                             deriv.arg = deriv.arg,
                             se = se,
                             varxij = varxij)

  xvars  <- result1$xvars
  xnames <- result1$xnames



  if (FALSE && !is.null(object@control$jix)) {




    myxij <- object@control$xij
    if (length(myxij)) {
    }
  }

  pred <- if (length(newdata)) {
    predict(object, newdata, type = "terms",
            raw = raw, se.fit = se,
            deriv.arg = deriv.arg)
  } else {
    predict(object, type = "terms",
            raw = raw, se.fit = se,
            deriv.arg = deriv.arg)
  }

  fits <- if (is.atomic(pred))
            NULL else pred$fit
  se.fit <- if (is.atomic(pred))
              NULL else pred$se.fit

  if (is.null(fits))
    fits <- pred
  fred <- attr(fits,
               "vterm.assign")  # NULL for M==1
  Constant <- attr(fits,
                   "constant")  # NULL if se = T

  gamplot <- xnames

  loop.var <- names(fred)
  for (term in loop.var) {
    .VGAM.x <- xvars[[term]]

      myylab <-
        if (all(substring(term, 1:nchar(term),
                          1:nchar(term)) != "("))
          paste("partial for", term) else term

    TT <- list(x = .VGAM.x,
      y = fits[, (if (is.null(fred)) term else
                  fred[[term]])],
          se.y = if (is.null(se.fit)) NULL else
          se.fit[, (if (is.null(fred)) term else
                            fred[[term]])],
          xlab = xnames[[term]],
          ylab = myylab)
    class(TT) <- "preplotvgam"
    gamplot[[term]] <- TT
  }
  attr(gamplot, "Constant") <- Constant
  invisible(gamplot)
}  # preplotvgam






 plotpreplotvgam <-
  function(x, y = NULL, residuals = NULL,
           rugplot = TRUE, se = FALSE,
           scale = 0,
           offset.arg = 0, deriv.arg = 0,
           overlay = FALSE,
           which.term = NULL, which.cf = NULL,
           control = NULL) {
  listof <- inherits(x[[1]], "preplotvgam")
  if (listof) {
    TT <- names(x)
    if (is.null(which.term))
      which.term <- TT  # Plot them all

    if (is.character(control$main))
      control.main.save <- rep(control$main,
                           length = length(TT))

    if (deriv.arg > 0 &&
        is.character(which.term)) {
        if (length(index.fun.call <-
            grep("[(]", which.term))) {
      terms2check <- which.term[index.fun.call]
      if (!any(substr(terms2check,
                      1, 2) == "s(")) {
warning("there appears to be no s() term, ",
        "so setting ",
  "argument 'deriv.arg' a positive value has ",
  "no effect. Setting its value to 0.")
        deriv.arg <- 0  # Replacing its value
      }
    }
  }

 

    plot.no <- 0
    for (ii in TT) {
      plot.no <- plot.no + 1
      control$main <- control.main.save[plot.no]
      if ((is.character(which.term) &&
           any(which.term == ii)) ||
          (is.numeric(which.term) &&
           any(which.term == plot.no)))
        plotpreplotvgam(x[[ii]], y = NULL,
                        residuals,
                        rugplot = rugplot,
                        se = se,
                        scale = scale,
                        offset.arg = offset.arg,
                        deriv.arg = deriv.arg,
                        overlay = overlay,
                        which.cf = which.cf,
                        control = control)
    }  #  ii
  } else {
    dummy <-
      function(residuals = NULL, rugplot = TRUE,
               se = FALSE, scale = 0,
               offset.arg = 0, deriv.arg = 0,
               overlay = FALSE,
               which.cf = NULL,
               control = plotvgam.control())
          c(list(residuals = residuals,
                 rugplot = rugplot,
            se = se, scale = scale,
            offset.arg = offset.arg,
            deriv.arg = deriv.arg,
            overlay = overlay,
            which.cf = which.cf), control)

      dd <- dummy(residuals = residuals,
                  rugplot = rugplot,
                se = se, scale = scale,
                offset.arg = offset.arg,
                deriv.arg = deriv.arg,
                overlay = overlay,
                which.cf = which.cf,
                control = control)

    uniq.comps <- unique(c(names(x), names(dd)))
    Call <- c(as.name("vplot"),
              c(dd, x)[uniq.comps])
    mode(Call) <- "call"
    invisible(eval(Call))
  }
}  # plotpreplotvgam






vplot.default <-
    function(x, y, se.y = NULL, xlab = "",
             ylab = "",
           residuals = NULL, rugplot = FALSE,
           scale = 0, se = FALSE,
           offset.arg = 0, deriv.arg = 0,
           overlay = FALSE,
           which.cf = NULL, ...) {
  switch(data.class(x)[1],
         logical = vplot.factor(factor(x),
                        y, se.y, xlab,
                        ylab,
                        residuals,
                        rugplot, scale,
                        se,
                        offset.arg = offset.arg,
                        overlay = overlay, ...),
         if (is.numeric(x)) {
             vplot.numeric(as.vector(x), y, se.y,
                           xlab, ylab,
                           residuals, rugplot,
                           scale, se,
                         offset.arg = offset.arg,
                         overlay = overlay, ...)
         } else {
           warning("The 'x' component of '",
             ylab, "' has class '",
             class(x),
             "'; no vplot() methods available")
         }
        )  # End of switch
}



vplot.list <-
  function(x, y, se.y = NULL, xlab, ylab,
           residuals = NULL, rugplot = FALSE,
           scale = 0, se = FALSE,
           offset.arg = 0, deriv.arg = 0,
           overlay = FALSE,
           which.cf = NULL, ...) {

  if (is.numeric(x[[1]])) {
    vplot.numeric(x[[1]], y, se.y, xlab, ylab,
                  residuals, rugplot, scale, se,
                  offset.arg = offset.arg,
                  deriv.arg = deriv.arg,
                  overlay = overlay, ...)
  } else {
    stop("this function has not been written yet")
  }
}  # vplot.list






 plotvgam.control <-
  function(which.cf = NULL,
           xlim = NULL, ylim = NULL,
           llty = par()$lty,
           slty = "dashed",
           pcex = par()$cex,
           pch = par()$pch,
           pcol = par()$col,
           lcol = par()$col,
           rcol = par()$col,
           scol = par()$col,
           llwd = par()$lwd,
           slwd = par()$lwd,
           add.arg = FALSE,
           one.at.a.time = FALSE,
           .include.dots = TRUE,
           noxmean = FALSE,
           shade = FALSE, shcol = "gray80",
           main = "",   # NULL,
           ...) {


  ans <-
  list(which.cf = which.cf,
       xlim = xlim, ylim = ylim,
       llty = llty, slty = slty,
       pcex = pcex, pch = pch,
       pcol = pcol, lcol = lcol,
       rcol = rcol, scol = scol,
       llwd = llwd, slwd = slwd,
       add.arg = add.arg,
       noxmean = noxmean,
       one.at.a.time = one.at.a.time,
       main = main,
       shade = shade, shcol = shcol)

  if (.include.dots) {
    c(list(...), ans)
  } else {
    default.vals <- plotvgam.control()
    return.list <- list()
    for (ii in names(default.vals)) {
    replace.val <-
        !((length(ans[[ii]]) ==
           length(default.vals[[ii]])) &&
        (length(default.vals[[ii]]) > 0) &&
        identical(ans[[ii]], default.vals[[ii]]))

      if (replace.val)
        return.list[[ii]] <- ans[[ii]]
    }
    if (length(return.list)) {
      names(return.list) <- names(return.list)
      return.list
    } else {
      NULL
    }
  }
}  # plotvgam.control






vplot.numeric <-
  function(x, y, se.y = NULL, xlab, ylab,
           residuals = NULL, rugplot = FALSE,
           se = FALSE, scale = 0,
           offset.arg = 0, deriv.arg = 0,
           overlay = FALSE,
           which.cf = NULL,
           xlim = NULL, ylim = NULL,
           llty = par()$lty,
           slty = "dashed",
           pcex = par()$cex,
           pch = par()$pch,
           pcol = par()$col,
           lcol = par()$col,
           rcol = par()$col,
           scol = par()$col,
           llwd = par()$lwd,
           slwd = par()$lwd,
           add.arg = FALSE,
           one.at.a.time = FALSE,
           noxmean = FALSE,
           separator = ":",
           shade = FALSE, shcol = "gray80",
           main = "",
           ...) {




    ylim0 <- ylim

    if (length(y) / length(x) !=
        round(length(y) / length(x)))
      stop("length of 'x' and 'y' do not ",
           "seem to match")
    y <- as.matrix(y)
    if (!length(which.cf))
      which.cf <- 1:ncol(y)  # Added 20040807

    if (!is.null(se.y))
      se.y <- as.matrix(se.y)
    if (!is.null(se.y) && anyNA(se.y))
      se.y <- NULL

    if (!is.null(residuals))  {
      residuals <- as.matrix(residuals)
      if (ncol(residuals) != ncol(y)) {
        warning("ncol(residuals) != ncol(y) so",
                " residuals are not plotted")
        residuals <- NULL
      }
    }

    offset.arg <- matrix(offset.arg, nrow(y),
                         ncol(y), byrow = TRUE)
    y <- y + offset.arg

    ylab <- add.hookey(ylab, deriv.arg)


    if (xmeanAdded <-
       (se && !is.null(se.y) && !noxmean &&
        all(substring(ylab, 1:nchar(ylab),
                      1:nchar(ylab)) != "("))) {
      x <- c(x, mean(x))
      y <- rbind(y, 0 * y[1, ])
      se.y <- rbind(se.y, 0 * se.y[1, ])
      if (!is.null(residuals))
        residuals <- rbind(residuals,
          NA * residuals[1, ])  # NAs not plotted
    }

    ux <- unique(sort(x))
    ooo <- match(ux, x)
    uy <- y[ooo, , drop = FALSE]


    xlim.orig <- xlim
    ylim.orig <- ylim
    xlim <- range(if (length(xlim))
                    NULL else ux,
                  xlim, na.rm = TRUE)
    ylim <- range(if (length(ylim))
                  NULL else uy[, which.cf],
                  ylim, na.rm = TRUE)


    if (rugplot) {
      usex <- if (xmeanAdded) x[-length(x)] else x
      jx <- jitter(usex[!is.na(usex)])
      xlim <- range(if (length(xlim.orig))
                    NULL else jx,
                    xlim.orig, na.rm = TRUE)
    }

    if (se && !is.null(se.y)) {
      se.upper <- uy + 2 * se.y[ooo, ,
                                drop = FALSE]
      se.lower <- uy - 2 * se.y[ooo, ,
                                drop = FALSE]

      ylim <- if (length(ylim.orig))
              range(ylim.orig) else
              range(c(ylim, se.upper[, which.cf],
                      se.lower[, which.cf]))
    }

    if (!is.null(residuals)) {
      if (length(residuals) == length(y)) {
        residuals <- as.matrix(y + residuals)
        ylim <- if (length(ylim.orig))
                    range(ylim.orig) else
           range(c(ylim, residuals[, which.cf]),
                      na.rm = TRUE)
      } else {
        residuals <- NULL
        warning("Residuals do not match 'x' in",
                " '", ylab, "' preplot object")
      }
    }


  all.missingy <- all(is.na(y))

  if (all.missingy)
    return()

  if (!length(ylim.orig))
    ylim <- ylim.scale(ylim, scale)

  if (overlay) {
    if (!length(which.cf))
      which.cf <- 1:ncol(uy)  # Added 20040807
    if (!add.arg) {
      matplot(ux, uy[, which.cf], type = "n",
              xlim = xlim, ylim = ylim,
              xlab = xlab, ylab = ylab,
              main = main, ...)
    }
    matlines(ux, uy[, which.cf],
             lwd = llwd, col = lcol, lty = llty)
    if (!is.null(residuals)) {
      if (ncol(y) == 1) {
          points(x, residuals, pch = pch,
                 col = pcol, cex = pcex)
      } else {
        matpoints(x, residuals[, which.cf],
                  pch = pch, col = pcol,
                  cex = pcex)  # add.arg = TRUE,
      }
    }
    if (rugplot)
      rug(jx, col = rcol)
    if (se && !is.null(se.y)) {
        matlines(ux, se.upper[, which.cf],
                 lty =  slty,
               lwd = slwd, col = scol)
        matlines(ux, se.lower[, which.cf],
                 lty =  slty,
               lwd = slwd, col = scol)
    }
  } else {
    YLAB <- ylab

    pcex <- rep_len(pcex,  ncol(uy))
    pch  <- rep_len(pch ,  ncol(uy))
    pcol <- rep_len(pcol,  ncol(uy))
    lcol <- rep_len(lcol,  ncol(uy))
    llty <- rep_len(llty,  ncol(uy))
    llwd <- rep_len(llwd,  ncol(uy))
    slty <- rep_len(slty,  ncol(uy))
    rcol <- rep_len(rcol,  ncol(uy))
    scol <- rep_len(scol,  ncol(uy))
    slwd <- rep_len(slwd,  ncol(uy))

    for (ii in 1:ncol(uy)) {
      if (!length(which.cf) ||
          ( length(which.cf) &&
            any(which.cf == ii))) {

        if (is.Numeric(ylim0, length.arg = 2)) {
          ylim <- ylim0
        } else {
          ylim <- range(ylim0, uy[, ii],
                        na.rm = TRUE)
          if (se && !is.null(se.y))
            ylim <- range(ylim0, se.lower[, ii],
                          se.upper[, ii],
                          na.rm = TRUE)
          if (!is.null(residuals))
              ylim <- range(c(ylim,
                              residuals[, ii]),
                          na.rm = TRUE)
          ylim <- ylim.scale(ylim, scale)
        }
        if (ncol(uy) > 1 && length(separator))
            YLAB <- paste(ylab, separator, ii,
                          sep = "")

        if (!add.arg) {
          if (one.at.a.time) {
     readline("Hit return for the next plot ")
          }
           plot(ux, uy[, ii], type = "n",
                xlim = xlim, ylim = ylim,
                xlab = xlab, ylab = YLAB,
                main = main, ...)
        }

        lines(ux, uy[, ii],
              lwd = llwd[ii], col = lcol[ii],
              lty = llty[ii])
        if (!is.null(residuals))
          points(x, residuals[, ii],
                 pch = pch[ii],
                 col = pcol[ii], cex = pcex[ii])
        if (rugplot)
          rug(jx, col = rcol[ii])

        if (se && !is.null(se.y)) {
          if (shade) {
            polygon(c(ux, rev(ux), ux[1]),
                    c(se.upper[, ii],
                      rev(se.lower[, ii]),
                      se.upper[1, ii]),
                    col = shcol, border = NA)
            lines(ux, uy[, ii],
                  lwd = llwd[ii], col = lcol[ii],
                  lty = llty[ii])
          } else {
            lines(ux, se.upper[, ii],
                  lty = slty[ii],
                  lwd = slwd[ii], col = scol[ii])
            lines(ux, se.lower[, ii],
                  lty = slty[ii],
                  lwd = slwd[ii], col = scol[ii])
          }  # !shade
        }  # se && !is.null(se.y))
      }
    }  # for()
  }  # overlay
}  # vplot.numeric()






vplot.matrix <-
  function(x, y, se.y = NULL, xlab, ylab,
           residuals = NULL, rugplot = FALSE,
           scale = 0, se = FALSE,
           offset.arg = 0, deriv.arg = 0,
           overlay = FALSE,
           which.cf = NULL, ...) {
  stop("You should not ever call this function!")
}






add.hookey <- function(ch, deriv.arg = 0) {

    if (!is.Numeric(deriv.arg,
                    integer.valued = TRUE,
                  length.arg = 1) ||
      deriv.arg < 0)
      stop("bad input for the 'deriv' argument")

  if (deriv.arg == 0)
    return(ch)

  hookey <- switch(deriv.arg, "'", "''",
                   "'''", "''''",
                   "'''''",
                   stop("too high a derivative"))
  nc <- nchar(ch)
  sub <- substring(ch, 1:nc, 1:nc)
  if (nc >= 2 && sub[1] == "s" && sub[2] == "(") {
    paste("s", hookey, substring(ch, 2, nc),
          sep = "", coll = "")
  } else {
    paste(ch, hookey, sep = "", collapse = "")
  }
}  # add.hookey






vplot.factor <-
  function(x, y, se.y = NULL, xlab, ylab,
           residuals = NULL,
           rugplot = FALSE, scale = 0,
           se = FALSE, xlim = NULL, ylim = NULL,
           offset.arg = 0, deriv.arg = 0,
           overlay = FALSE,
           which.cf = NULL, ...) {
  if (deriv.arg > 0)
    return(NULL)

  if (length(y)/length(x)  !=
      round(length(y)/length(x)))
    stop("length of 'x' and 'y' do not seem ",
         "to match")
  y <- as.matrix(y)

  if (!is.null(se.y))
    se.y <- as.matrix(se.y)
  if (!is.null(se.y) && anyNA(se.y))
    se.y <- NULL

  if (!is.null(residuals))  {
    residuals <- as.matrix(residuals)
    if (ncol(residuals) != ncol(y)) {
      warning("ncol(residuals) != ncol(y) so ",
              "residuals are not plotted")
      residuals <- NULL
    }
  }

    if (overlay) {
      warning("overlay = TRUE in vplot.factor: ",
              "assigning overlay <- FALSE")
      vvplot.factor(x, y,
                    se.y = if (is.null(se.y))
                               NULL else se.y,
                    xlab = xlab, ylab = ylab,
                    residuals = residuals,
                    rugplot = rugplot,
                    scale = scale,
                    se = se, xlim = xlim,
                    ylim = ylim, ...)
  } else {
    for (ii in 1:ncol(y)) {
      ylab <- rep_len(ylab, ncol(y))
      if (ncol(y) > 1)
        ylab <- dimnames(y)[[2]]
      vvplot.factor(x, y[, ii,drop = FALSE],
        se.y = if (is.null(se.y)) NULL else
                 se.y[, ii,drop = FALSE],
        xlab = xlab, ylab = ylab[ii],
        residuals = if (is.null(residuals))
          NULL else residuals[, ii,drop = FALSE],
        rugplot = rugplot, scale = scale,
        se = se, xlim = xlim, ylim = ylim, ...)
    }
  }
  invisible(NULL)
}  # vplot.factor






vvplot.factor <-
  function(x, y, se.y = NULL, xlab, ylab,
           residuals = NULL, rugplot = FALSE,
           scale = 0,
           se = FALSE, xlim = NULL, ylim = NULL,
           ...) {

  M <- ncol(y)
  nn <- as.numeric(table(x))
  codex <- as.numeric(x)
  ucodex <- seq(nn)[nn > 0]
  ooo <- match(ucodex, codex, 0)

  uy <- y[ooo, , drop = FALSE]
  ylim <- range(ylim, uy)
  xlim <- range(c(0, sum(nn), xlim))
  rightx <- cumsum(nn)
  leftx <- c(0, rightx[ -length(nn)])
  ux <- (leftx + rightx)/2
  delta <- (rightx - leftx)/8

  jx <- runif(length(codex),
             (ux - delta)[codex],
             (ux + delta)[codex])
  nnajx <- jx[!is.na(jx)]

  if (rugplot)
    xlim <- range(c(xlim, nnajx))
  if (se && !is.null(se.y)) {
    se.upper <- uy + 2 * se.y[ooo, , drop = FALSE]
    se.lower <- uy - 2 * se.y[ooo, , drop = FALSE]
    ylim <- range(c(ylim, se.upper, se.lower))
  }
  if (!is.null(residuals)) {
    if (length(residuals) == length(y)) {
      residuals <- y + residuals
      ylim <- range(c(ylim, residuals))
    } else {
      residuals <- NULL
      warning("Residuals do not match 'x' in \"",
              ylab, "\" preplot object")
    }
  }
  ylim <- ylim.scale(ylim, scale)
  Levels <- levels(x)
  if (any(nn == 0)) {
    keep <- nn > 0
    nn <- nn[keep]
    ux <- ux[keep]
    delta <- delta[keep]
    leftx <- leftx[keep]
    rightx <- rightx[keep]
    Levels <- Levels[keep]
  }


  about <- function(ux, M, Delta = 1 / M) {
    if (M == 1) return(cbind(ux))
    ans <- matrix(NA_real_, length(ux), M)
    grid <- seq(-Delta, Delta, len = M)
    for (ii in 1:M) {
      ans[, ii] <- ux + grid[ii]
    }
    ans
  }

  uxx <- about(ux, M, Delta = min(delta))
  xlim <- range(c(xlim, uxx))

      
  matplot(ux, uy, ylim = ylim, xlim = xlim,
          xlab = "", type = "n",
          ylab = ylab, axes = FALSE,
          frame.plot = TRUE)  # , ...
  mtext(xlab, 1, 2, adj = 0.5)
  axis(side = 2)
  lpos <- par("mar")[3]
  mtext(Levels, side = 3, line = lpos/2,
        at = ux, adj = 0.5, srt = 45)

  for (ii in 1:M)
    segments(uxx[, ii] - 1.0 * delta, uy[, ii],
             uxx[, ii] + 1.0 * delta, uy[, ii])
  if (!is.null(residuals)) {
    for (ii in 1:M) {
      jux <- uxx[, ii]
      jux <- jux[codex]
      jux <- jux + runif(length(jux),
                         -0.7*min(delta),
                         0.7*min(delta))
    if (M == 1) points(jux, residuals[, ii]) else
                points(jux, residuals[, ii],
                       pch = as.character(ii))
    }
  }
  if (rugplot)
    rug(nnajx)
  if (se) {
    for (ii in 1:M) {
        segments(uxx[, ii] + 0.5*delta,
                 se.upper[, ii],
                 uxx[, ii] - 0.5*delta,
                 se.upper[, ii])
        segments(uxx[, ii] + 0.5*delta,
                 se.lower[, ii],
                 uxx[, ii] - 0.5*delta,
                 se.lower[, ii])
      segments(uxx[, ii], se.lower[, ii],
               uxx[, ii], se.upper[, ii],
               lty = 2)
    }
  }
  invisible(diff(ylim))
}  # vvplot.factor






if (!isGeneric("vplot"))
    setGeneric("vplot", function(x, ...)
        standardGeneric("vplot"))


setMethod("vplot", "factor", function(x, ...)
           vplot.factor(x, ...))
setMethod("vplot", "list", function(x, ...)
           vplot.list(x, ...))
setMethod("vplot", "matrix", function(x, ...)
           vplot.matrix(x, ...))
setMethod("vplot", "numeric", function(x, ...)
           vplot.numeric(x, ...))




setMethod("plot", "vgam",
         function(x, y, ...) {
         if (!missing(y))
         stop("cannot process the 'y' argument")
         invisible(plot.vgam(x = x, y = y, ...))})









 plotqrrvglm <-
  function(object,
           rtype = c("response", "pearson",
                     "deviance", "working"),
           ask = FALSE,
           main = paste(Rtype,
             "residuals vs latent variable(s)"),
           xlab = "Latent Variable",
     I.tolerances = object@control$eq.tolerances,
           ...) {
  M <- object@misc$M
  n <- object@misc$n
  Rank <- object@control$Rank
  Coef.object <- Coef(object,
                      I.tolerances = I.tolerances)
  rtype <- match.arg(rtype,
                     c("response", "pearson",
                       "deviance", "working"))[1]
  res <- resid(object, type = rtype)

  my.ylab <- if (length(object@misc$ynames))
                 object@misc$ynames else
                 rep_len(" ", M)
  Rtype <- switch(rtype, pearson = "Pearson",
                  response = "Response",
                  deviance = "Deviance",
                  working = "Working")

  done <- 0
  for (rr in 1:Rank)
    for (ii in 1:M) {
      plot(Coef.object@latvar[, rr],
           res[, ii],
           xlab = paste0(xlab, if (Rank == 1)
                                   "" else rr),
           ylab = my.ylab[ii],
           main = main, ...)
      done <- done + 1
      if (done >= prod(par()$mfrow) && ask &&
          done != Rank*M) {
          done <- 0
     readline("Hit return for the next plot: ")
      }
    }
  object
}  # plotqrrvglm



setMethod("plot", "qrrvglm", function(x, y, ...)
         invisible(plotqrrvglm(object = x, ...)))









put.caption <-
    function(text.arg = "(a)",
             w.x = c(0.50, 0.50),
             w.y = c(0.07, 0.93), ...) {
  text(text.arg,
       x = weighted.mean(par()$usr[1:2], w = w.x),
       y = weighted.mean(par()$usr[3:4], w = w.y),
       ...)
}












setMethod("plot", "pvgam",
   function(x, y, ...) {
     if (!missing(y))
       stop("cannot process the 'y' argument")
       invisible(plot.vgam(x = x, y = y, ...))})