Extract hindcasts for a fitted mvgam
object
Arguments
- object
list
object of classmvgam
orjsdgam
. Seemvgam()
- ...
Ignored
- type
When this has the value
link
(default) the linear predictor is calculated on the link scale. Ifexpected
is used, predictions reflect the expectation of the response (the mean) but ignore uncertainty in the observation process. Whenresponse
is used, the predictions take uncertainty in the observation process into account to return predictions on the outcome scale. Whenvariance
is used, the variance of the response with respect to the mean (mean-variance relationship) is returned. Whentype = "terms"
, each component of the linear predictor is returned separately in the form of alist
(possibly with standard errors, ifsummary = TRUE
): this includes parametric model components, followed by each smooth component, but excludes any offset and any intercept. Two special cases are also allowed: typelatent_N
will return the estimated latent abundances from an N-mixture distribution, while typedetection
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,
silent = 2)
# Hindcasts on response scale
hc <- hindcast(mod)
str(hc)
#> List of 15
#> $ call :Class 'formula' language y ~ s(season, bs = "cc")
#> .. ..- attr(*, ".Environment")=<environment: 0x5594e13b2d98>
#> $ 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
#> ..$ series_1: int [1:75] 1 2 1 0 0 0 1 0 2 2 ...
#> ..$ series_2: int [1:75] 3 0 3 0 3 2 0 4 2 1 ...
#> ..$ series_3: int [1:75] 3 0 3 1 1 1 0 2 1 1 ...
#> $ train_times :List of 3
#> ..$ series_1: int [1:75] 1 2 3 4 5 6 7 8 9 10 ...
#> ..$ series_2: int [1:75] 1 2 3 4 5 6 7 8 9 10 ...
#> ..$ series_3: 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 3 1 1 1 0 1 3 1 3 ...
#> .. ..- 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] 2 6 2 3 4 3 10 2 5 3 ...
#> .. ..- 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] 7 5 7 6 3 6 3 6 2 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"
head(summary(hc), 12)
#> # A tibble: 12 × 7
#> series time predQ50 predQ2.5 predQ97.5 truth type
#> <fct> <int> <dbl> <dbl> <dbl> <int> <chr>
#> 1 series_1 1 2 0 7 1 response
#> 2 series_1 2 1 0 5 2 response
#> 3 series_1 3 1 0 4 1 response
#> 4 series_1 4 0 0 2 0 response
#> 5 series_1 5 0 0 2 0 response
#> 6 series_1 6 0 0 3 0 response
#> 7 series_1 7 0 0 3 1 response
#> 8 series_1 8 1 0 4 0 response
#> 9 series_1 9 1 0 5 2 response
#> 10 series_1 10 2 0 6 2 response
#> 11 series_1 11 2 0 6 2 response
#> 12 series_1 12 1 0 5 0 response
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')
head(summary(hc), 12)
#> # A tibble: 12 × 6
#> series time predQ50 predQ2.5 predQ97.5 type
#> <fct> <int> <dbl> <dbl> <dbl> <chr>
#> 1 series_1 1 2.28 0.758 5.18 expected
#> 2 series_1 2 1.51 0.497 3.87 expected
#> 3 series_1 3 0.829 0.218 2.35 expected
#> 4 series_1 4 0.313 0.0500 1.20 expected
#> 5 series_1 5 0.272 0.0367 1.11 expected
#> 6 series_1 6 0.466 0.0942 1.70 expected
#> 7 series_1 7 0.697 0.152 2.25 expected
#> 8 series_1 8 0.797 0.179 2.63 expected
#> 9 series_1 9 1.27 0.356 3.65 expected
#> 10 series_1 10 1.75 0.570 4.45 expected
#> 11 series_1 11 1.84 0.624 4.66 expected
#> 12 series_1 12 1.35 0.252 3.64 expected
plot(hc, series = 1)
plot(hc, series = 2)
plot(hc, series = 3)
# Estimated latent trends
hc <- hindcast(mod, type = 'trend')
head(summary(hc), 12)
#> # A tibble: 12 × 6
#> series time predQ50 predQ2.5 predQ97.5 type
#> <fct> <int> <dbl> <dbl> <dbl> <chr>
#> 1 series_1 1 -0.591 -1.82 0.282 trend
#> 2 series_1 2 -0.179 -1.35 0.777 trend
#> 3 series_1 3 0.0563 -1.37 1.13 trend
#> 4 series_1 4 -0.423 -2.09 0.979 trend
#> 5 series_1 5 -0.557 -2.54 0.724 trend
#> 6 series_1 6 -0.440 -2.14 0.864 trend
#> 7 series_1 7 -0.338 -1.98 0.826 trend
#> 8 series_1 8 -0.418 -1.98 0.672 trend
#> 9 series_1 9 0.0877 -1.16 1.14 trend
#> 10 series_1 10 0.198 -0.924 1.13 trend
#> 11 series_1 11 -0.452 -1.71 0.505 trend
#> 12 series_1 12 -1.14 -2.87 -0.00818 trend
plot(hc, series = 1)
plot(hc, series = 2)
plot(hc, series = 3)
# }