Skip to contents

Extract information on default prior distributions for an mvgam model

Source: R/get_mvgam_priors.R
get_mvgam_priors.Rd

This function lists the parameters that can have their prior distributions changed for a given model, as well listing their default distributions

Usage

get_mvgam_priors(
  formula,
  trend_formula,
  factor_formula,
  knots,
  trend_knots,
  trend_model = "None",
  family = poisson(),
  data,
  unit = time,
  species = series,
  use_lv = FALSE,
  n_lv,
  trend_map,
  ...
)

Arguments

formula

A formula object specifying the GAM observation model formula. These are exactly like the formula for a GLM except that smooth terms, s(), te(), ti(), t2(), as well as time-varying dynamic() terms, nonparametric gp() terms and offsets using offset(), can be added to the right hand side to specify that the linear predictor depends on smooth functions of predictors (or linear functionals of these).

In nmix() family models, the formula is used to set up a linear predictor for the detection probability. Details of the formula syntax used by mvgam can be found in mvgam_formulae

trend_formula

An optional formula object specifying the GAM process model formula. If supplied, a linear predictor will be modelled for the latent trends to capture process model evolution separately from the observation model.

Important notes:

  • Should not have a response variable specified on the left-hand side (e.g., ~ season + s(year))

  • Use trend instead of series for effects that vary across time series

  • Only available for RW(), AR() and VAR() trend models

  • In nmix() family models, sets up linear predictor for latent abundance

  • Consider dropping one intercept using - 1 convention to avoid estimation challenges

factor_formula

Can be supplied instead trend_formula to match syntax from jsdgam

knots

An optional list containing user specified knot values for basis construction. For most bases the user simply supplies the knots to be used, which must match up with the k value supplied. Different terms can use different numbers of knots, unless they share a covariate.

trend_knots

As for knots above, this is an optional list of knot values for smooth functions within the trend_formula.

trend_model

character or function specifying the time series dynamics for the latent trend.

Available options:

  • None: No latent trend component (GAM component only, like gam)

  • ZMVN or ZMVN(): Zero-Mean Multivariate Normal (Stan only)

  • 'RW' or RW(): Random Walk

  • 'AR1', 'AR2', 'AR3' or AR(p = 1, 2, 3): Autoregressive models

  • 'CAR1' or CAR(p = 1): Continuous-time AR (Ornstein–Uhlenbeck process)

  • 'VAR1' or VAR(): Vector Autoregressive (Stan only)

  • 'PWlogistic', 'PWlinear' or PW(): Piecewise trends (Stan only)

  • 'GP' or GP(): Gaussian Process with squared exponential kernel (Stan only)

Additional features:

  • Moving average and/or correlated process error terms available for most types (e.g., RW(cor = TRUE) for multivariate Random Walk)

  • Hierarchical correlations possible for structured data

  • See mvgam_trends for details and ZMVN() for examples

family

family specifying the exponential observation family for the series.

Supported families:

  • gaussian(): Real-valued data

  • betar(): Proportional data on (0,1)

  • lognormal(): Non-negative real-valued data

  • student_t(): Real-valued data

  • Gamma(): Non-negative real-valued data

  • bernoulli(): Binary data

  • poisson(): Count data (default)

  • nb(): Overdispersed count data

  • binomial(): Count data with imperfect detection when number of trials is known (use cbind() to bind observations and trials)

  • beta_binomial(): As binomial() but allows for overdispersion

  • nmix(): Count data with imperfect detection when number of trials is unknown (State-Space N-Mixture model with Poisson latent states and Binomial observations)

See mvgam_families for more details.

data

A dataframe or list containing the model response variable and covariates required by the GAM formula and optional trend_formula.

Required columns for most models:

  • series: A factor index of the series IDs (number of levels should equal number of unique series labels)

  • time: numeric or integer index of time points. For most dynamic trend types, time should be measured in discrete, regularly spaced intervals (i.e., c(1, 2, 3, ...)). Irregular spacing is allowed for trend_model = CAR(1), but zero intervals are adjusted to 1e-12 to prevent sampling errors.

Special cases:

  • Models with hierarchical temporal correlation (e.g., AR(gr = region, subgr = species)) should NOT include a series identifier

  • Models without temporal dynamics (trend_model = 'None' or trend_model = ZMVN()) don't require a time variable

unit

The unquoted name of the variable that represents the unit of analysis in data over which latent residuals should be correlated. This variable should be either a numeric or integer variable in the supplied data. Defaults to time to be consistent with other functionalities in mvgam, though note that the data need not be time series in this case. See examples below for further details and explanations

species

The unquoted name of the factor variable that indexes the different response units in data (usually 'species' in a JSDM). Defaults to series to be consistent with other mvgam models

use_lv

logical. If TRUE, use dynamic factors to estimate series' latent trends in a reduced dimension format. Only available for RW(), AR() and GP() trend models. Default is FALSE. See lv_correlations for examples.

n_lv

integer specifying the number of latent dynamic factors to use if use_lv == TRUE. Cannot exceed n_series. Default is min(2, floor(n_series / 2)).

trend_map

Optional data.frame specifying which series should depend on which latent trends. Enables multiple series to depend on the same latent trend process with different observation processes.

Required structure:

  • Column series: Single unique entry for each series (matching factor levels in data)

  • Column trend: Integer values indicating which trend each series depends on

Notes:

  • Sets up latent factor model by enabling use_lv = TRUE

  • Process model intercept is NOT automatically suppressed

  • Not yet supported for continuous time models (CAR())

...

Not currently used

Value

either a data.frame containing the prior definitions (if any suitable priors can be altered by the user) or NULL, indicating that no priors in the model can be modified

Details

