R/Interplot_mlmmi.R

In interplot: Plot the Effects of Variables in Interaction Terms

if(getRversion() >= "2.15.1") utils::globalVariables(c(".", "X.weights."))

#' Plot Conditional Coefficients in Mixed-Effects Models with Imputed Data and Interaction Terms
#'
#' \code{interplot.mlmmi} is a method to calculate conditional coefficient estimates from the results of multilevel (mixed-effects) regression models with interaction terms and multiply imputed data.
#'
#' @param m A model object including an interaction term, or, alternately, a data frame recording conditional coefficients.
#' @param var1 The name (as a string) of the variable of interest in the interaction term; its conditional coefficient estimates will be plotted.
#' @param var2 The name (as a string) of the other variable in the interaction term.
#' @param plot A logical value indicating whether the output is a plot or a dataframe including the conditional coefficient estimates of var1, their upper and lower bounds, and the corresponding values of var2.
#' @param steps Desired length of the sequence. A non-negative number, which for seq and seq.int will be rounded up if fractional. The default is 100 or the unique categories in the \code{var2} (when it is less than 100. Also see \code{\link{unique}}).
#' @param ci A numeric value defining the confidence intervals. The default value is 95\% (0.95).
#' @param adjCI Not working for `lmer` outputs yet.
#' @param hist A logical value indicating if there is a histogram of `var2` added at the bottom of the conditional effect plot.
#' @param var2_dt A numerical value indicating the frequency distribution of `var2`. It is only used when `hist == TRUE`. When the object is a model, the default is the distribution of `var2` of the model.
#' @param predPro A logical value with default of `FALSE`. When the `m` is an object of class `glmerMod` and the argument is set to `TRUE`, the function will plot predicted probabilities at the values given by `var2_vals`.
#' @param var2_vals A numerical value indicating the values the predicted probabilities are estimated, when `predPro` is `TRUE`.
#' @param point A logical value determining the format of plot. By default, the function produces a line plot when var2 takes on ten or more distinct values and a point (dot-and-whisker) plot otherwise; option TRUE forces a point plot.
#' @param sims Number of independent simulation draws used to calculate upper and lower bounds of coefficient estimates: lower values run faster; higher values produce smoother curves.
#' @param xmin A numerical value indicating the minimum value shown of x shown in the graph. Rarely used.
#' @param xmax A numerical value indicating the maximum value shown of x shown in the graph. Rarely used.
#' @param ercolor A character value indicating the outline color of the whisker or ribbon.
#' @param esize A numerical value indicating the size of the whisker or ribbon.
#' @param ralpha A numerical value indicating the transparency of the ribbon.
#' @param rfill A character value indicating the filling color of the ribbon.
#' @param stats_cp A character value indicating what statistics to present as the plot note. Three options are available: "none", "ci", and "ks". The default is "none". See the Details for more information.
#' @param txt_caption A character string to add a note for the plot, a value will sending to \code{ggplot2::labs(caption = txt_caption))}.
#' @param facet_labs An optional character vector of facet labels to be used when plotting an interaction with a factor variable.
#' @param var3 An optional name (as a string) of a third variable for a three-way interaction \code{var1 * var2 * var3}, pooled across imputations. When supplied, the conditional effect of \code{var1} across \code{var2} is shown at several values/levels of \code{var3}. The default \code{NULL} gives the standard two-way behavior. Requires continuous \code{var1} and \code{var2}.
#' @param var3_vals An optional numeric vector giving the values of a continuous \code{var3} to condition on. The default is the mean and the mean plus or minus one standard deviation. Ignored when \code{var3} is a factor.
#' @param facet A logical value, used only with \code{var3}. \code{TRUE} (default) draws one panel per value/level of \code{var3}; \code{FALSE} overlays the curves colored by \code{var3}.
#' @param ... Other ggplot aesthetics arguments for points in the dot-whisker plot or lines in the line-ribbon plots. Not currently used.
#'
#' @details \code{interplot.mlmmi} and \code{interplot.gmlmmi} are S3 methods from the \code{interplot}. It works on lists of mixed-effects objects with class \code{lmerMod} and \code{glmerMod} generated by \code{mitools} and \code{lme4}.
#'
#' Because the output function is based on \code{\link[ggplot2]{ggplot}}, any additional arguments and layers supported by \code{ggplot2} can be added with the \code{+}.
#'
#' \code{interplot} visualizes the conditional effect based on simulated marginal effects. The simulation provides a probabilistic distribution of moderation effect of the conditioning variable (\code{var2}) at every preset values (including the minimum and maximum values) of the conditioned variable (\code{var1}), denoted as Emin and Emax. This output allows the function to further examine the conditional effect statistically in two ways. One is to examine if the distribution of \eqn{Emax - Emin} covers zero. The other is to directly compare Emin and Emax through statistical tools for distributional comparisons. Users can choose either method by setting the argument \code{stats_cp} to "ci" or "ks".
#' \itemize{
#'   \item "ci" provides the confidence interval of the difference of \eqn{Emax - Emin}. An interval including 0 suggests no statistical difference before and after the conditional effect is applied, and vise versa.
#'   \item "ks" presents the result of a two-sample Kolmogorov-Smirnov test of the simulated distributions of Emin and Emax. The output includes a D statistics and a p-value of the null hypothesis that the two distributions come from the same distribution at the 0.05 level.
#' }
#'
#' See an illustration in the package vignette.
#'
#'
#' @return The function returns a \code{ggplot} object.
#'
#' @importFrom abind abind
#' @importFrom arm sim
#' @importFrom stats quantile median plogis model.matrix ks.test
#' @importFrom purrr map map2 list_c list_rbind
#' @import ggplot2
#' @import dplyr
#'
#'
#' @export

