## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  message = FALSE,
  warning = FALSE
)

## ----install, eval=FALSE------------------------------------------------------
# install.packages("remotes")
# remotes::install_github("qc-zhao/VCMoE")

## ----packages-----------------------------------------------------------------
library(VCMoE)
library(ggplot2)

## ----gaussian-simulate--------------------------------------------------------
set.seed(1)

sim <- simulate_vcmoe_gaussian(
  n = 180,
  k = 2,
  seed = 1,
  separation = 1.6,
  scenario = "well_separated"
)

head(sim$data)

## ----gaussian-fit-------------------------------------------------------------
fit <- vcmoe_fit(
  y ~ z1 | x1,
  data = sim$data,
  u = "u",
  family = "gaussian",
  k = 2,
  bandwidth = 0.35,
  u_grid = seq(0.15, 0.85, length.out = 4),
  control = list(maxit = 60, n_starts = 1, seed = 2, warn_ambiguous = FALSE)
)

fit

## ----gaussian-coefficients----------------------------------------------------
expert_coef <- coef(fit, "expert")
dim(expert_coef)
expert_coef[, , "z1"]

## ----gaussian-predictions-----------------------------------------------------
head(predict(fit, type = "mean"))
head(predict(fit, type = "posterior"))
head(predict(fit, type = "component"))

## ----gaussian-diagnostics-----------------------------------------------------
diagnostics <- vcmoe_diagnostics(fit)
diagnostics[, c("u", "converged", "ambiguous", "posterior_entropy", "effective_n")]

## ----gaussian-coefficient-plot------------------------------------------------
plot_coefficients(fit, "expert")

## ----gaussian-posterior-plot--------------------------------------------------
plot_posterior(fit)

## ----gaussian-scb-------------------------------------------------------------
band <- vcmoe_confband(
  fit,
  data = sim$data,
  level = 0.95,
  type = "simultaneous",
  coefficient_set = "expert",
  strict = FALSE
)

band
head(band$intervals[, c(
  "u", "component", "term", "block", "estimate",
  "lower", "upper", "status", "block_reason"
)])

## ----gaussian-scb-plot--------------------------------------------------------
scb_plot <- subset(
  band$intervals,
  coefficient_set == "expert" & block == "intercept" & status == "ok"
)

ggplot(scb_plot, aes(x = u, y = estimate, color = component, fill = component)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.18, color = NA) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~ term, scales = "free_y") +
  labs(
    x = "u",
    y = "coefficient",
    color = "component",
    fill = "component"
  ) +
  theme_minimal(base_size = 12)

## ----gaussian-bootstrap-------------------------------------------------------
boot <- vcmoe_bootstrap(
  fit,
  data = sim$data,
  B = 6,
  seed = 5,
  min_successful = 2,
  control = list(maxit = 40, n_starts = 1, warn_ambiguous = FALSE)
)

boot
head(confint(boot, parm = "expert", type = "simultaneous"))

## ----gaussian-bootstrap-plot--------------------------------------------------
plot_inference(
  boot,
  coefficient_set = "expert",
  type = "simultaneous",
  level = 0.95
)

## ----bandwidth, eval=FALSE----------------------------------------------------
# selection <- vcmoe_select_bandwidth(
#   y ~ z1 | x1,
#   data = sim$data,
#   u = "u",
#   family = "gaussian",
#   k = 2,
#   bandwidth_grid = c(0.25, 0.35, 0.45),
#   folds = 3,
#   u_grid = seq(0.15, 0.85, length.out = 4),
#   control = list(maxit = 60, n_starts = 1, seed = 3),
#   seed = 4
# )
# 
# selection
# selection$best_bandwidth