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

# library(MEMWAS)

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

## ----install-github, eval=FALSE-----------------------------------------------
# remotes::install_github("EnochKang/MEMWAS")
# library(MEMWAS)

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

## ----fit-interface, eval=FALSE------------------------------------------------
# fit <- fit_MEMWAS(
#   formula,
#   family = "gaussian",
#   data,
#   id,
#   time,
#   random = ~ 1,
#   autocor = "AR(1)",
#   serial = NULL,
#   predictor_autocor = NULL,
#   autocor_predictors = NULL,
#   L1_penalty = 0,
#   L2_penalty = 0,
#   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,
#   check_assumptions = FALSE,
#   bootstrap_inference = FALSE,
#   prediction_inference = FALSE,
#   verbose = TRUE
# )

## ----multi-serial, eval=FALSE-------------------------------------------------
# serial_spec <- list(
#   outcome = serial_component(
#     structure = "OU",
#     name = "outcome_persistence"
#   ),
#   exposure = serial_component(
#     structure = "expOU",
#     predictor = "exposure",
#     name = "exposure_persistence"
#   )
# )
# 
# fit_multi <- fit_MEMWAS(
#   outcome ~ time_c + exposure + treatment,
#   family = "gaussian",
#   data = dat,
#   id = "subject_id",
#   time = "visit_day",
#   random = ~(1 + time_c | subject_id),
#   autocor = NULL,
#   serial = serial_spec
# )

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

## ----education-example, eval=FALSE--------------------------------------------
# set.seed(2026)
# 
# m <- 80L
# n_time <- 4L
# rho <- 0.45
# sigma_u <- 3
# sigma_eps <- 2
# sigma_b <- 5
# 
# student_dat <- data.frame(
#   student_id = rep(seq_len(m), each = n_time),
#   semester = rep(seq_len(n_time), times = m)
# )
# 
# student_dat$study_minutes <- rnorm(m * n_time, 90, 25)
# student_dat$attendance_rate <- runif(m * n_time, 0.70, 1.00)
# b0 <- rnorm(m, 0, sigma_b)
# 
# simulate_ar1 <- function(n, rho, sigma) {
#   u <- numeric(n)
#   u[1L] <- rnorm(1L, sd = sigma)
#   innovation_sd <- sigma * sqrt(1 - rho^2)
#   if (n > 1L) {
#     for (tt in 2L:n) {
#       u[tt] <- rho * u[tt - 1L] + rnorm(1L, sd = innovation_sd)
#     }
#   }
#   u
# }
# 
# u <- unlist(
#   lapply(seq_len(m), function(i) simulate_ar1(n_time, rho, sigma_u)),
#   use.names = FALSE
# )
# 
# eps <- rnorm(m * n_time, 0, sigma_eps)
# 
# student_dat$reading_score <-
#   55 +
#   0.08 * student_dat$study_minutes +
#   18 * student_dat$attendance_rate +
#   b0[student_dat$student_id] +
#   u + eps
# 
# fit_student <- fit_MEMWAS(
#   reading_score ~ study_minutes + attendance_rate,
#   family = "gaussian",
#   data = student_dat,
#   id = "student_id",
#   time = "semester",
#   random = ~ 1,
#   autocor = "AR(1)",
#   method = "REML",
#   approximation = "laplace",
#   verbose = FALSE
# )
# 
# summary(fit_student)
# diagnose_approximation(fit_student)

## ----healthcare-example, eval=FALSE-------------------------------------------
# set.seed(2026)
# 
# m <- 120L
# visit_times <- c(0, 2, 5, 9, 14, 20)
# n_time <- length(visit_times)
# 
# clinic_dat <- data.frame(
#   patient_id = rep(seq_len(m), each = n_time),
#   visit_week = rep(visit_times, times = m)
# )
# 
# clinic_dat$visit_week_c <-
#   clinic_dat$visit_week - mean(visit_times)
# clinic_dat$therapy_sessions <- rpois(m * n_time, lambda = 2)
# clinic_dat$sleep_hours <- rnorm(m * n_time, mean = 6.5, sd = 1.1)
# 
# # Illustrative response only; a realistic simulation would also generate
# # the stated random-slope and OU latent processes.
# clinic_dat$depression_score <-
#   14 -
#   0.12 * clinic_dat$visit_week_c -
#   0.30 * clinic_dat$therapy_sessions -
#   0.60 * clinic_dat$sleep_hours +
#   rnorm(m * n_time, sd = 4)
# 
# fit_clinic <- fit_MEMWAS(
#   depression_score ~ visit_week_c + therapy_sessions + sleep_hours,
#   family = "gaussian",
#   data = clinic_dat,
#   id = "patient_id",
#   time = "visit_week",
#   random = ~(1 + visit_week_c | patient_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(2026)
# 
# m <- 100L
# n_time <- 4L
# count_dat <- data.frame(
#   patient_id = rep(seq_len(m), each = n_time),
#   visit = rep(seq_len(n_time), times = m)
# )
# 
# count_dat$visit_c <- count_dat$visit - mean(seq_len(n_time))
# count_dat$therapy_sessions <- rpois(m * n_time, 2)
# count_dat$sleep_hours <- rnorm(m * n_time, 6.5, 1.1)
# b0 <- rnorm(m, 0, 0.45)
# 
# eta <-
#   -0.4 -
#   0.08 * count_dat$visit_c +
#   0.10 * count_dat$therapy_sessions -
#   0.06 * count_dat$sleep_hours +
#   b0[count_dat$patient_id]
# 
# count_dat$emergency_visits <- rnbinom(
#   nrow(count_dat),
#   mu = exp(eta),
#   size = 2.5
# )
# 
# fit_count <- fit_MEMWAS(
#   emergency_visits ~ visit_c + therapy_sessions + sleep_hours,
#   family = "negative_binomial",
#   data = count_dat,
#   id = "patient_id",
#   time = "visit",
#   random = ~ 1,
#   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",
#     "adaptive_gaussian_quadrature",
#     "variational_inference",
#     "pql"
#   ),
#   formula = emergency_visits ~ visit_c + therapy_sessions + sleep_hours,
#   family = "negative_binomial",
#   data = count_dat,
#   id = "patient_id",
#   time = "visit",
#   random = ~ 1,
#   autocor = "NONE",
#   quadrature_points = 7L
# )
# 
# print(approx_comparison)

## ----predictor-serial-example, eval=FALSE-------------------------------------
# fit_exposure <- fit_MEMWAS(
#   outcome ~ time_c + exposure + treatment,
#   family = "gaussian",
#   data = dat,
#   id = "subject_id",
#   time = "visit_day",
#   random = ~(1 | subject_id),
#   autocor = NULL,
#   serial = list(
#     outcome_persistence = serial_component(
#       structure = "OU",
#       name = "outcome_persistence"
#     ),
#     exposure_effect = serial_component(
#       structure = "expOU",
#       predictor = "exposure",
#       name = "exposure_effect"
#     )
#   )
# )

