From 6c932eff1951e7f45bcb9e8a427cbd842a655acd Mon Sep 17 00:00:00 2001 From: Nicholas Clark Date: Wed, 8 Nov 2023 17:41:16 +1000 Subject: [PATCH] better documentation and mcmc_plot options --- R/conditional_effects.R | 1 + R/forecast.mvgam.R | 10 ++---- R/formula.mvgam.R | 6 ++-- R/get_mvgam_priors.R | 9 +---- R/hindcast.mvgam.R | 10 ++---- R/index-mvgam.R | 2 +- R/mcmc_plot.mvgam.R | 21 ++++++++--- R/mvgam.R | 1 + R/plot.mvgam.R | 3 +- R/plot_mvgam_factors.R | 2 +- R/plot_mvgam_fc.R | 2 +- R/plot_mvgam_pterms.R | 4 +-- R/plot_mvgam_randomeffects.R | 4 +-- R/plot_mvgam_resids.R | 2 +- R/plot_mvgam_smooth.R | 4 +-- R/plot_mvgam_trend.R | 2 +- R/plot_mvgam_uncertainty.R | 2 +- R/ppc.mvgam.R | 14 ++++++++ R/predict.mvgam.R | 2 +- R/score.mvgam_forecast.R | 12 ++++--- R/update.mvgam.R | 49 ++------------------------ man/conditional_effects.mvgam.Rd | 3 ++ man/forecast.mvgam.Rd | 11 +++--- man/formula.mvgam.Rd | 8 ++--- man/get_mvgam_priors.Rd | 15 ++++---- man/hindcast.mvgam.Rd | 9 ++--- man/index-mvgam.Rd | 2 +- man/mcmc_plot.mvgam.Rd | 12 +++++-- man/mvgam.Rd | 1 + man/plot.mvgam.Rd | 2 +- man/plot_mvgam_factors.Rd | 2 -- man/plot_mvgam_forecasts.Rd | 2 -- man/plot_mvgam_pterms.Rd | 2 -- man/plot_mvgam_randomeffects.Rd | 2 -- man/plot_mvgam_resids.Rd | 2 -- man/plot_mvgam_smooth.Rd | 2 -- man/plot_mvgam_trend.Rd | 2 -- man/plot_mvgam_uncertainty.Rd | 2 -- man/posterior_epred.mvgam.Rd | 2 +- man/posterior_linpred.mvgam.Rd | 2 +- man/posterior_predict.mvgam.Rd | 2 +- man/ppc.mvgam.Rd | 16 +++++++++ man/predict.mvgam.Rd | 2 +- man/score.mvgam_forecast.Rd | 14 +++++--- man/update.mvgam.Rd | 58 +++++++++++++++++++++---------- src/mvgam.dll | Bin 1064960 -> 1064960 bytes tests/testthat/Rplots.pdf | Bin 22882 -> 22851 bytes 47 files changed, 168 insertions(+), 169 deletions(-) diff --git a/R/conditional_effects.R b/R/conditional_effects.R index b0b759ec..116dd70b 100644 --- a/R/conditional_effects.R +++ b/R/conditional_effects.R @@ -35,6 +35,7 @@ #' and create more bespoke conditional effects plots. #' @name conditional_effects.mvgam #' @author Nicholas J Clark +#' @seealso \code{\link[marginaleffects]{plot_predictions}}, \code{\link[marginaleffects]{plot_slopes}} #' @examples #' \dontrun{ #' # Simulate some data diff --git a/R/forecast.mvgam.R b/R/forecast.mvgam.R index 59dfdc39..90b9047c 100644 --- a/R/forecast.mvgam.R +++ b/R/forecast.mvgam.R @@ -3,7 +3,7 @@ #'@importFrom parallel clusterExport stopCluster setDefaultCluster #'@importFrom stats predict #'@importFrom rlang missing_arg -#'@param object \code{list} object returned from \code{mvgam}. See [mvgam()] +#'@inheritParams predict.mvgam #'@param newdata Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' #'in addition to any other variables included in the linear predictor of the original \code{formula}. If included, the #'covariate information in \code{newdata} will be used to generate forecasts from the fitted model equations. If @@ -19,17 +19,11 @@ #'if the fitted model contained multivariate trends (either as a dynamic factor or \code{VAR} process), #'as it saves recomputing the full set of trends for each series individually #'@param n_cores \code{integer} specifying number of cores for generating forecasts in parallel -#'@param type When this has the value \code{link} the linear predictor is calculated on the link scale. -#'If \code{expected} is used, predictions reflect the expectation of the response (the mean) -#'but ignore uncertainty in the observation process. When \code{response} (default) is used, -#'the predictions take uncertainty in the observation process into account to return -#'predictions on the outcome scale. When \code{trend} is used, only the forecast distribution for the -#'latent trend is returned #'@param ... Ignored #'@details Posterior predictions are drawn from the fitted \code{mvgam} and used to simulate a forecast distribution #'@return An object of class \code{mvgam_forecast} containing hindcast and forecast distributions. #'See \code{\link{mvgam_forecast-class}} for details. -#' +#'@seealso \code{\link{hindcast}}, \code{\link{score}} #'@export forecast <- function(object, ...){ UseMethod("forecast", object) diff --git a/R/formula.mvgam.R b/R/formula.mvgam.R index 98814486..f99e7d85 100644 --- a/R/formula.mvgam.R +++ b/R/formula.mvgam.R @@ -1,13 +1,13 @@ -#'Extract model.frame from a fitted mvgam object +#'Extract formulae from mvgam objects #' #'@rdname formula.mvgam #'@param x `mvgam` or `mvgam_prefit` object -#'@param trend_effects \code{logical}, return the model.frame from the +#'@param trend_effects \code{logical}, return the formula from the #'observation model (if \code{FALSE}) or from the underlying process #'model (if\code{TRUE}) #'@param ... Ignored #'@author Nicholas J Clark -#'@return A \code{matrix} containing the fitted model frame +#'@return A \code{formula} object #'@export formula.mvgam = function(x, trend_effects = FALSE, ...){ # Check trend_effects diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R index d6d2c8a5..c058e90d 100644 --- a/R/get_mvgam_priors.R +++ b/R/get_mvgam_priors.R @@ -4,13 +4,6 @@ #'changed for a given `mvgam` model, as well listing their default distributions #' #' @inheritParams mvgam -#'@param use_stan Logical. If \code{TRUE} and if \code{rstan} is installed, the model will be compiled and sampled using -#'the Hamiltonian Monte Carlo with a call to \code{\link[cmdstanr]{cmdstan_model}} or, if `cmdstanr` is not available, -#'a call to \code{\link[rstan]{stan}}. Note that this functionality is still in development and -#'not all options that are available in \code{JAGS} can be used, including: no option for a Tweedie family and no option for -#'dynamic factor trends. However, as \code{Stan} can estimate Hilbert base approximate Gaussian Processes, which -#'are much more computationally tractable than full GPs for time series with `>100` observations, estimation -#'in \code{Stan} can support latent GP trends while estimation in \code{JAGS} cannot #'@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 #'\code{data.frame} to the `mvgam` function using the argument `priors`. If using `Stan` as the backend, @@ -26,7 +19,7 @@ #' ensure that the code is legal (i.e. to check that lower bounds are smaller than upper bounds, for #' example) #'@author Nicholas J Clark -#'@seealso \code{\link{mvgam}} +#'@seealso \code{\link{mvgam}} \code{\link[brms]{prior}} #'@return either a \code{data.frame} containing the prior definitions (if any suitable #'priors can be altered by the user) or \code{NULL}, indicating that no priors in the model #'can be modified through the `mvgam` interface diff --git a/R/hindcast.mvgam.R b/R/hindcast.mvgam.R index f6ecd968..cce003d0 100644 --- a/R/hindcast.mvgam.R +++ b/R/hindcast.mvgam.R @@ -1,23 +1,17 @@ #'@title Extract hindcasts for a fitted \code{mvgam} object #'@name hindcast.mvgam #'@importFrom stats predict -#'@param object \code{list} object returned from \code{mvgam}. See [mvgam()] +#'@inheritParams predict.mvgam #'@param series Either a \code{integer} specifying which series in the set is to be forecast, #'or the character string \code{'all'}, specifying that all series should be forecast. This is preferable #'if the fitted model contained multivariate trends (either as a dynamic factor or \code{VAR} process), #'as it saves recomputing the full set of trends for each series individually -#'@param type When this has the value \code{link} the linear predictor is calculated on the link scale. -#'If \code{expected} is used, predictions reflect the expectation of the response (the mean) -#'but ignore uncertainty in the observation process. When \code{response} (default) is used, -#'the predictions take uncertainty in the observation process into account to return -#'predictions on the outcome scale. When \code{trend} is used, only the hindcast distribution for the -#'latent trend is returned #'@param ... Ignored #'@details Posterior retrodictions are drawn from the fitted \code{mvgam} and #'organized into a convenient format #'@return An object of class \code{mvgam_forecast} containing hindcast distributions. #'See \code{\link{mvgam_forecast-class}} for details. -#' +#'#'@seealso \code{\link{forecast.mvgam}} #'@export hindcast <- function(object, ...){ UseMethod("hindcast", object) diff --git a/R/index-mvgam.R b/R/index-mvgam.R index 0d0b1469..0b2c77df 100644 --- a/R/index-mvgam.R +++ b/R/index-mvgam.R @@ -12,7 +12,7 @@ NULL #' @rdname index-mvgam #' @importFrom posterior variables -#' @param x \code{list} object of class `mvgam` +#' @param x \code{list} object returned from \code{mvgam}. See [mvgam()] #' @method variables mvgam #' @export #' @export variables diff --git a/R/mcmc_plot.mvgam.R b/R/mcmc_plot.mvgam.R index 6667bc82..3b5331a7 100644 --- a/R/mcmc_plot.mvgam.R +++ b/R/mcmc_plot.mvgam.R @@ -9,15 +9,18 @@ #' @param type The type of the plot. #' Supported types are (as names) \code{hist}, \code{dens}, #' \code{hist_by_chain}, \code{dens_overlay}, -#' \code{violin}, \code{intervals}, \code{areas}, \code{acf}, -#' \code{acf_bar},\code{trace}, \code{trace_highlight}, \code{scatter}, +#' \code{violin}, \code{intervals}, \code{areas}, +#' \code{areas_ridges}, \code{combo}, \code{acf}, +#' \code{acf_bar}, \code{trace}, \code{trace_highlight}, +#' \code{scatter}, \code{hex}, \code{pairs}, \code{violin}, #' \code{rhat}, \code{rhat_hist}, \code{neff}, \code{neff_hist} #' and \code{nuts_energy}. #' For an overview on the various plot types see #' \code{\link[bayesplot:MCMC-overview]{MCMC-overview}}. #' @return A \code{\link[ggplot2:ggplot]{ggplot}} object #' that can be further customized using the \pkg{ggplot2} package. -#' +#' @seealso \code{\link{mvgam_draws}} for an overview of some of the shortcut strings +#' that can be used for argument `variable` #' @examples #' \dontrun{ #' simdat <- sim_mvgam(n_series = 1, trend_model = 'AR1') @@ -26,6 +29,8 @@ #' data = simdat$data_train) #' mcmc_plot(mod) #' mcmc_plot(mod, type = 'neff_hist') +#' mcmc_plot(mod, variable = 'betas', type = 'areas') +#' mcmc_plot(mod, variable = 'trend_params', type = 'combo') #' } #' @export mcmc_plot.mvgam = function(object, @@ -75,13 +80,21 @@ mcmc_plot.mvgam = function(object, # x refers to a data.frame of draws draws <- as.array(object, variable = variable, regex = regex, use_alias = use_alias) - sel_variables <- dimnames(draws)[[2]] + sel_variables <- dimnames(draws)$variable if(type %in% c("scatter", "hex") && length(sel_variables) != 2L){ stop("Exactly 2 parameters must be selected for this type.", "\nParameters selected: ", paste0("'", sel_variables, "'", collapse = ", "), call. = FALSE) } + + if(type == 'pairs' && length(sel_variables) == 1L){ + stop("2 or more parameters must be selected for this type.", + "\nParameters selected: ", + paste0("'", sel_variables, "'", collapse = ", "), + call. = FALSE) + } + mcmc_args$x <- draws } } diff --git a/R/mvgam.R b/R/mvgam.R index 08ad7f97..4c9d920b 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -63,6 +63,7 @@ #' \item`student_t()` for real-valued data #' \item`Gamma()` for non-negative real-valued data} #'Note that only `nb()` and `poisson()` are available if using `JAGS` as the backend. +#'Default is `poisson()`. #'See \code{\link{mvgam_families}} for more details #'@param use_lv \code{logical}. If \code{TRUE}, use dynamic factors to estimate series' #'latent trends in a reduced dimension format. Defaults to \code{FALSE} diff --git a/R/plot.mvgam.R b/R/plot.mvgam.R index 7e6f913f..3195bbfc 100644 --- a/R/plot.mvgam.R +++ b/R/plot.mvgam.R @@ -2,7 +2,6 @@ #' #'This function takes a fitted \code{mvgam} object and produces plots of smooth functions, forecasts, trends and #'uncertainty components -#' #'@param x \code{list} object returned from \code{mvgam}. See [mvgam()] #'@param type \code{character} specifying which type of plot to return. Options are: #'series, @@ -14,7 +13,7 @@ #'trend, #'uncertainty, #'factors -#'@param residuals \code{logical}. If \code{TRUE} and `type = residuals`, posterior quantiles of partial residuals are added +#'@param residuals \code{logical}. If \code{TRUE} and `type = 'smooths'`, posterior quantiles of partial residuals are added #'to plots of 1-D smooths as a series of ribbon rectangles. #'Partial residuals for a smooth term are the median Dunn-Smyth residuals that would be obtained by dropping the term #'concerned from the model, while leaving all other estimates fixed (i.e. the diff --git a/R/plot_mvgam_factors.R b/R/plot_mvgam_factors.R index baa0061f..be035c72 100644 --- a/R/plot_mvgam_factors.R +++ b/R/plot_mvgam_factors.R @@ -3,7 +3,7 @@ #'This function takes a fitted \code{mvgam} object and returns plots and summary statistics for #'the latent dynamic factors #' -#'@param object \code{list} object returned from \code{mvgam} +#'@inheritParams plot.mvgam #'@param plot \code{logical} specifying whether factors should be plotted #'@author Nicholas J Clark #'@details If the model in \code{object} was estimated using dynamic factors, it is possible that not all factors diff --git a/R/plot_mvgam_fc.R b/R/plot_mvgam_fc.R index e8287cdd..45729449 100644 --- a/R/plot_mvgam_fc.R +++ b/R/plot_mvgam_fc.R @@ -1,6 +1,6 @@ #'Plot mvgam posterior predictions for a specified series #'@importFrom stats formula terms -#'@param object \code{list} object returned from \code{mvgam} +#'@inheritParams plot.mvgam #'@param series \code{integer} specifying which series in the set is to be plotted #'@param newdata Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' #'in addition to any other variables included in the linear predictor of the original \code{formula}. If included, the diff --git a/R/plot_mvgam_pterms.R b/R/plot_mvgam_pterms.R index b59716df..6feebec6 100644 --- a/R/plot_mvgam_pterms.R +++ b/R/plot_mvgam_pterms.R @@ -4,9 +4,7 @@ #' #'@importFrom graphics layout title rug bxp #'@importFrom stats coef predict -#'@param object \code{list} object returned from \code{mvgam} -#'@param trend_effects logical. If `TRUE` and a `trend_formula` was used in model -#'fitting, terms from the trend (i.e. process) model will be plotted +#'@inheritParams plot.mvgam #'@details Posterior empirical quantiles of each parametric term's partial effect estimates #'(on the link scale) are calculated and visualised as ribbon plots. These effects can #'be interpreted as the partial effect that a parametric term contributes when all other diff --git a/R/plot_mvgam_randomeffects.R b/R/plot_mvgam_randomeffects.R index 419a4da7..3f05d65f 100644 --- a/R/plot_mvgam_randomeffects.R +++ b/R/plot_mvgam_randomeffects.R @@ -3,9 +3,7 @@ #'This function plots posterior empirical quantiles for random effect smooths (bs = re) #' #'@importFrom graphics layout title -#'@param object \code{list} object returned from \code{mvgam} -#'@param trend_effects logical. If `TRUE` and a `trend_formula` was used in model -#'fitting, terms from the trend (i.e. process) model will be plotted +#'@inheritParams plot.mvgam #'@details Posterior empirical quantiles of random effect coefficient estimates #'(on the link scale) are calculated and visualised as ribbon plots. #'Labels for coefficients are taken from the levels of the original factor variable diff --git a/R/plot_mvgam_resids.R b/R/plot_mvgam_resids.R index eff551f2..0318da36 100644 --- a/R/plot_mvgam_resids.R +++ b/R/plot_mvgam_resids.R @@ -5,7 +5,7 @@ #'@importFrom graphics layout title #'@importFrom stats complete.cases qqnorm qqline acf pacf na.pass #'@importFrom mgcv bam -#'@param object \code{list} object returned from \code{mvgam} +#'@inheritParams plot.mvgam #'@param series \code{integer} specifying which series in the set is to be plotted #'@param newdata Optional \code{dataframe} or \code{list} of test data containing at least 'series', 'y', and 'time' #'in addition to any other variables included in the linear predictor of \code{formula}. If included, the diff --git a/R/plot_mvgam_smooth.R b/R/plot_mvgam_smooth.R index 89827be7..2ce98030 100644 --- a/R/plot_mvgam_smooth.R +++ b/R/plot_mvgam_smooth.R @@ -4,9 +4,7 @@ #' #'@importFrom grDevices hcl.colors #'@importFrom stats quantile predict -#'@param object \code{list} object returned from \code{mvgam} -#'@param trend_effects logical. If `TRUE` and a `trend_formula` was used in model -#'fitting, terms from the trend (i.e. process) model will be plotted +#'@inheritParams plot.mvgam #'@param series \code{integer} specifying which series in the set is to be plotted #'@param smooth either a \code{character} or \code{integer} specifying which smooth term to be plotted #'@param residuals \code{logical}. If \code{TRUE} then posterior quantiles of partial residuals are added diff --git a/R/plot_mvgam_trend.R b/R/plot_mvgam_trend.R index 0abc4653..5fadbb6b 100644 --- a/R/plot_mvgam_trend.R +++ b/R/plot_mvgam_trend.R @@ -1,7 +1,7 @@ #'Plot mvgam latent trend for a specified series #'@importFrom graphics par lines polygon box abline #'@importFrom stats sd quantile -#'@param object \code{list} object returned from \code{mvgam} +#'@inheritParams plot.mvgam #'@param series \code{integer} specifying which series in the set is to be plotted #'@param newdata Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' #'in addition to any other variables included in the linear predictor of the original \code{formula}. diff --git a/R/plot_mvgam_uncertainty.R b/R/plot_mvgam_uncertainty.R index 693cfaca..2a095449 100644 --- a/R/plot_mvgam_uncertainty.R +++ b/R/plot_mvgam_uncertainty.R @@ -1,7 +1,7 @@ #'Plot mvgam forecast uncertainty contributions for a specified series #'@importFrom graphics legend #'@importFrom stats predict -#'@param object \code{list} object returned from \code{mvgam} +#'@inheritParams plot.mvgam #'@param series \code{integer} specifying which series in the set is to be plotted #'@param newdata A \code{dataframe} or \code{list} containing at least 'series' and 'time' for the forecast horizon, in #'addition to any other variables included in the linear predictor of \code{formula} diff --git a/R/ppc.mvgam.R b/R/ppc.mvgam.R index d1cb954d..47ddf108 100644 --- a/R/ppc.mvgam.R +++ b/R/ppc.mvgam.R @@ -39,7 +39,21 @@ #'series (for \code{type == 'hist'}), kernel density or empirical CDF estimates for #'posterior predictions (for \code{type == 'density'} or \code{type == 'cdf'}) or a Probability #'Integral Transform histogram (for \code{type == 'pit'}). +#' @examples +#' \dontrun{ +#' # Simulate some smooth effects and fit a model +#' set.seed(0) +#' dat <- mgcv::gamSim(1, n = 200, scale = 2) +#' dat$time <- 1:NROW(dat) +#' mod <- mvgam(y ~ s(x0) + s(x1) + s(x2) + s(x3), +#' data = dat, +#' family = gaussian()) #' +#' # Posterior checks +#' ppc(mod, type = 'hist') +#' ppc(mod, type = 'density') +#' ppc(mod, type = 'cdf') +#' } #'@export ppc <- function(object, ...){ UseMethod("ppc", object) diff --git a/R/predict.mvgam.R b/R/predict.mvgam.R index 2313b769..1d88c064 100644 --- a/R/predict.mvgam.R +++ b/R/predict.mvgam.R @@ -1,6 +1,6 @@ #'Predict from the GAM component of an mvgam model #'@importFrom stats predict -#'@param object Object of class `mvgam` +#'@param object \code{list} object returned from \code{mvgam}. See [mvgam()] #'@param newdata Optional \code{dataframe} or \code{list} of test data containing the #'variables included in the linear predictor of \code{formula}. If not supplied, #'predictions are generated for the original observations used for the model fit. diff --git a/R/score.mvgam_forecast.R b/R/score.mvgam_forecast.R index 3752b094..081e6c65 100644 --- a/R/score.mvgam_forecast.R +++ b/R/score.mvgam_forecast.R @@ -32,21 +32,23 @@ #'will only contain the linear predictors #'@examples #'\dontrun{ -#'#Simulate observations for three count-valued time series +#'# Simulate observations for three count-valued time series #'data <- sim_mvgam() -#'#Fit a dynamic model using 'newdata' to automatically produce forecasts +#'# Fit a dynamic model using 'newdata' to automatically produce forecasts #'mod <- mvgam(y ~ 1, #' trend_model = 'RW', #' data = data$data_train, #' newdata = data$data_test) #' -#'#Extract forecasts into a 'mvgam_forecast' object +#'# Extract forecasts into a 'mvgam_forecast' object #'fc <- forecast(mod) #' -#'#Score forecasts -#'score(fc, score = 'drps') +#'# Compute Discrete Rank Probability Scores and 0.90 interval coverages +#'fc_scores <- score(fc, score = 'drps') +#'str(fc_scores) #'} #'@method score mvgam_forecast +#'@seealso \code{\link{forecast.mvgam}} #'@export score.mvgam_forecast = function(object, score = 'crps', log = FALSE, weights, diff --git a/R/update.mvgam.R b/R/update.mvgam.R index f8ce3e2c..b61b003f 100644 --- a/R/update.mvgam.R +++ b/R/update.mvgam.R @@ -4,56 +4,11 @@ #'@name update.mvgam #'@importFrom mgcv nb betar #'@importFrom rlang missing_arg -#'@param object A fitted `mvgam` model +#'@inheritParams mvgam +#'@param object \code{list} object returned from \code{mvgam}. See [mvgam()] #'@param formula Optional new `formula` object. Note, `mvgam` currently does not support dynamic formula #'updates such as removal of specific terms with `- term`. When updating, the entire formula needs #'to be supplied -#'@param trend_formula An optional \code{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)`). This feature is experimental, and is only -#'currently available for Random Walk trend models. -#'@param data A \code{dataframe} or \code{list} containing the model response variable and covariates -#'required by the GAM \code{formula}. Should include columns: -#''series' (character or factor index of the series IDs) -#''time' (numeric index of the time point for each observation). -#'Any other variables to be included in the linear predictor of \code{formula} must also be present -#'@param newdata Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' -#'in addition to any other variables included in the linear predictor of \code{formula}. If included, the -#'observations in variable \code{y} will be set to \code{NA} when fitting the model so that posterior -#'simulations can be obtained -#'@param trend_model \code{character} specifying the time series dynamics for the latent trend. Options are: -#'\itemize{ -#' \item `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 by \code{\link[mgcv]{gam}}) -#' \item `RW` (random walk with possible drift) -#' \item `AR1` (with possible drift) -#' \item `AR2` (with possible drift) -#' \item `AR3` (with possible drift) -#' \item `VAR1` (with possible drift; only available in \code{Stan}) -#' \item `GP` (Gaussian Process with squared exponential kernel; -#'only available in \code{Stan})} -#'@param 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 setting -#'`use_lv = TRUE` and using the mapping to set up the shared trends. Needs to have column names -#'`series` and `trend`, with integer values in the `trend` column to state which trend each series -#'should depend on. The `series` column should have a single unique entry for each series in the -#'data (names should perfectly match factor levels of the `series` variable in `data`). See examples -#'in \code{\link{mvgam}} for details -#'@param family \code{family} specifying the exponential observation family for the series. Currently supported -#'families are: `nb()`, `poisson()`, `tweedie()`, `gaussian()`, `betar()`, `lognormal()`, `student_t()` and `Gamma()` -#'@param use_lv \code{logical}. If \code{TRUE}, use dynamic factors to estimate series' -#'latent trends in a reduced dimension format. If \code{FALSE}, estimate independent latent trends for each series -#'@param n_lv \code{integer} the number of latent dynamic factors to use if \code{use_lv == TRUE}. -#'Cannot be \code{>n_series}. Defaults arbitrarily to \code{min(2, floor(n_series / 2))} -#'@param priors An optional \code{data.frame} with prior -#'definitions. See \code{\link{get_mvgam_priors}} and -#'[mvgam] for more information on changing default prior distributions -#'@param 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` #'@param ... Other arguments to be passed to \code{\link{mvgam}} #' @examples #' \dontrun{ diff --git a/man/conditional_effects.mvgam.Rd b/man/conditional_effects.mvgam.Rd index 9b491b48..2f4da81b 100644 --- a/man/conditional_effects.mvgam.Rd +++ b/man/conditional_effects.mvgam.Rd @@ -101,6 +101,9 @@ conditional_effects(mod) conditional_effects(mod, conf_level = 0.5, type = 'link') } } +\seealso{ +\code{\link[marginaleffects]{plot_predictions}}, \code{\link[marginaleffects]{plot_slopes}} +} \author{ Nicholas J Clark } diff --git a/man/forecast.mvgam.Rd b/man/forecast.mvgam.Rd index 6e209f1d..3051e2d7 100644 --- a/man/forecast.mvgam.Rd +++ b/man/forecast.mvgam.Rd @@ -41,12 +41,12 @@ as it saves recomputing the full set of trends for each series individually} \item{n_cores}{\code{integer} specifying number of cores for generating forecasts in parallel} -\item{type}{When this has the value \code{link} the linear predictor is calculated on the link scale. +\item{type}{When this has the value \code{link} (default) the linear predictor is calculated on the link scale. If \code{expected} is used, predictions reflect the expectation of the response (the mean) -but ignore uncertainty in the observation process. When \code{response} (default) is used, +but ignore uncertainty in the observation process. When \code{response} is used, the predictions take uncertainty in the observation process into account to return -predictions on the outcome scale. When \code{trend} is used, only the forecast distribution for the -latent trend is returned} +predictions on the outcome scale. When \code{variance} is used, the variance of the response +with respect to the mean (mean-variance relationship) is returned} } \value{ An object of class \code{mvgam_forecast} containing hindcast and forecast distributions. @@ -87,3 +87,6 @@ plot(fc, series = 3) } } +\seealso{ +\code{\link{hindcast}}, \code{\link{score}} +} diff --git a/man/formula.mvgam.Rd b/man/formula.mvgam.Rd index cc8d2cf2..497051fc 100644 --- a/man/formula.mvgam.Rd +++ b/man/formula.mvgam.Rd @@ -3,7 +3,7 @@ \name{formula.mvgam} \alias{formula.mvgam} \alias{formula.mvgam_prefit} -\title{Extract model.frame from a fitted mvgam object} +\title{Extract formulae from mvgam objects} \usage{ \method{formula}{mvgam}(x, trend_effects = FALSE, ...) @@ -12,17 +12,17 @@ \arguments{ \item{x}{\code{mvgam} or \code{mvgam_prefit} object} -\item{trend_effects}{\code{logical}, return the model.frame from the +\item{trend_effects}{\code{logical}, return the formula from the observation model (if \code{FALSE}) or from the underlying process model (if\code{TRUE})} \item{...}{Ignored} } \value{ -A \code{matrix} containing the fitted model frame +A \code{formula} object } \description{ -Extract model.frame from a fitted mvgam object +Extract formulae from mvgam objects } \author{ Nicholas J Clark diff --git a/man/get_mvgam_priors.Rd b/man/get_mvgam_priors.Rd index 0a30f373..c261a995 100644 --- a/man/get_mvgam_priors.Rd +++ b/man/get_mvgam_priors.Rd @@ -54,6 +54,7 @@ families are: \item\code{student_t()} for real-valued data \item\code{Gamma()} for non-negative real-valued data} Note that only \code{nb()} and \code{poisson()} are available if using \code{JAGS} as the backend. +Default is \code{poisson()}. See \code{\link{mvgam_families}} for more details} \item{use_lv}{\code{logical}. If \code{TRUE}, use dynamic factors to estimate series' @@ -62,13 +63,11 @@ latent trends in a reduced dimension format. Defaults to \code{FALSE}} \item{n_lv}{\code{integer} the number of latent dynamic factors to use if \code{use_lv == TRUE}. Cannot be \code{>n_series}. Defaults arbitrarily to \code{min(2, floor(n_series / 2))}} -\item{use_stan}{Logical. If \code{TRUE} and if \code{rstan} is installed, the model will be compiled and sampled using -the Hamiltonian Monte Carlo with a call to \code{\link[cmdstanr]{cmdstan_model}} or, if \code{cmdstanr} is not available, -a call to \code{\link[rstan]{stan}}. Note that this functionality is still in development and -not all options that are available in \code{JAGS} can be used, including: no option for a Tweedie family and no option for -dynamic factor trends. However, as \code{Stan} can estimate Hilbert base approximate Gaussian Processes, which -are much more computationally tractable than full GPs for time series with \verb{>100} observations, estimation -in \code{Stan} can support latent GP trends while estimation in \code{JAGS} cannot} +\item{use_stan}{Logical. If \code{TRUE}, the model will be compiled and sampled using +Hamiltonian Monte Carlo with a call to \code{\link[cmdstanr]{cmdstan_model}} or +a call to \code{\link[rstan]{stan}}. Note that +there are many more options when using \code{Stan} vs \code{JAGS} (the only "advantage" of \code{JAGS} is the ability +to use a Tweedie family).} \item{trend_model}{\code{character} or \code{function} specifying the time series dynamics for the latent trend. Options are: \itemize{ @@ -241,7 +240,7 @@ code(mod2) } \seealso{ -\code{\link{mvgam}} +\code{\link{mvgam}} \code{\link[brms]{prior}} } \author{ Nicholas J Clark diff --git a/man/hindcast.mvgam.Rd b/man/hindcast.mvgam.Rd index 10e551eb..588c00d5 100644 --- a/man/hindcast.mvgam.Rd +++ b/man/hindcast.mvgam.Rd @@ -19,16 +19,17 @@ or the character string \code{'all'}, specifying that all series should be forec if the fitted model contained multivariate trends (either as a dynamic factor or \code{VAR} process), as it saves recomputing the full set of trends for each series individually} -\item{type}{When this has the value \code{link} the linear predictor is calculated on the link scale. +\item{type}{When this has the value \code{link} (default) the linear predictor is calculated on the link scale. If \code{expected} is used, predictions reflect the expectation of the response (the mean) -but ignore uncertainty in the observation process. When \code{response} (default) is used, +but ignore uncertainty in the observation process. When \code{response} is used, the predictions take uncertainty in the observation process into account to return -predictions on the outcome scale. When \code{trend} is used, only the hindcast distribution for the -latent trend is returned} +predictions on the outcome scale. When \code{variance} is used, the variance of the response +with respect to the mean (mean-variance relationship) is returned} } \value{ An object of class \code{mvgam_forecast} containing hindcast distributions. See \code{\link{mvgam_forecast-class}} for details. +#'@seealso \code{\link{forecast.mvgam}} } \description{ Extract hindcasts for a fitted \code{mvgam} object diff --git a/man/index-mvgam.Rd b/man/index-mvgam.Rd index 6d38c66e..1f5d8473 100644 --- a/man/index-mvgam.Rd +++ b/man/index-mvgam.Rd @@ -15,7 +15,7 @@ \method{variables}{mvgam}(x, ...) } \arguments{ -\item{x}{\code{list} object of class \code{mvgam}} +\item{x}{\code{list} object returned from \code{mvgam}. See \code{\link[=mvgam]{mvgam()}}} \item{...}{Arguments passed to individual methods (if applicable).} } diff --git a/man/mcmc_plot.mvgam.Rd b/man/mcmc_plot.mvgam.Rd index 42196808..7e178f76 100644 --- a/man/mcmc_plot.mvgam.Rd +++ b/man/mcmc_plot.mvgam.Rd @@ -19,8 +19,10 @@ \item{type}{The type of the plot. Supported types are (as names) \code{hist}, \code{dens}, \code{hist_by_chain}, \code{dens_overlay}, -\code{violin}, \code{intervals}, \code{areas}, \code{acf}, -\code{acf_bar},\code{trace}, \code{trace_highlight}, \code{scatter}, +\code{violin}, \code{intervals}, \code{areas}, +\code{areas_ridges}, \code{combo}, \code{acf}, +\code{acf_bar}, \code{trace}, \code{trace_highlight}, +\code{scatter}, \code{violin}, \code{rhat}, \code{rhat_hist}, \code{neff}, \code{neff_hist} and \code{nuts_energy}. For an overview on the various plot types see @@ -57,5 +59,11 @@ mod <- mvgam(y ~ s(season, bs = 'cc'), data = simdat$data_train) mcmc_plot(mod) mcmc_plot(mod, type = 'neff_hist') +mcmc_plot(mod, variable = 'betas', type = 'areas') +mcmc_plot(mod, variable = 'trend_params', type = 'combo') } } +\seealso{ +\code{\link{mvgam_draws}} for an overview of some of the shortcut strings +that can be used for argument \code{variable} +} diff --git a/man/mvgam.Rd b/man/mvgam.Rd index 94a26ff6..3d4da0e7 100644 --- a/man/mvgam.Rd +++ b/man/mvgam.Rd @@ -106,6 +106,7 @@ families are: \item\code{student_t()} for real-valued data \item\code{Gamma()} for non-negative real-valued data} Note that only \code{nb()} and \code{poisson()} are available if using \code{JAGS} as the backend. +Default is \code{poisson()}. See \code{\link{mvgam_families}} for more details} \item{use_lv}{\code{logical}. If \code{TRUE}, use dynamic factors to estimate series' diff --git a/man/plot.mvgam.Rd b/man/plot.mvgam.Rd index d0e0d2be..6c34742d 100644 --- a/man/plot.mvgam.Rd +++ b/man/plot.mvgam.Rd @@ -32,7 +32,7 @@ factors} \item{series}{\code{integer} specifying which series in the set is to be plotted. This is ignored if \code{type == 're'}} -\item{residuals}{\code{logical}. If \code{TRUE} and \code{type = residuals}, posterior quantiles of partial residuals are added +\item{residuals}{\code{logical}. If \code{TRUE} and \code{type = 'smooths'}, posterior quantiles of partial residuals are added to plots of 1-D smooths as a series of ribbon rectangles. Partial residuals for a smooth term are the median Dunn-Smyth residuals that would be obtained by dropping the term concerned from the model, while leaving all other estimates fixed (i.e. the diff --git a/man/plot_mvgam_factors.Rd b/man/plot_mvgam_factors.Rd index 859211de..d9f81f1e 100644 --- a/man/plot_mvgam_factors.Rd +++ b/man/plot_mvgam_factors.Rd @@ -7,8 +7,6 @@ plot_mvgam_factors(object, plot = TRUE) } \arguments{ -\item{object}{\code{list} object returned from \code{mvgam}} - \item{plot}{\code{logical} specifying whether factors should be plotted} } \value{ diff --git a/man/plot_mvgam_forecasts.Rd b/man/plot_mvgam_forecasts.Rd index e389682c..81bfe154 100644 --- a/man/plot_mvgam_forecasts.Rd +++ b/man/plot_mvgam_forecasts.Rd @@ -37,8 +37,6 @@ plot_mvgam_fc( ) } \arguments{ -\item{object}{\code{list} object returned from \code{mvgam}} - \item{series}{\code{integer} specifying which series in the set is to be plotted} \item{newdata}{Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' diff --git a/man/plot_mvgam_pterms.Rd b/man/plot_mvgam_pterms.Rd index 73ba2f9a..1799f82d 100644 --- a/man/plot_mvgam_pterms.Rd +++ b/man/plot_mvgam_pterms.Rd @@ -7,8 +7,6 @@ plot_mvgam_pterms(object, trend_effects = FALSE) } \arguments{ -\item{object}{\code{list} object returned from \code{mvgam}} - \item{trend_effects}{logical. If \code{TRUE} and a \code{trend_formula} was used in model fitting, terms from the trend (i.e. process) model will be plotted} } diff --git a/man/plot_mvgam_randomeffects.Rd b/man/plot_mvgam_randomeffects.Rd index 9fc4dfc1..9210005f 100644 --- a/man/plot_mvgam_randomeffects.Rd +++ b/man/plot_mvgam_randomeffects.Rd @@ -7,8 +7,6 @@ plot_mvgam_randomeffects(object, trend_effects = FALSE) } \arguments{ -\item{object}{\code{list} object returned from \code{mvgam}} - \item{trend_effects}{logical. If \code{TRUE} and a \code{trend_formula} was used in model fitting, terms from the trend (i.e. process) model will be plotted} } diff --git a/man/plot_mvgam_resids.Rd b/man/plot_mvgam_resids.Rd index 0d2d2d83..22503ad3 100644 --- a/man/plot_mvgam_resids.Rd +++ b/man/plot_mvgam_resids.Rd @@ -7,8 +7,6 @@ plot_mvgam_resids(object, series = 1, newdata, data_test) } \arguments{ -\item{object}{\code{list} object returned from \code{mvgam}} - \item{series}{\code{integer} specifying which series in the set is to be plotted} \item{newdata}{Optional \code{dataframe} or \code{list} of test data containing at least 'series', 'y', and 'time' diff --git a/man/plot_mvgam_smooth.Rd b/man/plot_mvgam_smooth.Rd index 65900446..6bf28bf6 100644 --- a/man/plot_mvgam_smooth.Rd +++ b/man/plot_mvgam_smooth.Rd @@ -18,8 +18,6 @@ plot_mvgam_smooth( ) } \arguments{ -\item{object}{\code{list} object returned from \code{mvgam}} - \item{trend_effects}{logical. If \code{TRUE} and a \code{trend_formula} was used in model fitting, terms from the trend (i.e. process) model will be plotted} diff --git a/man/plot_mvgam_trend.Rd b/man/plot_mvgam_trend.Rd index 36b9012e..900b09ce 100644 --- a/man/plot_mvgam_trend.Rd +++ b/man/plot_mvgam_trend.Rd @@ -20,8 +20,6 @@ plot_mvgam_trend( ) } \arguments{ -\item{object}{\code{list} object returned from \code{mvgam}} - \item{series}{\code{integer} specifying which series in the set is to be plotted} \item{newdata}{Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' diff --git a/man/plot_mvgam_uncertainty.Rd b/man/plot_mvgam_uncertainty.Rd index d7824ba6..8e217ee7 100644 --- a/man/plot_mvgam_uncertainty.Rd +++ b/man/plot_mvgam_uncertainty.Rd @@ -14,8 +14,6 @@ plot_mvgam_uncertainty( ) } \arguments{ -\item{object}{\code{list} object returned from \code{mvgam}} - \item{series}{\code{integer} specifying which series in the set is to be plotted} \item{newdata}{A \code{dataframe} or \code{list} containing at least 'series' and 'time' for the forecast horizon, in diff --git a/man/posterior_epred.mvgam.Rd b/man/posterior_epred.mvgam.Rd index c93368c6..1911f035 100644 --- a/man/posterior_epred.mvgam.Rd +++ b/man/posterior_epred.mvgam.Rd @@ -8,7 +8,7 @@ \method{posterior_epred}{mvgam}(object, newdata, data_test, process_error = TRUE, ...) } \arguments{ -\item{object}{Object of class \code{mvgam}} +\item{object}{\code{list} object returned from \code{mvgam}. See \code{\link[=mvgam]{mvgam()}}} \item{newdata}{Optional \code{dataframe} or \code{list} of test data containing the variables included in the linear predictor of \code{formula}. If not supplied, diff --git a/man/posterior_linpred.mvgam.Rd b/man/posterior_linpred.mvgam.Rd index 672f734b..dd125dcb 100644 --- a/man/posterior_linpred.mvgam.Rd +++ b/man/posterior_linpred.mvgam.Rd @@ -14,7 +14,7 @@ ) } \arguments{ -\item{object}{Object of class \code{mvgam}} +\item{object}{\code{list} object returned from \code{mvgam}. See \code{\link[=mvgam]{mvgam()}}} \item{transform}{Logical; if \code{FALSE} (the default), draws of the linear predictor are returned. diff --git a/man/posterior_predict.mvgam.Rd b/man/posterior_predict.mvgam.Rd index b765b6fc..0c437971 100644 --- a/man/posterior_predict.mvgam.Rd +++ b/man/posterior_predict.mvgam.Rd @@ -7,7 +7,7 @@ \method{posterior_predict}{mvgam}(object, newdata, data_test, process_error = TRUE, ...) } \arguments{ -\item{object}{Object of class \code{mvgam}} +\item{object}{\code{list} object returned from \code{mvgam}. See \code{\link[=mvgam]{mvgam()}}} \item{newdata}{Optional \code{dataframe} or \code{list} of test data containing the variables included in the linear predictor of \code{formula}. If not supplied, diff --git a/man/ppc.mvgam.Rd b/man/ppc.mvgam.Rd index 139a8d84..ca29054c 100644 --- a/man/ppc.mvgam.Rd +++ b/man/ppc.mvgam.Rd @@ -74,3 +74,19 @@ can also be compared to out of sample observations as long as these observations 'data_test' in the original model fit and supplied here. Rootograms are currently only plotted using the 'hanging' style } +\examples{ +\dontrun{ +# Simulate some smooth effects and fit a model +set.seed(0) +dat <- mgcv::gamSim(1, n = 200, scale = 2) +dat$time <- 1:NROW(dat) +mod <- mvgam(y ~ s(x0) + s(x1) + s(x2) + s(x3), + data = dat, + family = gaussian()) + +# Posterior checks +ppc(mod, type = 'hist') +ppc(mod, type = 'density') +ppc(mod, type = 'cdf') +} +} diff --git a/man/predict.mvgam.Rd b/man/predict.mvgam.Rd index 45f7535b..6a5265c0 100644 --- a/man/predict.mvgam.Rd +++ b/man/predict.mvgam.Rd @@ -7,7 +7,7 @@ \method{predict}{mvgam}(object, newdata, data_test, type = "link", process_error = TRUE, ...) } \arguments{ -\item{object}{Object of class \code{mvgam}} +\item{object}{\code{list} object returned from \code{mvgam}. See \code{\link[=mvgam]{mvgam()}}} \item{newdata}{Optional \code{dataframe} or \code{list} of test data containing the variables included in the linear predictor of \code{formula}. If not supplied, diff --git a/man/score.mvgam_forecast.Rd b/man/score.mvgam_forecast.Rd index a3a22cbb..8dedb3fd 100644 --- a/man/score.mvgam_forecast.Rd +++ b/man/score.mvgam_forecast.Rd @@ -59,18 +59,22 @@ Compute probabilistic forecast scores for mvgam objects } \examples{ \dontrun{ -#Simulate observations for three count-valued time series +# Simulate observations for three count-valued time series data <- sim_mvgam() -#Fit a dynamic model using 'newdata' to automatically produce forecasts +# Fit a dynamic model using 'newdata' to automatically produce forecasts mod <- mvgam(y ~ 1, trend_model = 'RW', data = data$data_train, newdata = data$data_test) -#Extract forecasts into a 'mvgam_forecast' object +# Extract forecasts into a 'mvgam_forecast' object fc <- forecast(mod) -#Score forecasts -score(fc, score = 'drps') +# Compute Discrete Rank Probability Scores and 0.90 interval coverages +fc_scores <- score(fc, score = 'drps') +str(fc_scores) } } +\seealso{ +\code{\link{forecast.mvgam}} +} diff --git a/man/update.mvgam.Rd b/man/update.mvgam.Rd index 85dbed0f..532918da 100644 --- a/man/update.mvgam.Rd +++ b/man/update.mvgam.Rd @@ -21,7 +21,7 @@ ) } \arguments{ -\item{object}{A fitted \code{mvgam} model} +\item{object}{\code{list} object returned from \code{mvgam}. See \code{\link[=mvgam]{mvgam()}}} \item{formula}{Optional new \code{formula} object. Note, \code{mvgam} currently does not support dynamic formula updates such as removal of specific terms with \code{- term}. When updating, the entire formula needs @@ -30,31 +30,39 @@ to be supplied} \item{trend_formula}{An optional \code{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 \code{~ season + s(year)}). This feature is experimental, and is only -currently available for Random Walk trend models.} +of the formula (i.e. a valid option would be \code{~ season + s(year)}). Also note that you should not use +the identifier \code{series} in this formula to specify effects that vary across time series. Instead you should use +\code{trend}. This will ensure that models in which a \code{trend_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 experimental, and is only currently available for Random Walk and AR trend models.} \item{data}{A \code{dataframe} or \code{list} containing the model response variable and covariates required by the GAM \code{formula}. Should include columns: -'series' (character or factor index of the series IDs) -'time' (numeric index of the time point for each observation). +\code{series} (a \code{factor} index of the series IDs;the number of levels should be identical +to the number of unique series labels (i.e. \code{n_series = length(levels(data$series))})) +\code{time} (\code{numeric} or \code{integer} index of the time point for each observation). Any other variables to be included in the linear predictor of \code{formula} must also be present} -\item{newdata}{Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' +\item{newdata}{Optional \code{dataframe} or \code{list} of test data containing at least \code{series} and \code{time} in addition to any other variables included in the linear predictor of \code{formula}. If included, the observations in variable \code{y} will be set to \code{NA} when fitting the model so that posterior simulations can be obtained} -\item{trend_model}{\code{character} specifying the time series dynamics for the latent trend. Options are: +\item{trend_model}{\code{character} or \code{function} specifying the time series dynamics for the latent trend. Options are: \itemize{ \item \code{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 by \code{\link[mgcv]{gam}}) -\item \code{RW} (random walk with possible drift) -\item \code{AR1} (with possible drift) -\item \code{AR2} (with possible drift) -\item \code{AR3} (with possible drift) -\item \code{VAR1} (with possible drift; only available in \code{Stan}) -\item \code{GP} (Gaussian Process with squared exponential kernel; -only available in \code{Stan})}} +\item \code{'RW'} or \code{RW()} +\item \code{'AR1'} or \code{AR(p = 1)} +\item \code{'AR2'} or \code{AR(p = 2)} +\item \code{'AR3'} or \code{AR(p = 3)} +\item \code{'VAR1'} or \code{VAR()}(only available in \code{Stan}) +\item \code{'GP'} (Gaussian Process with squared exponential kernel; +only available in \code{Stan})} + +For all types apart from \code{'GP'}, moving average and/or correlated +process error terms can also be estimated (for example, \code{RW(cor = TRUE)} will set up a +multivariate Random Walk if \code{n_series > 1}). See \link{mvgam_trends} for more details} \item{trend_map}{Optional \code{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 @@ -63,20 +71,32 @@ different observation processes. If supplied, a latent factor model is set up by \code{series} and \code{trend}, with integer values in the \code{trend} column to state which trend each series should depend on. The \code{series} column should have a single unique entry for each series in the data (names should perfectly match factor levels of the \code{series} variable in \code{data}). See examples -in \code{\link{mvgam}} for details} +for details} \item{use_lv}{\code{logical}. If \code{TRUE}, use dynamic factors to estimate series' -latent trends in a reduced dimension format. If \code{FALSE}, estimate independent latent trends for each series} +latent trends in a reduced dimension format. Defaults to \code{FALSE}} \item{n_lv}{\code{integer} the number of latent dynamic factors to use if \code{use_lv == TRUE}. Cannot be \code{>n_series}. Defaults arbitrarily to \code{min(2, floor(n_series / 2))}} \item{family}{\code{family} specifying the exponential observation family for the series. Currently supported -families are: \code{nb()}, \code{poisson()}, \code{tweedie()}, \code{gaussian()}, \code{betar()}, \code{lognormal()}, \code{student_t()} and \code{Gamma()}} +families are: +\itemize{ +\item\code{nb()} for count data +\item\code{poisson()} for count data +\item\code{gaussian()} for real-valued data +\item\code{betar()} for proportional data on \verb{(0,1)} +\item\code{lognormal()} for non-negative real-valued data +\item\code{student_t()} for real-valued data +\item\code{Gamma()} for non-negative real-valued data} +Note that only \code{nb()} and \code{poisson()} are available if using \code{JAGS} as the backend. +Default is \code{poisson()}. +See \code{\link{mvgam_families}} for more details} \item{priors}{An optional \code{data.frame} with prior -definitions. See \code{\link{get_mvgam_priors}} and -\link{mvgam} for more information on changing default prior distributions} +definitions (in JAGS or Stan syntax). if using Stan, this can also be an object of +class \code{brmsprior} (see. \code{\link[brms]{prior}} for details). See \link{get_mvgam_priors} and +'Details' for more information on changing default prior distributions} \item{lfo}{Logical indicating whether this is part of a call to \link{lfo_cv.mvgam}. Returns a lighter version of the model with no residuals and fewer monitored parameters to speed up diff --git a/src/mvgam.dll b/src/mvgam.dll index b2557de620a145f8bddf6a472a568fe59d629dce..ab047ba9f0cce1a8c53283e6d7dbc94693dd2309 100644 GIT binary patch delta 75 zcmWN|yAeP@07cQUSbzG^fOs2_Yd{?%G fkuV|_iAaSNnaD*UN>PbgG@`xe`sbgweY-GyBD5mX delta 75 zcmV~$y$wJ>0EXc>oPW26Ep7q5ZBFqeB9X9$UNkK7_6(u&JR$fH{5qN2ZF6syeMCgU eh*%^d6;@;-7lkNAC2G-#_P}+Ix!h@A_YFT2=pujs diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index 8dff305e56f3cdcb39e7e760b284b565edfdc51f..072501ea0892d373cc90f7ccbefe9e6c26f755d5 100644 GIT binary patch delta 9606 zcmZXY2{@F0_x~*+N=il9CKVBpU4%+e){tdvBV*r2)-2RaI=wu{rNZoj=#Qx_G zv`Hf$jf6yuJoRNjI}&k2Yh*BI>zXZF-;dqSk&dHLDm2t1tYN1tyEkgzJbnI7<>vw% z0+Gk*D`aq@pP5&Gp4pOn;=(HLrR5SgSP1*INg1rE*&{BD!6pN0b{tE|qcH3SSs&&c zN7jJ>yQEEWq0+1}qGqeLaDE}R#*fnA4Ukba{(JMeWx&E-Wr5-Tc7Nf%vD)C%y%_>x zf5H>HKbKz}aP1)Az+fY25Qc@S{~BJ-t-*bnFrW4lEd-EQ+|I@=pRwE145X*2>o+Z| zeA7C16b)|P5?AON-Wy$h8ffJkXuM4eZH5Jem+_u$q-;%BE&D*mkn1nYV7vF?WcD@@ zIvZE@0boC2@6H3!Wk)rZRtosJ`)fP+)Y8-%-2NUR6$rNpc(*aUT$Sjt*WKg!hTN}F zxobB!Fzt6f3Kvk~DoWflJ{b|eP4WY3tOLi^e(VmGi5#4|aiyldZ5P6<+QawBh~dmt zlfl5xo14oxf6ApwzJs&*tqNmTRMusbC`n2{07bioGVOh^iU78jHVLUU!<#$!n!VM{ zkuqSn;%LNRX8@V#eP9%@&r(LLG-;}j$6DKgX2+b z-nf--!nww9Kbv8=e~j!CXcb_=j(o5h_;xPnWl#6ZjoO2sx9VqW9y0Yjt30eG0WT- zg`^Vn@K0@S400~-Q;)+b_&x@u_4O(dbG7xrMb5efc^)dFr?Ut}q~&zHyXa_~{MF>G z1BUsXuyPh6zo29X=aRK>3#q%;2&}i$KV2eM+=A1`SIs7!u7qnXoO3Bb3Ya474JBbDF=oFubni*FI^NWLp&%0 z_*|U+LX4#!&Qh|^=26r+W+Y`?+Jy7s2%B92GZdRv6&|v&e~0*WAY|i&d(uF~;wK@g zSrNTACFdiqx=ZvAS}%=IN(QY;D6%^xOPxuTnzl_y-`)=#Z`k|*f=XHU6v0lsa#9hY z4!df@GG>El;J-44F(CBV3DmriFdl{^IZG#M#*e?VElFx)1OMUtTpcOa;bfoDSl!{| zpgYE4i_wP)yC>s6_76FI!4J>x;aU^g&AhJ?+q@)AIvu)PlpzIdM7B`|L%kUM014`> zsbRWcE9Ky42q=Y>j((k=u^e!L?qLRUaoJYK5<5D4oaid0BCpLkWDKIEoZWZr!xD$Q zimh#O6Z1M!oD$Z=0Q-gOkR|xBq-u!K>EV{oHN)tvd}hDYl2ZNfURjLM%V5m*+WzR| zkXtef9AYx&k$45xWvyf|y?RzNv0?XvsQq(HEH0ZF7;N&_i78L(n9TiLe)EkC5*oMZ zh->#x)_PvW1sk$TRWL4Z;WxgGEN-s6f1P96re77ldeFL1)_gBaW^{gU=9f~8?SKnu zGND1LUf;;Q-!|9JCgFuel`R_)E%!r$aPeje)Yg{NR{t`8iqP}GBBgsYdpN*|f!yXk zwCK?Z2ow9K5*#p#@VPRlmsj#WKmXLEG0B>be}fNU6yvN z!yD&zdimO{GG7Dzv}1N=6Rg#(agixd`Zp;Ym^Thx~$2KoVLM5&I{D(;6&tK+xz5(EAG3B>^oTX`GN z{Ny$C>2#@Yj+HTp6EKH?l9rPCjz(m08!F(+*bmtnyL``cdD!OFp%RP@(Gl3O;bvJp zK6fn^K}?eQ#KClVEN@a1<~+%}w=A`7(+$Z@8r2~4Y`CD^ukrcV>J-13O3uPVBLJEH z?bl>wMDAob`8L9|fy3hFeVJ`K_u83$lX$Or$DW#j81Hy%xoAZ>W`mwG5$G6kBd0|i zZ!fCbUdezH z9A+@n+hAF-)4yBzZ97B96|@otKpj5niS`~?#rFD}ZNw|4-U~3Z$Bh*f<})~ZWo`CI zvf`WNTPZcY=H->xDc?tk;Qb2kfe+1UqfkrH_34cc@#9Z{rDW@KUP`2lmZjEV3BQbW zc-H;2@Yd;6zZpd$veG*qb`D1F5YKN(t}ZSyEK=wr6QaF&emM4~0P`+KtCYdln7!9J z=9wQOe#qx42sWpvOAVQUZfwY_46aWv`Plj}PXn{b4JkhaS^KufqP+{m$J9t2JI0P{ z9XQH|8a;~L+pKNneh4XQbUUL$(stND%6T@~I^yvcMIT!+d=K)|Pyqk3@VZaLN{5_T zH@B~CLi?A30!`rU=)i}3jo}0`hhZ_EIJ%u(wo=B67*xp(-b9b~%C}}1&p9tFmX;@G z@)n==Lx3lpyYA;NnRUL+hTPFlHh#~{QPwe){A({i!6rtdUZ|-(*Tm^(u}-*AXp*N? zU1fzg4k@a+eDjB2k9?~MrvKiua3?vXdObpR*dR#~4&)|`LPKpslQyLWz~uNoFJTug z)=@I1Jum&u#^qX{UD;H!6pL7lH&2^-xp(|XYn_&2l^4%vJ^K^;5xSjlXi|dd-M@ zt~Vgu#5pG!h+9hb+*Moa5D!;jN)G!drkD45v^;wVn$zAg3WX_$21~K@ElnlcdLlFs zw-y#hp@~5o{e>85|+YH90yBI&H0-g)L~W>9l^ zMqCi=y3>5J^_9ldX8A$~K-DzwVae!~-273xvz5}g}z@zsL@ZxpJ_K0nk%DAOZ9 zVB4B?%U0<5wU_f^a5>!{mxwovD+k*$1&il8(8U~CH#w(#Y|M?P6dV&Zao8uz{t~1W zyufB@399d24BR$+#INPd++uZ)!!5uq`vH6pq=r3Pf(pzzc#jP93lyVsuqBO+LVh#N z8p(dPVl-@=y&QdwBQ*I9t_x3^9!IRsCI6C&m`Xk+G}3IWKaCtol!|n=BYwap!)@`o zW6j3%6mGmcW!}xvPNOJoSiuhUzDlvDd5Suk}w_1BS@)5k^2+43jySXavf#brmdb?iBy;;fY z%j(Of6^bJu!AiDofga4SFFpQIHb5K24Oi0O&fbwzH&O{kuZNKz_diqo zP|t?5YV7?nj<&dM|G`ZZ+!hVIL*UG(kewb{ANEEBi}Ndrk~f-SP3v9L|r zY~Ea|;4EX1uNGtBlgD1c7QZu-x6Z_VN|DbM<7p$F^%i{!ap{~Vjz0A(pd2ag&7LS)* zy+Uzu73)H&ic?&yQm-YLp@hzjGM?BOUI--n1-TZ zURL;Lta*D1Sx{7HtOsh!${LtgoxXW3TvSj7iksS45zngeRO^pj>Jjvx0r;-**OZh@ zHHZ&|_$LR7 ztYl}}Bl8Wl0B@e`iU(|+ma0vX_@3}8{+K0^?`B)m_%Br}Nz@;!jx!ILxCye;p7Fj> zFM6YPfXa~{)vvIpj^gJ>15<$|tV81>I@!V8`#7(gR*uo93_%C4i=B_ejLH5?#bLqvGmm|DOlv{C_@o-?pm< zhOPEXns)*4_o}lKlb7Lk&XTRUsKS%k`)|~4k@(JDK=Dbbi?Ip<3%_nWIyv?|lhf$j z%*9f*8&@6*@mSyPel_|lN^8UT_H5=TKa(lvcLk5VX=0hN=(ybQs^@K_Xf3*GOPujl z(BpTa+83V(YOVBjU_DQib20XXszM)12K{nFMN|$l(L@$$PU>>Tv{okW6nQRizYfBt z;n>Gbbp{^?Rh(jRbpj$!4tIq^r+x(a@x@*WGF}xpmgw@{w7@~U*^w?t6`%d&Q+Ttr zndaWbV`*N(7sLc$Zq2%1O+{n}9Ter-E`3C-Z0OGL?B;LhMDYbz67A!U6NVmQa>c?c zZ{vHNlKkiI;qzoD7Y+1nQT5hIn+@tNeQCDqbg>$v*NJIe=C=W9nXRMRcjC{?V{)K0 zc2B;a7*W4BlD24%(Nx^ba2g;`id-+LpWfwS7fmR4%1nAga~_NT<-&MoIQJH{9PCuN zgjMpL5t5{OVz#X^X1+jip57Yym@uU$8iSl|$36Pyyz}Db0Uyf4v-n&PbSYvdr`E`2d6R}kHDh>xrAY_-$zbD&V3gtV@d_R4l&Js=l>kpOMwQY z9c?N3+9@UTpwJv!Y*$?7fLtKV5BY!#rI6*p{P9^}DOQZQ3d z^B#}D*4livTS$@+El_TLCddo2AHF$W^ufUhDq7LEQ|{DoTjN#$f?N)aV%9&gTNjTm zww5Gi6p0cKtyi->bA*>wZsuwb1Ji=;_=T$ASRh8Ve=B?;pQnfpF}y%;tzGQA_QyEP zDJDDLgDfiime&xHJ`_)*W*CTQSKmSf*9vVN)l(h-(*uUw$o0WVHIJw1tI+N-f|sIH zmJ=|Lgv1Y>Vjr1hSqMXUj3hdQs!rz_2vxnJddWOybz#O0vfs6N@^93kFp)4vUL9{n z&N=58ZLC1jCI%<)vtroumyMUOT4Km;ecIt?$jCf2m4gSYJ+xb|YJT@b(n8UP&xO+h zg1%m?MJOJp|12~CpwK+eVHhM^z>IqG>EGQ_c%{mG)9Z*ev2ItcrdF~lPwykAw%_X! zdoTPk{#r0*&wpc9rzYu87xt2_>pVNA7K(dVB~bhbktAPBA8Rl3*f6=q!rknTte_uOE zbz}tm694Ae2DUgIni$)0W17!Q!ewaatks~cAKeAEGs6o<9iKUmM5+lb|35gVZi!J-vM3`NiAUvS54AupSsrU@qX{TO=*D}|ywKl@Cv z@!c!#;yq@e(ueeMI^)bawkeK(3Cj8IkK=vSpj#_DrlMiuDY{h@EiY14gI zmp0H8#Y4|tI}ss%`vUfftS6wCtxg}0!U+xYr~2whFA$Mv?V0aFm51y|f2imAsM4@s zZ+&O&^qowQvrgyyzlHOvfA+b{&8NSP<~6zb(p)rW(H_8G8q_WBdbFn5vm{74*U*GAF%hM{h189aZg96m|Azc$p#Zo{trTttE{G6-2@+0gJn@!ggWDh_i80Cf;@3}Zry(GRUKVP3Re+|j;4L& z2^?1ams_JKfKGlyw7mn_%-6t@*x#sL3a|NDPY2>w_~AvsDD5sZ2P~E|pBHtr^`?cZ z%7_XSmghW^WFC7sUA~LlC)3z`t(jsqwdkGd!+6Z&U;)COH$$@PJs>}K!xchXqnRR4 z@tI|AMA0qrf}`m_h?)L0C>pH*0I`qG9oUO`2NT1>bedYU!tlw6)#}!0IfDUpku69I zo2f3TGQ21!-+2b?ej6-Sv{(GTtlmMi`pM1DjCZvdQ67`J>;jID*9`w^ynhi3zm1PDLb+p_d)we*@Xxw%}_ z>i(te#9c;_*x)=ob0KBuEFx)l`i9Z>L*sBGJ-rug>#B?eofu}-tOxg;<5QExn265b zDmrQ=E}IQ}n)lfk2s*_Cy>O7R)VjXgE&XSxk$$PAFV6;#BPY_KIR0X9($-W_o=E#Cg8&y)T42Tr!qq39m7tobK~<_?*!J_PvP66VUt= zz0m0FR~P86&b>ctw&`wLy^n~sw{05gC_C)}p!oC^a?bFfri1UWnSB39S$nVhF9rnJ zWYkOh-PCm47NW?OVr5QHg71y^et{B0VH@xLl^2))LQmcKdF$tvOY_n#DVV<$uz`I} zzGmqEZbDroD2^C6)8pFxTRIS(1okz*d) z?pqA%FpR*ds@~ zq8b=PGM3n+dkBAkmyc{F+t(~lrUcsC7YI-{83iKCZ!ewU;1`{tYLaEDvWu^)aCL}v zzewBkFfeRMgy9vDI)X?Qk*;~BTsq2cmZKKEe+rNtV3y@K=?V8+t@JT`VNBKZ>S3=s zr%?G;T4d{6*;+BS@yTQnh9fJ~5hK)G%1Mn`()=h-=JBEZqg={$UEUW=Ang+iB31O? z9Cv(i)T*9KU5w*ym&IL=(`j~|H$A~1VD9+7@9AGVg0eIE3|$lA^ucGrw1;>PjXzo8N>giO}}swSgnd1u$=m_u(nEUtkOQCAIJg`dzhOvCSuI8#I zWZXz;{98!wYRB@?b08xaT=UMkYNTHi2@9hGeFGPS|CwJFLdV`oxW3{0zI` zm2qYK-PhY25>yK-6LZS(eYCP6tZi;6zAxQDaE7h#e>)94>OSLP6|)DFrB8lua(=E& z=lUl7)w7yDK3T}t%P79MK)5h<*9EJp%HShhCJeG6nDAxaghl0LS}p>R3clD4mH6^; zl{xhYmm}`eqZL$#&&-$+iwNfy3#m+e{(_WS*5C<$Ob}dHY@_Sj zHNkdof$v1BrXO$(|AA{;%Ir4Zqt$A~;>R70culCXb_fy}=E+DRXm%|3-|WjnY_z`frINlg38AqH7~_$D=2EGoC$$5&kG z*Qg>k!E;p2`^#VOg6e`*XZ82HXtkPkaD)q}@;SHkcgOrkbe6hLrMnh$OBiJhGO(da z^{oC@oz}%)SouT$_@Z|On zu;o@%lssSb;e>0O`sdBKq6TSu!m_6nLj?NYlLFZl2UlE=^Ofi2YkiSUiL(iPcM|(sU4i_6{<2#N6mMBiDjWbmgxdkFosg$)L1( z#@{Uwk`c9V)bEOjU5{^N%L}CCnkF>o5%qvJ8m4_@!T%X)Y93j}Gee4&{K)BpY zawPCn(dvOhAiTjJHB66^DxcV9z0lhMPVSN03WGaoV3b7dTC{VmSaKG@vI_Y)p69nn zlBqh%%Db%(Ed`(<(h~U6#e)46merM$ec-+B@Zz3?tXK`|Iu2f-nN6%90{ycNU~Y8m zwAK0wuI!o907#x4Fd;GzBcJgzS<2H91xgy~%%0R` z;Mu`mBI^9<%5$&ktaHsqMFOr=^MI5Dq;U(@MJ!-PLEbF@H z21(1I3te>my}xc!MWClD?|5A#bJwDOt(KJ$ce`TWZDy9yHl0ibGa;Pdps=yK;XoBK zH?cz4!Ib23d>}brXi0`s&?6gz~J(4mc)50on|at)pu%o$Un=TRP-s zTlo0%qDoODRns(e1TTv+ZtfmHn3l0b6OS5n)y2<2uz9dweKi<3ciU|`~K-Tc4<)S@kG z(Ib(Rqck89CB z5PZ7WeVT}oL#Q|Z-k!jD17Zb!f(`O|F6FuYt@7tpUQ~^P8(hL$z3yM4l2e|8!+V2k z(Oi`c_AI;-c-eFsRC_!4QoAo04tkmv#FPA|(N}tLJvSD{3nX%6-s`Vi6I}X(k(K5Y zFtP3|u5Hdme=pUb?!X!b!H4*NB2)yPs?X#NZ-FR)jK2ft2u_J1rcml8kC*qL z6bcoi>~!EG4AI{~4WCy=1{$=Y4u{@W;}>07+l54w#&4?2=3J`buwy8?4Qk`X+X%{c z^L?Z+O%shp6NU)CQNV7M>iZ|q1&~cjE1tkRZ4~5U0V4L77Api7KL$mdkNrrEHs}By4TP@{wb)@A{ zOIrjQYOdUXjNj`E0^;ob2YDzPPBY9zVJw~NBI#%y!7JbvS)Y{tFS^iZ>d{2Z{+Y`j zU@#ZwzbhSqRJ@Mj1G_o>NCn`Ix^@c!PCgKl_!ECF(iTFk)R)KKp)_t8zx>6r_&gje z8CCt`ZH5_{s;Ikv1{md1yhs?u0>oLajYvd&&!Ygt7ygPC)r-@t~R-v!bNbh3khu7hEnd{r*Gj zhPZ^}|2&hB{oiMjVlpu>BTHQ|@&9_exVZRrng4kvas9ttiHl1}{rAU;i%VXYh&c^a PVU!YQ;^Vugt-|zwp0U11^?u*q|9$^1*X3NtdCv1Z_vgOv&;6WwJj*DaC0dnTQBht2_EJ;e z&jrr35yAqCRb_zyx`mL^hl4ei4-+~*_p}zC@^z4uzbA1q=1}z0`_}UDe8&P{d0f#o zX*svwhjlX|4n=%_mGyM!gRv~W_vCc2+F+?#z?M|hu;%vWVquh~X5)(5I=aM;5oNLN zM%79FsWLSd>LR|gy_RL@!(Rem zXWiM_P}l%CAw;!trj$OWq%DmdR zLdhF?w1XtTGi?>ClmEG0AuWK>hRhy^6t-$sxiiIIx-GI9(k+!luO+;-5}ozF;hkQ= zc_-hc8rWLp34!d$m?TX5rFv!7FsCd7)>g5@0I<1Tv$IUni(*s<6pbx(rv-TVt*?d_ zUU5~cUf-ta?HE*V&yFQqItT!S2-pTe-0x;#lm_(_?g_&$i?W1jg5oEp>9IGCV$Azr z8tnp*rP)((0Yv7C_^mtV&fhJQ7DZn@x*f6vrOrkFgm=>i4biDfZQ=Z=k6AW-ZGdxK z&8H!ZcZsaohQL0{X%yZS0aUrhs?^|JEj2C5$UcdvwHL<j#97_f9Fg>xoW_g*TH0 zVJM~IXC7~?_II`L^d+orhbH7s#z-d!PDop;+_-|XBeD8Vg-N5;Y-n6(o8?kIEmY~( znDv`Iz^}`_o4(+?b?fd^XIR~m7hu@sVjOg(Vm2K?bH_QtU1HKy4SxvfOqV$Le_fIVKI{$vF=7+^$k{H z6;UG&%C>&@gS2}?eC}c*bl&N34>gs?Js9<&Sd(Xh(pToVS&UXK(p<@x12AgxL$+i2 zr|DWnjbk6)>ym`p*jK-}o1K@Ebb|Tdn)1oE=w@qycXFo79hIt60_I!xk;Kl=cU~*D zILipyzO$M@vpIUJ_X9JS#-W66n zn6+6kuTh3-l=K|L@r=(EKHq%kv%cY;)RLHpYkwd}f@$kU zY!pYl((rWD6H@qT+07 zS!Vut`Kq&A$4(2y+kC%hi;4Y2@c0-7wVvi?3JA|0T(5LXU9!zz#&ub@CaE#-CB?Z` z#81NZyt&49MGOim^EzDFKftl=H*+&i?L6)BWpHJ~JLs z?=AcSL<)TldLc9p7B^;V5(zgPkZrf$a}A}oRW`wjnI_8|otMlhMzQzi;H!R9_;ZX_ z=BoYxk+6SEZPn*vM}Y27ImZ2VEH)l+vRQ9O2w2LQ75e<6(CK5*K=0|Pd{KIukvH`4 z-Mx})MII?l2=bCLqjl?4w;^rAAiU>g6=BAxCC%4PUM61X8-mBl&u4t>JobWbknMQM z>8TV$yzN7ZQ}p+7c`I|nO@-jMXqqi;>7K)Nrx z(}z&hF%M_J2vVNDGqyd!Z~hG-ULWE>DbmGKVaM>K$ZZvdY?OY898c5sDMnsY#uIGg z8Ce>Fj6tE2@vU4Klm=7!m-1#z=|tcV>IH}A`$Fw&MXkb}AX@S+It;&d zS{au5MY$nIsY;SaFsVtcB7E@pIhY;l(T4~>^m*wKBX8^)P>I5l?3)m2!IjUc7hdGd zdYbU0J_yC~PaEN}{F3(6h#}PNwxqmQ`)O5j`LFj2Z&em|n%SRZ45%q%9!5X3bvUy~ zIOX7bF2R+Ncher?aCT4v>&B2BQ{I|LQCi18ycDeae1)J$9%YF0^V}+q5#3$g`6v?}SoAeaGj@rHBM=0cor= z9*CStQR>NQLI{Z6VfuwDi}xWCMu&Q-;eZ!8HPY(Nj^vf#Vjp+XidTEUy3eZ>D;*WH zHX;9&qM_~=ML$w#ljZo<&|Xq}WM+~EU1nBU-G2IkX7lC=Ay4kKJVBBlO&_ z@V+NBHt!ay`0P`)Y=an(-vw_F%Kgf}C{~)Xrs-IXvYJ5z)(2UN#Rqt~8rPcwAiI$p ziJGPF7dN%94fXvhrJ+%pMkPJuEw)*v8U5S9kb`1pj3&IidJqw=DJ8=^yf8->pK$5Q z&;JosA^xTVgszFR4l|5g#Rlu3=MLKRtD_XaUORDF?v+132y$!Yr zLnyswOcF2ZPW<$vlJzG0Ac8dYIjI*$8?=kcz4(IPGt}3PX(#$X>8=mkFVJudDJdT9 zb_`iucXKXEEL`+RCqj?zIhV5aqYX>08~;AUGznOUz$ys3SNTn`G>WGdvF9YV9b);xOQ{c!JX0XPfRHV~qY&lSiRhbaK%s$allNMr-+&K14SFTS z!#Y87xrM7FaqZPywntLSyM!_^K(%AA=#*>m=>%i1lriY-Rb(qy;f*X&ozACnVgi*& z4KLxJP@4IdP&pQp@2|JSpEb$_=5D>3NR-C;Y;PP^rJr_h89oaBS^J_l|GY+F?ThgE zBK6765$kR6SCjj^=G)qAiHc@GDb-u;OUkw3q#VnN5iH9iQ|}9z&-7O}d4?zF)CiQ5G?M@L!rj#eW*GD%l+kU=Hy0Vb zBX|jm+bnZLYv#&h^v%gOC*a_|*{d$2g(IUP}|ns4aLzB@K)jqhNm z5Ps1xGEPw2r{< z4vyCrhRQ@<@N+8T)i4p_;Q$pX?vPQ@szx-v4%;~6+Ix3#+Iy+vHBkDDwnqhq3O`oYE^lA}O1 zO!lfdfvrTY;Q438Gf6jl)R+~+uUDDFk#z))`53_D9LB%> zR?=Jg=qS5|Qnx4WuX#g1D0oG zayGlaA<~)c-eEDV_G1{MlqfB7G~}q5M85>PDZRE#Ad*+lVAeDz#&i5q`i6GfoEi)T8y;00*%YIr2_rCSLcB-uB;oC85_!L@8@QtA;kmD+D#u%KY_~Zm_hns z%^1^R%3nE%qLPFJ8_@=jMjyDJ%&y>VSY;Z78B7#PzT97ZmNNm;jmCOkQ!u*Jm*b$( zUx2Hu0(oWR7S{pdb*uqPYVuy-8(^=UJ&Lf$M8P?x+J`@h^f2Rr?TVTC!LiCdB=<}X znH3}Oo6O0ZqlmqDl#4jINDhPK|VjPI`+XkvH`ir`p)SO4u)-CNs<>5(r+O4=G5tE zi2s2a;JT1}oOUEr_|Lq4Q&;-(lu3A5&CdF(5P&vB+lg%QuSi?hJlpKizh455@( zwcA|FaLRRp@WoGVzsn1;QThv%GSLh+J&P{CM>}Hy_11;>2;&3Y4Jtp#aYpI{KDd}s zlUrAyDlT^+G|zy#E2P>Dr~@#P9LO33Dqopd@DbE+2~N?s#%O%12zQw)vN$; zP()u{U)+H9NNe15XafIUM{&HN|GkdVSTA0eXz-XNvax?Ad2`Bu3?G`9MdMt?&Dg|y znjr6~9kDJEcjX59)_60(LP;p^*H~H5$|15J>$9%N8659C~2Tl3M;DMRaPfxHtN1%q5y(GNk?hM2a>m zc*0@^t|B)CC^Ngl!AeKI%gCe*C|Vwa1RlifdB~q9)CaI!iP-um5RxRN0)Ss%r2Z&b zk=VVz52K^*_@5mn-eU9it#IFAWQYlZDAEPLWGQnj%1fH%$U33)N^n$_y!CFyTKNXx zA+mEDSx)+f$$wn{SjFwh&nv9Tqu@U6Rf!tn{usgd^2l<#V!H0G!u+X}rBHvBF#OHk z(J@+js%D0k)#mLb1~C3<8ncUR?E>}*HW$(9@AjpBn%7++ohNNO;F+USjfL-CoS!{4 z`p@`Su*h-bnan3rELE!3-z&3bpYdh>kK=j+ZtVC;r381P^tWW;)^ZsNZM5M7JgRd; zxV10Eg;_=Z1KQeG?XRq3|2PTM=!gGpP#eawKNSNz9nQ2jALUpmQZnTYyC4Ys6hfR; zIR%_Z!=i3cVHn@%-iFQnHf~ktIGX!Qs#ElIrFSc3yF_RS-)=y<4gv*fda5%SN4lNY zcx0+1A9OzQ95>Zq;5*yy4L@VPSVqw-HhhU`-o8w^bRAl{qwI_@ zBtyWO;X^ZCzfZoX`OiHjGHwshbI?>Ra}by~xC4(e$WAdwBZi&J%;WDY#nwsH$(%TD z+69_-1DZII{MQt2$+~FFvvz@iOKN+UG(4JPoxfQXd24{`j=UATTPh)&-Bu&(ryYs+ zN%7OBX``HWSrYe~-iyx6c{ADT^N(x_s}dT+*A~uPUop!w0)`x7+Ym0zF{)|;OMF{- zO^5JF=EZ*-PwED(7N#=WIf+2dI~=v^uMvGeTh zIEm?fsEyvV(m{2B|Dni?GwL&(Ze`F>Wa#24H@FZK=wTR19?UxG;_qhH8jVcrNz>sK+M#0qxr`C zL)am_HfmSW$T6*i?vsvJ{&UCT;VCRPs=~lBXQy6Yxj)W%SdOGZ8IW35p4qHpdcwd& za8u{OWq?bAd>C)8rQAF!%OzVPw}je zdNaDTtU|qlzC}XUyv7eaZ736YzF$O(<=i;5XXSAK8_8!6iNX&ibKnuGhgSqK%tE?J zlR~ppbqm8De|Cqt<}G|Sslsi#%rBgH_%)`2g`{tFZcl=!kUeCix4%|?`A-y~)fuwq z%CM30-m+i{;tI&N*fwFfr{hRs0cLy*@-13v>VSw=Q-{DqGwOYF7Llyj<-Zm%R4cFk zx`vIU4_Fzku%~TsPoVaYnmHn09oq~Oz?)M)?iE3E%@-Y~bw@%*66JlI`jAMPxyqSPY*I{qn^rs6n5PYDf@f?sU6G}= z@@J264y{MCe{}F1@m1l1QvsOk4KPfLa%GWYD%eP2ETkh+S|}bGpX-ad(qFTFnCpFc zrNQ~K_bIjVG7sRE`q$mW`VK0vkgf;jpZxKyipG3D7y5gj=;xP~rE>R8RR^KTavoJK zs-*yuOC;(y*dRo-g@tEuWZ$*b43@tC@fNXlJS7(WK*068}O-gx0@cKX05Lt_W z>$SM@AhaFrw^Nj|k)qWH!ME1^I<)dctIa@QGdJ*e%`-#hG~89=!u1v_Pq?y*~#0WDBAp9&C}k1R0MI zN0N*mh+g@xKJdh_o8E=t*-`UgvAHZx)RH2S{l_PJM+!RSeCD0#<;<&sTJd*SBXM8S z8=($+R9$@8W_Am@I}QK`mdt_!pC9{`RSTUm-BBJf9zJ(DuSS534q`>PH>uCvvwV(63n|RpnO4%YOJqYAc z^Mbu1)57qQe){k8pt5UT6v%`)*&W%CexSjb*})MR!7FmRa>}n#!V`6i3*H!WTTWPu zfB1V?Nn0R`*WYubva2JJE%Iz|qLh8+Xz-1szOL6n=&FEkOIqM!T|SjVky9KX+wULG zZb}WAAN*aCXKR164QFJEXHPO%lw(5Tn8XjBxDjCkF1Vox057>`l;ZNf3W~(Ak*fBR zUK}@_^^*-mYm5*5i6#tCo+MMS)H5VkAy&jYDo5n8r`tjrcliC-#6a&AcFpTOSQCf0 z(Tf=U&tw5gMPrDSr2iOs&o8e1ToS>>RX}RF=OvY6IFJNIo?#(nUJgnGptNsKO9xr+ zjr+%nCtZR~|I4}=#naEXI_|qEyS6wj{?5w3qjGE@df?(xFCOlBvlh3@R%K=b_N&@% z@oSr}`;5&I6p!`0*aNb@HGCVV17qu$p~Y7cbnfMW(34CfUOM1VU4O!jB4}4?BD_8vX0f+sC=9M#-ImBymM}5=$5r5m+ zqwk5wL56qiNobdw2R@n~!GNCjopW<@-EV#YHyo{W2K?`1%l3zTVRCAmu8tw>cHki- z8$6REdUL$L@YTx4k5t~co=+dC|II7dI6dvaO|>>%P;)OZwlw@&I|0<2gDxX^lWstE zblG*S?E`QsFPwV1)>J6h;_$gpUXfgIm$(ApPpSp&?plxNKX@2Q7~Z&vScg?l+u)Un ze<{Zn+5VlJbQfHw@q&Qtyz1J3-;apnyL591xhbXdHI*#!uv4?CR=@zXo!?>Db3|~F zu&WO+3eWh5Sbgp06%jrlGR-d$`#X*98o-T12@OU-Mds|Gk<*NS<*&cMA%bqIt#I^o zO;=u2Vj~Hj#azu`l&Xkpaqn7A@3|4dnp^Tt5u%{}S@(S3ja5lNock_uwnSa?mLREz zNAfY{y{`k&^_EBK?S{wc*m{%Oy9VGvvXZf;*esA){yn_|Agmw*d^kl+J=x1<+*mPV zz69(+=uXFi>xC{B@r%`B%!p`@-^yb7^p-5f@j9HLDIyFH*>!?bzMWv%*DfsZf|`AI zC|>O5_g0kOJ9QRP4jai;@~EHLW$?(b+e-r2PFizH_z1&cqW6?H^DXB8x)zY_y(9|v zTx=<2t9qw`y+Yvn-M{EFlv2MTn`>Z68)_vnug&3iF>goLo_VtC*>)v{-&;EPnPS&K#2@mYf`xr(-C3o7lOw2<$%;y_*FL^2 z&F;$qc(O=dx_)QSnATymo z=3|DPJ6CI7PfUAULm$|jN^hTT1|d15ea8HWW4AH1zu?fc_plTfQT9XCZtL3Rf@xqo z*!3!GXQ#74uLDCxOKYjHbHeYIR~O!VJlrw}xLzwir1m@WfuHLB<-7Fe6Bn5DO$XZ% zL6_qd{<+Y#__oSC11xBR@gX))-+Mc!MtX~b(C}u+@6TP}gyUX3VEch$e18&)wp16l zLvZ!(XBfpzoS2i=yLH^ww&6t)kAbfBQ?Al?nDXbe4BV5cGzEKm{Ec8E)S^RL!qr+F z*ul%N62lj|emlA_he+q}Rgb=c4P_P*l+zUVZ=j%pscGMjp1Ykh7A8RWcaB|28`{oS z4tn4I6JP?ElGY2f`dtF3fWz1-4qi2*f()R74B*=(%gK?%gd11VC-vxmNF{GzKw{48cazqW#v}q3ug?|(8XT;Fn$X9KE%MH!y&?8bB)g(A9AiDo*K;KbK z`0ViV26Fg$u^v->{v+QO_hgjQFH6g9fitixZiKSy*LgLRVG`$deBtDJ`^ z{hxb^^8cTEityMA##Xl!75>-b=g*&4fd8L+3JU*ws@=!__iN9eS5#1pT{hO{f-7)C KAXoIXx&I3(>xUr#