## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE
)

# library(MEMWAS)

## ----install-tarball, eval=FALSE----------------------------------------------
# install.packages(
#   "MEMWAS.tar.gz",
#   repos = NULL,
#   type = "source"
# )
# #library(MEMWAS)

## ----capabilities, eval=FALSE-------------------------------------------------
# MEMWAS_capabilities()

## ----pspline-fit, eval=FALSE--------------------------------------------------
# fit_smooth <- fit_MEMWAS(
#   y ~ x2,
#   data = sim_data,
#   id = "id",
#   time = "time",
#   random = ~ (1 | id),
#   autocor = "AR(1)",
#   smooth = list(
#     list(
#       name = "x1_curve",
#       type = "pspline",
#       variable = "x1",
#       k = 10L,
#       degree = 3L,
#       difference_order = 2L,
#       lambda = "auto"
#     )
#   ),
#   smooth_control = list(
#     optimizer = "grouped_cv",
#     cv_folds = 5L,
#     log_lambda_range = c(-4, 4),
#     metric = "RMSE"
#   )
# )
# summary(fit_smooth)$smooth_summary

## ----fit-interface, eval=FALSE------------------------------------------------
# fit <- fit_MEMWAS(
#   formula,
#   family = "gaussian",
#   data,
#   id,
#   time,
#   random = ~ (1 | id),
#   autocor = "AR(1)",
#   serial = NULL,
#   predictor_autocor = NULL,
#   L1_penalty = 0,
#   L2_penalty = 0,
#   smooth = NULL,
#   smooth_control = list(),
#   control = list(),
#   method = "ML",
#   random_cov = "unstructured",
#   approximation = "laplace",
#   init_approximation = "variational_inference",
#   quadrature_points = 7L,
#   se_method = "hessian",
#   link = NULL,
#   offset = NULL,
#   weights = NULL,
#   theta = NULL,
#   shape = NULL,
#   screen_nonlinear = FALSE,
#   nonlinear_spline = "restricted_cubic",
#   check_assumptions = FALSE,
#   bootstrap_inference = FALSE,
#   prediction_inference = FALSE,
#   verbose = TRUE
# )

## ----multi-serial, eval=FALSE-------------------------------------------------
# set.seed(1L)
# sim_data <- MEMWAS:::.simulate_panel_data(
#   n_id = 30L, n_time = 4L, beta = c(x1 = 0.6, x2 = -0.3, x3 = 0.2),
#   cor_matrix = diag(3L), intercept = 0.5,
#   sigma_eps = 0.5, sigma_b = 0.4, autocor = "NONE"
# )
# 
# serial_spec <- list(
#   outcome = serial_component(
#     structure = "OU",
#     name = "outcome_persistence"
#   ),
#   exposure = serial_component(
#     structure = "expOU",
#     predictor = "x1",
#     name = "exposure_persistence"
#   )
# )
# 
# fit_multi <- fit_MEMWAS(
#   y ~ x1 + x2 + x3,
#   family = "gaussian",
#   data = sim_data,
#   id = "id",
#   time = "time",
#   random = ~ (1 + x3 | id),
#   autocor = NULL,
#   serial = serial_spec
# )

## ----diagnostics, eval=FALSE--------------------------------------------------
# diagnose_approximation(fit)

## ----education-example, eval=FALSE--------------------------------------------
# set.seed(1L)
# sim_data <- MEMWAS:::.simulate_panel_data(
#   n_id = 30L, n_time = 4L,
#   beta = c(x1 = 0.7, x2 = -0.3),
#   cor_matrix = diag(2L), intercept = 1,
#   sigma_eps = 0.4, sigma_b = 0.5,
#   autocor = "AR(1)", autocor_param = list(rho = 0.45)
# )
# 
# fit_student <- fit_MEMWAS(
#   y ~ x1 + x2,
#   family = "gaussian",
#   data = sim_data,
#   id = "id",
#   time = "time",
#   random = ~ (1 | id),
#   autocor = "AR(1)",
#   method = "REML",
#   approximation = "laplace",
#   control = list(n_starts = 1L, cold_start_verification = FALSE),
#   verbose = FALSE
# )
# 
# summary(fit_student)
# diagnose_approximation(fit_student)

## ----healthcare-example, eval=FALSE-------------------------------------------
# set.seed(2L)
# sim_data <- MEMWAS:::.simulate_panel_data(
#   n_id = 30L, n_time = 4L,
#   beta = c(x1 = -0.5, x2 = -0.3, x3 = 0.2),
#   cor_matrix = diag(3L), intercept = 1.5,
#   sigma_eps = 0.5, sigma_b = 0.4, autocor = "NONE"
# )
# 
# fit_clinic <- fit_MEMWAS(
#   y ~ x1 + x2 + x3,
#   family = "gaussian",
#   data = sim_data,
#   id = "id",
#   time = "time",
#   random = ~ (1 + x1 | id),
#   random_cov = "unstructured",
#   autocor = "OU",
#   control = list(serial_time_scale = 1),
#   method = "REML",
#   verbose = FALSE
# )
# 
# summary(fit_clinic)

## ----nongaussian-example, eval=FALSE------------------------------------------
# set.seed(3L)
# sim_data <- MEMWAS:::.simulate_panel_data(
#   n_id = 30L, n_time = 3L,
#   beta = c(x1 = 0.15, x2 = -0.10),
#   cor_matrix = diag(2L), intercept = -0.2,
#   sigma_eps = 0, sigma_b = 0.35,
#   family = "poisson", autocor = "NONE"
# )
# 
# fit_count <- fit_MEMWAS(
#   y ~ x1 + x2,
#   family = "poisson",
#   data = sim_data,
#   id = "id",
#   time = "time",
#   random = ~ (1 | id),
#   autocor = "NONE",
#   approximation = "laplace",
#   init_approximation = "variational_inference",
#   se_method = "hessian",
#   verbose = FALSE
# )
# 
# summary(fit_count)
# diagnose_approximation(fit_count)
# 
# approx_comparison <- compare_approximations(
#   approximations = c("laplace", "variational_inference"),
#   formula = y ~ x1 + x2,
#   family = "poisson",
#   data = sim_data,
#   id = "id",
#   time = "time",
#   random = ~ (1 | id),
#   autocor = "NONE",
#   quadrature_points = 7L
# )
# 
# print(approx_comparison)

## ----predictor-serial-example, eval=FALSE-------------------------------------
# set.seed(4L)
# sim_data <- MEMWAS:::.simulate_panel_data(
#   n_id = 30L, n_time = 4L, beta = c(x1 = 0.6, x2 = -0.3),
#   cor_matrix = diag(2L), intercept = 0.5,
#   sigma_eps = 0.4, sigma_b = 0.3, autocor = "NONE"
# )
# 
# fit_exposure <- fit_MEMWAS(
#   y ~ x1 + x2,
#   family = "gaussian",
#   data = sim_data,
#   id = "id",
#   time = "time",
#   random = ~ (1 | id),
#   autocor = NULL,
#   serial = list(
#     outcome_persistence = serial_component(
#       structure = "OU",
#       name = "outcome_persistence"
#     ),
#     exposure_effect = serial_component(
#       structure = "expOU",
#       predictor = "x1",
#       name = "exposure_effect"
#     )
#   )
# )