diff --git a/R/forecast.mvgam.R b/R/forecast.mvgam.R index bfe7984e..6b1ac9c0 100644 --- a/R/forecast.mvgam.R +++ b/R/forecast.mvgam.R @@ -28,13 +28,12 @@ forecast <- function(object, ...){ #'@rdname forecast.mvgam #'@method forecast mvgam #' @examples -#' \dontrun{ -#' simdat <- sim_mvgam(n_series = 3, trend_model = 'AR1') +#' \donttest{ +#' simdat <- sim_mvgam(n_series = 3, trend_model = AR()) #' mod <- mvgam(y ~ s(season, bs = 'cc', k = 6), #' trend_model = AR(), +#' noncentred = TRUE, #' data = simdat$data_train, -#' burnin = 300, -#' samples = 300, #' chains = 2) #' #' # Hindcasts on response scale diff --git a/R/hindcast.mvgam.R b/R/hindcast.mvgam.R index 0a8e49b2..7715fd78 100644 --- a/R/hindcast.mvgam.R +++ b/R/hindcast.mvgam.R @@ -7,7 +7,8 @@ #'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}} +#' +#'@seealso \code{\link{forecast.mvgam}} #'@export hindcast <- function(object, ...){ UseMethod("hindcast", object) @@ -16,13 +17,12 @@ hindcast <- function(object, ...){ #'@rdname hindcast.mvgam #'@method hindcast mvgam #' @examples -#' \dontrun{ -#' simdat <- sim_mvgam(n_series = 3, trend_model = 'AR1') +#' \donttest{ +#' simdat <- sim_mvgam(n_series = 3, trend_model = AR()) #' mod <- mvgam(y ~ s(season, bs = 'cc'), #' trend_model = AR(), +#' noncentred = TRUE, #' data = simdat$data_train, -#' burnin = 300, -#' samples = 300, #' chains = 2) #' #' # Hindcasts on response scale diff --git a/R/loo.mvgam.R b/R/loo.mvgam.R index 318eb1d8..799f3106 100644 --- a/R/loo.mvgam.R +++ b/R/loo.mvgam.R @@ -39,7 +39,7 @@ #' forecast evaluations for further scrutiny of models (see for example \code{\link{forecast.mvgam}}, #' \code{\link{score.mvgam_forecast}} and \code{\link{lfo_cv}}) #'@examples -#'\dontrun{ +#'\donttest{ #'# Simulate 4 time series with hierarchical seasonality #'# and independent AR1 dynamic processes #'set.seed(111) diff --git a/R/mcmc_plot.mvgam.R b/R/mcmc_plot.mvgam.R index 67a9b0db..66432bad 100644 --- a/R/mcmc_plot.mvgam.R +++ b/R/mcmc_plot.mvgam.R @@ -23,12 +23,11 @@ #' that can be used for argument `variable` #' @examples #' \dontrun{ -#' simdat <- sim_mvgam(n_series = 1, trend_model = 'AR1') +#' simdat <- sim_mvgam(n_series = 1, trend_model = AR()) #' mod <- mvgam(y ~ s(season, bs = 'cc', k = 6), #' trend_model = AR(), +#' noncentred = TRUE, #' data = simdat$data_train, -#' burnin = 300, -#' samples = 300, #' chains = 2) #' mcmc_plot(mod) #' mcmc_plot(mod, type = 'neff_hist') diff --git a/R/monotonic.R b/R/monotonic.R index 4da7cca8..ab9dc8ed 100644 --- a/R/monotonic.R +++ b/R/monotonic.R @@ -44,7 +44,7 @@ #' @author Nicholas J Clark #' @name monotonic #' @examples -#' \dontrun{ +#' \donttest{ #' # Simulate data from a monotonically increasing function #' set.seed(123123) #' x <- runif(80) * 4 - 1 @@ -72,8 +72,6 @@ #' mod2 <- mvgam(y ~ s(x, bs = 'moi', k = 18), #' data = mod_data, #' family = gaussian(), -#' burnin = 300, -#' samples = 300, #' chains = 2) #' #' plot_predictions(mod2, @@ -108,8 +106,6 @@ #' mod <- mvgam(y ~ s(x, bs = 'moi', by = fac, k = 8), #' data = mod_data, #' family = gaussian(), -#' burnin = 300, -#' samples = 300, #' chains = 2) #' #' # Visualise the different monotonic functions diff --git a/R/ti.R b/R/ti.R index 846e452c..086cb643 100644 --- a/R/ti.R +++ b/R/ti.R @@ -18,7 +18,7 @@ #' #' #' @examples -#' \dontrun{ +#' \donttest{ #' # Simulate some data #' dat <- mgcv::gamSim(1, n = 200, scale = 2) #' diff --git a/docs/index.html b/docs/index.html index 23a67632..fff8e221 100644 --- a/docs/index.html +++ b/docs/index.html @@ -114,14 +114,14 @@

Getting started Gamma() for non-negative real-valued data
  • -betar() for proportional data on (0,1) +betar() for proportional data on (0,1)
  • bernoulli() for binary data
  • poisson() for count data
  • -nb() for overdispersed count data
    • +nb() for overdispersed count data
  • binomial() for count data with known number of trials
  • @@ -131,9 +131,9 @@

    Getting started tweedie() for overdispersed count data

    • -

    Note that only poisson(), nb(), and tweedie() are available if using JAGS. All families, apart from tweedie(), are supported if using Stan. See ??mvgam_families for more information. Below is a simple example for simulating and modelling proportional data with Beta observations over a set of seasonal series with independent Gaussian Process dynamic trends:

    +

    Note that only poisson(), nb(), and tweedie() are available if using JAGS. All families, apart from tweedie(), are supported if using Stan. See ??mvgam_families for more information. Below is a simple example for simulating and modelling proportional data with Beta observations over a set of seasonal series with independent Gaussian Process dynamic trends:

    -data <- sim_mvgam(family = betar(),
+data <- sim_mvgam(family = betar(),
                   T = 80,
                   trend_model = GP(),
                   trend_rel = 0.5, 
@@ -150,7 +150,7 @@ 
    Getting started             trend_model = GP(),
              data = data$data_train,
              newdata = data$data_test,
-             family = betar())
    + family = betar())

    Plot the estimated posterior hindcast and forecast distributions for each series

     layout(matrix(1:4, nrow = 2, byrow = TRUE))
diff --git a/docs/reference/Rplot001.png b/docs/reference/Rplot001.png
index 17a35806..d15c0814 100644
Binary files a/docs/reference/Rplot001.png and b/docs/reference/Rplot001.png differ
diff --git a/docs/reference/Rplot002.png b/docs/reference/Rplot002.png
index f60122ed..d6904cf6 100644
Binary files a/docs/reference/Rplot002.png and b/docs/reference/Rplot002.png differ
diff --git a/docs/reference/Rplot003.png b/docs/reference/Rplot003.png
index 6408d480..68a7fb36 100644
Binary files a/docs/reference/Rplot003.png and b/docs/reference/Rplot003.png differ
diff --git a/docs/reference/Rplot004.png b/docs/reference/Rplot004.png
index f0f03fa2..d581bbe8 100644
Binary files a/docs/reference/Rplot004.png and b/docs/reference/Rplot004.png differ
diff --git a/docs/reference/Rplot005.png b/docs/reference/Rplot005.png
index 3680a107..97f7f6c1 100644
Binary files a/docs/reference/Rplot005.png and b/docs/reference/Rplot005.png differ
diff --git a/docs/reference/Rplot006.png b/docs/reference/Rplot006.png
index 4df17cb8..5d7e0baa 100644
Binary files a/docs/reference/Rplot006.png and b/docs/reference/Rplot006.png differ
diff --git a/docs/reference/Rplot007.png b/docs/reference/Rplot007.png
index 13a98490..c3cae136 100644
Binary files a/docs/reference/Rplot007.png and b/docs/reference/Rplot007.png differ
diff --git a/docs/reference/Rplot008.png b/docs/reference/Rplot008.png
index f8ee922c..d0851188 100644
Binary files a/docs/reference/Rplot008.png and b/docs/reference/Rplot008.png differ
diff --git a/docs/reference/Rplot009.png b/docs/reference/Rplot009.png
index d06943d1..5c33a2ac 100644
Binary files a/docs/reference/Rplot009.png and b/docs/reference/Rplot009.png differ
diff --git a/docs/reference/Rplot010.png b/docs/reference/Rplot010.png
index a277b510..64b403d1 100644
Binary files a/docs/reference/Rplot010.png and b/docs/reference/Rplot010.png differ
diff --git a/docs/reference/Rplot011.png b/docs/reference/Rplot011.png
index af9ef77f..5cf06ea5 100644
Binary files a/docs/reference/Rplot011.png and b/docs/reference/Rplot011.png differ
diff --git a/docs/reference/Rplot012.png b/docs/reference/Rplot012.png
index 19add981..9adea771 100644
Binary files a/docs/reference/Rplot012.png and b/docs/reference/Rplot012.png differ
diff --git a/docs/reference/Rplot013.png b/docs/reference/Rplot013.png
index 5b5b97bd..bf4f53c9 100644
Binary files a/docs/reference/Rplot013.png and b/docs/reference/Rplot013.png differ
diff --git a/docs/reference/Rplot014.png b/docs/reference/Rplot014.png
index c4806d27..acf3ee19 100644
Binary files a/docs/reference/Rplot014.png and b/docs/reference/Rplot014.png differ
diff --git a/docs/reference/Rplot015.png b/docs/reference/Rplot015.png
index 281e787f..53c008cd 100644
Binary files a/docs/reference/Rplot015.png and b/docs/reference/Rplot015.png differ
diff --git a/docs/reference/Rplot016.png b/docs/reference/Rplot016.png
index 94d24362..6fe4c62d 100644
Binary files a/docs/reference/Rplot016.png and b/docs/reference/Rplot016.png differ
diff --git a/docs/reference/conditional_effects.mvgam-1.png b/docs/reference/conditional_effects.mvgam-1.png
index 616192ec..79c60ce8 100644
Binary files a/docs/reference/conditional_effects.mvgam-1.png and b/docs/reference/conditional_effects.mvgam-1.png differ
diff --git a/docs/reference/conditional_effects.mvgam-10.png b/docs/reference/conditional_effects.mvgam-10.png
index 45178090..1bd8b816 100644
Binary files a/docs/reference/conditional_effects.mvgam-10.png and b/docs/reference/conditional_effects.mvgam-10.png differ
diff --git a/docs/reference/conditional_effects.mvgam-11.png b/docs/reference/conditional_effects.mvgam-11.png
index e5689162..d67a627a 100644
Binary files a/docs/reference/conditional_effects.mvgam-11.png and b/docs/reference/conditional_effects.mvgam-11.png differ
diff --git a/docs/reference/conditional_effects.mvgam-12.png b/docs/reference/conditional_effects.mvgam-12.png
index 4aaf533d..115b3890 100644
Binary files a/docs/reference/conditional_effects.mvgam-12.png and b/docs/reference/conditional_effects.mvgam-12.png differ
diff --git a/docs/reference/conditional_effects.mvgam-2.png b/docs/reference/conditional_effects.mvgam-2.png
index 29840bc7..c3583f7c 100644
Binary files a/docs/reference/conditional_effects.mvgam-2.png and b/docs/reference/conditional_effects.mvgam-2.png differ
diff --git a/docs/reference/conditional_effects.mvgam-3.png b/docs/reference/conditional_effects.mvgam-3.png
index c8e36719..c4534ce0 100644
Binary files a/docs/reference/conditional_effects.mvgam-3.png and b/docs/reference/conditional_effects.mvgam-3.png differ
diff --git a/docs/reference/conditional_effects.mvgam-4.png b/docs/reference/conditional_effects.mvgam-4.png
index c235c2e8..9be3e83b 100644
Binary files a/docs/reference/conditional_effects.mvgam-4.png and b/docs/reference/conditional_effects.mvgam-4.png differ
diff --git a/docs/reference/conditional_effects.mvgam-5.png b/docs/reference/conditional_effects.mvgam-5.png
index d6f3b4fc..ecf6000c 100644
Binary files a/docs/reference/conditional_effects.mvgam-5.png and b/docs/reference/conditional_effects.mvgam-5.png differ
diff --git a/docs/reference/conditional_effects.mvgam-6.png b/docs/reference/conditional_effects.mvgam-6.png
index 918dcfde..2ce75d3e 100644
Binary files a/docs/reference/conditional_effects.mvgam-6.png and b/docs/reference/conditional_effects.mvgam-6.png differ
diff --git a/docs/reference/conditional_effects.mvgam-7.png b/docs/reference/conditional_effects.mvgam-7.png
index 5bc7e7a7..18c43fd6 100644
Binary files a/docs/reference/conditional_effects.mvgam-7.png and b/docs/reference/conditional_effects.mvgam-7.png differ
diff --git a/docs/reference/conditional_effects.mvgam-8.png b/docs/reference/conditional_effects.mvgam-8.png
index 2d7268f3..9943ca38 100644
Binary files a/docs/reference/conditional_effects.mvgam-8.png and b/docs/reference/conditional_effects.mvgam-8.png differ
diff --git a/docs/reference/conditional_effects.mvgam-9.png b/docs/reference/conditional_effects.mvgam-9.png
index 2a11b0e5..61b3fb8b 100644
Binary files a/docs/reference/conditional_effects.mvgam-9.png and b/docs/reference/conditional_effects.mvgam-9.png differ
diff --git a/docs/reference/conditional_effects.mvgam.html b/docs/reference/conditional_effects.mvgam.html
index efbff23c..663477e7 100644
--- a/docs/reference/conditional_effects.mvgam.html
+++ b/docs/reference/conditional_effects.mvgam.html
@@ -168,7 +168,7 @@ 
    Author<
 
     
    
     
    Examples
    
-    
    if (FALSE) {
+    
    # \donttest{
 # Simulate some data
 simdat <- sim_mvgam(family = poisson(),
                     seasonality = 'hierarchical')
@@ -177,32 +177,139 @@ 
    Examplesmod <- mvgam(y ~ s(season, by = series, k = 5) + year:series,
              family = poisson(),
              data = simdat$data_train,
-             chains = 2,
-             burnin = 300,
-             samples = 300)
+             chains = 2)
+#> Compiling Stan program using cmdstanr
+#> 
+#> 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/RtmpU5rG2n/model-3fc011f623e7.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
+#> Start sampling
+#> Running MCMC with 2 parallel chains...
+#> 
+#> Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 1 finished in 0.7 seconds.
+#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 finished in 0.8 seconds.
+#> 
+#> Both chains finished successfully.
+#> Mean chain execution time: 0.8 seconds.
+#> Total execution time: 1.0 seconds.
+#> 
 
 # Plot all main effects on the response scale
 conditional_effects(mod)
+
+
 
 # Change the prediction interval to 70% using plot_predictions() argument
 # 'conf_level'
 conditional_effects(mod, conf_level = 0.7)
+
+
 
 # Plot all main effects on the link scale
 conditional_effects(mod, type = 'link')
+
+
 
 # Works the same for smooth terms, including smooth interactions
 set.seed(0)
 dat <- mgcv::gamSim(1, n = 200, scale = 2)
+#> Gu & Wahba 4 term additive model
 mod <- mvgam(y ~ te(x0, x1, k = 5) + s(x2, k = 6) + s(x3, k = 6),
             data = dat,
             family = gaussian(),
-            chains = 2,
-            burnin = 300,
-            samples = 300)
+            chains = 2)
+#> Compiling Stan program using cmdstanr
+#> 
+#> 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/RtmpU5rG2n/model-3fc072543cb.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
+#> Start sampling
+#> Running MCMC with 2 parallel chains...
+#> 
+#> Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 1 finished in 2.5 seconds.
+#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 finished in 2.9 seconds.
+#> 
+#> Both chains finished successfully.
+#> Mean chain execution time: 2.7 seconds.
+#> Total execution time: 3.0 seconds.
+#> 
 conditional_effects(mod)
