if (FALSE) {
# ========================================================================
# 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
)
# Inspect the model file with default mvgam priors
stancode(mod_default)
# 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
# ========================================================================
# 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
)
stancode(mod)
# 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
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
)
stancode(mod)
# ========================================================================
# 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
)
stancode(mod)
# ========================================================================
# 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
)
stancode(mod2)
# ========================================================================
# 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
mod2 <- mvgam(
y ~ cov + s(season),
data = simdat$data_train,
trend_model = AR(),
family = poisson(),
priors = brmsprior,
run_model = FALSE
)
stancode(mod2)
# ========================================================================
# 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)
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)
}
Run the code above in your browser using DataLab