Users can supply a model formula, prior to fitting the model, so that default priors can be inspected and altered. To make alterations, change the contents of the prior column and supplying this data.frame to the mvgam or jsdgam functions using the argument priors. If using Stan as the backend, users can also modify the parameter bounds by modifying the new_lowerbound and/or new_upperbound columns. This will be necessary if using restrictive distributions on some parameters, such as a Beta distribution for the trend sd parameters for example (Beta only has support on (0,1)), so the upperbound cannot be above 1. Another option is to make use of the prior modification functions in brms (i.e. prior) to change prior distributions and bounds (just use the name of the parameter that you'd like to change as the class argument; see examples below)

Note

Only the prior, new_lowerbound and/or new_upperbound columns of the output should be altered when defining the user-defined priors for the model. Use only if you are familiar with the underlying probabilistic programming language. There are no sanity checks done to ensure that the code is legal (i.e. to check that lower bounds are smaller than upper bounds, for example)

See also

mvgam, mvgam_formulae, prior

Author

Nicholas J Clark

Examples

# \donttest{
# ========================================================================
# Example 1: Simulate data and inspect default priors
# ========================================================================

dat <- sim_mvgam(trend_rel = 0.5)

# Get a model file that uses default mvgam priors for inspection (not
# always necessary, but this can be useful for testing whether your
# updated priors are written correctly)
mod_default <- mvgam(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2),
  run_model = FALSE
)
#> Your model may benefit from using "noncentred = TRUE"

# Inspect the model file with default mvgam priors
stancode(mod_default)
#> // Stan model code generated by package mvgam
#> functions {
#>   vector rep_each(vector x, int K) {
#>     int N = rows(x);
#>     vector[N * K] y;
#>     int pos = 1;
#>     for (n in 1 : N) {
#>       for (k in 1 : K) {
#>         y[pos] = x[n];
#>         pos += 1;
#>       }
#>     }
#>     return y;
#>   }
#> }
#> data {
#>   int<lower=0> total_obs; // total number of observations
#>   int<lower=0> n; // number of timepoints per series
#>   int<lower=0> n_sp; // number of smoothing parameters
#>   int<lower=0> n_series; // number of series
#>   int<lower=0> num_basis; // total number of basis coefficients
#>   vector[num_basis] zero; // prior locations for basis coefficients
#>   matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#>   array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#>   matrix[8, 8] S1; // mgcv smooth penalty matrix S1
#>   int<lower=0> n_nonmissing; // number of nonmissing observations
#>   array[n_nonmissing] int<lower=0> flat_ys; // flattened nonmissing observations
#>   matrix[n_nonmissing, num_basis] flat_xs; // X values for nonmissing observations
#>   array[n_nonmissing] int<lower=0> obs_ind; // indices of nonmissing observations
#> }
#> parameters {
#>   // raw basis coefficients
#>   vector[num_basis] b_raw;
#>   
#>   // random effect variances
#>   vector<lower=0>[1] sigma_raw;
#>   
#>   // random effect means
#>   vector[1] mu_raw;
#>   
#>   // negative binomial overdispersion
#>   vector<lower=0>[n_series] phi_inv;
#>   
#>   // latent trend AR1 terms
#>   vector<lower=-1, upper=1>[n_series] ar1;
#>   
#>   // latent trend AR2 terms
#>   vector<lower=-1, upper=1>[n_series] ar2;
#>   
#>   // latent trend variance parameters
#>   vector<lower=0>[n_series] sigma;
#>   
#>   // latent trends
#>   matrix[n, n_series] trend;
#>   
#>   // smoothing parameters
#>   vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#>   // basis coefficients
#>   vector[num_basis] b;
#>   b[1 : 8] = b_raw[1 : 8];
#>   b[9 : 11] = mu_raw[1] + b_raw[9 : 11] * sigma_raw[1];
#> }
#> model {
#>   // prior for random effect population variances
#>   sigma_raw ~ inv_gamma(1.418, 0.452);
#>   
#>   // prior for random effect population means
#>   mu_raw ~ std_normal();
#>   
#>   // prior for s(season)...
#>   b_raw[1 : 8] ~ multi_normal_prec(zero[1 : 8], S1[1 : 8, 1 : 8] * lambda[1]);
#>   
#>   // prior (non-centred) for s(series)...
#>   b_raw[9 : 11] ~ std_normal();
#>   
#>   // priors for AR parameters
#>   ar1 ~ std_normal();
#>   ar2 ~ std_normal();
#>   
#>   // priors for smoothing parameters
#>   lambda ~ normal(5, 30);
#>   
#>   // priors for overdispersion parameters
#>   phi_inv ~ student_t(3, 0, 0.1);
#>   
#>   // priors for latent trend variance parameters
#>   sigma ~ inv_gamma(1.418, 0.452);
#>   
#>   // trend estimates
#>   trend[1, 1 : n_series] ~ normal(0, sigma);
#>   trend[2, 1 : n_series] ~ normal(trend[1, 1 : n_series] * ar1, sigma);
#>   for (s in 1 : n_series) {
#>     trend[3 : n, s] ~ normal(ar1[s] * trend[2 : (n - 1), s]
#>                              + ar2[s] * trend[1 : (n - 2), s], sigma[s]);
#>   }
#>   {
#>     // likelihood functions
#>     vector[n_nonmissing] flat_trends;
#>     array[n_nonmissing] real flat_phis;
#>     flat_trends = to_vector(trend)[obs_ind];
#>     flat_phis = to_array_1d(rep_each(phi_inv, n)[obs_ind]);
#>     flat_ys ~ neg_binomial_2(exp(append_col(flat_xs, flat_trends)
#>                                  * append_row(b, 1.0)),
#>                              inv(flat_phis));
#>   }
#> }
#> generated quantities {
#>   vector[total_obs] eta;
#>   matrix[n, n_series] mus;
#>   vector[n_sp] rho;
#>   vector[n_series] tau;
#>   array[n, n_series] int ypred;
#>   matrix[n, n_series] phi_vec;
#>   vector[n_series] phi;
#>   phi = inv(phi_inv);
#>   for (s in 1 : n_series) {
#>     phi_vec[1 : n, s] = rep_vector(phi[s], n);
#>   }
#>   rho = log(lambda);
#>   for (s in 1 : n_series) {
#>     tau[s] = pow(sigma[s], -2.0);
#>   }
#>   
#>   // posterior predictions
#>   eta = X * b;
#>   for (s in 1 : n_series) {
#>     mus[1 : n, s] = eta[ytimes[1 : n, s]] + trend[1 : n, s];
#>     ypred[1 : n, s] = neg_binomial_2_rng(exp(mus[1 : n, s]), phi_vec[1 : n, s]);
#>   }
#> }
#> 
#> 

