Package 'rstanarm'

Description

Stan Development Team

The rstanarm package is an appendage to the rstan package that enables many of the most common applied regression models to be estimated using Markov Chain Monte Carlo, variational approximations to the posterior distribution, or optimization. The rstanarm package allows these models to be specified using the customary R modeling syntax (e.g., like that of glm with a formula and a data.frame).

The sections below provide an overview of the modeling functions and estimation algorithms used by rstanarm.

Details

The set of models supported by rstanarm is large (and will continue to grow), but also limited enough so that it is possible to integrate them tightly with the pp_check function for graphical posterior predictive checks with bayesplot and the posterior_predict function to easily estimate the effect of specific manipulations of predictor variables or to predict the outcome in a training set.

The objects returned by the rstanarm modeling functions are called stanreg objects. In addition to all of the typical methods defined for fitted model objects, stanreg objects can be passed to the loo function in the loo package for model comparison or to the launch_shinystan function in the shinystan package in order to visualize the posterior distribution using the ShinyStan graphical user interface. See the rstanarm vignettes for more details about the entire process.

Prior distributions

See priors help page and the vignette Prior Distributions for rstanarm Models for an overview of the various choices the user can make for prior distributions. The package vignettes for the modeling functions also provide examples of using many of the available priors as well as more detailed descriptions of some of the novel priors used by rstanarm.

Modeling functions

The model estimating functions are described in greater detail in their individual help pages and vignettes. Here we provide a very brief overview:

stan_lm, stan_aov, stan_biglm

Similar to lm or aov but with novel regularizing priors on the model parameters that are driven by prior beliefs about $R^2$, the proportion of variance in the outcome attributable to the predictors in a linear model.

stan_glm, stan_glm.nb

Similar to glm but with various possible prior distributions for the coefficients and, if applicable, a prior distribution for any auxiliary parameter in a Generalized Linear Model (GLM) that is characterized by a family object (e.g. the shape parameter in Gamma models). It is also possible to estimate a negative binomial model in a similar way to the glm.nb function in the MASS package.

stan_glmer, stan_glmer.nb, stan_lmer

Similar to the glmer, glmer.nb and lmer functions in the lme4 package in that GLMs are augmented to have group-specific terms that deviate from the common coefficients according to a mean-zero multivariate normal distribution with a highly-structured but unknown covariance matrix (for which rstanarm introduces an innovative prior distribution). MCMC provides more appropriate estimates of uncertainty for models that consist of a mix of common and group-specific parameters.

stan_mvmer

A multivariate form of stan_glmer, whereby the user can specify one or more submodels each consisting of a GLM with group-specific terms. If more than one submodel is specified (i.e. there is more than one outcome variable) then a dependence is induced by assuming that the group-specific terms for each grouping factor are correlated across submodels.

stan_nlmer

Similar to nlmer in the lme4 package for nonlinear "mixed-effects" models, but the group-specific coefficients have flexible priors on their unknown covariance matrices.

stan_gamm4

Similar to gamm4 in the gamm4 package, which augments a GLM (possibly with group-specific terms) with nonlinear smooth functions of the predictors to form a Generalized Additive Mixed Model (GAMM). Rather than calling glmer like gamm4 does, stan_gamm4 essentially calls stan_glmer, which avoids the optimization issues that often crop up with GAMMs and provides better estimates for the uncertainty of the parameter estimates.

stan_polr

Similar to polr in the MASS package in that it models an ordinal response, but the Bayesian model also implies a prior distribution on the unknown cutpoints. Can also be used to model binary outcomes, possibly while estimating an unknown exponent governing the probability of success.

stan_betareg

Similar to betareg in that it models an outcome that is a rate (proportion) but, rather than performing maximum likelihood estimation, full Bayesian estimation is performed by default, with customizable prior distributions for all parameters.

stan_clogit

Similar to clogit in that it models an binary outcome where the number of successes and failures is fixed within each stratum by the research design. There are some minor syntactical differences relative to clogit that allow stan_clogit to accept group-specific terms as in stan_glmer.

stan_surv

Fits models to survival (i.e. time-to-event) data using either a hazard scale or accelerated failure time (AFT) formulation. The user can choose between either a flexible spline-based approximation for the baseline hazard or a standard parametric distributional assumption for the baseline hazard. Covariate effects can then be accommodated under proportional or non-proportional hazards assumptions (for models on the hazard scale) or using time-fixed or time-varying acceleration factors (for models on the AFT scale).

stan_jm

Fits shared parameter joint models for longitudinal and survival (i.e. time-to-event) data. The joint model can be univariate (i.e. one longitudinal outcome) or multivariate (i.e. more than one longitudinal outcome). A variety of parameterisations are available for linking the longitudinal and event processes (i.e. a variety of association structures).

Estimation algorithms

The modeling functions in the rstanarm package take an algorithm argument that can be one of the following:

Sampling (algorithm="sampling")

Uses Markov Chain Monte Carlo (MCMC) — in particular, Hamiltonian Monte Carlo (HMC) with a tuned but diagonal mass matrix — to draw from the posterior distribution of the parameters. See sampling (rstan) for more details. This is the slowest but most reliable of the available estimation algorithms and it is the default and recommended algorithm for statistical inference.

Mean-field (algorithm="meanfield")

Uses mean-field variational inference to draw from an approximation to the posterior distribution. In particular, this algorithm finds the set of independent normal distributions in the unconstrained space that — when transformed into the constrained space — most closely approximate the posterior distribution. Then it draws repeatedly from these independent normal distributions and transforms them into the constrained space. The entire process is much faster than HMC and yields independent draws but is not recommended for final statistical inference. It can be useful to narrow the set of candidate models in large problems, particularly when specifying QR=TRUE in stan_glm, stan_glmer, and stan_gamm4, but is only an approximation to the posterior distribution.

Full-rank (algorithm="fullrank")

Uses full-rank variational inference to draw from an approximation to the posterior distribution by finding the multivariate normal distribution in the unconstrained space that — when transformed into the constrained space — most closely approximates the posterior distribution. Then it draws repeatedly from this multivariate normal distribution and transforms the draws into the constrained space. This process is slower than meanfield variational inference but is faster than HMC. Although still an approximation to the posterior distribution and thus not recommended for final statistical inference, the approximation is more realistic than that of mean-field variational inference because the parameters are not assumed to be independent in the unconstrained space. Nevertheless, fullrank variational inference is a more difficult optimization problem and the algorithm is more prone to non-convergence or convergence to a local optimum.

Optimizing (algorithm="optimizing")

Finds the posterior mode using a C++ implementation of the LBGFS algorithm. See optimizing for more details. If there is no prior information, then this is equivalent to maximum likelihood, in which case there is no great reason to use the functions in the rstanarm package over the emulated functions in other packages. However, if priors are specified, then the estimates are penalized maximum likelihood estimates, which may have some redeeming value. Currently, optimization is only supported for stan_glm.

Author(s)

Maintainer: Ben Goodrich [email protected]

Authors:

• Jonah Gabry [email protected]

Other contributors:

• Imad Ali [contributor]

• Sam Brilleman [contributor]

• Jacqueline Buros Novik (R/stan_jm.R) [contributor]

• AstraZeneca (R/stan_jm.R) [contributor]

• Eren Elci (R/stan_surv.R) [contributor]

• Trustees of Columbia University [copyright holder]

• Simon Wood (R/stan_gamm4.R) [copyright holder]

• R Core Deveopment Team (R/stan_aov.R) [copyright holder]

• Douglas Bates (R/pp_data.R) [copyright holder]

• Martin Maechler (R/pp_data.R) [copyright holder]

• Ben Bolker (R/pp_data.R) [copyright holder]

• Steve Walker (R/pp_data.R) [copyright holder]

• Brian Ripley (R/stan_aov.R, R/stan_polr.R) [copyright holder]

• William Venables (R/stan_polr.R) [copyright holder]

• Paul-Christian Burkner [email protected] (R/misc.R) [copyright holder]

References

Bates, D., Maechler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-Effects models using lme4. Journal of Statistical Software. 67(1), 1–48.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian Data Analysis. Chapman & Hall/CRC Press, London, third edition. https://stat.columbia.edu/~gelman/book/

Gelman, A. and Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, Cambridge, UK. https://stat.columbia.edu/~gelman/arm/

Stan Development Team. Stan Modeling Language Users Guide and Reference Manual. https://mc-stan.org/users/documentation/.

Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing. 27(5), 1413–1432. doi:10.1007/s11222-016-9696-4. arXiv preprint: https://arxiv.org/abs/1507.04544

Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018) Using stacking to average Bayesian predictive distributions. Bayesian Analysis, advance publication, doi:10.1214/17-BA1091.

Gabry, J. , Simpson, D. , Vehtari, A. , Betancourt, M. and Gelman, A. (2019), Visualization in Bayesian workflow. J. R. Stat. Soc. A, 182: 389-402. doi:10.1111/rssa.12378, arXiv preprint, code on GitHub)

Muth, C., Oravecz, Z., and Gabry, J. (2018) User-friendly Bayesian regression modeling: A tutorial with rstanarm and shinystan. The Quantitative Methods for Psychology. 14(2), 99–119. https://www.tqmp.org/RegularArticles/vol14-2/p099/p099.pdf

See Also

• https://mc-stan.org/ for more information on the Stan C++ package used by rstanarm for model fitting.

• https://github.com/stan-dev/rstanarm/issues/ to submit a bug report or feature request.

• https://discourse.mc-stan.org to ask a question about rstanarm on the Stan-users forum.

Description

Details about the adapt_delta argument to rstanarm's modeling functions.

Details

For the No-U-Turn Sampler (NUTS), the variant of Hamiltonian Monte Carlo used used by rstanarm, adapt_delta is the target average proposal acceptance probability during Stan's adaptation period. adapt_delta is ignored by rstanarm if the algorithm argument is not set to "sampling".

The default value of adapt_delta is 0.95, except when the prior for the regression coefficients is R2, hs, or hs_plus, in which case the default is 0.99.

These defaults are higher (more conservative) than the default of adapt_delta=0.8 used in the rstan package, which may result in slower sampling speeds but will be more robust to posterior distributions with high curvature.

In general you should not need to change adapt_delta unless you see a warning message about divergent transitions, in which case you can increase adapt_delta from the default to a value closer to 1 (e.g. from 0.95 to 0.99, or from 0.99 to 0.999, etc). The step size used by the numerical integrator is a function of adapt_delta in that increasing adapt_delta will result in a smaller step size and fewer divergences. Increasing adapt_delta will typically result in a slower sampler, but it will always lead to a more robust sampler.

References

Stan Development Team. Stan Modeling Language Users Guide and Reference Manual. https://mc-stan.org/users/documentation/.

Brief Guide to Stan's Warnings: https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

Description

