2 chains

nicholasjclark · Apr 30, 2024 · 226aadd · 226aadd
1 parent 47bde90
commit 226aadd
Show file tree

Hide file tree

Showing 298 changed files with 6,758 additions and 1,018 deletions.
diff --git a/R/RW.R b/R/RW.R
@@ -66,7 +66,9 @@
 #'             trend_model = CAR(),
 #'             data = dat,
 #'             family = gaussian(),
-#'             run_model = TRUE)
+#'             burnin = 300,
+#'             samples = 300,
+#'             chains = 2)
 #'
 #'# View usual summaries and plots
 #'summary(mod)

diff --git a/R/as.data.frame.mvgam.R b/R/as.data.frame.mvgam.R
@@ -37,7 +37,9 @@
 #'mod1 <- mvgam(y ~ s(season, bs = 'cc'),
 #'              trend_model = 'AR1',
 #'              data = sim$data_train,
-#'              family = Gamma())
+#'              family = Gamma(),
+#'              chains = 2,
+#'              samples = 300)
 #'beta_draws_df <- as.data.frame(mod1, variable = 'betas')
 #'head(beta_draws_df)
 #'str(beta_draws_df)

diff --git a/R/data_grids.R b/R/data_grids.R
@@ -0,0 +1,93 @@
+#' Get data objects into correct order in case it is not already
+#'@noRd
+sort_data = function(data, series_time = FALSE){
+  if(inherits(data, 'list')){
+    data_arranged <- data
+    if(series_time){
+      temp_dat = data.frame(time = data$index..time..index,
+                            series = data$series) %>%
+        dplyr::mutate(index = dplyr::row_number()) %>%
+        dplyr::arrange(series, time)
+    } else {
+      temp_dat = data.frame(time = data$index..time..index,
+                            series = data$series) %>%
+        dplyr::mutate(index = dplyr::row_number()) %>%
+        dplyr::arrange(time, series)
+    }
+
+    data_arranged <- lapply(data, function(x){
+      if(is.matrix(x)){
+        matrix(x[temp_dat$index,], ncol = NCOL(x))
+      } else {
+        x[temp_dat$index]
+      }
+    })
+    names(data_arranged) <- names(data)
+  } else {
+    if(series_time){
+      data_arranged <- data %>%
+        dplyr::arrange(series, index..time..index)
+    } else {
+      data_arranged <- data %>%
+        dplyr::arrange(index..time..index, series)
+    }
+  }
+
+  return(data_arranged)
+}
+
+#' Create prediction grids, mostly for simple plotting functions
+#'@noRd
+data_grid = function(..., newdata){
+
+  dots <- list(...)
+  vars <- names(dots)
+
+  # Validate that vars exist in supplied data
+  for(i in seq_along(vars)){
+    if(!exists(vars[i], newdata)){
+      stop(paste0('Variable ', vars[i], ' not found in newdata'),
+           call. = FALSE)
+    }
+  }
+
+  # Create sample dummy dataframe to get the prediction grid, ensuring
+  # factors are preserved
+  newdat_grid <- data.frame(do.call(cbind.data.frame,
+                                    lapply(vars, function(x){
+    newdata[[x]]
+  })))
+  colnames(newdat_grid) <- vars
+
+  # Use the supplied conditions for making the datagrid
+  newdat_grid <- marginaleffects::datagrid(..., newdata = newdat_grid)
+
+  # Now replicate the first observation for all other variables
+  if(inherits(newdata, 'list')){
+    newdat_full <- lapply(seq_along(newdata), function(x){
+      if(names(newdata)[x] %in% vars){
+        newdat_grid[[names(newdata)[x]]]
+      } else {
+        if(is.matrix(newdata[[x]])){
+          t(replicate(NROW(newdat_grid), newdata[[x]][1, ]))
+        } else {
+          if(is.factor(newdata[[x]])){
+            factor(rep(newdata[[x]][1], NROW(newdat_grid)),
+                   levels = levels(newdata[[x]]))
+          } else {
+            rep(newdata[[x]][1], NROW(newdat_grid))
+          }
+        }
+      }
+    })
+    names(newdat_full) <- names(newdata)
+  } else {
+    newdat_full <-
+      dplyr::bind_cols(newdat_grid,
+                       data.frame(newdata %>%
+                                    dplyr::select(!vars) %>%
+                                    dplyr::slice_head(n = 1)))
+  }
+
+  return(newdat_full)
+}
diff --git a/R/dynamic.R b/R/dynamic.R
@@ -84,11 +84,12 @@
 #'
 #'# 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)
+#'              # 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,
+#'             chains = 2)
 #'
 #'# Inspect the summary
 #'summary(mod)

