nearly there with CAR(1), just need some final testing; add pp_check(…

…) and more detailed examples

NAMESPACE

@@ -44,6 +44,7 @@ S3method(plot,mvgam_lfo)
 S3method(posterior_epred,mvgam)
 S3method(posterior_linpred,mvgam)
 S3method(posterior_predict,mvgam)
+S3method(pp_check,mvgam)
 S3method(ppc,mvgam)
 S3method(predict,mvgam)
 S3method(print,mvgam)
@@ -95,6 +96,7 @@ export(plot_mvgam_series)
 export(plot_mvgam_smooth)
 export(plot_mvgam_trend)
 export(plot_mvgam_uncertainty)
+export(pp_check)
 export(ppc)
 export(roll_eval_mvgam)
 export(score)
@@ -109,6 +111,7 @@ importFrom(bayesplot,color_scheme_set)
 importFrom(bayesplot,log_posterior)
 importFrom(bayesplot,neff_ratio)
 importFrom(bayesplot,nuts_params)
+importFrom(bayesplot,pp_check)
 importFrom(brms,bernoulli)
 importFrom(brms,conditional_effects)
 importFrom(brms,dbeta_binomial)

R/RW.R

@@ -16,6 +16,66 @@
 #' @return An object of class \code{mvgam_trend}, which contains a list of
 #' arguments to be interpreted by the parsing functions in \code{mvgam}
 #' @rdname RW
+#'@examples
+#'\dontrun{
+#'# A short example to illustrate CAR(1) models
+#'# Function to simulate CAR1 data with seasonality
+#'sim_corcar1 = function(n = 120,
+#'                       phi = 0.5,
+#'                       sigma = 1,
+#'                       sigma_obs = 0.75){
+#'# Sample irregularly spaced time intervals
+#'time_dis <- c(0, runif(n - 1, -0.1, 1))
+#'time_dis[time_dis < 0] <- 0; time_dis <- time_dis * 5
+#'
+#'# Set up the latent dynamic process
+#'x <- vector(length = n); x[1] <- -0.3
+#'for(i in 2:n){
+#'  # zero-distances will cause problems in sampling, so mvgam uses a
+#'  # minimum threshold; this simulation function emulates that process
+#'  if(time_dis[i] == 0){
+#'    x[i] <- rnorm(1, mean = (phi ^ 1e-12) * x[i - 1], sd = sigma)
+#'   } else {
+#'     x[i] <- rnorm(1, mean = (phi ^ time_dis[i]) * x[i - 1], sd = sigma)
+#'   }
+#' }
+#'
+#'# Add 12-month seasonality
+#' cov1 <- sin(2 * pi * (1 : n) / 12); cov2 <- cos(2 * pi * (1 : n) / 12)
+#' beta1 <- runif(1, 0.3, 0.7); beta2 <- runif(1, 0.2, 0.5)
+#' seasonality <- beta1 * cov1 + beta2 * cov2
+#'
+#'# Take Gaussian observations with error and return
+#' data.frame(y = rnorm(n, mean = x + seasonality, sd = sigma),
+#'            season = rep(1:12, 20)[1:n],
+#'            time = cumsum(time_dis))
+#'}
+#'
+#'# Sample two time series
+#'dat <- rbind(dplyr::bind_cols(sim_corcar1(phi = 0.65,
+#'                                          sigma_obs = 0.55),
+#'                              data.frame(series = 'series1')),
+#'             dplyr::bind_cols(sim_corcar1(phi = 0.8,
+#'                              sigma_obs = 0.35),
+#'                              data.frame(series = 'series2'))) %>%
+#'       dplyr::mutate(series = as.factor(series))
+#'
+#'# mvgam with CAR(1) trends and series-level seasonal smooths
+#'mod <- mvgam(formula = y ~ s(season, bs = 'cc',
+#'                             k = 5, by = series),
+#'             trend_model = CAR(),
+#'             data = dat,
+#'             family = gaussian(),
+#'             run_model = TRUE)
+#'
+#'# View usual summaries and plots
+#'summary(mod)
+#'conditional_effects(mod, type = 'expected')
+#'plot(mod, type = 'trend', series = 1)
+#'plot(mod, type = 'trend', series = 2)
+#'plot(mod, type = 'residuals', series = 1)
+#'plot(mod, type = 'residuals', series = 2)
+#'}
 #' @export
 RW = function(ma = FALSE, cor = FALSE){
   out <- structure(list(trend_model = 'RW',

R/compute_edf.R

@@ -80,7 +80,8 @@ compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names,
     # Calculate variance using family's mean-variance relationship
     mu_variance <- predict(object,
                            process_error = FALSE,
-                           type = 'variance')[best_draw, ]
+                           type = 'variance',
+                           summary = FALSE)[best_draw, ]
     if(length(mu_variance) > 1){
       mu_variance <- mu_variance[1:length(eta)]
     }

R/families.R

@@ -1631,7 +1631,7 @@ dsresids_vec = function(object){
   # Need to know which series each observation belongs to so we can
   # pull out appropriate family-level parameters (overdispersions, shapes, etc...)
   all_dat <- data.frame(series = object$obs_data$series,
-                        time = object$obs_data$time,
+                        time = object$obs_data$index..time..index,
                         y = object$obs_data$y) %>%
     dplyr::arrange(time, series)