# Look at which priors can be updated in mvgam
test_priors <- get_mvgam_priors(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2)
)
test_priors
#>                                param_name param_length
#> 1           vector<lower=0>[n_sp] lambda;            2
#> 2                       vector[1] mu_raw;            1
#> 3           vector<lower=0>[1] sigma_raw;            1
#> 4 vector<lower=-1,upper=1>[n_series] ar1;            3
#> 5 vector<lower=-1,upper=1>[n_series] ar2;            3
#> 6        vector<lower=0>[n_series] sigma;            3
#> 7      vector<lower=0>[n_series] phi_inv;            3
#>                    param_info                                prior
#> 1 s(season) smooth parameters              lambda ~ normal(5, 30);
#> 2          s(series) pop mean               mu_raw ~ std_normal();
#> 3            s(series) pop sd sigma_raw ~ inv_gamma(1.418, 0.452);
#> 4       trend AR1 coefficient                  ar1 ~ std_normal();
#> 5       trend AR2 coefficient                  ar2 ~ std_normal();
#> 6                    trend sd     sigma ~ inv_gamma(1.418, 0.452);
#> 7   inverse of NB dispsersion      phi_inv ~ student_t(3, 0, 0.1);
#>                   example_change new_lowerbound new_upperbound
#> 1     lambda ~ exponential(0.6);             NA             NA
#> 2  mu_raw ~ normal(-0.07, 0.61);             NA             NA
#> 3 sigma_raw ~ exponential(0.22);             NA             NA
#> 4     ar1 ~ normal(-0.92, 0.66);             NA             NA
#> 5         ar2 ~ normal(0, 0.23);             NA             NA
#> 6     sigma ~ exponential(0.47);             NA             NA
#> 7  phi_inv ~ normal(-0.24, 0.4);             NA             NA

# ========================================================================
# Example 2: Modify priors manually
# ========================================================================

# Make a few changes; first, change the population mean for the
# series-level random intercepts
test_priors$prior[2] <- "mu_raw ~ normal(0.2, 0.5);"

# Now use stronger regularisation for the series-level AR2 coefficients
test_priors$prior[5] <- "ar2 ~ normal(0, 0.25);"