For models fit using MCMC (algorithm="sampling"), the posterior sample —the post-warmup draws from the posterior distribution— can be extracted from a fitted model object as a matrix, data frame, or array. The as.matrix and as.data.frame methods merge all chains together, whereas the as.array method keeps the chains separate. For models fit using optimization ("optimizing") or variational inference ("meanfield" or "fullrank"), there is no posterior sample but rather a matrix (or data frame) of 1000 draws from either the asymptotic multivariate Gaussian sampling distribution of the parameters or the variational approximation to the posterior distribution.

Usage

## S3 method for class 'stanreg'
as.matrix(x, ..., pars = NULL, regex_pars = NULL)

## S3 method for class 'stanreg'
as.array(x, ..., pars = NULL, regex_pars = NULL)

## S3 method for class 'stanreg'
as.data.frame(x, ..., pars = NULL, regex_pars = NULL)

Arguments

Value

A matrix, data.frame, or array, the dimensions of which depend on pars and regex_pars, as well as the model and estimation algorithm (see the Description section above).

See Also

stanreg-draws-formats, stanreg-methods

Examples

if (.Platform$OS.type != "windows") {

if (!exists("example_model")) example(example_model)
# Extract posterior sample after MCMC
draws <- as.matrix(example_model)
print(dim(draws))

# For example, we can see that the median of the draws for the intercept
# is the same as the point estimate rstanarm uses
print(median(draws[, "(Intercept)"]))
print(example_model$coefficients[["(Intercept)"]])

# The as.array method keeps the chains separate
draws_array <- as.array(example_model)
print(dim(draws_array)) # iterations x chains x parameters

# Extract draws from asymptotic Gaussian sampling distribution
# after optimization
fit <- stan_glm(mpg ~ wt, data = mtcars, algorithm = "optimizing")
draws <- as.data.frame(fit)
print(colnames(draws))
print(nrow(draws)) # 1000 draws are taken

# Extract draws from variational approximation to the posterior distribution
fit2 <- update(fit, algorithm = "meanfield")
draws <- as.data.frame(fit2, pars = "wt")
print(colnames(draws))
print(nrow(draws)) # 1000 draws are taken

}

Description

Estimation algorithms available for rstanarm models

Estimation algorithms

The modeling functions in the rstanarm package take an algorithm argument that can be one of the following:

Sampling (algorithm="sampling")

Uses Markov Chain Monte Carlo (MCMC) — in particular, Hamiltonian Monte Carlo (HMC) with a tuned but diagonal mass matrix — to draw from the posterior distribution of the parameters. See sampling (rstan) for more details. This is the slowest but most reliable of the available estimation algorithms and it is the default and recommended algorithm for statistical inference.

Mean-field (algorithm="meanfield")

Uses mean-field variational inference to draw from an approximation to the posterior distribution. In particular, this algorithm finds the set of independent normal distributions in the unconstrained space that — when transformed into the constrained space — most closely approximate the posterior distribution. Then it draws repeatedly from these independent normal distributions and transforms them into the constrained space. The entire process is much faster than HMC and yields independent draws but is not recommended for final statistical inference. It can be useful to narrow the set of candidate models in large problems, particularly when specifying QR=TRUE in stan_glm, stan_glmer, and stan_gamm4, but is only an approximation to the posterior distribution.

Full-rank (algorithm="fullrank")

Uses full-rank variational inference to draw from an approximation to the posterior distribution by finding the multivariate normal distribution in the unconstrained space that — when transformed into the constrained space — most closely approximates the posterior distribution. Then it draws repeatedly from this multivariate normal distribution and transforms the draws into the constrained space. This process is slower than meanfield variational inference but is faster than HMC. Although still an approximation to the posterior distribution and thus not recommended for final statistical inference, the approximation is more realistic than that of mean-field variational inference because the parameters are not assumed to be independent in the unconstrained space. Nevertheless, fullrank variational inference is a more difficult optimization problem and the algorithm is more prone to non-convergence or convergence to a local optimum.

Optimizing (algorithm="optimizing")

Finds the posterior mode using a C++ implementation of the LBGFS algorithm. See optimizing for more details. If there is no prior information, then this is equivalent to maximum likelihood, in which case there is no great reason to use the functions in the rstanarm package over the emulated functions in other packages. However, if priors are specified, then the estimates are penalized maximum likelihood estimates, which may have some redeeming value. Currently, optimization is only supported for stan_glm.

See Also

https://mc-stan.org/rstanarm/

Description

Modeling functions available in rstanarm

Modeling functions

The model estimating functions are described in greater detail in their individual help pages and vignettes. Here we provide a very brief overview:

stan_lm, stan_aov, stan_biglm

Similar to lm or aov but with novel regularizing priors on the model parameters that are driven by prior beliefs about $R^2$, the proportion of variance in the outcome attributable to the predictors in a linear model.

stan_glm, stan_glm.nb

Similar to glm but with various possible prior distributions for the coefficients and, if applicable, a prior distribution for any auxiliary parameter in a Generalized Linear Model (GLM) that is characterized by a family object (e.g. the shape parameter in Gamma models). It is also possible to estimate a negative binomial model in a similar way to the glm.nb function in the MASS package.

stan_glmer, stan_glmer.nb, stan_lmer

Similar to the glmer, glmer.nb and lmer functions in the lme4 package in that GLMs are augmented to have group-specific terms that deviate from the common coefficients according to a mean-zero multivariate normal distribution with a highly-structured but unknown covariance matrix (for which rstanarm introduces an innovative prior distribution). MCMC provides more appropriate estimates of uncertainty for models that consist of a mix of common and group-specific parameters.

stan_mvmer

A multivariate form of stan_glmer, whereby the user can specify one or more submodels each consisting of a GLM with group-specific terms. If more than one submodel is specified (i.e. there is more than one outcome variable) then a dependence is induced by assuming that the group-specific terms for each grouping factor are correlated across submodels.

stan_nlmer

Similar to nlmer in the lme4 package for nonlinear "mixed-effects" models, but the group-specific coefficients have flexible priors on their unknown covariance matrices.

stan_gamm4

Similar to gamm4 in the gamm4 package, which augments a GLM (possibly with group-specific terms) with nonlinear smooth functions of the predictors to form a Generalized Additive Mixed Model (GAMM). Rather than calling glmer like gamm4 does, stan_gamm4 essentially calls stan_glmer, which avoids the optimization issues that often crop up with GAMMs and provides better estimates for the uncertainty of the parameter estimates.

stan_polr

Similar to polr in the MASS package in that it models an ordinal response, but the Bayesian model also implies a prior distribution on the unknown cutpoints. Can also be used to model binary outcomes, possibly while estimating an unknown exponent governing the probability of success.

stan_betareg

Similar to betareg in that it models an outcome that is a rate (proportion) but, rather than performing maximum likelihood estimation, full Bayesian estimation is performed by default, with customizable prior distributions for all parameters.

stan_clogit

Similar to clogit in that it models an binary outcome where the number of successes and failures is fixed within each stratum by the research design. There are some minor syntactical differences relative to clogit that allow stan_clogit to accept group-specific terms as in stan_glmer.

stan_surv

Fits models to survival (i.e. time-to-event) data using either a hazard scale or accelerated failure time (AFT) formulation. The user can choose between either a flexible spline-based approximation for the baseline hazard or a standard parametric distributional assumption for the baseline hazard. Covariate effects can then be accommodated under proportional or non-proportional hazards assumptions (for models on the hazard scale) or using time-fixed or time-varying acceleration factors (for models on the AFT scale).

stan_jm

Fits shared parameter joint models for longitudinal and survival (i.e. time-to-event) data. The joint model can be univariate (i.e. one longitudinal outcome) or multivariate (i.e. more than one longitudinal outcome). A variety of parameterisations are available for linking the longitudinal and event processes (i.e. a variety of association structures).

See Also

https://mc-stan.org/rstanarm/

Description

Compute a Bayesian version of R-squared or LOO-adjusted R-squared for regression models.

Usage

## S3 method for class 'stanreg'
bayes_R2(object, ..., re.form = NULL)

## S3 method for class 'stanreg'
loo_R2(object, ...)

Arguments

Value

A vector of R-squared values with length equal to the posterior sample size (the posterior distribution of R-squared).

References

Andrew Gelman, Ben Goodrich, Jonah Gabry, and Aki Vehtari (2018). R-squared for Bayesian regression models. The American Statistician, to appear. doi:10.1080/00031305.2018.1549100 (Preprint, Notebook)

Examples

if (.Platform$OS.type != "windows") {
fit <- stan_glm(
  mpg ~ wt + cyl,
  data = mtcars,
  QR = TRUE,
  chains = 2,
  refresh = 0
)
rsq <- bayes_R2(fit)
print(median(rsq))
hist(rsq)

loo_rsq <- loo_R2(fit)
print(median(loo_rsq))

# multilevel binomial model
if (!exists("example_model")) example(example_model)
print(example_model)
median(bayes_R2(example_model))
median(bayes_R2(example_model, re.form = NA)) # exclude group-level
}

Description

Collapse the linear predictor across the lower level units clustered an individual, using the function specified in the 'grp_assoc' argument

Usage

collapse_within_groups(eta, grp_idx, grp_assoc = "sum")

Arguments

Value

A vector or matrix, depending on the method called.

Description

Collapse the linear predictor across the lower level units clustered an individual, using the function specified in the 'grp_assoc' argument

Usage

## Default S3 method:
collapse_within_groups(eta, grp_idx, grp_assoc)

Arguments

Value

A vector or matrix, depending on the method called.

Description

Collapse the linear predictor across the lower level units clustered an individual, using the function specified in the 'grp_assoc' argument

Usage

## S3 method for class 'matrix'
collapse_within_groups(eta, grp_idx, grp_assoc)

Arguments

Value

A vector or matrix, depending on the method called.

Description

Evaluate the log baseline hazard at the specified times given the vector or matrix of MCMC draws for the baseline hazard coeffients / parameters

Usage

evaluate_log_survival(log_haz, qnodes, qwts)

Arguments

Value

A vector or matrix of log survival probabilities.

Description

Evaluate the log baseline hazard at the specified times given the vector or matrix of MCMC draws for the baseline hazard coeffients / parameters

Usage

## Default S3 method:
evaluate_log_survival(log_haz, qnodes, qwts)

Arguments

Value

A vector or matrix of log survival probabilities.

Description

Evaluate the log baseline hazard at the specified times given the vector or matrix of MCMC draws for the baseline hazard coeffients / parameters

Usage

## S3 method for class 'matrix'
evaluate_log_survival(log_haz, qnodes, qwts)

Arguments

Value

A vector or matrix of log survival probabilities.

Description

A model for use in the rstanarm examples related to stan_jm.

Format

Calling example("example_jm") will run the model in the Examples section, below, and the resulting stanmvreg object will then be available in the global environment. The chains and iter arguments are specified to make this example be small in size. In practice, we recommend that they be left unspecified in order to use the default values or increased if there are convergence problems. The cores argument is optional and on a multicore system, the user may well want to set that equal to the number of chains being executed.

Examples

