Skip to content

Conversation

@Kwan-Jenny
Copy link
Collaborator

This adds my Chapter 2 .qmd file.

It includes:

  • simulation of fake data with known covariance
  • fitting the independence model (baseline)
  • setting up the correlated model with Kronecker precision
  • wrapper function run_mod_kron for running the model

I am mainly stuck on the prep_priors_multiB() and run_mod_kron() parts, and would like feedback on whether the priors and wrapper are set up correctly.

@codecov
Copy link

codecov bot commented Sep 23, 2025

Codecov Report

❌ Patch coverage is 87.70053% with 23 lines in your changes missing coverage. Please review.

Files with missing lines Patch % Lines
R/Run_Mod.R 41.37% 17 Missing ⚠️
R/simulate_multi_b_long.R 93.75% 4 Missing ⚠️
R/prep_priors_multi_b.R 88.88% 2 Missing ⚠️
Files with missing lines Coverage Δ
R/clean_priors.R 100.00% <100.00%> (ø)
R/inits_kron.R 100.00% <100.00%> (ø)
R/write_model_ch2_kron.R 100.00% <100.00%> (ø)
R/prep_priors_multi_b.R 88.88% <88.88%> (ø)
R/simulate_multi_b_long.R 93.75% <93.75%> (ø)
R/Run_Mod.R 84.68% <41.37%> (-15.32%) ⬇️

@github-actions
Copy link
Contributor

github-actions bot commented Sep 23, 2025

📖 https://ucd-serg.github.io/serodynamics/preview/pr133
Preview documentation for this PR (at commit 9279040)

Copy link
Collaborator

@d-morrison d-morrison left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Kwan, sorry to push this review back another day; I'll try to take a detailed look tomorrow. It would really help me if you could move all these functions into the R/ folder, lint them, and document them with short examples, so I can read and test them individually. My attention span is badly limited right now, and it really helps to have the code in bite-size chunks.

@Kwan-Jenny
Copy link
Collaborator Author

Hi Kwan, sorry to push this review back another day; I'll try to take a detailed look tomorrow. It would really help me if you could move all these functions into the R/ folder, lint them, and document them with short examples, so I can read and test them individually. My attention span is badly limited right now, and it really helps to have the code in bite-size chunks.

Hi @d-morrison ,

I will do exactly as you suggested and let you know once it’s ready, so you can start the review then. I really appreciate you taking the time to review this!

- Added functions: two_phase_y, simulate_multiB_long, simulate_multiB_long2,
  make_ab_vec, prep_priors_multiB, inits_kron, clean_priors,
  write_model_ch2_kron, run_mod_kron
- Added example scripts under inst/examples for each function
- Refactored code from Chapter 2 QMD into R/ with docs + examples
@Kwan-Jenny
Copy link
Collaborator Author

Hi Kwan, sorry to push this review back another day; I'll try to take a detailed look tomorrow. It would really help me if you could move all these functions into the R/ folder, lint them, and document them with short examples, so I can read and test them individually. My attention span is badly limited right now, and it really helps to have the code in bite-size chunks.

Hi @d-morrison ,

I will do exactly as you suggested and let you know once it’s ready, so you can start the review then. I really appreciate you taking the time to review this!

Hi @d-morrison
I finished moving the functions into R/, added docs and short examples, and linted them. It’s ready for you to review now. Thanks again for your help!

Copy link
Collaborator

@d-morrison d-morrison left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please finish moving all functions out of the .qmd file and make sure all github checks pass.

- Deleted two_phase_y and make_ab_vec functions (and their examples/docs),
  since simulation now directly uses serodynamics:::ab().
- Kept single simulate_multi_b_long() function with serodynamics:::ab()
  as the unified approach.
@Kwan-Jenny
Copy link
Collaborator Author

Please finish moving all functions out of the .qmd file and make sure all github checks pass.

Hi, @d-morrison

Done. I moved the remaining code from the qmd into R/ files, added roxygen, and updated the pkgdown reference. Locally devtools::check() now finishes cleanly (0 errors | 0 warnings | 0 notes). The examples that run JAGS are wrapped in \dontrun{} so CI doesn’t time out, but the model compiles.

What changed:

  • New/updated R files:

    • R/write_model_ch2_kron.R – JAGS model writer. Uses a single scalar n_blocks everywhere (replaces old B).

    • R/run_mod_kron.R – wrapper modeled after run_mod(). Passes n_blocks into JAGS via list(n_blocks = longdata$n_antigen_isos).

    • R/prep_priors_multi_b.R – builds OmegaP/nuP and OmegaB/nuB (Wishart pieces) with input validation.

    • R/simulate_multi_b_long.R – simulator (now documented and exported).

  • Fixed the JAGS “redefine node” error by computing the growth-rate term inside the observation loop with a unique temp variable (no reused node names).

  • Updated _pkgdown.yml reference index and examples under inst/examples/.

