Add Stan support as alternative Bayesian modeling backend

[Copilot]

Convert JAGS Model to Stan ✅

This task successfully adds Stan support to the serodynamics package as an alternative to JAGS for Bayesian MCMC modeling. Both modeling backends can coexist, allowing users to choose their preferred framework.

Completed Work

• Research Stan ecosystem and best practices

  • Identified cmdstanr as the recommended Stan R interface (not rstan)
  • Reviewed Stan modeling language syntax and conversion patterns
  • Understood key differences: precision vs covariance, explicit blocks, inverse Wishart conversion

• Create Stan model files

  • Converted model.jags to inst/extdata/model.stan (117 lines)
  • Converted model.dobson.jags to inst/extdata/model.dobson.stan (19 lines)
  • Handled precision-to-covariance conversions for multivariate normal priors
  • Used inverse Wishart in Stan (equivalent to Wishart on precision in JAGS)
  • Properly structured Stan model with data, parameters, transformed parameters, and model blocks
  • Implemented ifelse logic for active vs recovery infection phases
  • Optimized to use multi_normal_prec() for better performance and numerical stability

• Create R infrastructure for Stan

  • Created prep_data_stan() function (75 lines) - converts data to Stan format with proper padding
  • Created prep_priors_stan() function (78 lines) - prepares Stan-compatible prior specifications
  • Created run_mod_stan() function (197 lines) - complete Stan equivalent of run_mod()
  • Created sample_predictive_stan() function - generates posterior predictive samples from fitted Stan models
  • Made cmdstanr optional (in Suggests) with clear runtime error messages
  • Stan functions integrate seamlessly with existing workflow
  • Extracted shared helper functions to reduce code duplication (~150 lines)

• Update dependencies and documentation

  • Added cmdstanr and badger to DESCRIPTION Suggests
  • Added stan-dev/cmdstanr to Remotes for easy installation
  • Generated comprehensive roxygen2 documentation
  • Created example scripts for all Stan functions
  • Updated NEWS.md with new features
  • Incremented version to 0.0.0.9048
  • Fixed all linting issues (0 warnings)
  • Updated WORDLIST (added: burnin, cmdstanr, stan, stanfit)
  • Passed all spelling checks

• Update README and documentation

  • Updated README.Rmd to document both JAGS and Stan support
  • Added Stan installation instructions
  • Implemented badger package for maintainable badge generation
  • Updated repository URLs from UCD-SERG to ucdavis organization
  • Rendered README.md
  • Updated copilot-instructions.md comprehensively
  • Documented Stan model files in repository structure
  • Clarified JAGS is required only for JAGS models

• Package validation

  • Successfully passed R CMD check (0 errors, 0 warnings, 0 notes)
  • All examples run without errors
  • All tests pass (63 total: 55 existing + 8 new Stan tests)
  • Package builds and installs correctly
  • Documentation builds correctly

• CI/Build improvements

  • Added *.knit.* pattern to .gitignore to exclude intermediate knit files
  • Removed README.knit.md from git tracking
  • Added cmdstanr installation to copilot-setup-steps.yml workflow
  • Updated verification step to include cmdstanr in key packages list
  • Fixed R-CMD-check workflow failure by wrapping Stan examples in \dontrun{}

• Test coverage

  • Added comprehensive tests for prep_priors_stan() (output structure, modifiability)
  • Added comprehensive tests for prep_data_stan() (NA validation, data padding)
  • Added comprehensive tests for run_mod_stan() (error handling, output structure, stratification support)
  • All tests use skip_if_not_installed() for CI portability
  • Tests skip gracefully when CmdStan not available

Implementation Details

Stan Model Conversion:

• Converted JAGS precision parameterization to Stan standard deviation
• Transformed dnorm(mu, tau) to normal(mu, sigma) where sigma = 1/sqrt(tau)
• Converted dwish(omega, df) to inv_wishart(df, omega_inv) where omega_inv = inverse(omega)
• Converted dgamma(a, b) parameters (same in both)
• Used explicit typed declarations for all variables
• Implemented conditional logic for infection phases using if/else instead of JAGS step function
• Uses multi_normal_prec() to accept precision matrices directly (better performance and numerical stability)
• Only inverts omega for inverse Wishart (required by Stan)

