fit.models: Compare Fitted Models

#' Compaison QQ Plot
#' 
#' Produces side-by-side QQ plots.  An optional simulated confidence envelope
#' can be included in each plot.
#' 
#' 
#' @param x a \code{fit.models} object.
#' @param fun a function to extract the desired quantity from \code{x}.
#' @param envelope a logical value.  If \code{TRUE} a \code{level} confidence
#' envelope is simulated for each QQ plot.
#' @param half.normal a logical value.  If \code{TRUE} the plot is drawn using
#' the absolute values.
#' @param n.samples a positive integer value giving the number of samples to
#' compute in the simulation of the confidence envelope.
#' @param level a numeric value between 0 and 1 specifying the confidence level
#' for the envelope.
#' @param id.n a non-negative integer value specifying the number of extreme
#' points to identify.
#' @param qqline a logical value.  If \code{TRUE}, a QQ line is included in the
#' plot.
#' @param \dots additional arguments are passed to
#' \code{\link[lattice]{xyplot}}.
#' @return the \code{trellis} object is invisibly returned.
#' @keywords hplot


#' @importFrom lattice xyplot panel.xyplot panel.abline llines strip.default
#' @importFrom stats rnorm quantile qnorm qqnorm


#' @export
sideBySideQQPlot <- function(x, fun, envelope = TRUE, half.normal = FALSE,
                             n.samples = 250, level = .95, id.n = 3, qqline = TRUE,
                             ...)
{
  confidence.envelope <- function(n, sd = 1, n.samples = 250, level = 0.95,
                                  half.normal = FALSE)
  {
    lower <- upper <- matrix(0.0, n, n.models)

    alphaOver2 <- (1.0 - level) / 2.0
    probs <- c(alphaOver2, 1.0 - alphaOver2)

    env <- matrix(rnorm(n * n.samples, sd = sd), n.samples, n)

    if(half.normal)
      env <- abs(env)

    env <- apply(env, 1, sort)
    env <- apply(env, 1, quantile, probs = probs)

    list(lower = env[1, ], upper = env[2, ])
  }

  n.models <- length(x)
  mod.names <- names(x)

  res <- lapply(x, fun)
  n.res <- sapply(res, length)

  px <- py <- list()
  for(i in 1:n.models) {
    tmp <- qqnorm(res[[i]], plot.it = FALSE)
    px[[i]] <- tmp$x
    py[[i]] <- tmp$y
  }

  if(half.normal) {
    py <- lapply(py, abs)
    px <- lapply(n.res, function (u) .5 + (0:(u-1)) / (2*u))
    px <- lapply(px, qnorm)
    for(i in 1:n.models)
      px[[i]][order(py[[i]])] <- px[[i]]
  }

  if(envelope) {
    sigma.hats <- numeric(n.models)
    for(i in 1:n.models) {
      if(!is.null(x[[i]]$scale))
        sigma.hats[i] <- x[[i]]$scale
      else {
        x.sum <- summary(x[[i]])
        if(!is.null(x.sum$dispersion))
          sigma.hats[i] <- sqrt(x.sum$dispersion)
        else if(!is.null(x.sum$sigma))
          sigma.hats[i] <- x.sum$sigma
        else
          stop("unable to determine residual scale")
      }
    }

    env <- list()
    for(i in 1:n.models)
      env[[i]] <- confidence.envelope(n.res[i], sigma.hats[i], n.samples,
                                      level, half.normal)
    lower <- lapply(env, function(u) u$lower)
    upper <- lapply(env, function(u) u$upper)

    den.range <- c(min(unlist(py), unlist(lower)),
                   max(unlist(py), unlist(upper)))
  }

  else
    den.range <- c(min(unlist(py)), max(unlist(py)))

  if(envelope && half.normal) {
    mod <- factor(rep(rep(mod.names, n.res), 3), levels = mod.names)

    tdf <- data.frame(py = c(unlist(py),
                             unlist(lower),
                             unlist(upper)),
                      px = rep(unlist(px), 3),
                      mod = mod)

    panel.special <- function(x, y, id.n, qqline, ...) {
      dat.idx <- 1:(length(x)/3)

      panel.xyplot(x[dat.idx], y[dat.idx], ...)

      if(qqline) {
        u <- quantile(x[!is.na(x)], c(0.25, 0.75))
        v <- quantile(y[!is.na(y)], c(0.25, 0.75))
        slope <- diff(v) / diff(u)
        int <- v[1] - slope * u[1]
        panel.abline(int, slope, ...)
      }

      panel.addons(x[dat.idx], y[dat.idx], id.n = id.n)

      dat.idx <- ((length(x)/3)+1):(2*length(x)/3)

      llines(sort(x[dat.idx]), sort(y[dat.idx]), col.line = "black", lty = 2)

      dat.idx <- (2*(length(x)/3)+1):(length(x))

      llines(sort(x[dat.idx]), sort(y[dat.idx]), col.line = "black", lty = 2)

      invisible()
    }
  }

  else if(envelope) {
    mod <- factor(rep(rep(mod.names, n.res), 3), levels = mod.names)

    tdf <- data.frame(py = c(unlist(py),
                             unlist(lower),
                             unlist(upper)),
                      px = rep(unlist(px), 3),
                      mod = mod)

    panel.special <- function(x, y, id.n, qqline, ...) {
      dat.idx <- 1:(length(x)/3)

      panel.xyplot(x[dat.idx], y[dat.idx], ...)

      panel.addons(x[dat.idx], y[dat.idx], id.n = id.n)

      if(qqline) {
        u <- quantile(x[!is.na(x)], c(0.25, 0.75))
        v <- quantile(y[!is.na(y)], c(0.25, 0.75))
        slope <- diff(v) / diff(u)
        int <- v[1] - slope * u[1]
        panel.abline(int, slope, ...)
      }

      dat.idx <- ((length(x)/3)+1):(2*length(x)/3)

      llines(sort(x[dat.idx]), sort(y[dat.idx]), col.line = "black", lty = 2)

      dat.idx <- (2*(length(x)/3)+1):(length(x))

      llines(sort(x[dat.idx]), sort(y[dat.idx]), col.line = "black", lty = 2)

      invisible()
    }
  }

  else {
    mod <- factor(rep(mod.names, n.res), levels = mod.names)
    tdf <- data.frame(px = unlist(px), py = unlist(py), mod = mod)

    panel.special <- function(x, y, id.n, qqline, ...) {
      panel.xyplot(x, y, ...)
      panel.addons(x, y, id.n = id.n)

      if(qqline) {
        u <- quantile(x[!is.na(x)], c(0.25, 0.75))
        v <- quantile(y[!is.na(y)], c(0.25, 0.75))
        slope <- diff(v) / diff(u)
        int <- v[1] - slope * u[1]
        panel.abline(int, slope, ...)
      }

      invisible()
    }
  }

  p <- xyplot(py ~ px | mod,
              data = tdf,
              id.n = id.n,
              qqline = qqline,
              panel = panel.special,
              strip = function(...) strip.default(..., style = 1),
              layout = c(n.models, 1, 1),
              ...)

  print(p)
  invisible(p)
}