+
+
+
 conditional_effects(mod, conf_level = 0.5, type = 'link')
-}
+
+
+
+# }
 
    
     
    
   
    
@@ -148,69 +148,174 @@ 
    References
     
    Examples
    
-    
    if (FALSE) {
-# 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,
-             chains = 2)
-summary(mod)
-
-# Plot the posterior hindcast
-plot(mod, type = 'forecast')
-
-# View the changepoints with ggplot2 utilities
-library(ggplot2)
-mcmc_plot(mod, variable = 'delta_trend',
-          regex = TRUE) +
-scale_y_discrete(labels = mod$trend_model$changepoints) +
-labs(y = 'Potential changepoint',
-     x = 'Rate change')
-}
    
+    
    # \donttest{
+# 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,
+             chains = 2)
+#> Compiling Stan program using cmdstanr
+#> 
+#> 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/RtmpU5rG2n/model-3fc0236235ea.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
+#> Start sampling
+#> Running MCMC with 2 parallel chains...
+#> 
+#> Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 finished in 7.1 seconds.
+#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 1 finished in 7.3 seconds.
+#> 
+#> Both chains finished successfully.
+#> Mean chain execution time: 7.2 seconds.
+#> Total execution time: 7.4 seconds.
+#> 
+summary(mod)
+#> GAM formula:
+#> y ~ 1
+#> <environment: 0x000001ad823a6de8>
+#> 
+#> 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 = 1000; warmup = 500; thin = 1 
+#> Total post-warmup draws = 1000
+#> 
+#> 
+#> 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.26 -0.15 -0.071 1.01   348
+#> 
+#> Latent trend offset estimates:
+#>            2.5%  50% 97.5% Rhat n_eff
+#> m_trend[1]  -14 -4.3 -0.36 1.01   347
+#> 
+#> Stan MCMC diagnostics:
+#> n_eff / iter looks reasonable for all parameters
+#> Rhat looks reasonable for all parameters
+#> 0 of 1000 iterations ended with a divergence (0%)
+#> 0 of 1000 iterations saturated the maximum tree depth of 12 (0%)
+#> E-FMI indicated no pathological behavior
+#> 
+#> Samples were drawn using NUTS(diag_e) at Mon Jul 01 8:47:11 AM 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)
+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.
+
+# }
+
    diff --git a/docs/reference/plot.mvgam-1.png b/docs/reference/plot.mvgam-1.png index 37c2f527..794dbd1d 100644 Binary files a/docs/reference/plot.mvgam-1.png and b/docs/reference/plot.mvgam-1.png differ diff --git a/docs/reference/plot.mvgam-10.png b/docs/reference/plot.mvgam-10.png index d6a7cc24..0a5e0fef 100644 Binary files a/docs/reference/plot.mvgam-10.png and b/docs/reference/plot.mvgam-10.png differ diff --git a/docs/reference/plot.mvgam-11.png b/docs/reference/plot.mvgam-11.png index cd140256..3cb3ec9f 100644 Binary files a/docs/reference/plot.mvgam-11.png and b/docs/reference/plot.mvgam-11.png differ diff --git a/docs/reference/plot.mvgam-12.png b/docs/reference/plot.mvgam-12.png index dddbf79b..73b452d0 100644 Binary files a/docs/reference/plot.mvgam-12.png and b/docs/reference/plot.mvgam-12.png differ diff --git a/docs/reference/plot.mvgam-13.png b/docs/reference/plot.mvgam-13.png index 12a5dc0a..890b0cc9 100644 Binary files a/docs/reference/plot.mvgam-13.png and b/docs/reference/plot.mvgam-13.png differ diff --git a/docs/reference/plot.mvgam-2.png b/docs/reference/plot.mvgam-2.png index fcc00405..f6cb46ac 100644 Binary files a/docs/reference/plot.mvgam-2.png and b/docs/reference/plot.mvgam-2.png differ diff --git a/docs/reference/plot.mvgam-3.png b/docs/reference/plot.mvgam-3.png index 788bd2a6..e8253a0d 100644 Binary files a/docs/reference/plot.mvgam-3.png and b/docs/reference/plot.mvgam-3.png differ diff --git a/docs/reference/plot.mvgam-4.png b/docs/reference/plot.mvgam-4.png index 2aaa018b..a1b450ce 100644 Binary files a/docs/reference/plot.mvgam-4.png and b/docs/reference/plot.mvgam-4.png differ diff --git a/docs/reference/plot.mvgam-5.png b/docs/reference/plot.mvgam-5.png index 314c0665..601865c1 100644 Binary files a/docs/reference/plot.mvgam-5.png and b/docs/reference/plot.mvgam-5.png differ diff --git a/docs/reference/plot.mvgam-6.png b/docs/reference/plot.mvgam-6.png index 7ab3b635..fffb211a 100644 Binary files a/docs/reference/plot.mvgam-6.png and b/docs/reference/plot.mvgam-6.png differ diff --git a/docs/reference/plot.mvgam-7.png b/docs/reference/plot.mvgam-7.png index c44589d9..fa9a479c 100644 Binary files a/docs/reference/plot.mvgam-7.png and b/docs/reference/plot.mvgam-7.png differ diff --git a/docs/reference/plot.mvgam-8.png b/docs/reference/plot.mvgam-8.png index 343496ed..f7919b7c 100644 Binary files a/docs/reference/plot.mvgam-8.png and b/docs/reference/plot.mvgam-8.png differ diff --git a/docs/reference/plot.mvgam-9.png b/docs/reference/plot.mvgam-9.png index cf09993a..5f8c3d51 100644 Binary files a/docs/reference/plot.mvgam-9.png and b/docs/reference/plot.mvgam-9.png differ diff --git a/docs/reference/plot.mvgam.html b/docs/reference/plot.mvgam.html index 708e0652..9f4ae945 100644 --- a/docs/reference/plot.mvgam.html +++ b/docs/reference/plot.mvgam.html @@ -163,52 +163,163 @@

    Author<

    Examples

    - 
    if (FALSE) {
-# Simulate some time series
-dat <- sim_mvgam(T = 80, n_series = 3)
-
-# Fit a basic model
-mod <- mvgam(y ~ s(season, bs = 'cc') + s(series, bs = 're'),
-            data = dat$data_train,
-            trend_model = RW(),
-            burnin = 300,
-            samples = 300,
-            chains = 2)
-
-# Plot predictions and residuals for each series
-plot(mod, type = 'forecast', series = 1)
-plot(mod, type = 'forecast', series = 2)
-plot(mod, type = 'forecast', series = 3)
-plot(mod, type = 'residuals', series = 1)
-plot(mod, type = 'residuals', series = 2)
-plot(mod, type = 'residuals', series = 3)
-
-# Plot model effects
-plot(mod, type = 'smooths')
-plot(mod, type = 're')
-
-# More flexible plots with 'marginaleffects' utilities
-plot_predictions(mod, condition = 'season', type = 'link')
-plot_predictions(mod,
-                condition = c('season', 'series', 'series'),
-                type = 'link')
-plot_predictions(mod, condition = 'series', type = 'link')
-
-# When using a State-Space model with predictors on the process
-# model, set trend_effects = TRUE to visualise process effects
-mod <- mvgam(y ~ -1,
-            trend_formula = ~ s(season, bs = 'cc'),
-            data = dat$data_train,
-            trend_model = RW(),
-            burnin = 300,
-            samples = 300,
-            chains = 2)
-plot(mod, type = 'smooths', trend_effects = TRUE)
-
-# But marginaleffects functions work without any modification
-plot_predictions(mod, condition = 'season', type = 'link')
-
-}
    + 
    # \donttest{
