Get marginal effects from model terms

Description

ggpredict() computes predicted (fitted) values for the response, at the margin of specific values from certain model terms, where additional model terms indicate the grouping structure. ggaverage() computes the average predicted values. The result is returned as tidy data frame.

Usage

ggaverage(model, terms, ci.lvl = 0.95, type = c("fe", "re"),
  typical = "mean", ppd = FALSE, x.as.factor = FALSE, pretty = TRUE,
  condition = NULL, ...)

ggpredict(model, terms, ci.lvl = 0.95, type = c("fe", "re"),
  full.data = FALSE, typical = "mean", ppd = FALSE, x.as.factor = FALSE,
  pretty = TRUE, condition = NULL, ...)

Arguments

Details

Supported Models

Currently supported model-objects are: lm, glm, glm.nb, lme, lmer, glmer, glmer.nb, nlmer, glmmTMB, gam, vgam, gamm, gamm4, multinom, betareg, gls, gee, plm, lrm, polr, clm, hurdle, zeroinfl, svyglm, svyglm.nb, truncreg, coxph, stanreg, brmsfit, lmRob, glmRob, brglm, rlm. Other models not listed here are passed to a generic predict-function and might work as well, or maybe with ggeffect(), which effectively does the same as ggpredict(). The main difference between ggpredict() and ggeffect() is how factors are held constant: ggpredict() uses the reference level, while ggeffect() computes a kind of "average" value, which represents the proportions of each factor's category.

Marginal Effects at Specific Values

Specific values of model terms can be specified via the terms-argument. Indicating levels in square brackets allows for selecting only specific groups or values resp. value ranges. Term name and levels in brackets must be separated by a whitespace character, e.g. terms = c("age", "education [1,3]"). Numeric ranges, separated with colon, are also allowed: terms = c("education", "age [30:60]").

The terms-argument also supports the same shortcuts as the mdrt.values-argument in gginteraction(). So terms = "age [meansd]" would return predictions for the values one standard deviation below the mean age, the mean age and one SD above the mean age. terms = "age [quart2]" would calculate predictions at the value of the lower, median and upper quartile of age.

Furthermore, it is possible to specify a function name. Values for predictions will then be transformed, e.g. terms = "income [exp]". This is useful when model predictors were transformed for fitting the model and should be back-transformed to the original scale for predictions.

A last option for numeric terms is pretty, e.g. terms = "age [pretty]". In this case, values are based on a "pretty" range of the variable, which also means you can selectively prettify certain terms (while the logical argument pretty, see above, usually prettifies all variables with more than 25 unique values). See package vignettes.

Holding covariates at constant values

For ggpredict(), if full.data = FALSE, expand.grid() is called on all unique combinations of model.frame(model)[, terms] and used as newdata-argument for predict(). In this case, all remaining covariates that are not specified in terms are held constant. Numeric values are set to the mean (unless changed with the condition or typical-argument), factors are set to their reference level (may also be changed with condition) and character vectors to their mode (most common element).

ggaverage() computes the average predicted values, by calling ggpredict() with full.data = TRUE, where argument newdata = model.frame(model) is used in predict(). Hence, predictions are made on the model data. In this case, all remaining covariates that are not specified in terms are not held constant, but vary between observations (and are kept as they happen to be). The predicted values are then averaged for each group (if any).

Thus, ggpredict() can be considered as calculating marginal effects at the mean, while ggaverage() computes average marginal effects.

Bayesian Regression Models

ggpredict() also works with Stan-models from the rstanarm or brms-package. The predicted values are the median value of all drawn posterior samples. The confidence intervals for Stan-models are actually high density intervals, computed by hdi, unless ppd = TRUE. If ppd = TRUE, predictions are based on draws of the posterior predictive distribution and the uncertainty interval is computed using predictive_interval. By default (i.e. ppd = FALSE), the predictions are based on posterior_linpred and hence have some limitations: the uncertainty of the error term is not taken into account. The recommendation is to use the posterior predictive distribution (posterior_predict). Note that for binomial models, the newdata-argument used in posterior_predict() must also contain the vector with the number of trials. In this case, a dummy-vector is used, where all values for the response are set to 1.

Value

A tibble (with ggeffects class attribute) with consistent data columns:

x

the values of the first term in terms, used as x-position in plots.

predicted

the predicted values of the response, used as y-position in plots.

conf.low

the lower bound of the confidence interval for the predicted values.

conf.high

the upper bound of the confidence interval for the predicted values.

observed

if full.data = TRUE, this columns contains the observed values (the response vector).

residuals

if full.data = TRUE, this columns contains residuals.

group

the grouping level from the second term in terms, used as grouping-aesthetics in plots.

facet

the grouping level from the third term in terms, used to indicate facets in plots.

For proportional odds logistic regression (see polr) resp. cumulative link models (e.g., see clm), an additional column response.level is returned, which indicates the grouping of predictions based on the level of the model's response.

