R/plotcollindlnm.R

In collin: Visualization the Effects of Collinearity in Distributed Lag Models and Other Linear Models

#' Visualizes effects of collinearity in distributed lag model under an
#'   hypothetical linear effect pattern.
#'
#' Visualize the results from a distributed lag model under an hypothetical
#'   linear effect pattern provided by the user, generated using the function
#'   \code{\link{collindlnm}}.
#'
#' @param x an object of class \code{"collindlnmlinear"}, which is generated by
#'   the function \code{\link{collindlnm}} with \code{shape = "linear"}.
#' @param lags a number or a numeric vector indicating at what lags the results
#'   will be visualized. Default (\code{NULL}) shows all lags.
#' @param ... other parameters to be passed through to plot function.
#' @return A plot showing a comparison between the results under the fitted
#'   model and the results under the hypothetical true effect.
#' @seealso \code{\link{collindlnm}}, \code{\link{plot.collindlnmnonlinear}}.
#' @examples
#' # For detailed examples:
#' browseVignettes("collin")
#' @export
#' @importFrom graphics abline legend lines par plot
plot.collindlnmlinear <- function(x, lags = NULL, ...) {
  ### control "x"
  if (!inherits(x, "collindlnmlinear")) {
    stop("argument 'x' must be of class 'collindlnmlinear'")
  }

  ### control "lags"
  if (!is.null(lags) & (!inherits(lags, c("integer", "numeric")) || any(lags < 0))) {
    stop("Values in 'lags' must be non-negative integers.")
  }

  sims <- x$sim
  classmod <- x$classmod

  if (classmod == "lme") {
    pred <- x$pred$matfit
    h0 <- 0
  }
  if (classmod == "glm") {
    pred <- x$pred$matRRfit
    h0 <- 1
  }

  ### all lags
  nlags <- dim(sims)[2]
  alllags <- 0:(nlags - 1)

  ### lags required by the user
  userlags <- alllags
  if (!is.null(lags)) {
    userlags <- intersect(alllags, round(lags))
  }

  ### x-axis values
  xs <- userlags
  xrange <- range(xs)

  ### y-axis values
  if (!is.null(lags)) {
    tosel <- 1 + userlags
    sims <- sims[, tosel]
    pred <- pred[tosel]
  }
  yrange <- range(range(sims, na.rm = TRUE), range(pred, na.rm = TRUE))

  plot(xrange, yrange, type = "n", col = "gray70", ...)
  for (i in 1:nrow(sims)) {
    lines(xs, sims[i, ], col = "gray70")
  }
  lines(xs, pred, lwd = 3, col = "red")
  abline(h = h0, lty = 2, lwd = 2)
}