# set.seed(123)
  if (.Platform$OS.type != "windows" || .Platform$r_arch !="i386")
  example_jm <- 
     stan_jm(formulaLong = logBili ~ year + (1 | id), 
             dataLong = pbcLong[1:101,],
             formulaEvent = survival::Surv(futimeYears, death) ~ sex + trt, 
             dataEvent = pbcSurv[1:15,],
             time_var = "year",
             # this next line is only to keep the example small in size!
             chains = 1, seed = 12345, iter = 100, refresh = 0)

Description

A model for use in rstanarm examples.

Format

Calling example("example_model") will run the model in the Examples section, below, and the resulting stanreg object will then be available in the global environment. The chains and iter arguments are specified to make this example be small in size. In practice, we recommend that they be left unspecified in order to use the default values (4 and 2000 respectively) or increased if there are convergence problems. The cores argument is optional and on a multicore system, the user may well want to set that equal to the number of chains being executed.

See Also

cbpp for a description of the data.

Examples

if (.Platform$OS.type != "windows") {
example_model <-
  stan_glmer(cbind(incidence, size - incidence) ~ size + period + (1|herd),
             data = lme4::cbpp, family = binomial, QR = TRUE,
             # this next line is only to keep the example small in size!
             chains = 2, cores = 1, seed = 12345, iter = 1000, refresh = 0)
example_model
}

Description

Extract parameters from stanmat and return as a list

Usage

extract_pars(object, ...)

Arguments

Value

A named list

Description

Extract parameters from stanmat and return as a list

Usage

## S3 method for class 'stanmvreg'
extract_pars(object, stanmat = NULL, means = FALSE, ...)

Arguments

Value

A named list

Description

Extract parameters from stanmat and return as a list

Usage

## S3 method for class 'stansurv'
extract_pars(object, stanmat = NULL, means = FALSE, ...)

Arguments

Value

A named list

Description

It is necessary to drop unneeded variables though so that errors are not encountered if the original data contained NA values for variables unrelated to the model formula. We generate a data frame here for in-sample predictions rather than using a model frame, since some quantities will need to be recalculated at quadrature points etc, for example in posterior_survfit.

Usage

get_model_data(object, ...)

Arguments

Value

A data frame or list of data frames with all the (unevaluated) variables required for predictions.

Description

It is necessary to drop unneeded variables though so that errors are not encountered if the original data contained NA values for variables unrelated to the model formula. We generate a data frame here for in-sample predictions rather than using a model frame, since some quantities will need to be recalculated at quadrature points etc, for example in posterior_survfit.

Usage

## S3 method for class 'stanmvreg'
get_model_data(object, m = NULL, ...)

Arguments

Value

A data frame or list of data frames with all the (unevaluated) variables required for predictions.

Description

It is necessary to drop unneeded variables though so that errors are not encountered if the original data contained NA values for variables unrelated to the model formula. We generate a data frame here for in-sample predictions rather than using a model frame, since some quantities will need to be recalculated at quadrature points etc, for example in posterior_survfit.

Usage

## S3 method for class 'stansurv'
get_model_data(object, ...)

Arguments

Value

A data frame or list of data frames with all the (unevaluated) variables required for predictions.

Description

The kfold method performs exact $K$-fold cross-validation. First the data are randomly partitioned into $K$ subsets of equal size (or as close to equal as possible), or the user can specify the folds argument to determine the partitioning. Then the model is refit $K$ times, each time leaving out one of the $K$ subsets. If $K$ is equal to the total number of observations in the data then $K$-fold cross-validation is equivalent to exact leave-one-out cross-validation (to which loo is an efficient approximation).

Usage

## S3 method for class 'stanreg'
kfold(
  x,
  K = 10,
  ...,
  folds = NULL,
  save_fits = FALSE,
  cores = getOption("mc.cores", 1)
)

Arguments

Value

An object with classes 'kfold' and 'loo' that has a similar structure as the objects returned by the loo and waic methods and is compatible with the loo_compare function for comparing models.

References

Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing. 27(5), 1413–1432. doi:10.1007/s11222-016-9696-4. arXiv preprint: https://arxiv.org/abs/1507.04544

Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018) Using stacking to average Bayesian predictive distributions. Bayesian Analysis, advance publication, doi:10.1214/17-BA1091.

Examples

if (.Platform$OS.type != "windows") {

fit1 <- stan_glm(mpg ~ wt, data = mtcars, refresh = 0)
fit2 <- stan_glm(mpg ~ wt + cyl, data = mtcars, refresh = 0)
fit3 <- stan_glm(mpg ~ disp * as.factor(cyl), data = mtcars, refresh = 0)

# 10-fold cross-validation
# (if possible also specify the 'cores' argument to use multiple cores)
(kfold1 <- kfold(fit1, K = 10))
kfold2 <- kfold(fit2, K = 10)
kfold3 <- kfold(fit3, K = 10)
loo_compare(kfold1, kfold2, kfold3)

# stratifying by a grouping variable
# (note: might get some divergences warnings with this model but
# this is just intended as a quick example of how to code this)
fit4 <- stan_lmer(mpg ~ disp + (1|cyl), data = mtcars, refresh = 0)
table(mtcars$cyl)
folds_cyl <- loo::kfold_split_stratified(K = 3, x = mtcars$cyl)
table(cyl = mtcars$cyl, fold = folds_cyl)
kfold4 <- kfold(fit4, folds = folds_cyl, cores = 2)
print(kfold4)

}
# Example code demonstrating the different ways to specify the number
# of cores and how the cores are used
#
# options(mc.cores = NULL)
#
# # spread the K models over N_CORES cores (method 1)
# kfold(fit, K, cores = N_CORES)
#
# # spread the K models over N_CORES cores (method 2)
# options(mc.cores = N_CORES)
# kfold(fit, K)
#
# # fit K models sequentially using N_CORES cores for the Markov chains each time
# options(mc.cores = N_CORES)
# kfold(fit, K, cores = 1)

Description

The ShinyStan interface provides visual and numerical summaries of model parameters and convergence diagnostics.

Usage

## S3 method for class 'stanreg'
launch_shinystan(
  object,
  ppd = TRUE,
  seed = 1234,
  model_name = NULL,
  note = NULL,
  rstudio = getOption("shinystan.rstudio"),
  ...
)

Arguments

Details

The launch_shinystan function will accept a stanreg object as input. Currently, almost any model fit using one of rstanarm's model-fitting functions can be used with ShinyStan. The only exception is that ShinyStan does not currently support rstanarm models fit using algorithm='optimizing'. See the shinystan package documentation for more information.

Faster launch times

For some rstanarm models ShinyStan may take a very long time to launch. If this is the case with one of your models you may be able to speed up launch_shinystan in one of several ways:

Prevent ShinyStan from preparing graphical posterior predictive checks:

When used with a stanreg object (rstanarm model object) ShinyStan will draw from the posterior predictive distribution and prepare graphical posterior predictive checks before launching. That way when you go to the PPcheck page the plots are immediately available. This can be time consuming for models fit to very large datasets and you can prevent this behavior by creating a shinystan object before calling launch_shinystan. To do this use as.shinystan with optional argument ppd set to FALSE (see the Examples section below). When you then launch ShinyStan and go to the PPcheck page the plots will no longer be automatically generated and you will be presented with the standard interface requiring you to first specify the appropriate $y$ and $yrep$, which can be done for many but not all rstanarm models.

Use a shinystan object:

Even if you don't want to prevent ShinyStan from preparing graphical posterior predictive checks, first creating a shinystan object using as.shinystan can reduce future launch times. That is, launch_shinystan(sso) will be faster than launch_shinystan(fit), where sso is a shinystan object and fit is a stanreg object. It still may take some time for as.shinystan to create sso initially, but each time you subsequently call launch_shinystan(sso) it will reuse sso instead of internally creating a shinystan object every time. See the Examples section below.

References

Gabry, J. , Simpson, D. , Vehtari, A. , Betancourt, M. and Gelman, A. (2019), Visualization in Bayesian workflow. J. R. Stat. Soc. A, 182: 389-402. doi:10.1111/rssa.12378, arXiv preprint, code on GitHub)

Muth, C., Oravecz, Z., and Gabry, J. (2018) User-friendly Bayesian regression modeling: A tutorial with rstanarm and shinystan. The Quantitative Methods for Psychology. 14(2), 99–119. https://www.tqmp.org/RegularArticles/vol14-2/p099/p099.pdf

Examples

if (.Platform$OS.type != "windows") {
## Not run: 
if (!exists("example_model")) example(example_model)

# Launch the ShinyStan app without saving the resulting shinystan object
if (interactive()) launch_shinystan(example_model)

# Launch the ShinyStan app (saving resulting shinystan object as sso)
if (interactive()) sso <- launch_shinystan(example_model)

# First create shinystan object then call launch_shinystan
sso <- shinystan::as.shinystan(example_model)
if (interactive()) launch_shinystan(sso)

# Prevent ShinyStan from preparing graphical posterior predictive checks that
# can be time consuming. example_model is small enough that it won't matter
# much here but in general this can help speed up launch_shinystan
sso <- shinystan::as.shinystan(example_model, ppd = FALSE)
if (interactive()) launch_shinystan(sso)

## End(Not run)
}

Description

Make linear predictor vector from x and point estimates for beta, or linear predictor matrix from x and full posterior sample of beta.

Usage

linear_predictor(beta, x, offset = NULL)

Arguments

Value

A vector or matrix.

Description

Make linear predictor vector from x and point estimates for beta, or linear predictor matrix from x and full posterior sample of beta.

Usage

## Default S3 method:
linear_predictor(beta, x, offset = NULL)

Arguments

Value

A vector or matrix.

Description

Make linear predictor vector from x and point estimates for beta, or linear predictor matrix from x and full posterior sample of beta.

Usage

## S3 method for class 'matrix'
linear_predictor(beta, x, offset = NULL)

Arguments

Value

A vector or matrix.

Description

Get inverse link function

Usage

linkinv(x, ...)

Arguments

Value

The inverse link function associated with x.

Description

Get inverse link function

Usage

## S3 method for class 'character'
linkinv(x, ...)

Arguments

Value

The inverse link function associated with x.

Description

Get inverse link function

Usage

## S3 method for class 'family'
linkinv(x, ...)

Arguments

Value

The inverse link function associated with x.

Description

Get inverse link function

Usage

## S3 method for class 'stanmvreg'
linkinv(x, m = NULL, ...)

Arguments

Value

The inverse link function associated with x.

Description

Get inverse link function

Usage

## S3 method for class 'stanreg'
linkinv(x, ...)

Arguments

Value

The inverse link function associated with x.

Description

get arguments needed for ll_fun

Usage

ll_args(object, ...)

Arguments

Value

a named list with elements data, draws, S (posterior sample size) and N = number of observations

Description

Alternative ll_args method for stanjm objects that allows data and pars to be passed directly, rather than constructed using pp_data within the ll_args method. This can be much faster when used in the MH algorithm within posterior_survfit, since it doesn't require repeated calls to pp_data.

Usage

## S3 method for class 'stanjm'
ll_args(object, data, pars, m = 1, reloo_or_kfold = FALSE, ...)

