Extract or compute hindcasts and forecasts for a fitted mvgam object

Extract or compute hindcasts and forecasts for a fitted mvgam object

Usage

forecast(object, ...)

# S3 method for mvgam
forecast(object, newdata, data_test, n_cores = 1, type = "response", ...)

Arguments

object: list object returned from mvgam. See mvgam()
...: Ignored
newdata: Optional dataframe or list of test data containing the same variables that were included in the original data used to fit the model. If included, the covariate information in newdata will be used to generate forecasts from the fitted model equations. If this same newdata was originally included in the call to mvgam, then forecasts have already been produced by the generative model and these will simply be extracted and plotted. However if no newdata was supplied to the original model call, an assumption is made that the newdata supplied here comes sequentially after the data supplied in the original model (i.e. we assume there is no time gap between the last observation of series 1 in the original data and the first observation for series 1 in newdata)
data_test: Deprecated. Still works in place of newdata but users are recommended to use newdata instead for more seamless integration into R workflows
n_cores: integer specifying number of cores for generating forecasts in parallel
type: When this has the value link (default) the linear predictor is calculated on the link scale. If expected is used, predictions reflect the expectation of the response (the mean) but ignore uncertainty in the observation process. When response is used, the predictions take uncertainty in the observation process into account to return predictions on the outcome scale. When variance is used, the variance of the response with respect to the mean (mean-variance relationship) is returned. When type = "terms", each component of the linear predictor is returned separately in the form of a list (possibly with standard errors, if summary = TRUE): this includes parametric model components, followed by each smooth component, but excludes any offset and any intercept. Two special cases are also allowed: type latent_N will return the estimated latent abundances from an N-mixture distribution, while type detection will return the estimated detection probability from an N-mixture distribution

Value

An object of class mvgam_forecast containing hindcast and forecast distributions. See mvgam_forecast-class for details.

Details

Posterior predictions are drawn from the fitted mvgam and used to simulate a forecast distribution

Examples

# \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,
            chains = 2)
#> Compiling Stan program using cmdstanr
#> 
<<<<<<< HEAD
#> 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/Rtmp2bnpq5/model-3cd012a31411.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
=======
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f
#> 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) 
<<<<<<< HEAD
=======
#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f
#> 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 2 Iteration: 400 / 1000 [ 40%]  (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: 500 / 1000 [ 50%]  (Warmup) 
#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
<<<<<<< HEAD
#> 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 1.1 seconds.
#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 1 finished in 1.2 seconds.
#> 
#> Both chains finished successfully.
#> Mean chain execution time: 1.1 seconds.
#> Total execution time: 1.3 seconds.
=======
#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 1 finished in 2.7 seconds.
#> Chain 2 finished in 2.6 seconds.
#> 
#> Both chains finished successfully.
#> Mean chain execution time: 2.7 seconds.
#> Total execution time: 2.8 seconds.
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f
#> 

