Add RStan backend with multivariate observation model (Model B)

[Copilot]

Adds a Stan-based MCMC backend alongside the existing JAGS backend. Model B replaces K independent univariate normal likelihoods with a single K-variate normal per time point, capturing residual correlations across antigen-isotype pairs via an LKJ-prior covariance matrix. JAGS backend is fully unchanged.

New API

# Requires: install.packages("rstan")
fit <- run_mod_stan(
  data    = sim_data,   # same as_case_data() input as run_mod()
  model   = "model_b",
  chains  = 4,
  iter    = 2000,
  warmup  = 1000
)
# Returns sr_model tibble — identical columns to run_mod() output
# plot_predicted_curve(), post_summ(), etc. work without modification

attr(fit, "Omega_eps")  # K×K posterior-mean residual correlation matrix
attr(fit, "Sigma_eps")  # K×K posterior-mean residual covariance matrix

New files

File	Purpose
inst/stan/model_b.stan	Two-phase kinetics curve in functions{} block; K-variate normal likelihood; LKJ(2) prior on Omega_eps; half-Cauchy on tau_eps; natural-scale parameters in generated quantities
R/prep_data_stan.R	Converts prep_data() output → Stan list; strips NA padding (replaced with 0, masked by n_obs[i]) and dimension names
R/postprocess_stan_output.R	Extracts y0_nat/y1_nat/t1_nat/alpha_nat/shape_nat from generated quantities, pivots to long sr_model tibble
R/run_mod_stan.R	Orchestrates pipeline; calls prep_data(..., add_newperson=FALSE); compiles Stan model once per session; attaches Omega_eps, Sigma_eps, fitted_residuals, stan.fit attributes
vignettes/articles/stan-multivariate.qmd	Demonstrates Model B workflow; compares to JAGS backend; documents data flow
tests/testthat/test-run_mod_stan.R	Data-prep tests run unconditionally; sampling/output tests skip when rstan is absent

Modified files

• DESCRIPTION: rstan added to Suggests (optional — package installs without Stan); version bumped to 0.0.0.9048
• R/serodynamics-package.R: Stan availability note
• .Rbuildignore: Stan build artifact patterns (.o, .so, .d)
• NEWS.md, inst/WORDLIST: Updated for new symbols and terminology

Design notes

• rstan in Suggests, not Imports — avoids mandatory Stan/C++ toolchain for all users; run_mod_stan() errors informatively if rstan is absent
• No rstantools pre-compilation — Stan model compiled at first use per session via rstan::stan_model(); eliminates src/ Makevars infrastructure
• add_newperson = FALSE — Stan handles out-of-sample prediction via generated quantities, not a dummy row with NA observations

Original prompt

---

This section details on the original issue you should resolve

<issue_title>Add Stan backend with multivariate observation model (Model B)</issue_title>
<issue_description>## Summary

Add an RStan backend to serodynamics alongside the existing JAGS backend, starting with Model B: a multivariate observation model that captures residual correlations across antigen-isotype pairs within each time point.

Currently, run_mod() fits each biomarker independently via JAGS (model.jags). Model B replaces the K independent univariate normal likelihoods with a single K-variate normal likelihood per time point, introducing a residual covariance matrix Σ_ε. This is a pathogen-agnostic extension.

Motivation

• The Serodynamics primer (Hay et al., Epidemics 2024) identifies multi-biomarker kinetic modeling as a key future direction: "extensions to describe multiple biomarker kinetics simultaneously can improve TSI estimate precision"
• Ignoring within-individual residual correlations between IgG and IgA discards information and can bias downstream seroincidence estimates
• Stan offers LKJ priors for correlation matrices and non-centered parameterization, which are better suited for covariance estimation than JAGS's Wishart prior

Why Stan? (not extending JAGS)

1. JAGS cannot recognize Kronecker product conjugacy — confirmed on the JAGS forums. This blocks the Model C extension.
2. No LKJ prior in JAGS — the Wishart prior tangles scale and shape, especially problematic in small samples (Gelman 2006)
3. Non-centered parameterization — essential for hierarchical covariance models, not available in JAGS
4. C++ compilation is handled automatically — rstan::stan_model() compiles Stan → C++ transparently. Users just need rstan installed.

Package integration strategy

Use rstantools to properly integrate Stan into the R package:

# One-time setup (run by developer):
rstantools::use_rstan(pkgdir = ".", auto_config = TRUE)

This creates src/ directory with Makevars, inst/stan/ for .stan files, and configures the package to pre-compile Stan models at install time. Users then get fast model access without runtime compilation.

Key: rstan goes in Imports (not Suggests) if we use rstantools, because the compiled Stan models are part of the package binary. But we keep the JAGS backend fully functional — run_mod() is unchanged.

Scope of changes

New files

File	Purpose
inst/stan/model_b.stan	Stan model specification for Model B
R/run_mod_stan.R	Exported function: Stan analogue of run_mod()
R/prep_data_stan.R	Convert prepped_jags_data → Stan-ready list
R/postprocess_stan_output.R	Convert stanfit → sr_model class tibble
vignettes/articles/stan-multivariate.qmd	Vignette demonstrating Model B
tests/testthat/test-run_mod_stan.R	Unit tests

Files to modify

File	Change
DESCRIPTION	Add rstan, StanHeaders, rstantools to Imports; add SystemRequirements: GNU make
R/serodynamics-package.R	Document Stan availability
.Rbuildignore	Add Stan-related build artifacts

Files NOT modified

run_mod(), prep_data(), prep_priors(), calc_fit_mod(), model.jags — all stay exactly as is. Full backward compatibility.

API design: run_mod_stan()

Mirrors the existing run_mod() interface:

#' @title Run Stan Model (Multivariate)
#' @description Fits antibody kinetics model using RStan.
#'   Model B: multivariate observation model with correlated residuals.
#'   Model C: adds Kronecker product covariance for cross-biomarker
#'            parameter correlations.
#' @param data A data.frame from as_case_data()
#' @param model Character: "model_b" or "model_c"
#' @param chains Integer: number of MCMC chains (default 4)
#' @param iter Integer: total iterations per chain (default 2000)
#' @param warmup Integer: warmup iterations (default 1000)
#' @param adapt_delta Numeric: target acceptance rate (default 0.95)
#' @param max_treedepth Integer: max tree depth for NUTS (default 12)
#' @param strat Character: stratification variable (default NA)
#' @param ... Additional arguments passed to prep_priors()
#' @returns An sr_model class tibble (same format as run_mod())
#' @export
run_mod_stan <- function(data, model = c("model_b", "model_c"),
                          chains = 4, iter = 2000, warmup = 1000,
                          adapt_delta = 0.95, max_treedepth = 12,
                          strat = NA, ...) { ... }

Input: Same data from as_case_data() — no new data format needed.

Output: Same sr_model tibble with columns: Iteration, Chain, Parameter, Iso_type, Stratification, Subject, value. All existing plot/su...

• Fixes Add Stan backend with multivariate observation model (Model B) #190

---

⚡ Quickly spin up Copilot coding agent tasks from anywhere on your macOS or Windows machine with Raycast.