params <-
list(EVAL = TRUE)

## ---- echo = FALSE------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  eval = if (isTRUE(exists("params"))) params$EVAL else FALSE
)


## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE,   
  dpi = 100,
  fig.asp = 0.8,
  fig.width = 6,
  out.width = "60%",
  fig.align = "center")
library(mvgam)
library(ggplot2)
theme_set(theme_bw(base_size = 12, base_family = 'serif'))


## -----------------------------------------------------------------------------
set.seed(2345)
simdat <- sim_mvgam(T = 100, 
                    n_series = 3, 
                    trend_model = GP(),
                    prop_trend = 0.75,
                    family = poisson(),
                    prop_missing = 0.10)


## -----------------------------------------------------------------------------
str(simdat)


## ---- fig.alt = "Plotting time series features for GAM models in mvgam"-------
plot_mvgam_series(data = simdat$data_train, 
                  series = 'all')


## ---- fig.alt = "Plotting time series features for GAM models in mvgam"-------
plot_mvgam_series(data = simdat$data_train, 
                  newdata = simdat$data_test,
                  series = 1)


## ----include=FALSE------------------------------------------------------------
mod1 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
                s(time, by = series, bs = 'cr', k = 20),
              knots = list(season = c(0.5, 12.5)),
              trend_model = 'None',
              data = simdat$data_train)


## ----eval=FALSE---------------------------------------------------------------
## mod1 <- mvgam(y ~ s(season, bs = 'cc', k = 8) +
##                 s(time, by = series, bs = 'cr', k = 20),
##               knots = list(season = c(0.5, 12.5)),
##               trend_model = 'None',
##               data = simdat$data_train,
##               silent = 2)


## -----------------------------------------------------------------------------
summary(mod1, include_betas = FALSE)


## ---- fig.alt = "Plotting GAM smooth functions using mvgam"-------------------
conditional_effects(mod1, type = 'link')


## ----include=FALSE------------------------------------------------------------
mod2 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
                gp(time, by = series, c = 5/4, k = 20,
                   scale = FALSE),
              knots = list(season = c(0.5, 12.5)),
              trend_model = 'None',
              data = simdat$data_train,
              adapt_delta = 0.98,
              silent = 2)


## ----eval=FALSE---------------------------------------------------------------
## mod2 <- mvgam(y ~ s(season, bs = 'cc', k = 8) +
##                 gp(time, by = series, c = 5/4, k = 20,
##                    scale = FALSE),
##               knots = list(season = c(0.5, 12.5)),
##               trend_model = 'None',
##               data = simdat$data_train,
##               silent = 2)


## -----------------------------------------------------------------------------
summary(mod2, include_betas = FALSE)


## ---- fig.alt = "Summarising latent Gaussian Process parameters in mvgam"-----
mcmc_plot(mod2, variable = c('alpha_gp'), regex = TRUE, type = 'areas')


## ---- fig.alt = "Summarising latent Gaussian Process parameters in mvgam"-----
mcmc_plot(mod2, variable = c('rho_gp'), regex = TRUE, type = 'areas')


## ---- fig.alt = "Plotting latent Gaussian Process effects in mvgam and marginaleffects"----
conditional_effects(mod2, type = 'link')


## -----------------------------------------------------------------------------
fc_mod1 <- forecast(mod1, newdata = simdat$data_test)
fc_mod2 <- forecast(mod2, newdata = simdat$data_test)


## -----------------------------------------------------------------------------
str(fc_mod1)


## -----------------------------------------------------------------------------
plot(fc_mod1, series = 1)
plot(fc_mod2, series = 1)

plot(fc_mod1, series = 2)
plot(fc_mod2, series = 2)


## ----include=FALSE------------------------------------------------------------
mod2 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
                gp(time, by = series, c = 5/4, k = 20,
                   scale = FALSE),
              knots = list(season = c(0.5, 12.5)),
              trend_model = 'None',
              data = simdat$data_train,
              newdata = simdat$data_test,
              adapt_delta = 0.98,
              silent = 2)