+# Simulate some time series
+dat <- sim_mvgam(T = 80, n_series = 3)
+
+# Fit a basic model
+mod <- mvgam(y ~ s(season, bs = 'cc') + s(series, bs = 're'),
+            data = dat$data_train,
+            trend_model = RW(),
+            chains = 2)
+#> Your model may benefit from using "noncentred = TRUE"
+#> Compiling Stan program using cmdstanr
+#> 
+#> 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/RtmpU5rG2n/model-3fc04a6c5d12.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
+#> Start sampling
+#> Running MCMC with 2 parallel chains...
+#> 
+#> Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 1 finished in 9.0 seconds.
+#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 finished in 14.2 seconds.
+#> 
+#> Both chains finished successfully.
+#> Mean chain execution time: 11.6 seconds.
+#> Total execution time: 14.3 seconds.
+#> 
+
+# Plot predictions and residuals for each series
+plot(mod, type = 'forecast', series = 1)
+
+plot(mod, type = 'forecast', series = 2)
+
+plot(mod, type = 'forecast', series = 3)
+
+plot(mod, type = 'residuals', series = 1)
+
+plot(mod, type = 'residuals', series = 2)
+
+plot(mod, type = 'residuals', series = 3)
+
+
+# Plot model effects
+plot(mod, type = 'smooths')
+
+plot(mod, type = 're')
+
+
+# More flexible plots with 'marginaleffects' utilities
+library(marginaleffects)
+plot_predictions(mod, condition = 'season', type = 'link')
+
+plot_predictions(mod,
+                condition = c('season', 'series', 'series'),
+                type = 'link')
+
+plot_predictions(mod, condition = 'series', type = 'link')
+
+
+# When using a State-Space model with predictors on the process
+# model, set trend_effects = TRUE to visualise process effects
+mod <- mvgam(y ~ -1,
+            trend_formula = ~ s(season, bs = 'cc'),
+            data = dat$data_train,
+            trend_model = RW(),
+            chains = 2)
+#> Your model may benefit from using "noncentred = TRUE"
+#> Compiling Stan program using cmdstanr
+#> 
+#> 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/RtmpU5rG2n/model-3fc059bc5efa.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
+#> Start sampling
+#> Running MCMC with 2 parallel chains...
+#> 
+#> Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 finished in 3.6 seconds.
+#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 1 finished in 3.6 seconds.
+#> 
+#> Both chains finished successfully.
+#> Mean chain execution time: 3.6 seconds.
+#> Total execution time: 3.7 seconds.
+#> 
+plot(mod, type = 'smooths', trend_effects = TRUE)
+
+
+# But marginaleffects functions work without any modification
+plot_predictions(mod, condition = 'season', type = 'link')
+
+
+# }
+
    diff --git a/docs/reference/plot_mvgam_factors-1.png b/docs/reference/plot_mvgam_factors-1.png index 894ac9f4..540fa862 100644 Binary files a/docs/reference/plot_mvgam_factors-1.png and b/docs/reference/plot_mvgam_factors-1.png differ diff --git a/docs/reference/plot_mvgam_factors.html b/docs/reference/plot_mvgam_factors.html index ebb54ec5..a109be3a 100644 --- a/docs/reference/plot_mvgam_factors.html +++ b/docs/reference/plot_mvgam_factors.html @@ -1,137 +1,207 @@ - -Latent factor summaries for a fitted mvgam object — plot_mvgam_factors • mvgam - Skip to contents - - -
    -
    -
    -

    Latent factor summaries for a fitted mvgam object

    - Source: R/plot_mvgam_factors.R -
    plot_mvgam_factors.Rd
    -
    - -
    -

    This function takes a fitted mvgam object and returns plots and summary statistics for -the latent dynamic factors

    -
    - -
    -

    Usage

    - 
    plot_mvgam_factors(object, plot = TRUE)
    -
    - -
    -

    Arguments

    -
    object
    -

    list object returned from mvgam. See mvgam()

    - - -
    plot
    -

    logical specifying whether factors should be plotted

    - -
    -
    -

    Value

    - - -

    A dataframe of factor contributions and, -optionally, a series of base R plots

    -
    -
    -

    Details

    -

    If the model in object was estimated using dynamic factors, it is possible that not all factors -contributed to the estimated trends. This is due to the regularisation penalty that acts independently on each -factor's Gaussian precision, which will squeeze un-needed factors to a white noise process (effectively dropping -that factor from the model). In this function, each factor is tested against a null hypothesis of white noise by -calculating the sum of the factor's 2nd derivatives. A factor that has a larger contribution will have a larger -sum due to the weaker penalty on the factor's precision. If -plot == TRUE, the factors are also plotted.

    -
    -
    -

    Author

    -

    Nicholas J Clark

    -
    - -
    -

    Examples

    - 
    if (FALSE) {
-simdat <- sim_mvgam()
-mod <- mvgam(y ~ s(season, bs = 'cc',
-                  k = 6),
-            trend_model = AR(),
-            use_lv = TRUE,
-            n_lv = 2,
-            data = simdat$data_train,
-            chains = 2)
-plot_mvgam_factors(mod)
-}
    -
    -
    - - -
    - - - - - - - + +Latent factor summaries for a fitted mvgam object — plot_mvgam_factors • mvgam + Skip to contents + + +
    +
    +
    +

    Latent factor summaries for a fitted mvgam object

    + Source: R/plot_mvgam_factors.R +
    plot_mvgam_factors.Rd
    +
    + +
    +

    This function takes a fitted mvgam object and returns plots and summary statistics for +the latent dynamic factors

    +
    + +
    +

    Usage

    + 
    plot_mvgam_factors(object, plot = TRUE)
    +
    + +
    +

    Arguments

    +
    object
    +

    list object returned from mvgam. See mvgam()

    + + +
    plot
    +

    logical specifying whether factors should be plotted

    + +
    +
    +

    Value

    + + +

    A dataframe of factor contributions and, +optionally, a series of base R plots

    +
    +
    +

    Details

    +

    If the model in object was estimated using dynamic factors, it is possible that not all factors +contributed to the estimated trends. This is due to the regularisation penalty that acts independently on each +factor's Gaussian precision, which will squeeze un-needed factors to a white noise process (effectively dropping +that factor from the model). In this function, each factor is tested against a null hypothesis of white noise by +calculating the sum of the factor's 2nd derivatives. A factor that has a larger contribution will have a larger +sum due to the weaker penalty on the factor's precision. If +plot == TRUE, the factors are also plotted.

    +
    +
    +

    Author

    +

    Nicholas J Clark

    +
    + +
    +

    Examples

    + 
    # \donttest{