# Coding function for mlm, mi objects
interplot.mlmmi <- function(m, var1, var2, plot = TRUE, steps = NULL,
                            ci = .95, adjCI = FALSE, hist = FALSE,
                            var2_dt = NA, predPro = FALSE, var2_vals = NULL,
                            point = FALSE, sims = 5000, xmin = NA, xmax = NA,
                            ercolor = NA, esize = 0.5, ralpha = 0.5,
                            rfill = "grey70", stats_cp = "none",
                            txt_caption = NULL, facet_labs = NULL,
                            var3 = NULL, var3_vals = NULL, facet = TRUE, ...) {

  if (predPro == TRUE)
    stop("Predicted probability is estimated only for general linear models.")

  # Combine simulations from multiple imputed models
  m.list <- m
  m <- m.list[[1]]
  class(m.list) <- class(m)
  m.sims.list <- map(m.list, \(i) arm::sim(i, sims))
  m.sims <- m.sims.list[[1]]

  for (i in 2:length(m.sims.list)) {
    m.sims@fixef <- rbind(m.sims@fixef, m.sims.list[[i]]@fixef)
    m.sims@ranef[[1]] <- abind::abind(m.sims@ranef[[1]], m.sims.list[[i]]@ranef[[1]],
                                       along = 1)
  }

  coef_names <- unlist(dimnames(m@pp$X)[2])

  ## For factor base terms ####
  orig_var1 <- var1
  orig_var2 <- var2

  factor_v1 <- is.factor(m@frame[[var1]])
  factor_v2 <- is.factor(m@frame[[var2]])

  if (factor_v1 & factor_v2)
    stop("The function does not support interactions between two factors.")

  # Three-way and nonlinear paths (pooled across imputations) ####
  fixef_draws <- m.sims@fixef
  colnames(fixef_draws) <- coef_names
  out <- interplot_mi_extra(
    m = m, m.list = m.list, draws_coef = fixef_draws, frame = m@frame,
    tt_fixed = stats::terms(m, fixed.only = TRUE), coef_fun = lme4::fixef,
    vcov_fun = stats::vcov, df_resid = NULL, var1 = var1, var2 = var2,
    var3 = var3, var3_vals = var3_vals, facet = facet, factor_v1 = factor_v1,
    factor_v2 = factor_v2, predPro = predPro, plot = plot, steps = steps,
    ci = ci, adjCI = adjCI, hist = hist, var2_dt = var2_dt, point = point,
    sims = sims, xmin = xmin, xmax = xmax, ercolor = ercolor, esize = esize,
    ralpha = ralpha, rfill = rfill, stats_cp = stats_cp,
    txt_caption = txt_caption, orig_var1 = orig_var1, orig_var2 = orig_var2, ...
  )
  if (!is.null(out)) return(out$value)

  if (factor_v1) {
    var1 <- paste0(var1, levels(m@frame[[var1]]))
  } else if (factor_v2) {
    var2 <- paste0(var2, levels(m@frame[[var2]]))
  }

  # Build interaction term names
  if (factor_v1 | factor_v2) {
    var12 <- paste0(var2, ":", var1)[-1]
    for (i in seq(var12)) {
      if (!var12[i] %in% coef_names)
        var12[i] <- paste0(var1, ":", var2)[-1][i]
      if (!var12[i] %in% coef_names)
        stop(paste("Model does not include the interaction of", var1, "and", var2, "."))
    }
  } else {
    if (var1 == var2) {
      var12 <- paste0("I(", var1, "^2)")
    } else {
      var12 <- paste0(var2, ":", var1)
    }
    if (!var12 %in% coef_names)
      var12 <- paste0(var1, ":", var2)
    if (!var12 %in% coef_names)
      stop(paste("Model does not include the interaction of", var1, "and", var2, "."))
  }

  # Set the min, max, and steps ####
  if (factor_v2) {
    xmin <- 0
    xmax <- 1
    steps <- 2
  } else {
    if (is.na(xmin)) xmin <- min(m@frame[[var2]], na.rm = TRUE)
    if (is.na(xmax)) xmax <- max(m@frame[[var2]], na.rm = TRUE)
    if (is.null(steps)) {
      steps <- length(unique(na.omit(m@frame[[var2]])))
    }
    if (steps > 100) steps <- 100
  }

  fake_seq <- seq(xmin, xmax, length.out = steps)
  ci_diff <- NULL
  ks_diff <- NULL

  # Calculate the effects ####
  if (factor_v1) {
    results <- map(seq(var1)[-length(var1)], \(j) {
      sim_v1 <- m.sims@fixef[, match(var1[j + 1], coef_names)]
      sim_v12 <- m.sims@fixef[, match(var12[j], coef_names)]
      coef <- sim_coef_vec(sim_v1, sim_v12, fake_seq, ci)
      d <- sim_diff_stats(sim_v1, sim_v12, xmin, xmax, ci)
      coef$value <- var1[j + 1]
      list(coef = coef, ci_diff = d$ci_diff)
    })

    coef_df <- map(results, "coef") |> list_rbind()
    ci_diff <- map(results, "ci_diff") |> setNames(var12)

    if (plot == TRUE) {
      if (hist == TRUE & all(is.na(var2_dt))) var2_dt <- m@frame[[var2]]
      if (is.null(facet_labs)) facet_labs <- unique(coef_df$value)

      if (stats_cp == "ci") {
        facet_labs <- map2(ci_diff, facet_labs, \(aCI, aLevel) {
          paste0(aLevel, "\n CI(Max - Min): [",
                 round(aCI[1], digits = 3), ", ",
                 round(aCI[2], digits = 3), "]")
        }) |> list_c()
      }

      coef_df$value <- factor(coef_df$value, labels = facet_labs)
      interplot.plot(m = coef_df, var1 = orig_var1, var2 = orig_var2,
                     hist = hist, steps = steps, var2_dt = var2_dt,
                     point = point, ercolor = ercolor, esize = esize,
                     ralpha = ralpha, rfill = rfill, ci_diff = ci_diff,
                     ks_diff = ks_diff, stats_cp = stats_cp,
                     txt_caption = txt_caption, ...) +
        facet_grid(. ~ value)
    } else {
      names(coef_df)[1:4] <- c(var2, "coef", "ub", "lb")
      return(list(df_coef = coef_df, stats_ci = ci_diff))
    }

  } else if (factor_v2) {
    results <- map(seq(var2)[-length(var2)], \(j) {
      sim_v1 <- m.sims@fixef[, match(var1, coef_names)]
      sim_v12 <- m.sims@fixef[, match(var12[j], coef_names)]
      coef <- sim_coef_vec(sim_v1, sim_v12, fake_seq, ci)
      d <- sim_diff_stats(sim_v1, sim_v12, xmin, xmax, ci)
      coef$value <- var2[j + 1]
      list(coef = coef, ci_diff = d$ci_diff)
    })

    coef_df <- map(results, "coef") |> list_rbind()
    ci_diff <- map(results, "ci_diff") |> setNames(var12)

    if (plot == TRUE) {
      if (hist == TRUE & all(is.na(var2_dt))) var2_dt <- m@frame[[var2]]
      if (is.null(facet_labs)) facet_labs <- unique(coef_df$value)

      if (stats_cp == "ci") {
        facet_labs <- map2(ci_diff, facet_labs, \(aCI, aLevel) {
          paste0(aLevel, "\n CI(Max - Min): [",
                 round(aCI[1], digits = 3), ", ",
                 round(aCI[2], digits = 3), "]")
        }) |> list_c()
      }

      coef_df$value <- factor(coef_df$value, labels = facet_labs)
      interplot.plot(m = coef_df, var1 = orig_var1, var2 = orig_var2,
                     steps = steps, hist = hist, var2_dt = var2_dt,
                     point = point, ercolor = ercolor, esize = esize,
                     ralpha = ralpha, rfill = rfill, ci_diff = ci_diff,
                     ks_diff = ks_diff, stats_cp = stats_cp,
                     txt_caption = txt_caption, ...) +
        facet_grid(. ~ value)
    } else {
      names(coef_df)[1:4] <- c(var2, "coef", "ub", "lb")
      return(list(df_coef = coef_df, stats_ci = ci_diff))
    }

  } else {
    ## Numeric case
    multiplier <- if (var1 == var2) 2 else 1

    sim_v1 <- m.sims@fixef[, match(var1, coef_names)]
    sim_v12 <- m.sims@fixef[, match(var12, coef_names)]

    coef <- sim_coef_vec(sim_v1, sim_v12, fake_seq, ci, multiplier)

    diff_stats <- sim_diff_stats(sim_v1, sim_v12, xmin, xmax, ci, multiplier)
    ci_diff <- diff_stats$ci_diff
    ks_diff <- diff_stats$ks_diff

    if (plot == TRUE) {
      if (hist == TRUE & all(is.na(var2_dt))) var2_dt <- m@frame[[var2]]
      interplot.plot(m = coef, var1 = orig_var1, var2 = orig_var2,
                     hist = hist, var2_dt = var2_dt, point = point,
                     ercolor = ercolor, esize = esize, ralpha = ralpha,
                     rfill = rfill, ci_diff = ci_diff, ks_diff = ks_diff,
                     stats_cp = stats_cp, txt_caption = txt_caption, ...)
    } else {
      names(coef) <- c(var2, "coef", "ub", "lb")
      return(list(df_coef = coef, stats_ci = ci_diff))
    }
  }
}