# Check that the changes are made to the model file without any warnings
# by setting 'run_model = FALSE'
mod <- mvgam(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2),
  priors = test_priors,
  run_model = FALSE
)
#> Your model may benefit from using "noncentred = TRUE"
stancode(mod)
#> // Stan model code generated by package mvgam
#> functions {
#>   vector rep_each(vector x, int K) {
#>     int N = rows(x);
#>     vector[N * K] y;
#>     int pos = 1;
#>     for (n in 1 : N) {
#>       for (k in 1 : K) {
#>         y[pos] = x[n];
#>         pos += 1;
#>       }
#>     }
#>     return y;
#>   }
#> }
#> data {
#>   int<lower=0> total_obs; // total number of observations
#>   int<lower=0> n; // number of timepoints per series
#>   int<lower=0> n_sp; // number of smoothing parameters
#>   int<lower=0> n_series; // number of series
#>   int<lower=0> num_basis; // total number of basis coefficients
#>   vector[num_basis] zero; // prior locations for basis coefficients
#>   matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#>   array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#>   matrix[8, 8] S1; // mgcv smooth penalty matrix S1
#>   int<lower=0> n_nonmissing; // number of nonmissing observations
#>   array[n_nonmissing] int<lower=0> flat_ys; // flattened nonmissing observations
#>   matrix[n_nonmissing, num_basis] flat_xs; // X values for nonmissing observations
#>   array[n_nonmissing] int<lower=0> obs_ind; // indices of nonmissing observations
#> }
#> parameters {
#>   // raw basis coefficients
#>   vector[num_basis] b_raw;
#>   
#>   // random effect variances
#>   vector<lower=0>[1] sigma_raw;
#>   
#>   // random effect means
#>   vector[1] mu_raw;
#>   
#>   // negative binomial overdispersion
#>   vector<lower=0>[n_series] phi_inv;
#>   
#>   // latent trend AR1 terms
#>   vector<lower=-1, upper=1>[n_series] ar1;
#>   
#>   // latent trend AR2 terms
#>   vector<lower=-1, upper=1>[n_series] ar2;
#>   
#>   // latent trend variance parameters
#>   vector<lower=0>[n_series] sigma;
#>   
#>   // latent trends
#>   matrix[n, n_series] trend;
#>   
#>   // smoothing parameters
#>   vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#>   // basis coefficients
#>   vector[num_basis] b;
#>   b[1 : 8] = b_raw[1 : 8];
#>   b[9 : 11] = mu_raw[1] + b_raw[9 : 11] * sigma_raw[1];
#> }
#> model {
#>   // prior for random effect population variances
#>   sigma_raw ~ inv_gamma(1.418, 0.452);
#>   
#>   // prior for random effect population means
#>   mu_raw ~ normal(0.2, 0.5);
#>   
#>   // prior for s(season)...
#>   b_raw[1 : 8] ~ multi_normal_prec(zero[1 : 8], S1[1 : 8, 1 : 8] * lambda[1]);
#>   
#>   // prior (non-centred) for s(series)...
#>   b_raw[9 : 11] ~ std_normal();
#>   
#>   // priors for AR parameters
#>   ar1 ~ std_normal();
#>   ar2 ~ normal(0, 0.25);
#>   
#>   // priors for smoothing parameters
#>   lambda ~ normal(5, 30);
#>   
#>   // priors for overdispersion parameters
#>   phi_inv ~ student_t(3, 0, 0.1);
#>   
#>   // priors for latent trend variance parameters
#>   sigma ~ inv_gamma(1.418, 0.452);
#>   
#>   // trend estimates
#>   trend[1, 1 : n_series] ~ normal(0, sigma);
#>   trend[2, 1 : n_series] ~ normal(trend[1, 1 : n_series] * ar1, sigma);
#>   for (s in 1 : n_series) {
#>     trend[3 : n, s] ~ normal(ar1[s] * trend[2 : (n - 1), s]
#>                              + ar2[s] * trend[1 : (n - 2), s], sigma[s]);
#>   }
#>   {
#>     // likelihood functions
#>     vector[n_nonmissing] flat_trends;
#>     array[n_nonmissing] real flat_phis;
#>     flat_trends = to_vector(trend)[obs_ind];
#>     flat_phis = to_array_1d(rep_each(phi_inv, n)[obs_ind]);
#>     flat_ys ~ neg_binomial_2(exp(append_col(flat_xs, flat_trends)
#>                                  * append_row(b, 1.0)),
#>                              inv(flat_phis));
#>   }
#> }
#> generated quantities {
#>   vector[total_obs] eta;
#>   matrix[n, n_series] mus;
#>   vector[n_sp] rho;
#>   vector[n_series] tau;
#>   array[n, n_series] int ypred;
#>   matrix[n, n_series] phi_vec;
#>   vector[n_series] phi;
#>   phi = inv(phi_inv);
#>   for (s in 1 : n_series) {
#>     phi_vec[1 : n, s] = rep_vector(phi[s], n);
#>   }
#>   rho = log(lambda);
#>   for (s in 1 : n_series) {
#>     tau[s] = pow(sigma[s], -2.0);
#>   }
#>   
#>   // posterior predictions
#>   eta = X * b;
#>   for (s in 1 : n_series) {
#>     mus[1 : n, s] = eta[ytimes[1 : n, s]] + trend[1 : n, s];
#>     ypred[1 : n, s] = neg_binomial_2_rng(exp(mus[1 : n, s]), phi_vec[1 : n, s]);
#>   }
#> }
#> 
#> 

# No warnings, the model is ready for fitting now in the usual way with
# the addition of the 'priors' argument

# ========================================================================
# Example 3: Use brms syntax for prior modification
# ========================================================================

# The same can be done using 'brms' functions; here we will also change
# the ar1 prior and put some bounds on the ar coefficients to enforce
# stationarity; we set the prior using the 'class' argument in all brms
# prior functions
brmsprior <- c(
  prior(normal(0.2, 0.5), class = mu_raw),
  prior(normal(0, 0.25), class = ar1, lb = -1, ub = 1),
  prior(normal(0, 0.25), class = ar2, lb = -1, ub = 1)
)
brmsprior
#>             prior  class coef group resp dpar nlpar   lb   ub tag source
#>  normal(0.2, 0.5) mu_raw                            <NA> <NA>       user
#>   normal(0, 0.25)    ar1                              -1    1       user
#>   normal(0, 0.25)    ar2                              -1    1       user

