R/plot.R

Defines functions intDotplot plot.RRglmerMod plot.RRglm

Documented in intDotplot plot.RRglm plot.RRglmerMod

#' Plot diagnostics for a RRglm object
#'
#' Six plots (selectable by \code{which}) are currently available: (1) a plot of estimated population prevalence per RR model,
#' (2) a plot of estimated population prevalence per protection level,
#' (3) a plot of ungrouped residuals against predicted honest response,
#' (4) a plot of grouped (on covariates) residuals against predicted honest response,
#' (5) a plot of grouped Hosmer-Lemeshow residuals against predicted response,
#' and (6) a Normal Q-Q plot of grouped (on covariates) residuals. By default, plots 1, 3, 4 and 6 are provided.
#' Reference: Fox, J-P, Klotzke, K. and Veen, D. (2016). \emph{Generalized Linear Mixed Models for Randomized
#' Responses.} Manuscript submitted for publication.
#'
#' @param x
#' an object of class RRglm.
#' @param which
#' if a subset of the plots is required, specify a subset of the numbers 1:6 (default: 1, 3, 4, 6).
#' @param type
#' the type of residuals which should be used to be used for plots 3, 4 and 6. The alternatives are: "deviance" (default) and "pearson".
#' @param ngroups
#' the number of groups to compute the Hosmer-Lemeshow residuals for (default: 10).
#' @param ...
#' further arguments passed to or from other methods.
#' @method plot RRglm
#' @export
#' @examples
#' out <- RRglm(response ~ Gender + RR + pp + age, link="RRlink.logit", RRmodel=RRmodel,
#'          p1=RRp1, p2=RRp2, data=Plagiarism, etastart=rep(0.01, nrow(Plagiarism)))
#' plot(out, which = 1:6, type = "deviance", ngroups = 50)
plot.RRglm <- function(x, which = c(1, 3, 4, 6), type = c("deviance", "pearson"), ngroups = 10, ...)
{
  type <- match.arg(type)
  show <- rep(FALSE, 6)
  show[which] <- TRUE

  eta <- x$linear.predictors
  pi <- exp(eta) / (1 + exp(eta))

  if(any(show[1], show[2]))
  {
    s <- summary(x)
  }
  if(show[3])
  {
    if (type == "deviance")
      ylab.res <- "Deviance residuals"
    else if (type == "pearson")
      ylab.res <- "Pearson residuals"
  }
  if(any(show[4], show[6]))
  {
    vars <- all.vars(x$formula)
    unique.groups <- getUniqueGroups(x$model[, vars[2:length(vars)]])

    if (type == "deviance")
    {
      grouped.res <- residuals(x, type = "deviance.grouped")
      ylab.res <- "Deviance residuals"
      ylab.grouped.res <- "Deviance residuals (grouped)"
    }
    else if (type == "pearson")
    {
      grouped.res <- residuals(x, type = "pearson.grouped")
      ylab.res <- "Pearson residuals"
      ylab.grouped.res <- "Pearson residuals (grouped)"
    }
  }

  oask <- devAskNewPage(TRUE)
  on.exit(devAskNewPage(oask))

  if(show[1])
  {
    p <- s$prevalence.weighted
    tmp <- list()
    tmp$Prevalence <- data.frame("Prevalence.Weighted" = p$estimate.weighted)
    attr(tmp$Prevalence, "postVar") = array((p$se.weighted^2), dim=c(1,1,length(p$se.weighted)))
    rownames(tmp$Prevalence) <- paste(p$Item, " (", p$RRmodel, ")", sep="")

    dev.hold()
    ignore <- capture.output(print(intDotplot(tmp, mtit = "Prevalence (95%-CI)")))
    dev.flush()
  }
  if(show[2])
  {
    p <- s$prevalence
    tmp <- list()
    tmp$Prevalence <- data.frame("Prevalence" = p$estimate)
    attr(tmp$Prevalence, "postVar") = array(p$se^2, dim=c(1,1,length(p$se)))
    rownames(tmp$Prevalence) <- paste(p$Item, " (", p$model.level, ")", sep="")

    dev.hold()
    ignore <- capture.output(print(intDotplot(tmp, mtit = "Prevalence (95%-CI)")))
    dev.flush()
  }
  if(show[3])
  {
    res <- residuals(x, type = type)
    dev.hold()
    print(lattice::xyplot(res ~ pi, main = "Residuals vs. Predicted honest response", ylab = ylab.res,
         xlab = "Predicted affirmative honest response probability"))
    dev.flush()
  }
  if(show[4])
  {
    grouped.pi <- getCellMeans(y = pi, factor.groups = unique.groups)
    dev.hold()
    print(lattice::xyplot(grouped.res ~ grouped.pi, main = "Residuals vs. Predicted honest response (grouped)", ylab = ylab.grouped.res,
         xlab = "Predicted affirmative honest response probability"))
    dev.flush()
  }
  if(show[5])
  {
    grouped.res.hl <- residuals(x, type = "hosmer-lemeshow", ngroups = ngroups)
    dev.hold()
    print(lattice::xyplot(grouped.res.hl ~ 1:ngroups, main = "Residuals vs. Predicted honest response (HL)", ylab = "Hosmer-Lemeshow residuals",
         xlab = "Groups by predicted affirmative response probability"))
    dev.flush()
  }
  if(show[6])
  {
    dev.hold()
    print(lattice::qqmath( ~ grouped.res, xlab = "Theoretical Quantiles", ylab = ylab.grouped.res,
                           prepanel = lattice::prepanel.qqmathline,
                           panel = function(x, y) {
                             lattice::panel.qqmathline(x, distribution = qnorm)
                             lattice::panel.qqmath(x)
                           }))
    dev.flush()
  }

  invisible()
}

