Partial dependence on one or two predictors

Averages the fitted surface over the sample at each value of one or two predictors, so that what is left is how the prediction moves with them. partial_dependence() returns the values and plot() draws them.

Usage

partial_dependence(
  object,
  variables,
  newdata = NULL,
  grid = 26L,
  values = NULL,
  level = 0.95,
  interval = "eti",
  type = "response",
  ...
)

# S3 method for class 'bartisan_partial'
print(x, digits = 3L, n_print = 10L, ...)

# S3 method for class 'bartisan_partial'
plot(x, ...)

# S3 method for class 'bartisan_fit'
plot(x, y, ...)

Arguments

object: a <bartisan_fit> object; the output of a call to bartisan().
variables: a one-sided formula naming one or two predictors, as in ~ age or ~ age + sex, or a character vector of their names.
newdata: optional; a data frame to average over. Default is the data the model was fit to.
grid: integer; how many values of a numeric predictor to evaluate. Default is 26. A factor is evaluated at each of its levels whatever this is, and a numeric predictor named second at three values rather than grid of them; see Details.
values: optional; a named list giving the values to evaluate a predictor at, which overrides grid for the predictors it names. An entry may be a function of the predictor rather than the values themselves, so that values = list(age = unique) evaluates age at every value it takes; NA is dropped from what such a function returns, while a vector written out by hand is used exactly as written. A numeric second predictor this does not name is held at three of its values near its quartiles, and a message reports which.
level: numeric; the level of the credible interval. Default is .95.
interval: string; "eti" (the default) for an equal-tailed interval from the quantiles of the draws, or "hpdi" for the highest posterior density interval, which is the shortest interval containing level of the posterior mass.
type: string; the prediction scale, passed to predict.bartisan_fit(). Default is "response".
...: for plot.bartisan_fit(), further arguments passed to partial_dependence(); for partial_dependence(), further arguments passed to predict.bartisan_fit().
x: a <bartisan_partial> object; the output of a call to partial_dependence(). For plot.bartisan_fit(), a <bartisan_fit>.
digits: integer; for print(), the number of significant digits to print the estimates and their interval to. Default is 3.
n_print: integer; for print(), the total number of rows to show, taken half from the top of the grid and half from the bottom, with the odd row going to the top. A line between the two halves counts what was left out. Default is 10, and Inf shows every row.
y: for plot.bartisan_fit(), the predictors to plot, as variables above.

Value

A <bartisan_partial> object, a data frame with one row per grid point, columns for the predictors, and estimate, lower and upper. plot() returns a ggplot2 object.

Details

At each grid value every unit is assigned that value, the prediction is taken for all of them, and the average over units is taken within each posterior draw. The interval is then computed from those averages, so it is an interval on the average prediction and not on any one unit's.

The second predictor groups the curves rather than adding an axis, which is readable for a factor and for a numeric predictor with a few values, and not for a continuous one: grid values of it would give grid curves, each with a ribbon of its own. A numeric second predictor with more than three distinct values is therefore held at three of them and a message says which, with the values argument there to choose others and values = list(z = unique) the short way to ask for all of them, which usually suits a predictor with four or five values. The three are the values nearest its quartiles, and they are distinct even when the quartiles are not, since a predictor with a large mass at one value takes that value for two or three of them; each quartile in turn takes the nearest value the predictor has that an earlier one did not take. They are values the predictor takes rather than points on an even grid, which also keeps the legend readable.

The usual caveat on a partial dependence plot applies. Averaging over the other predictors evaluates the model at covariate combinations that may not occur, so a curve over a region where the predictor's values are sparse says more about the prior than about the data, and the result summarizes the fitted function rather than supporting a causal claim; estimate_effect() is for a contrast that is meant causally.

Examples

data("rhc")
set.seed(123)

# Note: we use small sampler settings that would not be
# appropriate for analysis; this is just for the
# examples
fit <- bartisan(death ~ age + sex + meanbp + aps,
                data = rhc, num_trees = 10, num_burn = 50,
                num_draws = 50, verbose = FALSE)
#> ℹ Using `family = binomial()`.
#> ℹ Set `family` explicitly to silence this message.

# How the fitted risk moves with mean blood pressure
pd <- partial_dependence(fit, ~ meanbp)
pd
#> Partial dependence
#> 
#> Predictor: "meanbp"
#> Averaged over 1500 units, on the "response" scale
#> 
#>  meanbp estimate lower upper
#>    0.00    0.657 0.615 0.709
#>    8.88    0.657 0.615 0.709
#>   17.76    0.657 0.615 0.709
#>   26.64    0.657 0.615 0.708
#>   35.52    0.657 0.615 0.707
#>   --- 16 rows omitted ---
#>  186.48    0.637 0.544 0.688
#>  195.36    0.637 0.543 0.688
#>  204.24    0.637 0.542 0.688
#>  213.12    0.636 0.542 0.688
#>  222.00    0.636 0.542 0.688
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ `n_print` in `print()` (`?bartisan::print.bartisan_partial()`) sets how many
#>   rows are shown, half from each end; `print(., n_print = Inf)` shows all of
#>   them.

plot(pd)


# The same thing from the fit, using the `plot()` method
plot(fit, ~ meanbp)


# \donttest{
# Two predictors, one of them a factor, which gives a
# curve per level
plot(fit, ~ meanbp + sex)


# Two numeric predictors, where the second is held at
# three values and a message says which
plot(fit, ~ meanbp + aps)
#> ℹ Grouping by `aps` at 41, 54, and 68, three of its values near its quartiles.
#> ℹ Set `values` to choose them yourself.


# An entry of `values` may be a function of the
# predictor, which is how to ask for values of your own
# without naming them
three_values <- function(x) quantile(x, c(.1, .5, .9))
plot(fit, ~ meanbp + aps, values = list(aps = three_values))

# }

Usage

Arguments

Value

Details

See also

Examples