more examples

.github/workflows/R-CMD-check.yaml

@@ -52,6 +52,7 @@ jobs:
               knitr
               collapse
               rmarkdown
+              ggplot2
               rjags
               coda
               runjags

.github/workflows/test-coverage.yaml

@@ -34,6 +34,7 @@ jobs:
               knitr
               collapse
               rmarkdown
+              ggplot2
               rjags
               coda
               runjags

DESCRIPTION

@@ -46,6 +46,7 @@ Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
 Suggests: 
     cmdstanr (>= 0.4.0),
+    ggplot2 (>= 2.0.0),
     knitr, 
     collapse,
     rmarkdown, 

NAMESPACE

@@ -128,6 +128,7 @@ importFrom(stats,logLik)
 importFrom(stats,make.link)
 importFrom(stats,median)
 importFrom(stats,model.frame)
+importFrom(stats,na.fail)
 importFrom(stats,na.pass)
 importFrom(stats,pacf)
 importFrom(stats,pbeta)
@@ -165,4 +166,5 @@ importFrom(stats,update.formula)
 importFrom(stats,var)
 importFrom(utils,getFromNamespace)
 importFrom(utils,head)
+importFrom(utils,lsf.str)
 importFrom(utils,tail)

R/as.data.frame.mvgam.R

@@ -38,9 +38,12 @@
 #'head(beta_draws_df)
 #'str(beta_draws_df)
 #'
-#'beta_draws_mat <- as.data.frame(mod1, variable = 'betas')
+#'beta_draws_mat <- as.matrix(mod1, variable = 'betas')
 #'head(beta_draws_mat)
-#'str(beta_draws_mat)}
+#'str(beta_draws_mat)
+#'
+#'shape_pars <- as.matrix(mod1, variable = 'shape', regex = TRUE)
+#'head(shape_pars)}
 NULL
 
 #'@rdname mvgam_draws
@@ -228,7 +231,6 @@ as.data.frame.mvgam = function(x,
   }
 
   return(post)
-
 }
 
 #'@rdname mvgam_draws
@@ -244,6 +246,5 @@ as.matrix.mvgam = function(x,
   colnames(post) <- varnames
 
   return(post)
-
 }
 

R/dynamic.R

@@ -27,6 +27,61 @@
 #' set automatically to ensure enough basis functions are used to approximate the expected
 #' wiggliness of the underlying dynamic function (\code{k} will increase as \code{rho} decreases)
 #' @rdname dynamic