#' Visualize effects of collinearity in distributed lag model under an
#'   hypothetical non-linear effect pattern.
#'
#' Visualize the results from a distributed lag model under an hypothetical
#'   non-linear effect pattern provided by the user, generated using the function
#'   \code{\link{collindlnm}}. The number of plots shown is equal to the number
#'   of values passed by \code{at} in the function \code{\link{collindlnm}}.
#'   The way in which these plots are displayed is controlled by the user through
#'   the argument \code{show}.
#'
#' @param x an object of class \code{"collindlnmnonlinear"}, which is generated
#'   by the function \code{\link{collindlnm}} with \code{shape = "nonlinear"}.
#' @param lags a number or a numeric vector indicating at what lags the results
#'   will be visualized. Default (\code{NULL}) shows all lags.
#' @param show character indicating how the multiple plots will be shown. If
#'   \code{show = "manual"}, then it is expected that the user have previously
#'   set the graphical parameters to arrange them (i.e. setting \code{mfrow})
#'   with \code{\link{par}}. If \code{show = "auto"}, then the arrangement of
#'   the plots (i.e. the value of \code{mfrow}) is automatically set. If
#'   \code{show = "sequence"}, then the plots are sequentially overlaid.
#' @param addlegend logical indicating whether a legend indicating at what
#'   value (of \code{at} passed in \code{\link{collindlnm}}) the results
#'   correspond to.
#' @param varlegend character indicating the label for the explored variable to
#'   be shown in the legend.
#' @param ... other parameters to be passed through to plot function.
#' @return A plot showing a comparison between results under the fitted model
#'   and the results under the hypothetical true effect, for each of the
#'   different values of the variable of interest where effects were explored.
#' @seealso \code{\link{collindlnm}}, \code{\link{plot.collindlnmlinear}}.
#' @export
#' @importFrom graphics abline legend lines par plot
#' @importFrom grDevices dev.off devAskNewPage
plot.collindlnmnonlinear <- function(x,
                                     lags = NULL,
                                     show = c("manual", "auto", "sequence"),
                                     addlegend = TRUE,
                                     varlegend = NULL, ...) {
  ### control "x"
  if (!inherits(x, "collindlnmnonlinear")) {
    stop("argument 'x' must be of class 'collindlnmnonlinear'")
  }

  ### control "lags"
  if (!is.null(lags) & (!inherits(lags, c("integer", "numeric")) || any(lags < 0))) {
    stop("Values in 'lags' must be non-negative integers.")
  }

  ### control "show"
  show <- match.arg(show)

  ### control "addlegend"
  if (!is.logical(addlegend) || length(addlegend) != 1) {
    stop("'addlegend' must be TRUE or FALSE.")
  }

  ### control "varlegend"
  if (!is.null(varlegend) & (!inherits(varlegend, "character") || length(varlegend) != 1)) {
    stop("'varlegend' must be a character of length 1.")
  }

  sims <- x$sim
  classmod <- x$classmod
  at <- x$at

  if (classmod == "glm") {
    pred <- x$pred$matRRfit
    h0 <- 1
  }

  ### all lags
  nlags <-  dim(sims[[1]])[2]
  alllags <- 0:(nlags - 1)
  ### lags required by the user
  userlags <- alllags
  if (!is.null(lags)) {
    userlags <- intersect(alllags, round(lags))
  }

  ### x-axis values
  xs <- userlags
  xrange <- range(xs)

  ### legend
  if (addlegend) {
    if (is.null(varlegend))
      varlegend <- "at"
    leg <- paste(varlegend, "=", at)
  }

  ### plot
  if (show == "sequence") {
    dev.off()
  }

  nplots <- length(at)
  if (show == "auto") {
    oldpar <- par(no.readonly = TRUE) # all par settings which could be changed
    on.exit(par(oldpar))
    op <- par(mfrow = grDevices::n2mfrow(nplots))
  }

  nsim <- length(sims)

  for (j in 1:nplots) {
    simsatj <- t(sapply(sims, FUN = function(s) s[j, ]))
    predatj <- pred[j, ]
    # only those selected by the user
    if (!is.null(lags)) {
      tosel <- 1 + userlags
      simsatj <- simsatj[, tosel]
      predatj <- predatj[tosel]
    }
    yrange <- range(range(simsatj, na.rm = TRUE), range(predatj, na.rm = TRUE))

    if (show == "sequence") {
      devAskNewPage(TRUE)
    }

    plot(xrange, yrange, xlab = "Lag", ylab = "RR", type = "n", col = "gray70", ...)
    for (i in 1:nsim) {
      lines(xs, simsatj[i, ], col = "gray70")
    }
    lines(xs, predatj, lwd = 2, col = "red")
    abline(h = h0, lty = 2)
    #mtext("Lag", side = 1, line = 2, cex = 0.7)
    if (addlegend) {
      legend("topright", legend = leg[j])
    }
    devAskNewPage(FALSE)
  }
  devAskNewPage(options("device.ask.default")[[1]])
  if (show == "auto") {
    invisible() #-- now,  par(oldpar)  will be executed
  }
}