#' @export
interplot.gmlmmi <- function(m, var1, var2, plot = TRUE, steps = NULL,
                             ci = .95, adjCI = FALSE, hist = FALSE,
                             var2_dt = NA, predPro = FALSE, var2_vals = NULL,
                             point = FALSE, sims = 5000, xmin = NA, xmax = NA,
                             ercolor = NA, esize = 0.5, ralpha = 0.5,
                             rfill = "grey70", stats_cp = "none",
                             txt_caption = NULL, facet_labs = NULL,
                             var3 = NULL, var3_vals = NULL, facet = TRUE, ...) {

  # Combine simulations from multiple imputed models
  m.list <- m
  m <- m.list[[1]]
  class(m.list) <- class(m)
  m.sims.list <- map(m.list, \(i) arm::sim(i, sims))
  m.sims <- m.sims.list[[1]]

  for (i in 2:length(m.sims.list)) {
    m.sims@fixef <- rbind(m.sims@fixef, m.sims.list[[i]]@fixef)
    m.sims@ranef[[1]] <- abind::abind(m.sims@ranef[[1]], m.sims.list[[i]]@ranef[[1]],
                                       along = 1)
  }

  coef_names <- unlist(dimnames(m@pp$X)[2])

  ## For factor base terms ####
  orig_var1 <- var1
  orig_var2 <- var2

  factor_v1 <- is.factor(m@frame[[var1]])
  factor_v2 <- is.factor(m@frame[[var2]])

  if (factor_v1 & factor_v2)
    stop("The function does not support interactions between two factors.")
  if ((factor_v1 | factor_v2) & predPro == TRUE)
    stop("The current version does not support estimating predicted probabilities for factor base terms.")

  # Three-way and nonlinear paths (pooled across imputations) ####
  fixef_draws <- m.sims@fixef
  colnames(fixef_draws) <- coef_names
  out <- interplot_mi_extra(
    m = m, m.list = m.list, draws_coef = fixef_draws, frame = m@frame,
    tt_fixed = stats::terms(m, fixed.only = TRUE), coef_fun = lme4::fixef,
    vcov_fun = stats::vcov, df_resid = NULL, var1 = var1, var2 = var2,
    var3 = var3, var3_vals = var3_vals, facet = facet, factor_v1 = factor_v1,
    factor_v2 = factor_v2, predPro = predPro, plot = plot, steps = steps,
    ci = ci, adjCI = adjCI, hist = hist, var2_dt = var2_dt, point = point,
    sims = sims, xmin = xmin, xmax = xmax, ercolor = ercolor, esize = esize,
    ralpha = ralpha, rfill = rfill, stats_cp = stats_cp,
    txt_caption = txt_caption, orig_var1 = orig_var1, orig_var2 = orig_var2, ...
  )
  if (!is.null(out)) return(out$value)

  if (factor_v1) {
    var1 <- paste0(var1, levels(m@frame[[var1]]))
  } else if (factor_v2) {
    var2 <- paste0(var2, levels(m@frame[[var2]]))
  }

  # Build interaction term names
  if (factor_v1 | factor_v2) {
    var12 <- paste0(var2, ":", var1)[-1]
    for (i in seq(var12)) {
      if (!var12[i] %in% coef_names)
        var12[i] <- paste0(var1, ":", var2)[-1][i]
      if (!var12[i] %in% coef_names)
        stop(paste("Model does not include the interaction of", var1, "and", var2, "."))
    }
  } else {
    if (var1 == var2) {
      var12 <- paste0("I(", var1, "^2)")
    } else {
      var12 <- paste0(var2, ":", var1)
    }
    if (!var12 %in% coef_names)
      var12 <- paste0(var1, ":", var2)
    if (!var12 %in% coef_names)
      stop(paste("Model does not include the interaction of", var1, "and", var2, "."))
  }

  # Set the min, max, and steps ####
  if (factor_v2) {
    xmin <- 0
    xmax <- 1
    steps <- 2
  } else {
    if (is.na(xmin)) xmin <- min(m@frame[[var2]], na.rm = TRUE)
    if (is.na(xmax)) xmax <- max(m@frame[[var2]], na.rm = TRUE)
    if (is.null(steps)) {
      steps <- length(unique(na.omit(m@frame[[var2]])))
    }
    if (steps > 100) steps <- 100
  }

  fake_seq <- seq(xmin, xmax, length.out = steps)
  ci_diff <- NULL
  ks_diff <- NULL

  # Calculate the effects ####
  if (factor_v1) {
    results <- map(seq(var1)[-length(var1)], \(j) {
      sim_v1 <- m.sims@fixef[, match(var1[j + 1], coef_names)]
      sim_v12 <- m.sims@fixef[, match(var12[j], coef_names)]
      coef <- sim_coef_vec(sim_v1, sim_v12, fake_seq, ci)
      d <- sim_diff_stats(sim_v1, sim_v12, xmin, xmax, ci)
      coef$value <- var1[j + 1]
      list(coef = coef, ci_diff = d$ci_diff)
    })

    coef_df <- map(results, "coef") |> list_rbind()
    ci_diff <- map(results, "ci_diff") |> setNames(var12)

    if (plot == TRUE) {
      if (hist == TRUE & all(is.na(var2_dt))) var2_dt <- m@frame[[var2]]
      if (is.null(facet_labs)) facet_labs <- unique(coef_df$value)

      if (stats_cp == "ci") {
        facet_labs <- map2(ci_diff, facet_labs, \(aCI, aLevel) {
          paste0(aLevel, "\n CI(Max - Min): [",
                 round(aCI[1], digits = 3), ", ",
                 round(aCI[2], digits = 3), "]")
        }) |> list_c()
      }

      coef_df$value <- factor(coef_df$value, labels = facet_labs)
      interplot.plot(m = coef_df, var1 = orig_var1, var2 = orig_var2,
                     steps = steps, hist = hist, var2_dt = var2_dt,
                     point = point, ercolor = ercolor, esize = esize,
                     ralpha = ralpha, rfill = rfill, ci_diff = ci_diff,
                     ks_diff = ks_diff, stats_cp = stats_cp,
                     txt_caption = txt_caption, ...) +
        facet_grid(. ~ value)
    } else {
      names(coef_df)[1:4] <- c(var2, "coef", "ub", "lb")
      return(list(df_coef = coef_df, stats_ci = ci_diff))
    }

  } else if (factor_v2) {
    results <- map(seq(var2)[-length(var2)], \(j) {
      sim_v1 <- m.sims@fixef[, match(var1, coef_names)]
      sim_v12 <- m.sims@fixef[, match(var12[j], coef_names)]
      coef <- sim_coef_vec(sim_v1, sim_v12, fake_seq, ci)
      d <- sim_diff_stats(sim_v1, sim_v12, xmin, xmax, ci)
      coef$value <- var2[j + 1]
      list(coef = coef, ci_diff = d$ci_diff)
    })

    coef_df <- map(results, "coef") |> list_rbind()
    ci_diff <- map(results, "ci_diff") |> setNames(var12)

    if (plot == TRUE) {
      if (hist == TRUE & all(is.na(var2_dt))) var2_dt <- m@frame[[var2]]
      if (is.null(facet_labs)) facet_labs <- unique(coef_df$value)

      if (stats_cp == "ci") {
        facet_labs <- map2(ci_diff, facet_labs, \(aCI, aLevel) {
          paste0(aLevel, "\n CI(Max - Min): [",
                 round(aCI[1], digits = 3), ", ",
                 round(aCI[2], digits = 3), "]")
        }) |> list_c()
      }

      coef_df$value <- factor(coef_df$value, labels = facet_labs)
      interplot.plot(m = coef_df, var1 = orig_var1, var2 = orig_var2,
                     steps = steps, hist = hist, var2_dt = var2_dt,
                     point = point, ercolor = ercolor, esize = esize,
                     ralpha = ralpha, rfill = rfill, ci_diff = ci_diff,
                     ks_diff = ks_diff, stats_cp = stats_cp,
                     txt_caption = txt_caption, ...) +
        facet_grid(. ~ value)
    } else {
      names(coef_df)[1:4] <- c(var2, "coef", "ub", "lb")
      return(list(df_coef = coef_df, stats_ci = ci_diff))
    }

  } else {
    ## Numeric case
    multiplier <- if (var1 == var2) 2 else 1

    if (predPro == TRUE) {
      if (is.null(var2_vals))
        stop("The predicted probabilities cannot be estimated without defining 'var2_vals'.")

      df <- data.frame(m@frame)
      df[[names(m@flist)]] <- NULL
      if (any(grepl("X.weights.", names(df))))
        df <- select(df, -X.weights.)
      df_temp <- select(df, 1)
      df <- df[-1] |>
        map(\(v) {
          if (is.factor(v)) {
            as.data.frame(model.matrix(~ v - 1)[, -1, drop = FALSE])
          } else {
            as.numeric(v)
          }
        })

      for (i in seq(df)) {
        if (!is.data.frame(df[[i]])) {
          namesUpdate <- c(names(df_temp), names(df)[[i]])
          df_temp <- cbind(df_temp, df[[i]])
          names(df_temp) <- namesUpdate
        } else {
          df_temp <- cbind(df_temp, df[[i]])
        }
      }
      df <- df_temp

      names(df)[1] <- "(Intercept)"
      df$`(Intercept)` <- 1

      if (var1 == var2) {
        names(df) <- sub("I\\.(.*)\\.2\\.", "I\\(\\1\\^2\\)", names(df))
      }

      iv_medians <- summarize(df, across(everything(), \(x) median(x, na.rm = TRUE)))
      fake_data <- iv_medians[rep(1, each = steps * length(var2_vals)), ]
      fake_data[[var1]] <- rep(
        seq(min(df[[var1]]), max(df[[var1]]), length.out = steps),
        times = length(var2_vals)
      )
      fake_data[[var2]] <- rep(var2_vals, each = steps)
      fake_data[[var12]] <- fake_data[[var1]] * fake_data[[var2]]

      linpred <- data.matrix(fake_data) %*% t(data.matrix(m.sims@fixef))
      pp_mat <- plogis(linpred)
      probs <- c((1 - ci) / 2, 1 - (1 - ci) / 2)

      coef <- data.frame(
        fake = fake_data[[var1]],
        value = as.factor(fake_data[[var2]]),
        coef1 = rowMeans(pp_mat) * 100,
        lb = apply(pp_mat, 1, quantile, probs = probs[1]) * 100,
        ub = apply(pp_mat, 1, quantile, probs = probs[2]) * 100
      )
    } else {
      sim_v1 <- m.sims@fixef[, match(var1, coef_names)]
      sim_v12 <- m.sims@fixef[, match(var12, coef_names)]

      coef <- sim_coef_vec(sim_v1, sim_v12, fake_seq, ci, multiplier)

      diff_stats <- sim_diff_stats(sim_v1, sim_v12, xmin, xmax, ci, multiplier)
      ci_diff <- diff_stats$ci_diff
      ks_diff <- diff_stats$ks_diff
    }

    if (plot == TRUE) {
      if (hist == TRUE & all(is.na(var2_dt))) var2_dt <- m@frame[[var2]]
      interplot.plot(m = coef, var1 = orig_var1, var2 = orig_var2,
                     steps = steps, hist = hist, predPro = predPro,
                     var2_vals = var2_vals, var2_dt = var2_dt, point = point,
                     ercolor = ercolor, esize = esize, ralpha = ralpha,
                     rfill = rfill, ci_diff = ci_diff, ks_diff = ks_diff,
                     stats_cp = stats_cp, txt_caption = txt_caption, ...)
    } else {
      if (predPro == TRUE) {
        names(coef) <- c(var2, paste0("values_in_", var1), "coef", "ub", "lb")
      } else {
        names(coef) <- c(var2, "coef", "ub", "lb")
      }
      return(list(df_coef = coef, stats_ci = ci_diff))
    }
  }
}