Bayesian inference for GAMMs with flexible priors.
stan_gamm4(
formula,
random = NULL,
family = gaussian(),
data,
weights = NULL,
subset = NULL,
na.action,
knots = NULL,
drop.unused.levels = TRUE,
...,
prior = default_prior_coef(family),
prior_intercept = default_prior_intercept(family),
prior_smooth = exponential(autoscale = FALSE),
prior_aux = exponential(autoscale = TRUE),
prior_covariance = decov(),
prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
adapt_delta = NULL,
QR = FALSE,
sparse = FALSE
)plot_nonlinear(
x,
smooths,
...,
prob = 0.9,
facet_args = list(),
alpha = 1,
size = 0.75
)
A stanreg object is returned
for stan_gamm4
.
plot_nonlinear
returns a ggplot object.
Same as for
gamm4
. We strongly advise against
omitting the data
argument. Unless data
is specified (and is
a data frame) many post-estimation functions (including update
,
loo
, kfold
) are not guaranteed to work properly.
Same as glm
,
but rarely specified.
Further arguments passed to sampling
(e.g.
iter
, chains
, cores
, etc.) or to
vb
(if algorithm
is "meanfield"
or
"fullrank"
).
The prior distribution for the (non-hierarchical) regression coefficients.
The default priors are described in the vignette
Prior
Distributions for rstanarm Models.
If not using the default, prior
should be a call to one of the
various functions provided by rstanarm for specifying priors. The
subset of these functions that can be used for the prior on the
coefficients can be grouped into several "families":
Family | Functions |
Student t family | normal , student_t , cauchy |
Hierarchical shrinkage family | hs , hs_plus |
Laplace family | laplace , lasso |
Product normal family | product_normal |
See the priors help page for details on the families and
how to specify the arguments for all of the functions in the table above.
To omit a prior ---i.e., to use a flat (improper) uniform prior---
prior
can be set to NULL
, although this is rarely a good
idea.
Note: Unless QR=TRUE
, if prior
is from the Student t
family or Laplace family, and if the autoscale
argument to the
function used to specify the prior (e.g. normal
) is left at
its default and recommended value of TRUE
, then the default or
user-specified prior scale(s) may be adjusted internally based on the
scales of the predictors. See the priors help page and the
Prior Distributions vignette for details on the rescaling and the
prior_summary
function for a summary of the priors used for a
particular model.
The prior distribution for the intercept (after centering all predictors, see note below).
The default prior is described in the vignette
Prior
Distributions for rstanarm Models.
If not using the default, prior_intercept
can be a call to
normal
, student_t
or cauchy
. See the
priors help page for details on these functions. To omit a
prior on the intercept ---i.e., to use a flat (improper) uniform prior---
prior_intercept
can be set to NULL
.
Note: If using a dense representation of the design matrix
---i.e., if the sparse
argument is left at its default value of
FALSE
--- then the prior distribution for the intercept is set so it
applies to the value when all predictors are centered (you don't
need to manually center them). This is explained further in
[Prior Distributions for rstanarm Models](https://mc-stan.org/rstanarm/articles/priors.html)
If you prefer to specify a prior on the intercept without the predictors
being auto-centered, then you have to omit the intercept from the
formula
and include a column of ones as a predictor,
in which case some element of prior
specifies the prior on it,
rather than prior_intercept
. Regardless of how
prior_intercept
is specified, the reported estimates of the
intercept always correspond to a parameterization without centered
predictors (i.e., same as in glm
).
The prior distribution for the hyperparameters in GAMs, with lower values yielding less flexible smooth functions.
prior_smooth
can be a call to exponential
to
use an exponential distribution, or normal
, student_t
or
cauchy
, which results in a half-normal, half-t, or half-Cauchy
prior. See priors
for details on these functions. To omit a
prior ---i.e., to use a flat (improper) uniform prior--- set
prior_smooth
to NULL
. The number of hyperparameters depends
on the model specification but a scalar prior will be recylced as necessary
to the appropriate length.
The prior distribution for the "auxiliary" parameter (if
applicable). The "auxiliary" parameter refers to a different parameter
depending on the family
. For Gaussian models prior_aux
controls "sigma"
, the error
standard deviation. For negative binomial models prior_aux
controls
"reciprocal_dispersion"
, which is similar to the
"size"
parameter of rnbinom
:
smaller values of "reciprocal_dispersion"
correspond to
greater dispersion. For gamma models prior_aux
sets the prior on
to the "shape"
parameter (see e.g.,
rgamma
), and for inverse-Gaussian models it is the
so-called "lambda"
parameter (which is essentially the reciprocal of
a scale parameter). Binomial and Poisson models do not have auxiliary
parameters.
The default prior is described in the vignette
Prior
Distributions for rstanarm Models.
If not using the default, prior_aux
can be a call to
exponential
to use an exponential distribution, or normal
,
student_t
or cauchy
, which results in a half-normal, half-t,
or half-Cauchy prior. See priors
for details on these
functions. To omit a prior ---i.e., to use a flat (improper) uniform
prior--- set prior_aux
to NULL
.
Cannot be NULL
; see decov
for
more information about the default arguments.
A logical scalar (defaulting to FALSE
) indicating
whether to draw from the prior predictive distribution instead of
conditioning on the outcome.
A string (possibly abbreviated) indicating the
estimation approach to use. Can be "sampling"
for MCMC (the
default), "optimizing"
for optimization, "meanfield"
for
variational inference with independent normal distributions, or
"fullrank"
for variational inference with a multivariate normal
distribution. See rstanarm-package
for more details on the
estimation algorithms. NOTE: not all fitting functions support all four
algorithms.
Only relevant if algorithm="sampling"
. See
the adapt_delta help page for details.
A logical scalar defaulting to FALSE
, but if TRUE
applies a scaled qr
decomposition to the design matrix. The
transformation does not change the likelihood of the data but is
recommended for computational reasons when there are multiple predictors.
See the QR-argument documentation page for details on how
rstanarm does the transformation and important information about how
to interpret the prior distributions of the model parameters when using
QR=TRUE
.
A logical scalar (defaulting to FALSE
) indicating
whether to use a sparse representation of the design (X) matrix.
If TRUE
, the the design matrix is not centered (since that would
destroy the sparsity) and likewise it is not possible to specify both
QR = TRUE
and sparse = TRUE
. Depending on how many zeros
there are in the design matrix, setting sparse = TRUE
may make
the code run faster and can consume much less RAM.
An object produced by stan_gamm4
.
An optional character vector specifying a subset of the smooth
functions specified in the call to stan_gamm4
. The default is
include all smooth terms.
For univarite smooths, a scalar between 0 and 1 governing the width of the uncertainty interval.
An optional named list of arguments passed to
facet_wrap
(other than the facets
argument).
For univariate smooths, passed to
geom_ribbon
. For bivariate smooths, size/2
is
passed to geom_contour
.
The stan_gamm4
function is similar in syntax to
gamm4
in the gamm4 package. But rather than performing
(restricted) maximum likelihood estimation with the lme4 package,
the stan_gamm4
function utilizes MCMC to perform Bayesian
estimation. The Bayesian model adds priors on the common regression
coefficients (in the same way as stan_glm
), priors on the
standard deviations of the smooth terms, and a prior on the decomposition
of the covariance matrices of any group-specific parameters (as in
stan_glmer
). Estimating these models via MCMC avoids
the optimization issues that often crop up with GAMMs and provides better
estimates for the uncertainty in the parameter estimates.
See gamm4
for more information about the model
specicification and priors
for more information about the
priors on the main coefficients. The formula
should include at least
one smooth term, which can be specified in any way that is supported by the
jagam
function in the mgcv package. The
prior_smooth
argument should be used to specify a prior on the unknown
standard deviations that govern how smooth the smooth function is. The
prior_covariance
argument can be used to specify the prior on the
components of the covariance matrix for any (optional) group-specific terms.
The gamm4
function in the gamm4 package uses
group-specific terms to implement the departure from linearity in the smooth
terms, but that is not the case for stan_gamm4
where the group-specific
terms are exactly the same as in stan_glmer
.
The plot_nonlinear
function creates a ggplot object with one facet for
each smooth function specified in the call to stan_gamm4
in the case
where all smooths are univariate. A subset of the smooth functions can be
specified using the smooths
argument, which is necessary to plot a
bivariate smooth or to exclude the bivariate smooth and plot the univariate
ones. In the bivariate case, a plot is produced using
geom_contour
. In the univariate case, the resulting
plot is conceptually similar to plot.gam
except the
outer lines here demark the edges of posterior uncertainty intervals
(credible intervals) rather than confidence intervals and the inner line
is the posterior median of the function rather than the function implied
by a point estimate. To change the colors used in the plot see
color_scheme_set
.
Crainiceanu, C., Ruppert D., and Wand, M. (2005). Bayesian analysis for penalized spline regression using WinBUGS. Journal of Statistical Software. 14(14), 1--22. https://www.jstatsoft.org/article/view/v014i14
stanreg-methods
and
gamm4
.
The vignette for stan_glmer
, which also discusses
stan_gamm4
. https://mc-stan.org/rstanarm/articles/
if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") {
# from example(gamm4, package = "gamm4"), prefixing gamm4() call with stan_
# \donttest{
dat <- mgcv::gamSim(1, n = 400, scale = 2) ## simulate 4 term additive truth
## Now add 20 level random effect `fac'...
dat$fac <- fac <- as.factor(sample(1:20, 400, replace = TRUE))
dat$y <- dat$y + model.matrix(~ fac - 1) %*% rnorm(20) * .5
br <- stan_gamm4(y ~ s(x0) + x1 + s(x2), data = dat, random = ~ (1 | fac),
chains = 1, iter = 500) # for example speed
print(br)
plot_nonlinear(br)
plot_nonlinear(br, smooths = "s(x0)", alpha = 2/3)
# }
}
Run the code above in your browser using DataLab