R/merPlots.R

In merTools: Tools for Analyzing Mixed Effect Regression Models

#' @title Plot the results of a simulation of the random effects
#' @name plotREsim
#' @description Plot the simulated random effects on a ggplot2 chart. Points that
#' are distinguishable from zero (i.e. the confidence band based on \code{level}
#' does not cross the red line) are highlighted. Currently, the plots are ordered
#' according to the grouping factor.
#' @param data a data.frame generated by \code{\link{REsim}} with simulations of
#' the random effects of a \code{\link{merMod}}
#' @param level the width of the confidence interval
#' @param stat a character value indicating the variable name in data of the
#' midpoint of the estimated interval, e.g. "mean" or "median"
#' @param sd a logical indicating whether or not to plot error bars around
#' the estimates (default is TRUE). Calculates the width of the error bars
#' based on \code{level} and the variable named "sd" in \code{data}
#' @param sigmaScale a numeric value to divide the estimate and the standard
#' deviation by in the case of doing an effect size calculation
#' @param oddsRatio logical, should the parameters be converted to odds ratios
#' before plotting
#' @param labs logical, include the labels of the groups on the x-axis
#' @param facet Accepts either logical (\code{TRUE}) or \code{list} to specify which
#' random effects to plot. If \code{TRUE}, facets by both \code{groupFctr} and \code{term}.
#' If \code{list} selects the panel specified by the named elements of the list
#' @return a ggplot2 plot of the coefficient effects
#' @examples
#' \donttest{
#'  fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
#'  (p1 <- plotREsim(REsim(fm1)))
#'  #Plot just the random effects for the Days slope
#'  (p2 <- plotREsim(REsim(fm1), facet= list(groupFctr= "Subject", term= "Days")))
#'  }
#' @export
#' @import ggplot2
plotREsim <- function(data, level = 0.95, stat = "median", sd = TRUE,
                      sigmaScale = NULL, oddsRatio = FALSE, labs = FALSE,
                      facet= TRUE){
  # error checking
  plot_sim_error_chks(type= "RE", level= level, stat= stat, sd= sd, sigmaScale= sigmaScale,
                      oddsRatio= oddsRatio, labs= labs, facet= facet)
  # check for faceting
  facet_logical <- is.logical(facet)
  if (!facet_logical) {
    data <- data[data$groupFctr == facet[[1]] & data$term == facet[[2]], ]
  }

  if(!missing(sigmaScale)){
    data[, "sd"] <- data[, "sd"] / sigmaScale
    data[, stat] <- data[, stat] / sigmaScale
  }
  data[, "sd"] <- data[, "sd"] * qnorm(1-((1-level)/2))
  data[, "ymax"] <- data[, stat] + data[, "sd"]
  data[, "ymin"] <- data[, stat] - data[, "sd"]
  data[, "sig"] <- data[, "ymin"] > 0 | data[, "ymax"] < 0
  hlineInt <- 0
  if(oddsRatio == TRUE){
    data[, "ymax"] <- exp(data[, "ymax"])
    data[, stat] <- exp(data[, stat])
    data[, "ymin"] <- exp(data[, "ymin"])
    hlineInt <- 1
  }
  data <- data[order(data[,"groupFctr"], data[,"term"], data[,stat]),]
  rownames(data) <- 1:nrow(data)
  data[,"xvar"] <- factor(paste(data$groupFctr, data$groupID, sep=""),
                          levels=unique(paste(data$groupFctr,data$groupID, sep="")),
                          ordered=TRUE)
  if(labs == TRUE){
    xlabs.tmp <- element_text(face = "bold", angle=90, vjust=.5)
  } else {
    data[,"xvar"] <- as.numeric(data[,"xvar"])
    xlabs.tmp <- element_blank()
  }

  p <- ggplot(data, aes(x = .data[["xvar"]], y = .data[[stat]], ymax = .data[["ymax"]], ymin = .data[["ymin"]])) +
         geom_hline(yintercept = hlineInt, color = I("red"), linewidth = I(1.1)) +
         geom_point(color="gray75", alpha=1/(nrow(data)^.33), size=I(0.5)) +
         geom_point(data=subset(data, sig==TRUE), size=I(3)) +
         labs(x = "Group", y = "Effect Range", title = "Effect Ranges") +
         theme_bw() +
         theme(panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               axis.text.x = xlabs.tmp,
               axis.ticks.x = element_blank())
  if (sd) {
    p <- p +
      geom_pointrange(alpha = 1/(nrow(data)^.33)) +
      geom_pointrange(data=subset(data, sig==TRUE), alpha = 0.25)
  }
  # check facet
  if (facet_logical) {
    return(p + facet_grid(term ~ groupFctr, scales = "free_x"))
  } else {
    return(p)
  }
}

#' @title Plot the results of a simulation of the fixed effects
#' @name plotFEsim
#' @description Plot the simulated fixed effects on a ggplot2 chart
#' @param data a data.frame generated by \code{\link{FEsim}} with simulations of
#' the fixed effects of a \code{\link{merMod}}
#' @param level the width of the confidence interval
#' @param stat a character value indicating the variable name in data of the
#' midpoint of the estimated interval, e.g. "mean" or "median"
#' @param sd logical, indicating whether or not to plot error bars around
#' the estimates (default is TRUE). Calculates the width of the error bars
#' based on \code{level} and the variable named "sd" in \code{data}
#' @param intercept logical, should the intercept be included, default is FALSE
#' @param sigmaScale a numeric value to divide the estimate and the standard
#' deviation by in the case of doing an effect size calculation
#' @param oddsRatio logical, should the parameters be converted to odds ratios
#' before plotting
#' @return a ggplot2 plot of the coefficient effects
#' @examples
#'  fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
#'  (p1 <- plotFEsim(FEsim(fm1)))
#' @export
#' @import ggplot2
plotFEsim <- function(data, level=0.95, stat = "median", sd = TRUE,
                      intercept = FALSE, sigmaScale = NULL, oddsRatio = FALSE){
  # error checking
  plot_sim_error_chks(type= "FE", level= level, stat= stat, sd= sd, sigmaScale= sigmaScale,
                      oddsRatio= oddsRatio, labs= TRUE, facet= TRUE)

  if(!missing(sigmaScale)){
    data[, "sd"] <- data[, "sd"] / sigmaScale
    data[, stat] <- data[, stat] / sigmaScale
  }
  if(intercept == FALSE){
    data <- data[data$term != "(Intercept)", ]
  }
  data[, "sd"] <- data[, "sd"] * qnorm(1-((1-level)/2))
  data[, "ymax"] <- data[, stat] + data[, "sd"]
  data[, "ymin"] <- data[, stat] - data[, "sd"]
  hlineInt <- 0
  if(oddsRatio == TRUE){
    data[, "ymax"] <- exp(data[, "ymax"])
    data[, stat] <- exp(data[, stat])
    data[, "ymin"] <- exp(data[, "ymin"])
    hlineInt <- 1
  }
  xvar <- "term"
  data$term <- as.character(data$term)
  data$term <- factor(data$term , levels = data[order(data[, stat]), 1])
  p <- ggplot(aes(x = .data[[xvar]], y = .data[[stat]], ymax = .data[["ymax"]], ymin = .data[["ymin"]]), data = data) +
    geom_hline(yintercept = hlineInt, color = I("red")) +
    geom_point(size=I(3)) +
    coord_flip() +
    theme_bw()
  if (sd) {
    p <- p + geom_errorbar(width = 0.2)
  }
  p
}