diff --git a/R/evaluate_mvgams.R b/R/evaluate_mvgams.R
@@ -76,14 +76,16 @@
 #'                 trend_model = AR(p = 2),
 #'                 family = poisson(),
 #'                 data = dat$data_train,
-#'                 newdata = dat$data_test)
+#'                 newdata = dat$data_test,
+#'                 chains = 2)
 #'
 #'# 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)
+#'                newdata = dat$data_test,
+#'                chains = 2)
 #'
 #'# Compare Discrete Ranked Probability Scores for the testing period
 #'fc_ar2 <- forecast(mod_ar2)

diff --git a/R/families.R b/R/families.R
@@ -161,7 +161,8 @@ student_t = function(link = 'identity'){
 #'             priors = c(prior(std_normal(), class = b),
 #'                        prior(normal(1, 1.5), class = Intercept_trend)),
 #'             burnin = 300,
-#'             samples = 300)
+#'             samples = 300,
+#'             chains = 2)
 #'
 #' # The usual diagnostics
 #' summary(mod)

diff --git a/R/forecast.mvgam.R b/R/forecast.mvgam.R
@@ -34,7 +34,8 @@ forecast <- function(object, ...){
 #'             trend_model = AR(),
 #'             data = simdat$data_train,
 #'             burnin = 300,
-#'             samples = 300)
+#'             samples = 300,
+#'             chains = 2)
 #'
 #' # Hindcasts on response scale
 #' hc <- hindcast(mod)

diff --git a/R/hindcast.mvgam.R b/R/hindcast.mvgam.R
@@ -22,7 +22,8 @@ hindcast <- function(object, ...){
 #'             trend_model = AR(),
 #'             data = simdat$data_train,
 #'             burnin = 300,
-#'             samples = 300)
+#'             samples = 300,
+#'             chains = 2)
 #'
 #' # Hindcasts on response scale
 #' hc <- hindcast(mod)

diff --git a/R/index-mvgam.R b/R/index-mvgam.R
@@ -23,7 +23,8 @@ NULL
 #'              trend_model = AR(),
 #'              data = simdat$data_train,
 #'             burnin = 300,
-#'             samples = 300)
+#'             samples = 300,
+#'             chains = 2)
 #' variables(mod)
 #' }
 #' @export

diff --git a/R/lfo_cv.mvgam.R b/R/lfo_cv.mvgam.R
@@ -60,7 +60,8 @@
 #'                data = dat$data_train,
 #'                newdata = dat$data_test,
 #'                burnin = 300,
-#'                samples = 300)
+#'                samples = 300,
+#'                chains = 2)
 #'
 #'# Fit a less appropriate model
 #'mod_rw <- mvgam(y ~ s(season, bs = 'cc', k = 6),
@@ -69,7 +70,8 @@
 #'               data = dat$data_train,
 #'               newdata = dat$data_test,
 #'               burnin = 300,
-#'               samples = 300)
+#'               samples = 300,
+#'               chains = 2)
 #'
 #'# Compare Discrete Ranked Probability Scores for the testing period
 #'fc_ar2 <- forecast(mod_ar2)

diff --git a/R/logLik.mvgam.R b/R/logLik.mvgam.R
@@ -24,10 +24,11 @@
 #' # Simulate some data and fit a model
 #' simdat <- sim_mvgam(n_series = 1, trend_model = 'AR1')
 #' mod <- mvgam(y ~ s(season, bs = 'cc', k = 6),
-#'             trend_model = AR(),
-#'             data = simdat$data_train,
-#'             burnin = 300,
-#'             samples = 300)
+#'              trend_model = AR(),
+#'              data = simdat$data_train,
+#'              burnin = 300,
+#'              samples = 300,
+#'              chains = 2)
 #'
 #'# Extract logLikelihood values
 #'lls <- logLik(mod)

diff --git a/R/loo.mvgam.R b/R/loo.mvgam.R
@@ -24,7 +24,8 @@
 #'              simdat$data_test),
 #'              family = gaussian(),
 #'              burnin = 300,
-#'              samples = 300)
+#'              samples = 300,
+#'              chains = 2)
 #'plot(mod1, type = 'smooths')
 #'loo(mod1)
 #'
@@ -34,15 +35,17 @@
 #'               s(season, series, bs = 'fs',
 #'               xt = list(bs = 'cc'), k = 4),
 #'               burnin = 300,
-#'               samples = 300)
+#'               samples = 300,
+#'               chains = 2)
 #'plot(mod2, type = 'smooths')
 #'loo(mod2)
 #'
 #'# Now add AR1 dynamic errors to mod2
 #'mod3 <- update(mod2,
 #'               trend_model = 'AR1',
 #'               burnin = 300,
-#'               samples = 300)
+#'               samples = 300,
+#'               chains = 2)
 #'plot(mod3, type = 'smooths')
 #'plot(mod3, type = 'trend')
 #'loo(mod3)

diff --git a/R/lv_correlations.R b/R/lv_correlations.R
@@ -17,7 +17,8 @@
 #'             n_lv = 2,
 #'             data = simdat$data_train,
 #'             burnin = 300,
-#'             samples = 300)
+#'             samples = 300,
+#'             chains = 2)
 #'lvcors <- lv_correlations(mod)
 #'names(lvcors)
 #'lapply(lvcors, class)

diff --git a/R/mcmc_plot.mvgam.R b/R/mcmc_plot.mvgam.R
@@ -28,7 +28,8 @@
 #'              trend_model = AR(),
 #'              data = simdat$data_train,
 #'              burnin = 300,
-#'              samples = 300)
+#'              samples = 300,
+#'              chains = 2)
 #' mcmc_plot(mod)
 #' mcmc_plot(mod, type = 'neff_hist')
 #' mcmc_plot(mod, variable = 'betas', type = 'areas')

diff --git a/R/monotonic.R b/R/monotonic.R
@@ -73,7 +73,8 @@
 #'              data = mod_data,
 #'              family = gaussian(),
 #'              burnin = 300,
-#'              samples = 300)
+#'              samples = 300,
+#'              chains = 2)
 #'
 #' plot_predictions(mod2,
 #'                  by = 'x',
@@ -108,7 +109,8 @@
 #'              data = mod_data,
 #'              family = gaussian(),
 #'              burnin = 300,
-#'              samples = 300)
+#'              samples = 300,
+#'              chains = 2)
 #'
 #' # Visualise the different monotonic functions
 #' plot_predictions(mod, condition = c('x', 'fac', 'fac'),

diff --git a/R/mvgam.R b/R/mvgam.R
@@ -335,7 +335,8 @@
 #'               trend_model = RW(),
 #'               family = poisson(),
 #'               burnin = 300,
-#'               samples = 300)
+#'               samples = 300,
+#'               chains = 2)
 #'
 #' # Extract the model summary
 #' summary(mod1)
@@ -399,7 +400,8 @@
 #'               data = mod_data,
 #'               return_model_data = TRUE,
 #'               burnin = 300,
-#'               samples = 300)
+#'               samples = 300,
+#'               chains = 2)
 #'
 #' # The mapping matrix is now supplied as data to the model in the 'Z' element
 #' mod1$model_data$Z
@@ -443,7 +445,8 @@
 #'              data = data_train,
 #'              newdata = data_test,
 #'              burnin = 300,
-#'              samples = 300)
+#'              samples = 300,
+#'              chains = 2)
 #'
 #' # Inspect the model summary, forecast and time-varying coefficient distribution
 #' summary(mod)
@@ -475,7 +478,8 @@
 #'              data = dat$data_train,
 #'              trend_model = 'None',
 #'              burnin = 300,
-#'              samples = 300)
+#'              samples = 300,
+#'              chains = 2)
 #'
 #' # Inspect the model file to see the modification to the linear predictor
 #' # (eta)
@@ -569,7 +573,8 @@
 #'              family = poisson(),
 #'              data = mod_data,
 #'              burnin = 300,
-#'              samples = 300)
+#'              samples = 300,
+#'              chains = 2)
 #' summary(mod)
 #'
 #' # Plot the posterior hindcast
@@ -610,7 +615,8 @@
 #'              family = binomial(),
 #'              data = dat,
 #'              burnin = 300,
-#'              samples = 300)
+#'              samples = 300,
+#'              chains = 2)
 #' summary(mod)
 #' pp_check(mod, type = "bars_grouped",
 #'          group = "series", ndraws = 50)

diff --git a/R/mvgam_diagnostics.R b/R/mvgam_diagnostics.R
@@ -21,7 +21,8 @@
 #'             trend_model = AR(),
 #'             data = simdat$data_train,
 #'             burnin = 300,
-#'             samples = 300)
+#'             samples = 300,
+#'             chains = 2)
 #' np <- nuts_params(mod)
 #' head(np)
 #'

diff --git a/R/pairs.mvgam.R b/R/pairs.mvgam.R
@@ -17,8 +17,9 @@
 #' \donttest{
 #' simdat <- sim_mvgam(n_series = 1, trend_model = 'AR1')
 #' mod <- mvgam(y ~ s(season, bs = 'cc'),
-#'             trend_model = AR(),
-#'             data = simdat$data_train)
+#'              trend_model = AR(),
+#'              data = simdat$data_train,
+#'              chains = 2)
 #' pairs(mod)
 #' pairs(mod, variable = c('ar1', 'sigma'), regex = TRUE)
 #' }