#' Plot diagnostics for a RRglmerMod object
#'
#' Five plots (selectable by \code{which}) are currently available: (1) a plot of estimated population prevalence per RR model,
#' (2) a plot of estimated population prevalence per protection level,
#' (3) a plot of random effects and their conditional variance (95%-cI),
#' (4) a plot of conditional pearson residuals against predicted honest response,
#' and (5) a plot of unconditional pearson residuals against predicted honest response.
#' By default, plots 1, 3, 4 and 5 are provided.
#' Reference: Fox, J-P, Klotzke, K. and Veen, D. (2016). \emph{Generalized Linear Mixed Models for Randomized
#' Responses.} Manuscript submitted for publication.
#'
#' @param x
#' an object of class RRglmerMod.
#' @param which
#' if a subset of the plots is required, specify a subset of the numbers 1:5 (default: 1, 3, 4, 5).
#' @param ...
#' further arguments passed to or from other methods.
#' @method plot RRglmerMod
#' @export
#' @examples
#' out <- RRglmer(response ~ Gender + RR + pp + (1+pp|age), link="RRlink.logit", RRmodel=RRmodel,
#'          p1=RRp1, p2=RRp2, data=Plagiarism, na.action = "na.omit",
#'          etastart = rep(0.01, nrow(Plagiarism)),
#'          control = glmerControl(optimizer = "Nelder_Mead", tolPwrss = 1e-03), nAGQ = 1)
#' plot(out, which = 1:5)
plot.RRglmerMod <- function(x, which = c(1, 3 , 4, 5), ...)
{
  #3: random effects
  #4: conditional pearson
  #5: unconditional pearson

  show <- rep(FALSE, 3)
  show[which] <- TRUE

  eta <- predict(x, type = "link")
  pi <- exp(eta) / (1 + exp(eta))

  if(any(show[1], show[2]))
  {
    s <- summary(x)
  }

  oask <- devAskNewPage(TRUE)
  on.exit(devAskNewPage(oask))
  if(show[1])
  {
    p <- s$prevalence.weighted
    tmp <- list()
    tmp$Prevalence <- data.frame("Prevalence.Weighted" = p$estimate.weighted)
    attr(tmp$Prevalence, "postVar") = array((p$se.weighted^2), dim=c(1,1,length(p$se.weighted)))
    rownames(tmp$Prevalence) <- paste(p$Item, " (", p$RRmodel, ")", sep="")

    dev.hold()
    ignore <- capture.output(print(intDotplot(tmp, mtit = "Prevalence (95%-CI)")))
    dev.flush()
  }
  if(show[2])
  {
    p <- s$prevalence
    tmp <- list()
    tmp$Prevalence <- data.frame("Prevalence" = p$estimate)
    attr(tmp$Prevalence, "postVar") = array(p$se^2, dim=c(1,1,length(p$se)))
    rownames(tmp$Prevalence) <- paste(p$Item, " (", p$model.level, ")", sep="")

    dev.hold()
    ignore <- capture.output(print(intDotplot(tmp, mtit = "Prevalence (95%-CI)")))
    dev.flush()
  }
  if(show[3])
  {
    random.effects <- ranef(x, condVar=TRUE)
    dev.hold()
    ignore <- capture.output(print(intDotplot(random.effects, mtit = "Random effects (95%-CI)")))
    dev.flush()
  }
  if(show[4])
  {
    res.conditional <- residuals(x, type = "pearson")
    ylab <- "Conditional pearson residuals"
    xlab <- "Predicted affirmative honest response probability"
    dev.hold()
    print(lattice::xyplot(res.conditional ~ pi, xlab = xlab, ylab = ylab))
    dev.flush()
  }
  if(show[5])
  {
    res.unconditional <- residuals(x, type = "unconditional.pearson")
    ylab <- "Unconditional pearson residuals"
    xlab <- "Predicted affirmative honest response probability"
    dev.hold()
    print(lattice::xyplot(res.unconditional ~ pi, xlab = xlab, ylab = ylab))
    dev.flush()
  }

  invisible()
}