+simdat <- sim_mvgam()
+mod <- mvgam(y ~ s(season, bs = 'cc',
+                  k = 6),
+            trend_model = AR(),
+            use_lv = TRUE,
+            n_lv = 2,
+            data = simdat$data_train,
+            chains = 2)
+#> Warning in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpU5rG2n/model_b6696a4e4d5e19ee8abb2233ebec3436.stan', line 23, column 31: Found

+#>     int division:

+#>       n_lv * (n_lv - 1) / 2

+#>     Values will be rounded towards zero. If rounding is not desired you can

+#>     write

+#>     the division as

+#>       n_lv * (n_lv - 1) / 2.0

+#>     If rounding is intended please use the integer division operator %/%.

+#> Compiling Stan program using cmdstanr
+#> 
+#> Warning in 'C:/Users/uqnclar2/AppData/Local/Temp/RtmpU5rG2n/model-3fc01824595c.stan', line 23, column 33: Found

+#>     int division:

+#>       n_lv * (n_lv - 1) / 2

+#>     Values will be rounded towards zero. If rounding is not desired you can

+#>     write

+#>     the division as

+#>       n_lv * (n_lv - 1) / 2.0

+#>     If rounding is intended please use the integer division operator %/%.

+#> 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/RtmpU5rG2n/model-3fc01824595c.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
+#> Start sampling
+#> Running MCMC with 2 parallel chains...
+#> 
+#> Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 1 finished in 1.3 seconds.
+#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 finished in 1.3 seconds.
+#> 
+#> Both chains finished successfully.
+#> Mean chain execution time: 1.3 seconds.
+#> Total execution time: 1.4 seconds.
+#> 
+plot_mvgam_factors(mod)
+
+#>         Contribution
+#> Factor1    0.5061927
+#> Factor2    0.4938073
+# }
+
    +
    +
    + + +
    + + + + + + + diff --git a/docs/reference/plot_mvgam_series-1.png b/docs/reference/plot_mvgam_series-1.png new file mode 100644 index 00000000..c051038d Binary files /dev/null and b/docs/reference/plot_mvgam_series-1.png differ diff --git a/docs/reference/plot_mvgam_series-2.png b/docs/reference/plot_mvgam_series-2.png new file mode 100644 index 00000000..9e4ab458 Binary files /dev/null and b/docs/reference/plot_mvgam_series-2.png differ diff --git a/docs/reference/plot_mvgam_series-3.png b/docs/reference/plot_mvgam_series-3.png new file mode 100644 index 00000000..8eb2b0b9 Binary files /dev/null and b/docs/reference/plot_mvgam_series-3.png differ diff --git a/docs/reference/plot_mvgam_series.html b/docs/reference/plot_mvgam_series.html index 64365ba3..309bdfdf 100644 --- a/docs/reference/plot_mvgam_series.html +++ b/docs/reference/plot_mvgam_series.html @@ -153,14 +153,21 @@

    Author<

    Examples

    - 
    # Simulate and plot series with observations bounded at 0 and 1 (Beta responses)