mod <- mvgam(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2),
  priors = brmsprior,
  run_model = FALSE
)
#> Your model may benefit from using "noncentred = TRUE"
stancode(mod)
#> // Stan model code generated by package mvgam
#> functions {
#>   vector rep_each(vector x, int K) {
#>     int N = rows(x);
#>     vector[N * K] y;
#>     int pos = 1;
#>     for (n in 1 : N) {
#>       for (k in 1 : K) {
#>         y[pos] = x[n];
#>         pos += 1;
#>       }
#>     }
#>     return y;
#>   }
#> }
#> data {
#>   int<lower=0> total_obs; // total number of observations
#>   int<lower=0> n; // number of timepoints per series
#>   int<lower=0> n_sp; // number of smoothing parameters
#>   int<lower=0> n_series; // number of series
#>   int<lower=0> num_basis; // total number of basis coefficients
#>   vector[num_basis] zero; // prior locations for basis coefficients
#>   matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#>   array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#>   matrix[8, 8] S1; // mgcv smooth penalty matrix S1
#>   int<lower=0> n_nonmissing; // number of nonmissing observations
#>   array[n_nonmissing] int<lower=0> flat_ys; // flattened nonmissing observations
#>   matrix[n_nonmissing, num_basis] flat_xs; // X values for nonmissing observations
#>   array[n_nonmissing] int<lower=0> obs_ind; // indices of nonmissing observations
#> }
#> parameters {
#>   // raw basis coefficients
#>   vector[num_basis] b_raw;
#>   
#>   // random effect variances
#>   vector<lower=0>[1] sigma_raw;
#>   
#>   // random effect means
#>   vector[1] mu_raw;
#>   
#>   // negative binomial overdispersion
#>   vector<lower=0>[n_series] phi_inv;
#>   
#>   // latent trend AR1 terms
#>   vector<lower=-1, upper=1>[n_series] ar1;
#>   
#>   // latent trend AR2 terms
#>   vector<lower=-1, upper=1>[n_series] ar2;
#>   
#>   // latent trend variance parameters
#>   vector<lower=0>[n_series] sigma;
#>   
#>   // latent trends
#>   matrix[n, n_series] trend;
#>   
#>   // smoothing parameters
#>   vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#>   // basis coefficients
#>   vector[num_basis] b;
#>   b[1 : 8] = b_raw[1 : 8];
#>   b[9 : 11] = mu_raw[1] + b_raw[9 : 11] * sigma_raw[1];
#> }
#> model {
#>   // prior for random effect population variances
#>   sigma_raw ~ inv_gamma(1.418, 0.452);
#>   
#>   // prior for random effect population means
#>   mu_raw ~ normal(0.2, 0.5);
#>   
#>   // prior for s(season)...
#>   b_raw[1 : 8] ~ multi_normal_prec(zero[1 : 8], S1[1 : 8, 1 : 8] * lambda[1]);
#>   
#>   // prior (non-centred) for s(series)...
#>   b_raw[9 : 11] ~ std_normal();
#>   
#>   // priors for AR parameters
#>   ar1 ~ normal(0, 0.25);
#>   ar2 ~ normal(0, 0.25);
#>   
#>   // priors for smoothing parameters
#>   lambda ~ normal(5, 30);
#>   
#>   // priors for overdispersion parameters
#>   phi_inv ~ student_t(3, 0, 0.1);
#>   
#>   // priors for latent trend variance parameters
#>   sigma ~ inv_gamma(1.418, 0.452);
#>   
#>   // trend estimates
#>   trend[1, 1 : n_series] ~ normal(0, sigma);
#>   trend[2, 1 : n_series] ~ normal(trend[1, 1 : n_series] * ar1, sigma);
#>   for (s in 1 : n_series) {
#>     trend[3 : n, s] ~ normal(ar1[s] * trend[2 : (n - 1), s]
#>                              + ar2[s] * trend[1 : (n - 2), s], sigma[s]);
#>   }
#>   {
#>     // likelihood functions
#>     vector[n_nonmissing] flat_trends;
#>     array[n_nonmissing] real flat_phis;
#>     flat_trends = to_vector(trend)[obs_ind];
#>     flat_phis = to_array_1d(rep_each(phi_inv, n)[obs_ind]);
#>     flat_ys ~ neg_binomial_2(exp(append_col(flat_xs, flat_trends)
#>                                  * append_row(b, 1.0)),
#>                              inv(flat_phis));
#>   }
#> }
#> generated quantities {
#>   vector[total_obs] eta;
#>   matrix[n, n_series] mus;
#>   vector[n_sp] rho;
#>   vector[n_series] tau;
#>   array[n, n_series] int ypred;
#>   matrix[n, n_series] phi_vec;
#>   vector[n_series] phi;
#>   phi = inv(phi_inv);
#>   for (s in 1 : n_series) {
#>     phi_vec[1 : n, s] = rep_vector(phi[s], n);
#>   }
#>   rho = log(lambda);
#>   for (s in 1 : n_series) {
#>     tau[s] = pow(sigma[s], -2.0);
#>   }
#>   
#>   // posterior predictions
#>   eta = X * b;
#>   for (s in 1 : n_series) {
#>     mus[1 : n, s] = eta[ytimes[1 : n, s]] + trend[1 : n, s];
#>     ypred[1 : n, s] = neg_binomial_2_rng(exp(mus[1 : n, s]), phi_vec[1 : n, s]);
#>   }
#> }
#> 
#> 

# ========================================================================
# Example 4: Error handling example
# ========================================================================