# Hindcasts on response scale
hc <- hindcast(mod)
str(hc)
#> List of 15
#>  $ call              :Class 'formula'  language y ~ s(season, bs = "cc", k = 6)
<<<<<<< HEAD
#>   .. ..- attr(*, ".Environment")=<environment: 0x000001b48c2257f8> 
=======
#>   .. ..- attr(*, ".Environment")=<environment: 0x000001f1650b5778> 
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f
#>  $ trend_call        : NULL
#>  $ family            : chr "poisson"
#>  $ trend_model       :List of 7
#>   ..$ trend_model: chr "AR1"
#>   ..$ ma         : logi FALSE
#>   ..$ cor        : logi FALSE
#>   ..$ unit       : chr "time"
#>   ..$ gr         : chr "NA"
#>   ..$ subgr      : chr "series"
#>   ..$ label      : language AR()
#>   ..- attr(*, "class")= chr "mvgam_trend"
#>  $ drift             : logi FALSE
#>  $ use_lv            : logi FALSE
#>  $ fit_engine        : chr "stan"
#>  $ type              : chr "response"
#>  $ series_names      : chr [1:3] "series_1" "series_2" "series_3"
#>  $ train_observations:List of 3
<<<<<<< HEAD
#>   ..$ series_1: int [1:75] 1 0 0 1 0 0 1 1 1 4 ...
#>   ..$ series_2: int [1:75] 1 1 2 0 1 0 2 4 0 6 ...
#>   ..$ series_3: int [1:75] 1 2 2 0 0 0 3 0 0 3 ...
=======
#>   ..$ series_1: int [1:75] 0 1 2 4 0 1 13 6 1 3 ...
#>   ..$ series_2: int [1:75] 0 0 0 0 0 0 0 1 2 0 ...
#>   ..$ series_3: int [1:75] 0 0 0 0 0 1 2 2 1 0 ...
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f
#>  $ train_times       : int [1:75] 1 2 3 4 5 6 7 8 9 10 ...
#>  $ test_observations : NULL
#>  $ test_times        : NULL
#>  $ hindcasts         :List of 3
<<<<<<< HEAD
#>   ..$ series_1: num [1:1000, 1:75] 1 0 0 0 0 2 3 1 3 2 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,1]" "ypred[2,1]" "ypred[3,1]" "ypred[4,1]" ...
#>   ..$ series_2: num [1:1000, 1:75] 0 3 1 0 5 1 1 0 0 2 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,2]" "ypred[2,2]" "ypred[3,2]" "ypred[4,2]" ...
#>   ..$ series_3: num [1:1000, 1:75] 1 0 1 2 1 5 0 1 3 1 ...
=======
#>   ..$ series_1: num [1:1000, 1:75] 1 0 0 0 0 0 0 0 1 0 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,1]" "ypred[2,1]" "ypred[3,1]" "ypred[4,1]" ...
#>   ..$ series_2: num [1:1000, 1:75] 1 2 0 1 0 0 0 0 1 0 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,2]" "ypred[2,2]" "ypred[3,2]" "ypred[4,2]" ...
#>   ..$ series_3: num [1:1000, 1:75] 0 2 1 0 1 2 1 1 1 0 ...
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,3]" "ypred[2,3]" "ypred[3,3]" "ypred[4,3]" ...
#>  $ forecasts         : NULL
#>  - attr(*, "class")= chr "mvgam_forecast"
plot(hc, series = 1)
#> No non-missing values in test_observations; cannot calculate forecast score

plot(hc, series = 2)
#> No non-missing values in test_observations; cannot calculate forecast score

plot(hc, series = 3)
#> No non-missing values in test_observations; cannot calculate forecast score