Note

Since data for ggaverage() comes from the model frame, not all possible combinations of values in terms might be present in the data, thus lines or confidence bands from plot() might not span over the complete x-axis-range.

There are some limitations for certain model objects. For example, it is currently only possible to compute predicted risk scores for coxph-models, but not expected number of events nor survival probabilities.

polr- or clm-models have an additional column response.level, which indicates with which level of the response variable the predicted values are associated.

There is a summary()-method, which gives a cleaner output (especially for predictions by groups), and which indicates at which values covariates were held constant.

Examples

data(efc)
fit <- lm(barthtot ~ c12hour + neg_c_7 + c161sex + c172code, data = efc)

ggpredict(fit, terms = "c12hour")
ggpredict(fit, terms = "c12hour", full.data = TRUE)
ggpredict(fit, terms = c("c12hour", "c172code"))
ggpredict(fit, terms = c("c12hour", "c172code", "c161sex"))

# only range of 40 to 60 for variable 'c12hour'
ggpredict(fit, terms = "c12hour [40:60]")

# using "summary()" shows that covariate "neg_c_7" is held
# constant at a value of 11.84 (its mean value). To use a
# different value, use "condition"
ggpredict(fit, terms = "c12hour [40:60]", condition = c(neg_c_7 = 20))

# to plot ggeffects-objects, you can use the 'plot()'-function.
# the following examples show how to build your ggplot by hand.

# plot predicted values, remaining covariates held constant
library(ggplot2)
mydf <- ggpredict(fit, terms = "c12hour")
ggplot(mydf, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)

# with "full.data = TRUE", remaining covariates vary between
# observations, so fitted values can be plotted
mydf <- ggpredict(fit, terms = "c12hour", full.data = TRUE)
ggplot(mydf, aes(x, predicted)) + geom_point()

# you can add a smoothing-geom to show the linear trend of fitted values
ggplot(mydf, aes(x, predicted)) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_point()

# three variables, so we can use facets and groups
mydf <- ggpredict(
  fit,
  terms = c("c12hour", "c161sex", "c172code"),
  full.data = TRUE
)
ggplot(mydf, aes(x = x, y = predicted, colour = group)) +
  stat_smooth(method = "lm", se = FALSE) +
  facet_wrap(~facet, ncol = 2)

# average marginal effects
mydf <- ggaverage(fit, terms = c("c12hour", "c172code"))
ggplot(mydf, aes(x = x, y = predicted, colour = group)) +
  stat_smooth(method = "lm", se = FALSE)

# select specific levels for grouping terms
mydf <- ggpredict(fit, terms = c("c12hour", "c172code [1,3]", "c161sex"))
ggplot(mydf, aes(x = x, y = predicted, colour = group)) +
  stat_smooth(method = "lm", se = FALSE) +
  facet_wrap(~facet) +
  labs(
    y = get_y_title(mydf),
    x = get_x_title(mydf),
    colour = get_legend_title(mydf)
  )

# level indication also works for factors with non-numeric levels
# and in combination with numeric levels for other variables
library(sjlabelled)
data(efc)
efc$c172code <- as_label(efc$c172code)
fit <- lm(barthtot ~ c12hour + neg_c_7 + c161sex + c172code, data = efc)
ggpredict(fit, terms = c("c12hour",
  "c172code [low level of education, high level of education]",
  "c161sex [1]"))

# use categorical value on x-axis, use axis-labels, add error bars
dat <- ggpredict(fit, terms = c("c172code", "c161sex"))
ggplot(dat, aes(x, predicted, colour = group)) +
  geom_point(position = position_dodge(.1)) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(.1)
  ) +
  scale_x_continuous(breaks = 1:3, labels = get_x_labels(dat))

# 3-way-interaction with 2 continuous variables
data(efc)
# make categorical
efc$c161sex <- as_factor(efc$c161sex)
fit <- lm(neg_c_7 ~ c12hour * barthtot * c161sex, data = efc)
# select only levels 30, 50 and 70 from continuous variable Barthel-Index
dat <- ggpredict(fit, terms = c("c12hour", "barthtot [30,50,70]", "c161sex"))
ggplot(dat, aes(x = x, y = predicted, colour = group)) +
  stat_smooth(method = "lm", se = FALSE, fullrange = TRUE) +
  facet_wrap(~facet) +
  labs(
    colour = get_legend_title(dat),
    x = get_x_title(dat),
    y = get_y_title(dat),
    title = get_title(dat)
  )

# or with ggeffects' plot-method
## Not run: 
plot(dat, ci = FALSE)
## End(Not run)

# use factor levels as x-column in returned data frame
data(efc)
efc$c161sex <- as_label(efc$c161sex)
fit <- lm(neg_c_7 ~ c12hour + c161sex, data = efc)
ggpredict(fit, terms = "c161sex", x.as.factor = TRUE)