# Look at what is returned when an incorrect spelling is used
test_priors$prior[5] <- "ar2_bananas ~ normal(0, 0.25);"
mod <- mvgam(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2),
  priors = test_priors,
  run_model = FALSE
)
#> Warning: no match found in model_file for parameter: ar2_bananas
#> Your model may benefit from using "noncentred = TRUE"
stancode(mod)
#> // Stan model code generated by package mvgam
#> functions {
#>   vector rep_each(vector x, int K) {
#>     int N = rows(x);
#>     vector[N * K] y;
#>     int pos = 1;
#>     for (n in 1 : N) {
#>       for (k in 1 : K) {
#>         y[pos] = x[n];
#>         pos += 1;
#>       }
#>     }
#>     return y;
#>   }
#> }
#> data {
#>   int<lower=0> total_obs; // total number of observations
#>   int<lower=0> n; // number of timepoints per series
#>   int<lower=0> n_sp; // number of smoothing parameters
#>   int<lower=0> n_series; // number of series
#>   int<lower=0> num_basis; // total number of basis coefficients
#>   vector[num_basis] zero; // prior locations for basis coefficients
#>   matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#>   array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#>   matrix[8, 8] S1; // mgcv smooth penalty matrix S1
#>   int<lower=0> n_nonmissing; // number of nonmissing observations
#>   array[n_nonmissing] int<lower=0> flat_ys; // flattened nonmissing observations
#>   matrix[n_nonmissing, num_basis] flat_xs; // X values for nonmissing observations
#>   array[n_nonmissing] int<lower=0> obs_ind; // indices of nonmissing observations
#> }
#> parameters {
#>   // raw basis coefficients
#>   vector[num_basis] b_raw;
#>   
#>   // random effect variances
#>   vector<lower=0>[1] sigma_raw;
#>   
#>   // random effect means
#>   vector[1] mu_raw;
#>   
#>   // negative binomial overdispersion
#>   vector<lower=0>[n_series] phi_inv;
#>   
#>   // latent trend AR1 terms
#>   vector<lower=-1, upper=1>[n_series] ar1;
#>   
#>   // latent trend AR2 terms
#>   vector<lower=-1, upper=1>[n_series] ar2;
#>   
#>   // latent trend variance parameters
#>   vector<lower=0>[n_series] sigma;
#>   
#>   // latent trends
#>   matrix[n, n_series] trend;
#>   
#>   // smoothing parameters
#>   vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#>   // basis coefficients
#>   vector[num_basis] b;
#>   b[1 : 8] = b_raw[1 : 8];
#>   b[9 : 11] = mu_raw[1] + b_raw[9 : 11] * sigma_raw[1];
#> }
#> model {
#>   // prior for random effect population variances
#>   sigma_raw ~ inv_gamma(1.418, 0.452);
#>   
#>   // prior for random effect population means
#>   mu_raw ~ normal(0.2, 0.5);
#>   
#>   // prior for s(season)...
#>   b_raw[1 : 8] ~ multi_normal_prec(zero[1 : 8], S1[1 : 8, 1 : 8] * lambda[1]);
#>   
#>   // prior (non-centred) for s(series)...
#>   b_raw[9 : 11] ~ std_normal();
#>   
#>   // priors for AR parameters
#>   ar1 ~ std_normal();
#>   ar2 ~ std_normal();
#>   
#>   // priors for smoothing parameters
#>   lambda ~ normal(5, 30);
#>   
#>   // priors for overdispersion parameters
#>   phi_inv ~ student_t(3, 0, 0.1);
#>   
#>   // priors for latent trend variance parameters
#>   sigma ~ inv_gamma(1.418, 0.452);
#>   
#>   // trend estimates
#>   trend[1, 1 : n_series] ~ normal(0, sigma);
#>   trend[2, 1 : n_series] ~ normal(trend[1, 1 : n_series] * ar1, sigma);
#>   for (s in 1 : n_series) {
#>     trend[3 : n, s] ~ normal(ar1[s] * trend[2 : (n - 1), s]
#>                              + ar2[s] * trend[1 : (n - 2), s], sigma[s]);
#>   }
#>   {
#>     // likelihood functions
#>     vector[n_nonmissing] flat_trends;
#>     array[n_nonmissing] real flat_phis;
#>     flat_trends = to_vector(trend)[obs_ind];
#>     flat_phis = to_array_1d(rep_each(phi_inv, n)[obs_ind]);
#>     flat_ys ~ neg_binomial_2(exp(append_col(flat_xs, flat_trends)
#>                                  * append_row(b, 1.0)),
#>                              inv(flat_phis));
#>   }
#> }
#> generated quantities {
#>   vector[total_obs] eta;
#>   matrix[n, n_series] mus;
#>   vector[n_sp] rho;
#>   vector[n_series] tau;
#>   array[n, n_series] int ypred;
#>   matrix[n, n_series] phi_vec;
#>   vector[n_series] phi;
#>   phi = inv(phi_inv);
#>   for (s in 1 : n_series) {
#>     phi_vec[1 : n, s] = rep_vector(phi[s], n);
#>   }
#>   rho = log(lambda);
#>   for (s in 1 : n_series) {
#>     tau[s] = pow(sigma[s], -2.0);
#>   }
#>   
#>   // posterior predictions
#>   eta = X * b;
#>   for (s in 1 : n_series) {
#>     mus[1 : n, s] = eta[ytimes[1 : n, s]] + trend[1 : n, s];
#>     ypred[1 : n, s] = neg_binomial_2_rng(exp(mus[1 : n, s]), phi_vec[1 : n, s]);
#>   }
#> }
#> 
#> 

# ========================================================================
# Example 5: Parametric (fixed effect) priors
# ========================================================================

simdat <- sim_mvgam()

# Add a fake covariate
simdat$data_train$cov <- rnorm(NROW(simdat$data_train))

priors <- get_mvgam_priors(
  y ~ cov + s(season),
  data = simdat$data_train,
  family = poisson(),
  trend_model = AR()
)

# Change priors for the intercept and fake covariate effects
priors$prior[1] <- "(Intercept) ~ normal(0, 1);"
priors$prior[2] <- "cov ~ normal(0, 0.1);"

