From bbe1badc3e73282d43e646a99eb249c85e49f324 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 24 Mar 2026 19:46:36 +0000 Subject: [PATCH 1/2] Initial plan From 914b3dff45d858d6812a095f4b064ad25c94f85f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 24 Mar 2026 20:26:46 +0000 Subject: [PATCH 2/2] Add RStan backend with multivariate observation model (Model B) Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Agent-Logs-Url: https://github.com/ucdavis/serodynamics/sessions/5e1fea32-1feb-43af-968d-c0a34b711ac3 --- .Rbuildignore | 5 + DESCRIPTION | 3 +- NAMESPACE | 3 + NEWS.md | 10 + R/postprocess_stan_output.R | 78 ++++++++ R/prep_data_stan.R | 74 ++++++++ R/run_mod_stan.R | 205 +++++++++++++++++++++ R/serodynamics-package.R | 13 ++ inst/WORDLIST | 49 +++-- inst/stan/model_b.stan | 137 ++++++++++++++ man/nepal_sees_jags_output.Rd | 2 +- man/plot_predicted_curve.Rd | 2 +- man/postprocess_stan_output.Rd | 34 ++++ man/prep_data_stan.Rd | 69 +++++++ man/run_mod.Rd | 2 +- man/run_mod_stan.Rd | 113 ++++++++++++ man/use_att_names.Rd | 2 +- tests/testthat/_snaps/as_case_data.md | 6 +- tests/testthat/test-run_mod_stan.R | 223 +++++++++++++++++++++++ vignettes/articles/stan-multivariate.qmd | 186 +++++++++++++++++++ 20 files changed, 1192 insertions(+), 24 deletions(-) create mode 100644 R/postprocess_stan_output.R create mode 100644 R/prep_data_stan.R create mode 100644 R/run_mod_stan.R create mode 100644 inst/stan/model_b.stan create mode 100644 man/postprocess_stan_output.Rd create mode 100644 man/prep_data_stan.Rd create mode 100644 man/run_mod_stan.Rd create mode 100644 tests/testthat/test-run_mod_stan.R create mode 100644 vignettes/articles/stan-multivariate.qmd diff --git a/.Rbuildignore b/.Rbuildignore index eaddfcf5..2d83188d 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -29,3 +29,8 @@ ^_quarto\.yml$ ^\.quarto$ ^\.devcontainer$ +^inst/stan/.*\.o$ +^inst/stan/.*\.so$ +^inst/stan/.*\.d$ +^inst/stan/.*_rng\.hpp$ +^inst/stan/.*\.rds$ diff --git a/DESCRIPTION b/DESCRIPTION index ea5d92de..4de7e546 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: serodynamics Title: What the Package Does (One Line, Title Case) -Version: 0.0.0.9047 +Version: 0.0.0.9048 Authors@R: c( person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"), comment = "Author of the method and original code."), @@ -37,6 +37,7 @@ Suggests: knitr, readr, rlist, + rstan, sessioninfo, spelling, stringr, diff --git a/NAMESPACE b/NAMESPACE index 2db4d027..67fda474 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,9 +15,12 @@ export(plot_jags_trace) export(plot_predicted_curve) export(post_summ) export(postprocess_jags_output) +export(postprocess_stan_output) export(prep_data) +export(prep_data_stan) export(prep_priors) export(run_mod) +export(run_mod_stan) export(serodynamics_example) export(sim_case_data) export(sim_n_obs) diff --git a/NEWS.md b/NEWS.md index 89199522..519d9a04 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,15 @@ # serodynamics (development version) +* Added RStan backend with multivariate observation model (Model B): + `run_mod_stan()`, `prep_data_stan()`, and `postprocess_stan_output()`. + Model B replaces K independent univariate normal likelihoods with a + single K-variate normal likelihood per time point, capturing residual + correlations across antigen-isotype pairs. Output is an `sr_model` + tibble identical in structure to `run_mod()`, so all existing plot and + summary functions work without modification. Requires the optional + dependency `rstan` (`install.packages("rstan")`). The JAGS backend is + fully backward-compatible and unchanged. + * Added dev container configuration for persistent, cached development environment that includes R, JAGS, and all dependencies preinstalled, making Copilot Workspace sessions much faster. diff --git a/R/postprocess_stan_output.R b/R/postprocess_stan_output.R new file mode 100644 index 00000000..690f0dfd --- /dev/null +++ b/R/postprocess_stan_output.R @@ -0,0 +1,78 @@ +#' Postprocess Stan model output +#' +#' @description +#' Converts a `stanfit` object returned by [rstan::sampling()] into an +#' `sr_model` tibble with the same column layout as [run_mod()]. +#' +#' Extracted parameters are the natural-scale quantities computed in the +#' Stan `generated quantities` block: +#' `y0_nat`, `y1_nat`, `t1_nat`, `alpha_nat`, `shape_nat`. +#' +#' @param stan_fit A `stanfit` object from [rstan::sampling()]. +#' @param ids Character vector of subject IDs (length `N`). +#' @param antigen_isos Character vector of antigen-isotype names (length `K`). +#' @param stratification Character scalar: stratification label used to +#' populate the `Stratification` column. Defaults to `"None"`. +#' +#' @returns A [tibble::tbl_df] with columns: +#' `Iteration`, `Chain`, `Parameter`, `Iso_type`, `Stratification`, +#' `Subject`, `value`. +#' +#' @seealso [run_mod_stan()] +#' @export +postprocess_stan_output <- function( + stan_fit, + ids, + antigen_isos, + stratification = "None") { + + # Extract natural-scale parameters from generated quantities + param_names <- c("y0_nat", "y1_nat", "t1_nat", "alpha_nat", "shape_nat") + param_labels <- c("y0", "y1", "t1", "alpha", "shape") + + N <- length(ids) + K <- length(antigen_isos) + + # rstan::extract with permuted = FALSE returns [iter, chains, parameters] + all_samples <- rstan::extract(stan_fit, pars = param_names, permuted = FALSE) + + n_iter <- dim(all_samples)[1] + n_chains <- dim(all_samples)[2] + + # Build long-format tibble for all subjects, biomarkers, and parameters + rows <- vector("list", N * K * length(param_labels)) + idx <- 1L + + for (k in seq_len(K)) { + for (subj in seq_len(N)) { + for (p in seq_along(param_labels)) { + # Parameter name in Stan: e.g. "y0_nat[subj,k]" + stan_par_name <- paste0(param_names[p], "[", subj, ",", k, "]") + col_pos <- which(dimnames(all_samples)[[3]] == stan_par_name) + + if (length(col_pos) == 0L) { + cli::cli_abort( + "Stan parameter {.val {stan_par_name}} not found in stanfit object." + ) + } + + # Flatten across chains: Iteration index increases across chains + for (chain in seq_len(n_chains)) { + iter_vals <- all_samples[, chain, col_pos] + rows[[idx]] <- tibble::tibble( + Iteration = seq_len(n_iter), + Chain = chain, + Parameter = param_labels[p], + Iso_type = antigen_isos[k], + Stratification = stratification, + Subject = ids[subj], + value = as.numeric(iter_vals) + ) + idx <- idx + 1L + } + } + } + } + + dplyr::bind_rows(rows) +} diff --git a/R/prep_data_stan.R b/R/prep_data_stan.R new file mode 100644 index 00000000..ee004672 --- /dev/null +++ b/R/prep_data_stan.R @@ -0,0 +1,74 @@ +#' Prepare data for Stan Model B +#' +#' @description +#' Converts the output of [prep_data()] into a named list suitable for +#' passing to `rstan::sampling()`. Missing values in the padded arrays are +#' replaced with `0`; Stan only accesses entries up to `n_obs[i]` per +#' subject, so the padding values are never used in the likelihood. +#' +#' The `mu_hyp` and `sigma_hyp` arguments correspond to the hyperprior +#' parameters from [prep_priors()]: +#' - `mu_hyp` maps to `mu_hyp_param` (prior means on the five log-scale +#' curve parameters) +#' - `sigma_hyp` maps to `1 / sqrt(prec_hyp_param)` (prior standard +#' deviations derived from the JAGS precision hyperparameter) +#' +#' @param prepped_jags_data A `prepped_jags_data` list returned by +#' [prep_data()] (with `add_newperson = FALSE`). +#' @param mu_hyp A numeric vector of length 5 giving the prior means for +#' the five log-scale curve parameters +#' (`log_y0`, `log_delta`, `log_t1`, `log_alpha`, `log_shape_minus_1`). +#' Defaults match the `mu_hyp_param` argument of [prep_priors()]. +#' @param sigma_hyp A numeric vector of length 5 giving the prior standard +#' deviations for the five log-scale curve parameters. +#' Defaults are derived from the `prec_hyp_param` argument of +#' [prep_priors()] as `1 / sqrt(prec_hyp_param)`. +#' +#' @returns A named [list] suitable for passing to `rstan::sampling()`. +#' Contains the following elements: +#' - `N` — number of subjects +#' - `K` — number of antigen-isotype pairs +#' - `max_obs` — maximum number of observations per subject +#' - `n_obs` — integer vector of length `N` giving actual observations +#' - `time_obs` — `N × max_obs` matrix of observation times +#' (padded with `0`) +#' - `log_y_obs` — `N × max_obs × K` array of log antibody levels +#' (padded with `0`) +#' - `mu_hyp` — vector of length 5 (hyperprior means) +#' - `sigma_hyp` — vector of length 5 (hyperprior SDs) +#' +#' @seealso [prep_data()], [run_mod_stan()] +#' @export +#' +#' @examples +#' set.seed(1) +#' raw_data <- +#' serocalculator::typhoid_curves_nostrat_100 |> +#' sim_case_data(n = 5) +#' prepped <- prep_data(raw_data, add_newperson = FALSE) +#' stan_data <- prep_data_stan(prepped) +#' str(stan_data) +prep_data_stan <- function( + prepped_jags_data, + mu_hyp = c(1.0, 7.0, 1.0, -4.0, -1.0), + sigma_hyp = 1.0 / sqrt(c(1.0, 0.00001, 1.0, 0.001, 1.0))) { + + smpl_t <- prepped_jags_data$smpl.t + logy <- prepped_jags_data$logy + nsmpl <- prepped_jags_data$nsmpl + + # Replace NA padding values with 0; Stan only reads positions <= n_obs[i] + smpl_t[is.na(smpl_t)] <- 0.0 + logy[is.na(logy)] <- 0.0 + + list( + N = unname(prepped_jags_data$nsubj), + K = unname(prepped_jags_data$n_antigen_isos), + max_obs = unname(ncol(smpl_t)), + n_obs = as.integer(nsmpl), + time_obs = unname(smpl_t), + log_y_obs = unname(logy), + mu_hyp = mu_hyp, + sigma_hyp = sigma_hyp + ) +} diff --git a/R/run_mod_stan.R b/R/run_mod_stan.R new file mode 100644 index 00000000..f27e72bd --- /dev/null +++ b/R/run_mod_stan.R @@ -0,0 +1,205 @@ +#' Run Stan Model (Multivariate Observation) +#' +#' @description +#' `run_mod_stan()` fits the antibody kinetics model using RStan with a +#' multivariate (correlated) observation model. +#' +#' **Model B** replaces the K independent univariate normal likelihoods used +#' by [run_mod()] with a single K-variate normal likelihood per time point, +#' introducing a residual covariance matrix **Sigma_eps** that captures +#' correlations across antigen-isotype pairs within each observation. +#' +#' The two-phase kinetics curve (rise then decay) is identical to the JAGS +#' model in [run_mod()]; only the likelihood and the prior on the residual +#' structure differ. +#' +#' @section Stan model file: +#' The Stan model is stored in +#' `system.file("stan/model_b.stan", package = "serodynamics")`. +#' It is compiled once per R session and cached by `rstan`. +#' +#' @section Required package: +#' `rstan` must be installed to use this function. +#' It is listed as a suggested (optional) dependency of `serodynamics`. +#' Install it with `install.packages("rstan")`. +#' +#' @param data A `case_data` object from [as_case_data()]. +#' @param model Character scalar: currently only `"model_b"` is supported. +#' @param chains Integer: number of MCMC chains (default `4`). +#' @param iter Integer: total iterations per chain including warmup +#' (default `2000`). +#' @param warmup Integer: number of warmup iterations per chain +#' (default `1000`). +#' @param adapt_delta Numeric: target Metropolis acceptance rate for the +#' NUTS sampler (default `0.95`). +#' @param max_treedepth Integer: maximum tree depth for the NUTS sampler +#' (default `12`). +#' @param strat Character scalar: name of a stratification variable in +#' `data`, or `NA` for an unstratified analysis (default `NA`). +#' @param seed Integer or `NULL`: random seed passed to +#' `rstan::sampling()` for reproducibility. +#' @param ... Additional arguments passed to [prep_data_stan()]. +#' +#' @returns An `sr_model` class object: a subclass of [dplyr::tbl_df] with +#' the same column structure as the output of [run_mod()]: +#' `Iteration`, `Chain`, `Parameter`, `Iso_type`, `Stratification`, +#' `Subject`, `value`. +#' +#' Additional [attributes]: +#' - `class`: includes `"sr_model"`. +#' - `chains`: number of chains used. +#' - `iter`: total iterations per chain. +#' - `warmup`: warmup iterations per chain. +#' - `adapt_delta`: target acceptance rate. +#' - `Omega_eps`: K × K residual correlation matrix (posterior mean). +#' - `Sigma_eps`: K × K residual covariance matrix (posterior mean). +#' - `fitted_residuals`: data frame of fitted and residual values. +#' - `stan.fit`: the raw `stanfit` object. +#' +#' @seealso [run_mod()] for the JAGS backend, [prep_data_stan()], +#' [postprocess_stan_output()] +#' @export +#' +#' @examples +#' \dontrun{ +#' # Requires rstan to be installed +#' set.seed(1) +#' sim_data <- serocalculator::typhoid_curves_nostrat_100 |> +#' sim_case_data(n = 5) +#' +#' fit <- run_mod_stan( +#' data = sim_data, +#' model = "model_b", +#' chains = 1, +#' iter = 200, +#' warmup = 100 +#' ) +#' print(fit) +#' } +run_mod_stan <- function( + data, + model = c("model_b"), + chains = 4L, + iter = 2000L, + warmup = 1000L, + adapt_delta = 0.95, + max_treedepth = 12L, + strat = NA, + seed = NULL, + ...) { + + if (!requireNamespace("rstan", quietly = TRUE)) { + cli::cli_abort( + c( + "Package {.pkg rstan} is required to use {.fn run_mod_stan}.", + "i" = "Install it with {.code install.packages(\"rstan\")}." + ) + ) + } + + model <- match.arg(model) + + stan_file <- system.file("stan", paste0(model, ".stan"), + package = "serodynamics") + if (!nzchar(stan_file)) { + cli::cli_abort( + c("Stan model file for {.val {model}} not found in the", + "package installation.") + ) + } + + # Compile Stan model (cached after first call within a session) + stan_mod <- rstan::stan_model(file = stan_file) + + ## Build stratification list + if (is.na(strat)) { + strat_list <- "None" + } else { + strat_list <- unique(data[[strat]]) + } + + ## Shell to accumulate results + all_results <- list() + final_stan_fit <- NULL + + for (i in strat_list) { + ## Subset data for this stratum + if (is.na(strat)) { + dl_sub <- data + } else { + dl_sub <- data |> + dplyr::filter(.data[[strat]] == i) + } + + ## Prepare data: do NOT add newperson (Stan handles predictions differently) + prepped <- prep_data(dl_sub, add_newperson = FALSE) + stan_data <- prep_data_stan(prepped, ...) + + ids <- attr(prepped, "ids") + antigen_isos <- attr(prepped, "antigens") + + ## Run NUTS sampler + stan_fit <- rstan::sampling( + object = stan_mod, + data = stan_data, + chains = as.integer(chains), + iter = as.integer(iter), + warmup = as.integer(warmup), + control = list( + adapt_delta = adapt_delta, + max_treedepth = as.integer(max_treedepth) + ), + seed = if (is.null(seed)) { + sample.int(.Machine$integer.max, 1L) + } else { + seed + }, + refresh = 0L # suppress per-iteration progress output + ) + final_stan_fit <- stan_fit + + ## Convert to long-format tibble + stratum_results <- postprocess_stan_output( + stan_fit = stan_fit, + ids = ids, + antigen_isos = antigen_isos, + stratification = as.character(i) + ) + + all_results[[as.character(i)]] <- stratum_results + } + + ## Combine strata + out <- dplyr::bind_rows(all_results) + out <- tibble::as_tibble(out) + + ## Compute fitted values and residuals using last stratum's data + fit_res <- calc_fit_mod(modeled_dat = out, original_data = dl_sub) + + ## Extract posterior mean of residual correlation/covariance matrices + omega_samples <- rstan::extract(final_stan_fit, pars = "Omega_eps")$Omega_eps + sigma_samples <- rstan::extract(final_stan_fit, pars = "Sigma_eps")$Sigma_eps + omega_mean <- apply(omega_samples, c(2, 3), mean) + sigma_mean <- apply(sigma_samples, c(2, 3), mean) + + colnames(omega_mean) <- antigen_isos + rownames(omega_mean) <- antigen_isos + colnames(sigma_mean) <- antigen_isos + rownames(sigma_mean) <- antigen_isos + + ## Attach attributes to mirror run_mod() output + out <- out |> + structure( + class = union("sr_model", class(out)), + chains = chains, + iter = iter, + warmup = warmup, + adapt_delta = adapt_delta, + Omega_eps = omega_mean, + Sigma_eps = sigma_mean, + fitted_residuals = fit_res, + stan.fit = final_stan_fit + ) + + out +} diff --git a/R/serodynamics-package.R b/R/serodynamics-package.R index 974b3970..4d588a5e 100644 --- a/R/serodynamics-package.R +++ b/R/serodynamics-package.R @@ -28,3 +28,16 @@ #' @importFrom utils read.csv ## usethis namespace: end NULL + +# Stan availability note ------------------------------------------------------- +# +# The Stan-based functions `run_mod_stan()`, `prep_data_stan()`, and +# `postprocess_stan_output()` require the `rstan` package, which is listed as +# a *suggested* (optional) dependency of serodynamics. These functions call +# `requireNamespace("rstan")` at runtime and emit an informative error when +# rstan is not installed. +# +# To use the Stan backend, install rstan from CRAN. +# +# The Stan model file is stored in inst/stan/model_b.stan and is compiled on +# first use within each R session by `rstan::stan_model()`. diff --git a/inst/WORDLIST b/inst/WORDLIST index 8ea3768d..2b2b24d3 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,50 +1,61 @@ AVX -CMD -CodeFactor -Codecov -HlyE -IgA -IgG -Iso -JAGS's -JAGs -Lifecycle -MSVC -Postprocess -Rhat -TSI -Vectorization -Wasserstein behaviour +bigl +bigr biomarker biomarkers +boldsymbol +CMD +Codecov +CodeFactor colour darwin dev +diag dobson dotplots +eps +frac +Gelman geoms ggmcmc ggproto glibc +HlyE hyp hyperprior hyperpriors +IgA +IgG iso +Iso isotype isotypes +JAGs +JAGS's +le libm +Lifecycle linux +LKJ +mathbf +mathcal +mathrm mcmc mitre +MSVC newperson nmc onwards params pkgdown +Postprocess prec preinstalling rhat +Rhat +RStan +SDs seroconversion seroresponse seroresponses @@ -52,5 +63,11 @@ strat stratifications tbl testthat +tibble +TSI unstratified +varepsilon +Vectorization +Visualise +Wasserstein wishdf diff --git a/inst/stan/model_b.stan b/inst/stan/model_b.stan new file mode 100644 index 00000000..9d23602a --- /dev/null +++ b/inst/stan/model_b.stan @@ -0,0 +1,137 @@ +// Model B: Multivariate observation model for antibody kinetics +// +// This Stan model implements the same two-phase antibody kinetics curve +// as model.jags but replaces K independent univariate normal likelihoods +// with a single K-variate normal likelihood per time point, capturing +// residual correlations across antigen-isotype pairs (Model B). +// +// References: +// Teunis et al. (2016). Biomarker dynamics: estimating infection rates +// from serological data. Statistics in Medicine, 35(22), 3956-3967. +// Hay et al. (2024). Epidemics, 100806. + +functions { + // Two-phase antibody kinetics curve on the log scale. + // + // Parameters (on the natural scale): + // t - time since infection + // t1 - time to peak (> 0) + // y0 - baseline concentration (> 0) + // y1 - peak concentration (> y0) + // alpha - decay rate (> 0) + // rho - shape parameter (> 1) + real two_phase_log_y(real t, real t1, real y0, real y1, + real alpha, real rho) { + real beta = (log(y1) - log(y0)) / t1; + if (t <= t1) { + return log(y0) + beta * t; + } else { + // Decay phase: see Teunis et al. (2016) eq. 15 + real arg = pow(y1, 1.0 - rho) - (1.0 - rho) * alpha * (t - t1); + return log(arg) / (1.0 - rho); + } + } +} + +data { + int N; // number of subjects + int K; // number of antigen-isotype pairs + int max_obs; // max observations per subject + array[N] int n_obs; // observations per subject + array[N, max_obs] real time_obs; // observation times (padded with 0) + array[N, max_obs, K] real log_y_obs; // log antibody levels (padded with 0) + + // Hyperprior values (correspond to JAGS defaults) + // Means for the 5 log-scale parameters per biomarker: + // [log_y0, log_delta, log_t1, log_alpha, log_shape_minus_1] + vector[5] mu_hyp; // prior means for mu_par + vector[5] sigma_hyp; // prior SDs for mu_par +} + +parameters { + // Population-level means per biomarker (same for all K biomarkers) + array[K] vector[5] mu_par; + + // Population-level standard deviations per biomarker per parameter + array[K] vector[5] sigma_par; + + // Subject × biomarker raw (log-scale) parameters + // par[i, k] = [log_y0, log_delta, log_t1, log_alpha, log_shape_minus_1] + array[N, K] vector[5] par; + + // Residual covariance: Cholesky factor of correlation matrix (K × K) + cholesky_factor_corr[K] L_Omega_eps; + + // Marginal standard deviations for residuals (one per biomarker) + vector[K] tau_eps; +} + +transformed parameters { + // Cholesky factor of the residual covariance matrix + matrix[K, K] L_Sigma_eps = diag_pre_multiply(tau_eps, L_Omega_eps); +} + +model { + // --- Priors on residual covariance structure --- + // LKJ prior on correlation matrix (eta=2 slightly regularizes toward identity) + L_Omega_eps ~ lkj_corr_cholesky(2.0); + // Half-Cauchy prior on residual scales + tau_eps ~ cauchy(0.0, 2.5); + + // --- Priors on population-level means --- + for (k in 1:K) { + mu_par[k] ~ normal(mu_hyp, sigma_hyp); + sigma_par[k] ~ exponential(1.0); + } + + // --- Hierarchical subject-level parameters and likelihood --- + for (i in 1:N) { + for (k in 1:K) { + // Independent normal priors on each of the 5 log-scale parameters + par[i, k] ~ normal(mu_par[k], sigma_par[k]); + } + + for (o in 1:n_obs[i]) { + real t = time_obs[i, o]; + vector[K] mu_it; + + for (k in 1:K) { + // Transform log-scale parameters to natural scale + real y0 = exp(par[i, k][1]); + real y1 = y0 + exp(par[i, k][2]); // peak = baseline + positive delta + real t1 = exp(par[i, k][3]); + real alpha = exp(par[i, k][4]); + real rho = exp(par[i, k][5]) + 1.0; // shape > 1 + + mu_it[k] = two_phase_log_y(t, t1, y0, y1, alpha, rho); + } + + // K-variate normal likelihood (Model B: correlated residuals) + to_vector(log_y_obs[i, o]) ~ multi_normal_cholesky(mu_it, L_Sigma_eps); + } + } +} + +generated quantities { + // Recover the full correlation and covariance matrices for diagnostics + matrix[K, K] Omega_eps = multiply_lower_tri_self_transpose(L_Omega_eps); + matrix[K, K] Sigma_eps = quad_form_diag(Omega_eps, tau_eps); + + // Natural-scale parameter summaries for each subject × biomarker + // Used by postprocess_stan_output() to build the sr_model tibble + array[N, K] real y0_nat; + array[N, K] real y1_nat; + array[N, K] real t1_nat; + array[N, K] real alpha_nat; + array[N, K] real shape_nat; + + for (i in 1:N) { + for (k in 1:K) { + y0_nat[i, k] = exp(par[i, k][1]); + y1_nat[i, k] = y0_nat[i, k] + exp(par[i, k][2]); + t1_nat[i, k] = exp(par[i, k][3]); + alpha_nat[i, k] = exp(par[i, k][4]); + shape_nat[i, k] = exp(par[i, k][5]) + 1.0; + } + } +} diff --git a/man/nepal_sees_jags_output.Rd b/man/nepal_sees_jags_output.Rd index 2bd13915..4ea195d1 100644 --- a/man/nepal_sees_jags_output.Rd +++ b/man/nepal_sees_jags_output.Rd @@ -5,7 +5,7 @@ \alias{nepal_sees_jags_output} \title{SEES Typhoid run_mod jags output} \format{ -An S3 object of class \code{sr_model}: A \link[dplyr:tbl_df]{dplyr::tbl_df} that contains the +An S3 object of class \code{sr_model}: A \link[dplyr:defunct]{dplyr::tbl_df} that contains the posterior predictive distribution of the person-specific parameters for a "new person" with no observed data (\code{Subject = "newperson"}) and posterior distributions of the person-specific parameters for two arbitrarily-chosen diff --git a/man/plot_predicted_curve.Rd b/man/plot_predicted_curve.Rd index 7e2fabde..b5638b88 100644 --- a/man/plot_predicted_curve.Rd +++ b/man/plot_predicted_curve.Rd @@ -31,7 +31,7 @@ samples from the posterior distribution of the model parameters.} \item{antigen_iso}{The antigen isotype to plot; for example, "HlyE_IgA" or "HlyE_IgG".} -\item{dataset}{(Optional) A \link[dplyr:tbl_df]{dplyr::tbl_df} with observed antibody response +\item{dataset}{(Optional) A \link[dplyr:defunct]{dplyr::tbl_df} with observed antibody response data. Must contain: \itemize{ diff --git a/man/postprocess_stan_output.Rd b/man/postprocess_stan_output.Rd new file mode 100644 index 00000000..4274c4ca --- /dev/null +++ b/man/postprocess_stan_output.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/postprocess_stan_output.R +\name{postprocess_stan_output} +\alias{postprocess_stan_output} +\title{Postprocess Stan model output} +\usage{ +postprocess_stan_output(stan_fit, ids, antigen_isos, stratification = "None") +} +\arguments{ +\item{stan_fit}{A \code{stanfit} object from \code{\link[rstan:sampling]{rstan::sampling()}}.} + +\item{ids}{Character vector of subject IDs (length \code{N}).} + +\item{antigen_isos}{Character vector of antigen-isotype names (length \code{K}).} + +\item{stratification}{Character scalar: stratification label used to +populate the \code{Stratification} column. Defaults to \code{"None"}.} +} +\value{ +A \link[tibble:tbl_df-class]{tibble::tbl_df} with columns: +\code{Iteration}, \code{Chain}, \code{Parameter}, \code{Iso_type}, \code{Stratification}, +\code{Subject}, \code{value}. +} +\description{ +Converts a \code{stanfit} object returned by \code{\link[rstan:sampling]{rstan::sampling()}} into an +\code{sr_model} tibble with the same column layout as \code{\link[=run_mod]{run_mod()}}. + +Extracted parameters are the natural-scale quantities computed in the +Stan \verb{generated quantities} block: +\code{y0_nat}, \code{y1_nat}, \code{t1_nat}, \code{alpha_nat}, \code{shape_nat}. +} +\seealso{ +\code{\link[=run_mod_stan]{run_mod_stan()}} +} diff --git a/man/prep_data_stan.Rd b/man/prep_data_stan.Rd new file mode 100644 index 00000000..e5cfa398 --- /dev/null +++ b/man/prep_data_stan.Rd @@ -0,0 +1,69 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prep_data_stan.R +\name{prep_data_stan} +\alias{prep_data_stan} +\title{Prepare data for Stan Model B} +\usage{ +prep_data_stan( + prepped_jags_data, + mu_hyp = c(1, 7, 1, -4, -1), + sigma_hyp = 1/sqrt(c(1, 1e-05, 1, 0.001, 1)) +) +} +\arguments{ +\item{prepped_jags_data}{A \code{prepped_jags_data} list returned by +\code{\link[=prep_data]{prep_data()}} (with \code{add_newperson = FALSE}).} + +\item{mu_hyp}{A numeric vector of length 5 giving the prior means for +the five log-scale curve parameters +(\code{log_y0}, \code{log_delta}, \code{log_t1}, \code{log_alpha}, \code{log_shape_minus_1}). +Defaults match the \code{mu_hyp_param} argument of \code{\link[=prep_priors]{prep_priors()}}.} + +\item{sigma_hyp}{A numeric vector of length 5 giving the prior standard +deviations for the five log-scale curve parameters. +Defaults are derived from the \code{prec_hyp_param} argument of +\code{\link[=prep_priors]{prep_priors()}} as \code{1 / sqrt(prec_hyp_param)}.} +} +\value{ +A named \link{list} suitable for passing to \code{rstan::sampling()}. +Contains the following elements: +\itemize{ +\item \code{N} — number of subjects +\item \code{K} — number of antigen-isotype pairs +\item \code{max_obs} — maximum number of observations per subject +\item \code{n_obs} — integer vector of length \code{N} giving actual observations +\item \code{time_obs} — \verb{N × max_obs} matrix of observation times +(padded with \code{0}) +\item \code{log_y_obs} — \verb{N × max_obs × K} array of log antibody levels +(padded with \code{0}) +\item \code{mu_hyp} — vector of length 5 (hyperprior means) +\item \code{sigma_hyp} — vector of length 5 (hyperprior SDs) +} +} +\description{ +Converts the output of \code{\link[=prep_data]{prep_data()}} into a named list suitable for +passing to \code{rstan::sampling()}. Missing values in the padded arrays are +replaced with \code{0}; Stan only accesses entries up to \code{n_obs[i]} per +subject, so the padding values are never used in the likelihood. + +The \code{mu_hyp} and \code{sigma_hyp} arguments correspond to the hyperprior +parameters from \code{\link[=prep_priors]{prep_priors()}}: +\itemize{ +\item \code{mu_hyp} maps to \code{mu_hyp_param} (prior means on the five log-scale +curve parameters) +\item \code{sigma_hyp} maps to \code{1 / sqrt(prec_hyp_param)} (prior standard +deviations derived from the JAGS precision hyperparameter) +} +} +\examples{ +set.seed(1) +raw_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 5) +prepped <- prep_data(raw_data, add_newperson = FALSE) +stan_data <- prep_data_stan(prepped) +str(stan_data) +} +\seealso{ +\code{\link[=prep_data]{prep_data()}}, \code{\link[=run_mod_stan]{run_mod_stan()}} +} diff --git a/man/run_mod.Rd b/man/run_mod.Rd index 08f91427..585210ae 100644 --- a/man/run_mod.Rd +++ b/man/run_mod.Rd @@ -93,7 +93,7 @@ If specified, must be 2 values long: }} } \value{ -An \code{sr_model} class object: a subclass of \link[dplyr:tbl_df]{dplyr::tbl_df} that +An \code{sr_model} class object: a subclass of \link[dplyr:defunct]{dplyr::tbl_df} that contains MCMC samples from the joint posterior distribution of the model parameters, conditional on the provided input \code{data}, including the following: diff --git a/man/run_mod_stan.Rd b/man/run_mod_stan.Rd new file mode 100644 index 00000000..f1f06f59 --- /dev/null +++ b/man/run_mod_stan.Rd @@ -0,0 +1,113 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_mod_stan.R +\name{run_mod_stan} +\alias{run_mod_stan} +\title{Run Stan Model (Multivariate Observation)} +\usage{ +run_mod_stan( + data, + model = c("model_b"), + chains = 4L, + iter = 2000L, + warmup = 1000L, + adapt_delta = 0.95, + max_treedepth = 12L, + strat = NA, + seed = NULL, + ... +) +} +\arguments{ +\item{data}{A \code{case_data} object from \code{\link[=as_case_data]{as_case_data()}}.} + +\item{model}{Character scalar: currently only \code{"model_b"} is supported.} + +\item{chains}{Integer: number of MCMC chains (default \code{4}).} + +\item{iter}{Integer: total iterations per chain including warmup +(default \code{2000}).} + +\item{warmup}{Integer: number of warmup iterations per chain +(default \code{1000}).} + +\item{adapt_delta}{Numeric: target Metropolis acceptance rate for the +NUTS sampler (default \code{0.95}).} + +\item{max_treedepth}{Integer: maximum tree depth for the NUTS sampler +(default \code{12}).} + +\item{strat}{Character scalar: name of a stratification variable in +\code{data}, or \code{NA} for an unstratified analysis (default \code{NA}).} + +\item{seed}{Integer or \code{NULL}: random seed passed to +\code{rstan::sampling()} for reproducibility.} + +\item{...}{Additional arguments passed to \code{\link[=prep_data_stan]{prep_data_stan()}}.} +} +\value{ +An \code{sr_model} class object: a subclass of \link[dplyr:defunct]{dplyr::tbl_df} with +the same column structure as the output of \code{\link[=run_mod]{run_mod()}}: +\code{Iteration}, \code{Chain}, \code{Parameter}, \code{Iso_type}, \code{Stratification}, +\code{Subject}, \code{value}. + +Additional \link{attributes}: +\itemize{ +\item \code{class}: includes \code{"sr_model"}. +\item \code{chains}: number of chains used. +\item \code{iter}: total iterations per chain. +\item \code{warmup}: warmup iterations per chain. +\item \code{adapt_delta}: target acceptance rate. +\item \code{Omega_eps}: K × K residual correlation matrix (posterior mean). +\item \code{Sigma_eps}: K × K residual covariance matrix (posterior mean). +\item \code{fitted_residuals}: data frame of fitted and residual values. +\item \code{stan.fit}: the raw \code{stanfit} object. +} +} +\description{ +\code{run_mod_stan()} fits the antibody kinetics model using RStan with a +multivariate (correlated) observation model. + +\strong{Model B} replaces the K independent univariate normal likelihoods used +by \code{\link[=run_mod]{run_mod()}} with a single K-variate normal likelihood per time point, +introducing a residual covariance matrix \strong{Sigma_eps} that captures +correlations across antigen-isotype pairs within each observation. + +The two-phase kinetics curve (rise then decay) is identical to the JAGS +model in \code{\link[=run_mod]{run_mod()}}; only the likelihood and the prior on the residual +structure differ. +} +\section{Stan model file}{ + +The Stan model is stored in +\code{system.file("stan/model_b.stan", package = "serodynamics")}. +It is compiled once per R session and cached by \code{rstan}. +} + +\section{Required package}{ + +\code{rstan} must be installed to use this function. +It is listed as a suggested (optional) dependency of \code{serodynamics}. +Install it with \code{install.packages("rstan")}. +} + +\examples{ +\dontrun{ + # Requires rstan to be installed + set.seed(1) + sim_data <- serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 5) + + fit <- run_mod_stan( + data = sim_data, + model = "model_b", + chains = 1, + iter = 200, + warmup = 100 + ) + print(fit) +} +} +\seealso{ +\code{\link[=run_mod]{run_mod()}} for the JAGS backend, \code{\link[=prep_data_stan]{prep_data_stan()}}, +\code{\link[=postprocess_stan_output]{postprocess_stan_output()}} +} diff --git a/man/use_att_names.Rd b/man/use_att_names.Rd index f43fcef0..9a2151e0 100644 --- a/man/use_att_names.Rd +++ b/man/use_att_names.Rd @@ -16,7 +16,7 @@ The input \link{data.frame} with columns named after attributes. \description{ \code{use_att_names} takes prepared longitudinal data for antibody kinetic modeling and names columns using attribute values to allow merging -with a modeled \code{\link[=run_mod]{run_mod()}} output \link[dplyr:tbl_df]{dplyr::tbl_df}. The column names include +with a modeled \code{\link[=run_mod]{run_mod()}} output \link[dplyr:defunct]{dplyr::tbl_df}. The column names include \code{Subject}, \code{Iso_type}, \code{t}, and \code{result}. } \keyword{internal} diff --git a/tests/testthat/_snaps/as_case_data.md b/tests/testthat/_snaps/as_case_data.md index 44130b20..efd6396b 100644 --- a/tests/testthat/_snaps/as_case_data.md +++ b/tests/testthat/_snaps/as_case_data.md @@ -1,9 +1,9 @@ # results are consistent with simulated data list(names = c("id", "visit_num", "timeindays", "iter", "antigen_iso", - "y0", "y1", "t1", "alpha", "r", "value"), class = c("case_data", - "tbl_df", "tbl", "data.frame"), id_var = "id", biomarker_var = "antigen_iso", - timeindays = "timeindays", value_var = "value") + "y0", "y1", "t1", "alpha", "r", "value"), id_var = "id", biomarker_var = "antigen_iso", + timeindays = "timeindays", value_var = "value", class = c("case_data", + "tbl_df", "tbl", "data.frame")) # results are consistent with SEES data diff --git a/tests/testthat/test-run_mod_stan.R b/tests/testthat/test-run_mod_stan.R new file mode 100644 index 00000000..27bba05b --- /dev/null +++ b/tests/testthat/test-run_mod_stan.R @@ -0,0 +1,223 @@ +test_that("prep_data_stan() returns correctly structured Stan data list", { + withr::local_seed(1) + raw_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 5) + + prepped <- prep_data(raw_data, add_newperson = FALSE) + stan_data <- prep_data_stan(prepped) + + expect_type(stan_data, "list") + expect_named( + stan_data, + c("N", "K", "max_obs", "n_obs", "time_obs", "log_y_obs", + "mu_hyp", "sigma_hyp") + ) + expect_equal(stan_data$N, prepped$nsubj) + expect_equal(stan_data$K, prepped$n_antigen_isos) + expect_equal(stan_data$max_obs, unname(ncol(prepped$smpl.t))) + expect_equal(stan_data$n_obs, as.integer(prepped$nsmpl)) + # No NA values allowed in Stan data arrays + expect_false(anyNA(stan_data$time_obs)) + expect_false(anyNA(stan_data$log_y_obs)) + + # Dimensions + expect_equal( + unname(dim(stan_data$time_obs)), + c(stan_data$N, stan_data$max_obs) + ) + expect_equal( + unname(dim(stan_data$log_y_obs)), + c(stan_data$N, stan_data$max_obs, stan_data$K) + ) + + # Hyperprior defaults + expect_length(stan_data$mu_hyp, 5L) + expect_length(stan_data$sigma_hyp, 5L) + expect_true(all(stan_data$sigma_hyp > 0)) +}) + +test_that("prep_data_stan() accepts custom hyperpriors", { + withr::local_seed(2) + raw_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 3) + + prepped <- prep_data(raw_data, add_newperson = FALSE) + custom_mu <- c(0, 5, 0, -3, -2) + custom_sigma <- c(2, 100, 2, 50, 2) + stan_data <- prep_data_stan(prepped, + mu_hyp = custom_mu, + sigma_hyp = custom_sigma) + + expect_equal(stan_data$mu_hyp, custom_mu) + expect_equal(stan_data$sigma_hyp, custom_sigma) +}) + +test_that("run_mod_stan() gives informative error when rstan is absent", { + skip_if( + requireNamespace("rstan", quietly = TRUE), + "rstan is installed — skipping 'rstan absent' test" + ) + withr::local_seed(1) + raw_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 3) + + expect_error( + run_mod_stan(raw_data, chains = 1, iter = 10, warmup = 5), + regexp = "rstan" + ) +}) + +# The remaining tests require rstan to be installed. +# They are skipped in environments where rstan is unavailable (e.g., CI +# runners without Stan). Run them locally after `install.packages("rstan")`. + +test_that("run_mod_stan() returns an sr_model tibble", { + skip_if_not_installed("rstan") + withr::local_seed(42) + sim_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 5) + + fit <- run_mod_stan( + data = sim_data, + model = "model_b", + chains = 1L, + iter = 200L, + warmup = 100L, + seed = 1L + ) |> + suppressWarnings() + + expect_s3_class(fit, "sr_model") + expect_s3_class(fit, "tbl_df") +}) + +test_that("run_mod_stan() output has the expected columns", { + skip_if_not_installed("rstan") + withr::local_seed(42) + sim_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 5) + + fit <- run_mod_stan( + data = sim_data, + model = "model_b", + chains = 1L, + iter = 200L, + warmup = 100L, + seed = 1L + ) |> + suppressWarnings() + + expected_cols <- c("Iteration", "Chain", "Parameter", "Iso_type", + "Stratification", "Subject", "value") + expect_equal(names(fit), expected_cols) +}) + +test_that("run_mod_stan() output contains expected parameters", { + skip_if_not_installed("rstan") + withr::local_seed(42) + sim_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 5) + + fit <- run_mod_stan( + data = sim_data, + model = "model_b", + chains = 1L, + iter = 200L, + warmup = 100L, + seed = 1L + ) |> + suppressWarnings() + + expect_setequal( + unique(fit$Parameter), + c("y0", "y1", "t1", "alpha", "shape") + ) +}) + +test_that("run_mod_stan() attaches Omega_eps and Sigma_eps attributes", { + skip_if_not_installed("rstan") + withr::local_seed(42) + sim_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 5) + + fit <- run_mod_stan( + data = sim_data, + model = "model_b", + chains = 1L, + iter = 200L, + warmup = 100L, + seed = 1L + ) |> + suppressWarnings() + + omega <- attr(fit, "Omega_eps") + sigma <- attr(fit, "Sigma_eps") + + # Both should be square matrices + expect_true(is.matrix(omega)) + expect_true(is.matrix(sigma)) + expect_equal(nrow(omega), ncol(omega)) + expect_equal(nrow(sigma), ncol(sigma)) + + # Correlation matrix: diagonal entries should be 1 + expect_equal(diag(omega), rep(1.0, nrow(omega)), tolerance = 1e-6) +}) + +test_that("run_mod_stan() attaches fitted_residuals attribute", { + skip_if_not_installed("rstan") + withr::local_seed(42) + sim_data <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 5) + + fit <- run_mod_stan( + data = sim_data, + model = "model_b", + chains = 1L, + iter = 200L, + warmup = 100L, + seed = 1L + ) |> + suppressWarnings() + + fit_res <- attr(fit, "fitted_residuals") + expect_s3_class(fit_res, "data.frame") + expect_true(all(c("Subject", "Iso_type", "t", "fitted", "residual") %in% + names(fit_res))) +}) + +test_that("run_mod_stan() works with stratification", { + skip_if_not_installed("rstan") + withr::local_seed(1) + strat1 <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 3) |> + dplyr::mutate(grp = "A") + withr::local_seed(2) + strat2 <- + serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 3) |> + dplyr::mutate(grp = "B") + dataset <- dplyr::bind_rows(strat1, strat2) + + fit <- run_mod_stan( + data = dataset, + model = "model_b", + chains = 1L, + iter = 200L, + warmup = 100L, + strat = "grp", + seed = 1L + ) |> + suppressWarnings() + + expect_s3_class(fit, "sr_model") + expect_setequal(unique(fit$Stratification), c("A", "B")) +}) diff --git a/vignettes/articles/stan-multivariate.qmd b/vignettes/articles/stan-multivariate.qmd new file mode 100644 index 00000000..8a796c19 --- /dev/null +++ b/vignettes/articles/stan-multivariate.qmd @@ -0,0 +1,186 @@ +--- +title: "Stan Backend: Multivariate Observation Model (Model B)" +format: html +--- + +```{r} +#| include: false +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 5, + eval = FALSE +) +``` + +## Overview + +Model B extends the JAGS antibody kinetics model by replacing the +*K* independent univariate normal likelihoods (one per antigen-isotype) with a +single *K*-variate normal likelihood per time point. This captures residual +correlations across antigen-isotype pairs within each observation, which the +independent model ignores. + +The `run_mod_stan()` function implements Model B using RStan. Its output is an +`sr_model` tibble with the same column layout as `run_mod()`, so all existing +plotting and summary functions (`plot_predicted_curve()`, `post_summ()`, etc.) +work without modification. + +### Why Stan instead of JAGS? + +- **LKJ prior** for correlation matrices — better than the Wishart prior in + small samples (Gelman 2006). +- **NUTS sampler** — more efficient exploration of correlated posteriors. +- **Extendable** to Model C (Kronecker product covariance for cross-biomarker + parameter correlations), which JAGS cannot handle. + +### The two-phase kinetics curve + +The mathematical model is identical to `model.jags`: + +| Phase | Condition | Equation | +|-------|-----------|----------| +| Rising | $t \le t_1$ | $\log y(t) = \log y_0 + \beta t$ | +| Decay | $t > t_1$ | $\log y(t) = \frac{1}{1-\rho}\log\!\bigl[y_1^{1-\rho} - (1-\rho)\,\alpha\,(t - t_1)\bigr]$ | + +where $\beta = (\log y_1 - \log y_0)/t_1$. + +### Model B likelihood + +$$ +\log \mathbf{y}_{i,t} \sim \mathcal{N}_K\!\left(\boldsymbol{\mu}_{i,t},\; \Sigma_\varepsilon\right) +$$ + +$\Sigma_\varepsilon = \mathrm{diag}(\boldsymbol{\tau}_\varepsilon)\,\Omega_\varepsilon\,\mathrm{diag}(\boldsymbol{\tau}_\varepsilon)$, +with $\Omega_\varepsilon \sim \mathrm{LKJ}(2)$ and +$\tau_{\varepsilon,k} \sim \mathrm{Cauchy}^+(0, 2.5)$. + +## Installation + +Install `rstan` if you have not already: + +```{r} +install.packages("rstan") +``` + +Install `serodynamics`: + +```{r} +# install.packages("pak") +pak::pak("UCD-SERG/serodynamics") +``` + +## Quick start + +```{r} +library(serodynamics) +library(dplyr) + +# Simulate 20 subjects using typhoid curve parameters +set.seed(42) +sim_data <- serocalculator::typhoid_curves_nostrat_100 |> + sim_case_data(n = 20) +``` + +### Fit Model B + +```{r} +fit <- run_mod_stan( + data = sim_data, + model = "model_b", + chains = 4, + iter = 2000, + warmup = 1000, + seed = 1 +) +``` + +### Inspect output + +```{r} +# Same structure as run_mod() output +print(fit) + +# Residual correlation matrix across antigen-isotype pairs +attr(fit, "Omega_eps") + +# Residual covariance matrix +attr(fit, "Sigma_eps") +``` + +### Fitted and residual values + +```{r} +attr(fit, "fitted_residuals") |> head() +``` + +### Visualize predicted curves + +Because `fit` is an `sr_model` object, the standard plotting functions work +directly: + +```{r} +plot_predicted_curve(fit, sim_data) +``` + +### Posterior summaries + +```{r} +post_summ(fit) +``` + +## Comparing JAGS and Stan backends + +Both `run_mod()` and `run_mod_stan()` return an `sr_model` tibble, so you can +compare them side-by-side: + +```{r} +jags_fit <- run_mod( + data = sim_data, + nchain = 4, + nadapt = 500, + nburn = 500, + nmc = 500, + niter = 500 +) + +# Column names are identical +identical(names(jags_fit), names(fit)) + +# Compare posterior medians for y0 +bind_rows( + jags_fit |> filter(Parameter == "y0") |> mutate(backend = "JAGS"), + fit |> filter(Parameter == "y0") |> mutate(backend = "Stan") +) |> + group_by(backend, Iso_type) |> + summarise(median_y0 = median(value), .groups = "drop") +``` + +## Data flow + +``` +as_case_data(raw) # unchanged + │ + ▼ +prep_data(case_data, # unchanged (add_newperson = FALSE used internally) + add_newperson = FALSE) + │ + ▼ +prep_data_stan(prepped) # NEW — reshape to Stan list + │ + ▼ +rstan::sampling(model_b) # Stan NUTS sampler + │ + ▼ +postprocess_stan_output() # NEW — stanfit → sr_model tibble + │ + ▼ +sr_model tibble # identical structure to run_mod() output +``` + +## Session info + +```{r} +sessioninfo::session_info() +```