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(1111)
N <- 200
beta_temp <- mvgam:::sim_gp(rnorm(1), alpha_gp = 0.75, rho_gp = 10, h = N) + 0.5


## ----fig.alt = "Simulating time-varying effects in mvgam and R"------------------------------
plot(
  beta_temp,
  type = "l",
  lwd = 3,
  bty = "l",
  xlab = "Time",
  ylab = "Coefficient",
  col = "darkred"
)
box(bty = "l", lwd = 2)


## --------------------------------------------------------------------------------------------
temp <- rnorm(N, sd = 1)


## ----fig.alt = "Simulating time-varying effects in mvgam and R"------------------------------
out <- rnorm(N, mean = 4 + beta_temp * temp, sd = 0.25)
time <- seq_along(temp)
plot(
  out,
  type = "l",
  lwd = 3,
  bty = "l",
  xlab = "Time",
  ylab = "Outcome",
  col = "darkred"
)
box(bty = "l", lwd = 2)


## --------------------------------------------------------------------------------------------
data <- data.frame(out, temp, time)
data_train <- data[1:190, ]
data_test <- data[191:200, ]


## ----include=FALSE---------------------------------------------------------------------------
mod <- mvgam(
  out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 40),
  family = gaussian(),
  data = data_train,
  silent = 2
)


## ----eval=FALSE------------------------------------------------------------------------------
# mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 40),
#   family = gaussian(),
#   data = data_train,
#   silent = 2
# )

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


## --------------------------------------------------------------------------------------------
plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 190, lty = "dashed", lwd = 2)
lines(beta_temp, lwd = 2.5, col = "white")
lines(beta_temp, lwd = 2)


## --------------------------------------------------------------------------------------------
require(marginaleffects)
range_round <- function(x) {
  round(range(x, na.rm = TRUE), 2)
}
plot_predictions(
  mod,
  newdata = datagrid(
    time = unique,
    temp = range_round
  ),
  by = c("time", "temp", "temp"),
  type = "link"
)


## --------------------------------------------------------------------------------------------
fc <- forecast(mod, newdata = data_test)
plot(fc)


## ----include=FALSE---------------------------------------------------------------------------
mod <- mvgam(
  out ~ dynamic(temp, k = 40),
  family = gaussian(),
  data = data_train,
  silent = 2
)


## ----eval=FALSE------------------------------------------------------------------------------
# mod <- mvgam(out ~ dynamic(temp, k = 40),
#   family = gaussian(),
#   data = data_train,
#   silent = 2
# )

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


## --------------------------------------------------------------------------------------------
plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 190, lty = "dashed", lwd = 2)
lines(beta_temp, lwd = 2.5, col = "white")
lines(beta_temp, lwd = 2)


## --------------------------------------------------------------------------------------------
load(url("https://github.com/atsa-es/MARSS/raw/master/data/SalmonSurvCUI.rda"))
dplyr::glimpse(SalmonSurvCUI)


## --------------------------------------------------------------------------------------------
SalmonSurvCUI %>%
  # create a time variable
  dplyr::mutate(time = dplyr::row_number()) %>%
  # create a series variable
  dplyr::mutate(series = as.factor("salmon")) %>%
  # z-score the covariate CUI.apr
  dplyr::mutate(CUI.apr = as.vector(scale(CUI.apr))) %>%
  # convert logit-transformed survival back to proportional
  dplyr::mutate(survival = plogis(logit.s)) -> model_data


## --------------------------------------------------------------------------------------------
dplyr::glimpse(model_data)


## --------------------------------------------------------------------------------------------
plot_mvgam_series(data = model_data, y = "survival")


## ----include = FALSE-------------------------------------------------------------------------
mod0 <- mvgam(
  formula = survival ~ 1,
  trend_model = AR(),
  noncentred = TRUE,
  priors = prior(normal(-3.5, 0.5), class = Intercept),
  family = betar(),
  data = model_data,
  silent = 2
)


## ----eval = FALSE----------------------------------------------------------------------------
# mod0 <- mvgam(
#   formula = survival ~ 1,
#   trend_model = AR(),
#   noncentred = TRUE,
#   priors = prior(normal(-3.5, 0.5), class = Intercept),
#   family = betar(),
#   data = model_data,
#   silent = 2
# )

## --------------------------------------------------------------------------------------------
summary(mod0)


## --------------------------------------------------------------------------------------------
plot(mod0, type = "trend")


## ----include=FALSE---------------------------------------------------------------------------
mod1 <- mvgam(
  formula = survival ~ 1,
  trend_formula = ~ dynamic(CUI.apr, k = 25, scale = FALSE) - 1,
  trend_model = AR(),
  noncentred = TRUE,
  priors = prior(normal(-3.5, 0.5), class = Intercept),
  family = betar(),
  data = model_data,
  control = list(adapt_delta = 0.99),
  silent = 2
)


## ----eval=FALSE------------------------------------------------------------------------------
# mod1 <- mvgam(
#   formula = survival ~ 1,
#   trend_formula = ~ dynamic(CUI.apr, k = 25, scale = FALSE) - 1,
#   trend_model = AR(),
#   noncentred = TRUE,
#   priors = prior(normal(-3.5, 0.5), class = Intercept),
#   family = betar(),
#   data = model_data,
#   silent = 2
# )

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


## --------------------------------------------------------------------------------------------
plot(mod1, type = "trend")


## --------------------------------------------------------------------------------------------
plot(mod1, type = "forecast")


## --------------------------------------------------------------------------------------------
# Extract estimates of the process error 'sigma' for each model
mod0_sigma <- as.data.frame(mod0, variable = "sigma", regex = TRUE) %>%
  dplyr::mutate(model = "Mod0")
mod1_sigma <- as.data.frame(mod1, variable = "sigma", regex = TRUE) %>%
  dplyr::mutate(model = "Mod1")
sigmas <- rbind(mod0_sigma, mod1_sigma)

# Plot using ggplot2
require(ggplot2)
ggplot(sigmas, aes(y = `sigma[1]`, fill = model)) +
  geom_density(alpha = 0.3, colour = NA) +
  coord_flip()


## --------------------------------------------------------------------------------------------
plot(mod1, type = "smooths", trend_effects = TRUE)


## --------------------------------------------------------------------------------------------
loo_compare(mod0, mod1)


## ----include=FALSE---------------------------------------------------------------------------
lfo_mod0 <- lfo_cv(mod0, min_t = 30)
lfo_mod1 <- lfo_cv(mod1, min_t = 30)


## ----eval=FALSE------------------------------------------------------------------------------
# lfo_mod0 <- lfo_cv(mod0, min_t = 30)
# lfo_mod1 <- lfo_cv(mod1, min_t = 30)

## --------------------------------------------------------------------------------------------
sum(lfo_mod0$elpds)
sum(lfo_mod1$elpds)


## ----fig.alt = "Comparing forecast skill for dynamic beta regression models in mvgam and R"----
plot(
  x = 1:length(lfo_mod0$elpds) + 30,
  y = lfo_mod0$elpds - lfo_mod1$elpds,
  ylab = "ELPDmod0 - ELPDmod1",
  xlab = "Evaluation time point",
  pch = 16,
  col = "darkred",
  bty = "l"
)
abline(h = 0, lty = "dashed")