Fit a Bayesian dynamic GAM to a univariate or multivariate set of time series
Source:R/mvgam.R
mvgam.Rd
This function estimates the posterior distribution for Generalised Additive Models (GAMs) that can include
smooth spline functions, specified in the GAM formula, as well as latent temporal processes,
specified by trend_model
. Further modelling options include State-Space representations to allow covariates
and dynamic processes to occur on the latent 'State' level while also capturing observation-level effects.
Prior specifications are flexible and explicitly encourage users to apply
prior distributions that actually reflect their beliefs. In addition, model fits can easily be assessed and
compared with posterior predictive checks, forecast comparisons and leave-one-out / leave-future-out cross-validation.
Usage
mvgam(
formula,
trend_formula,
knots,
trend_knots,
data,
data_train,
newdata,
data_test,
run_model = TRUE,
prior_simulation = FALSE,
return_model_data = FALSE,
family = "poisson",
share_obs_params = FALSE,
use_lv = FALSE,
n_lv,
trend_map,
trend_model = "None",
drift = FALSE,
chains = 4,
burnin = 500,
samples = 500,
thin = 1,
parallel = TRUE,
threads = 1,
priors,
refit = FALSE,
lfo = FALSE,
residuals = TRUE,
use_stan = TRUE,
backend = getOption("brms.backend", "cmdstanr"),
algorithm = getOption("brms.algorithm", "sampling"),
autoformat = TRUE,
save_all_pars = FALSE,
max_treedepth,
adapt_delta,
jags_path,
...
)
Arguments
- formula
A
character
string 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-varyingdynamic()
terms, 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). Innmix()
family models, theformula
is used to set up a linear predictor for the detection probability. Details of the formula syntax used by mvgam can be found inmvgam_formulae
- trend_formula
An optional
character
string 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. Should not have a response variable specified on the left-hand side of the formula (i.e. a valid option would be~ season + s(year)
). Also note that you should not use the identifierseries
in this formula to specify effects that vary across time series. Instead you should usetrend
. This will ensure that models in which atrend_map
is supplied will still work consistently (i.e. by allowing effects to vary across process models, even when some time series share the same underlying process model). This feature is only currently available forRW()
,AR()
andVAR()
trend models. Innmix()
family models, thetrend_formula
is used to set up a linear predictor for the underlying latent abundance- knots
An optional
list
containing user specified knot values to be used for basis construction. For most bases the user simply supplies the knots to be used, which must match up with the k value supplied (note that the number of knots is not always justk
). Different terms can use different numbers of knots, unless they share a covariate- trend_knots
As for
knots
above, this is an optionallist
of knot values for smooth functions within thetrend_formula
- data
A
dataframe
orlist
containing the model response variable and covariates required by the GAMformula
and optionaltrend_formula
. Should include columns: #'series
(afactor
index of the series IDs; the number of levels should be identical to the number of unique series labels (i.e.n_series = length(levels(data$series))
))time
(numeric
orinteger
index of the time point for each observation). For most dynamic trend types available inmvgam
(see argumenttrend_model
), time should be measured in discrete, regularly spaced intervals (i.e.c(1, 2, 3, ...)
). However you can use irregularly spaced intervals if usingtrend_model = CAR(1)
, though note that any temporal intervals that are exactly0
will be adjusted to a very small number (1e-12
) to prevent sampling errors. See an example ofCAR()
trends inCAR
Should also include any other variables to be included in the linear predictor of
formula
- data_train
Deprecated. Still works in place of
data
but users are recommended to usedata
instead for more seamless integration intoR
workflows- newdata
Optional
dataframe
orlist
of test data containing at leastseries
andtime
in addition to any other variables included in the linear predictor offormula
. If included, the observations in variabley
will be set toNA
when fitting the model so that posterior simulations can be obtained- data_test
Deprecated. Still works in place of
newdata
but users are recommended to usenewdata
instead for more seamless integration intoR
workflows- run_model
logical
. IfFALSE
, the model is not fitted but instead the function will return the model file and the data / initial values that are needed to fit the model outside ofmvgam
- prior_simulation
logical
. IfTRUE
, no observations are fed to the model, and instead simulations from prior distributions are returned- return_model_data
logical
. IfTRUE
, the list of data that is needed to fit the model is returned, along with the initial values for smooth and AR parameters, once the model is fitted. This will be helpful if users wish to modify the model file to add other stochastic elements that are not currently available inmvgam
. Default isFALSE
to reduce the size of the returned object, unlessrun_model == FALSE
- family
family
specifying the exponential observation family for the series. Currently supported families are:gaussian()
for real-valued databetar()
for proportional data on(0,1)
lognormal()
for non-negative real-valued datastudent_t()
for real-valued dataGamma()
for non-negative real-valued databernoulli()
for binary datanb()
for count datapoisson()
for count databinomial()
for count data with imperfect detection when the number of trials is known; note that thecbind()
function must be used to bind the discrete observations and the number of trialsnmix()
for count data with imperfect detection when the number of trials is unknown and should be modeled via a State-Space N-Mixture model. The latent states are Poisson, capturing the 'true' latent abundance, while the observation process is Binomial to account for imperfect detection. Seemvgam_families
for an example of how to use this family
Note that only
nb()
andpoisson()
are available if usingJAGS
as the backend. Default ispoisson()
. Seemvgam_families
for more details- share_obs_params
logical
. IfTRUE
and thefamily
has additional family-specific observation parameters (e.g. variance components instudent_t()
orgaussian()
, or dispersion parameters innb()
orbetar()
), these parameters will be shared across all series. This is handy if you have multiple time series that you believe share some properties, such as being from the same species over different spatial units. Default isFALSE
.- use_lv
logical
. IfTRUE
, use dynamic factors to estimate series' latent trends in a reduced dimension format. Only available forRW()
,AR()
andGP()
trend models. Defaults toFALSE
- n_lv
integer
the number of latent dynamic factors to use ifuse_lv == TRUE
. Cannot be> n_series
. Defaults arbitrarily tomin(2, floor(n_series / 2))
- trend_map
Optional
data.frame
specifying which series should depend on which latent trends. Useful for allowing multiple series to depend on the same latent trend process, but with different observation processes. If supplied, a latent factor model is set up by settinguse_lv = TRUE
and using the mapping to set up the shared trends. Needs to have column namesseries
andtrend
, with integer values in thetrend
column to state which trend each series should depend on. Theseries
column should have a single unique entry for each series in the data (names should perfectly match factor levels of theseries
variable indata
). See examples for details- trend_model
character
orfunction
specifying the time series dynamics for the latent trend. Options are:None
(no latent trend component; i.e. the GAM component is all that contributes to the linear predictor, and the observation process is the only source of error; similarly to what is estimated bygam
)'RW'
orRW()
'AR1'
orAR(p = 1)
'AR2'
orAR(p = 2)
'AR3'
orAR(p = 3)
'CAR1'
orCAR(p = 1)
'VAR1'
orVAR()
(only available inStan
)'PWlogistic
,'PWlinear'
orPW()
(only available inStan
)'GP'
orGP()
(Gaussian Process with squared exponential kernel; only available inStan
)
For all trend types apart from
GP()
,CAR()
andPW()
, moving average and/or correlated process error terms can also be estimated (for example,RW(cor = TRUE)
will set up a multivariate Random Walk ifn_series > 1
). See mvgam_trends for more details- drift
logical
estimate a drift parameter in the latent trend components. Useful if the latent trend is expected to broadly follow a non-zero slope. Only available forRW()
andAR()
trend models. Note that if the latent trend is more or less stationary, the drift parameter can become unidentifiable, especially if an intercept term is included in the GAM linear predictor (which it is by default when callingjagam
). Drift parameters will also likely be unidentifiable if using dynamic factor models. Therefore this defaults toFALSE
- chains
integer
specifying the number of parallel chains for the model. Ignored ifalgorithm %in% c('meanfield', 'fullrank', 'pathfinder', 'laplace')
- burnin
integer
specifying the number of warmup iterations of the Markov chain to run to tune sampling algorithms. Ignored ifalgorithm %in% c('meanfield', 'fullrank', 'pathfinder', 'laplace')
- samples
integer
specifying the number of post-warmup iterations of the Markov chain to run for sampling the posterior distribution- thin
Thinning interval for monitors. Ignored if
algorithm %in% c('meanfield', 'fullrank', 'pathfinder', 'laplace')
- parallel
logical
specifying whether multiple cores should be used for generating MCMC simulations in parallel. IfTRUE
, the number of cores to use will bemin(c(chains, parallel::detectCores() - 1))
- threads
integer
Experimental option to use multithreading for within-chain parallelisation inStan
. We recommend its use only if you are experienced withStan
'sreduce_sum
function and have a slow running model that cannot be sped up by any other means. Only available for some families(poisson()
,nb()
,gaussian()
) and when usingCmdstan
as the backend- priors
An optional
data.frame
with prior definitions (in JAGS or Stan syntax). if using Stan, this can also be an object of classbrmsprior
(see.prior
for details). See get_mvgam_priors and 'Details' for more information on changing default prior distributions- refit
Logical indicating whether this is a refit, called using update.mvgam. Users should leave as
FALSE
- lfo
Logical indicating whether this is part of a call to lfo_cv.mvgam. Returns a lighter version of the model with no residuals and fewer monitored parameters to speed up post-processing. But other downstream functions will not work properly, so users should always leave this set as
FALSE
- residuals
Logical indicating whether to compute series-level randomized quantile residuals and include them as part of the returned object. Defaults to
TRUE
, but you can set toFALSE
to save computational time and reduce the size of the returned object (users can always add residuals to an object of classmvgam
using add_residuals)- use_stan
Logical. If
TRUE
, the model will be compiled and sampled using Hamiltonian Monte Carlo with a call tocmdstan_model
or a call tostan
. Note that there are many more options when usingStan
vsJAGS
- backend
Character string naming the package to use as the backend for fitting the Stan model (if
use_stan = TRUE
). Options are "cmdstanr" (the default) or "rstan". Can be set globally for the current R session via the"brms.backend"
option (seeoptions
). Details on the rstan and cmdstanr packages are available at https://mc-stan.org/rstan/ and https://mc-stan.org/cmdstanr/, respectively- algorithm
Character string naming the estimation approach to use. Options are
"sampling"
for MCMC (the default),"meanfield"
for variational inference with factorized normal distributions,"fullrank"
for variational inference with a multivariate normal distribution,"laplace"
for a Laplace approximation (only available when using cmdstanr as the backend) or"pathfinder"
for the pathfinder algorithm (only currently available when using cmdstanr as the backend). Can be set globally for the current R session via the"brms.algorithm"
option (seeoptions
). Limited testing suggests that"meanfield"
performs best out of the non-MCMC approximations for dynamic GAMs, possibly because of the difficulties estimating covariances among the many spline parameters and latent trend parameters. But rigorous testing has not been carried out- autoformat
Logical
. Use thestanc
parser to automatically format theStan
code and check for deprecations. Defaults toTRUE
- save_all_pars
Logical
flag to indicate if draws from all variables defined in Stan'sparameters
block should be saved (default isFALSE
).- max_treedepth
positive integer placing a cap on the number of simulation steps evaluated during each iteration when
use_stan == TRUE
. Default is12
. Increasing this value can sometimes help with exploration of complex posterior geometries, but it is rarely fruitful to go above amax_treedepth
of14
- adapt_delta
positive numeric between
0
and1
defining the target average proposal acceptance probability during Stan's adaptation period, ifuse_stan == TRUE
. Default is0.8
. In general you should not need to change adapt_delta unless you see a warning message about divergent transitions, in which case you can increase adapt_delta from the default to a value closer to1
(e.g. from0.95
to0.99
, or from0.99
to0.999
, etc). The step size used by the numerical integrator is a function ofadapt_delta
in that increasingadapt_delta
will result in a smaller step size and fewer divergences. Increasingadapt_delta
will typically result in a slower sampler, but it will always lead to a more robust sampler- jags_path
Optional character vector specifying the path to the location of the
JAGS
executable (.exe) to use for modelling ifuse_stan == FALSE
. If missing, the path will be recovered from a call tofindjags
- ...
Further arguments passed to Stan. For
backend = "rstan"
the arguments are passed tosampling
orvb
. Forbackend = "cmdstanr"
the arguments are passed to thecmdstanr::sample
,cmdstanr::variational
,cmdstanr::laplace
orcmdstanr::pathfinder
method
Value
A list
object of class mvgam
containing model output, the text representation of the model file,
the mgcv model output (for easily generating simulations at
unsampled covariate values), Dunn-Smyth residuals for each series and key information needed
for other functions in the package. See mvgam-class
for details.
Use methods(class = "mvgam")
for an overview on available methods.
Details
Dynamic GAMs are useful when we wish to predict future values from time series that show temporal dependence
but we do not want to rely on extrapolating from a smooth term (which can sometimes lead to unpredictable and unrealistic behaviours).
In addition, smooths can often try to wiggle excessively to capture any autocorrelation that is present in a time series,
which exacerbates the problem of forecasting ahead. As GAMs are very naturally viewed through a Bayesian lens, and we often
must model time series that show complex distributional features and missing data, parameters for mvgam
models are estimated
in a Bayesian framework using Markov Chain Monte Carlo by default. A general overview is provided
in the primary vignettes: vignette("mvgam_overview")
and vignette("data_in_mvgam")
.
For a full list of available vignettes see vignette(package = "mvgam")
Formula syntax: Details of the formula syntax used by mvgam can be found in
mvgam_formulae
. Note that it is possible to supply an empty formula where
there are no predictors or intercepts in the observation model (i.e. y ~ 0
or y ~ -1
).
In this case, an intercept-only observation model will be set up but the intercept coefficient
will be fixed at zero. This can be handy if you wish to fit pure State-Space models where
the variation in the dynamic trend controls the average expectation, and/or where intercepts
are non-identifiable (as in piecewise trends, see examples below)
Families and link functions: Details of families supported by mvgam can be found in
mvgam_families
.
Trend models: Details of latent trend dynamic models supported by mvgam can be found in
mvgam_trends
.
Priors: A jagam
model file is generated from formula
and
modified to include any latent
temporal processes. Default priors for intercepts and any scale parameters are generated
using the same practice as brms. Prior distributions for most important model parameters can be
altered by the user to inspect model
sensitivities to given priors (see get_mvgam_priors
for details).
Note that latent trends are estimated on the link scale so choose priors
accordingly. However more control over the model specification can be accomplished by first using mvgam
as a
baseline, then editing the returned model accordingly. The model file can be edited and run outside
of mvgam
by setting run_model = FALSE
and this is encouraged for complex
modelling tasks. Note, no priors are
formally checked to ensure they are in the right syntax for the respective
probabilistic modelling framework, so it is
up to the user to ensure these are correct (i.e. use dnorm
for normal
densities in JAGS
, with the mean and precision
parameterisation; but use normal
for normal densities in Stan
, with the mean
and standard deviation parameterisation)
Random effects: For any smooth terms using the random effect basis (smooth.construct.re.smooth.spec
),
a non-centred parameterisation is automatically employed to avoid degeneracies that are common in hierarchical models.
Note however that centred versions may perform better for series that are particularly informative, so as with any
foray into Bayesian modelling, it is worth building an understanding of the model's assumptions and limitations by following a
principled workflow. Also note that models are parameterised using drop.unused.levels = FALSE
in jagam
to ensure predictions can be made for all levels of the supplied factor variable
Observation level parameters: When more than one series is included in data
and an
observation family that contains more than one parameter is used, additional observation family parameters
(i.e. phi
for nb()
or sigma
for gaussian()
) are
by default estimated independently for each series. But if you wish for the series to share
the same observation parameters, set share_obs_params = TRUE
Factor regularisation: When using a dynamic factor model for the trends with JAGS
factor precisions are given
regularized penalty priors to theoretically allow some factors to be dropped from the model by squeezing increasing
factors' variances to zero. This is done to help protect against selecting too many latent factors than are needed to
capture dependencies in the data, so it can often be advantageous to set n_lv
to a slightly larger number. However
larger numbers of factors do come with additional computational costs so these should be balanced as well. When using
Stan
, all factors are parameterised with fixed variance parameters
Residuals: For each series, randomized quantile (i.e. Dunn-Smyth) residuals are calculated for inspecting model diagnostics
If the fitted model is appropriate then Dunn-Smyth residuals will be standard normal in distribution and no
autocorrelation will be evident. When a particular observation is missing, the residual is calculated by comparing independent
draws from the model's posterior distribution
Using Stan: mvgam
is primarily designed to use Hamiltonian Monte Carlo for parameter estimation
via the software Stan
(using either the cmdstanr
or rstan
interface).
There are great advantages when using Stan
over Gibbs / Metropolis Hastings samplers, which includes the option
to estimate smooth latent trends via Hilbert space approximate Gaussian Processes.
This often makes sense for ecological series, which we expect to change smoothly. In mvgam
, latent squared
exponential GP trends are approximated using by default 20
basis functions, which saves computational
costs compared to fitting full GPs while adequately estimating
GP alpha
and rho
parameters. Because of the many advantages of Stan
over JAGS
,
further development of the package will only be applied to Stan
. This includes the planned addition
of more response distributions, plans to handle zero-inflation, and plans to incorporate a greater
variety of trend models. Users are strongly encouraged to opt for Stan
over JAGS
in any proceeding workflows
References
Nicholas J Clark & Konstans Wells (2020). Dynamic generalised additive models (DGAMs) for forecasting discrete ecological time series. Methods in Ecology and Evolution. 14:3, 771-784.
See also
jagam
, gam
, gam.models
,
Examples
# \donttest{
# Simulate a collection of three time series that have shared seasonal dynamics
dat <- sim_mvgam(T = 80, n_series = 3, prop_missing = 0.1,
prop_trend = 0.6)
# Plot key summary statistics for a single series
plot_mvgam_series(data = dat$data_train, series = 1)
# Plot all series together
plot_mvgam_series(data = dat$data_train, series = 'all')
# Formulate a model using Stan where series share a cyclic smooth for
# seasonality and each series has an independent random walk temporal process;
# Set run_model = FALSE to inspect the returned objects
mod1 <- mvgam(formula = y ~ s(season, bs = 'cc', k = 6),
data = dat$data_train,
trend_model = RW(),
family = poisson(),
use_stan = TRUE,
run_model = FALSE)
# View the model code in Stan language
code(mod1)
#> // 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[4, 4] 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 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] ~ student_t(3, 0, 2.5);
#>
#> // prior for s(season)...
#> b_raw[2 : 5] ~ multi_normal_prec(zero[2 : 5], S1[1 : 4, 1 : 4] * lambda[1]);
#>
#> // priors for smoothing parameters
#> lambda ~ normal(5, 30);
#>
#> // priors for latent trend variance parameters
#> sigma ~ student_t(3, 0, 2.5);
#>
#> // trend estimates
#> trend[1, 1 : n_series] ~ normal(0, sigma);
#> for (s in 1 : n_series) {
#> trend[2 : n, s] ~ normal(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]);
#> }
#> }
#>
#>
# Now fit the model
mod1 <- mvgam(formula = y ~ s(season, bs = 'cc', k = 6),
data = dat$data_train,
trend_model = RW(),
family = poisson(),
burnin = 300,
samples = 300,
chains = 2)
#> Using cmdstanr as the backend
#>
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:359,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:16,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486dd2316.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
#> Init values were only set for a subset of parameters.
#> Missing init values for the following parameters:
#> - chain 1: sigma, trend, lambda
#> - chain 2: sigma, trend, lambda
#> Running MCMC with 2 parallel chains...
#>
#> Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: normal_lpdf: Scale parameter[2] is 0, but must be positive! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486dd2316.stan', line 49, column 2 to column 44)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 1 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 2 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 2 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 1 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 2 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 2 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 2 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 1 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 1 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 1 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 2 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 1 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 2 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 2 finished in 1.5 seconds.
#> Chain 1 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 1 finished in 2.0 seconds.
#>
#> Both chains finished successfully.
#> Mean chain execution time: 1.8 seconds.
#> Total execution time: 2.2 seconds.
#>
#> Warning: 2 of 2 chains had an E-BFMI less than 0.2.
#> See https://mc-stan.org/misc/warnings for details.
# Extract the model summary
summary(mod1)
#> GAM formula:
#> y ~ s(season, bs = "cc", k = 6)
#> <environment: 0x0000020e354f6b38>
#>
#> Family:
#> poisson
#>
#> Link function:
#> log
#>
#> Trend model:
#> RW()
#>
#> N series:
#> 3
#>
#> N timepoints:
#> 60
#>
#> Status:
#> Fitted using Stan
#> 2 chains, each with iter = 600; warmup = 300; thin = 1
#> Total post-warmup draws = 600
#>
#>
#> GAM coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> (Intercept) -0.1400 0.270 0.790 1.02 83
#> s(season).1 -0.6200 -0.320 -0.049 1.00 326
#> s(season).2 -1.1000 -0.670 -0.350 1.00 625
#> s(season).3 -0.2900 -0.014 0.310 1.00 523
#> s(season).4 -0.0039 0.290 0.540 1.01 474
#>
#> Approximate significance of GAM smooths:
#> edf Ref.df Chi.sq p-value
#> s(season) 2.94 4 27.6 0.00062 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Latent trend variance estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> sigma[1] 0.031 0.13 0.38 1.15 15
#> sigma[2] 0.140 0.35 0.69 1.05 30
#> sigma[3] 0.082 0.18 0.41 1.08 36
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
#> Rhats above 1.05 found for 12 parameters
#> *Diagnose further to investigate why the chains have not mixed
#> 0 of 600 iterations ended with a divergence (0%)
#> 0 of 600 iterations saturated the maximum tree depth of 12 (0%)
#> Chain 1: E-FMI = 0.0691
#> Chain 2: E-FMI = 0.1523
#> *E-FMI below 0.2 indicates you may need to reparameterize your model
#>
#> Samples were drawn using NUTS(diag_e) at Wed May 01 1:24:35 PM 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split MCMC chains
#> (at convergence, Rhat = 1)
# Plot the estimated historical trend and forecast for one series
plot(mod1, type = 'trend', series = 1)
plot(mod1, type = 'forecast', series = 1)
# Residual diagnostics
plot(mod1, type = 'residuals', series = 1)
resids <- residuals(mod1)
str(resids)
#> num [1:180, 1:4] -0.539 0.784 0.255 -0.307 0.523 ...
#> - attr(*, "dimnames")=List of 2
#> ..$ : NULL
#> ..$ : chr [1:4] "Estimate" "Est.Error" "Q2.5" "Q97.5"
# Compute the forecast using covariate information in data_test
fc <- forecast(mod1, newdata = dat$data_test)
str(fc)
#> List of 16
#> $ call :Class 'formula' language y ~ s(season, bs = "cc", k = 6)
#> .. ..- attr(*, ".Environment")=<environment: 0x0000020e354f6b38>
#> $ trend_call : NULL
#> $ family : chr "poisson"
#> $ family_pars : NULL
#> $ trend_model :List of 4
#> ..$ trend_model: chr "RW"
#> ..$ ma : logi FALSE
#> ..$ cor : logi FALSE
#> ..$ label : language RW()
#> ..- attr(*, "class")= chr "mvgam_trend"
#> $ drift : logi FALSE
#> $ use_lv : logi FALSE
#> $ fit_engine : chr "stan"
#> $ type : chr "response"
#> $ series_names : Factor w/ 3 levels "series_1","series_2",..: 1 2 3
#> $ train_observations:List of 3
#> ..$ series_1: int [1:60] 1 1 2 1 3 1 2 3 2 2 ...
#> ..$ series_2: int [1:60] 3 2 NA NA 0 0 1 0 0 1 ...
#> ..$ series_3: int [1:60] 2 1 NA 0 1 0 0 1 NA 3 ...
#> $ train_times : int [1:60] 1 2 3 4 5 6 7 8 9 10 ...
#> $ test_observations :List of 3
#> ..$ series_1: int [1:20] 2 1 NA 1 0 1 2 5 4 NA ...
#> ..$ series_2: int [1:20] 0 1 0 1 2 1 2 4 8 3 ...
#> ..$ series_3: int [1:20] 0 0 0 0 0 0 NA 0 NA 0 ...
#> $ test_times : int [1:20] 61 62 63 64 65 66 67 68 69 70 ...
#> $ hindcasts :List of 3
#> ..$ series_1: num [1:600, 1:60] 2 3 6 2 1 0 0 2 2 2 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:60] "ypred[1,1]" "ypred[2,1]" "ypred[3,1]" "ypred[4,1]" ...
#> ..$ series_2: num [1:600, 1:60] 2 4 3 4 2 3 0 1 3 4 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:60] "ypred[1,2]" "ypred[2,2]" "ypred[3,2]" "ypred[4,2]" ...
#> ..$ series_3: num [1:600, 1:60] 3 4 2 1 0 1 0 2 2 0 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:60] "ypred[1,3]" "ypred[2,3]" "ypred[3,3]" "ypred[4,3]" ...
#> $ forecasts :List of 3
#> ..$ series_1: int [1:600, 1:20] 3 1 2 1 0 1 2 3 2 2 ...
#> ..$ series_2: int [1:600, 1:20] 5 6 10 4 8 2 1 3 1 8 ...
#> ..$ series_3: int [1:600, 1:20] 1 0 3 0 1 1 1 1 2 0 ...
#> - attr(*, "class")= chr "mvgam_forecast"
plot(fc)
#> Out of sample DRPS:
#> 35.5156611111111
# Plot the estimated seasonal smooth function
plot(mod1, type = 'smooths')
# Plot estimated first derivatives of the smooth
plot(mod1, type = 'smooths', derivatives = TRUE)
# Plot partial residuals of the smooth
plot(mod1, type = 'smooths', residuals = TRUE)
# Plot posterior realisations for the smooth
plot(mod1, type = 'smooths', realisations = TRUE)
# Plot conditional response predictions using marginaleffects
conditional_effects(mod1)
#> Warning: Removed 16 rows containing missing values or values outside the scale range
#> (`geom_point()`).
plot_predictions(mod1, condition = 'season', points = 0.5)
#> Warning: Removed 16 rows containing missing values or values outside the scale range
#> (`geom_point()`).
# Generate posterior predictive checks using bayesplot
pp_check(mod1)
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
#> Warning: NA responses are not shown in 'pp_check'.
# Extract observation model beta coefficient draws as a data.frame
beta_draws_df <- as.data.frame(mod1, variable = 'betas')
head(beta_draws_df)
#> (Intercept) s(season).1 s(season).2 s(season).3 s(season).4
#> 1 0.508462 -0.571868 -0.682598 -0.00488394 0.319060
#> 2 0.453657 -0.480477 -0.648053 0.06537980 0.275158
#> 3 0.177792 -0.404498 -0.483719 -0.06786070 0.307115
#> 4 0.200146 -0.366601 -1.023480 0.03332720 0.463580
#> 5 0.248963 -0.182699 -0.314912 -0.19939700 0.260020
#> 6 0.123240 -0.368114 -1.002160 0.06694740 0.249951
str(beta_draws_df)
#> 'data.frame': 600 obs. of 5 variables:
#> $ (Intercept): num 0.508 0.454 0.178 0.2 0.249 ...
#> $ s(season).1: num -0.572 -0.48 -0.404 -0.367 -0.183 ...
#> $ s(season).2: num -0.683 -0.648 -0.484 -1.023 -0.315 ...
#> $ s(season).3: num -0.00488 0.06538 -0.06786 0.03333 -0.1994 ...
#> $ s(season).4: num 0.319 0.275 0.307 0.464 0.26 ...
# Investigate model fit
options(mc.cores = 1)
loo(mod1)
#> Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
#>
#> Computed from 600 by 164 log-likelihood matrix
#>
#> Estimate SE
#> elpd_loo -246.8 9.9
#> p_loo 23.9 2.6
#> looic 493.5 19.7
#> ------
#> Monte Carlo SE of elpd_loo is NA.
#>
#> Pareto k diagnostic values:
#> Count Pct. Min. n_eff
#> (-Inf, 0.5] (good) 148 90.2% 18
#> (0.5, 0.7] (ok) 14 8.5% 19
#> (0.7, 1] (bad) 2 1.2% 12
#> (1, Inf) (very bad) 0 0.0% <NA>
#> See help('pareto-k-diagnostic') for details.
# Example of supplying a trend_map so that some series can share
# latent trend processes
sim <- sim_mvgam(n_series = 3)
mod_data <- sim$data_train
# Here, we specify only two latent trends; series 1 and 2 share a trend,
# while series 3 has it's own unique latent trend
trend_map <- data.frame(series = unique(mod_data$series),
trend = c(1, 1, 2))
# Fit the model using AR1 trends
mod <- mvgam(y ~ s(season, bs = 'cc', k = 6),
trend_map = trend_map,
trend_model = AR(),
data = mod_data,
return_model_data = TRUE,
burnin = 300,
samples = 300,
chains = 2)
#> Using cmdstanr as the backend
#>
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:359,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:16,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e48369921f5.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
#> Init values were only set for a subset of parameters.
#> Missing init values for the following parameters:
#> - chain 1: sigma, ar1, LV, lambda
#> - chain 2: sigma, ar1, LV, lambda
#> Running MCMC with 2 parallel chains...
#>
#> Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: normal_lpdf: Scale parameter[1] is 0, but must be positive! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e48369921f5.stan', line 70, column 2 to column 37)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 1 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 1 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 2 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 2 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 1 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 2 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 1 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 2 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 1 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 2 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 1 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 2 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 2 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 2 finished in 1.1 seconds.
#> Chain 1 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 1 finished in 1.2 seconds.
#>
#> Both chains finished successfully.
#> Mean chain execution time: 1.2 seconds.
#> Total execution time: 1.3 seconds.
#>
#> Warning: 2 of 2 chains had an E-BFMI less than 0.2.
#> See https://mc-stan.org/misc/warnings for details.
# The mapping matrix is now supplied as data to the model in the 'Z' element
mod1$model_data$Z
#> NULL
code(mod)
#> // 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_lv; // number of dynamic factors
#> int<lower=0> n_sp; // number of smoothing parameters
#> int<lower=0> n_series; // number of series
#> matrix[n_series, n_lv] Z; // matrix mapping series to latent trends
#> 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[4, 4] 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
#> }
#> transformed data {
#>
#> }
#> parameters {
#> // raw basis coefficients
#> vector[num_basis] b_raw;
#>
#> // latent factor SD terms
#> vector<lower=0>[n_lv] sigma;
#>
#> // latent factor AR1 terms
#> vector<lower=-1.5, upper=1.5>[n_lv] ar1;
#>
#> // dynamic factors
#> matrix[n, n_lv] LV;
#>
#> // smoothing parameters
#> vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#> // trends and dynamic factor loading matrix
#> matrix[n, n_series] trend;
#>
#> // basis coefficients
#> vector[num_basis] b;
#> b[1 : num_basis] = b_raw[1 : num_basis];
#>
#> // derived latent trends
#> for (i in 1 : n) {
#> for (s in 1 : n_series) {
#> trend[i, s] = dot_product(Z[s, : ], LV[i, : ]);
#> }
#> }
#> }
#> model {
#> // prior for (Intercept)...
#> b_raw[1] ~ student_t(3, 0, 2.5);
#>
#> // prior for s(season)...
#> b_raw[2 : 5] ~ multi_normal_prec(zero[2 : 5], S1[1 : 4, 1 : 4] * lambda[1]);
#>
#> // priors for AR parameters
#> ar1 ~ std_normal();
#>
#> // priors for smoothing parameters
#> lambda ~ normal(5, 30);
#>
#> // priors for factor SD parameters
#> sigma ~ student_t(3, 0, 2.5);
#>
#> // dynamic factor estimates
#> LV[1, 1 : n_lv] ~ normal(0, sigma);
#> for (j in 1 : n_lv) {
#> LV[2 : n, j] ~ normal(ar1[j] * LV[1 : (n - 1), j], sigma[j]);
#> }
#> {
#> // 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_lv] penalty;
#> array[n, n_series] int ypred;
#> rho = log(lambda);
#> penalty = 1.0 / (sigma .* sigma);
#>
#> matrix[n_series, n_lv] lv_coefs = Z;
#> // 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]);
#> }
#> }
#>
#>
# The first two series share an identical latent trend; the third is different
plot(mod, type = 'trend', series = 1)
plot(mod, type = 'trend', series = 2)
plot(mod, type = 'trend', series = 3)
# Example of how to use dynamic coefficients
# Simulate a time-varying coefficient for the effect of temperature
set.seed(123)
N <- 200
beta_temp <- vector(length = N)
beta_temp[1] <- 0.4
for(i in 2:N){
beta_temp[i] <- rnorm(1, mean = beta_temp[i - 1] - 0.0025, sd = 0.05)
}
plot(beta_temp)
# Simulate a covariate called 'temp'
temp <- rnorm(N, sd = 1)
# Simulate the Gaussian observation process
out <- rnorm(N, mean = 4 + beta_temp * temp,
sd = 0.5)
# Gather necessary data into a data.frame; split into training / testing
data = data.frame(out, temp, time = seq_along(temp))
data_train <- data[1:180,]
data_test <- data[181:200,]
# Fit the model using the dynamic() formula helper
mod <- mvgam(out ~
dynamic(temp,
scale = FALSE,
k = 40),
family = gaussian(),
data = data_train,
newdata = data_test,
burnin = 300,
samples = 300,
chains = 2)
#> Using cmdstanr as the backend
#>
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:359,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:16,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
#> Init values were only set for a subset of parameters.
#> Missing init values for the following parameters:
#> - chain 1: alpha_gp_time_bytemp, rho_gp_time_bytemp, z_gp_time_bytemp, sigma_obs, lambda
#> - chain 2: alpha_gp_time_bytemp, rho_gp_time_bytemp, z_gp_time_bytemp, sigma_obs, lambda
#> Running MCMC with 2 parallel chains...
#>
#> Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: normal_id_glm_lpdf: Scale vector[1] is inf, but must be positive finite! (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4865131b09.stan', line 99, column 4 to column 61)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 1 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 2 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 2 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 1 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 1 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 1 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 2 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 2 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 1 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 2 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 1 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 2 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 1 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 2 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 1 finished in 0.6 seconds.
#> Chain 2 finished in 0.6 seconds.
#>
#> Both chains finished successfully.
#> Mean chain execution time: 0.6 seconds.
#> Total execution time: 0.8 seconds.
#>
#> Warning: 30 of 600 (5.0%) transitions ended with a divergence.
#> See https://mc-stan.org/misc/warnings for details.
# Inspect the model summary, forecast and time-varying coefficient distribution
summary(mod)
#> GAM formula:
#> out ~ gp(time, by = temp, c = 5/4, k = 40, scale = FALSE)
#> <environment: 0x0000020e354f6b38>
#>
#> Family:
#> gaussian
#>
#> Link function:
#> identity
#>
#> Trend model:
#> None
#>
#> N series:
#> 1
#>
#> N timepoints:
#> 200
#>
#> Status:
#> Fitted using Stan
#> 2 chains, each with iter = 600; warmup = 300; thin = 1
#> Total post-warmup draws = 600
#>
#>
#> Observation error parameter estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> sigma_obs[1] 0.44 0.49 0.54 1 580
#>
#> GAM coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> (Intercept) 4.0000 4.0e+00 4.100 1.00 985
#> gp(time):temp.1 1.0000 3.4e+00 7.500 1.00 508
#> gp(time):temp.2 -2.8000 1.6e+00 5.600 1.01 296
#> gp(time):temp.3 -5.8000 -1.4e+00 5.400 1.00 519
#> gp(time):temp.4 -5.1000 -1.1e+00 2.100 1.01 368
#> gp(time):temp.5 -2.7000 4.6e-01 4.100 1.00 554
#> gp(time):temp.6 -2.4000 1.3e-01 3.200 1.00 576
#> gp(time):temp.7 -3.7000 -3.1e-01 2.000 1.00 589
#> gp(time):temp.8 -1.5000 2.1e-01 2.600 1.00 505
#> gp(time):temp.9 -1.0000 2.7e-01 2.800 1.00 497
#> gp(time):temp.10 -2.4000 -2.4e-01 0.720 1.00 444
#> gp(time):temp.11 -1.7000 -4.4e-03 1.400 1.00 762
#> gp(time):temp.12 -0.7600 6.6e-02 1.800 1.00 424
#> gp(time):temp.13 -1.1000 3.6e-09 1.100 1.00 593
#> gp(time):temp.14 -1.7000 -7.4e-04 0.740 1.00 325
#> gp(time):temp.15 -0.8900 2.4e-12 1.200 1.00 630
#> gp(time):temp.16 -0.9200 -2.3e-17 0.990 1.00 641
#> gp(time):temp.17 -0.7000 1.0e-10 0.780 1.00 616
#> gp(time):temp.18 -0.7500 -2.7e-14 0.550 1.00 780
#> gp(time):temp.19 -0.7900 -1.3e-16 0.330 1.00 565
#> gp(time):temp.20 -0.3900 -1.9e-33 0.530 1.00 529
#> gp(time):temp.21 -0.2100 6.4e-12 0.650 1.00 543
#> gp(time):temp.22 -0.4400 -3.0e-32 0.270 1.01 528
#> gp(time):temp.23 -0.5900 -4.9e-23 0.180 1.00 352
#> gp(time):temp.24 -0.2100 5.6e-24 0.470 1.00 801
#> gp(time):temp.25 -0.3200 -3.5e-43 0.140 1.00 401
#> gp(time):temp.26 -0.4600 -2.6e-37 0.068 1.01 246
#> gp(time):temp.27 -0.0840 5.4e-49 0.250 1.00 339
#> gp(time):temp.28 -0.1200 8.7e-58 0.250 1.01 289
#> gp(time):temp.29 -0.1200 -5.5e-71 0.130 1.00 969
#> gp(time):temp.30 -0.1400 3.0e-52 0.110 1.00 387
#> gp(time):temp.31 -0.0690 1.3e-69 0.082 1.00 456
#> gp(time):temp.32 -0.1200 1.4e-78 0.051 1.00 670
#> gp(time):temp.33 -0.0940 6.8e-73 0.072 1.00 570
#> gp(time):temp.34 -0.0650 -1.5e-75 0.062 1.00 678
#> gp(time):temp.35 -0.0320 -8.6e-54 0.050 1.00 714
#> gp(time):temp.36 -0.0360 1.2e-71 0.045 1.00 369
#> gp(time):temp.37 -0.0930 -2.3e-48 0.012 1.01 267
#> gp(time):temp.38 -0.0130 -1.1e-104 0.035 1.00 664
#> gp(time):temp.39 -0.0082 7.5e-116 0.063 1.01 198
#> gp(time):temp.40 -0.0120 6.6e-86 0.014 1.01 255
#> gp(time):temp.41 -0.0100 -1.8e-92 0.031 1.00 723
#>
#> GAM gp term marginal deviation (alpha) and length scale (rho) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> alpha_gp(time):temp 0.18 0.34 0.78 1.00 302
#> rho_gp(time):temp 13.00 33.00 100.00 1.01 106
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
#> Rhat looks reasonable for all parameters
#> 30 of 600 iterations ended with a divergence (5%)
#> *Try running with larger adapt_delta to remove the divergences
#> 0 of 600 iterations saturated the maximum tree depth of 12 (0%)
#> E-FMI indicated no pathological behavior
#>
#> Samples were drawn using NUTS(diag_e) at Wed May 01 1:26:11 PM 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split MCMC chains
#> (at convergence, Rhat = 1)
plot(mod, type = 'smooths')
fc <- forecast(mod, newdata = data_test)
plot(fc)
#> Out of sample CRPS:
#> 6.33971271101183
# Propagating the smooth term shows how the coefficient is expected to evolve
plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 180, lty = 'dashed', lwd = 2)
points(beta_temp, pch = 16)
# Example showing how to incorporate an offset; simulate some count data
# with different means per series
set.seed(100)
dat <- sim_mvgam(prop_trend = 0, mu = c(0, 2, 2), seasonality = 'hierarchical')
# Add offset terms to the training and testing data
dat$data_train$offset <- 0.5 * as.numeric(dat$data_train$series)
dat$data_test$offset <- 0.5 * as.numeric(dat$data_test$series)
# Fit a model that includes the offset in the linear predictor as well as
# hierarchical seasonal smooths
mod <- mvgam(formula = y ~ offset(offset) +
s(series, bs = 're') +
s(season, bs = 'cc') +
s(season, by = series, m = 1, k = 5),
data = dat$data_train,
trend_model = 'None',
burnin = 300,
samples = 300,
chains = 2)
#> Using cmdstanr as the backend
#>
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:359,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:16,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e482d266aa8.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
#> Init values were only set for a subset of parameters.
#> Missing init values for the following parameters:
#> - chain 1: sigma_raw, mu_raw, lambda
#> - chain 2: sigma_raw, mu_raw, lambda
#> Running MCMC with 2 parallel chains...
#>
#> Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: multi_normal_prec_lpdf: Precision matrix is not symmetric. Precision matrix[3,6] = 1.48544e+109, but Precision matrix[6,3] = 1.48544e+109 (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e482d266aa8.stan', line 51, column 2 to column 78)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 1 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 2 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 2 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 1 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 2 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 2 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 1 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 1 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 2 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 2 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 1 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 2 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 2 finished in 3.1 seconds.
#> Chain 1 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 1 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 1 finished in 4.1 seconds.
#>
#> Both chains finished successfully.
#> Mean chain execution time: 3.6 seconds.
#> Total execution time: 4.2 seconds.
#>
#> Warning: 2 of 600 (0.0%) transitions ended with a divergence.
#> See https://mc-stan.org/misc/warnings for details.
# Inspect the model file to see the modification to the linear predictor
# (eta)
code(mod)
#> // 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
#> vector[total_obs] off_set; // offset vector
#> 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
#> matrix[4, 4] S2; // mgcv smooth penalty matrix S2
#> matrix[4, 4] S3; // mgcv smooth penalty matrix S3
#> matrix[4, 4] S4; // mgcv smooth penalty matrix S4
#> 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;
#>
#> // smoothing parameters
#> vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#> // basis coefficients
#> vector[num_basis] b;
#> b[1 : 21] = b_raw[1 : 21];
#> b[22 : 24] = mu_raw[1] + b_raw[22 : 24] * sigma_raw[1];
#> }
#> model {
#> // prior for random effect population variances
#> sigma_raw ~ student_t(3, 0, 2.5);
#>
#> // prior for random effect population means
#> mu_raw ~ std_normal();
#>
#> // prior for (Intercept)...
#> b_raw[1] ~ student_t(3, 1.4, 2.5);
#>
#> // prior for s(season)...
#> b_raw[2 : 9] ~ multi_normal_prec(zero[2 : 9], S1[1 : 8, 1 : 8] * lambda[1]);
#>
#> // prior for s(season):seriesseries_1...
#> b_raw[10 : 13] ~ multi_normal_prec(zero[10 : 13],
#> S2[1 : 4, 1 : 4] * lambda[2]);
#>
#> // prior for s(season):seriesseries_2...
#> b_raw[14 : 17] ~ multi_normal_prec(zero[14 : 17],
#> S3[1 : 4, 1 : 4] * lambda[3]);
#>
#> // prior for s(season):seriesseries_3...
#> b_raw[18 : 21] ~ multi_normal_prec(zero[18 : 21],
#> S4[1 : 4, 1 : 4] * lambda[4]);
#>
#> // prior (non-centred) for s(series)...
#> b_raw[22 : 24] ~ std_normal();
#>
#> // priors for smoothing parameters
#> lambda ~ normal(5, 30);
#> {
#> // likelihood functions
#> flat_ys ~ poisson_log_glm(flat_xs, off_set[obs_ind], b);
#> }
#> }
#> generated quantities {
#> vector[total_obs] eta;
#> matrix[n, n_series] mus;
#> vector[n_sp] rho;
#> array[n, n_series] int ypred;
#> rho = log(lambda);
#>
#> // posterior predictions
#> eta = X * b + off_set;
#> for (s in 1 : n_series) {
#> mus[1 : n, s] = eta[ytimes[1 : n, s]];
#> ypred[1 : n, s] = poisson_log_rng(mus[1 : n, s]);
#> }
#> }
#>
#>
# Forecasts for the first two series will differ in magnitude
fc <- forecast(mod, newdata = dat$data_test)
layout(matrix(1:2, ncol = 2))
plot(fc, series = 1, ylim = c(0, 75))
#> Out of sample DRPS:
#> 22.7671555555556
plot(fc, series = 2, ylim = c(0, 75))
#> Out of sample DRPS:
#> 78.9129777777778
layout(1)
# Changing the offset for the testing data should lead to changes in
# the forecast
dat$data_test$offset <- dat$data_test$offset - 2
fc <- forecast(mod, newdata = dat$data_test)
plot(fc)
#> Out of sample DRPS:
#> 38.7073305555556
# Relative Risks can be computed by fixing the offset to the same value
# for each series
dat$data_test$offset <- rep(1, NROW(dat$data_test))
preds_rr <- predict(mod, type = 'link', newdata = dat$data_test,
summary = FALSE)
series1_inds <- which(dat$data_test$series == 'series_1')
series2_inds <- which(dat$data_test$series == 'series_2')
# Relative Risks are now more comparable among series
layout(matrix(1:2, ncol = 2))
plot(preds_rr[1, series1_inds], type = 'l', col = 'grey75',
ylim = range(preds_rr),
ylab = 'Series1 Relative Risk', xlab = 'Time')
for(i in 2:50){
lines(preds_rr[i, series1_inds], col = 'grey75')
}
plot(preds_rr[1, series2_inds], type = 'l', col = 'darkred',
ylim = range(preds_rr),
ylab = 'Series2 Relative Risk', xlab = 'Time')
for(i in 2:50){
lines(preds_rr[i, series2_inds], col = 'darkred')
}
layout(1)
# Example of logistic growth with possible changepoints
# Simple logistic growth model
dNt = function(r, N, k){
r * N * (k - N)
}
# Iterate growth through time
Nt = function(r, N, t, k) {
for (i in 1:(t - 1)) {
# population at next time step is current population + growth,
# but we introduce several 'shocks' as changepoints
if(i %in% c(5, 15, 25, 41, 45, 60, 80)){
N[i + 1] <- max(1, N[i] + dNt(r + runif(1, -0.1, 0.1),
N[i], k))
} else {
N[i + 1] <- max(1, N[i] + dNt(r, N[i], k))
}
}
N
}
# Simulate expected values
set.seed(11)
expected <- Nt(0.004, 2, 100, 30)
plot(expected, xlab = 'Time')
# Take Poisson draws
y <- rpois(100, expected)
plot(y, xlab = 'Time')
# Assemble data into dataframe and model. We set a
# fixed carrying capacity of 35 for this example, but note that
# this value is not required to be fixed at each timepoint
mod_data <- data.frame(y = y,
time = 1:100,
cap = 35,
series = as.factor('series_1'))
plot_mvgam_series(data = mod_data)
# The intercept is nonidentifiable when using piecewise
# trends because the trend functions have their own offset
# parameters 'm'; it is recommended to always drop intercepts
# when using these trend models
mod <- mvgam(y ~ 0,
trend_model = PW(growth = 'logistic'),
family = poisson(),
data = mod_data,
burnin = 300,
samples = 300,
chains = 2)
#> Using cmdstanr as the backend
#>
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:359,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:16,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e4836da3d16.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
#> Init values were only set for a subset of parameters.
#> Missing init values for the following parameters:
#> - chain 1: k_trend, m_trend, delta_trend
#> - chain 2: k_trend, m_trend, delta_trend
#> Running MCMC with 2 parallel chains...
#>
#> Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 1 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 2 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 1 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 2 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 1 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 1 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 1 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 2 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 2 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 1 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 2 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 1 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 1 finished in 6.0 seconds.
#> Chain 2 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 2 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 2 finished in 7.3 seconds.
#>
#> Both chains finished successfully.
#> Mean chain execution time: 6.7 seconds.
#> Total execution time: 7.4 seconds.
#>
summary(mod)
#> GAM formula:
#> y ~ 1
#> <environment: 0x0000020e354f6b38>
#>
#> Family:
#> poisson
#>
#> Link function:
#> log
#>
#> Trend model:
#> PW(growth = "logistic")
#>
#> N series:
#> 1
#>
#> N timepoints:
#> 100
#>
#> Status:
#> Fitted using Stan
#> 2 chains, each with iter = 600; warmup = 300; thin = 1
#> Total post-warmup draws = 600
#>
#>
#> GAM coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> (Intercept) 0 0 0 NaN NaN
#>
#> Latent trend growth rate estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> k_trend[1] -0.25 -0.14 -0.07 1.01 220
#>
#> Latent trend offset estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> m_trend[1] -16 -4.4 -0.21 1.01 235
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
#> Rhat looks reasonable for all parameters
#> 0 of 600 iterations ended with a divergence (0%)
#> 0 of 600 iterations saturated the maximum tree depth of 12 (0%)
#> E-FMI indicated no pathological behavior
#>
#> Samples were drawn using NUTS(diag_e) at Wed May 01 1:28:32 PM 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split MCMC chains
#> (at convergence, Rhat = 1)
# Plot the posterior hindcast
plot(mod, type = 'forecast')
# View the changepoints with ggplot2 utilities
library(ggplot2)
#> Warning: package ‘ggplot2’ was built under R version 4.3.3
mcmc_plot(mod, variable = 'delta_trend',
regex = TRUE) +
scale_y_discrete(labels = mod$trend_model$changepoints) +
labs(y = 'Potential changepoint',
x = 'Rate change')
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.
# Example showcasing how cbind() is needed for Binomial observations
# Simulate two time series of Binomial trials
trials <- sample(c(20:25), 50, replace = TRUE)
x <- rnorm(50)
detprob1 <- plogis(-0.5 + 0.9*x)
detprob2 <- plogis(-0.1 -0.7*x)
dat <- rbind(data.frame(y = rbinom(n = 50, size = trials, prob = detprob1),
time = 1:50,
series = 'series1',
x = x,
ntrials = trials),
data.frame(y = rbinom(n = 50, size = trials, prob = detprob2),
time = 1:50,
series = 'series2',
x = x,
ntrials = trials))
dat <- dplyr::mutate(dat, series = as.factor(series))
dat <- dplyr::arrange(dat, time, series)
plot_mvgam_series(data = dat, series = 'all')
# Fit a model using the binomial() family; must specify observations
# and number of trials in the cbind() wrapper
mod <- mvgam(cbind(y, ntrials) ~ series + s(x, by = series),
family = binomial(),
data = dat,
burnin = 300,
samples = 300,
chains = 2)
#> Warning: Binomial and Beta-binomial families require cbind(n_successes, n_trials)
#> in the formula left-hand side. Do not use cbind(n_successes, n_failures)!
#> This warning is displayed once per session.
#> Using cmdstanr as the backend
#>
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:359,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:16,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486b806119.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
#> Init values were only set for a subset of parameters.
#> Missing init values for the following parameters:
#> - chain 1: lambda
#> - chain 2: lambda
#> Running MCMC with 2 parallel chains...
#>
#> Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: multi_normal_prec_lpdf: LDLT_Factor of precision parameter is not positive definite. last conditional variance is 2.00674e-247. (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486b806119.stan', line 40, column 2 to line 42, column 70)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: multi_normal_prec_lpdf: LDLT_Factor of precision parameter is not positive definite. last conditional variance is 7.60059e-247. (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486b806119.stan', line 40, column 2 to line 42, column 70)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: multi_normal_prec_lpdf: LDLT_Factor of precision parameter is not positive definite. last conditional variance is 4.55016e-195. (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486b806119.stan', line 45, column 2 to line 47, column 71)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: multi_normal_prec_lpdf: LDLT_Factor of precision parameter is not positive definite. last conditional variance is 8.5849e-308. (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486b806119.stan', line 40, column 2 to line 42, column 70)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: multi_normal_prec_lpdf: LDLT_Factor of precision parameter is not positive definite. last conditional variance is 8.5849e-308. (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486b806119.stan', line 40, column 2 to line 42, column 70)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: multi_normal_prec_lpdf: LDLT_Factor of precision parameter is not positive definite. last conditional variance is 8.5849e-308. (in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpO0S52W/model-2e486b806119.stan', line 40, column 2 to line 42, column 70)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 1 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 1 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 1 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 1 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 2 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 1 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 2 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 2 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 2 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 1 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 1 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 2 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 2 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 1 finished in 0.7 seconds.
#> Chain 2 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 2 finished in 0.8 seconds.
#>
#> Both chains finished successfully.
#> Mean chain execution time: 0.8 seconds.
#> Total execution time: 0.9 seconds.
#>
summary(mod)
#> GAM formula:
#> cbind(y, ntrials) ~ series + s(x, by = series)
#> <environment: 0x0000020e354f6b38>
#>
#> Family:
#> binomial
#>
#> Link function:
#> logit
#>
#> Trend model:
#> None
#>
#> N series:
#> 2
#>
#> N timepoints:
#> 50
#>
#> Status:
#> Fitted using Stan
#> 2 chains, each with iter = 600; warmup = 300; thin = 1
#> Total post-warmup draws = 600
#>
#>
#> GAM coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> (Intercept) -0.590 -0.46000 -0.330 1.01 417
#> seriesseries2 0.130 0.33000 0.500 1.00 387
#> s(x):seriesseries1.1 -0.130 0.00031 0.140 1.00 329
#> s(x):seriesseries1.2 -0.200 -0.00210 0.140 1.02 249
#> s(x):seriesseries1.3 -0.055 -0.00180 0.055 1.01 273
#> s(x):seriesseries1.4 -0.100 0.00130 0.082 1.02 228
#> s(x):seriesseries1.5 -0.050 0.00070 0.047 1.02 245
#> s(x):seriesseries1.6 -0.084 0.00007 0.093 1.02 230
#> s(x):seriesseries1.7 -0.055 -0.00086 0.067 1.03 223
#> s(x):seriesseries1.8 -0.400 -0.00310 0.290 1.02 217
#> s(x):seriesseries1.9 0.700 0.92000 1.100 1.00 339
#> s(x):seriesseries2.1 -0.230 -0.02300 0.098 1.01 144
#> s(x):seriesseries2.2 -0.180 0.00630 0.390 1.03 59
#> s(x):seriesseries2.3 -0.053 0.00450 0.096 1.00 206
#> s(x):seriesseries2.4 -0.130 -0.00320 0.130 1.02 81
#> s(x):seriesseries2.5 -0.059 -0.00310 0.064 1.02 75
#> s(x):seriesseries2.6 -0.140 0.00530 0.110 1.02 63
#> s(x):seriesseries2.7 -0.076 0.00390 0.079 1.02 69
#> s(x):seriesseries2.8 -0.420 -0.01900 0.470 1.02 65
#> s(x):seriesseries2.9 -0.710 -0.52000 -0.240 1.02 85
#>
#> Approximate significance of GAM smooths:
#> edf Ref.df Chi.sq p-value
#> s(x):seriesseries1 0.995 9 236 <2e-16 ***
#> s(x):seriesseries2 0.974 9 107 0.012 *
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
#> Rhats above 1.05 found for 2 parameters
#> *Diagnose further to investigate why the chains have not mixed
#> 0 of 600 iterations ended with a divergence (0%)
#> 0 of 600 iterations saturated the maximum tree depth of 12 (0%)
#> E-FMI indicated no pathological behavior
#>
#> Samples were drawn using NUTS(diag_e) at Wed May 01 1:29:08 PM 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split MCMC chains
#> (at convergence, Rhat = 1)
pp_check(mod, type = "bars_grouped",
group = "series", ndraws = 50)
pp_check(mod, type = "ecdf_overlay_grouped",
group = "series", ndraws = 50)
conditional_effects(mod, type = 'link')
# }