What I would like reviewed specifically:

  1. write_model_ch2_kron() – data list that JAGS expects (n_blocks, mu.hyp, prec.hyp, OmegaP/nuP, OmegaB/nuB, measurement pieces) and the Tau construction (TauB ⊗ TauP).

  2. run_mod_kron() – how the priors are combined (prep_priors() + prep_priors_multi_b()) and the monitor set (c("y0","y1","t1","alpha","shape","TauB","TauP")).

Thanks!
Kwan

@Kwan-Jenny
Copy link
Collaborator Author

Please finish moving all functions out of the .qmd file and make sure all github checks pass.

Hi @d-morrison ,
quick clarification before I proceed:

Do you want me to get all GitHub Actions checks (R-CMD-check on macOS/Windows/Linux, docs, coverage) green before review/approval, or is a clean local devtools::check() sufficient?

Right now devtools::check() passes locally, but the PR shows some CI failures. If your expectation is CI-green, I will keep tightening the examples/tests (e.g., wrap JAGS bits in \dontrun{}/skip on CI, fix pkgdown) until everything passes.

@d-morrison
Copy link
Collaborator

Please finish moving all functions out of the .qmd file and make sure all github checks pass.

Hi @d-morrison , quick clarification before I proceed:

Do you want me to get all GitHub Actions checks (R-CMD-check on macOS/Windows/Linux, docs, coverage) green before review/approval, or is a clean local devtools::check() sufficient?

Right now devtools::check() passes locally, but the PR shows some CI failures. If your expectation is CI-green, I will keep tightening the examples/tests (e.g., wrap JAGS bits in \dontrun{}/skip on CI, fix pkgdown) until everything passes.

@Kwan-Jenny please get all the CIs to pass. It's ok to use \dontrun{} on the JAGS-dependent examples. For the tests, they should probably pass on at least one operating system, but if you look at the run_mod() tests, you'll see that we skip on windows and linux because the results are OS-dependent.

Replace crossing(visit_num,time_days) with visit_num + lookup; rows now n_id*n_blocks*length(time_grid).
- Add `.runjags_fun` arg (defaults to runjags::run.jags).

- Replace direct `runjags::run.jags(...)` call with `.runjags_fun(...)`.

- No runtime change; enables clean mocking in testthat/CI.
- Add `.runjags_fun` (default `runjags::run.jags`) so tests can fake JAGS.
- Replace brittle joins on attr names with clean lookup tables
  (Subnum→Iso_type, Subject index→Subject_label), then rename to Subject.
- Use `select(-any_of(...))` / `rename(... = any_of(...))` to avoid NSE notes.
- Output schema unchanged; enables CI-safe unit test of run_mod_kron.
- Keep `used_priors` attribute while dropping legacy fields (`omega`, `wishdf`, `Omega`, `WishDF`, `prec.par`).

- Adjust test to assert attribute preservation and retain existing checks.
…od_kron

- Add correlated flag to run_mod to toggle Kronecker (Chapter 2) model.
- New args: file_mod_kron (path to CH2 JAGS file).
- When correlated=TRUE:
Clean base priors via clean_priors().
Add Kronecker priors from prep_priors_multi_b() and n_blocks.
Switch inits to inits_kron() and monitor TauB, TauP.
Auto-write fallback model_ch2_kron.jags if missing.

- Preserve original behavior when correlated=FALSE.
- Remove deprecated run_mod_kron(); functionality folded into run_mod.
@Kwan-Jenny
Copy link
Collaborator Author

Kwan-Jenny commented Oct 9, 2025

Hello @d-morrison ,

Question: Kronecker model in JAGS—choose Option A (fix TauP) or Option B (switch backend)?

I implemented the Chapter-2 “correlated biomarkers” model using a Kronecker precision
Tau = TauB ⊗ TauP, with both TauP ~ Wishart(OmegaP, nuP) and TauB ~ Wishart(OmegaB, nuB).

When I try to run it in JAGS, I get:

Error in node TauB: Unable to find appropriate sampler

This happens because JAGS needs conjugacy to sample Wishart nodes. Once Tau is the product of two stochastic matrices, neither TauP nor TauB is conjugate, so JAGS can’t assign a sampler.

Two paths forward:

  1. Option 1 (stay in JAGS, minimal change):
    Treat TauP as fixed data (TauP_data) and only estimate TauB ~ Wishart. That preserves cross-biomarker correlation (our Chapter-2 goal), while keeping within-biomarker covariance fixed. I can expose TauP_data as an argument (defaulting to diag(5) or a weakly informative precision). This will run in JAGS.

  2. Option 2 (full model, different backend):
    Move the Kronecker model (both TauP and TauB stochastic) to STAN (LKJ/Cholesky) or NIMBLE (custom samplers). This is the statistically cleaner version but requires a new backend and some refactoring.