# Forecasts on response scale
fc <- forecast(mod, newdata = simdat$data_test)
str(fc)
#> List of 16
#>  $ call              :Class 'formula'  language y ~ s(season, bs = "cc", k = 6)
<<<<<<< HEAD
#>   .. ..- attr(*, ".Environment")=<environment: 0x000001b48c2257f8> 
=======
#>   .. ..- attr(*, ".Environment")=<environment: 0x000001f1650b5778> 
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f
#>  $ trend_call        : NULL
#>  $ family            : chr "poisson"
#>  $ family_pars       : NULL
#>  $ trend_model       :List of 7
#>   ..$ trend_model: chr "AR1"
#>   ..$ ma         : logi FALSE
#>   ..$ cor        : logi FALSE
#>   ..$ unit       : chr "time"
#>   ..$ gr         : chr "NA"
#>   ..$ subgr      : chr "series"
#>   ..$ label      : language AR()
#>   ..- attr(*, "class")= chr "mvgam_trend"
#>  $ drift             : logi FALSE
#>  $ use_lv            : logi FALSE
#>  $ fit_engine        : chr "stan"
#>  $ type              : chr "response"
#>  $ series_names      : Factor w/ 3 levels "series_1","series_2",..: 1 2 3
#>  $ train_observations:List of 3
<<<<<<< HEAD
#>   ..$ series_1: int [1:75] 1 0 0 1 0 0 1 1 1 4 ...
#>   ..$ series_2: int [1:75] 1 1 2 0 1 0 2 4 0 6 ...
#>   ..$ series_3: int [1:75] 1 2 2 0 0 0 3 0 0 3 ...
#>  $ train_times       : int [1:75] 1 2 3 4 5 6 7 8 9 10 ...
#>  $ test_observations :List of 3
#>   ..$ series_1: int [1:25] 0 0 0 0 0 1 1 4 2 3 ...
#>   ..$ series_2: int [1:25] 1 0 0 2 0 3 2 9 8 1 ...
#>   ..$ series_3: int [1:25] 0 0 1 1 0 0 0 4 1 1 ...
#>  $ test_times        : int [1:25] 76 77 78 79 80 81 82 83 84 85 ...
#>  $ hindcasts         :List of 3
#>   ..$ series_1: num [1:1000, 1:75] 1 0 0 0 0 2 3 1 3 2 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,1]" "ypred[2,1]" "ypred[3,1]" "ypred[4,1]" ...
#>   ..$ series_2: num [1:1000, 1:75] 0 3 1 0 5 1 1 0 0 2 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,2]" "ypred[2,2]" "ypred[3,2]" "ypred[4,2]" ...
#>   ..$ series_3: num [1:1000, 1:75] 1 0 1 2 1 5 0 1 3 1 ...
=======
#>   ..$ series_1: int [1:75] 0 1 2 4 0 1 13 6 1 3 ...
#>   ..$ series_2: int [1:75] 0 0 0 0 0 0 0 1 2 0 ...
#>   ..$ series_3: int [1:75] 0 0 0 0 0 1 2 2 1 0 ...
#>  $ train_times       : int [1:75] 1 2 3 4 5 6 7 8 9 10 ...
#>  $ test_observations :List of 3
#>   ..$ series_1: int [1:25] 2 1 1 12 2 6 0 5 0 0 ...
#>   ..$ series_2: int [1:25] 0 0 0 10 0 0 1 1 0 0 ...
#>   ..$ series_3: int [1:25] 1 1 1 4 1 3 0 0 1 0 ...
#>  $ test_times        : int [1:25] 76 77 78 79 80 81 82 83 84 85 ...
#>  $ hindcasts         :List of 3
#>   ..$ series_1: num [1:1000, 1:75] 1 0 0 0 0 0 0 0 1 0 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,1]" "ypred[2,1]" "ypred[3,1]" "ypred[4,1]" ...
#>   ..$ series_2: num [1:1000, 1:75] 1 2 0 1 0 0 0 0 1 0 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,2]" "ypred[2,2]" "ypred[3,2]" "ypred[4,2]" ...
#>   ..$ series_3: num [1:1000, 1:75] 0 2 1 0 1 2 1 1 1 0 ...
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "ypred[1,3]" "ypred[2,3]" "ypred[3,3]" "ypred[4,3]" ...
#>  $ forecasts         :List of 3
<<<<<<< HEAD
#>   ..$ series_1: int [1:1000, 1:25] 0 2 0 0 1 0 0 1 0 2 ...
#>   ..$ series_2: int [1:1000, 1:25] 0 1 0 1 0 0 0 0 1 0 ...
#>   ..$ series_3: int [1:1000, 1:25] 0 1 0 1 0 0 1 1 0 0 ...
#>  - attr(*, "class")= chr "mvgam_forecast"
plot(fc, series = 1)
#> Out of sample DRPS:
#> 11.465246

plot(fc, series = 2)
#> Out of sample DRPS:
#> 25.378196

plot(fc, series = 3)
#> Out of sample DRPS:
#> 17.418698
=======
#>   ..$ series_1: int [1:1000, 1:25] 2 1 1 0 6 0 0 0 0 0 ...
#>   ..$ series_2: int [1:1000, 1:25] 1 2 1 0 0 2 0 0 0 1 ...
#>   ..$ series_3: int [1:1000, 1:25] 1 0 0 1 1 0 0 2 1 0 ...
#>  - attr(*, "class")= chr "mvgam_forecast"
plot(fc, series = 1)
#> Out of sample DRPS:
#> 40.733034

plot(fc, series = 2)
#> Out of sample DRPS:
#> 20.222428

plot(fc, series = 3)
#> Out of sample DRPS:
#> 13.041625
>>>>>>> 79e31b766e09ac39e22b93cb199f3eb4fda5e08f


# Forecasts as expectations
fc <- forecast(mod, newdata = simdat$data_test, type = 'expected')
plot(fc, series = 1)

plot(fc, series = 2)

plot(fc, series = 3)


# Dynamic trend extrapolations
fc <- forecast(mod, newdata = simdat$data_test, type = 'trend')
plot(fc, series = 1)

plot(fc, series = 2)

plot(fc, series = 3)

# }

Extract or compute hindcasts and forecasts for a fitted `mvgam` object

Usage

Arguments

Value

Details

See also

Examples