mod2 <- mvgam(
  y ~ cov + s(season),
  data = simdat$data_train,
  trend_model = AR(),
  family = poisson(),
  priors = priors,
  run_model = FALSE
)
#> Your model may benefit from using "noncentred = TRUE"
stancode(mod2)
#> // Stan model code generated by package mvgam
#> data {
#>   int<lower=0> total_obs; // total number of observations
#>   int<lower=0> n; // number of timepoints per series
#>   int<lower=0> n_sp; // number of smoothing parameters
#>   int<lower=0> n_series; // number of series
#>   int<lower=0> num_basis; // total number of basis coefficients
#>   vector[num_basis] zero; // prior locations for basis coefficients
#>   matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#>   array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#>   matrix[9, 18] S1; // mgcv smooth penalty matrix S1
#>   int<lower=0> n_nonmissing; // number of nonmissing observations
#>   array[n_nonmissing] int<lower=0> flat_ys; // flattened nonmissing observations
#>   matrix[n_nonmissing, num_basis] flat_xs; // X values for nonmissing observations
#>   array[n_nonmissing] int<lower=0> obs_ind; // indices of nonmissing observations
#> }
#> parameters {
#>   // raw basis coefficients
#>   vector[num_basis] b_raw;
#>   
#>   // latent trend AR1 terms
#>   vector<lower=-1, upper=1>[n_series] ar1;
#>   
#>   // latent trend variance parameters
#>   vector<lower=0>[n_series] sigma;
#>   
#>   // latent trends
#>   matrix[n, n_series] trend;
#>   
#>   // smoothing parameters
#>   vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#>   // basis coefficients
#>   vector[num_basis] b;
#>   b[1 : num_basis] = b_raw[1 : num_basis];
#> }
#> model {
#>   // prior for (Intercept)...
#>   b_raw[1] ~ normal(0, 1);
#>   
#>   // prior for cov...
#>   b_raw[2] ~ normal(0, 0.1);
#>   
#>   // prior for s(season)...
#>   b_raw[3 : 11] ~ multi_normal_prec(zero[3 : 11],
#>                                     S1[1 : 9, 1 : 9] * lambda[1]
#>                                     + S1[1 : 9, 10 : 18] * lambda[2]);
#>   
#>   // priors for AR parameters
#>   ar1 ~ std_normal();
#>   
#>   // priors for smoothing parameters
#>   lambda ~ normal(5, 30);
#>   
#>   // priors for latent trend variance parameters
#>   sigma ~ inv_gamma(1.418, 0.452);
#>   
#>   // trend estimates
#>   trend[1, 1 : n_series] ~ normal(0, sigma);
#>   for (s in 1 : n_series) {
#>     trend[2 : n, s] ~ normal(ar1[s] * trend[1 : (n - 1), s], sigma[s]);
#>   }
#>   {
#>     // likelihood functions
#>     vector[n_nonmissing] flat_trends;
#>     flat_trends = to_vector(trend)[obs_ind];
#>     flat_ys ~ poisson_log_glm(append_col(flat_xs, flat_trends), 0.0,
#>                               append_row(b, 1.0));
#>   }
#> }
#> generated quantities {
#>   vector[total_obs] eta;
#>   matrix[n, n_series] mus;
#>   vector[n_sp] rho;
#>   vector[n_series] tau;
#>   array[n, n_series] int ypred;
#>   rho = log(lambda);
#>   for (s in 1 : n_series) {
#>     tau[s] = pow(sigma[s], -2.0);
#>   }
#>   
#>   // posterior predictions
#>   eta = X * b;
#>   for (s in 1 : n_series) {
#>     mus[1 : n, s] = eta[ytimes[1 : n, s]] + trend[1 : n, s];
#>     ypred[1 : n, s] = poisson_log_rng(mus[1 : n, s]);
#>   }
#> }
#> 
#> 

# ========================================================================
# Example 6: Alternative brms syntax for fixed effects
# ========================================================================

# Likewise using 'brms' utilities (note that you can use Intercept rather
# than `(Intercept)`) to change priors on the intercept
brmsprior <- c(
  prior(normal(0.2, 0.5), class = cov),
  prior(normal(0, 0.25), class = Intercept)
)
brmsprior
#>             prior     class coef group resp dpar nlpar   lb   ub tag source
#>  normal(0.2, 0.5)       cov                            <NA> <NA>       user
#>   normal(0, 0.25) Intercept                            <NA> <NA>       user

mod2 <- mvgam(
  y ~ cov + s(season),
  data = simdat$data_train,
  trend_model = AR(),
  family = poisson(),
  priors = brmsprior,
  run_model = FALSE
)
#> Your model may benefit from using "noncentred = TRUE"
stancode(mod2)
#> // Stan model code generated by package mvgam
#> data {
#>   int<lower=0> total_obs; // total number of observations
#>   int<lower=0> n; // number of timepoints per series
#>   int<lower=0> n_sp; // number of smoothing parameters
#>   int<lower=0> n_series; // number of series
#>   int<lower=0> num_basis; // total number of basis coefficients
#>   vector[num_basis] zero; // prior locations for basis coefficients
#>   matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#>   array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#>   matrix[9, 18] S1; // mgcv smooth penalty matrix S1
#>   int<lower=0> n_nonmissing; // number of nonmissing observations
#>   array[n_nonmissing] int<lower=0> flat_ys; // flattened nonmissing observations
#>   matrix[n_nonmissing, num_basis] flat_xs; // X values for nonmissing observations
#>   array[n_nonmissing] int<lower=0> obs_ind; // indices of nonmissing observations
#> }
#> parameters {
#>   // raw basis coefficients
#>   vector[num_basis] b_raw;
#>   
#>   // latent trend AR1 terms
#>   vector<lower=-1, upper=1>[n_series] ar1;
#>   
#>   // latent trend variance parameters
#>   vector<lower=0>[n_series] sigma;
#>   
#>   // latent trends
#>   matrix[n, n_series] trend;
#>   
#>   // smoothing parameters
#>   vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#>   // basis coefficients
#>   vector[num_basis] b;
#>   b[1 : num_basis] = b_raw[1 : num_basis];
#> }
#> model {
#>   // prior for (Intercept)...
#>   b_raw[1] ~ normal(0, 0.25);
#>   
#>   // prior for cov...
#>   b_raw[2] ~ normal(0.2, 0.5);
#>   
#>   // prior for s(season)...
#>   b_raw[3 : 11] ~ multi_normal_prec(zero[3 : 11],
#>                                     S1[1 : 9, 1 : 9] * lambda[1]
#>                                     + S1[1 : 9, 10 : 18] * lambda[2]);
#>   
#>   // priors for AR parameters
#>   ar1 ~ std_normal();
#>   
#>   // priors for smoothing parameters
#>   lambda ~ normal(5, 30);
#>   
#>   // priors for latent trend variance parameters
#>   sigma ~ inv_gamma(1.418, 0.452);
#>   
#>   // trend estimates
#>   trend[1, 1 : n_series] ~ normal(0, sigma);
#>   for (s in 1 : n_series) {
#>     trend[2 : n, s] ~ normal(ar1[s] * trend[1 : (n - 1), s], sigma[s]);
#>   }
#>   {
#>     // likelihood functions
#>     vector[n_nonmissing] flat_trends;
#>     flat_trends = to_vector(trend)[obs_ind];
#>     flat_ys ~ poisson_log_glm(append_col(flat_xs, flat_trends), 0.0,
#>                               append_row(b, 1.0));
#>   }
#> }
#> generated quantities {
#>   vector[total_obs] eta;
#>   matrix[n, n_series] mus;
#>   vector[n_sp] rho;
#>   vector[n_series] tau;
#>   array[n, n_series] int ypred;
#>   rho = log(lambda);
#>   for (s in 1 : n_series) {
#>     tau[s] = pow(sigma[s], -2.0);
#>   }
#>   
#>   // posterior predictions
#>   eta = X * b;
#>   for (s in 1 : n_series) {
#>     mus[1 : n, s] = eta[ytimes[1 : n, s]] + trend[1 : n, s];
#>     ypred[1 : n, s] = poisson_log_rng(mus[1 : n, s]);
#>   }
#> }
#> 
#> 

