-
Notifications
You must be signed in to change notification settings - Fork 0
Add chapter2.qmd: fake-data simulation and correlated model setup #133
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Conversation
Codecov Report❌ Patch coverage is
|
|
📖 https://ucd-serg.github.io/serodynamics/preview/pr133 |
There was a problem hiding this 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.
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
Hi @d-morrison |
d-morrison
left a comment
There was a problem hiding this 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.
Hi, @d-morrison Done. I moved the remaining code from the qmd into What changed:
What I would like reviewed specifically:
Thanks! |
Hi @d-morrison , 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 Right now |
@Kwan-Jenny please get all the CIs to pass. It's ok to use |
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.
|
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 When I try to run it in JAGS, I get:
This happens because JAGS needs conjugacy to sample Wishart nodes. Once Two paths forward:
What do you think is the best?
I added an example in the docs (wrapped in |
|
Hi @d-morrison , What I changed (short version)
New files (what they do and why)
Role: make fake data with multiple biomarkers that are correlated.
Role: provide priors for the Kronecker precision pieces.
Role: remove legacy fields from the usual priors when we switch to the Kronecker prior.
Role: safe initial values for the correlated model.
Role: write the JAGS model file for Chapter 2 to disk. Note: JAGS struggles to find a sampler when we set both Updated file
Goal: let users choose the old model (independent) or the new model (correlated) with one function. What’s new (kept backward-compatible):
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
|
There was a problem hiding this 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 acorrelatedparameter 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 |
Copilot
AI
Oct 13, 2025
There was a problem hiding this comment.
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'.
| #' @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 |
There was a problem hiding this comment.
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?
| cli::cli_abort("`omega_b_scale` must be a numeric vector of length | ||
| `n_blocks`.") | ||
| } |
Copilot
AI
Oct 13, 2025
There was a problem hiding this comment.
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.
| 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`.") | |
| } | |
| } |
| # 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 |
Copilot
AI
Oct 13, 2025
There was a problem hiding this comment.
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.
| # Step 3: run the unified runner in correlated (Chapter 2) mode | |
| # Step 3: run the model in correlated (Chapter 2) mode using run_mod() |
| #' - 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 |
There was a problem hiding this comment.
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:
| #' @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, ... |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| max_antigens = longdata$n_antigen_isos, ... | |
| max_antigens = longdata$n_antigen_isos, | |
| ... |
| 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) |
There was a problem hiding this comment.
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?
| 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, ... |
There was a problem hiding this comment.
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?
This adds my Chapter 2 .qmd file.
It includes:
run_mod_kronfor running the modelI am mainly stuck on the
prep_priors_multiB()andrun_mod_kron()parts, and would like feedback on whether the priors and wrapper are set up correctly.