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"))


## -----------------------------------------------------------------------
simdat <- sim_mvgam(
  n_series = 4,
  T = 24,
  prop_missing = 0.2
)
head(simdat$data_train, 16)


## -----------------------------------------------------------------------
class(simdat$data_train$series)
levels(simdat$data_train$series)


## -----------------------------------------------------------------------
all(
  levels(simdat$data_train$series) %in%
    unique(simdat$data_train$series)
)


## -----------------------------------------------------------------------
summary(glm(
  y ~ series + time,
  data = simdat$data_train,
  family = poisson()
))


## -----------------------------------------------------------------------
summary(mgcv::gam(
  y ~ series + s(time, by = series),
  data = simdat$data_train,
  family = poisson()
))


## -----------------------------------------------------------------------
gauss_dat <- data.frame(
  outcome = rnorm(10),
  series = factor("series1", levels = "series1"),
  time = 1:10
)
gauss_dat


## -----------------------------------------------------------------------
mgcv::gam(outcome ~ time, family = betar(), data = gauss_dat)


## ----error=TRUE---------------------------------------------------------
try({
  mvgam(outcome ~ time, family = betar(), data = gauss_dat)
})


## -----------------------------------------------------------------------
# A function to ensure all timepoints within a sequence are identical
all_times_avail <- function(time, min_time, max_time) {
  identical(
    as.numeric(sort(time)),
    as.numeric(seq.int(from = min_time, to = max_time))
  )
}

# Get min and max times from the data
min_time <- min(simdat$data_train$time)
max_time <- max(simdat$data_train$time)

# Check that all times are recorded for each series
data.frame(
  series = simdat$data_train$series,
  time = simdat$data_train$time
) %>%
  dplyr::group_by(series) %>%
  dplyr::summarise(
    all_there = all_times_avail(
      time,
      min_time,
      max_time
    )
  ) -> checked_times
if (any(checked_times$all_there == FALSE)) {
  warning(
    "One or more series in is missing observations for one or more timepoints"
  )
} else {
  cat("All series have observations at all timepoints :)")
}


## -----------------------------------------------------------------------
bad_times <- data.frame(
  time = seq(1, 16, by = 2),
  series = factor("series_1"),
  outcome = rnorm(8)
)
bad_times


## ----error = TRUE-------------------------------------------------------
try({
  get_mvgam_priors(outcome ~ 1, data = bad_times, family = gaussian())
})


## -----------------------------------------------------------------------
bad_times %>%
  dplyr::right_join(expand.grid(
    time = seq(
      min(bad_times$time),
      max(bad_times$time)
    ),
    series = factor(unique(bad_times$series), levels = levels(bad_times$series))
  )) %>%
  dplyr::arrange(time) -> good_times
good_times


## ----error = TRUE-------------------------------------------------------
try({
  get_mvgam_priors(outcome ~ 1, data = good_times, family = gaussian())
})


## -----------------------------------------------------------------------
bad_levels <- data.frame(
  time = 1:8,
  series = factor(
    "series_1",
    levels = c(
      "series_1",
      "series_2"
    )
  ),
  outcome = rnorm(8)
)

levels(bad_levels$series)


## ----error = TRUE-------------------------------------------------------
try({
  get_mvgam_priors(outcome ~ 1, data = bad_levels, family = gaussian())
})


## -----------------------------------------------------------------------
setdiff(levels(bad_levels$series), unique(bad_levels$series))


## -----------------------------------------------------------------------
bad_levels %>%
  dplyr::mutate(series = droplevels(series)) -> good_levels
levels(good_levels$series)


## ----error = TRUE-------------------------------------------------------
try({
  get_mvgam_priors(
    outcome ~ 1,
    data = good_levels,
    family = gaussian()
  )
})


## -----------------------------------------------------------------------
miss_dat <- data.frame(
  outcome = rnorm(10),
  cov = c(NA, rnorm(9)),
  series = factor("series1", levels = "series1"),
  time = 1:10
)
miss_dat


## ----error = TRUE-------------------------------------------------------
try({
  get_mvgam_priors(
    outcome ~ cov,
    data = miss_dat,
    family = gaussian()
  )
})


## -----------------------------------------------------------------------
miss_dat <- list(
  outcome = rnorm(10),
  series = factor("series1", levels = "series1"),
  time = 1:10
)
miss_dat$cov <- matrix(rnorm(50), ncol = 5, nrow = 10)
miss_dat$cov[2, 3] <- NA


## ----error=TRUE---------------------------------------------------------
try({
  get_mvgam_priors(
    outcome ~ cov,
    data = miss_dat,
    family = gaussian()
  )
})


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


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


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


## -----------------------------------------------------------------------
data("all_neon_tick_data")
str(dplyr::ungroup(all_neon_tick_data))


## -----------------------------------------------------------------------
plotIDs <- c(
  "SCBI_013",
  "SCBI_002",
  "SERC_001",
  "SERC_005",
  "SERC_006",
  "SERC_012",
  "BLAN_012",
  "BLAN_005"
)


## -----------------------------------------------------------------------
model_dat <- all_neon_tick_data %>%
  dplyr::ungroup() %>%
  dplyr::mutate(target = ixodes_scapularis) %>%
  dplyr::filter(plotID %in% plotIDs) %>%
  dplyr::select(Year, epiWeek, plotID, target) %>%
  dplyr::mutate(epiWeek = as.numeric(epiWeek))


## -----------------------------------------------------------------------
model_dat %>%
  # Create all possible combos of plotID, Year and epiWeek;
  # missing outcomes will be filled in as NA
  dplyr::full_join(expand.grid(
    plotID = unique(model_dat$plotID),
    Year = unique(model_dat$Year),
    epiWeek = seq(1, 52)
  )) %>%
  # left_join back to original data so plotID and siteID will
  # match up, in case you need the siteID for anything else later on
  dplyr::left_join(
    all_neon_tick_data %>%
      dplyr::select(siteID, plotID) %>%
      dplyr::distinct()
  ) -> model_dat


## -----------------------------------------------------------------------
model_dat %>%
  dplyr::mutate(
    series = plotID,
    y = target
  ) %>%
  dplyr::mutate(
    siteID = factor(siteID),
    series = factor(series)
  ) %>%
  dplyr::select(-target, -plotID) %>%
  dplyr::arrange(Year, epiWeek, series) -> model_dat


## -----------------------------------------------------------------------
model_dat %>%
  dplyr::ungroup() %>%
  dplyr::group_by(series) %>%
  dplyr::arrange(Year, epiWeek) %>%
  dplyr::mutate(time = seq(1, dplyr::n())) %>%
  dplyr::ungroup() -> model_dat


## -----------------------------------------------------------------------
levels(model_dat$series)


## ----error=TRUE---------------------------------------------------------
try({
  get_mvgam_priors(
    y ~ 1,
    data = model_dat,
    family = poisson()
  )
})


## -----------------------------------------------------------------------
testmod <- mvgam(
  y ~
    s(epiWeek, by = series, bs = "cc") +
      s(series, bs = "re"),
  trend_model = AR(),
  data = model_dat,
  backend = "cmdstanr",
  run_model = FALSE
)


## -----------------------------------------------------------------------
str(testmod$model_data)


## -----------------------------------------------------------------------
stancode(testmod)