## ----eval=FALSE---------------------------------------------------------------
## mod2 <- mvgam(y ~ s(season, bs = 'cc', k = 8) +
##                 gp(time, by = series, c = 5/4, k = 20,
##                    scale = FALSE),
##               knots = list(season = c(0.5, 12.5)),
##               trend_model = 'None',
##               data = simdat$data_train,
##               newdata = simdat$data_test,
##               silent = 2)


## -----------------------------------------------------------------------------
fc_mod2 <- forecast(mod2)


## ----warning=FALSE, fig.alt = "Plotting posterior forecast distributions using mvgam and R"----
plot(fc_mod2, series = 1)


## ----warning=FALSE------------------------------------------------------------
crps_mod1 <- score(fc_mod1, score = 'crps')
str(crps_mod1)
crps_mod1$series_1


## ----warning=FALSE------------------------------------------------------------
crps_mod1 <- score(fc_mod1, score = 'crps', interval_width = 0.6)
crps_mod1$series_1


## -----------------------------------------------------------------------------
link_mod1 <- forecast(mod1, newdata = simdat$data_test, type = 'link')
score(link_mod1, score = 'elpd')$series_1


## -----------------------------------------------------------------------------
energy_mod2 <- score(fc_mod2, score = 'energy')
str(energy_mod2)


## -----------------------------------------------------------------------------
energy_mod2$all_series


## -----------------------------------------------------------------------------
crps_mod1 <- score(fc_mod1, score = 'crps')
crps_mod2 <- score(fc_mod2, score = 'crps')

diff_scores <- crps_mod2$series_1$score -
  crps_mod1$series_1$score
plot(diff_scores, pch = 16, cex = 1.25, col = 'darkred', 
     ylim = c(-1*max(abs(diff_scores), na.rm = TRUE),
              max(abs(diff_scores), na.rm = TRUE)),
     bty = 'l',
     xlab = 'Forecast horizon',
     ylab = expression(CRPS[GP]~-~CRPS[spline]))
abline(h = 0, lty = 'dashed', lwd = 2)
gp_better <- length(which(diff_scores < 0))
title(main = paste0('GP better in ', gp_better, ' of 25 evaluations',
                    '\nMean difference = ', 
                    round(mean(diff_scores, na.rm = TRUE), 2)))


diff_scores <- crps_mod2$series_2$score -
  crps_mod1$series_2$score
plot(diff_scores, pch = 16, cex = 1.25, col = 'darkred', 
     ylim = c(-1*max(abs(diff_scores), na.rm = TRUE),
              max(abs(diff_scores), na.rm = TRUE)),
     bty = 'l',
     xlab = 'Forecast horizon',
     ylab = expression(CRPS[GP]~-~CRPS[spline]))
abline(h = 0, lty = 'dashed', lwd = 2)
gp_better <- length(which(diff_scores < 0))
title(main = paste0('GP better in ', gp_better, ' of 25 evaluations',
                    '\nMean difference = ', 
                    round(mean(diff_scores, na.rm = TRUE), 2)))

diff_scores <- crps_mod2$series_3$score -
  crps_mod1$series_3$score
plot(diff_scores, pch = 16, cex = 1.25, col = 'darkred', 
     ylim = c(-1*max(abs(diff_scores), na.rm = TRUE),
              max(abs(diff_scores), na.rm = TRUE)),
     bty = 'l',
     xlab = 'Forecast horizon',
     ylab = expression(CRPS[GP]~-~CRPS[spline]))
abline(h = 0, lty = 'dashed', lwd = 2)
gp_better <- length(which(diff_scores < 0))
title(main = paste0('GP better in ', gp_better, ' of 25 evaluations',
                    '\nMean difference = ', 
                    round(mean(diff_scores, na.rm = TRUE), 2)))