Arguments

Description

get arguments needed for ll_fun

Usage

## S3 method for class 'stanreg'
ll_args(
  object,
  newdata = NULL,
  offset = NULL,
  m = NULL,
  reloo_or_kfold = FALSE,
  ...
)

Arguments

Value

a named list with elements data, draws, S (posterior sample size) and N = number of observations

Description

get arguments needed for ll_fun

Usage

## S3 method for class 'stansurv'
ll_args(object, newdata = NULL, ...)

Arguments

Value

a named list with elements data, draws, S (posterior sample size) and N = number of observations

Description

For models fit using MCMC only, the log_lik method returns the $S$ by $N$ pointwise log-likelihood matrix, where $S$ is the size of the posterior sample and $N$ is the number of data points, or in the case of the stanmvreg method (when called on stan_jm model objects) an $S$ by $Npat$ matrix where $Npat$ is the number of individuals.

Usage

## S3 method for class 'stanreg'
log_lik(object, newdata = NULL, offset = NULL, ...)

## S3 method for class 'stanmvreg'
log_lik(object, m = 1, newdata = NULL, ...)

## S3 method for class 'stanjm'
log_lik(object, newdataLong = NULL, newdataEvent = NULL, ...)

Arguments

Value

For the stanreg and stanmvreg methods an $S$ by $N$ matrix, where $S$ is the size of the posterior sample and $N$ is the number of data points. For the stanjm method an $S$ by $Npat$ matrix where $Npat$ is the number of individuals.

Examples

if (.Platform$OS.type != "windows") {

 roaches$roach100 <- roaches$roach1 / 100
 fit <- stan_glm(
    y ~ roach100 + treatment + senior,
    offset = log(exposure2),
    data = roaches,
    family = poisson(link = "log"),
    prior = normal(0, 2.5),
    prior_intercept = normal(0, 10),
    iter = 500, # just to speed up example,
    refresh = 0
 )
 ll <- log_lik(fit)
 dim(ll)
 all.equal(ncol(ll), nobs(fit))

 # using newdata argument
 nd <- roaches[1:2, ]
 nd$treatment[1:2] <- c(0, 1)
 ll2 <- log_lik(fit, newdata = nd, offset = c(0, 0))
 head(ll2)
 dim(ll2)
 all.equal(ncol(ll2), nrow(nd))

}

Description

Logit and inverse logit

Usage

logit(x)

invlogit(x)

Arguments

Value

A numeric vector the same length as x.

Description

These functions are wrappers around the E_loo function (loo package) that provide compatibility for rstanarm models.

Usage

## S3 method for class 'stanreg'
loo_predict(
  object,
  type = c("mean", "var", "quantile"),
  probs = 0.5,
  ...,
  psis_object = NULL
)

## S3 method for class 'stanreg'
loo_linpred(
  object,
  type = c("mean", "var", "quantile"),
  probs = 0.5,
  transform = FALSE,
  ...,
  psis_object = NULL
)

## S3 method for class 'stanreg'
loo_predictive_interval(object, prob = 0.9, ..., psis_object = NULL)

Arguments

Value

A list with elements value and pareto_k.

For loo_predict and loo_linpred the value component is a vector with one element per observation.

For loo_predictive_interval the value component is a matrix with one row per observation and two columns (like predictive_interval). loo_predictive_interval(..., prob = p) is equivalent to loo_predict(..., type = "quantile", probs = c(a, 1-a)) with a = (1 - p)/2, except it transposes the result and adds informative column names.

See E_loo and pareto-k-diagnostic for details on the pareto_k diagnostic.

References

Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing. 27(5), 1413–1432. doi:10.1007/s11222-016-9696-4. arXiv preprint: https://arxiv.org/abs/1507.04544

Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018) Using stacking to average Bayesian predictive distributions. Bayesian Analysis, advance publication, doi:10.1214/17-BA1091.

Gabry, J. , Simpson, D. , Vehtari, A. , Betancourt, M. and Gelman, A. (2019), Visualization in Bayesian workflow. J. R. Stat. Soc. A, 182: 389-402. doi:10.1111/rssa.12378, arXiv preprint, code on GitHub)

Examples

if (.Platform$OS.type != "windows") {
## Not run: 
if (!exists("example_model")) example(example_model)

# optionally, log-weights can be pre-computed and reused
psis_result <- loo::psis(log_ratios = -log_lik(example_model))

loo_probs <- loo_linpred(example_model, type = "mean", transform = TRUE, psis_object = psis_result)
str(loo_probs)

loo_pred_var <- loo_predict(example_model, type = "var", psis_object = psis_result)
str(loo_pred_var)

loo_pred_ints <- loo_predictive_interval(example_model, prob = 0.8, psis_object = psis_result)
str(loo_pred_ints)

## End(Not run)
}

Description

For models fit using MCMC, compute approximate leave-one-out cross-validation (LOO, LOOIC) or, less preferably, the Widely Applicable Information Criterion (WAIC) using the loo package. (For $K$-fold cross-validation see kfold.stanreg.) Functions for model comparison, and model weighting/averaging are also provided.

Note: these functions are not guaranteed to work properly unless the data argument was specified when the model was fit. Also, as of loo version 2.0.0 the default number of cores is now only 1, but we recommend using as many (or close to as many) cores as possible by setting the cores argument or using options(mc.cores = VALUE) to set it for an entire session.

Usage

## S3 method for class 'stanreg'
loo(
  x,
  ...,
  cores = getOption("mc.cores", 1),
  save_psis = FALSE,
  k_threshold = NULL
)

## S3 method for class 'stanreg'
waic(x, ...)

## S3 method for class 'stanreg'
loo_compare(x, ..., criterion = c("loo", "kfold", "waic"), detail = FALSE)

## S3 method for class 'stanreg_list'
loo_compare(x, ..., criterion = c("loo", "kfold", "waic"), detail = FALSE)

## S3 method for class 'stanreg_list'
loo_model_weights(x, ..., cores = getOption("mc.cores", 1), k_threshold = NULL)

compare_models(..., loos = list(), detail = FALSE)

Arguments

Value

The structure of the objects returned by loo and waic methods are documented in detail in the Value section in loo and waic (from the loo package).

loo_compare returns a matrix with class 'compare.loo'. See the Comparing models section below for more details.

Approximate LOO CV

The loo method for stanreg objects provides an interface to the loo package for approximate leave-one-out cross-validation (LOO). The LOO Information Criterion (LOOIC) has the same purpose as the Akaike Information Criterion (AIC) that is used by frequentists. Both are intended to estimate the expected log predictive density (ELPD) for a new dataset. However, the AIC ignores priors and assumes that the posterior distribution is multivariate normal, whereas the functions from the loo package do not make this distributional assumption and integrate over uncertainty in the parameters. This only assumes that any one observation can be omitted without having a major effect on the posterior distribution, which can be judged using the diagnostic plot provided by the plot.loo method and the warnings provided by the print.loo method (see the How to Use the rstanarm Package vignette for an example of this process).

How to proceed when loo gives warnings (k_threshold)

The k_threshold argument to the loo method for rstanarm models is provided as a possible remedy when the diagnostics reveal problems stemming from the posterior's sensitivity to particular observations. Warnings about Pareto $k$ estimates indicate observations for which the approximation to LOO is problematic (this is described in detail in Vehtari, Gelman, and Gabry (2017) and the loo package documentation). The k_threshold argument can be used to set the $k$ value above which an observation is flagged. If k_threshold is not NULL and there are $J$ observations with $k$ estimates above k_threshold then when loo is called it will refit the original model $J$ times, each time leaving out one of the $J$ problematic observations. The pointwise contributions of these observations to the total ELPD are then computed directly and substituted for the previous estimates from these $J$ observations that are stored in the object created by loo. Another option to consider is K-fold cross-validation, which is documented on a separate page (see kfold).

Note: in the warning messages issued by loo about large Pareto $k$ estimates we recommend setting k_threshold to at least $0.7$. There is a theoretical reason, explained in Vehtari, Gelman, and Gabry (2017), for setting the threshold to the stricter value of $0.5$, but in practice they find that errors in the LOO approximation start to increase non-negligibly when $k > 0.7$.

Comparing models

"loo" (or "waic" or "kfold") objects can be passed to the loo_compare function in the loo package to perform model comparison. rstanarm also provides a loo_compare.stanreg method that can be used if the "loo" (or "waic" or "kfold") object has been added to the fitted model object (see the Examples section below for how to do this). This second method allows rstanarm to perform some extra checks that can't be done by the loo package itself (e.g., verifying that all models to be compared were fit using the same outcome variable).

loo_compare will return a matrix with one row per model and columns containing the ELPD difference and the standard error of the difference. In the first row of the matrix will be the model with the largest ELPD (smallest LOOIC) and will contain zeros (there is no difference between this model and itself). For each of the remaining models the ELPD difference and SE are reported relative to the model with the best ELPD (the first row). See the Details section at the loo_compare page in the loo package for more information.

Model weights

The loo_model_weights method can be used to compute model weights for a "stanreg_list" object, which is a list of fitted model objects made with stanreg_list. The end of the Examples section has a demonstration. For details see the loo_model_weights documentation in the loo package.

References

Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing. 27(5), 1413–1432. doi:10.1007/s11222-016-9696-4. arXiv preprint: https://arxiv.org/abs/1507.04544

Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018) Using stacking to average Bayesian predictive distributions. Bayesian Analysis, advance publication, doi:10.1214/17-BA1091.

Gabry, J. , Simpson, D. , Vehtari, A. , Betancourt, M. and Gelman, A. (2019), Visualization in Bayesian workflow. J. R. Stat. Soc. A, 182: 389-402. doi:10.1111/rssa.12378, arXiv preprint, code on GitHub)

See Also

• The new loo package vignettes and various rstanarm vignettes for more examples using loo and related functions with rstanarm models.

• pareto-k-diagnostic in the loo package for more on Pareto $k$ diagnostics.

• log_lik.stanreg to directly access the pointwise log-likelihood matrix.

Examples