Data Preparation:

• Stan requires padded arrays (no ragged arrays) - prep_data_stan() handles this
• Maximum sample size determined and arrays padded with zeros
• Loop limits controlled by nsmpl[subj] to ignore padding
• Validates original input data for NA values (not padded JAGS arrays)
• Stan cannot handle NA values - no add_newperson parameter in prep_data_stan()

Posterior Predictive Sampling:

• New sample_predictive_stan() function for generating predictions from fitted Stan models
• Samples from the marginal posterior distribution of model parameters
• Generates predictions at user-specified time points
• Returns both raw samples and summary statistics (mean, median, 95% credible intervals)
• Separate from model fitting to maintain clean separation of concerns

Code Quality:

• Extracted shared logic into helper files (prep_priors_helpers.R, run_mod_helpers.R)
• Backend-neutral error messages in shared validation functions
• Proper handling of NA stratification levels
• Preserves stratification column types (Date, POSIXct, factor, etc.)
• Factor-safe list indexing with as.character()
• Robust attribute selection using named selection instead of positional indexing

Integration:

• Stan functions return same output format as JAGS functions (sr_model class)
• Works with existing post-processing and plotting functions
• Uses same prior specification interface
• Maintains backward compatibility - no breaking changes

CI Environment:

• Copilot setup workflow now installs cmdstanr and CmdStan automatically
• Ensures consistent development environment for Stan-based work
• Verifies cmdstanr installation in environment checks

Example Code:

• Example files use plain R code with runtime checks (requireNamespace + cmdstan_version)
• Roxygen wraps examples with \dontrun{} for R CMD check compatibility
• Examples require CmdStan to be installed
• Users can run examples interactively after installing cmdstanr and CmdStan

What Users Can Do Now

1. Use JAGS (existing functionality):

   run_mod(data, file_mod = "model.jags", nchain = 4, ...)

2. Use Stan (new functionality):

   # Fit the model
   model <- run_mod_stan(data, file_mod = "model.stan", nchain = 4, ...)
   
   # Generate posterior predictive samples
   predictions <- sample_predictive_stan(model, time_points = c(5, 30, 90))

Both JAGS and Stan produce compatible outputs that work with the same analysis and visualization functions.

Note: Stan requires complete data (no NA values). For posterior predictive sampling with Stan models, use the sample_predictive_stan() helper function which samples from the marginal posterior distribution of model parameters.

Original prompt

---

This section details on the original issue you should resolve

<issue_title>convert jags model to stan</issue_title>
<issue_description></issue_description>

Comments on the Issue (you are @copilot in this section)

• Fixes convert jags model to stan #174

---

💬 We'd love your input! Share your thoughts on Copilot coding agent in our 2 minute survey.

[d-morrison]

fix failing workflow checks

[d-morrison]

address workflow check failures

[d-morrison]

fix failing workflows

[claude]

Broken reference — hallucination risk: model_stan_transform_from_jags.md does not exist in this repository. Searching the whole tree finds no such file. The companion derivation lives in vignettes/articles/jags-to-stan.qmd. Please update both this line and line 76 (same reference) to point to the real document.

Suggested change

	// reparameterizations in the companion note model_stan_transform_from_jags.md.
	// reparameterizations in the companion vignette vignettes/articles/jags-to-stan.qmd.

[claude]

Same broken reference — model_stan_transform_from_jags.md does not exist. Update to reference vignettes/articles/jags-to-stan.qmd.

Suggested change

	// lognormal factorization (derivation in model_stan_transform_from_jags.md).
	// Moment-match the JAGS inverse-Wishart covariance prior to an LKJ +
	// lognormal factorization (derivation in vignettes/articles/jags-to-stan.qmd).

[claude]

Numerical instability near shape ≈ 1: this formula is algebraically equivalent to the Stan model's log_two_phase function but less numerically stable. When shape_pop is close to 1 (i.e. one_minus_shape ≈ 0), log_arg → 1 and log(1)/one_minus_shape → 0/0, producing NaN for those posterior draws. The pmax() guard prevents log(negative) but not the 0/0 singularity.

