R/stan_gamm4.R
stan_gamm4.Rd
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 )
formula, random, family, data, knots, drop.unused.levels  Same as for



subset, weights, na.action  Same as 

...  Further arguments passed to 

prior  The prior distribution for the (nonhierarchical) regression coefficients. The default priors are described in the vignette
Prior
Distributions for rstanarm Models.
If not using the default,
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
Note: Unless 

prior_intercept  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, Note: If using a dense representation of the design matrix
i.e., if the 

prior_smooth  The prior distribution for the hyperparameters in GAMs, with lower values yielding less flexible smooth functions.


prior_aux  The prior distribution for the "auxiliary" parameter (if
applicable). The "auxiliary" parameter refers to a different parameter
depending on the The default prior is described in the vignette
Prior
Distributions for rstanarm Models.
If not using the default, 

prior_covariance  Cannot be 

prior_PD  A logical scalar (defaulting to 

algorithm  A string (possibly abbreviated) indicating the
estimation approach to use. Can be 

adapt_delta  Only relevant if 

QR  A logical scalar defaulting to 

sparse  A logical scalar (defaulting to 

x  An object produced by 

smooths  An optional character vector specifying a subset of the smooth
functions specified in the call to 

prob  For univarite smooths, a scalar between 0 and 1 governing the width of the uncertainty interval. 

facet_args  An optional named list of arguments passed to


alpha, size  For univariate smooths, passed to

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. 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 groupspecific 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) groupspecific terms.
The gamm4
function in the gamm4 package uses
groupspecific terms to implement the departure from linearity in the smooth
terms, but that is not the case for stan_gamm4
where the groupspecific
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), 122. https://www.jstatsoft.org/article/view/v014i14
stanregmethods
and
gamm4
.
The vignette for stan_glmer
, which also discusses
stan_gamm4
. http://mcstan.org/rstanarm/articles/
# from example(gamm4, package = "gamm4"), prefixing gamm4() call with stan_ # \donttest{ dat < mgcv::gamSim(1, n = 400, scale = 2) ## simulate 4 term additive truth#> Gu & Wahba 4 term additive model## 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#> #> SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 9.8e05 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 500 [ 0%] (Warmup) #> Chain 1: Iteration: 50 / 500 [ 10%] (Warmup) #> Chain 1: Iteration: 100 / 500 [ 20%] (Warmup) #> Chain 1: Iteration: 150 / 500 [ 30%] (Warmup) #> Chain 1: Iteration: 200 / 500 [ 40%] (Warmup) #> Chain 1: Iteration: 250 / 500 [ 50%] (Warmup) #> Chain 1: Iteration: 251 / 500 [ 50%] (Sampling) #> Chain 1: Iteration: 300 / 500 [ 60%] (Sampling) #> Chain 1: Iteration: 350 / 500 [ 70%] (Sampling) #> Chain 1: Iteration: 400 / 500 [ 80%] (Sampling) #> Chain 1: Iteration: 450 / 500 [ 90%] (Sampling) #> Chain 1: Iteration: 500 / 500 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 3.96087 seconds (Warmup) #> Chain 1: 1.7028 seconds (Sampling) #> Chain 1: 5.66367 seconds (Total) #> Chain 1:#> Warning: There were 8 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See #> http://mcstan.org/misc/warnings.html#divergenttransitionsafterwarmup#> Warning: Examine the pairs() plot to diagnose sampling problems#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#bulkess#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#tailessprint(br)#> stan_gamm4 #> family: gaussian [identity] #> formula: y ~ s(x0) + x1 + s(x2) #> observations: 400 #>  #> Median MAD_SD #> (Intercept) 4.3 0.2 #> x1 6.4 0.4 #> s(x0).1 0.9 2.2 #> s(x0).2 0.5 1.9 #> s(x0).3 0.5 2.2 #> s(x0).4 0.6 1.8 #> s(x0).5 1.3 1.9 #> s(x0).6 1.9 1.2 #> s(x0).7 0.7 0.8 #> s(x0).8 2.7 1.6 #> s(x0).9 0.1 1.1 #> s(x2).1 42.9 12.0 #> s(x2).2 5.7 8.1 #> s(x2).3 39.5 8.8 #> s(x2).4 40.3 6.1 #> s(x2).5 4.2 4.5 #> s(x2).6 10.4 2.0 #> s(x2).7 8.9 1.6 #> s(x2).8 9.8 5.1 #> s(x2).9 3.5 5.4 #> #> Auxiliary parameter(s): #> Median MAD_SD #> sigma 2.0 0.1 #> #> Smoothing terms: #> Median MAD_SD #> smooth_sd[s(x0)1] 2.3 0.8 #> smooth_sd[s(x0)2] 1.5 1.4 #> smooth_sd[s(x2)1] 18.6 2.9 #> smooth_sd[s(x2)2] 3.6 3.8 #> #> Error terms: #> Groups Name Std.Dev. #> fac (Intercept) 0.45 #> Residual 2.02 #> Num. levels: fac 20 #> #>  #> * For help interpreting the printed output see ?print.stanreg #> * For info on the priors used see ?prior_summary.stanregplot_nonlinear(br)plot_nonlinear(br, smooths = "s(x0)", alpha = 2/3)# }