-sim_data <- sim_mvgam(family = betar(), trend_model = RW(), prop_trend = 0.6)
-plot_mvgam_series(data = sim_data$data_train, series = 'all')
-plot_mvgam_series(data = sim_data$data_train, newdata = sim_data$data_test, series = 1)
-
-# Now simulate series with overdispersed discrete observations
-sim_data <- sim_mvgam(family = nb(), trend_model = RW(), prop_trend = 0.6, phi = 10)
-plot_mvgam_series(data = sim_data$data_train, series = 'all')
    + 
    # Simulate and plot series with observations bounded at 0 and 1 (Beta responses)
+sim_data <- sim_mvgam(family = betar(),
+                     trend_model = RW(), prop_trend = 0.6)
+plot_mvgam_series(data = sim_data$data_train, series = 'all')
+
+plot_mvgam_series(data = sim_data$data_train,
+                 newdata = sim_data$data_test, series = 1)
+
+
+# Now simulate series with overdispersed discrete observations
+sim_data <- sim_mvgam(family = nb(), trend_model = RW(),
+                     prop_trend = 0.6, phi = 10)
+plot_mvgam_series(data = sim_data$data_train, series = 'all')
+
+
    diff --git a/docs/reference/plot_mvgam_trend-1.png b/docs/reference/plot_mvgam_trend-1.png index 382a6df8..10cfe97a 100644 Binary files a/docs/reference/plot_mvgam_trend-1.png and b/docs/reference/plot_mvgam_trend-1.png differ diff --git a/docs/reference/plot_mvgam_trend-2.png b/docs/reference/plot_mvgam_trend-2.png index e095608d..048246d8 100644 Binary files a/docs/reference/plot_mvgam_trend-2.png and b/docs/reference/plot_mvgam_trend-2.png differ diff --git a/docs/reference/plot_mvgam_trend-3.png b/docs/reference/plot_mvgam_trend-3.png index 974097e7..29ebf54a 100644 Binary files a/docs/reference/plot_mvgam_trend-3.png and b/docs/reference/plot_mvgam_trend-3.png differ diff --git a/docs/reference/plot_mvgam_trend-4.png b/docs/reference/plot_mvgam_trend-4.png index 812c5189..26ced11b 100644 Binary files a/docs/reference/plot_mvgam_trend-4.png and b/docs/reference/plot_mvgam_trend-4.png differ diff --git a/docs/reference/plot_mvgam_trend-5.png b/docs/reference/plot_mvgam_trend-5.png index a5807e28..550b7454 100644 Binary files a/docs/reference/plot_mvgam_trend-5.png and b/docs/reference/plot_mvgam_trend-5.png differ diff --git a/docs/reference/plot_mvgam_trend-6.png b/docs/reference/plot_mvgam_trend-6.png index b49a315c..5c16cdc3 100644 Binary files a/docs/reference/plot_mvgam_trend-6.png and b/docs/reference/plot_mvgam_trend-6.png differ diff --git a/docs/reference/plot_mvgam_trend.html b/docs/reference/plot_mvgam_trend.html index 201bbef6..8cab1961 100644 --- a/docs/reference/plot_mvgam_trend.html +++ b/docs/reference/plot_mvgam_trend.html @@ -145,29 +145,75 @@

    Value

    Examples

    - 
    # \donttest{
