R/marginalplot.R

In eirm: Explanatory Item Response Modeling for Dichotomous and Polytomous Items

#' @title Plots of Estimated Marginal Effects in Explanatory IRT Models
#' @name marginalplot
#' @importFrom ggeffects ggpredict get_x_title get_legend_title
#' @importFrom ggplot2 ggplot aes_string geom_point position_dodge geom_errorbar scale_y_continuous labs theme_bw facet_wrap
#' @description
#' This function uses \code{\link[ggeffects]{ggpredict}} to calculate marginal effects for explanatory variables in
#' an explanatory IRT model estimated with the \code{\link{eirm}} function. It returns a plot of estimated probabilities
#' generated by the explanatory IRT model while holding some predictors constant.
#'
#' @param x An eirm object returned from the \code{\link{eirm}} function.
#' @param predictors Character vector with the names of up to three categorical predictors from the eirm model. The first predictor
#' is plotted on the x-axis; the second predictor is used as a group variable; the third predictor is used as
#' a facet in the plot.
#' @param conf.int Confidence interval to be used in the plot (default = 0.95 for 95\% confidence intervals).
#' @param plot.title A title to be used in the plot.
#' @return A ggplot2 object.
#'
#' @examples
#' data("VerbAgg")
#' mod <- eirm(formula = "r2 ~ -1 + situ + btype + mode + (1|id)", data = VerbAgg)
#'
#' # Only one predictor
#' p1 <- marginalplot(mod, predictors = c("situ"))
#'
#' # Two predictors
#' p2 <- marginalplot(mod, predictors = c("situ", "btype"))
#'
#' # All three predictors
#' p3 <- marginalplot(mod, predictors = c("situ", "btype", "mode"))
#' @export

marginalplot <- function(x, predictors, conf.int = 0.95, plot.title = NULL) {

  mod <- ggpredict(x$model, terms = predictors, ci.lvl = conf.int)

  if(length(predictors) == 1) {
    # Only one variable
    p <- ggplot(mod,
                aes_string(x="x", y="predicted")) +
      geom_point(position = position_dodge(width = 0.25), size = 2) +
      geom_errorbar(aes_string(ymin="conf.low", ymax="conf.high"),
                    width=.1, position = position_dodge(width = 0.25)) +
      scale_y_continuous(limits = c(0, 1)) +
      labs(title = plot.title,
           x = get_x_title(mod),
           y = "Predicted Probability",
           shape = get_legend_title(mod)) +
      theme_bw()
  } else if(length(predictors) == 2) {
    # Two variables
    p <- ggplot(mod,
                aes_string(x="x", y="predicted", group="group")) +
      geom_point(aes_string(shape = "group"),
                 position = position_dodge(width = 0.25), size = 2) +
      geom_errorbar(aes_string(ymin="conf.low", ymax="conf.high"),
                    width=.1, position = position_dodge(width = 0.25)) +
      scale_y_continuous(limits = c(0, 1)) +
      labs(title = plot.title,
           x = get_x_title(mod),
           y = "Predicted Probability",
           shape = get_legend_title(mod)) +
      theme_bw()
  } else {
    # Three variables
    p <- ggplot(mod,
                aes_string(x="x", y="predicted", group="group")) +
      geom_point(aes_string(shape = "group"),
                 position = position_dodge(width = 0.25), size = 2) +
      geom_errorbar(aes_string(ymin="conf.low", ymax="conf.high"),
                    width=.1, position = position_dodge(width = 0.25)) +
      scale_y_continuous(limits = c(0, 1)) +
      labs(title = plot.title,
           x = get_x_title(mod),
           y = "Predicted Probability",
           shape = get_legend_title(mod)) +
      facet_wrap(~facet, scales = "free_x") +
      theme_bw()
  }

  return(p)
}