if (.Platform$OS.type != "windows") {

fit1 <- stan_glm(mpg ~ wt, data = mtcars, refresh = 0)
fit2 <- stan_glm(mpg ~ wt + cyl, data = mtcars, refresh = 0)

# (for bigger models use as many cores as possible)
loo1 <- loo(fit1, cores = 1)
print(loo1)
loo2 <- loo(fit2, cores = 1)
print(loo2)

# when comparing models the loo objects can be passed to loo_compare
# as individual arguments or as a list of loo objects
loo_compare(loo1, loo2)
loo_compare(list(loo1, loo2))

# if the fitted model objects contain a loo object in the component "loo"
# then the model objects can be passed directly or as a stanreg_list
fit1$loo <- loo1
fit2$loo <- loo2
loo_compare(fit1, fit2)

# if the fitted model objects contain a loo object _and_ a waic or kfold
# object, then the criterion argument determines which of them the comparison
# is based on
fit1$waic <- waic(fit1)
fit2$waic <- waic(fit2)
loo_compare(fit1, fit2, criterion = "waic")

# the models can also be combined into a stanreg_list object, and more
# informative model names can be provided to use when printing
model_list <- stanreg_list(fit1, fit2, model_names = c("Fewer predictors", "More predictors"))
loo_compare(model_list)

fit3 <- stan_glm(mpg ~ disp * as.factor(cyl), data = mtcars, refresh = 0)
loo3 <- loo(fit3, cores = 2, k_threshold = 0.7)
loo_compare(loo1, loo2, loo3)

# setting detail=TRUE will also print model formulas if used with
# loo_compare.stanreg or loo_compare.stanreg_list
fit3$loo <- loo3
model_list <- stanreg_list(fit1, fit2, fit3)
loo_compare(model_list, detail=TRUE)

# Computing model weights
#
# if the objects in model_list already have 'loo' components then those
# will be used. otherwise loo will be computed for each model internally
# (in which case the 'cores' argument may also be used and is passed to loo())
loo_model_weights(model_list)  # defaults to method="stacking"
loo_model_weights(model_list,  method = "pseudobma")
loo_model_weights(model_list,  method = "pseudobma", BB = FALSE)

# you can also pass precomputed loo objects directly to loo_model_weights
loo_list <- list(A = loo1, B = loo2, C = loo3) # names optional (affects printing)
loo_model_weights(loo_list)

}

Description

Specifies the information required to fit a Negative Binomial GLM in a similar way to negative.binomial. However, here the overdispersion parameter theta is not specified by the user and always estimated (really the reciprocal of the dispersion parameter is estimated). A call to this function can be passed to the family argument of stan_glm or stan_glmer to estimate a Negative Binomial model. Alternatively, the stan_glm.nb and stan_glmer.nb wrapper functions may be used, which call neg_binomial_2 internally.

Usage

neg_binomial_2(link = "log")

Arguments

Value

An object of class family very similar to that of poisson but with a different family name.

Examples

if (.Platform$OS.type != "windows")
stan_glm(Days ~ Sex/(Age + Eth*Lrn), data = MASS::quine, seed = 123,
         family = neg_binomial_2, QR = TRUE, algorithm = "optimizing")

# or, equivalently, call stan_glm.nb() without specifying the family

Description

The methods documented on this page are actually some of the least important methods defined for stanreg objects. The most important methods are documented separately, each with its own page. Links to those pages are provided in the See Also section, below.

Usage

## S3 method for class 'stanmvreg'
nobs(object, ...)

## S3 method for class 'stanreg'
coef(object, ...)

## S3 method for class 'stanreg'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'stanreg'
fitted(object, ...)

## S3 method for class 'stanreg'
nobs(object, ...)

## S3 method for class 'stanreg'
residuals(object, ...)

## S3 method for class 'stanreg'
se(object, ...)

## S3 method for class 'stanreg'
update(object, formula., ..., evaluate = TRUE)

## S3 method for class 'stanreg'
vcov(object, correlation = FALSE, ...)

## S3 method for class 'stanreg'
fixef(object, ...)

## S3 method for class 'stanreg'
ngrps(object, ...)

## S3 method for class 'stanreg'
nsamples(object, ...)

## S3 method for class 'stanreg'
ranef(object, ...)

## S3 method for class 'stanreg'
sigma(object, ...)

## S3 method for class 'stanreg'
VarCorr(x, sigma = 1, ...)

Arguments

Details

The methods documented on this page are similar to the methods defined for objects of class 'lm', 'glm', 'glmer', etc. However there are a few key differences:

residuals

Residuals are always of type "response" (not "deviance" residuals or any other type). However, in the case of stan_polr with more than two response categories, the residuals are the difference between the latent utility and its linear predictor.

coef

Medians are used for point estimates. See the Point estimates section in print.stanreg for more details.

se

The se function returns standard errors based on mad. See the Uncertainty estimates section in print.stanreg for more details.

confint

For models fit using optimization, confidence intervals are returned via a call to confint.default. If algorithm is "sampling", "meanfield", or "fullrank", the confint will throw an error because the posterior_interval function should be used to compute Bayesian uncertainty intervals.

nsamples

The number of draws from the posterior distribution obtained

See Also

• The print, summary, and prior_summary methods for stanreg objects for information on the fitted model.

• launch_shinystan to use the ShinyStan GUI to explore a fitted rstanarm model.

• The plot method to plot estimates and diagnostics.

• The pp_check method for graphical posterior predictive checking.

• The posterior_predict and predictive_error methods for predictions and predictive errors.

• The posterior_interval and predictive_interval methods for uncertainty intervals for model parameters and predictions.

• The loo, kfold, and log_lik methods for leave-one-out or K-fold cross-validation, model comparison, and computing the log-likelihood of (possibly new) data.

• The as.matrix, as.data.frame, and as.array methods to access posterior draws.

Description

Interface to bayesplot's mcmc_pairs function for use with rstanarm models. Be careful not to specify too many parameters to include or the plot will be both hard to read and slow to render.

Usage

## S3 method for class 'stanreg'
pairs(
  x,
  pars = NULL,
  regex_pars = NULL,
  condition = pairs_condition(nuts = "accept_stat__"),
  ...
)

Arguments

Examples

if (.Platform$OS.type != "windows") {

if (!exists("example_model")) example(example_model)

bayesplot::color_scheme_set("purple")

# see 'condition' argument above for details on the plots below and
# above the diagonal. default is to split by accept_stat__.
pairs(example_model, pars = c("(Intercept)", "log-posterior"))

# for demonstration purposes, intentionally fit a model that
# will (almost certainly) have some divergences
fit <- stan_glm(
  mpg ~ ., data = mtcars,
  iter = 1000,
  # this combo of prior and adapt_delta should lead to some divergences
  prior = hs(),
  adapt_delta = 0.9,
  refresh = 0
)

pairs(fit, pars = c("wt", "sigma", "log-posterior"))

# requires hexbin package
# pairs(
# fit,
#   pars = c("wt", "sigma", "log-posterior"),
#   transformations = list(sigma = "log"), # show log(sigma) instead of sigma
#  off_diag_fun = "hex" # use hexagonal heatmaps instead of scatterplots
# )

bayesplot::color_scheme_set("brightblue")
pairs(
  fit,
  pars = c("(Intercept)", "wt", "sigma", "log-posterior"),
  transformations = list(sigma = "log"),
  off_diag_args = list(size = 3/4, alpha = 1/3), # size and transparency of scatterplot points
  np_style = pairs_style_np(div_color = "black", div_shape = 2) # color and shape of the divergences
)

# Using the condition argument to show divergences above the diagonal
pairs(
  fit,
  pars = c("(Intercept)", "wt", "log-posterior"),
  condition = pairs_condition(nuts = "divergent__")
)


}

Description

This generic plot method for predict.stanjm objects will plot the estimated subject-specific or marginal longitudinal trajectory using the data frame returned by a call to posterior_traj. To ensure that enough data points are available to plot the longitudinal trajectory, it is assumed that the call to posterior_traj would have used the default interpolate = TRUE, and perhaps also extrapolate = TRUE (the latter being optional, depending on whether or not the user wants to see extrapolation of the longitudinal trajectory beyond the last observation time).

Usage

## S3 method for class 'predict.stanjm'
plot(
  x,
  ids = NULL,
  limits = c("ci", "pi", "none"),
  xlab = NULL,
  ylab = NULL,
  vline = FALSE,
  plot_observed = FALSE,
  facet_scales = "free_x",
  ci_geom_args = NULL,
  grp_overlay = FALSE,
  ...
)

Arguments

Value

A ggplot object, also of class plot.predict.stanjm. This object can be further customised using the ggplot2 package. It can also be passed to the function plot_stack_jm.

See Also

posterior_traj, plot_stack_jm, posterior_survfit, plot.survfit.stanjm

Examples

# Run example model if not already loaded
  if (!exists("example_jm")) example(example_jm)
  
  # For a subset of individuals in the estimation dataset we will
  # obtain subject-specific predictions for the longitudinal submodel 
  # at evenly spaced times between 0 and their event or censoring time.
  pt1 <- posterior_traj(example_jm, ids = c(7,13,15), interpolate = TRUE)
  plot(pt1)                  # credible interval for mean response
  plot(pt1, limits = "pi")   # prediction interval for raw response
  plot(pt1, limits = "none") # no uncertainty interval
  
  # We can also extrapolate the longitudinal trajectories.
  pt2 <- posterior_traj(example_jm, ids = c(7,13,15), interpolate = TRUE,
                        extrapolate = TRUE)
  plot(pt2)
  plot(pt2, vline = TRUE)    # add line indicating event or censoring time
  plot(pt2, vline = TRUE, plot_observed = TRUE)  # overlay observed longitudinal data
 
  # We can change or add attributes to the plot
  plot1 <- plot(pt2, ids = c(7,13,15), xlab = "Follow up time",
                     vline = TRUE, plot_observed = TRUE, 
                     facet_scales = "fixed", color = "blue", linetype = 2,
                     ci_geom_args = list(fill = "red"))
  plot1
       
  # Since the returned plot is also a ggplot object, we can
  # modify some of its attributes after it has been returned
  plot1 + 
    ggplot2::theme(strip.background = ggplot2::element_blank()) +
    ggplot2::labs(title = "Some plotted longitudinal trajectories")

Description

The plot method for stanreg-objects provides a convenient interface to the MCMC module in the bayesplot package for plotting MCMC draws and diagnostics. It is also straightforward to use the functions from the bayesplot package directly rather than via the plot method. Examples of both methods of plotting are given below.

Usage

## S3 method for class 'stanreg'
plot(x, plotfun = "intervals", pars = NULL, regex_pars = NULL, ...)

## S3 method for class 'stansurv'
plot(
  x,
  plotfun = "basehaz",
  pars = NULL,
  regex_pars = NULL,
  ...,
  prob = 0.95,
  limits = c("ci", "none"),
  ci_geom_args = NULL,
  n = 1000
)

Arguments

Value

Either a ggplot object that can be further customized using the ggplot2 package, or an object created from multiple ggplot objects (e.g. a gtable object created by arrangeGrob).

References

Gabry, J. , Simpson, D. , Vehtari, A. , Betancourt, M. and Gelman, A. (2019), Visualization in Bayesian workflow. J. R. Stat. Soc. A, 182: 389-402. doi:10.1111/rssa.12378, arXiv preprint, code on GitHub)

See Also

• The vignettes in the bayesplot package for many examples.

• MCMC-overview (bayesplot) for links to the documentation for all the available plotting functions.

• color_scheme_set (bayesplot) to change the color scheme used for plotting.

• pp_check for graphical posterior predictive checks.

• plot_nonlinear for models with nonlinear smooth functions fit using stan_gamm4.

Examples

