Extract hindcasts for a fitted mvgam object
Arguments
- object
listobject of classmvgamorjsdgam. Seemvgam()- ...
Ignored
- type
When this has the value
link(default) the linear predictor is calculated on the link scale. Ifexpectedis used, predictions reflect the expectation of the response (the mean) but ignore uncertainty in the observation process. Whenresponseis used, the predictions take uncertainty in the observation process into account to return predictions on the outcome scale. Whenvarianceis 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_Nwill return the estimated latent abundances from an N-mixture distribution, while typedetectionwill 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 hindcasts (i.e. retrodictions) are drawn from the fitted mvgam and
organized into a convenient format for plotting
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: 0x55c5ca5ddb48>
#> $ 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] 2 2 0 0 2 4 6 6 4 1 ...
#> ..$ series_2: int [1:75] 2 1 3 1 0 6 6 3 3 7 ...
#> ..$ series_3: int [1:75] 1 0 0 0 0 1 1 0 1 0 ...
#> $ 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 1 2 2 1 2 0 1 0 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] 4 0 3 2 1 1 2 0 4 1 ...
#> .. ..- 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 1 0 1 4 0 4 0 1 0 ...
#> .. ..- 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 1 0 4 2 response
#> 2 series_1 2 0 0 3 2 response
#> 3 series_1 3 0 0 2 0 response
#> 4 series_1 4 0 0 3 0 response
#> 5 series_1 5 1 0 4 2 response
#> 6 series_1 6 3 0 8 4 response
#> 7 series_1 7 4 0 10 6 response
#> 8 series_1 8 4 0 9 6 response
#> 9 series_1 9 3 0 8 4 response
#> 10 series_1 10 3 0 8 1 response
#> 11 series_1 11 2 0 7 3 response
#> 12 series_1 12 1 0 4 3 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 0.983 0.495 1.99 expected
#> 2 series_1 2 0.637 0.325 1.38 expected
#> 3 series_1 3 0.438 0.200 0.863 expected
#> 4 series_1 4 0.523 0.220 1.11 expected
#> 5 series_1 5 1.14 0.593 2.42 expected
#> 6 series_1 6 2.78 1.60 5.04 expected
#> 7 series_1 7 4.39 2.46 7.63 expected
#> 8 series_1 8 3.61 1.99 6.77 expected
#> 9 series_1 9 3.16 1.81 5.76 expected
#> 10 series_1 10 2.97 1.34 5.32 expected
#> 11 series_1 11 2.41 1.27 4.31 expected
#> 12 series_1 12 1.08 0.584 2.44 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.0738 -0.461 0.773 trend
#> 2 series_1 2 0.121 -0.434 0.915 trend
#> 3 series_1 3 -0.0539 -0.762 0.558 trend
#> 4 series_1 4 -0.0480 -0.798 0.557 trend
#> 5 series_1 5 0.0514 -0.485 0.756 trend
#> 6 series_1 6 0.0780 -0.407 0.727 trend
#> 7 series_1 7 0.127 -0.367 0.760 trend
#> 8 series_1 8 0.176 -0.308 0.898 trend
#> 9 series_1 9 0.0423 -0.507 0.660 trend
#> 10 series_1 10 -0.175 -0.893 0.320 trend
#> 11 series_1 11 0.0326 -0.535 0.644 trend
#> 12 series_1 12 0.145 -0.374 1.02 trend
plot(hc, series = 1)
plot(hc, series = 2)
plot(hc, series = 3)
# }