What do you think is the best?

  • If Option 1, I will modify write_model_ch2_kron() to accept TauP_data and remove TauP ~ Wishart from the JAGS file, then pass TauP_data from run_mod(correlated=TRUE, ...).

  • If Option 2, I will prototype a STAN/NIMBLE version and keep the JAGS path for Chapter-1 (independence) only.

I added an example in the docs (wrapped in \dontrun{}) that reproduces the exact JAGS error so reviewers can see it locally.

@Kwan-Jenny
Copy link
Collaborator Author

Hi @d-morrison ,

What I changed (short version)

  • I kept the old independent model as the default.

  • I added an option correlated = TRUE to run_mod() to fit a Chapter 2 model that allows correlation across biomarkers using a Kronecker structure.

  • To support that, I added small helper functions and one JAGS model writer. These are only used when correlated = TRUE. They do not break the old code.


New files (what they do and why)

  1. simulate_multi_b_long()

Role: make fake data with multiple biomarkers that are correlated.
Why: so we can quickly test the Chapter 2 ideas end-to-end (simulate → fit).
What it does:

  • Draws subject-level latent parameters from a Kronecker covariance: Σ_total = Σ_B ⊗ Σ_P (5 params per biomarker).
  • Generates noisy observations on a time grid using our existing ab() curve.
    Output: a long tibble (Subject, visit_num, antigen_iso, time_days, value) plus a truth list with the matrices used.
  1. prep_priors_multi_b()

Role: provide priors for the Kronecker precision pieces.
Why: the correlated model needs two Wisharts: one within-biomarker (TauP, 5x5) and one across-biomarker (TauB, BxB).
What it returns: a small list {OmegaP, nuP, OmegaB, nuB} with weakly-informative defaults (proper but loose).

  1. clean_priors()

Role: remove legacy fields from the usual priors when we switch to the Kronecker prior.
Why: the independent model uses per-biomarker prec.par and its Wishart hyper-prior (omega, wishdf). Those conflict with the Kronecker prior.
What it does: drop omega, wishdf, prec.par, etc., but keep the used_priors attribute.

  1. inits_kron()

Role: safe initial values for the correlated model.
Why: if TauB, TauP, or prec.par appear in user inits, JAGS can complain.
What it does: returns RNG seeds and explicitly NULLs TauB, TauP, and prec.par so JAGS can initialize them itself.

  1. write_model_ch2_kron()

Role: write the JAGS model file for Chapter 2 to disk.
Why: the prior is different: draw all biomarkers at once per subject using Kronecker precision: Tau = TauB ⊗ TauP.
What is unchanged: transforms to natural scale and the likelihood for log-antibody values. Only the prior on the stacked parameters changed.

Note: JAGS struggles to find a sampler when we set both TauB ~ Wishart and TauP ~ Wishart (error: “Unable to find appropriate sampler for TauB”). I documented this in the examples and tests. We can discuss two paths: (1) fix by treating one precision as fixed/data or re-parameterize, or (2) move the correlated prior to STAN/NIMBLE.


Updated file

run_mod() (the main user function)

Goal: let users choose the old model (independent) or the new model (correlated) with one function.

What’s new (kept backward-compatible):

  • New arguments:

    • correlated = FALSE (default = old behavior)

    • file_mod_kron = "model_ch2_kron.jags" (path for Chapter 2 JAGS file)

  • Internal switch:

    • If correlated = FALSE: same as before (use prep_priors(), file_mod, initsfunction, monitor core params).

    • If correlated = TRUE:

      1. Build base priors with prep_priors() and then run clean_priors() to drop legacy fields.

      2. Add Kronecker priors with prep_priors_multi_b(n_blocks = n_antigen_isos) and pass scalar n_blocks.

      3. Choose the Chapter 2 model file: use file_mod_kron if it exists; otherwise auto-write one via write_model_ch2_kron(tempfile).

      4. Use inits_kron() for safe initialization.

      5. Monitor c("y0","y1","t1","alpha","shape","TauB","TauP").

What did not change: the output shape/class, attributes, and the fitted/residuals calculation stay the same, so downstream code should still work.


Why these helpers live as separate functions

  • They keep run_mod() simple and maintain backward compatibility.

  • They make the correlated path explicit and testable (each helper has its own example/tests).

  • If we later move to STAN/NIMBLE, these helpers can be reused (data prep and hyper-prior assembly stay the same).

@Kwan-Jenny Kwan-Jenny requested a review from d-morrison October 9, 2025 21:48
@d-morrison d-morrison requested a review from Copilot October 13, 2025 19:41
Copy link
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull Request Overview