if (.Platform$OS.type != "windows") {

# Use rstanarm example model
if (!exists("example_model")) example(example_model)
fit <- example_model

#####################################
### Intervals and point estimates ###
#####################################
plot(fit) # same as plot(fit, "intervals"), plot(fit, "mcmc_intervals")

p <- plot(fit, pars = "size", regex_pars = "period",
          prob = 0.5, prob_outer = 0.9)
p + ggplot2::ggtitle("Posterior medians \n with 50% and 90% intervals")

# Shaded areas under densities
bayesplot::color_scheme_set("brightblue")
plot(fit, "areas", regex_pars = "period",
     prob = 0.5, prob_outer = 0.9)

# Make the same plot by extracting posterior draws and calling
# bayesplot::mcmc_areas directly
x <- as.array(fit, regex_pars = "period")
bayesplot::mcmc_areas(x, prob = 0.5, prob_outer = 0.9)

# Ridgelines version of the areas plot
bayesplot::mcmc_areas_ridges(x, regex_pars = "period", prob = 0.9)


##################################
### Histograms & density plots ###
##################################
plot_title <- ggplot2::ggtitle("Posterior Distributions")
plot(fit, "hist", regex_pars = "period") + plot_title
plot(fit, "dens_overlay", pars = "(Intercept)",
     regex_pars = "period") + plot_title

####################
### Scatterplots ###
####################
bayesplot::color_scheme_set("teal")
plot(fit, "scatter", pars = paste0("period", 2:3))
plot(fit, "scatter", pars = c("(Intercept)", "size"),
     size = 3, alpha = 0.5) +
     ggplot2::stat_ellipse(level = 0.9)


####################################################
### Rhat, effective sample size, autocorrelation ###
####################################################
bayesplot::color_scheme_set("red")

# rhat
plot(fit, "rhat")
plot(fit, "rhat_hist")

# ratio of effective sample size to total posterior sample size
plot(fit, "neff")
plot(fit, "neff_hist")

# autocorrelation by chain
plot(fit, "acf", pars = "(Intercept)", regex_pars = "period")
plot(fit, "acf_bar", pars = "(Intercept)", regex_pars = "period")


##################
### Traceplots ###
##################
# NOTE: rstanarm doesn't store the warmup draws (to save space because they
# are not so essential for diagnosing the particular models implemented in
# rstanarm) so the iterations in the traceplot are post-warmup iterations

bayesplot::color_scheme_set("pink")
(trace <- plot(fit, "trace", pars = "(Intercept)"))

# change traceplot colors to ggplot defaults or custom values
trace + ggplot2::scale_color_discrete()
trace + ggplot2::scale_color_manual(values = c("maroon", "skyblue2"))

# changing facet layout
plot(fit, "trace", pars = c("(Intercept)", "period2"),
     facet_args = list(nrow = 2))
# same plot by calling bayesplot::mcmc_trace directly
x <- as.array(fit, pars = c("(Intercept)", "period2"))
bayesplot::mcmc_trace(x, facet_args = list(nrow = 2))


############
### More ###
############

# regex_pars examples
plot(fit, regex_pars = "herd:1\\]")
plot(fit, regex_pars = "herd:[279]")
plot(fit, regex_pars = "herd:[279]|period2")
plot(fit, regex_pars = c("herd:[279]", "period2"))


# For graphical posterior predictive checks see
# help("pp_check.stanreg")
}

Description

This generic plot method for survfit.stanjm objects will plot the estimated subject-specific or marginal survival function using the data frame returned by a call to posterior_survfit. The call to posterior_survfit should ideally have included an "extrapolation" of the survival function, obtained by setting the extrapolate argument to TRUE.

The plot_stack_jm function takes arguments containing the plots of the estimated subject-specific longitudinal trajectory (or trajectories if a multivariate joint model was estimated) and the plot of the estimated subject-specific survival function and combines them into a single figure. This is most easily understood by running the Examples below.

Usage

## S3 method for class 'survfit.stanjm'
plot(
  x,
  ids = NULL,
  limits = c("ci", "none"),
  xlab = NULL,
  ylab = NULL,
  facet_scales = "free",
  ci_geom_args = NULL,
  ...
)

## S3 method for class 'survfit.stansurv'
plot(
  x,
  ids = NULL,
  limits = c("ci", "none"),
  xlab = NULL,
  ylab = NULL,
  facet_scales = "free",
  ci_geom_args = NULL,
  ...
)

plot_stack_jm(yplot, survplot)

Arguments

Value

The plot method returns a ggplot object, also of class plot.survfit.stanjm. This object can be further customised using the ggplot2 package. It can also be passed to the function plot_stack_jm.

plot_stack_jm returns an object of class bayesplot_grid that includes plots of the estimated subject-specific longitudinal trajectories stacked on top of the associated subject-specific survival curve.

See Also

posterior_survfit, plot_stack_jm, posterior_traj, plot.predict.stanjm

plot.predict.stanjm, plot.survfit.stanjm, posterior_predict, posterior_survfit

Examples

# Run example model if not already loaded
  if (!exists("example_jm")) example(example_jm)

  # Obtain subject-specific conditional survival probabilities
  # for all individuals in the estimation dataset.
  ps1 <- posterior_survfit(example_jm, extrapolate = TRUE)

  # We then plot the conditional survival probabilities for
  # a subset of individuals
  plot(ps1, ids = c(7,13,15))
  # We can change or add attributes to the plot
  plot(ps1, ids = c(7,13,15), limits = "none")
  plot(ps1, ids = c(7,13,15), xlab = "Follow up time")
  plot(ps1, ids = c(7,13,15), ci_geom_args = list(fill = "red"),
       color = "blue", linetype = 2)
  plot(ps1, ids = c(7,13,15), facet_scales = "fixed")

  # Since the returned plot is also a ggplot object, we can
  # modify some of its attributes after it has been returned
  plot1 <- plot(ps1, ids = c(7,13,15))
  plot1 +
    ggplot2::theme(strip.background = ggplot2::element_blank()) +
    ggplot2::coord_cartesian(xlim = c(0, 15)) +
    ggplot2::labs(title = "Some plotted survival functions")

  # We can also combine the plot(s) of the estimated
  # subject-specific survival functions, with plot(s)
  # of the estimated longitudinal trajectories for the
  # same individuals
  ps1 <- posterior_survfit(example_jm, ids = c(7,13,15))
  pt1 <- posterior_traj(example_jm, , ids = c(7,13,15))
  plot_surv <- plot(ps1)
  plot_traj <- plot(pt1, vline = TRUE, plot_observed = TRUE)
  plot_stack_jm(plot_traj, plot_surv)

  # Lastly, let us plot the standardised survival function
  # based on all individuals in our estimation dataset
  ps2 <- posterior_survfit(example_jm, standardise = TRUE, times = 0,
                          control = list(epoints = 20))
  plot(ps2)

if (.Platform$OS.type != "windows") {

  if (!exists("example_jm")) example(example_jm)
  ps1 <- posterior_survfit(example_jm, ids = c(7,13,15))
  pt1 <- posterior_traj(example_jm, ids = c(7,13,15), extrapolate = TRUE)
  plot_surv <- plot(ps1)
  plot_traj <- plot(pt1, vline = TRUE, plot_observed = TRUE)
  plot_stack_jm(plot_traj, plot_surv)

}

Description

For models fit using MCMC (algorithm="sampling") or one of the variational approximations ("meanfield" or "fullrank"), the posterior_interval function computes Bayesian posterior uncertainty intervals. These intervals are often referred to as credible intervals, but we use the term uncertainty intervals to highlight the fact that wider intervals correspond to greater uncertainty.

Usage

## S3 method for class 'stanreg'
posterior_interval(
  object,
  prob = 0.9,
  type = "central",
  pars = NULL,
  regex_pars = NULL,
  ...
)

Arguments

Details

Interpretation

Unlike for a frenquentist confidence interval, it is valid to say that, conditional on the data and model, we believe that with probability $p$ the value of a parameter is in its $100p$% posterior interval. This intuitive interpretation of Bayesian intervals is often erroneously applied to frequentist confidence intervals. See Morey et al. (2015) for more details on this issue and the advantages of using Bayesian posterior uncertainty intervals (also known as credible intervals).

Default 90% intervals

We default to reporting $90$% intervals rather than $95$% intervals for several reasons:

• Computational stability: $90$% intervals are more stable than $95$% intervals (for which each end relies on only $2.5$% of the posterior draws).

• Relation to Type-S errors (Gelman and Carlin, 2014): $95$% of the mass in a $90$% central interval is above the lower value (and $95$% is below the upper value). For a parameter $\theta$, it is therefore easy to see if the posterior probability that $\theta > 0$ (or $\theta < 0$) is larger or smaller than $95$%.

Of course, if $95$% intervals are desired they can be computed by specifying prob=0.95.

Types of intervals

Currently posterior_interval only computes central intervals because other types of intervals are rarely useful for the models that rstanarm can estimate. Additional possibilities may be provided in future releases as more models become available.

Value

A matrix with two columns and as many rows as model parameters (or the subset of parameters specified by pars and/or regex_pars). For a given value of prob, $p$, the columns correspond to the lower and upper $100p$% interval limits and have the names $100\alpha/2$% and $100(1 - \alpha/2)$%, where $\alpha = 1-p$. For example, if prob=0.9 is specified (a $90$% interval), then the column names will be "5%" and "95%", respectively.

References

Gelman, A. and Carlin, J. (2014). Beyond power calculations: assessing Type S (sign) and Type M (magnitude) errors. Perspectives on Psychological Science. 9(6), 641–51.

Morey, R. D., Hoekstra, R., Rouder, J., Lee, M. D., and Wagenmakers, E. (2016). The fallacy of placing confidence in confidence intervals. Psychonomic Bulletin & Review. 23(1), 103–123.

See Also

confint.stanreg, which, for models fit using optimization, can be used to compute traditional confidence intervals.

predictive_interval for predictive intervals.

Examples

if (.Platform$OS.type != "windows") {
if (!exists("example_model")) example(example_model)
posterior_interval(example_model)
posterior_interval(example_model, regex_pars = "herd")
posterior_interval(example_model, pars = "period2", prob = 0.5)
}

Description

Extract the posterior draws of the linear predictor, possibly transformed by the inverse-link function. This function is occasionally useful, but it should be used sparingly: inference and model checking should generally be carried out using the posterior predictive distribution (i.e., using posterior_predict).

Usage

## S3 method for class 'stanreg'
posterior_linpred(
  object,
  transform = FALSE,
  newdata = NULL,
  draws = NULL,
  re.form = NULL,
  offset = NULL,
  XZ = FALSE,
  ...
)

## S3 method for class 'stanreg'
posterior_epred(
  object,
  newdata = NULL,
  draws = NULL,
  re.form = NULL,
  offset = NULL,
  XZ = FALSE,
  ...
)

Arguments

Details

The posterior_linpred function returns the posterior distribution of the linear predictor, while the posterior_epred function returns the posterior distribution of the conditional expectation. In the special case of a Gaussian likelihood with an identity link function, these two concepts are the same. The posterior_epred function is a less noisy way to obtain expectations over the output of posterior_predict.

Value