#' Internal function to draw a dotplot
#'
#' Sligthly modified function dotplot.ranef.mer from lme4 package.
#'
#' @param x
#' data.
#' @param mtit
#' main title.
#' @keywords internal
#' @export
intDotplot <- function (x, mtit = "95%-CI", ...)
{
  prepanel.ci <- function(x, y, se, subscripts, ...) {
    if (is.null(se))
      return(list())
    x <- as.numeric(x)
    hw <- 1.96 * as.numeric(se[subscripts])
    list(xlim = range(x - hw, x + hw, finite = TRUE))
  }
  panel.ci <- function(x, y, se, subscripts, pch = 16, horizontal = TRUE,
                       col = dot.symbol$col, lty = dot.line$lty, lwd = dot.line$lwd,
                       col.line = dot.line$col, levels.fos = unique(y), groups = NULL,
                       ...) {
    x <- as.numeric(x)
    y <- as.numeric(y)
    dot.line <- lattice::trellis.par.get("dot.line")
    dot.symbol <- lattice::trellis.par.get("dot.symbol")
    sup.symbol <- lattice::trellis.par.get("superpose.symbol")
    lattice::panel.abline(h = levels.fos, col = col.line, lty = lty,
                 lwd = lwd)
    lattice::panel.abline(v = 0, col = col.line, lty = lty, lwd = lwd)
    if (!is.null(se)) {
      se <- as.numeric(se[subscripts])
      lattice::panel.segments(x - 1.96 * se, y, x + 1.96 * se, y,
                     col = "black")
    }
    lattice::panel.xyplot(x, y, pch = pch, ...)
  }
  f <- function(nx, ...) {
    xt <- x[[nx]]
    ss <- stack(xt)
    ss$ind <- factor(as.character(ss$ind), levels = colnames(xt))
    ss$.nn <- factor(rownames(xt))# rep.int(reorder(factor(rownames(xt)), xt[[1]],
                              #FUN = mean, sort = sort), ncol(xt))
    se <- NULL
    if (!is.null(pv <- attr(xt, "postVar")))
      se <- unlist(lapply(1:(dim(pv)[1]), function(i) sqrt(pv[i,
                                                              i, ])))
    lattice::dotplot(.nn ~ values | ind, ss, se = se, prepanel = prepanel.ci,
            panel = panel.ci, xlab = NULL, main = mtit, ...)
  }
  setNames(lapply(names(x), f, ...), names(x))
}

Try the GLMMRR package in your browser

Any scripts or data that you put into this service are public.

GLMMRR documentation built on May 30, 2017, 3:31 a.m.