-simdat <- sim_mvgam(n_series = 3, trend_model = 'AR1')
-mod <- mvgam(y ~ s(season, bs = 'cc', k = 6),
-            trend_model = AR(),
-            data = simdat$data_train,
-            burnin = 300,
-            samples = 300,
-            chains = 2)
-
-# Plot estimated trends for some series
-plot_mvgam_trend(mod)
-plot_mvgam_trend(mod, series = 2)
-
-# Extrapolate trends forward in time and plot on response scale
-plot_mvgam_trend(mod, newdata = simdat$data_test)
-plot_mvgam_trend(mod, newdata = simdat$data_test, series = 2)
-
-# But it is recommended to compute extrapolations for all series
-# first and then plot
-trend_fc <- forecast(mod, newdata = simdat$data_test)
-plot(trend_fc, series = 1)
-plot(trend_fc, series = 2)
-# }
    + 
    # \donttest{
+simdat <- sim_mvgam(n_series = 3, trend_model = 'AR1')
+mod <- mvgam(y ~ s(season, bs = 'cc', k = 6),
+            trend_model = AR(),
+            noncentred = TRUE,
+            data = simdat$data_train,
+            chains = 2)
+#> Compiling Stan program using cmdstanr
+#> 
+#> Start sampling
+#> Running MCMC with 2 parallel chains...
+#> 
+#> Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
+#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
+#> Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
+#> Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
+#> Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
+#> Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
+#> Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
+#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
+#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
+#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
+#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
+#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 2 finished in 1.0 seconds.
+#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
+#> Chain 1 finished in 1.1 seconds.
+#> 
+#> Both chains finished successfully.
+#> Mean chain execution time: 1.0 seconds.
+#> Total execution time: 1.2 seconds.
+#> 
+
+# Plot estimated trends for some series
+plot_mvgam_trend(mod)
+
+plot_mvgam_trend(mod, series = 2)
+
+
+# Extrapolate trends forward in time and plot on response scale
+plot_mvgam_trend(mod, newdata = simdat$data_test)
+
+plot_mvgam_trend(mod, newdata = simdat$data_test, series = 2)
+
+
+# But it is recommended to compute extrapolations for all series
+# first and then plot
+trend_fc <- forecast(mod, newdata = simdat$data_test)
+plot(trend_fc, series = 1)
+#> Out of sample DRPS:
+#> 20.849655
+
+plot(trend_fc, series = 2)
+#> Out of sample DRPS:
+#> 25.91199
+
+# }
+
    diff --git a/docs/reference/posterior_epred.mvgam.html b/docs/reference/posterior_epred.mvgam.html index 18951486..6e4fc6a6 100644 --- a/docs/reference/posterior_epred.mvgam.html +++ b/docs/reference/posterior_epred.mvgam.html @@ -28,7 +28,7 @@ mvgam - 1.1.0 + 1.1.2