Add serodynamics Stan model for Chapter 2 (Model 2a, Chapter1+α)

[Kwan-Jenny]

Summary

Adds the Chapter 2 Stan model Model 2a ("Chapter1+α") to serodynamics, plus its runner and supporting functions. Model 2a is the honest generalization of the Chapter 1 model: it keeps Chapter 1's two free 5×5 within-biomarker
covariance blocks unchanged and adds only the 5 same-parameter cross-biomarker covariances (a diagonal cross-block). Setting those 5 terms to zero recovers Chapter 1 exactly, so Model 2a strictly nests Chapter 1 (35 vs 30 covariance
parameters).

This PR is built on main and is self-contained: it does not depend on the (unmerged) Chapter 1 Stan PR.

What this enables

A subject's IgA and IgG curve parameters can now be correlated (e.g. a person with fast IgA decay also tends to have fast IgG decay), via the same-parameter cross-biomarker terms c = (c_y0, c_y1, c_t1, c_alpha, c_r). Chapter 1 forces these to zero (biomarkers independent).

The covariance structure

        IgG block      IgA block
      ┌────────────┬────────────┐
 IgG  │  Σ_G full  │  diag(c)   │   ← only the diagonal cross-block is new
      ├────────────┼────────────┤
 IgA  │  diag(c)   │  Σ_A full  │
      └────────────┴────────────┘
        30 (within)  +  5 (cross)  =  35 params      (C = 0  ⇒  Chapter 1)

Positive-definiteness is guaranteed for any c via a conditional (Schur) construction: θ_IgG ~ MVN(μ_G, Σ_G), θ_IgA | θ_IgG ~ MVN(μ_A + B(θ_IgG − μ_G), Ψ) with B = C Σ_G⁻¹. This yields a diagonal cross-block exactly and reduces to the block-diagonal Chapter 1 model when C = 0.

Why this structure (and not the alternatives)

• Unstructured (55 params): too many for n ≈ 48; discards interpretability.
• Kronecker Σ_B ⊗ Σ_P (17 params): ties the two biomarkers to one shared within-block correlation, and is not nested with Chapter 1 — so the "Chapter 1 converges, therefore this should too" argument does not apply.
• Chapter1+α (35 params): keeps all of Chapter 1's within-biomarker flexibility and adds the minimum cross-biomarker term. The only option here that strictly nests Chapter 1.

Changes (committed)

file	description
inst/extdata/model_2a.stan	Model 2a (Chapter1+α): conditional/Schur construction, non-centered, LKJ + lognormal priors matched to Chapter 1, Normal(0, ·) prior on the 5 cross-couplings
R/run_serodynamics_stan_2a.R	runner; fits model_2a.stan. file_mod lets calling code point at an alternative model with the same data interface. Requires exactly 2 biomarkers.
R/prep_data_serodynamics_stan_2a.R	case data → padded Stan list (reuses prep_data() from main)
R/prep_priors_serodynamics_stan_2a.R	per-antigen prior inputs (same as the Chapter 1 Stan model)
R/serodynamics_stan_2a_helpers.R	summarise_pop_2a(), log_two_phase_r(), predict_newperson_2a() (conditional new-person prediction)

Naming follows the lab rename run_mod() → run_serodynamics(); the Chapter 2 runner uses the new convention (run_serodynamics_stan_2a).

Priors

Matched to Chapter 1, with Chapter 1 at the prior mode:

• Σ_G and Ψ (the conditional/Schur covariance of the IgA block) each get the Chapter 1 per-antigen LKJ + lognormal factorization. At C = 0, Ψ = Σ_A, so the prior is exactly Chapter 1's at the nesting point.
• the 5 cross-couplings c_p ~ Normal(0, scale_p), centred at 0 (so the model leans toward "no cross-correlation" until the data say otherwise). scale_p is derived from the matched marginal scales inside the model.

Sampler defaults copy the Chapter 1 values validated on this likelihood (adapt_delta = 0.95, max_treedepth = 12, default init).

Dependencies

• cmdstanr (as in the Chapter 1 Stan work).
• Reuses prep_data() and sim_case_data() from main.

Validation (not in this PR)

Validation is done locally and is intentionally not committed (validation/ is in .gitignore / .Rbuildignore). It generates a no-correlation control and two correlated datasets (same marginals, same truth) and fits both the Chapter 1 baseline and Model 2a on each, checking:

• no-correlation → the two models agree (Model 2a shows no advantage), and
• with correlation → Model 2a recovers the IgA–IgG decay correlation and gives a tighter, closer-to-truth new-person prediction.

