R/get_predictions_merMod.R

In javifar/ggeffects: Create Tidy Data Frames of Marginal Effects for 'ggplot' from Model Outputs

get_predictions_merMod <- function(model, fitfram, ci.lvl, linv, type, terms, value_adjustment, condition, ...) {
  # does user want standard errors?
  se <- !is.null(ci.lvl) && !is.na(ci.lvl)

  # compute ci, two-ways
  if (!is.null(ci.lvl) && !is.na(ci.lvl))
    ci <- (1 + ci.lvl) / 2
  else
    ci <- .975

  # check whether predictions should be conditioned
  # on random effects (grouping level) or not.
  if (type == "fe")
    ref <- NA
  else
    ref <- NULL

  if (type %in% c("sim", "sim_re")) {

    # simulate predictions
    fitfram <- .do_simulate(model, terms, ci, type, ...)

  } else {

    fitfram$predicted <- suppressWarnings(stats::predict(
      model,
      newdata = fitfram,
      type = "response",
      re.form = ref,
      allow.new.levels = TRUE,
      ...
    ))

    if (se) {
      # get standard errors from variance-covariance matrix
      se.pred <-
        .standard_error_predictions(
          model = model,
          prediction_data = fitfram,
          value_adjustment = value_adjustment,
          terms = terms,
          type = type,
          condition = condition
        )

      if (.check_returned_se(se.pred)) {
        se.fit <- se.pred$se.fit
        fitfram <- se.pred$prediction_data

        if (is.null(linv)) {
          # calculate CI for linear mixed models
          fitfram$conf.low <- fitfram$predicted - stats::qnorm(ci) * se.fit
          fitfram$conf.high <- fitfram$predicted + stats::qnorm(ci) * se.fit
        } else {
          # get link-function and back-transform fitted values
          # to original scale, so we compute proper CI
          lf <- insight::link_function(model)

          # calculate CI for glmm
          fitfram$conf.low <- linv(lf(fitfram$predicted) - stats::qnorm(ci) * se.fit)
          fitfram$conf.high <- linv(lf(fitfram$predicted) + stats::qnorm(ci) * se.fit)
        }

        # copy standard errors
        attr(fitfram, "std.error") <- se.fit
        attr(fitfram, "prediction.interval") <- attr(se.pred, "prediction_interval")
      } else {
        fitfram$conf.low <- NA
        fitfram$conf.high <- NA
      }

    } else {
      fitfram$conf.low <- NA
      fitfram$conf.high <- NA
    }

  }

  fitfram
}