The default is to return a draws by nrow(newdata) matrix of simulations from the posterior distribution of the (possibly transformed) linear predictor. The exception is if the argument XZ is set to TRUE (see the XZ argument description above).

Note

For models estimated with stan_clogit, the number of successes per stratum is ostensibly fixed by the research design. Thus, when calling posterior_linpred with new data and transform = TRUE, the data.frame passed to the newdata argument must contain an outcome variable and a stratifying factor, both with the same name as in the original data.frame. Then, the probabilities will condition on this outcome in the new data.

See Also

posterior_predict to draw from the posterior predictive distribution of the outcome, which is typically preferable.

Examples

if (.Platform$OS.type != "windows") {
if (!exists("example_model")) example(example_model)
print(family(example_model))

# linear predictor on log-odds scale
linpred <- posterior_linpred(example_model)
colMeans(linpred)

# probabilities
# same as posterior_linpred(example_model, transform = TRUE)
probs <- posterior_epred(example_model)
colMeans(probs)

# not conditioning on any group-level parameters
probs2 <- posterior_epred(example_model, re.form = NA)
apply(probs2, 2, median)
}

Description

The posterior predictive distribution is the distribution of the outcome implied by the model after using the observed data to update our beliefs about the unknown parameters in the model. Simulating data from the posterior predictive distribution using the observed predictors is useful for checking the fit of the model. Drawing from the posterior predictive distribution at interesting values of the predictors also lets us visualize how a manipulation of a predictor affects (a function of) the outcome(s). With new observations of predictor variables we can use the posterior predictive distribution to generate predicted outcomes.

Usage

## S3 method for class 'stanreg'
posterior_predict(
  object,
  newdata = NULL,
  draws = NULL,
  re.form = NULL,
  fun = NULL,
  seed = NULL,
  offset = NULL,
  ...
)

## S3 method for class 'stanmvreg'
posterior_predict(
  object,
  m = 1,
  newdata = NULL,
  draws = NULL,
  re.form = NULL,
  fun = NULL,
  seed = NULL,
  offset = NULL,
  ...
)

Arguments

Value

A draws by nrow(newdata) matrix of simulations from the posterior predictive distribution. Each row of the matrix is a vector of predictions generated using a single draw of the model parameters from the posterior distribution.

Note

For binomial models with a number of trials greater than one (i.e., not Bernoulli models), if newdata is specified then it must include all variables needed for computing the number of binomial trials to use for the predictions. For example if the left-hand side of the model formula is cbind(successes, failures) then both successes and failures must be in newdata. The particular values of successes and failures in newdata do not matter so long as their sum is the desired number of trials. If the left-hand side of the model formula were cbind(successes, trials - successes) then both trials and successes would need to be in newdata, probably with successes set to 0 and trials specifying the number of trials. See the Examples section below and the How to Use the rstanarm Package for examples.

For models estimated with stan_clogit, the number of successes per stratum is ostensibly fixed by the research design. Thus, when doing posterior prediction with new data, the data.frame passed to the newdata argument must contain an outcome variable and a stratifying factor, both with the same name as in the original data.frame. Then, the posterior predictions will condition on this outcome in the new data.

See Also

pp_check for graphical posterior predictive checks. Examples of posterior predictive checking can also be found in the rstanarm vignettes and demos.

predictive_error and predictive_interval.

Examples

if (.Platform$OS.type != "windows") {
if (!exists("example_model")) example(example_model)
yrep <- posterior_predict(example_model)
table(yrep)


# Using newdata
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
dat <- data.frame(counts, treatment, outcome)
fit3 <- stan_glm(
  counts ~ outcome + treatment,
  data = dat,
  family = poisson(link="log"),
  prior = normal(0, 1, autoscale = FALSE),
  prior_intercept = normal(0, 5, autoscale = FALSE),
  refresh = 0
)
nd <- data.frame(treatment = factor(rep(1,3)), outcome = factor(1:3))
ytilde <- posterior_predict(fit3, nd, draws = 500)
print(dim(ytilde))  # 500 by 3 matrix (draws by nrow(nd))

ytilde <- data.frame(
  count = c(ytilde),
  outcome = rep(nd$outcome, each = 500)
)
ggplot2::ggplot(ytilde, ggplot2::aes(x=outcome, y=count)) +
  ggplot2::geom_boxplot() +
  ggplot2::ylab("predicted count")


# Using newdata with a binomial model.
# example_model is binomial so we need to set
# the number of trials to use for prediction.
# This could be a different number for each
# row of newdata or the same for all rows.
# Here we'll use the same value for all.
nd <- lme4::cbpp
print(formula(example_model))  # cbind(incidence, size - incidence) ~ ...
nd$size <- max(nd$size) + 1L   # number of trials
nd$incidence <- 0  # set to 0 so size - incidence = number of trials
ytilde <- posterior_predict(example_model, newdata = nd)


# Using fun argument to transform predictions
mtcars2 <- mtcars
mtcars2$log_mpg <- log(mtcars2$mpg)
fit <- stan_glm(log_mpg ~ wt, data = mtcars2, refresh = 0)
ytilde <- posterior_predict(fit, fun = exp)

}

Description

This function allows us to generate predicted quantities for survival models at specified times. These quantities include the hazard rate, cumulative hazard, survival probability, or failure probability (i.e. CDF). Note that the cumulative hazard, survival probability, or failure probability may be conditional on a last known survival time (see the condition argument discussed below). Predictions are obtained using unique draws from the posterior distribution of each of the model parameters and then summarised into a median and posterior uncertainty interval. For stan_jm models "dynamic" predictions are allowed and are in fact the default when new data is provided (see the dynamic argument discussed below).

Usage

posterior_survfit(object, ...)

## S3 method for class 'stansurv'
posterior_survfit(
  object,
  newdata = NULL,
  type = "surv",
  extrapolate = TRUE,
  control = list(),
  condition = FALSE,
  last_time = NULL,
  prob = 0.95,
  times = NULL,
  standardise = FALSE,
  draws = NULL,
  seed = NULL,
  return_matrix = FALSE,
  ...
)

## S3 method for class 'stanjm'
posterior_survfit(
  object,
  newdataLong = NULL,
  newdataEvent = NULL,
  type = "surv",
  extrapolate = TRUE,
  control = list(),
  condition = NULL,
  last_time = NULL,
  prob = 0.95,
  ids,
  times = NULL,
  standardise = FALSE,
  dynamic = TRUE,
  scale = 1.5,
  draws = NULL,
  seed = NULL,
  return_matrix = FALSE,
  ...
)

Arguments

Details

By default, the predicted quantities are evaluated conditional on observed values of the fixed effect covariates. That is, predictions will be obtained using either:

• the design matrices used in the original stan_surv or stan_jm model call, or

• the covariate values provided in the newdata argument (or newdataLong and newdataEvent arugments for the stanjm method).

However, if you wish to average over the observed distribution of the fixed effect covariates then this is possible – such predictions are sometimes referred to as standardised survival probabilties – see the standardise argument.

For stansurv objects, the predicted quantities are calculated for each row of the prediction data, at the specified times as well as any times generated through extrapolation (when extrapolate = TRUE).

For stanjm objects, the predicted quantities are calculated for each individual, at the specified times as well as any times generated through extrapolation (when extrapolate = TRUE).

Dynamic versus marginalised predictions

The following also applies for stanjm objects. By default the survival probabilities are conditional on an individual's group-specific coefficients (i.e. their individual-level random effects). If prediction data is provided via the newdataLong and newdataEvent arguments, then the default behaviour is to sample new group-specific coefficients for the individuals in the new data using a Monte Carlo scheme that conditions on their longitudinal outcome data provided in newdataLong (sometimes referred to as "dynamic predictions", see Rizopoulos (2011)). This default behaviour can be stopped by specifying dynamic = FALSE, in which case the predicted survival probabilities will be marginalised over the distribution of the group-specific coefficients. This has the benefit that the user does not need to provide longitudinal outcome measurements for the new individuals, however, it does mean that the predictions will incorporate all the uncertainty associated with between-individual variation in the biomarker (longitudinal outcome) values since the predictions aren't conditional on any observed biomarker (longitudinal outcome) data for the individual.

Value

When return_matrix = FALSE (the default), a data frame of class survfit.stansurv or survfit.stanjm. The data frame includes columns for each of the following: (i) the median of the posterior predictions (median); (ii) each of the lower and upper limits of the corresponding uncertainty interval for the posterior predictions (ci_lb and ci_ub); (iii) an observation identifier (for stan_surv models) or an individual identifier (for stan_jm models), unless standardised predictions were requested; (iv) the time that the prediction corresponds to (time). (v) the last known survival time on which the prediction is conditional (cond_time); this will be set to NA if not relevant. The returned object also includes a number of additional attributes.

When return_matrix = TRUE a list of matrices is returned. Each matrix contains the predictions evaluated at one step of the extrapolation time sequence (note that if extrapolate = FALSE then the list will be of length one, i.e. the predictions are only evaluated at times which corresponds to just one time point for each individual). Each matrix will have draws rows and nrow(newdata) columns, such that each row contains a vector of predictions generated using a single draw of the model parameters from the posterior distribution. The returned list also includes a number of additional attributes.

Note

Note that if any variables were transformed (e.g. rescaled) in the data used to fit the model, then these variables must also be transformed in newdataLong and newdataEvent. This only applies if variables were transformed before passing the data to one of the modeling functions and not if transformations were specified inside the model formula.

References

Rizopoulos, D. (2011). Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data. Biometrics 67, 819.

See Also

plot.survfit.stanjm for plotting the estimated survival probabilities
ps_check for for graphical checks of the estimated survival function
posterior_traj for estimating the marginal or subject-specific longitudinal trajectories
plot_stack_jm for combining plots of the estimated subject-specific longitudinal trajectory and survival function

Examples