A write-up will follow as an articles/*.qmd once the results are finalized (separate from this PR).

Notes for review

• This is built on main; once the Chapter 1 Stan PR merges, the temporary duplicates (prep_*_serodynamics_stan_2a, and the validation-only Chapter 1 baseline) can be replaced by the shared prep_*_stan / run_mod_stan. The runner already takes file_mod and isolates its prep, so that switch is a small change.
• model_2a.stan mirrors the validated Chapter 1 model.stan as closely as possible — only the covariance construction changes.

[claude]

Claude finished @Kwan-Jenny's task in 4m 32s —— View job

---

Code Review: PR #259 — Model 2a Stan (Chapter1+α)

Todo

• Load tools
• Read changed R files
• Read Stan model
• Read documentation files
• Analyze for correctness, idioms, deprecations, hallucinations
• Post inline findings
• Post summary

---

Overall verdict

The statistical construction is sound — the conditional/Schur factorization is correctly implemented and the B = C Σ₁⁻¹ formula is implemented correctly (mdivide_right_spd, quad_form_sym with the transposed argument). The nesting claim holds. Two issues need fixing before merge; the rest are lower-priority.

---

Must-fix

1. cmdstanr not declared in DESCRIPTION (R/run_serodynamics_stan_2a.R:86)

cmdstanr::cmdstan_model() is the runner's hard dependency but cmdstanr appears nowhere in DESCRIPTION. R CMD check will flag an undeclared namespace. Add to Suggests: and guard with rlang::check_installed("cmdstanr") before the call.

2. log_y1 vs log_y1y0 naming inconsistency (R/serodynamics_stan_2a_helpers.R:21)

summarise_pop_2a labels Stan parameter par[k][2] as "log_y1y0" (= log(y1 − y0), the log-differential). log_two_phase_r names its second curve parameter log_y1, expecting the log of the absolute peak log(y0 + exp(par[2])). These are different quantities — feeding summarise_pop_2a output directly into log_two_phase_r produces silently wrong curves. Fix by either aligning the parameter names or prominently documenting the required conversion (log_y1 = log(exp(log_y0) + exp(log_y1y0))).

---

Should-fix

3. Dead code / misleading "strip attributes" comment (R/run_serodynamics_stan_2a.R:83)

c(longdata, priorspec_clean) already drops all non-names attributes (class, antigens, n_antigens, ids) — R's c() on lists does not propagate extra attributes. The subsequent filter checks for element names that don't exist, making it a no-op. The misleading comment is the real problem: a future reader may rely on it when debugging attribute-related regressions.

4. given length not validated in predict_newperson_2a (R/serodynamics_stan_2a_helpers.R:85)

A wrong-length given fails inside the matrix multiply with an opaque dimension-mismatch error. A simple length guard at the top of the function would give a much clearer message.

5. No tests

The new data-prep, prior-prep, and curve helpers have no unit tests. Given that the test suite runs with lifecycle_verbosity = "error", even a small set of smoke tests (data round-trip, curve value at t=0 and t=t1, invalid-input errors) would increase confidence and catch regressions when the Chapter 1 Stan functions eventually replace the duplicated prep_* code.

---

Minor / informational

6. chol() result named L (upper, not lower triangular) (R/serodynamics_stan_2a_helpers.R:132)
R's chol() returns the upper factor U; t(L) %*% z is mathematically correct, but the variable name L implies lower triangular (the Stan/LAPACK convention). Rename to U or add a comment.

7. log_two_phase_r docstring says "identical to Stan" — signatures differ (R/serodynamics_stan_2a_helpers.R:51)
The Stan function takes a pre-computed y1 argument for efficiency; the R function does not. The formula is the same; the interface is not. "Implements the same formula as" is more accurate than "identical to."

8. cross_cor not clamped to [-1, 1] (inst/extdata/model_2a.stan:213)
Theoretically bounded, but floating-point accumulation can push values slightly outside the range. A fmax(-1, fmin(1, ...)) clamp would prevent downstream surprises.

9. quad_form_sym(Sigma_1, B') orientation confirmed correct (inst/extdata/model_2a.stan:197)
quad_form_sym(A, B) = B' A B, so quad_form_sym(Sigma_1, B') = B Σ₁ B' ✓ — a brief inline comment would help future readers.

---

No hallucinations found

mdivide_right_spd, quad_form_sym, multiply_lower_tri_self_transpose, lkj_corr_cholesky, diag_pre_multiply — all verified as real Stan Math functions. cmdstanr::cmdstan_model, prep_data(), get_timeindays_var, serocalculator::get_values_var — all verified as real functions in the dependency chain.

[codecov]

Codecov Report

❌ Patch coverage is 0% with 178 lines in your changes missing coverage. Please review.
✅ All tests successful. No failed tests found.

Files with missing lines	Patch %	Lines
R/serodynamics_stan_2a_helpers.R	0.00%	55 Missing ⚠️
R/prep_data_serodynamics_stan_2a.R	0.00%	46 Missing ⚠️
R/prep_priors_serodynamics_stan_2a.R	0.00%	39 Missing ⚠️
R/run_serodynamics_stan_2a.R	0.00%	38 Missing ⚠️

Files with missing lines	Coverage Δ
R/run_serodynamics_stan_2a.R	0.00% <0.00%> (ø)
R/prep_priors_serodynamics_stan_2a.R	0.00% <0.00%> (ø)
R/prep_data_serodynamics_stan_2a.R	0.00% <0.00%> (ø)
R/serodynamics_stan_2a_helpers.R	0.00% <0.00%> (ø)

[claude]

cmdstanr is not declared in DESCRIPTION.

cmdstanr::cmdstan_model() is called here (and the whole runner fails without it), but cmdstanr appears in neither Imports nor Suggests in DESCRIPTION. R CMD check will flag an undeclared namespace, and users without cmdstanr will get an opaque error.

Add to Suggests in DESCRIPTION and guard with rlang::check_installed("cmdstanr") before this call (or earlier in the function), following the pattern the package uses for other optional dependencies.

[claude]

Dead code — attribute stripping is already done by c().

c(longdata, priorspec_clean) returns a plain list: R's c() on list objects drops all non-names attributes (class, antigens, n_antigens, ids) from the inputs. By the time this filter runs, stan_data is already a plain named list and none of the filtered names appear as element names, so the filter is a no-op.

The comment "strip R-side attributes" is therefore misleading; a future reader may trust it and assume it works, masking a regression if the data-prep function ever returns those fields as proper list elements. Either remove the filter (with a comment explaining why c() already strips attributes), or make the stripping explicit with attributes(stan_data) <- list(names = names(stan_data)) if belt-and-suspenders safety is desired.

[claude]

Naming inconsistency with log_two_phase_r: log_y1y0 ≠ log_y1.

summarise_pop_2a labels parameter j = 2 as "log_y1y0", which is Stan's par[k][2] = log(y1 − y0) (the log of the differential).

log_two_phase_r (line 51) takes a parameter named log_y1, expecting the log of the absolute peak, i.e. log(y1) = log(y0 + exp(par[2])).

These are different quantities. A user who feeds summarise_pop_2a output directly into log_two_phase_r (a natural workflow) will get silently wrong curves. The required conversion is:

log_y1 <- log(exp(log_y0) + exp(log_y1y0))

Fix by either:

• Renaming the log_two_phase_r parameter from log_y1 to log_y1y0 and updating its body (y1 <- exp(log_y0) + exp(log_y1y0)), or
• Documenting the conversion prominently in both functions' @param docs.

[claude]

Docstring claims "identical to the Stan log_two_phase" but the signatures differ.

The Stan function signature is log_two_phase(t, log_y0, log_y1, y1, t1, alpha, shape) — it takes y1 (the actual peak value) as a pre-computed argument. The R function omits that argument and computes y1 <- exp(log_y1) internally.

The functions compute the same formula, but the signatures are intentionally different (Stan takes the pre-computed y1 for efficiency). The docstring should say "implements the same formula" rather than "identical to" to avoid confusion.

[claude]

given length is not validated.

If given is not NULL and not exactly length n_params, the error surfaces inside the matrix multiply as an opaque dimension mismatch. A guard here would give a much clearer message:

if (!is.null(given) && length(given) != n_params) {
  stop("`given` must be NULL or a length-", n_params, " vector; got length ", length(given), ".")
}

[claude]

chol() returns the upper Cholesky factor; t(L) is correct but the variable name L implies lower triangular.

R's chol(S) returns upper triangular U such that S = U'U. The sampling step mu + t(L) %*% z = mu + U' z is mathematically correct for drawing from N(mu, S). But naming the result L (conventional for the lower Cholesky factor) is misleading. Rename to U or add an inline comment explaining the convention.

[claude]

quad_form_sym orientation is correct — confirming for reviewers.

quad_form_sym(A, B) computes B' A B, so quad_form_sym(Sigma_1, B') = (B')' Sigma_1 B' = B Sigma_1 B', which matches the Schur formula Sigma_2 = B Sigma_1 B' + Psi. ✓

One note: quad_form_sym requires A to be symmetric. Sigma_1 = multiply_lower_tri_self_transpose(L_1) is PSD by construction, so this is fine. A future reader may want a brief inline comment explaining the B' argument.

[claude]

cross_cor is not clamped to [-1, 1].

The Schur construction guarantees the theoretical joint covariance is PD, so |cross_cor[p]| ≤ 1 in exact arithmetic. In practice, floating-point accumulation (especially when c_cross[p] is large relative to the geometric mean of the diagonal variances) can push a computed value slightly outside [-1, 1]. Downstream R code that interprets cross_cor as a correlation may be surprised.

Consider adding a post-hoc clamp or at least noting this in the generated quantities block:

cross_cor[p] = fmax(-1.0, fmin(1.0, c_cross[p] / sqrt(Sigma_1[p, p] * Sigma_2[p, p])));