plot.vsel: Plot summary statistics of a variable selection

View source: R/methods.R

plot.vselR Documentation

Plot summary statistics of a variable selection


This is the plot() method for vsel objects (returned by varsel() or cv_varsel()).


## S3 method for class 'vsel'
  nterms_max = NULL,
  stats = "elpd",
  deltas = FALSE,
  alpha = 0.32,
  baseline = if (!inherits(x$refmodel, "datafit")) "ref" else "best",
  thres_elpd = NA,



An object of class vsel (returned by varsel() or cv_varsel()).


Maximum submodel size for which the statistics are calculated. Using NULL is effectively the same as using length(solution_terms(object)). Note that nterms_max does not count the intercept, so use nterms_max = 0 for the intercept-only model. For plot.vsel(), nterms_max must be at least 1.


One or more character strings determining which performance statistics (i.e., utilities or losses) to calculate. Available statistics are:

  • "elpd": (expected) sum of log predictive densities.

  • "mlpd": mean log predictive density, that is, "elpd" divided by the number of observations.

  • "mse": mean squared error.

  • "rmse": root mean squared error. For the corresponding standard error and lower and upper confidence interval bounds, bootstrapping is used.

  • "acc" (or its alias, "pctcorr"): classification accuracy (binomial() family only).

  • "auc": area under the ROC curve (binomial() family only). For the corresponding standard error and lower and upper confidence interval bounds, bootstrapping is used.


If TRUE, the submodel statistics are estimated as differences from the baseline model (see argument baseline) instead of estimating the actual values of the statistics.


A number determining the (nominal) coverage 1 - alpha of the normal-approximation (or bootstrap; see argument stats) confidence intervals. For example, in case of the normal approximation, alpha = 0.32 corresponds to one-standard-error intervals (because of the coverage of 68%).


For summary.vsel(): Only relevant if deltas is TRUE. For plot.vsel(): Always relevant. Either "ref" or "best", indicating whether the baseline is the reference model or the best submodel found (in terms of stats[1]), respectively.


Only relevant if any(stats %in% c("elpd", "mlpd")). The threshold for the ELPD difference (taking the submodel's ELPD minus the baseline model's ELPD) above which the submodel's ELPD is considered to be close enough to the baseline model's ELPD. An equivalent rule is applied in case of the MLPD. See suggest_size() for a formalization. Supplying NA deactivates this.


Arguments passed to the internal function which is used for bootstrapping (if applicable; see argument stats). Currently, relevant arguments are B (the number of bootstrap samples, defaulting to 2000) and seed (see set.seed(), defaulting to$integer.max, 1), but can also be NA to not call set.seed() at all).


As long as the reference model's performance is computable, it is always shown in the plot as a dashed red horizontal line. If baseline = "best", the baseline model's performance is shown as a dotted black horizontal line. If ! and any(stats %in% c("elpd", "mlpd")), the value supplied to thres_elpd (which is automatically adapted internally in case of the MLPD or deltas = FALSE) is shown as a dot-dashed gray horizontal line for the reference model and, if baseline = "best", as a long-dashed green horizontal line for the baseline model.


if (requireNamespace("rstanarm", quietly = TRUE)) {
  # Data:
  dat_gauss <- data.frame(y = df_gaussian$y, df_gaussian$x)

  # The "stanreg" fit which will be used as the reference model (with small
  # values for `chains` and `iter`, but only for technical reasons in this
  # example; this is not recommended in general):
  fit <- rstanarm::stan_glm(
    y ~ X1 + X2 + X3 + X4 + X5, family = gaussian(), data = dat_gauss,
    QR = TRUE, chains = 2, iter = 500, refresh = 0, seed = 9876

  # Variable selection (here without cross-validation and with small values
  # for `nterms_max`, `nclusters`, and `nclusters_pred`, but only for the
  # sake of speed in this example; this is not recommended in general):
  vs <- varsel(fit, nterms_max = 3, nclusters = 5, nclusters_pred = 10,
               seed = 5555)

projpred documentation built on Nov. 10, 2022, 5:15 p.m.