The Stan model avoids this with log1p, which handles the a → 0 removable singularity smoothly (see the functions {} block in model.stan). Use the same formula here:

Suggested change

	{
	one_minus_shape <- 1 - shape_pop
	# Compute the argument to log, ensuring it stays positive
	log_arg <- y1_pop^one_minus_shape -
	one_minus_shape * alpha_pop * (t - t1_pop)
	# Use pmax to ensure log_arg > 0 (avoid NaN from log of negative)
	(1 / one_minus_shape) * log(pmax(log_arg, .Machine$double.eps))
	}
	{
	# Use the numerically stable log1p form that matches model.stan's
	# log_two_phase(): log_y1 - log1p(a*alpha*(t-t1)*y1^a) / a
	# Stable as a = shape-1 -> 0 (no 0/0 singularity).
	a <- shape_pop - 1 # > 0 since shape > 1
	log(y1_pop) - log1p(a * alpha_pop * (t - t1_pop) * y1_pop^a) / a
	}

[claude]

Misleading ids attribute in stratified runs: longdata here is always from the last stratum processed, so attr(stan_out, "ids") will only contain the last stratum's subject IDs, not the full set. While sample_predictive_stan() currently ignores this attribute (it only reads antigens), anyone who inspects attr(result, "ids") after a stratified run will get an incomplete list.

Consider either (a) accumulating all IDs across strata:

all_ids <- c(all_ids, attr(longdata, "ids"))  # inside the loop
# ...
attr(stan_out, "ids") <- all_ids

or (b) setting the attribute to NULL / documenting clearly that it only reflects the last stratum, to prevent silent misuse.

[claude]

Claude finished @Kwan-Jenny's task in 5m 16s —— View job

---

PR #175 Code Review: Add Stan support as alternative Bayesian modeling backend

• Check what changed since last review (commit 6a60ad4)
• Read key changed files: model.stan, sample_predictive_stan.R, run_mod_stan.R, run_mod_helpers.R, prep_data_stan.R, prep_priors_helpers.R
• Verify prior findings from last review
• Post inline findings and final verdict

---

Status: 3 prior findings still unresolved — no new findings

The previous Claude review (2026-06-27) identified three findings. All three are confirmed still present in the current HEAD after reviewing the commits since 6a60ad4 (7dde102, c6c4ed9, e59d994, 9fa3c18, a4bd849, 02af72b, c850326, 63cb294, 3d7a25e). All three are re-posted as inline comments.

---

Unresolved Inline Findings

1. inst/extdata/model.stan:11,76 — Broken reference (hallucination risk)

Both lines say "companion note model_stan_transform_from_jags.md". This file does not exist in the repository. It was introduced in c6c4ed9 when the header block was rewritten — but the same commit also added the real derivation document at vignettes/articles/jags-to-stan.qmd. Fix: update both comments to reference vignettes/articles/jags-to-stan.qmd.

2. R/sample_predictive_stan.R:204–211 — Numerically unstable recovery-phase formula

The R code uses:

one_minus_shape <- 1 - shape_pop
log_arg <- y1_pop^one_minus_shape - one_minus_shape * alpha_pop * (t - t1_pop)
(1 / one_minus_shape) * log(pmax(log_arg, .Machine$double.eps))

Near shape ≈ 1, log_arg → 1 and log(1)/0 → 0/0 → NaN. The pmax guard protects against log(negative) only, not this singularity. Worse, if pmax does clamp, dividing log(.Machine$double.eps) ≈ -36 by a small negative one_minus_shape gives a large positive result — the wrong sign. The Stan model (log_two_phase, line 53) uses log1p to handle both issues. The R code should match:

a <- shape_pop - 1
log(y1_pop) - log1p(a * alpha_pop * (t - t1_pop) * y1_pop^a) / a

3. R/run_mod_stan.R:274 — ids attribute reflects last stratum only

In stratified runs longdata is always the last iteration's value, so attr(stan_out, "ids") contains only the last stratum's subject IDs. sample_predictive_stan() doesn't use this attribute, so no current functional bug, but the attribute is misleading for callers. Either accumulate IDs across strata inside the loop, set to NULL, or document the limitation explicitly.

---

Changes Since Prior Review — All Confirmed Correct