# Run example model if not already loaded
  if (!exists("example_jm")) example(example_jm)

  # Obtain subject-specific survival probabilities for a few
  # selected individuals in the estimation dataset who were
  # known to survive up until their censoring time. By default
  # the posterior_survfit function will estimate the conditional
  # survival probabilities, that is, conditional on having survived
  # until the event or censoring time, and then by default will
  # extrapolate the survival predictions forward from there.
  ps1 <- posterior_survfit(example_jm, ids = c(7,13,15))
  # We can plot the estimated survival probabilities using the
  # associated plot function
  plot(ps1)

  # If we wanted to estimate the survival probabilities for the
  # same three individuals as the previous example, but this time
  # we won't condition on them having survived up until their
  # censoring time. Instead, we will estimate their probability
  # of having survived between 0 and 5 years given their covariates
  # and their estimated random effects.
  # The easiest way to achieve the time scale we want (ie, 0 to 5 years)
  # is to specify that we want the survival time estimated at time 0
  # and then extrapolated forward 5 years. We also specify that we
  # do not want to condition on their last known survival time.
  ps2 <- posterior_survfit(example_jm, ids = c(7,13,15), times = 0,
    extrapolate = TRUE, condition = FALSE, control = list(edist = 5))

  # Instead we may want to estimate subject-specific survival probabilities
  # for a set of new individuals. To demonstrate this, we will simply take
  # the first two individuals in the estimation dataset, but pass their data
  # via the newdata arguments so that posterior_survfit will assume we are
  # predicting survival for new individuals and draw new random effects
  # under a Monte Carlo scheme (see Rizopoulos (2011)).
  ndL <- pbcLong[pbcLong$id %in% c(1,2),]
  ndE <- pbcSurv[pbcSurv$id %in% c(1,2),]
  ps3 <- posterior_survfit(example_jm,
    newdataLong = ndL, newdataEvent = ndE,
    last_time = "futimeYears", seed = 12345)
  head(ps3)
  # We can then compare the estimated random effects for these
  # individuals based on the fitted model and the Monte Carlo scheme
  ranef(example_jm)$Long1$id[1:2,,drop=FALSE] # from fitted model
  colMeans(attr(ps3, "b_new"))                # from Monte Carlo scheme

  # Lastly, if we wanted to obtain "standardised" survival probabilities,
  # (by averaging over the observed distribution of the fixed effect
  # covariates, as well as averaging over the estimated random effects
  # for individuals in our estimation sample or new data) then we can
  # specify 'standardise = TRUE'. We can then plot the resulting
  # standardised survival curve.
  ps4 <- posterior_survfit(example_jm, standardise = TRUE,
                           times = 0, extrapolate = TRUE)
  plot(ps4)

Description

This function allows us to generate an estimated longitudinal trajectory (either subject-specific, or by marginalising over the distribution of the group-specific parameters) based on draws from the posterior predictive distribution.

Usage

posterior_traj(
  object,
  m = 1,
  newdata = NULL,
  newdataLong = NULL,
  newdataEvent = NULL,
  interpolate = TRUE,
  extrapolate = FALSE,
  control = list(),
  last_time = NULL,
  prob = 0.95,
  ids,
  dynamic = TRUE,
  scale = 1.5,
  draws = NULL,
  seed = NULL,
  return_matrix = FALSE,
  ...
)

Arguments

Details

The posterior_traj function acts as a wrapper to the posterior_predict function, but allows predictions to be easily generated at time points that are interpolated and/or extrapolated between time zero (baseline) and the last known survival time for the individual, thereby providing predictions that correspond to a smooth estimate of the longitudinal trajectory (useful for the plotting via the associated plot.predict.stanjm method). In addition it returns a data frame by default, whereas the posterior_predict function returns a matrix; see the Value section below for details. Also, posterior_traj allows predictions to only be generated for a subset of individuals, via the ids argument.

Value

When return_matrix = FALSE, a data frame of class predict.stanjm. The data frame includes a column for the median of the posterior predictions of the mean longitudinal response (yfit), a column for each of the lower and upper limits of the uncertainty interval corresponding to the posterior predictions of the mean longitudinal response (ci_lb and ci_ub), and a column for each of the lower and upper limits of the prediction interval corresponding to the posterior predictions of the (raw) longitudinal response. The data frame also includes columns for the subject ID variable, and each of the predictor variables. The returned object also includes a number of attributes.

When return_matrix = TRUE, the returned object is the same as that described for posterior_predict.

See Also

plot.predict.stanjm, posterior_predict, posterior_survfit

Examples

# Run example model if not already loaded
  if (!exists("example_jm")) example(example_jm)
  
  # Obtain subject-specific predictions for all individuals 
  # in the estimation dataset
  pt1 <- posterior_traj(example_jm, interpolate = FALSE, extrapolate = FALSE)
  head(pt1)
  
  # Obtain subject-specific predictions only for a few selected individuals
  pt2 <- posterior_traj(example_jm, ids = c(1,3,8))
  
  # If we wanted to obtain subject-specific predictions in order to plot the 
  # longitudinal trajectories, then we might want to ensure a full trajectory 
  # is obtained by interpolating and extrapolating time. We can then use the 
  # generic plot function to plot the subject-specific predicted trajectories
  # for the first three individuals. Interpolation and extrapolation is 
  # carried out by default.
  pt3 <- posterior_traj(example_jm)
  head(pt3) # predictions at additional time points compared with pt1 
  plot(pt3, ids = 1:3)
  
  # If we wanted to extrapolate further in time, but decrease the number of 
  # discrete time points at which we obtain predictions for each individual, 
  # then we could specify a named list in the 'control' argument
  pt4 <- posterior_traj(example_jm, control = list(ipoints = 10, epoints = 10, eprop = 0.5))
  
  # If we have prediction data for a new individual, and we want to
  # estimate the longitudinal trajectory for that individual conditional
  # on this new data (perhaps extrapolating forward from our last
  # longitudinal measurement) then we can do that. It requires drawing
  # new individual-specific parameters, based on the full likelihood,
  # so we must supply new data for both the longitudinal and event 
  # submodels. These are sometimes known as dynamic predictions.
  ndL <- pbcLong[pbcLong$id == 8, , drop = FALSE]
  ndE <- pbcSurv[pbcSurv$id == 8, , drop = FALSE]
  ndL$id <- "new_subject" # new id can't match one used in training data
  ndE$id <- "new_subject"
  pt5 <- posterior_traj(example_jm, 
                        newdataLong = ndL,
                        newdataEvent = ndE)
                        
  # By default it is assumed that the last known survival time for 
  # the individual is the time of their last biomarker measurement,
  # but if we know they survived to some later time then we can
  # condition on that information using the last_time argument
  pt6 <- posterior_traj(example_jm, 
                        newdataLong = ndL,
                        newdataEvent = ndE, 
                        last_time = "futimeYears")
  
  # Alternatively we may want to estimate the marginal longitudinal
  # trajectory for a given set of covariates. To do this, we can pass
  # the desired covariate values in a new data frame (however the only
  # covariate in our fitted model was the time variable, year). To make sure  
  # that we marginalise over the random effects, we need to specify an ID value
  # which does not correspond to any of the individuals who were used in the
  # model estimation and specify the argument dynamic=FALSE.
  # The marginal prediction is obtained by generating subject-specific 
  # predictions using a series of random draws from the random 
  # effects distribution, and then integrating (ie, averaging) over these. 
  # Our marginal prediction will therefore capture the between-individual 
  # variation associated with the random effects.
  
  nd <- data.frame(id = rep("new1", 11), year = (0:10 / 2))
  pt7 <- posterior_traj(example_jm, newdataLong = nd, dynamic = FALSE)
  head(pt7)  # note the greater width of the uncertainty interval compared 
             # with the subject-specific predictions in pt1, pt2, etc
  
  # Alternatively, we could have estimated the "marginal" trajectory by 
  # ignoring the random effects (ie, assuming the random effects were set 
  # to zero). This will generate a predicted longitudinal trajectory only
  # based on the fixed effect component of the model. In essence, for a 
  # linear mixed effects model (ie, a model that uses an identity link 
  # function), we should obtain a similar point estimate ("yfit") to the
  # estimates obtained in pt5 (since the mean of the estimated random effects
  # distribution will be approximately 0). However, it is important to note that
  # the uncertainty interval will be much more narrow, since it completely
  # ignores the between-individual variability captured by the random effects.
  # Further, if the model uses a non-identity link function, then the point
  # estimate ("yfit") obtained only using the fixed effect component of the
  # model will actually provide a biased estimate of the marginal prediction.
  # Nonetheless, to demonstrate how we can obtain the predictions only using 
  # the fixed effect component of the model, we simply specify 're.form = NA'. 
  # (We will use the same covariate values as used in the prediction for 
  # example for pt5).
  
  pt8 <- posterior_traj(example_jm, newdataLong = nd, dynamic = FALSE, 
                        re.form = NA)
  head(pt8)  # note the much narrower ci, compared with pt5

Description

Plot medians and central intervals comparing parameter draws from the prior and posterior distributions. If the plotted priors look different than the priors you think you specified it is likely either because of internal rescaling or the use of the QR argument (see the documentation for the prior_summary method for details on these special cases).

Usage

posterior_vs_prior(object, ...)

## S3 method for class 'stanreg'
posterior_vs_prior(
  object,
  pars = NULL,
  regex_pars = NULL,
  prob = 0.9,
  color_by = c("parameter", "vs", "none"),
  group_by_parameter = FALSE,
  facet_args = list(),
  ...
)

Arguments

Value

A ggplot object that can be further customized using the ggplot2 package.

References

Gabry, J. , Simpson, D. , Vehtari, A. , Betancourt, M. and Gelman, A. (2019), Visualization in Bayesian workflow. J. R. Stat. Soc. A, 182: 389-402. doi:10.1111/rssa.12378, arXiv preprint, code on GitHub)

Examples

if (.Platform$OS.type != "windows") {
## Not run: 
if (!exists("example_model")) example(example_model)
# display non-varying (i.e. not group-level) coefficients
posterior_vs_prior(example_model, pars = "beta")

# show group-level (varying) parameters and group by parameter
posterior_vs_prior(example_model, pars = "varying",
                   group_by_parameter = TRUE, color_by = "vs")

# group by parameter and allow axis scales to vary across facets
posterior_vs_prior(example_model, regex_pars = "period",
                   group_by_parameter = TRUE, color_by = "none",
                   facet_args = list(scales = "free"))

# assign to object and customize with functions from ggplot2
(gg <- posterior_vs_prior(example_model, pars = c("beta", "varying"), prob = 0.8))

gg +
 ggplot2::geom_hline(yintercept = 0, size = 0.3, linetype = 3) +
 ggplot2::coord_flip() +
 ggplot2::ggtitle("Comparing the prior and posterior")

# compare very wide and very narrow priors using roaches example
# (see help(roaches, "rstanarm") for info on the dataset)
roaches$roach100 <- roaches$roach1 / 100
wide_prior <- normal(0, 10)
narrow_prior <- normal(0, 0.1)
fit_pois_wide_prior <- stan_glm(y ~ treatment + roach100 + senior,
                                offset = log(exposure2),
                                family = "poisson", data = roaches,
                                prior = wide_prior)
posterior_vs_prior(fit_pois_wide_prior, pars = "beta", prob = 0.5,
                   group_by_parameter = TRUE, color_by = "vs",
                   facet_args = list(scales = "free"))

fit_pois_narrow_prior <- update(fit_pois_wide_prior, prior = narrow_prior)
posterior_vs_prior(fit_pois_narrow_prior, pars = "beta", prob = 0.5,
                   group_by_parameter = TRUE, color_by = "vs",
                   facet_args = list(scales = "free"))


# look at cutpoints for ordinal model
fit_polr <- stan_polr(tobgp ~ agegp, data = esoph, method = "probit",
                      prior = R2(0.2, "mean"), init_r = 0.1)
(gg_polr <- posterior_vs_prior(fit_polr, regex_pars = "\\|", color_by = "vs",
                               group_by_parameter = TRUE))
# flip the x and y axes
gg_polr + ggplot2::coord_flip()

## End(Not run)
}