+#'@examples
+#'\dontrun{
+#'# Simulate a time-varying coefficient \
+#'#(as a Gaussian Process with length scale = 10)
+#'set.seed(1111)
+#'N <- 200
+#'beta <- mvgam:::sim_gp(rnorm(1),
+#'                       alpha_gp = 0.75,
+#'                       rho_gp = 10,
+#'                       h = N) + 0.5
+#'plot(beta, type = 'l', lwd = 3,
+#'     bty = 'l', xlab = 'Time',
+#'     ylab = 'Coefficient',
+#'     col = 'darkred')
+#'
+#'# Simulate the predictor as a standard normal
+#'predictor <- rnorm(N, sd = 1)
+#'
+#'# Simulate a Gaussian outcome variable
+#'out <- rnorm(N, mean = 4 + beta * predictor,
+#'             sd = 0.25)
+#'time <- seq_along(predictor)
+#'plot(out,  type = 'l', lwd = 3,
+#'     bty = 'l', xlab = 'Time', ylab = 'Outcome',
+#'     col = 'darkred')
+#'
+#'# Gather into a data.frame and fit a dynamic coefficient mmodel
+#'data <- data.frame(out, predictor, time)
+#'
+#'# Split into training and testing
+#'data_train <- data[1:190,]
+#'data_test <- data[191:200,]
+#'
+#'# Fit a model using the dynamic function
+#'mod <- mvgam(out ~
+#'             # mis-specify the length scale slightly as this
+#'             # won't be known in practice
+#'             dynamic(predictor, rho = 8, stationary = TRUE),
+#'            family = gaussian(),
+#'            data = data_train)
+#'
+#'# Inspect the summary
+#'summary(mod)
+#'
+#'# Plot the time-varying coefficient estimates
+#'plot(mod, type = 'smooths')
+#'
+#'# Extrapolate the coefficient forward in time
+#'plot_mvgam_smooth(mod, smooth = 1, newdata = data)
+#'abline(v = 190, lty = 'dashed', lwd = 2)
+#'
+#'# Overlay the true simulated time-varying coefficient
+#'lines(beta, lwd = 2.5, col = 'white')
+#'lines(beta, lwd = 2)
+#'}
 #' @author Nicholas J Clark
 #' @export
 dynamic = function(variable, rho = 5, stationary = TRUE){

R/evaluate_mvgams.R

@@ -4,6 +4,7 @@
 #'@importFrom stats quantile ecdf median predict
 #'@importFrom parallel clusterExport stopCluster setDefaultCluster clusterEvalQ
 #'@importFrom grDevices devAskNewPage
+#'@importFrom utils lsf.str
 #'@param object \code{list} object returned from \code{mvgam}
 #'@param n_samples \code{integer} specifying the number of samples to generate from the model's
 #'posterior distribution
@@ -56,7 +57,75 @@
 #'extrapolation of smooths and let the latent trends in the \code{mvgam} model capture any
 #'temporal dependencies in the data. These trends are time series models and so will provide much more
 #'stable forecasts
+#'@examples
+#'\dontrun{
+#'# Simulate from a Poisson-AR2 model with a seasonal smooth
+#'set.seed(100)
+#'dat <- sim_mvgam(T = 75,
+#'                 n_series = 1,
+#'                 prop_trend = 0.75,
+#'                  trend_model = 'AR2',
+#'                  family = poisson())
 #'
+#'# Plot the time series
+#'plot_mvgam_series(data = dat$data_train,
+#'                  newdata = dat$data_test,
+#'                  series = 1)
+#'
+#'# Fit an appropriate model
+#'mod_ar2 <- mvgam(y ~ s(season, bs = 'cc'),
+#'                trend_model = 'AR2',
+#'                family = poisson(),
+#'                data = dat$data_train,
+#'                newdata = dat$data_test)
+#'
+#'# Fit a less appropriate model
+#'mod_rw <- mvgam(y ~ s(season, bs = 'cc'),
+#'               trend_model = 'RW',
+#'               family = poisson(),
+#'               data = dat$data_train,
+#'               newdata = dat$data_test)
+#'
+#'# Compare Discrete Ranked Probability Scores for the testing period
+#'fc_ar2 <- forecast(mod_ar2)
+#'fc_rw <- forecast(mod_rw)
+#'score_ar2 <- score(fc_ar2, score = 'drps')
+#'score_rw <- score(fc_rw, score = 'drps')
+#'sum(score_ar2$series_1$score)
+#'sum(score_rw$series_1$score)
+#'
+#'# Use rolling evaluation for approximate comparisons of 3-step ahead
+#'# forecasts across the training period
+#'compare_mvgams(mod_ar2,
+#'               mod_rw,
+#'               fc_horizon = 3,
+#'               n_samples = 1000,
+#'               n_evaluations = 5)
+#'
+#'# Now use approximate leave-future-out CV to compare
+#'# rolling forecasts; start at time point 40 to reduce
+#'# computational time and to ensure enough data is available
+#'# for estimating model parameters
+#'lfo_ar2 <- lfo_cv(mod_ar2,
+#'                  min_t = 40,
+#'                  fc_horizon = 3)
+#'lfo_rw <- lfo_cv(mod_rw,
+#'                 min_t = 40,
+#'                 fc_horizon = 3)
+#'
+#'# Plot Pareto-K values and ELPD estimates
+#'plot(lfo_ar2)
+#'plot(lfo_rw)
+#'
+#'# Proportion of timepoints in which AR2 model gives
+#'# better forecasts
+#'length(which((lfo_ar2$elpds - lfo_rw$elpds) > 0)) /
+#'       length(lfo_ar2$elpds)
+#'
+#'# A higher total ELPD is preferred
+#'lfo_ar2$sum_ELPD
+#'lfo_rw$sum_ELPD
+#'}
 #' @name evaluate_mvgams
 NULL
 

R/lfo_cv.mvgam.R

@@ -37,6 +37,67 @@
 #'the Pareto-k shape values and 'the specified `pareto_k_threshold`
 #'@references Paul-Christian Bürkner, Jonah Gabry & Aki Vehtari (2020). Approximate leave-future-out cross-validation for Bayesian time series models
 #'Journal of Statistical Computation and Simulation. 90:14, 2499-2523.
+#'@examples
+#'\dontrun{
+#'# Simulate from a Poisson-AR2 model with a seasonal smooth
+#'set.seed(100)
+#'dat <- sim_mvgam(T = 75,
+#'                 n_series = 1,
+#'                 prop_trend = 0.75,
+#'                  trend_model = 'AR2',
+#'                  family = poisson())
+#'
+#'# Plot the time series
+#'plot_mvgam_series(data = dat$data_train,
+#'                  newdata = dat$data_test,
+#'                  series = 1)
+#'
+#'# Fit an appropriate model
+#'mod_ar2 <- mvgam(y ~ s(season, bs = 'cc'),
+#'                trend_model = 'AR2',
+#'                family = poisson(),
+#'                data = dat$data_train,
+#'                newdata = dat$data_test)
+#'
+#'# Fit a less appropriate model
+#'mod_rw <- mvgam(y ~ s(season, bs = 'cc'),
+#'               trend_model = 'RW',
+#'               family = poisson(),
+#'               data = dat$data_train,
+#'               newdata = dat$data_test)
+#'
+#'# Compare Discrete Ranked Probability Scores for the testing period
+#'fc_ar2 <- forecast(mod_ar2)
+#'fc_rw <- forecast(mod_rw)
+#'score_ar2 <- score(fc_ar2, score = 'drps')
+#'score_rw <- score(fc_rw, score = 'drps')
+#'sum(score_ar2$series_1$score)
+#'sum(score_rw$series_1$score)
+#'
+#'# Now use approximate leave-future-out CV to compare
+#'# rolling forecasts; start at time point 40 to reduce
+#'# computational time and to ensure enough data is available
+#'# for estimating model parameters
+#'lfo_ar2 <- lfo_cv(mod_ar2,
+#'                  min_t = 40,
+#'                  fc_horizon = 3)
+#'lfo_rw <- lfo_cv(mod_rw,
+#'                 min_t = 40,
+#'                 fc_horizon = 3)
+#'
+#'# Plot Pareto-K values and ELPD estimates
+#'plot(lfo_ar2)
+#'plot(lfo_rw)
+#'
+#'# Proportion of timepoints in which AR2 model gives
+#'# better forecasts
+#'length(which((lfo_ar2$elpds - lfo_rw$elpds) > 0)) /
+#'       length(lfo_ar2$elpds)
+#'
+#'# A higher total ELPD is preferred
+#'lfo_ar2$sum_ELPD
+#'lfo_rw$sum_ELPD
+#'}
 #'@author Nicholas J Clark
 #'@export
 lfo_cv <- function(object, ...){