Extract hindcasts for a fitted mvgam object

Extract hindcasts for a fitted mvgam object

Usage

hindcast(object, ...)

# S3 method for mvgam
hindcast(object, type = "response", ...)

Arguments

object: list object returned from mvgam. See mvgam()
...: Ignored
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 distributions. See mvgam_forecast-class for details.

Details

Posterior retrodictions are drawn from the fitted mvgam and organized into a convenient format

Examples

# \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,
            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 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
#> Chain 2 Iteration: 100 / 1000 [ 10%]  (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 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 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
#> Chain 1 finished in 2.7 seconds.
#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 2 finished in 3.0 seconds.
#> 
#> Both chains finished successfully.
#> Mean chain execution time: 2.9 seconds.
#> Total execution time: 3.2 seconds.
#> 

# Hindcasts on response scale
hc <- hindcast(mod)
str(hc)
#> List of 15
#>  $ call              :Class 'formula'  language y ~ s(season, bs = "cc")
#>   .. ..- attr(*, ".Environment")=<environment: 0x000001341428cdf8> 
#>  $ trend_call        : NULL
#>  $ family            : chr "poisson"
#>  $ trend_model       :List of 4
#>   ..$ trend_model: chr "AR1"
#>   ..$ ma         : logi FALSE
#>   ..$ cor        : logi FALSE
#>   ..$ 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
#>   ..$ series_1: int [1:75] 1 0 1 1 1 0 2 1 5 0 ...
#>   ..$ series_2: int [1:75] 0 0 1 1 1 0 3 5 5 1 ...
#>   ..$ series_3: int [1:75] 1 1 2 2 0 2 9 9 2 2 ...
#>  $ 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
#>   ..$ series_1: num [1:1000, 1:75] 1 5 0 1 1 5 1 1 2 1 ...
#>   .. ..- 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] 5 0 5 1 0 0 1 0 2 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] 4 1 2 1 0 2 2 1 3 3 ...
#>   .. ..- 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


# Hindcasts as expectations
hc <- hindcast(mod, type = 'expected')
str(hc)
#> List of 15
#>  $ call              :Class 'formula'  language y ~ s(season, bs = "cc")
#>   .. ..- attr(*, ".Environment")=<environment: 0x000001341428cdf8> 
#>  $ trend_call        : NULL
#>  $ family            : chr "poisson"
#>  $ trend_model       :List of 4
#>   ..$ trend_model: chr "AR1"
#>   ..$ ma         : logi FALSE
#>   ..$ cor        : logi FALSE
#>   ..$ label      : language AR()
#>   ..- attr(*, "class")= chr "mvgam_trend"
#>  $ drift             : logi FALSE
#>  $ use_lv            : logi FALSE
#>  $ fit_engine        : chr "stan"
#>  $ type              : chr "expected"
#>  $ series_names      : chr [1:3] "series_1" "series_2" "series_3"
#>  $ train_observations:List of 3
#>   ..$ series_1: int [1:75] 1 0 1 1 1 0 2 1 5 0 ...
#>   ..$ series_2: int [1:75] 0 0 1 1 1 0 3 5 5 1 ...
#>   ..$ series_3: int [1:75] 1 1 2 2 0 2 9 9 2 2 ...
#>  $ 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
#>   ..$ series_1: num [1:1000, 1:75] 1.61 3.77 0.77 2.87 1.85 ...
#>   ..$ series_2: num [1:1000, 1:75] 1.801 1.27 2.474 0.896 1.998 ...
#>   ..$ series_3: num [1:1000, 1:75] 2.12 1.54 1.97 2.21 1.42 ...
#>  $ forecasts         : NULL
#>  - attr(*, "class")= chr "mvgam_forecast"
plot(hc, series = 1)

plot(hc, series = 2)

plot(hc, series = 3)


# Estimated latent trends
hc <- hindcast(mod, type = 'trend')
str(hc)
#> List of 15
#>  $ call              :Class 'formula'  language y ~ s(season, bs = "cc")
#>   .. ..- attr(*, ".Environment")=<environment: 0x000001341428cdf8> 
#>  $ trend_call        : NULL
#>  $ family            : chr "poisson"
#>  $ trend_model       :List of 4
#>   ..$ trend_model: chr "AR1"
#>   ..$ ma         : logi FALSE
#>   ..$ cor        : logi FALSE
#>   ..$ label      : language AR()
#>   ..- attr(*, "class")= chr "mvgam_trend"
#>  $ drift             : logi FALSE
#>  $ use_lv            : logi FALSE
#>  $ fit_engine        : chr "stan"
#>  $ type              : chr "trend"
#>  $ series_names      : chr [1:3] "series_1" "series_2" "series_3"
#>  $ train_observations:List of 3
#>   ..$ series_1: int [1:75] 1 0 1 1 1 0 2 1 5 0 ...
#>   ..$ series_2: int [1:75] 0 0 1 1 1 0 3 5 5 1 ...
#>   ..$ series_3: int [1:75] 1 1 2 2 0 2 9 9 2 2 ...
#>  $ 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
#>   ..$ series_1: num [1:1000, 1:75] -0.139 0.589 -0.978 0.23 -0.11 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "trend[1,1]" "trend[2,1]" "trend[3,1]" "trend[4,1]" ...
#>   ..$ series_2: num [1:1000, 1:75] -0.0278 -0.499 0.1888 -0.9333 -0.0321 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "trend[1,2]" "trend[2,2]" "trend[3,2]" "trend[4,2]" ...
#>   ..$ series_3: num [1:1000, 1:75] 0.1364 -0.3046 -0.0367 -0.0314 -0.3705 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:75] "trend[1,3]" "trend[2,3]" "trend[3,3]" "trend[4,3]" ...
#>  $ forecasts         : NULL
#>  - attr(*, "class")= chr "mvgam_forecast"
plot(hc, series = 1)

plot(hc, series = 2)

plot(hc, series = 3)

# }

Extract hindcasts for a fitted `mvgam` object

Usage

Arguments

Value

Details

See also

Examples