## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 10,
  fig.height = 6,
  dpi = 120,
  out.width = "100%",
  warning = FALSE,
  message = FALSE
)

## -----------------------------------------------------------------------------
library(VARcheck)
set.seed(35032)

## ----simulate-----------------------------------------------------------------
# Simulate a VAR(1) process
T  <- 150
p  <- 4
phi <- matrix(c(0.5, 0.1, 0.0, 0.0,
                0.2, 0.4, 0.0, 0.1,
                0.0, 0.0, 0.6, 0.0,
                0.1, 0.0, 0.2, 0.3), nrow = p, byrow = TRUE)

emp <- matrix(0, T, p)
for (t in 2:T) {
  emp[t, ] <- phi %*% emp[t - 1, ] + rnorm(p, sd = 0.5)
}

# Fit a simple lag-1 linear model per variable (stand-in for a real VAR fit)
pred <- matrix(NA, T, p)
res  <- matrix(NA, T, p)
for (j in seq_len(p)) {
  fit        <- lm(emp[-1, j] ~ emp[-T, j])
  pred[-1, j] <- predict(fit)
  res[-1, j]  <- residuals(fit)
}

# Simulate from the fitted model for the posterior predictive check
sim <- matrix(0, T, p)
for (t in 2:T) {
  sim[t, ] <- phi %*% sim[t - 1, ] + rnorm(p, sd = 0.5)
}

## ----vd-----------------------------------------------------------------------
vd <- new_var_data(
  empirical  = emp,
  predicted  = pred,
  residuals  = res,
  simulated  = sim,
  var_names  = c("Mood", "Energy", "Fatigue", "Anxiety")
)

vd

## ----basic-plot, fig.height = 9-----------------------------------------------
plot_var_check(vd)

## ----vars, fig.height = 5-----------------------------------------------------
plot_var_check(vd, vars = c("Mood", "Fatigue"))

## ----panels, fig.height = 9---------------------------------------------------
plot_var_check(vd, panels = c("data", "residuals"))

## ----multi-subject, fig.height = 9--------------------------------------------
emp2  <- emp + matrix(rnorm(T * p, sd = 0.2), T, p)
pred2 <- pred + matrix(rnorm(T * p, sd = 0.1), T, p)
res2  <- emp2 - pred2

vd_multi <- new_var_data(
  empirical = list(emp, emp2),
  predicted = list(pred, pred2),
  residuals = list(res, res2),
  simulated = list(sim, sim),
  var_names = c("Mood", "Energy", "Fatigue", "Anxiety")
)

plot_var_check(vd_multi, subject = 2)

## ----colours, fig.height = 9--------------------------------------------------
plot_var_check(vd, colors = list(predicted = "steelblue4"))

## ----theme, fig.height = 9----------------------------------------------------
plot_var_check(
  vd,
  theme = ggplot2::theme(
    text      = ggplot2::element_text(size = 9, family = "sans"),
    panel.grid.minor = ggplot2::element_blank()
  )
)

## ----ylim, fig.height = 9-----------------------------------------------------
plot_var_check(vd, ylim_data = c(-3, 3), ylim_res = c(-2, 2))

## ----mlvar-fit, eval = FALSE--------------------------------------------------
# library(mlVAR)
# 
# vars <- c("A", "B", "C", "D")
# 
# mlVAR_out <- mlVAR(data, vars = c("A", "B", "C", "D"),
#                    idvar = "id", lags = 1,
#                    dayvar = "day", beepvar = "beep")
# 
# pred_df <- predict(mlVAR_out)
# res_df  <- residuals(mlVAR_out)
# sim_out <- mlVAR:::resimulate(mlVAR_out, keep_missing = TRUE,
#                               variance = "empirical")

## ----mlvar-single, eval = FALSE-----------------------------------------------
# i <- 1  # subject index
# unique_ids <- unique(data$id)
# 
# check_df <- new_var_data(
#   empirical = as.matrix(data[data$id == unique_ids[i], vars]),
#   predicted = as.matrix(pred_df[pred_df$id == unique_ids[i], vars]),
#   residuals = as.matrix(res_df[res_df$id == unique_ids[i], vars]),
#   simulated = as.matrix(sim_out[sim_out$id == unique_ids[i], vars])
# )
# 
# plot_var_check(check_df)

## ----mlvar-multi, eval = FALSE------------------------------------------------
# n_subj   <- length(unique_ids)
# 
# check_df_all <- new_var_data(
#   empirical = lapply(unique_ids, \(id) as.matrix(data[data$id == id, vars])),
#   predicted = lapply(unique_ids, \(id) as.matrix(pred_df[pred_df$id == id, vars])),
#   residuals = lapply(unique_ids, \(id) as.matrix(res_df[res_df$id == id, vars])),
#   simulated = lapply(unique_ids, \(id) as.matrix(sim_out[sim_out$id == id, vars]))
# )
# 
# plot_var_check(check_df_all, subject = 3)