# ========================================================================
# Example 7: Bulk prior assignment
# ========================================================================

# The "class = 'b'" shortcut can be used to put the same prior on all
# 'fixed' effect coefficients (apart from any intercepts)
set.seed(0)
dat <- mgcv::gamSim(1, n = 200, scale = 2)
#> Gu & Wahba 4 term additive model
dat$time <- 1:NROW(dat)
mod <- mvgam(
  y ~ x0 + x1 + s(x2) + s(x3),
  priors = prior(normal(0, 0.75), class = "b"),
  data = dat,
  family = gaussian(),
  run_model = FALSE
)
stancode(mod)
#> // Stan model code generated by package mvgam
#> functions {
#>   vector rep_each(vector x, int K) {
#>     int N = rows(x);
#>     vector[N * K] y;
#>     int pos = 1;
#>     for (n in 1 : N) {
#>       for (k in 1 : K) {
#>         y[pos] = x[n];
#>         pos += 1;
#>       }
#>     }
#>     return y;
#>   }
#> }
#> data {
#>   int<lower=0> total_obs; // total number of observations
#>   int<lower=0> n; // number of timepoints per series
#>   int<lower=0> n_sp; // number of smoothing parameters
#>   int<lower=0> n_series; // number of series
#>   int<lower=0> num_basis; // total number of basis coefficients
#>   vector[num_basis] zero; // prior locations for basis coefficients
#>   matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#>   array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#>   matrix[9, 18] S1; // mgcv smooth penalty matrix S1
#>   matrix[9, 18] S2; // mgcv smooth penalty matrix S2
#>   int<lower=0> n_nonmissing; // number of nonmissing observations
#>   vector[n_nonmissing] flat_ys; // flattened nonmissing observations
#>   matrix[n_nonmissing, num_basis] flat_xs; // X values for nonmissing observations
#>   array[n_nonmissing] int<lower=0> obs_ind; // indices of nonmissing observations
#> }
#> parameters {
#>   // raw basis coefficients
#>   vector[num_basis] b_raw;
#>   
#>   // gaussian observation error
#>   vector<lower=0>[n_series] sigma_obs;
#>   
#>   // smoothing parameters
#>   vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#>   // basis coefficients
#>   vector[num_basis] b;
#>   b[1 : num_basis] = b_raw[1 : num_basis];
#> }
#> model {
#>   // prior for (Intercept)...
#>   b_raw[1] ~ student_t(3, 7.4, 3.7);
#>   
#>   // prior for x0...
#>   b_raw[2] ~ normal(0, 0.75);
#>   
#>   // prior for x1...
#>   b_raw[3] ~ normal(0, 0.75);
#>   
#>   // prior for s(x2)...
#>   b_raw[4 : 12] ~ multi_normal_prec(zero[4 : 12],
#>                                     S1[1 : 9, 1 : 9] * lambda[1]
#>                                     + S1[1 : 9, 10 : 18] * lambda[2]);
#>   
#>   // prior for s(x3)...
#>   b_raw[13 : 21] ~ multi_normal_prec(zero[13 : 21],
#>                                      S2[1 : 9, 1 : 9] * lambda[3]
#>                                      + S2[1 : 9, 10 : 18] * lambda[4]);
#>   
#>   // priors for smoothing parameters
#>   lambda ~ normal(5, 30);
#>   
#>   // priors for observation error parameters
#>   sigma_obs ~ inv_gamma(1.418, 0.452);
#>   {
#>     // likelihood functions
#>     vector[n_nonmissing] flat_sigma_obs;
#>     flat_sigma_obs = rep_each(sigma_obs, n)[obs_ind];
#>     flat_ys ~ normal_id_glm(flat_xs, 0.0, b, flat_sigma_obs);
#>   }
#> }
#> generated quantities {
#>   vector[total_obs] eta;
#>   matrix[n, n_series] sigma_obs_vec;
#>   matrix[n, n_series] mus;
#>   vector[n_sp] rho;
#>   array[n, n_series] real ypred;
#>   rho = log(lambda);
#>   
#>   // posterior predictions
#>   eta = X * b;
#>   for (s in 1 : n_series) {
#>     sigma_obs_vec[1 : n, s] = rep_vector(sigma_obs[s], n);
#>   }
#>   for (s in 1 : n_series) {
#>     mus[1 : n, s] = eta[ytimes[1 : n, s]];
#>     ypred[1 : n, s] = normal_rng(mus[1 : n, s], sigma_obs_vec[1 : n, s]);
#>   }
#> }
#> 
#> 
# }