Commit	Change	Verdict
c6c4ed9	Full Stan model rewrite: LKJ + lognormal replaces inverse-Wishart; log1p in log_two_phase	✓ Math correct (moment-matching derivation verified)
7dde102	nBurnin ← nadapt; removed noisy cli_inform; seq_len(nchain); improved @returns	✓ All prior P1 findings resolved
9fa3c18	Added vignettes/articles/jags-to-stan.qmd	✓ Correct
e59d994	Removed unused copilot scaffolding	✓
02af72b–3d7a25e	Version bump, WORDLIST, pkgdown index	✓

CLAUDE.md Requirements — All Met

Requirement	Status
dplyr::*_join() must set relationship	✓ Both joins in run_mod_helpers.R use relationship = "many-to-one"
seq_len() not 1:n	✓ Fixed in 7dde102; initialize_prior_arrays also uses seq_len
No .data in selection contexts	✓ .data[[strat]] only inside dplyr::filter() (data-masking verb)
tibble::as_tibble() not dplyr::as_tibble()	✓ Correct in format_model_output()
No base::merge()	✓
No deprecated lifecycle calls	✓

[claude]

Still unresolved from prior review — broken reference (hallucination risk): model_stan_transform_from_jags.md does not exist anywhere in the repository. Searching the full tree finds no such file. This reference was introduced in c6c4ed9 when the header was rewritten. The companion derivation was added simultaneously as vignettes/articles/jags-to-stan.qmd — update this line (and line 76 below) to point to that real document.

Suggested change

	// reparameterizations in the companion note model_stan_transform_from_jags.md.
	// reparameterizations in the companion vignette vignettes/articles/jags-to-stan.qmd.

[claude]

Still unresolved from prior review — same broken reference: model_stan_transform_from_jags.md does not exist. The derivation lives in vignettes/articles/jags-to-stan.qmd.

Suggested change

	// lognormal factorization (derivation in model_stan_transform_from_jags.md).
	// Moment-match the JAGS inverse-Wishart covariance prior to an LKJ +
	// lognormal factorization (derivation in vignettes/articles/jags-to-stan.qmd).

[claude]

Still unresolved from prior review — numerically unstable recovery formula: This implements the same expression the Stan model (log_two_phase) explicitly avoids. Near shape ≈ 1 (i.e., one_minus_shape → 0), log_arg → 1, so log(log_arg) → 0, and you get (1/0) * 0 → NaN. The pmax() guard protects against log(negative) only, not this 0/0 removable singularity.

Additionally, when pmax does clamp (e.g. large alpha*(t-t1)), dividing log(.Machine$double.eps) ≈ -36 by a small negative one_minus_shape gives a large positive result — the wrong sign.

The Stan model uses log1p to handle both issues cleanly (see functions {} block, line 53). Use the same formula here:

Suggested change

	{
	one_minus_shape <- 1 - shape_pop
	# Compute the argument to log, ensuring it stays positive
	log_arg <- y1_pop^one_minus_shape -
	one_minus_shape * alpha_pop * (t - t1_pop)
	# Use pmax to ensure log_arg > 0 (avoid NaN from log of negative)
	(1 / one_minus_shape) * log(pmax(log_arg, .Machine$double.eps))
	}
	{
	# Numerically stable log1p form matching model.stan's log_two_phase().
	# Stable as a = shape-1 -> 0 (removes the 0/0 singularity).
	a <- shape_pop - 1 # > 0 since shape_pop = exp(...) + 1 > 1
	log(y1_pop) -
	log1p(a * alpha_pop * (t - t1_pop) * y1_pop^a) / a
	}

[claude]

Still unresolved from prior review — ids attribute reflects last stratum only: longdata here is always the result from the final loop iteration, so attr(stan_out, "ids") only contains that stratum's subject IDs. sample_predictive_stan() doesn't currently use this attribute, but any caller inspecting attr(result, "ids") after a stratified fit gets an incomplete and misleading list.

Options:

• Accumulate IDs inside the loop: all_ids <- c(all_ids, attr(longdata, "ids")) and attach all_ids here.
• Set to NULL and document clearly that it is not populated for stratified models.
• Document the limitation with an explicit comment so the next maintainer doesn't rely on it.