This PR adds Chapter 2 functionality for fitting a correlated model using Kronecker precision matrices across biomarkers. The main purpose is to extend the existing independence model to allow correlations both within and across biomarkers using a Kronecker structure.

  • Implements fake data simulation with known covariance structure for testing
  • Adds Kronecker model infrastructure with new JAGS model file generation
  • Extends run_mod() with a correlated parameter to switch between independence and Kronecker models

Reviewed Changes

Copilot reviewed 30 out of 30 changed files in this pull request and generated 3 comments.

Show a summary per file
File Description
vignettes/articles/Chapter2.qmd Documentation and examples for Chapter 2 Kronecker model implementation
R/write_model_ch2_kron.R Generates JAGS model file with Kronecker precision structure
R/simulate_multi_b_long.R Simulates longitudinal data with Kronecker covariance
R/prep_priors_multi_b.R Sets up Wishart hyperparameters for Kronecker prior
R/inits_kron.R Provides safe initialization for Kronecker model
R/clean_priors.R Removes legacy prior fields to avoid conflicts
R/Run_Mod.R Extends run_mod() to support correlated model fitting
tests/testthat/*.R Unit tests for new functions
man/*.Rd Documentation files for new functions
inst/examples/*.R Example code for new functions

Tip: Customize your code reviews with copilot-instructions.md. Create the file or learn how to get started.

#' Note: These objects can be large.
#' @param correlated Logical; use Chapter-2 Kronecker prior across biomarkers.
#' Default FALSE (independence).
#' @param file_mod_kron Path to a JAGS file for the Kronecker model.If
Copy link

Copilot AI Oct 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing space after period in 'model.If' - should be 'model. If'.

Suggested change
#' @param file_mod_kron Path to a JAGS file for the Kronecker model.If
#' @param file_mod_kron Path to a JAGS file for the Kronecker model. If

Copilot uses AI. Check for mistakes.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Kwan-Jenny please consider this suggestion and explain your decision one way or the other?

Comment on lines +37 to +39
cli::cli_abort("`omega_b_scale` must be a numeric vector of length
`n_blocks`.")
}
Copy link

Copilot AI Oct 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The error message is split across lines with inconsistent indentation. Consider using a single line or proper string concatenation for better readability.

Suggested change
cli::cli_abort("`omega_b_scale` must be a numeric vector of length
`n_blocks`.")
}
cli::cli_abort("`omega_b_scale` must be a numeric vector of length `n_blocks`.")
}
}

Copilot uses AI. Check for mistakes.
# Step 2: write the new Kronecker model file (once per session/project)
write_model_ch2_kron()

# Step 3: run the unified runner in correlated (Chapter 2) mode
Copy link

Copilot AI Oct 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The comment mentions 'unified runner' but the function being called is run_mod(). This could be confusing since there's no mention of what makes it 'unified'. Consider clarifying or removing this term.

Suggested change
# Step 3: run the unified runner in correlated (Chapter 2) mode
# Step 3: run the model in correlated (Chapter 2) mode using run_mod()

Copilot uses AI. Check for mistakes.
#' - An optional `"jags.post"` attribute, included when argument
#' `with_post` = TRUE.
#' @inheritDotParams prep_priors
#' @seealso clean_priors, prep_priors_multi_b, inits_kron, write_model_ch2_kron
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please cross-reference wherever possible:

Suggested change
#' @seealso clean_priors, prep_priors_multi_b, inits_kron, write_model_ch2_kron
#' @seealso [clean_priors], [prep_priors_multi_b], [inits_kron], [write_model_ch2_kron]

if (!correlated) {
# original (independence) behavior
priorspec <- prep_priors(
max_antigens = longdata$n_antigen_isos, ...
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comment on lines +149 to +156
base_priors <- serodynamics::clean_priors(base_priors)

kron_priors <- serodynamics::prep_priors_multi_b(
n_blocks = longdata$n_antigen_isos
)
B_scalar <- list(n_blocks = longdata$n_antigen_isos)

priorspec <- c(base_priors, kron_priors, B_scalar)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could we decompose this section into a single function that takes base_priors and longdata as arguments and returns priorspec?

Comment on lines +138 to +147
priorspec <- prep_priors(
max_antigens = longdata$n_antigen_isos, ...
)
model_path <- file_mod # UNCHANGED for independence
init_fun <- initsfunction
to_monitor <- c("y0", "y1", "t1", "alpha", "shape")
} else {
# CH2 behavior: Kronecker prior across biomarkers
base_priors <- prep_priors(
max_antigens = longdata$n_antigen_isos, ...
Copy link
Collaborator

@d-morrison d-morrison Nov 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

since

priorspec <- prep_priors(
        max_antigens = longdata$n_antigen_isos, ...
)

and

base_priors <- prep_priors(
        max_antigens = longdata$n_antigen_isos, ...
)

are the same function call, let's move them out of the if-else and combine them into a single call?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants