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 = normal(),
prior_intercept = normal(), prior_aux = cauchy(0, 5),
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)
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.
The prior distribution for the regression coefficients.
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 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.
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. 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.
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.
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
adapt_delta for details.
A logical scalar defaulting to FALSE, but if TRUE
applies a scaled qr decomposition to the design matrix,
\(X = Q^\ast R^\ast\), where \(Q^\ast = Q \sqrt{n-1}\) and \(R^\ast = \frac{1}{\sqrt{n-1}} R\). The coefficients relative to \(Q^\ast\) are obtained and then
premultiplied by the inverse of \(R^{\ast}\) to obtain coefficients
relative to the original predictors, \(X\). These transformations do not
change the likelihood of the data but are recommended for computational
reasons when there are multiple predictors. Importantly, while the columns
of \(X\) are almost always correlated, the columns of \(Q^\ast\)
are uncorrelated by design, which often makes sampling from the posterior
easier. However, because when QR is TRUE the prior
argument applies to the coefficients relative to \(Q^\ast\) (and
those are not very interpretable), setting QR=TRUE is only
recommended if you do not have an informative prior for the regression
coefficients.
A logical scalar (defaulting to FALSE) indicating
whether to use a sparse representation of the design (X) matrix.
Setting this to TRUE will likely be twice as slow, even if the
design matrix has a considerable number of zeros, but it may allow the
model to be estimated when the computer has too little RAM to
utilize a dense design matrix. If TRUE, the the design matrix
is not centered (since that would destroy the sparsity) and it is
not possible to specify both QR = TRUE and sparse = TRUE.
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.
A stanreg object is returned
for stan_gamm4.
plot_nonlinear returns a ggplot object.
The stan_gamm4 function is similar in syntax to
gamm4 in the gamm4 package, which accepts a syntax
that is similar to (but not quite as extensive as) that for
gamm in the mgcv package and converts
it internally into the syntax accepted by glmer in the
lme4 package. But rather than performing (restricted) maximum likelihood
estimation, 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) and priors on the
terms of a decomposition of the covariance matrices of the group-specific
parameters, including the smooths. 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. If random = NULL, the output is a subset of that produced by
gam in the sense that there are several estimated components
for each smooth term. However, the parameterization used to estimate the model
is different and corresponds to the parameterization in
gamm4 where is smooth term is decomposed into a linear
and a non-linear part. If prior is not NULL, then the number
of parameters to place priors on is equal to the number of linear terms in
the formula. The prior on the non-linear part of each smooth term is
handled by the decov function. If random is not NULL,
then there are additional group-specific terms whose priors are also handled
by the decov function and whose posterior medians can be extracted
by calling ranef.
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.
# NOT RUN {
# from example(gamm4, package = "gamm4"), prefixing gamm4() call with stan_
# }
# NOT RUN {
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),
QR = TRUE, chains = 1, iter = 200) # for example speed
print(br)
plot_nonlinear(br)
plot_nonlinear(br, smooths = "s(x0)", alpha = 2/3)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab