diff --git a/.Rbuildignore b/.Rbuildignore index 1478b93b..4de8f509 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -11,6 +11,9 @@ ^data-raw$ ^README\.Rmd$ ^\.lintr\.R$ +^\.lintr$ +^scripts$ +^slurm$ ^_quarto\.yml$ ^.*\.qmd$ ^.*\.png$ diff --git a/.do_commit.sh b/.do_commit.sh new file mode 100644 index 00000000..306a71b7 --- /dev/null +++ b/.do_commit.sh @@ -0,0 +1,12 @@ +#!/bin/bash +git -c user.name="Claude" -c user.email="claude[bot]@users.noreply.github.com" \ + commit --author="claude[bot] " \ + -m "fix(sampler): increase max_treedepth default from 12 to 15 + +50% treedepth saturation on n=5 phase0 run (500/1000 transitions hit +the ceiling) is the dominant driver of ESS_bulk=16 and R-hat=1.766. +The LKJ Cholesky geometry requires longer HMC trajectories than +max_treedepth=12 allows; raising to 15 gives 8x more trajectory +length budget before truncation. + +Co-Authored-By: Claude Sonnet 4.6 " diff --git a/.gitconfig_tmp b/.gitconfig_tmp new file mode 100644 index 00000000..0ed23cd6 --- /dev/null +++ b/.gitconfig_tmp @@ -0,0 +1,3 @@ +[user] + name = Claude + email = claude[bot]@users.noreply.github.com diff --git a/.github/workflows/phase0-debug.yaml b/.github/workflows/phase0-debug.yaml index dbf525e0..ed656e8d 100644 --- a/.github/workflows/phase0-debug.yaml +++ b/.github/workflows/phase0-debug.yaml @@ -8,6 +8,16 @@ name: Phase 0 Debug Loop # to the branch under outputs/ci/ so Claude can read them in the # next turn. # +# Architecture note: +# This workflow uses the standard r-lib/actions/setup-r-dependencies@v2 +# pattern, matching every other R workflow in this repo (R-CMD-check, +# test-coverage, pkgdown, etc.). That action handles DESCRIPTION +# parsing, Remotes: field resolution, system dependencies, caching, +# and (critically) GitHub authentication correctly. Earlier versions +# of this workflow tried to call remotes::install_deps() directly +# and hit 401 errors that turned out to be due to token plumbing +# that setup-r-dependencies handles transparently. +# # Trigger: # Manual only (workflow_dispatch). Use either the "Run workflow" # button in the Actions tab, or: @@ -49,67 +59,71 @@ permissions: contents: write concurrency: - # one phase0 run at a time per branch — prevents conflicting commits group: phase0-debug-${{ github.ref }} cancel-in-progress: false jobs: phase0: runs-on: ubuntu-latest - timeout-minutes: 350 # 350 min is safely under the 360 free-tier cap + timeout-minutes: 350 env: - N_SUBJECTS: ${{ inputs.n }} - ITER_WARMUP: ${{ inputs.iter_warmup }} - ITER_SAMPLING: ${{ inputs.iter_sampling }} - CHAINS: ${{ inputs.chains }} - CMDSTAN_VERSION: "2.38.0" + # Matches every other R workflow in this repo. Required for setup-r-dependencies + # to authenticate package downloads from GitHub (including UCD-SERG/serodynamics + # via the Remotes: field in DESCRIPTION). + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + N_SUBJECTS: ${{ inputs.n }} + ITER_WARMUP: ${{ inputs.iter_warmup }} + ITER_SAMPLING: ${{ inputs.iter_sampling }} + CHAINS: ${{ inputs.chains }} + CMDSTAN_VERSION: "2.38.0" steps: - name: Checkout PR branch uses: actions/checkout@v4 with: ref: ${{ github.ref }} - # persist-credentials default true; needed for the later push. - name: Setup R uses: r-lib/actions/setup-r@v2 with: - r-version: "release" use-public-rspm: true + - name: Install package dependencies + # This step handles: + # - DESCRIPTION Imports / Depends / LinkingTo / Suggests + # - Remotes: field (stan-dev/cmdstanr, UCD-SERG/serodynamics) + # - System libraries (apt-get install) + # - Caching (transparent) + # The same action and same auth setup are used by R-CMD-check.yaml. + uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: | + any::devtools + any::posterior + stan-dev/cmdstanr + needs: check + - name: Cache cmdstan id: cache-cmdstan uses: actions/cache@v4 with: path: ~/.cmdstan key: cmdstan-${{ runner.os }}-${{ env.CMDSTAN_VERSION }} - - - name: System dependencies for R packages - run: | - sudo apt-get update - sudo apt-get install -y libcurl4-openssl-dev libssl-dev libxml2-dev \ - libfontconfig1-dev libharfbuzz-dev libfribidi-dev \ - libfreetype6-dev libpng-dev libtiff5-dev libjpeg-dev - - - name: Install R packages (CRAN + r-universe) - run: | - install.packages(c("devtools", "remotes", "posterior", "rlang")) - install.packages("cmdstanr", - repos = c("https://stan-dev.r-universe.dev", - "https://cloud.r-project.org")) - shell: Rscript {0} + restore-keys: | + cmdstan-${{ runner.os }}- - name: Install cmdstan (only if cache miss) if: steps.cache-cmdstan.outputs.cache-hit != 'true' run: | cmdstanr::check_cmdstan_toolchain(fix = TRUE) - cmdstanr::install_cmdstan(version = Sys.getenv("CMDSTAN_VERSION"), - cores = 2) + cmdstanr::install_cmdstan( + version = Sys.getenv("CMDSTAN_VERSION"), + cores = 2 + ) shell: Rscript {0} - name: Register cmdstan path run: | - # Whether cached or freshly installed, point cmdstanr at it. paths <- list.files("~/.cmdstan", pattern = "^cmdstan-", full.names = TRUE) stopifnot(length(paths) >= 1) @@ -118,14 +132,14 @@ jobs: cat("cmdstan version:", cmdstanr::cmdstan_version(), "\n") shell: Rscript {0} - - name: Install shigella package + dependencies - run: | - # remotes handles the Remotes: field in DESCRIPTION - # (UCD-SERG/serodynamics, stan-dev/cmdstanr). - remotes::install_deps(".", dependencies = TRUE, upgrade = "never") - devtools::install(".", dependencies = FALSE, upgrade = "never", - quick = TRUE, build = FALSE) - shell: Rscript {0} + - name: Install local shigella package + # Use R CMD INSTALL (base R, shell-level) instead of devtools::install + # or remotes::install_local. setup-r-dependencies@v2 doesn't guarantee + # remotes/devtools remain in the user library after it finishes, and + # devtools::install's `upgrade` arg expectations vary by version. + # R CMD INSTALL has no such dependencies — it's the lowest-level + # installer and always works. + run: R CMD INSTALL --no-multiarch --with-keep.source . - name: Run Phase 0 diagnostic id: run @@ -137,7 +151,6 @@ jobs: mkdir -p "${OUT_DIR}" Rscript -e ' - # Explicit %||% for R < 4.4 compatibility on GitHub runners `%||%` <- function(a, b) if (is.null(a)) b else a out_dir <- Sys.getenv("OUT_DIR") @@ -165,14 +178,13 @@ jobs: NULL }) - # Write a structured SUMMARY.txt that Claude can grep. summary_path <- file.path(out_dir, "SUMMARY.txt") lines <- c( - sprintf("RUN_ID: %s", Sys.getenv("GITHUB_RUN_ID")), - sprintf("N: %d", n), + sprintf("RUN_ID: %s", Sys.getenv("GITHUB_RUN_ID")), + sprintf("N: %d", n), sprintf("ITER_WARMUP: %d", iter_warmup), sprintf("ITER_SAMPLING: %d", iter_samp), - sprintf("CHAINS: %d", chains) + sprintf("CHAINS: %d", chains) ) if (is.null(res)) { @@ -195,7 +207,6 @@ jobs: sprintf("TREEDEPTH: %d / %d", sum(d$num_max_treedepth %||% 0L), total_iters) ) - # Diagnostic verdict — Claude can read this directly. ess <- s$ess_bulk %||% 0 rhat <- s$rhat %||% Inf divg <- sum(d$num_divergent %||% 0L) / total_iters @@ -216,16 +227,40 @@ jobs: cat("\n=== End SUMMARY.txt ===\n") ' 2>&1 | tee "${OUT_DIR}/run.log" - # Always succeed at the workflow level — the SUMMARY.txt - # carries the pass/fail signal. exit 0 env: OUT_DIR: outputs/ci/phase0_n${{ inputs.n }}_run${{ github.run_id }} + + - name: Run bimodality diagnostic + # Runs after the fit so the new RDS is available. Reads its own + # OUT_DIR from the positional argument we pass. Skips with a + # warning (not a fail) if the script or RDS is missing. + if: always() + env: + OUT_DIR: outputs/ci/phase0_n${{ inputs.n }}_run${{ github.run_id }} + run: | + set +e + if [ ! -f scripts/diagnostic_bimodality.R ]; then + echo "diagnostic_bimodality.R missing; skipping." + exit 0 + fi + if [ ! -f "${OUT_DIR}/one_fit_n${N_SUBJECTS}_ci.rds" ]; then + echo "Fit RDS missing at ${OUT_DIR}; skipping diagnostic." + exit 0 + fi + echo "Running diagnostic_bimodality.R on ${OUT_DIR}..." + Rscript scripts/diagnostic_bimodality.R "${OUT_DIR}" 2>&1 \ + | tee "${OUT_DIR}/diagnostic_bimodality.log" + exit 0 - name: Commit outputs back to the branch if: always() run: | set -euo pipefail + if [ ! -d outputs/ci ] || [ -z "$(ls -A outputs/ci 2>/dev/null)" ]; then + echo "outputs/ci/ is empty or missing; nothing to commit." + exit 0 + fi git config user.name "github-actions[bot]" git config user.email "github-actions[bot]@users.noreply.github.com" git add outputs/ci/ @@ -233,9 +268,7 @@ jobs: echo "No new outputs to commit." exit 0 fi - # [skip ci] prevents the push from triggering other workflows git commit -m "ci(phase0): n=${N_SUBJECTS} run ${GITHUB_RUN_ID} [skip ci]" - # Push to the branch we ran on BRANCH="${GITHUB_REF#refs/heads/}" git push origin "HEAD:${BRANCH}" diff --git a/.gitignore b/.gitignore index daf8d693..eb3d1bbe 100644 --- a/.gitignore +++ b/.gitignore @@ -18,3 +18,5 @@ README_files shigella.Rcheck/ shigella*.tar.gz shigella*.tgz +commit.sh +gitconfig diff --git a/.lintr.R b/.lintr.R index 500ecde0..14d9f3c3 100644 --- a/.lintr.R +++ b/.lintr.R @@ -1,9 +1,8 @@ - -undesirable_functions <- - lintr::default_undesirable_functions |> +undesirable_functions <- + lintr::default_undesirable_functions |> lintr::modify_defaults( - + # following https://github.com/r-lib/devtools/blob/2aa51ef/.lintr.R: # Base messaging "message" = "use cli::cli_inform()", @@ -18,32 +17,33 @@ undesirable_functions <- "cli_alert_info" = "use cli::cli_inform()", "cli_alert_success" = "use cli::cli_inform()", "cli_alert_warning" = "use cli::cli_inform()", - + library = paste( "\nuse `::`, `usethis::use_import_from()`, or `withr::local_package()`", "instead of modifying the global search path.", "\nSee:\n", " and\n", "", - "\nfor more details" + "\nfor more details." ), - - structure = NULL - # see https://github.com/r-lib/lintr/pull/2227 and + + structure = NULL, + browser = NULL + # see https://github.com/r-lib/lintr/pull/2227 and # rebuttal https://github.com/r-lib/lintr/pull/2227#issuecomment-1800302675 - + ) # define snake_case with uppercase acronyms allowed; # see https://github.com/r-lib/lintr/issues/2844 for details: withr::local_package("rex") -snake_case_ACRO = rex::rex( +snake_case_ACROs1 <- rex::rex( start, maybe("."), - some_of(lower, digit) %or% some_of(upper, digit), + list(some_of(upper), maybe("s"), zero_or_more(digit)) %or% list(some_of(lower), zero_or_more(digit)), zero_or_more( "_", - some_of(lower, digit) %or% some_of(upper, digit) + list(some_of(upper), maybe("s"), zero_or_more(digit)) %or% list(some_of(lower), zero_or_more(digit)) ), end ) @@ -54,7 +54,7 @@ linters <- lintr::linters_with_defaults( lintr::redundant_equals_linter(), lintr::pipe_consistency_linter(pipe = "|>"), lintr::object_name_linter( - regexes = c(snake_case_ACRO = snake_case_ACRO) + regexes = c(snake_case_ACROs1 = snake_case_ACROs1) ), lintr::undesirable_function_linter( fun = undesirable_functions, @@ -64,9 +64,22 @@ linters <- lintr::linters_with_defaults( # prevent warnings from lintr::read_settings: rm(undesirable_functions) -rm(snake_case_ACRO) +rm(snake_case_ACROs1) exclusions <- list( `data-raw` = list( - pipe_consistency_linter = Inf - ) + pipe_consistency_linter = Inf, + undesirable_function_linter = Inf + ), + vignettes = list( + undesirable_function_linter = Inf, + object_name_linter = Inf + ), + "inst/examples" = list( + undesirable_function_linter = Inf + ), + "tests/testthat.R" = list( + undesirable_function_linter = Inf + ), + "quarto/mermaid-diagrams.qmd" = Inf + ) diff --git a/DESCRIPTION b/DESCRIPTION index 8287647d..5af11f59 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,23 +1,34 @@ Package: shigella -Title: What the Package Does (One Line, Title Case) -Version: 0.0.0.9008 +Title: Bayesian Modeling of Shigella Antibody Kinetics +Version: 0.0.0.9009 Authors@R: c( person("Kwan Ho", "Lee", , "ksjlee@ucdavis.edu", role = c("aut", "cre")), person("Douglas Ezra", "Morrison", , "demorrison@ucdavis.edu", role = c("aut"), comment = c(ORCID = "0000-0002-7195-830X"))) -Description: What the package does (one paragraph). +Description: Tools for multivariate Bayesian hierarchical modeling of + antibody response trajectories following confirmed Shigella infection, + supporting kinetic parameter estimation and serosurveillance + applications via Stan (cmdstanr) backends. License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) URL: https://ucd-serg.github.io/shigella/ Remotes: + stan-dev/cmdstanr, UCD-SERG/serodynamics Suggests: + cmdstanr, knitr, + posterior, rmarkdown, serodynamics (>= 0.0.0.9011), spelling, testthat (>= 3.0.0) +Imports: + cli, + dplyr, + MASS, + tibble VignetteBuilder: knitr Depends: R (>= 3.5) diff --git a/NAMESPACE b/NAMESPACE index 6ae92683..4453a68a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,2 +1,10 @@ # Generated by roxygen2: do not edit by hand +export(postprocess_stan_output) +export(prep_data_stan) +export(prep_priors_stan) +export(run_mod_stan) +export(sim_correlated_case_data) +export(write_status) +importFrom(stats,median) +importFrom(stats,rnorm) diff --git a/R/build_sigma_matrices.R b/R/build_sigma_matrices.R new file mode 100644 index 00000000..4e69d304 --- /dev/null +++ b/R/build_sigma_matrices.R @@ -0,0 +1,28 @@ +# Helper: build Sigma_P, Sigma_B, Sigma_eps, Sigma_full, and mu_vec. +# Returns a list with sigma_eps, sigma_full, mu_vec. +#' @keywords internal +#' @noRd +.build_sigma_matrices <- function(mu, tau_P, tau_B, tau_eps, + omega_P, omega_B, omega_eps) { + n_param <- length(mu) + n_biomarker <- length(tau_B) + + sigma_P <- diag(tau_P) %*% omega_P %*% diag(tau_P) + sigma_B <- diag(tau_B) %*% omega_B %*% diag(tau_B) + sigma_eps <- diag(tau_eps) %*% omega_eps %*% diag(tau_eps) + + # Sigma_full = Sigma_B kron Sigma_P, dimension PK x PK + sigma_full <- kronecker(sigma_B, sigma_P) + + # vec(M) where M is P x K (columns = biomarkers) + # Assume same mu for all biomarkers (can be extended) + mean_matrix <- matrix( + mu, + nrow = n_param, + ncol = n_biomarker, + byrow = FALSE + ) + mu_vec <- as.vector(mean_matrix) + + list(sigma_eps = sigma_eps, sigma_full = sigma_full, mu_vec = mu_vec) +} diff --git a/R/case_data_to_prepped_jags.R b/R/case_data_to_prepped_jags.R new file mode 100644 index 00000000..3ac6b2da --- /dev/null +++ b/R/case_data_to_prepped_jags.R @@ -0,0 +1,15 @@ +# Helper: convert a case_data object to prepped_jags_data via serodynamics. +#' @keywords internal +#' @noRd +.case_data_to_prepped_jags <- function(data) { + if (!requireNamespace("serodynamics", quietly = TRUE)) { + cli::cli_abort(c( + "Package {.pkg serodynamics} is required.", + "i" = "Install it before using {.fn prep_data_stan}." + )) + } + serodynamics::prep_data( + data, + add_newperson = FALSE + ) +} diff --git a/R/compute_kinetics_at_time.R b/R/compute_kinetics_at_time.R new file mode 100644 index 00000000..d8ea59c8 --- /dev/null +++ b/R/compute_kinetics_at_time.R @@ -0,0 +1,42 @@ +# Helper: compute log-scale mean for one biomarker at one time point. +# Implements the two-phase power-law kinetics model. +#' @keywords internal +#' @noRd +.compute_kinetics_at_time <- function( + log_y0, log_y1m0, # nolint: object_name_linter. Stan param name. + log_t1, log_alpha, log_rm1, tt) { + y0 <- exp(log_y0) + y1 <- y0 + exp(log_y1m0) + t1_j <- exp(log_t1) + alpha <- exp(log_alpha) + shape <- exp(log_rm1) + 1 + + if (tt <= t1_j) { + beta_growth <- (log(y1) - log(y0)) / t1_j + return(log(y0) + beta_growth * tt) + } + + # Tolerance is sqrt(.Machine$double.eps) (~1.5e-8): well outside the + # production prior on shape (typical |1 - shape| > 0.3), but excludes + # the numerically-unstable region where log(term) / (1 - shape) loses + # meaningful precision. + if (abs(1 - shape) < sqrt(.Machine$double.eps)) { + cli::cli_abort(c( + "shape ~= 1 is degenerate for the two-phase model", + "i" = "log_rm1 = {log_rm1} produces shape = {shape}, |1 - shape| = + {abs(1 - shape)}", + "i" = paste0("the decay-phase formula log(term) / (1 - shape) is", + " undefined or numerically unstable in this region.") + )) + } + + term <- y1^(1 - shape) - (1 - shape) * alpha * (tt - t1_j) + if (term <= 0) { + cli::cli_abort(c( + "Trajectory infeasibility detected at t = {tt}", + "i" = paste0("term = {term} <= 0; parameter combination", + " is invalid for the two-phase model") + )) + } + log(term) / (1 - shape) +} diff --git a/R/compute_log_mu_k.R b/R/compute_log_mu_k.R new file mode 100644 index 00000000..14b6e201 --- /dev/null +++ b/R/compute_log_mu_k.R @@ -0,0 +1,20 @@ +# Helper: compute log mu for each biomarker for subject i at time tt. +# Returns log_mu_k (length-K numeric vector). +#' @keywords internal +#' @noRd +.compute_log_mu_k <- function(theta_arr, i, n_biomarker, tt) { + vapply( + seq_len(n_biomarker), + function(j) { + .compute_kinetics_at_time( # nolint: object_usage_linter + theta_arr[i, 1, j], + theta_arr[i, 2, j], + theta_arr[i, 3, j], + theta_arr[i, 4, j], + theta_arr[i, 5, j], + tt + ) + }, + numeric(1L) + ) +} diff --git a/R/draw_subject_params.R b/R/draw_subject_params.R new file mode 100644 index 00000000..b8d7d82e --- /dev/null +++ b/R/draw_subject_params.R @@ -0,0 +1,30 @@ +# Helper: draw subject parameters from the Kronecker covariance. +# Returns theta_arr (N x P x K). +#' @keywords internal +#' @noRd +.draw_subject_params <- function(n, mu_vec, sigma_full, + n_param, n_biomarker, antigen_isos) { + theta_vec <- MASS::mvrnorm(n = n, mu = mu_vec, Sigma = sigma_full) + if (is.null(dim(theta_vec))) { + theta_vec <- matrix(theta_vec, nrow = 1L) + } + + # dim n x PK; reshape to N x P x K + theta_arr <- array(NA_real_, dim = c(n, n_param, n_biomarker)) + + for (i in seq_len(n)) { + theta_arr[i, , ] <- matrix( + theta_vec[i, ], + nrow = n_param, + ncol = n_biomarker + ) + } + + dimnames(theta_arr) <- list( + subject = as.character(seq_len(n)), + param = c("log_y0", "log_y1m0", "log_t1", "log_alpha", "log_rm1"), + biomarker = antigen_isos + ) + + theta_arr +} diff --git a/R/generate_obs_for_subject.R b/R/generate_obs_for_subject.R new file mode 100644 index 00000000..d5dc4303 --- /dev/null +++ b/R/generate_obs_for_subject.R @@ -0,0 +1,37 @@ +# Helper: generate all observation rows for one subject. +# Returns a flat list of data frames (one per biomarker x time point). +#' @keywords internal +#' @noRd +.generate_obs_for_subject <- function(i, n_obs_per_subject, time_grid, + n_biomarker, theta_arr, antigen_isos, + u_eps) { + obs_times <- sort( + sample(time_grid, size = n_obs_per_subject, replace = FALSE) + ) + + rows <- list() + row_counter <- 1L + + for (tt_idx in seq_along(obs_times)) { + tt <- obs_times[tt_idx] + log_mu_k <- .compute_log_mu_k(theta_arr, i, n_biomarker, tt) + z <- rnorm(n_biomarker) + # u_eps is upper-triangular Cholesky from R's chol(); + # t(u_eps) %*% z gives MVN(0, sigma_eps) draws. + log_y_obs <- log_mu_k + as.vector(t(u_eps) %*% z) + + for (j in seq_len(n_biomarker)) { + rows[[row_counter]] <- data.frame( + id = as.character(i), + visit_num = tt_idx, + timeindays = tt, + antigen_iso = antigen_isos[j], + value = exp(log_y_obs[j]), + stringsAsFactors = FALSE + ) + row_counter <- row_counter + 1L + } + } + + rows +} diff --git a/R/generate_obs_rows.R b/R/generate_obs_rows.R new file mode 100644 index 00000000..d29dca6e --- /dev/null +++ b/R/generate_obs_rows.R @@ -0,0 +1,15 @@ +# Helper: generate all observation rows for all subjects. +# Returns a list of data frames (to be dplyr::bind_rows'd). +#' @keywords internal +#' @noRd +.generate_obs_rows <- function(n, n_obs_per_subject, time_grid, + n_biomarker, theta_arr, antigen_isos, u_eps) { + rows <- list() + for (i in seq_len(n)) { + rows <- c(rows, .generate_obs_for_subject( + i, n_obs_per_subject, time_grid, + n_biomarker, theta_arr, antigen_isos, u_eps + )) + } + rows +} diff --git a/R/locate_stan_file.R b/R/locate_stan_file.R new file mode 100644 index 00000000..498ef1ca --- /dev/null +++ b/R/locate_stan_file.R @@ -0,0 +1,35 @@ +# Helper: locate the Stan source file for the given model. +#' @keywords internal +#' @noRd +.locate_stan_file <- function(model, stan_dir) { + stan_basename <- paste0(model, ".stan") + + if (is.null(stan_dir)) { + stan_file <- system.file( + "stan", + stan_basename, + package = "shigella", + mustWork = FALSE + ) + + # Fallback for interactive development before the package is installed. + if (identical(stan_file, "") || !file.exists(stan_file)) { + stan_file <- file.path("inst", "stan", stan_basename) + } + } else { + stan_file <- file.path(stan_dir, stan_basename) + } + + if (!file.exists(stan_file)) { + cli::cli_abort(c( + "Cannot locate Stan file: {.file {stan_file}}", + "i" = "Working directory is: {.path {getwd()}}", + "i" = "If running interactively, check that + {.file inst/stan/{stan_basename}} exists.", + "i" = "If running from an installed package, use system.file('stan', + '{stan_basename}', package = 'shigella')." + )) + } + + stan_file +} diff --git a/R/postprocess_stan_output.R b/R/postprocess_stan_output.R new file mode 100644 index 00000000..f599e163 --- /dev/null +++ b/R/postprocess_stan_output.R @@ -0,0 +1,244 @@ +# ─── Internal extraction helpers ───────────────────────────────────────────── +# Consolidated here (from separate extract_*.R files) to eliminate cross-file +# object_usage_linter false positives on dotted-prefix internal helper calls. + +# Helper: extract all chain draws for one parameter cell `pname[subj, iso]`. +#' @keywords internal +#' @noRd +.extract_draws_for_cell <- function(col_name, draws_df, pname, + ids, antigens, stratification, n_chain) { + m <- regmatches(col_name, regexec("\\[(\\d+),(\\d+)\\]", + col_name))[[1]] + subj_idx <- as.integer(m[2]) + iso_idx <- as.integer(m[3]) + sub_df <- draws_df[, c(".chain", ".iteration", col_name)] + + lapply(seq_len(n_chain), function(ch) { + chain_data <- sub_df[sub_df$.chain == ch, ] + tibble::tibble( + Iteration = chain_data$.iteration, + Chain = ch, + Parameter = pname, + Iso_type = antigens[iso_idx], + Stratification = stratification, + Subject = ids[subj_idx], + value = chain_data[[col_name]] + ) + }) +} + +# Helper: extract draws for all param_names into a single tibble. +#' @keywords internal +#' @noRd +.extract_param_draws <- function(param_names, + draws_df, + N, + K, + ids, + antigens, + stratification, + n_chain) { + draws_per_param <- lapply(param_names, function(pname) { + matching_cols <- grep(paste0("^", pname, "\\["), colnames(draws_df), + value = TRUE) + if (length(matching_cols) != N * K) { + cli::cli_abort( + "Expected {N * K} {.var {pname}} draws; got {length(matching_cols)}." + ) + } + draws_per_col <- lapply( + matching_cols, .extract_draws_for_cell, + draws_df = draws_df, + pname = pname, + ids = ids, + antigens = antigens, + stratification = stratification, + n_chain = n_chain + ) + dplyr::bind_rows(unlist(draws_per_col, recursive = FALSE)) + }) + + dplyr::bind_rows(draws_per_param) +} + +# Helper: extract Omega_eps and Sigma_eps from a cmdstanr fit. Model 2 only. +#' @keywords internal +#' @noRd +.extract_residual_cov_stan <- function(stan_fit, K, antigens) { + tryCatch({ + omega_eps_arr <- posterior::as_draws_array( + stan_fit$draws(variables = "Omega_eps") + ) + sigma_eps_arr <- posterior::as_draws_array( + stan_fit$draws(variables = "Sigma_eps") + ) + omega_eps <- .summarize_matrix_draws(omega_eps_arr, "Omega_eps", K, K) + sigma_eps <- .summarize_matrix_draws(sigma_eps_arr, "Sigma_eps", K, K) + dimnames(omega_eps) <- list(antigens, antigens) + dimnames(sigma_eps) <- list(antigens, antigens) + list(Omega_eps = omega_eps, Sigma_eps = sigma_eps) + }, error = function(e) { + cli::cli_warn("Omega_eps/Sigma_eps not extracted: {e$message}") + list() + }) +} + +# Helper: extract Omega_B, Sigma_B, Omega_P, Sigma_P from a cmdstanr fit. +# Model 2 only. +#' @keywords internal +#' @noRd +.extract_kron_matrices <- function(stan_fit, K, param_names, antigens) { + tryCatch({ + omega_B_arr <- posterior::as_draws_array( + stan_fit$draws(variables = "Omega_B") + ) + sigma_B_arr <- posterior::as_draws_array( + stan_fit$draws(variables = "Sigma_B") + ) + omega_P_arr <- posterior::as_draws_array( + stan_fit$draws(variables = "Omega_P") + ) + sigma_P_arr <- posterior::as_draws_array( + stan_fit$draws(variables = "Sigma_P") + ) + + omega_B <- .summarize_matrix_draws(omega_B_arr, "Omega_B", K, K) + sigma_B <- .summarize_matrix_draws(sigma_B_arr, "Sigma_B", K, K) + P <- length(param_names) + omega_P <- .summarize_matrix_draws(omega_P_arr, "Omega_P", P, P) + sigma_P <- .summarize_matrix_draws(sigma_P_arr, "Sigma_P", P, P) + + dimnames(omega_B) <- list(antigens, antigens) + dimnames(sigma_B) <- list(antigens, antigens) + dimnames(omega_P) <- list(param_names, param_names) + dimnames(sigma_P) <- list(param_names, param_names) + + list(Omega_B = omega_B, Sigma_B = sigma_B, + Omega_P = omega_P, Sigma_P = sigma_P) + }, error = function(e) { + cli::cli_warn("Kronecker matrices not extracted: {e$message}") + list() + }) +} + +# Helper: extract per-biomarker Omega_P list from a model_1 cmdstanr fit. +# model_1 generates array[K] corr_matrix[P] Omega_P; cmdstanr names +# cells Omega_P[k,p,q]. Returns a named list of K matrices. +#' @keywords internal +#' @noRd +.extract_model1_omega_p_stan <- function(stan_fit, K, param_names, antigens) { + P <- length(param_names) + tryCatch({ + omega_P_arr <- posterior::as_draws_array( + stan_fit$draws(variables = "Omega_P") + ) + omega_P_list <- .summarize_matrix_array(omega_P_arr, "Omega_P", K, P, P) + for (k in seq_len(K)) { + dimnames(omega_P_list[[k]]) <- list(param_names, param_names) + } + names(omega_P_list) <- antigens + list(Omega_P = omega_P_list) + }, error = function(e) { + cli::cli_warn("model_1 Omega_P not extracted: {e$message}") + list() + }) +} + +# Helper: extract the log_lik draws matrix from a cmdstanr fit. +#' @keywords internal +#' @noRd +.extract_log_lik_stan <- function(stan_fit) { + tryCatch({ + list(log_lik = posterior::as_draws_matrix( + stan_fit$draws(variables = "log_lik") + )) + }, error = function(e) { + cli::cli_warn("log_lik not extracted: {e$message}") + list() + }) +} + +# ─── Main function ──────────────────────────────────────────────────────────── + +#' @title Post-process Stan output to sr_model format (cmdstanr version) +#' @description +#' Converts a CmdStanMCMC object (from cmdstanr's mod$sample()) into the +#' long-format tibble produced by run_mod(), so downstream plotting/summary +#' functions work without modification. +#' +#' @param stan_fit CmdStanMCMC object from cmdstanr (Not rstan stanfit) +#' @param ids subject IDs from attr(stan_data, "ids") +#' @param antigens biomarker names from attr(stan_data, "antigens") +#' @param model "model_1", "model_2" +#' @param stratification label for this stratum +#' @param param_names Character vector of parameter names to extract from the +#' Stan model's \code{generated quantities} block. Must match the variable +#' names exactly as declared in both \file{inst/stan/model_1.stan} and +#' \file{inst/stan/model_2.stan}. Defaults to +#' \code{c("y0", "y1", "t1", "alpha", "shape")}. +#' @return list with sr_tibble and cov_summaries +#' @example inst/examples/postprocess_stan_output-examples.R +#' @export +postprocess_stan_output <- function(stan_fit, + ids, + antigens, + model = c("model_2", "model_1"), + stratification = "None", + param_names = c("y0", "y1", "t1", + "alpha", "shape")) { + + model <- match.arg(model) + has_kron <- identical(model, "model_2") + + if (!requireNamespace("posterior", quietly = TRUE)) { + cli::cli_abort("Package {.pkg posterior} required for cmdstanr + postprocessing.") + } + + avail <- stan_fit$metadata()$stan_variables + missing_vars <- setdiff(param_names, avail) + if (length(missing_vars) > 0) { + cli::cli_abort(c( + "{.arg param_names} contains variables absent from the Stan model.", + "x" = "Missing: {.val {missing_vars}}", + "i" = "Available Stan variables: {.val {avail}}" + )) + } + + N <- length(ids) + K <- length(antigens) + + draws_df <- tibble::as_tibble(posterior::as_draws_df( + stan_fit$draws(variables = param_names) + )) + n_chain <- max(draws_df$.chain) + + sr_tibble <- .extract_param_draws( + param_names = param_names, + draws_df = draws_df, + N = N, + K = K, + ids = ids, + antigens = antigens, + stratification = stratification, + n_chain = n_chain + ) + + if (has_kron) { + cov_summaries <- c( + .extract_residual_cov_stan(stan_fit, K, antigens), + .extract_kron_matrices(stan_fit, K, param_names, antigens) + ) + } else { + cov_summaries <- .extract_model1_omega_p_stan( + stan_fit, K, param_names, antigens + ) + } + + cov_summaries <- c(cov_summaries, .extract_log_lik_stan(stan_fit)) + + return(list( + sr_tibble = sr_tibble, + cov_summaries = cov_summaries + )) +} diff --git a/R/prep_data_stan.R b/R/prep_data_stan.R new file mode 100644 index 00000000..e3807760 --- /dev/null +++ b/R/prep_data_stan.R @@ -0,0 +1,107 @@ +#' @title Prepare data for the Stan backend +#' @description +#' Converts case data into the structured list that the Stan models +#' (`model_1.stan`, `model_2.stan`) expect. Accepts either a raw +#' `case_data` object (the typical entry point) or a +#' `prepped_jags_data` object already produced by +#' [serodynamics::prep_data()] (the JAGS-side prep step). +#' +#' When given a `case_data` object, this function internally calls +#' [serodynamics::prep_data()] with `add_newperson = FALSE` (Stan +#' handles posterior prediction in the `generated quantities` block +#' rather than via a dummy missing-data subject). +#' +#' The Stan models expect: +#' - `N`: number of subjects +#' - `K`: number of antigen-isotype biomarkers +#' - `P`: number of kinetic parameters (always 5) +#' - `max_obs`: max number of observations per subject +#' - `n_obs[N]`: actual number of observations per subject +#' - `time_obs[N, max_obs]`: observation times (NA -> 0, ignored by +#' the likelihood via the `n_obs[i]` guard) +#' - `log_y[N, max_obs, K]`: log-transformed antibody observations +#' +#' @param data either a `case_data` object (output of +#' [sim_correlated_case_data()] or [serodynamics::as_case_data()]) +#' or a `prepped_jags_data` object (output of +#' [serodynamics::prep_data()]). +#' @param drop_newperson [logical] whether to drop the JAGS dummy +#' "newperson" row if it is present in a `prepped_jags_data` +#' input. Default `TRUE`. Has no effect when `data` is a +#' `case_data` object because the internal `prep_data()` call uses +#' `add_newperson = FALSE`. +#' +#' @returns a named [list] with attributes `ids` and `antigens`, +#' ready to pass to a compiled Stan model via +#' `cmdstanr::cmdstan_model()$sample()`. +#' @export +#' @example inst/examples/prep_data_stan-examples.R +prep_data_stan <- function(data, + drop_newperson = TRUE) { + + # Route prepped_jags_data directly; convert case_data via helper + if (inherits(data, "prepped_jags_data")) { + prepped_jags_data <- data + } else if (inherits(data, "case_data")) { + prepped_jags_data <- .case_data_to_prepped_jags(data) + } else { + cli::cli_abort(c( + "{.arg data} must be a {.cls case_data} or {.cls prepped_jags_data} + object", "i" = "Got an object of class {.cls {class(data)}}." + )) + } + + # Extract arrays + smpl_t <- prepped_jags_data$smpl.t # [nsubj, max_visits] + logy <- prepped_jags_data$logy # [nsubj, max_visits, K] + nsmpl <- as.integer(prepped_jags_data$nsmpl) + K <- prepped_jags_data$n_antigen_isos + N_full <- prepped_jags_data$nsubj + ids_all <- attr(prepped_jags_data, "ids") + + # Drop the "newperson" dummy row if present + if (drop_newperson && "newperson" %in% ids_all) { + keep_idx <- which(ids_all != "newperson") + smpl_t <- smpl_t[keep_idx, , drop = FALSE] + logy <- logy[keep_idx, , , drop = FALSE] + nsmpl <- nsmpl[keep_idx] + ids_kept <- ids_all[keep_idx] + N <- length(keep_idx) + } else { + ids_kept <- ids_all + N <- N_full + } + + max_obs <- ncol(smpl_t) + # P = 5 is dictated by the kinetic parameter layout in inst/stan/model_2.stan + # (data block): log_y0, log_y1m0, log_t1, log_alpha, log_rm1. + # It cannot be derived from prepped_jags_data dimensions and must stay in sync + # with the Stan model manually. + P <- 5L + + # Replace NA with 0; Stan ignores these via the n_obs[i] guard in + # the likelihood loop (for (t_idx in 1:n_obs[i])). + time_obs <- smpl_t + time_obs[is.na(time_obs)] <- 0 + log_y <- logy + log_y[is.na(log_y)] <- 0 + + # Sanity checks + .validate_stan_arrays(nsmpl, max_obs) + + antigens <- attr(prepped_jags_data, "antigens") + stan_data <- list( + N = N, + K = as.integer(K), + P = P, + max_obs = as.integer(max_obs), + n_obs = nsmpl, + time_obs = time_obs, + log_y = log_y + ) + + # Attach metadata for postprocessing + attr(stan_data, "ids") <- ids_kept + attr(stan_data, "antigens") <- antigens + return(stan_data) +} diff --git a/R/prep_priors_stan.R b/R/prep_priors_stan.R new file mode 100644 index 00000000..6d43f620 --- /dev/null +++ b/R/prep_priors_stan.R @@ -0,0 +1,74 @@ +#' @title Prepare priors for Stan backend +#' @description +#' Translates the JAGS prior specification into Stan's LKJ + half-Cauchy +#' decomposition. +#' +#' Defaults match the JAGS prior specification, with two adjustments for +#' Stan compatibility: +#' - mu_hyp_sd capped at ~5 (was 316 in JAGS) - Stan's HMC sampler +#' handles weakly-informative priors better with more reasonable scales. +#' JAGS Gibbs sampling tolerates wider priors, but Stan's gradient-based +#' HMC explores the tails too aggressively when sd is huge. +#' - tau scales = 1.0 (was 2.5) - keeps initial steps reasonable +#' +#' +#' @param mu_hyp_mean [numeric] length-5 prior mean for population params +#' @param mu_hyp_sd [numeric] length-5 prior SD for population params. +#' Weakly informative, Stan-friendly. JAGS original was c(1, 316, 1, 32, 1). +#' 5.0 on log-scale params covers ~5 orders of magnitude - plenty wide. +#' @param tau_P_scale half-Cauchy scale for parameter SDs +#' @param tau_B_scale half-Cauchy scale for biomarker SDs (Model 2 only) +#' @param tau_eps_scale half-Cauchy scale for residual SDs +#' @param lkj_P_eta LKJ shape for parameter correlation +#' @param lkj_B_eta LKJ shape for biomarker correlation (Model 2 only) +#' @param lkj_eps_eta LKJ shape for residual correlation +#' @param model character: "model_1", "model_2" +#' +#' @returns named list with priors for the Stan data block +#' @example inst/examples/prep_priors_stan-examples.R +#' @export +prep_priors_stan <- function( + mu_hyp_mean = c(1.0, 7.0, 1.0, -4.0, -1.0), + mu_hyp_sd = c(5.0, 5.0, 5.0, 5.0, 5.0), + tau_P_scale = 1.0, + tau_B_scale = 1.0, + tau_eps_scale = 1.0, + lkj_P_eta = 2.0, + lkj_B_eta = 4.0, + lkj_eps_eta = 2.0, + model = c("model_2", "model_1")) { + + model <- match.arg(model) + + has_kron <- identical(model, "model_2") + + if (length(mu_hyp_mean) != 5) { + cli::cli_abort("{.arg mu_hyp_mean} must be length 5.") + } + if (length(mu_hyp_sd) != 5) { + cli::cli_abort("{.arg mu_hyp_sd} must be length 5.") + } + + priors <- list( + mu_hyp_mean = mu_hyp_mean, + mu_hyp_sd = mu_hyp_sd, + tau_P_scale = tau_P_scale, + tau_eps_scale = tau_eps_scale, + lkj_P_eta = lkj_P_eta + ) + + if (has_kron) { + # model_2 uses an LKJ prior on the epsilon correlation matrix; + # model_1 does not. + priors$tau_B_scale <- tau_B_scale + priors$lkj_B_eta <- lkj_B_eta + priors$lkj_eps_eta <- lkj_eps_eta + } + + attr(priors, "model") <- model + # Snapshot of the list itself so run_mod_stan can attach it to the output + # without re-calling prep_priors_stan (self-referential by design). + attr(priors, "stan_input_snapshot") <- priors + class(priors) <- c("curve_params_priors_stan", "list") + return(priors) +} diff --git a/R/run_mod_stan.R b/R/run_mod_stan.R new file mode 100644 index 00000000..353689a0 --- /dev/null +++ b/R/run_mod_stan.R @@ -0,0 +1,172 @@ +#' @title Run Stan model using the cmdstanr backend +#' @description +#' Fits the two-phase antibody kinetics model using **cmdstanr**. +#' The `compile_dir` argument allows compiled Stan binaries to be written to +#' a writable directory (default: `/tmp`), which is useful on HPC systems +#' where the home directory is mounted non-executable. +#' +#' Output: an `sr_model` tibble with the same column schema as `run_mod()`, +#' so all existing plot / summary functions work unchanged. Stan-specific +#' attributes are also attached: +#' - `Omega_eps`, `Sigma_eps`: residual covariance (Model 2 only) +#' - `Omega_B`, `Sigma_B`: biomarker covariance (Model 2 only) +#' - `Omega_P`: parameter corr matrix (Model 2: single matrix; +#' Model 1: named list of K matrices, one per +#' biomarker) +#' - `Sigma_P`: parameter covariance (Model 2 only) +#' - `stan_fit`: raw CmdStanMCMC object (when with_post = TRUE) +#' +#' @param data case_data object (from sim_correlated_case_data() or +#' as_case_data()) +#' @param model character: "model_1", "model_2" +#' @param chains Number of chains to run. +#' @param strat optional stratification variable (default NA) +#' @param with_post return raw CmdStanMCMC object as attribute (default FALSE) +#' @param stan_dir Optional directory containing `model_*.stan` files. +#' If `NULL`, the function first looks for Stan files installed with the +#' package using `system.file("stan", ..., package = "shigella")`, then falls +#' back to `inst/stan` for interactive development. +#' @param compile_dir directory where cmdstanr writes compiled binaries. +#' Default uses STAN_COMPILE_DIR env var, or +#' /tmp//cmdstan_bin. +#' @param init initial value strategy. Numeric value scales down random init +#' (default 0.1 to avoid -inf in multi_normal_cholesky_lpdf) +#' @param ... additional priors passed to prep_priors_stan() +#' @param iter_sampling Number of post-warmup iterations per chain. +#' @param iter_warmup Number of warmup iterations per chain. +#' @param adapt_delta Target average acceptance probability for Stan sampling. +#' @param max_treedepth Maximum tree depth for Stan NUTS sampling. +#' @param seed Random seed passed to Stan. +#' @param parallel_chains Number of chains to run in parallel. +#' @param refresh Stan progress refresh interval. +#' @param show_messages Logical; whether to show CmdStan messages. +#' +#' @returns sr_model tibble +#' @example inst/examples/run_mod_stan-examples.R +#' @export +run_mod_stan <- function(data, + model = c("model_2", "model_1"), + chains = 4, + iter_sampling = 1000, + iter_warmup = 1000, + adapt_delta = 0.95, + max_treedepth = 15, + seed = sample.int(.Machine$integer.max, 1), + strat = NA, + parallel_chains = chains, + with_post = FALSE, + stan_dir = NULL, + compile_dir = NULL, + init = 0.1, + refresh = 200, + show_messages = TRUE, + ...) { + + if (!requireNamespace("cmdstanr", quietly = TRUE)) { + cli::cli_abort(c( + "Package {.pkg cmdstanr} is required.", + "i" = "Install it with: install.packages('cmdstanr', + repos = 'https://mc-stan.org/r-packages/')" + )) + } + if (!requireNamespace("serodynamics", quietly = TRUE)) { + cli::cli_abort("Package {.pkg serodynamics} is required for prep_data().") + } + + model <- match.arg(model) + + # ---- Locate Stan source file ---- + stan_file <- .locate_stan_file(model, stan_dir) + cli::cli_inform(c("i" = "Using Stan file: {.file {stan_file}}")) + + # ---- Determine compile output directory ---- + compile_dir <- .setup_compile_dir(compile_dir) + cli::cli_inform(c("i" = "Compile output directory: {.path {compile_dir}}")) + + # ---- Stratification ---- + if (is.na(strat)) { + strat_list <- "None" + } else { + strat_list <- as.character(unique(data[[strat]])) + } + + combined_out <- list() + stanfit_list <- list() + cov_list <- list() + priors <- NULL + + # ---- Compile model once(cmdstanr caches and avoids repeated filesystem hits) + cli::cli_inform(c("i" = "Compiling {.strong {model}} (or using cache)...")) + mod <- cmdstanr::cmdstan_model( + stan_file = stan_file, + dir = compile_dir, + compile = TRUE + ) + + for (i in strat_list) { + i <- as.character(i) + result <- .run_single_stratum( + stratum = i, data = data, strat = strat, + mod = mod, model = model, chains = chains, + parallel_chains = parallel_chains, + iter_warmup = iter_warmup, iter_sampling = iter_sampling, + seed = seed, adapt_delta = adapt_delta, + max_treedepth = max_treedepth, init = init, + refresh = refresh, show_messages = show_messages, + ... + ) + combined_out[[i]] <- result$sr_tibble + cov_list[[i]] <- result$cov_summaries + stanfit_list[[i]] <- result$stan_fit + priors <- result$priors + } + + if (is.null(priors)) { + cli::cli_abort(c( + "No strata were fitted.", + "i" = paste0("{.code strat_list} appears to be empty;", + " provide at least one stratum.") + )) + } + + sr_out <- dplyr::bind_rows(combined_out) + + sr_out <- sr_out |> + structure( + nChain = chains, + nParameters = 5L, + nIterations = iter_sampling + iter_warmup, + nWarmup = iter_warmup, + model_type = model, + priors = attr(priors, "stan_input_snapshot") + ) + + if (length(cov_list) == 1) { + for (nm in names(cov_list[[1]])) { + attr(sr_out, nm) <- cov_list[[1]][[nm]] + } + } else { + attr(sr_out, "cov_by_stratum") <- cov_list + } + + # Calculate fitted/residuals + # calc_fit_mod is not exported by serodynamics; ::: is a deliberate, + # tolerated compromise until the upstream package exports it. + fit_res <- tryCatch( + serodynamics:::calc_fit_mod(modeled_dat = sr_out, original_data = data), # nolint: namespace_linter + error = function(e) { + cli::cli_warn("calc_fit_mod failed: {e$message}") + NULL + } + ) + if (!is.null(fit_res)) { + attr(sr_out, "fitted_residuals") <- fit_res + } + + if (with_post) { + attr(sr_out, "stan_fit") <- stanfit_list + } + + class(sr_out) <- union("sr_model", class(sr_out)) + return(sr_out) +} diff --git a/R/run_phase0_diagnostic.R b/R/run_phase0_diagnostic.R new file mode 100644 index 00000000..698309ab --- /dev/null +++ b/R/run_phase0_diagnostic.R @@ -0,0 +1,239 @@ +#' Run Phase 0 interactive SLURM reproducibility diagnostic +#' +#' Simulates correlated case data, fits model_2 under an interactive SLURM +#' allocation (`salloc`), extracts diagnostics, and saves a result bundle. +#' Use this to establish a Phase 0 baseline for comparing with Phase 1 +#' sbatch results and confirming determinism across allocation modes. +#' +#' @param n Number of subjects to simulate. +#' @param iter_warmup Number of Stan warmup iterations per chain. +#' @param iter_sampling Number of Stan sampling iterations per chain. +#' @param tag File-naming tag, e.g. `"n5"` or `"n48"`. +#' @param output_dir Directory for output files (default `"outputs/phase0"`). +#' @param true_rho_B True Kronecker biomarker correlation (default `0.6`). +#' @param seed Random seed for simulation and Stan (default `20260513`). +#' @param chains Number of MCMC chains (default `2`). +#' @param adapt_delta Stan `adapt_delta` (default `0.95`). +#' @param max_treedepth Stan `max_treedepth` (default `12`). +#' @param compile_dir Directory for compiled Stan binaries. If `NULL`, +#' defaults to `/tmp//cmdstan_bin_phase0_`. +#' @return Invisibly returns the result bundle list, or `NULL` if the fit +#' crashed. +#' @example inst/examples/run_phase0_diagnostic-examples.R +#' @keywords internal +run_phase0_diagnostic <- function(n, + iter_warmup, + iter_sampling, + tag, + output_dir = "outputs/phase0", + true_rho_B = 0.6, + seed = 20260513L, + chains = 2L, + adapt_delta = 0.95, + max_treedepth = 15L, + compile_dir = NULL) { + cli::cli_h1("PHASE 0: INTERACTIVE SLURM REPRODUCIBILITY TEST ({tag})") + cli::cli_inform(c( + "Purpose: fit via salloc to compare determinism with Phase 1 sbatch", + "Started at: {format(Sys.time())}", + "Host: {Sys.info()[['nodename']]}" + )) + + dir.create("logs/phase0", recursive = TRUE, showWarnings = FALSE) + dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) + + status_file <- file.path(output_dir, "PHASE0_STATUS.txt") + unlink(status_file) + write_status(status_file, "INIT", "Phase 0 started") + + # ----- 1. Package versions ----- + write_status(status_file, "LOAD_PACKAGES", "logging") + pkg_versions <- c( + R = R.version.string, + cmdstanr = .safe_pkg_version("cmdstanr"), + posterior = .safe_pkg_version("posterior"), + shigella = .safe_pkg_version("shigella"), + serodynamics = .safe_pkg_version("serodynamics") + ) + cmdstan_ver <- tryCatch( + cmdstanr::cmdstan_version(), error = function(e) "UNKNOWN" + ) + for (pkg in names(pkg_versions)) { + cat(sprintf(" %-15s %s\n", pkg, pkg_versions[pkg])) # nolint: undesirable_function_linter + } + cat(sprintf(" %-15s %s\n\n", "cmdstan", cmdstan_ver)) # nolint: undesirable_function_linter + saveRDS( + c(pkg_versions, cmdstan = cmdstan_ver), + file.path(output_dir, "env_versions.rds") + ) + write_status(status_file, "LOAD_PACKAGES", "OK") + + # ----- 2. Compile dir ----- + write_status(status_file, "COMPILE_DIR", "setting up") + if (is.null(compile_dir)) { + user <- Sys.getenv("USER", unset = "unknown") + compile_dir <- file.path("/tmp", user, + paste0("cmdstan_bin_phase0_", tag)) + } + if (!dir.exists(compile_dir)) { + dir.create(compile_dir, recursive = TRUE, mode = "0755") + } + cli::cli_inform(c( + "compile_dir: {compile_dir}", + "existing files: {length(list.files(compile_dir))}" + )) + write_status(status_file, "COMPILE_DIR", "OK") + + # ----- 3. Simulate ----- + write_status(status_file, "SIMULATE", "running") + set.seed(seed) + omega_B_true <- matrix(c(1, true_rho_B, true_rho_B, 1), 2, 2) + sim_data <- sim_correlated_case_data( + n = n, + omega_B = omega_B_true, + antigen_isos = c("IgG", "IgA"), + n_obs_per_subject = 5L + ) + cli::cli_inform(sprintf( + "n_subjects: %d, rows: %d, true rho_B: %.3f", + length(unique(sim_data$id)), nrow(sim_data), true_rho_B + )) + saveRDS(sim_data, file.path(output_dir, sprintf("sim_data_%s.rds", tag))) + write_status(status_file, "SIMULATE", "OK") + + # ----- 4. Fit ----- + write_status(status_file, "FIT", "running") + out_file <- file.path(output_dir, sprintf("one_fit_%s.rds", tag)) + t_start <- Sys.time() + saveRDS( + list(scenario = "phase0_interactive", status = "FIT_STARTED", + started_at = format(t_start), true_rho_B = true_rho_B, n = n), + out_file + ) + fit <- tryCatch({ + run_mod_stan( + data = sim_data, + model = "model_2", + chains = chains, + iter_warmup = iter_warmup, + iter_sampling = iter_sampling, + parallel_chains = chains, + adapt_delta = adapt_delta, + max_treedepth = max_treedepth, + init = 0.1, + with_post = TRUE, + compile_dir = compile_dir, + refresh = 100, + show_messages = TRUE + ) + }, error = function(e) { + cli::cli_inform("[FIT ERROR]: {conditionMessage(e)}") + write_status(status_file, "FIT", paste("CRASHED:", conditionMessage(e))) + saveRDS( + list(scenario = "phase0_interactive", status = "FIT_FAILED", + error = conditionMessage(e), crashed_at = format(Sys.time()), + true_rho_B = true_rho_B, n = n), + out_file + ) + NULL + }) + + elapsed <- as.numeric(Sys.time() - t_start, units = "mins") + cli::cli_inform(sprintf("Fit elapsed: %.2f min", elapsed)) + + if (is.null(fit)) { + # write status before signaling — ensures cleanup under options(warn = 2) + write_status(status_file, "DONE", "Phase 0 FAILED") + cli::cli_warn("PHASE 0 RESULT: FIT CRASHED") + return(invisible(NULL)) + } + write_status(status_file, "FIT", "OK") + + # ----- 5. Diagnostics ----- + write_status(status_file, "DIAG", "extracting") + sf <- attr(fit, "stan_fit")[[1]] + diag <- sf$diagnostic_summary( + diagnostics = c("divergences", "treedepth", "ebfmi") + ) + total_iters <- chains * iter_sampling + draws_summary <- tryCatch( + posterior::summarise_draws( + sf$draws(variables = "Omega_B[1,2]"), + "median", "mean", "sd", + ~ quantile(.x, c(0.025, 0.975), na.rm = TRUE), + "ess_bulk", "rhat" + ), + error = function(e) NULL + ) + if (!is.null(draws_summary)) { + cli::cli_inform("Omega_B[1,2] posterior summary:") + print(draws_summary) + } + write_status(status_file, "DIAG", "OK") + + # ----- 6. Save bundle ----- + write_status(status_file, "SAVE", "writing rds") + result_bundle <- list( + scenario = "phase0_interactive", + status = "OK", + elapsed_min = elapsed, + started_at = format(t_start), + completed_at = format(Sys.time()), + host = Sys.info()[["nodename"]], + pkg_versions = pkg_versions, + cmdstan_version = cmdstan_ver, + true_rho_B = true_rho_B, + n_subjects = n, + fit_settings = list( + chains = chains, warmup = iter_warmup, sampling = iter_sampling, + adapt_delta = adapt_delta, max_treedepth = max_treedepth + ), + diagnostic_summary = diag, + omega_B_summary = draws_summary, + rho_B_posterior = tryCatch( + as.vector(posterior::as_draws_array(sf$draws("Omega_B[1,2]"))), + error = function(e) NULL + ) + ) + saveRDS(result_bundle, out_file) + saveRDS( + result_bundle$diagnostic_summary, + file.path(output_dir, sprintf("one_fit_%s_diag.rds", tag)) + ) + write_status(status_file, "SAVE", "OK") + cli::cli_inform("saved -> {out_file}") + + # ----- 7. Summary ----- + cli::cli_h1("PHASE 0 RESULT SUMMARY") + cli::cli_inform(sprintf(" Status: OK")) + cli::cli_inform(sprintf(" Elapsed: %.2f min", elapsed)) + cli::cli_inform(sprintf(" True rho_B: %+.3f", true_rho_B)) + if (!is.null(draws_summary)) { + cli::cli_inform(sprintf(" Recovered median: %+.3f [%.3f, %.3f]", + draws_summary$median, + draws_summary$`2.5%`, + draws_summary$`97.5%`)) + cli::cli_inform( + sprintf(" ESS_bulk: %.0f", draws_summary$ess_bulk) + ) + cli::cli_inform(sprintf(" R-hat: %.3f", draws_summary$rhat)) + } + cli::cli_inform(sprintf(" Divergent: %d / %d", + sum(diag$num_divergent), total_iters)) + cli::cli_inform(sprintf(" Max-treedepth hits: %d / %d", + sum(diag$num_max_treedepth), total_iters)) + + cli::cli_h2("NEXT STEP") + cli::cli_inform(sprintf(" 1. Inspect %s + logs/phase0/*.log", out_file)) + cli::cli_inform(c( + " 2. If divergent rate <= 5% AND R-hat <= 1.01:", + " -> Proceed to Phase 1 (sbatch slurm/phase1_single.sbatch)" + )) + cli::cli_inform(c( + " 3. If divergent rate > 10% OR R-hat > 1.02:", + " -> Skip Phase 1-3, jump to Phase 4 diagnosis." + )) + + write_status(status_file, "DONE", "Phase 0 completed successfully") + invisible(result_bundle) +} diff --git a/R/run_phase1_diagnostic.R b/R/run_phase1_diagnostic.R new file mode 100644 index 00000000..6b2be28a --- /dev/null +++ b/R/run_phase1_diagnostic.R @@ -0,0 +1,301 @@ +#' Run Phase 1 SLURM single-job reproducibility diagnostic +#' +#' Simulates correlated case data (identical seed and parameters to Phase 0), +#' fits model_2 inside a SLURM sbatch job, extracts diagnostics, and saves a +#' result bundle. Optionally compares output with a Phase 0 baseline to +#' isolate SLURM-vs-code attribution of any sampling issues. +#' +#' @param n Number of subjects to simulate (must match Phase 0). +#' @param iter_warmup Number of Stan warmup iterations per chain. +#' @param iter_sampling Number of Stan sampling iterations per chain. +#' @param tag File-naming tag, e.g. `"n5"` or `"n48"`. +#' @param output_dir Directory for output files (default `"outputs/phase1"`). +#' @param phase0_dir Directory containing the Phase 0 result bundle +#' (default `"outputs/phase0"`). Used for the Phase 0 vs Phase 1 +#' comparison table. +#' @param true_rho_B True Kronecker biomarker correlation (default `0.6`). +#' @param seed Random seed - must match Phase 0 (default `20260513`). +#' @param chains Number of MCMC chains (default `2`). +#' @param adapt_delta Stan `adapt_delta` (default `0.95`). +#' @param max_treedepth Stan `max_treedepth` (default `12`). +#' @param compile_dir Directory for compiled Stan binaries. If `NULL`, +#' falls back to the `STAN_COMPILE_DIR` environment variable, then +#' `/tmp//cmdstan_bin_phase1__`. +#' @return Invisibly returns the result bundle list, or `NULL` if the fit +#' crashed. +#' @example inst/examples/run_phase1_diagnostic-examples.R +#' @keywords internal +run_phase1_diagnostic <- function(n, + iter_warmup, + iter_sampling, + tag, + output_dir = "outputs/phase1", + phase0_dir = "outputs/phase0", + true_rho_B = 0.6, + seed = 20260513L, + chains = 2L, + adapt_delta = 0.95, + max_treedepth = 12L, + compile_dir = NULL) { + cli::cli_h1("PHASE 1: SLURM SINGLE JOB DIAGNOSTIC ({tag})") + cli::cli_inform( + "Purpose: fit inside Slurm; compare with Phase 0 to isolate attribution" + ) + + # ----- 0. SLURM environment ----- + slurm_env <- c( + SLURM_JOB_ID = Sys.getenv("SLURM_JOB_ID"), + SLURM_JOB_NAME = Sys.getenv("SLURM_JOB_NAME"), + SLURM_NODELIST = Sys.getenv("SLURM_NODELIST"), + SLURM_CPUS_PER_TASK = Sys.getenv("SLURM_CPUS_PER_TASK"), + SLURM_MEM_PER_NODE = Sys.getenv("SLURM_MEM_PER_NODE"), + SLURM_SUBMIT_DIR = Sys.getenv("SLURM_SUBMIT_DIR"), + USER = Sys.getenv("USER"), + HOSTNAME = Sys.info()[["nodename"]], + TMPDIR = Sys.getenv("TMPDIR") + ) + cli::cli_h2("SLURM environment") + for (env_name in names(slurm_env)) { + cat(sprintf(" %-22s = %s\n", env_name, slurm_env[env_name])) # nolint: undesirable_function_linter + } + cli::cli_inform(c( + "Started at: {format(Sys.time())}", + "R version: {R.version.string}" + )) + + dir.create("logs/phase1", recursive = TRUE, showWarnings = FALSE) + dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) + + job_id <- slurm_env["SLURM_JOB_ID"] + if (job_id == "") job_id <- format(Sys.time(), "%Y%m%d_%H%M%S") + status_file <- file.path(output_dir, + sprintf("PHASE1_STATUS_%s.txt", job_id)) + write_status(status_file, "INIT", + sprintf("Phase 1 started, jobid=%s", job_id)) + + # ----- 1. Package versions ----- + write_status(status_file, "LOAD_PACKAGES", "logging") + pkg_versions <- c( + R = R.version.string, + cmdstanr = .safe_pkg_version("cmdstanr"), + posterior = .safe_pkg_version("posterior"), + shigella = .safe_pkg_version("shigella"), + serodynamics = .safe_pkg_version("serodynamics") + ) + cmdstan_ver <- tryCatch( + cmdstanr::cmdstan_version(), error = function(e) "UNKNOWN" + ) + cli::cli_h2("Package versions") + for (pkg in names(pkg_versions)) { + cat(sprintf(" %-15s %s\n", pkg, pkg_versions[pkg])) # nolint: undesirable_function_linter + } + cat(sprintf(" %-15s %s\n\n", "cmdstan", cmdstan_ver)) # nolint: undesirable_function_linter + write_status(status_file, "LOAD_PACKAGES", "OK") + + # ----- 2. Compile dir (SLURM-specific, per-task subdir) ----- + write_status(status_file, "COMPILE_DIR", "setting up") + if (is.null(compile_dir)) { + compile_dir <- Sys.getenv("STAN_COMPILE_DIR", unset = "") + if (compile_dir == "") { + user <- Sys.getenv("USER", unset = "unknown") + compile_dir <- file.path( + "/tmp", user, sprintf("cmdstan_bin_phase1_%s_%s", tag, job_id) + ) + } + } + if (!dir.exists(compile_dir)) { + dir.create(compile_dir, recursive = TRUE, mode = "0755") + } + cli::cli_inform(c( + "compile_dir: {compile_dir}", + "existing files: {length(list.files(compile_dir))}" + )) + write_status(status_file, "COMPILE_DIR", "OK") + + # ----- 3. Simulate (identical to Phase 0 - same seed, same n) ----- + write_status(status_file, "SIMULATE", "running") + set.seed(seed) + omega_B_true <- matrix(c(1, true_rho_B, true_rho_B, 1), 2, 2) + sim_data <- sim_correlated_case_data( + n = n, + omega_B = omega_B_true, + antigen_isos = c("IgG", "IgA"), + n_obs_per_subject = 5L + ) + cli::cli_inform(sprintf( + "Simulated: n=%d, rho_B=%.1f, total rows=%d", + length(unique(sim_data$id)), true_rho_B, nrow(sim_data) + )) + write_status(status_file, "SIMULATE", "OK") + + # ----- 4. Save STARTED placeholder ----- + scenario <- sprintf("phase1_slurm_single_%s", tag) + out_file <- file.path(output_dir, + sprintf("one_fit_%s_jobid_%s.rds", tag, job_id)) + saveRDS( + list(scenario = scenario, status = "FIT_STARTED", job_id = job_id, + started_at = format(Sys.time()), slurm_env = slurm_env, + pkg_versions = pkg_versions, cmdstan_version = cmdstan_ver, + true_rho_B = true_rho_B, n = n), + out_file + ) + + # ----- 5. Fit ----- + write_status(status_file, "FIT", "running") + cli::cli_h2("Fitting model_2 (Kronecker)") + t_start <- Sys.time() + fit <- tryCatch({ + run_mod_stan( + data = sim_data, + model = "model_2", + chains = chains, + iter_warmup = iter_warmup, + iter_sampling = iter_sampling, + parallel_chains = chains, + adapt_delta = adapt_delta, + max_treedepth = max_treedepth, + init = 0.1, + with_post = TRUE, + compile_dir = compile_dir, + refresh = 100, + show_messages = TRUE + ) + }, error = function(e) { + cli::cli_inform("[FIT ERROR]: {conditionMessage(e)}") + write_status(status_file, "FIT", paste("CRASHED:", conditionMessage(e))) + saveRDS( + list(scenario = scenario, status = "FIT_FAILED", job_id = job_id, + error = conditionMessage(e), crashed_at = format(Sys.time()), + slurm_env = slurm_env, pkg_versions = pkg_versions), + out_file + ) + NULL + }) + + elapsed <- as.numeric(Sys.time() - t_start, units = "mins") + cli::cli_inform(sprintf("Fit elapsed: %.2f min", elapsed)) + + if (is.null(fit)) { + phase0_file <- file.path(phase0_dir, sprintf("one_fit_%s.rds", tag)) + # write status before signaling — ensures cleanup under options(warn = 2) + write_status(status_file, "DONE", "Phase 1 FAILED") + cli::cli_warn("PHASE 1 RESULT: FIT CRASHED INSIDE SLURM") + cli::cli_inform("Compare with {phase0_file} to determine:") + cli::cli_inform(c( + " " = "If Phase 0 OK but Phase 1 FAIL -> Slurm env issue", + " " = "If both fail -> code / model identifiability issue" + )) + return(invisible(NULL)) + } + write_status(status_file, "FIT", "OK") + + # ----- 6. Diagnostics ----- + write_status(status_file, "DIAG", "extracting") + sf <- attr(fit, "stan_fit")[[1]] + diag <- sf$diagnostic_summary( + diagnostics = c("divergences", "treedepth", "ebfmi") + ) + total_iters <- chains * iter_sampling + draws_summary <- tryCatch( + posterior::summarise_draws( + sf$draws(variables = "Omega_B[1,2]"), + "median", "mean", "sd", + ~ quantile(.x, c(0.025, 0.975), na.rm = TRUE), + "ess_bulk", "rhat" + ), + error = function(e) NULL + ) + cli::cli_h2("Diagnostics") + cli::cli_inform(sprintf(" divergent: %d / %d (%.2f%%)", + sum(diag$num_divergent), total_iters, + 100 * sum(diag$num_divergent) / total_iters)) + cli::cli_inform(sprintf(" max-treedepth: %d / %d (%.2f%%)", + sum(diag$num_max_treedepth), total_iters, + 100 * sum(diag$num_max_treedepth) / total_iters)) + cli::cli_inform(sprintf(" E-BFMI: %s", + paste(sprintf("%.3f", diag$ebfmi), collapse = ", "))) + if (!is.null(draws_summary)) { + cli::cli_inform("Omega_B[1,2] summary:") + print(draws_summary) + } else { + cli::cli_inform("[INFO] Omega_B[1,2] not available in this fit.") + } + write_status(status_file, "DIAG", "OK") + + # ----- 7. Save bundle ----- + write_status(status_file, "SAVE", "writing rds") + result_bundle <- list( + scenario = scenario, + status = "OK", + job_id = job_id, + elapsed_min = elapsed, + started_at = format(t_start), + completed_at = format(Sys.time()), + slurm_env = slurm_env, + pkg_versions = pkg_versions, + cmdstan_version = cmdstan_ver, + true_rho_B = true_rho_B, + n_subjects = n, + fit_settings = list( + chains = chains, warmup = iter_warmup, sampling = iter_sampling, + adapt_delta = adapt_delta, max_treedepth = max_treedepth + ), + diagnostic_summary = diag, + omega_B_summary = draws_summary, + rho_B_posterior = tryCatch( + as.vector(posterior::as_draws_array(sf$draws("Omega_B[1,2]"))), + error = function(e) NULL + ) + ) + saveRDS(result_bundle, out_file) + + # ----- 8. Compare with Phase 0 if available ----- + phase0_file <- file.path(phase0_dir, sprintf("one_fit_%s.rds", tag)) + if (file.exists(phase0_file)) { + ph0 <- readRDS(phase0_file) + if (!is.null(ph0$omega_B_summary) && !is.null(draws_summary)) { + cli::cli_h2("Phase 0 vs Phase 1 comparison") + cmp <- data.frame( + metric = c("status", "elapsed_min", "post_median", + "post_lo_2.5", "post_hi_97.5", + "ess_bulk", "rhat", "n_divergent", "n_treedepth"), + phase0 = c( + ph0$status, + round(ph0$elapsed_min, 2), + round(ph0$omega_B_summary$median, 3), + round(ph0$omega_B_summary$`2.5%`, 3), + round(ph0$omega_B_summary$`97.5%`, 3), + round(ph0$omega_B_summary$ess_bulk, 0), + round(ph0$omega_B_summary$rhat, 3), + sum(ph0$diagnostic_summary$num_divergent), + sum(ph0$diagnostic_summary$num_max_treedepth) + ), + phase1 = c( + "OK", + round(elapsed, 2), + round(draws_summary$median, 3), + round(draws_summary$`2.5%`, 3), + round(draws_summary$`97.5%`, 3), + round(draws_summary$ess_bulk, 0), + round(draws_summary$rhat, 3), + sum(diag$num_divergent), + sum(diag$num_max_treedepth) + ) + ) + print(cmp, row.names = FALSE) + saveRDS(cmp, file.path(output_dir, + sprintf("p0_vs_p1_comparison_%s.rds", job_id))) + } + } else { + cli::cli_inform(c( + "[INFO] Phase 0 result not found - comparison skipped.", + paste0(" Run phase0_interactive_reproducibility_{tag}.R", + " first for direct comparison.") + )) + } + + write_status(status_file, "SAVE", "OK") + cli::cli_inform("Phase 1 complete. Results: {out_file}") + write_status(status_file, "DONE", "Phase 1 OK") + invisible(result_bundle) +} diff --git a/R/run_single_stratum.R b/R/run_single_stratum.R new file mode 100644 index 00000000..24a0bf39 --- /dev/null +++ b/R/run_single_stratum.R @@ -0,0 +1,60 @@ +# Helper: run prep + sample + postprocess for one stratum. +# Accepts the full dataset plus strat/stratum identifiers and slices +# internally, so the caller loop body only needs one function call. +# Returns a list with sr_tibble, cov_summaries, stan_fit, and priors. +#' @keywords internal +#' @noRd +.run_single_stratum <- function(stratum, data, strat, + mod, model, chains, parallel_chains, + iter_warmup, iter_sampling, seed, + adapt_delta, max_treedepth, init, + refresh, show_messages, ...) { + dl_sub <- if (is.na(strat)) { + data + } else { + sub <- data[data[[strat]] == stratum, , drop = FALSE] + # Restore all non-structural attributes (class, id_var, biomarker_var, + # time_in_days, value_var, etc.) dropped by [.data.frame subsetting. + standard_attrs <- c("names", "class", "row.names") + for (a in setdiff(names(attributes(data)), standard_attrs)) { + attr(sub, a) <- attr(data, a) + } + class(sub) <- class(data) + sub + } + + prepped <- serodynamics::prep_data(dl_sub, add_newperson = FALSE) + stan_data <- prep_data_stan(prepped) + priors <- prep_priors_stan(model = model, ...) + full_data <- c(stan_data, priors) + + cli::cli_inform(c("i" = "Sampling {.strong {model}} with {chains} chains...")) + fit <- mod$sample( + data = full_data, + chains = chains, + parallel_chains = parallel_chains, + iter_warmup = iter_warmup, + iter_sampling = iter_sampling, + seed = seed, + adapt_delta = adapt_delta, + max_treedepth = max_treedepth, + init = init, + refresh = refresh, + show_messages = show_messages + ) + + processed <- postprocess_stan_output( + stan_fit = fit, + ids = attr(stan_data, "ids"), + antigens = attr(stan_data, "antigens"), + model = model, + stratification = stratum + ) + + list( + sr_tibble = processed$sr_tibble, + cov_summaries = processed$cov_summaries, + stan_fit = fit, + priors = priors + ) +} diff --git a/R/safe_pkg_version.R b/R/safe_pkg_version.R new file mode 100644 index 00000000..a4ab9066 --- /dev/null +++ b/R/safe_pkg_version.R @@ -0,0 +1,11 @@ +# Returns the installed version of a package as a string, or "not installed" +# if the package is absent. Used for robust version-logging in diagnostic +# functions where suggested packages may not be present. +#' @keywords internal +#' @noRd +.safe_pkg_version <- function(pkg) { + tryCatch( + as.character(utils::packageVersion(pkg)), + error = function(e) "not installed" + ) +} diff --git a/R/setup_compile_dir.R b/R/setup_compile_dir.R new file mode 100644 index 00000000..0d862931 --- /dev/null +++ b/R/setup_compile_dir.R @@ -0,0 +1,17 @@ +# Helper: resolve and create the compile output directory. +# Priority: argument > environment variable > /tmp fallback. +#' @keywords internal +#' @noRd +.setup_compile_dir <- function(compile_dir) { + if (is.null(compile_dir)) { + compile_dir <- Sys.getenv("STAN_COMPILE_DIR", unset = "") + if (compile_dir == "") { + user <- Sys.getenv("USER", unset = "default") + compile_dir <- file.path("/tmp", user, "cmdstan_bin") + } + } + if (!dir.exists(compile_dir)) { + dir.create(compile_dir, recursive = TRUE, mode = "0755") + } + compile_dir +} diff --git a/R/shigella-package.R b/R/shigella-package.R index a65cf643..13fc44bd 100644 --- a/R/shigella-package.R +++ b/R/shigella-package.R @@ -1,4 +1,5 @@ #' @keywords internal +#' @importFrom stats median rnorm "_PACKAGE" ## usethis namespace: start diff --git a/R/sim_correlated_case_data.R b/R/sim_correlated_case_data.R new file mode 100644 index 00000000..15231e0d --- /dev/null +++ b/R/sim_correlated_case_data.R @@ -0,0 +1,162 @@ +#' @title Simulate correlated longitudinal case data +#' @description +#' Extends [serodynamics::sim_case_data()] to inject known correlation +#' structure at two levels: +#' +#' 1. **Parameter-level (between-biomarker) correlation** via a Kronecker +#' covariance on the vectorised per-subject parameter matrix. +#' 2. **Residual-level correlation** via multivariate log-normal +#' observation noise. +#' +#' **Data-generating process** +#' +#' *Subject parameters.* Let \eqn{\Theta_i} be the \eqn{P \times K} +#' matrix of log-scale kinetic parameters for subject \eqn{i} +#' (rows = parameters, columns = biomarkers). Draw +#' \deqn{ +#' \mathrm{vec}(\Theta_i) \sim +#' \mathcal{N}\!\bigl(\mathrm{vec}(M),\; +#' \Sigma_B \otimes \Sigma_P\bigr), +#' } +#' where \eqn{M} is a \eqn{P \times K} population-mean matrix, +#' \eqn{\Sigma_P = \mathrm{diag}(\tau_P)\,\Omega_P\,\mathrm{diag}(\tau_P)} +#' is the \eqn{P \times P} within-biomarker parameter covariance, and +#' \eqn{\Sigma_B = \mathrm{diag}(\tau_B)\,\Omega_B\,\mathrm{diag}(\tau_B)} +#' is the \eqn{K \times K} between-biomarker covariance. +#' +#' *Observation model.* For each subject \eqn{i}, time \eqn{t}, and +#' biomarker \eqn{k}, +#' \deqn{ +#' \log y_{i,t,k} = \log \mu_{i,t,k} + \varepsilon_{i,t,k}, +#' \quad +#' \boldsymbol{\varepsilon}_{i,t} \sim +#' \mathcal{N}(\mathbf{0},\, \Sigma_\varepsilon), +#' } +#' where +#' \eqn{\Sigma_\varepsilon = +#' \mathrm{diag}(\tau_\varepsilon)\,\Omega_\varepsilon\, +#' \mathrm{diag}(\tau_\varepsilon)}. +#' +#' *Two-phase kinetics.* The deterministic log-mean follows +#' \deqn{ +#' \log \mu_{i,t,k} = +#' \begin{cases} +#' \log(y0_{ik}) + \beta_{ik}\,t, & t \le t1_{ik},\\[4pt] +#' \dfrac{1}{1-s_{ik}} +#' \log\!\bigl(y1_{ik}^{1-s_{ik}} +#' - (1-s_{ik})\,\alpha_{ik}(t - t1_{ik})\bigr), +#' & t > t1_{ik}, +#' \end{cases} +#' } +#' with growth rate +#' \eqn{\beta_{ik} = (\log y1_{ik} - \log y0_{ik})\,/\,t1_{ik}} +#' and shape \eqn{s_{ik} = \exp(\mathtt{log\_rm1}_{ik}) + 1 > 1}. +#' +#' This is the data-generating process for a Kronecker-correlated simulation +#' study. +#' +#' @param n [integer] number of individuals to simulate +#' @param mu [numeric] length-P vector of population means on log scale +#' (defaults match JAGS prep_priors) +#' @param tau_P [numeric] length-P vector of SDs across kinetic +#' parameters +#' @param tau_B [numeric] length-K vector of SDs across biomarkers +#' @param tau_eps [numeric] length-K vector of residual SDs +#' @param omega_P [matrix] P x P parameter correlation matrix +#' (default: identity - no within-biomarker parameter correlation) +#' @param omega_B [matrix] K x K biomarker correlation matrix +#' (default: identity - Scenario 2, residual correlation only) +#' @param omega_eps [matrix] K x K residual correlation matrix +#' (default: identity - no residual correlation) +#' @param antigen_isos [character] names for the K biomarkers +#' @param n_obs_per_subject [integer] number of observations per subject +#' (default 5, matching the Shigella SOSAR cohort) +#' @param time_grid [numeric] follow-up times in days +#' (default c(2, 7, 30, 90, 180)) +#' @param seed [integer] RNG seed +#' +#' @returns a `case_data` object plus attributes recording the truth: +#' - `"truth"` - list with mu, tau_P, tau_B, tau_eps, omega_P, +#' omega_B, omega_eps, sigma_P, sigma_B, sigma_eps +#' - `"theta_true"` - N x P x K array of true subject parameters in +#' **Stan's internal log-scale parameterisation**: +#' `log_y0`, `log_y1m0`, `log_t1`, `log_alpha`, `log_rm1`. +#' Fitted posterior summaries use natural-scale names +#' (`y0`, `y1`, `t1`, `alpha`, `shape`); apply `exp()` before comparing. +#' @export +#' @example inst/examples/sim_correlated_case_data-examples.R +sim_correlated_case_data <- function( + n = 48, + mu = c(1.0, 7.0, 1.0, -4.0, -1.0), + tau_P = c(0.5, 0.7, 0.3, 1.0, 0.4), + tau_B = c(0.8, 0.8), + tau_eps = c(0.3, 0.3), + omega_P = diag(5), + omega_B = diag(2), + omega_eps = diag(2), + antigen_isos = c("biomarker_1", "biomarker_2"), + n_obs_per_subject = 5L, + time_grid = c(2, 7, 30, 90, 180), + seed = NULL) { + + if (!is.null(seed)) { + set.seed(seed) + } + + n_param <- length(mu) + n_biomarker <- length(antigen_isos) + + .validate_sim_inputs( + n_param, n_biomarker, + tau_B, tau_eps, tau_P, + omega_P, omega_B, omega_eps, + time_grid, n_obs_per_subject + ) + + .validate_corr_matrix(omega_P, "omega_P") + .validate_corr_matrix(omega_B, "omega_B") + .validate_corr_matrix(omega_eps, "omega_eps") + + mats <- .build_sigma_matrices(mu, tau_P, tau_B, tau_eps, + omega_P, omega_B, omega_eps) + sigma_eps <- mats$sigma_eps + sigma_full <- mats$sigma_full + mu_vec <- mats$mu_vec + + theta_arr <- .draw_subject_params(n, mu_vec, sigma_full, + n_param, n_biomarker, antigen_isos) + + # --- Generate observations --- + u_eps <- chol(sigma_eps) + + rows <- .generate_obs_rows(n, n_obs_per_subject, time_grid, + n_biomarker, theta_arr, antigen_isos, u_eps) + + sim_df <- dplyr::bind_rows(rows) + + # Convert to case_data + case <- sim_df |> + serodynamics::as_case_data( + id_var = "id", + biomarker_var = "antigen_iso", + time_in_days = "timeindays", + value_var = "value" + ) + + # Attach ground truth + attr(case, "truth") <- list( + mu = mu, + tau_P = tau_P, + tau_B = tau_B, + tau_eps = tau_eps, + omega_P = omega_P, + omega_B = omega_B, + omega_eps = omega_eps, + sigma_P = diag(tau_P) %*% omega_P %*% diag(tau_P), + sigma_B = diag(tau_B) %*% omega_B %*% diag(tau_B), + sigma_eps = sigma_eps + ) + attr(case, "theta_true") <- theta_arr + + return(case) +} diff --git a/R/summarize_matrix_array.R b/R/summarize_matrix_array.R new file mode 100644 index 00000000..56d5bd7f --- /dev/null +++ b/R/summarize_matrix_array.R @@ -0,0 +1,32 @@ +#' Summarize an array-of-matrices variable from posterior draws +#' +#' For Stan variables declared as `array[n_arr] matrix[n_row, n_col]`, +#' cmdstanr names cells `var[k,i,j]`. Returns a list of `n_arr` matrices. +#' +#' @details Returns element-wise posterior medians only. Credible intervals +#' are not computed. For full posterior summaries use +#' [posterior::summarise_draws()] directly. +#' +#' @keywords internal +#' @noRd +.summarize_matrix_array <- function(draws_arr, var_name, + n_arr, n_row, n_col) { + var_dim <- dimnames(draws_arr)$variable + + lapply(seq_len(n_arr), function(k) { + mat <- matrix(NA_real_, nrow = n_row, ncol = n_col) + + for (i in seq_len(n_row)) { + for (j in seq_len(n_col)) { + cell_name <- sprintf("%s[%d,%d,%d]", var_name, k, i, j) + + if (cell_name %in% var_dim) { + cell_draws <- as.numeric(draws_arr[, , cell_name]) + mat[i, j] <- median(cell_draws, na.rm = TRUE) + } + } + } + + mat + }) +} diff --git a/R/summarize_matrix_draws.R b/R/summarize_matrix_draws.R new file mode 100644 index 00000000..e19fdaec --- /dev/null +++ b/R/summarize_matrix_draws.R @@ -0,0 +1,22 @@ +#' Summarize a matrix variable from posterior draws +#' +#' @details Returns element-wise posterior medians only. Credible intervals +#' are not computed. For full posterior summaries use +#' [posterior::summarise_draws()] directly. +#' +#' @keywords internal +#' @noRd +.summarize_matrix_draws <- function(draws_arr, var_name, nrow, ncol) { + result <- matrix(NA_real_, nrow = nrow, ncol = ncol) + var_dim <- dimnames(draws_arr)$variable + for (i in seq_len(nrow)) { + for (j in seq_len(ncol)) { + cell_name <- sprintf("%s[%d,%d]", var_name, i, j) + if (cell_name %in% var_dim) { + cell_draws <- as.numeric(draws_arr[, , cell_name]) + result[i, j] <- median(cell_draws, na.rm = TRUE) + } + } + } + result +} diff --git a/R/validate_corr_matrix.R b/R/validate_corr_matrix.R new file mode 100644 index 00000000..ca9795e6 --- /dev/null +++ b/R/validate_corr_matrix.R @@ -0,0 +1,29 @@ +#' Validate that a matrix is a correlation matrix +#' +#' @param M Numeric matrix to validate. +#' @param name Character; name to use in error messages. +#' @param tol Numeric tolerance for symmetry and unit-diagonal checks. +#' @keywords internal +#' @noRd +.validate_corr_matrix <- function(M, name, tol = 1e-8) { + if (!is.matrix(M) || !is.numeric(M)) { + cli::cli_abort("{.arg {name}} must be a numeric matrix.") + } + if (nrow(M) != ncol(M)) { + cli::cli_abort("{.arg {name}} must be square; got {nrow(M)} x {ncol(M)}.") + } + if (max(abs(M - t(M))) > tol) { + cli::cli_abort("{.arg {name}} must be symmetric.") + } + if (max(abs(diag(M) - 1)) > tol) { + cli::cli_abort("{.arg {name}} must have unit diagonal.") + } + eig <- eigen(M, symmetric = TRUE, only.values = TRUE)$values + if (min(eig) < -tol) { + cli::cli_abort(c( + "{.arg {name}} must be positive semi-definite.", + "i" = "Smallest eigenvalue: {min(eig)}" + )) + } + invisible(M) +} diff --git a/R/validate_sim_inputs.R b/R/validate_sim_inputs.R new file mode 100644 index 00000000..8187c01f --- /dev/null +++ b/R/validate_sim_inputs.R @@ -0,0 +1,39 @@ +# Helper: validate simulation inputs. +#' @keywords internal +#' @noRd +.validate_sim_inputs <- function(n_param, n_biomarker, + tau_B, tau_eps, tau_P, + omega_P, omega_B, omega_eps, + time_grid, n_obs_per_subject) { + if (n_biomarker != length(tau_B)) { + cli::cli_abort("{.arg tau_B} must have length K.") + } + + if (n_biomarker != length(tau_eps)) { + cli::cli_abort("{.arg tau_eps} must have length K.") + } + + if (n_param != length(tau_P)) { + cli::cli_abort("{.arg tau_P} must have length P.") + } + + if (!is.matrix(omega_P) || !identical(dim(omega_P), c(n_param, n_param))) { + cli::cli_abort("{.arg omega_P} must be a P x P matrix.") + } + + if (!is.matrix(omega_B) || !identical(dim(omega_B), c(n_biomarker, + n_biomarker))) { + cli::cli_abort("{.arg omega_B} must be a K x K matrix.") + } + + if (!is.matrix(omega_eps) || !identical(dim(omega_eps), c(n_biomarker, + n_biomarker))) { + cli::cli_abort("{.arg omega_eps} must be a K x K matrix.") + } + + if (length(time_grid) < n_obs_per_subject) { + cli::cli_abort( + "{.arg time_grid} must have at least {.arg n_obs_per_subject} entries." + ) + } +} diff --git a/R/validate_stan_arrays.R b/R/validate_stan_arrays.R new file mode 100644 index 00000000..8adb557b --- /dev/null +++ b/R/validate_stan_arrays.R @@ -0,0 +1,15 @@ +# Helper: sanity-check array sizes and zero-observation subjects. +#' @keywords internal +#' @noRd +.validate_stan_arrays <- function(nsmpl, max_obs) { + if (any(nsmpl > max_obs)) { + cli::cli_abort( + "n_obs[{which(nsmpl > max_obs)}] > max_obs. Array sizes inconsistent." + ) + } + if (any(nsmpl == 0)) { + cli::cli_warn( + "Subject(s) with 0 observations detected; these contribute no likelihood." + ) + } +} diff --git a/R/write_status.R b/R/write_status.R new file mode 100644 index 00000000..2407509b --- /dev/null +++ b/R/write_status.R @@ -0,0 +1,19 @@ +#' Write a step-status line to a diagnostic log file +#' +#' Appends a timestamped status entry to a log file. Intended for use in +#' long-running HPC diagnostic scripts so crash location is visible even +#' when the script dies mid-way. +#' +#' @param status_file Path to the status log file (character). +#' @param step Character label for the current step. +#' @param msg Optional message string (default `""`). +#' @examples +#' log_file <- tempfile(fileext = ".txt") +#' write_status(log_file, "INIT", "starting") +#' write_status(log_file, "DONE") +#' readLines(log_file) +#' @export +write_status <- function(status_file, step, msg = "") { + cat(sprintf("[%s] STEP=%s | %s\n", format(Sys.time()), step, msg), + file = status_file, append = TRUE) +} diff --git a/do_commit.sh b/do_commit.sh new file mode 100644 index 00000000..e69de29b diff --git a/inst/WORDLIST b/inst/WORDLIST index d92948f0..d3d393e2 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,23 +1,95 @@ +Betancourt +Bonačić +Bordetella +Bürkner CMD +CRC +CmdStan +CmdStanMCMC Codecov -GitHub +Cov +CrI +Crowther +Diekmann +Dunson +Eijkeren +Gelman +Girolami +Graaf +Guikema +HMC +HPC IgA IgG -Iso -Isotype -MFI +IpaB +Kretzschmar +Kurowicka +LKJ +Lewandowski +Marinović +Mattoo ORCID -Postprocess +SDs +SLURM SOSAR -Seroresponse -df -ggplot -ipab +Schrader +Stavnezer +Talts +Teunis +Treedepth +UC +Vehtari +al +attr +bigl +bigr +biomarker +biomarkers +boldsymbol +cdot +cholesky +cmdstan +cmdstanr +cmdstanr's +conda +cov +de +diag +env +eps +estimand +et +hyp isotype -isotypes +iter +ldots +le +lpdf +mathbb +mathbf +mathcal +mathrm +mucosal +multimodality +neq newperson -pre -repo +observationally +otimes +params +pathogenesis +qquad +rstan +sbatch +sd serodynamics -serotype +seroincidence +serological +seroresponse +serosurveillance +sr +stan +stanfit tibble +tmp +treedepth +widehat diff --git a/inst/examples/postprocess_stan_output-examples.R b/inst/examples/postprocess_stan_output-examples.R new file mode 100644 index 00000000..0526a69e --- /dev/null +++ b/inst/examples/postprocess_stan_output-examples.R @@ -0,0 +1,46 @@ +## Example: postprocess_stan_output() +## +## Convert a raw cmdstanr fit object into a list with sr_tibble +## (tidy parameter draws) and cov_summaries (covariance matrices). +## Requires a compiled cmdstan installation. + +if (interactive()) { +if (requireNamespace("cmdstanr", quietly = TRUE)) { + + set.seed(2026) + + sim_data <- sim_correlated_case_data( + n = 5, + omega_B = diag(2), + antigen_isos = c("IgG", "IgA"), + n_obs_per_subject = 5L + ) + + stan_data <- prep_data_stan(sim_data) + priors <- prep_priors_stan(model = "model_2") + + ## Compile and sample directly (without run_mod_stan wrapper) + mod <- cmdstanr::cmdstan_model( + system.file("stan", "model_2.stan", package = "shigella") + ) + raw_fit <- mod$sample( + data = c(stan_data, priors), + chains = 1, + iter_warmup = 200, + iter_sampling = 200, + refresh = 0, + show_messages = FALSE + ) + + ## Post-process the raw fit into tidy sr_model format + processed <- postprocess_stan_output( + stan_fit = raw_fit, + ids = attr(stan_data, "ids"), + antigens = attr(stan_data, "antigens"), + model = "model_2" + ) + + class(processed$sr_tibble) + names(processed$cov_summaries) +} +} diff --git a/inst/examples/prep_data_stan-examples.R b/inst/examples/prep_data_stan-examples.R new file mode 100644 index 00000000..0a09613d --- /dev/null +++ b/inst/examples/prep_data_stan-examples.R @@ -0,0 +1,11 @@ +## Example: prep_data_stan() +## +## Convert a case_data object directly into the list format that +## Stan models expect. Internally calls serodynamics::prep_data() +## with add_newperson = FALSE. + +sim <- sim_correlated_case_data(n = 5, seed = 2026) + +stan_data <- prep_data_stan(sim) + +str(stan_data) diff --git a/inst/examples/prep_priors_stan-examples.R b/inst/examples/prep_priors_stan-examples.R new file mode 100644 index 00000000..fc2fc6a0 --- /dev/null +++ b/inst/examples/prep_priors_stan-examples.R @@ -0,0 +1,7 @@ +## Example: prep_priors_stan() +## +## Return the prior hyperparameters for the Kronecker correlated Stan model. + +priors <- prep_priors_stan(model = "model_2") + +str(priors) diff --git a/inst/examples/run_mod_stan-examples.R b/inst/examples/run_mod_stan-examples.R new file mode 100644 index 00000000..036b3f3b --- /dev/null +++ b/inst/examples/run_mod_stan-examples.R @@ -0,0 +1,33 @@ +## Example: run_mod_stan() +## +## Fit the Kronecker Stan model on a small synthetic dataset. +## Uses minimal MCMC settings so the example completes quickly. +## Requires a compiled cmdstan installation. + +if (interactive()) { + +if (requireNamespace("cmdstanr", quietly = TRUE)) { + + set.seed(2026) + + sim_data <- sim_correlated_case_data( + n = 5, + omega_B = matrix(c(1, 0.6, 0.6, 1), nrow = 2), + antigen_isos = c("IgG", "IgA"), + n_obs_per_subject = 5L + ) + + fit <- run_mod_stan( + data = sim_data, + model = "model_2", + chains = 1, + iter_warmup = 200, + iter_sampling = 200, + refresh = 0, + show_messages = FALSE + ) + + class(fit) + head(fit) +} +} diff --git a/inst/examples/run_phase0_diagnostic-examples.R b/inst/examples/run_phase0_diagnostic-examples.R new file mode 100644 index 00000000..e7d86207 --- /dev/null +++ b/inst/examples/run_phase0_diagnostic-examples.R @@ -0,0 +1,22 @@ +## Example: run_phase0_diagnostic() +## +## Runs the Phase 0 interactive SLURM reproducibility diagnostic. +## Requires a compiled cmdstan installation and is intended for use +## under an interactive SLURM allocation (salloc), not CI. + +if (interactive()) { + + if (requireNamespace("cmdstanr", quietly = TRUE)) { + + result <- run_phase0_diagnostic( + n = 5, + iter_warmup = 200, + iter_sampling = 200, + tag = "n5", + output_dir = file.path(tempdir(), "phase0"), + chains = 1L + ) + + names(result) # status, elapsed_min, diagnostic_summary, omega_B_summary, ... + } +} diff --git a/inst/examples/run_phase1_diagnostic-examples.R b/inst/examples/run_phase1_diagnostic-examples.R new file mode 100644 index 00000000..8259e14b --- /dev/null +++ b/inst/examples/run_phase1_diagnostic-examples.R @@ -0,0 +1,23 @@ +## Example: run_phase1_diagnostic() +## +## Runs the Phase 1 SLURM single-job reproducibility diagnostic. +## Intended for use inside a SLURM sbatch job. Requires a compiled +## cmdstan installation and matching Phase 0 outputs for comparison. + +if (interactive()) { + + if (requireNamespace("cmdstanr", quietly = TRUE)) { + + result <- run_phase1_diagnostic( + n = 5, + iter_warmup = 200, + iter_sampling = 200, + tag = "n5", + output_dir = file.path(tempdir(), "phase1"), + phase0_dir = file.path(tempdir(), "phase0"), + chains = 1L + ) + + names(result) # status, elapsed_min, diagnostic_summary, omega_B_summary, ... + } +} diff --git a/inst/examples/sim_correlated_case_data-examples.R b/inst/examples/sim_correlated_case_data-examples.R new file mode 100644 index 00000000..66c60d1b --- /dev/null +++ b/inst/examples/sim_correlated_case_data-examples.R @@ -0,0 +1,18 @@ +## Example +## +## Generate synthetic Shigella antibody-kinetics data with a known +## IgG-IgA correlation rho_B. + +set.seed(2026) + +omega_B <- matrix(c(1.0, 0.6, + 0.6, 1.0), nrow = 2) + +sim_data <- sim_correlated_case_data( + n = 5, + omega_B = omega_B, + antigen_isos = c("IgG", "IgA"), + n_obs_per_subject = 5L +) + +head(sim_data) diff --git a/inst/examples/write_status-examples.R b/inst/examples/write_status-examples.R new file mode 100644 index 00000000..f32c2cd0 --- /dev/null +++ b/inst/examples/write_status-examples.R @@ -0,0 +1,17 @@ +## Example: write_status() +## +## Appends a timestamped step-status line to a log file. + +if (interactive()) { + + log_file <- tempfile(fileext = ".txt") + + write_status(log_file, "INIT", "starting workflow") + write_status(log_file, "STEP_1", "data loaded") + write_status(log_file, "STEP_2") # msg defaults to "" + write_status(log_file, "DONE", "workflow complete") + + readLines(log_file) + unlink(log_file) + +} diff --git a/inst/stan/model_1.stan b/inst/stan/model_1.stan new file mode 100644 index 00000000..796793a3 --- /dev/null +++ b/inst/stan/model_1.stan @@ -0,0 +1,155 @@ +// ============================================================ +// model_1.stan +// +// Model 1 (formerly Model A): Independent biomarker antibody kinetics. +// Each biomarker fit independently — no cross-biomarker correlation. +// +// Log-space kinetics: uses log_two_phase_curve() ported from model_2.stan. +// This removes the discontinuous fallback (no `if base <= 0`) and aligns +// the likelihood implementation with the validated model_2.stan. +// +// Compatible with serodynamics::prep_data_stan() output: +// N, K, P, max_obs, n_obs[N], time_obs[N, max_obs], log_y[N, max_obs, K] +// +// Compatible with prep_priors_stan(model = "model_1") output. +// ============================================================ + +functions { + // Compute log(y(t)) DIRECTLY (matching JAGS reference and model_2.stan) + // Uses log-space parameters where natural: + // - log_y0, log_y1 are log of baseline / peak antibody + // - t1, alpha, shape are in natural scale + real log_two_phase_curve(real t, + real log_y0, real log_y1, + real t1, real alpha, real shape) { + if (t <= t1) { + // Active phase: log(y(t)) = log(y0) + beta * t + // where beta = (log(y1) - log(y0)) / t1 + real beta = (log_y1 - log_y0) / t1; + return log_y0 + beta * t; + } else { + // Recovery phase: log(y(t)) = 1/(1-shape) * log(inside) + // inside = y1^(1-shape) - (1-shape)*alpha*(t-t1) + // + // Since shape > 1, (1 - shape) < 0: + // y1^(1-shape) = exp((1-shape) * log_y1) (small positive) + // -(1-shape)*alpha*(t-t1) = (shape-1)*alpha*(t-t1) (positive) + // So inside = small_positive + positive = positive ✓ + // No need for fallback because both terms are guaranteed positive. + real one_minus_shape = 1 - shape; + real first_term = exp(one_minus_shape * log_y1); + real second_term = (shape - 1) * alpha * (t - t1); + real inside = first_term + second_term; + return log(inside) / one_minus_shape; + } + } +} + +data { + int N; + int K; + int P; + int max_obs; + array[N] int n_obs; + array[N, max_obs] real time_obs; + array[N, max_obs, K] real log_y; + + vector[P] mu_hyp_mean; + vector[P] mu_hyp_sd; + real tau_P_scale; + real tau_eps_scale; + real lkj_P_eta; + // lkj_eps_eta is NOT used here: model_1 uses independent per-biomarker + // residual scales (tau_eps[k]), not an LKJ prior on epsilon correlation. + // For the LKJ epsilon prior, see model_2.stan. +} + +parameters { + matrix[K, P] M; + array[K] cholesky_factor_corr[P] L_Omega_P; + array[K] vector[P] tau_P; + vector[K] tau_eps; + array[K] matrix[N, P] Z; +} + +transformed parameters { + array[N, K] vector[P] theta; + for (k in 1:K) { + matrix[P, P] L_Sigma_P_k = diag_pre_multiply(tau_P[k], L_Omega_P[k]); + for (i in 1:N) { + theta[i, k] = M[k]' + L_Sigma_P_k * Z[k][i]'; + } + } +} + +model { + for (k in 1:K) { + for (p in 1:P) { + M[k, p] ~ normal(mu_hyp_mean[p], mu_hyp_sd[p]); + } + L_Omega_P[k] ~ lkj_corr_cholesky(lkj_P_eta); + tau_P[k] ~ cauchy(0, tau_P_scale); + to_vector(Z[k]) ~ std_normal(); + } + tau_eps ~ cauchy(0, tau_eps_scale); + + for (i in 1:N) { + if (n_obs[i] > 0) { + for (o in 1:n_obs[i]) { + for (k in 1:K) { + real log_y0_o = theta[i, k][1]; + real log_y1_o = log_sum_exp(log_y0_o, theta[i, k][2]); + real t1_o = exp(theta[i, k][3]); + real alpha_o = exp(theta[i, k][4]); + real shape_o = exp(theta[i, k][5]) + 1; + real mu_log = log_two_phase_curve(time_obs[i, o], log_y0_o, log_y1_o, + t1_o, alpha_o, shape_o); + log_y[i, o, k] ~ normal(mu_log, tau_eps[k]); + } + } + } + } +} + +generated quantities { + // NOTE: Kinetics recomputed here intentionally — Stan scoping requires + // local variables to be redefined; this is not duplication that can be eliminated. + array[K] corr_matrix[P] Omega_P; + for (k in 1:K) { + Omega_P[k] = multiply_lower_tri_self_transpose(L_Omega_P[k]); + } + + array[N, K] real y0; + array[N, K] real y1; + array[N, K] real t1; + array[N, K] real alpha; + array[N, K] real shape; + for (i in 1:N) { + for (k in 1:K) { + y0[i, k] = exp(theta[i, k][1]); + y1[i, k] = y0[i, k] + exp(theta[i, k][2]); + t1[i, k] = exp(theta[i, k][3]); + alpha[i, k] = exp(theta[i, k][4]); + shape[i, k] = exp(theta[i, k][5]) + 1; + } + } + + vector[N] log_lik; + for (i in 1:N) { + log_lik[i] = 0; + if (n_obs[i] > 0) { + for (o in 1:n_obs[i]) { + for (k in 1:K) { + real log_y0_o = theta[i, k][1]; + real log_y1_o = log_sum_exp(log_y0_o, theta[i, k][2]); + real t1_o = exp(theta[i, k][3]); + real alpha_o = exp(theta[i, k][4]); + real shape_o = exp(theta[i, k][5]) + 1; + real mu_log = log_two_phase_curve(time_obs[i, o], log_y0_o, log_y1_o, + t1_o, alpha_o, shape_o); + log_lik[i] += normal_lpdf(log_y[i, o, k] | mu_log, tau_eps[k]); + } + } + } + } +} diff --git a/inst/stan/model_2.stan b/inst/stan/model_2.stan new file mode 100644 index 00000000..355d06d6 --- /dev/null +++ b/inst/stan/model_2.stan @@ -0,0 +1,228 @@ +// ============================================================ +// model_2.stan — JAGS-ALIGNED version (v3) +// +// Key change from previous version: compute log(y(t)) DIRECTLY +// (matching the JAGS reference model.jags from Chapter 1), +// rather than computing y(t) then taking log(). +// +// This avoids the exp -> arithmetic -> log round-trip, which: +// 1. Avoids overflow when y1 is large +// 2. Removes the need for a discontinuous fallback (no `if base <= 0`) +// 3. Aligns numerically with the JAGS model that works for Chapter 1 +// +// Also: y1 = y0 + exp(Theta[k,2]), so log(y1) = log_sum_exp(log_y0, Theta[k,2]) +// This is more stable than log(exp(.) + exp(.)). +// +// Model 2: Kronecker correlated antibody kinetics. +// vec(Theta_i) ~ MVN_KP(vec(M), Sigma_B kron Sigma_P) +// log y[i,o,1:K] ~ MVN_K(mu_log[i,o,1:K], Sigma_eps) +// ============================================================ + +functions { + // Compute log(y(t)) DIRECTLY (matching JAGS reference) + // Uses log-space parameters where natural: + // - log_y0, log_y1 are log of baseline / peak antibody + // - t1, alpha, shape are in natural scale + real log_two_phase_curve(real t, + real log_y0, real log_y1, + real t1, real alpha, real shape) { + if (t <= t1) { + // Active phase: log(y(t)) = log(y0) + beta * t + // where beta = (log(y1) - log(y0)) / t1 + real beta = (log_y1 - log_y0) / t1; + return log_y0 + beta * t; + } else { + // Recovery phase: log(y(t)) = 1/(1-shape) * log(inside) + // inside = y1^(1-shape) - (1-shape)*alpha*(t-t1) + // + // Since shape > 1, (1 - shape) < 0: + // y1^(1-shape) = exp((1-shape) * log_y1) (small positive) + // -(1-shape)*alpha*(t-t1) = (shape-1)*alpha*(t-t1) (positive) + // So inside = small_positive + positive = positive ✓ + // No need for fallback because both terms are guaranteed positive. + real one_minus_shape = 1 - shape; + real first_term = exp(one_minus_shape * log_y1); + real second_term = (shape - 1) * alpha * (t - t1); + real inside = first_term + second_term; + return log(inside) / one_minus_shape; + } + } + + matrix kron_chol(matrix L_B, matrix L_P) { + int K = rows(L_B); + int P = rows(L_P); + matrix[K * P, K * P] L_out = rep_matrix(0, K * P, K * P); + for (i in 1:K) { + for (j in 1:i) { + for (p in 1:P) { + for (q in 1:p) { + L_out[(i - 1) * P + p, (j - 1) * P + q] = L_B[i, j] * L_P[p, q]; + } + } + } + } + return L_out; + } +} + +data { + int N; + int K; + int P; + int max_obs; + array[N] int n_obs; + array[N, max_obs] real time_obs; + array[N, max_obs, K] real log_y; + + vector[P] mu_hyp_mean; + vector[P] mu_hyp_sd; + real tau_P_scale; + real tau_B_scale; + real tau_eps_scale; + real lkj_P_eta; + real lkj_B_eta; + real lkj_eps_eta; +} + +parameters { + matrix[K, P] M; + + cholesky_factor_corr[K] L_Omega_B; + cholesky_factor_corr[P] L_Omega_P; + vector[K] tau_B; + vector[P] tau_P; + + cholesky_factor_corr[K] L_Omega_eps; + vector[K] tau_eps; + + matrix[N, K * P] Z; +} + +transformed parameters { + array[N] matrix[K, P] Theta; + + matrix[K, K] L_Sigma_B = diag_pre_multiply(tau_B, L_Omega_B); + matrix[P, P] L_Sigma_P = diag_pre_multiply(tau_P, L_Omega_P); + matrix[K * P, K * P] L_kron = kron_chol(L_Sigma_B, L_Sigma_P); + + vector[K * P] mu_vec; + for (k in 1:K) { + for (p in 1:P) { + mu_vec[(k - 1) * P + p] = M[k, p]; + } + } + + for (i in 1:N) { + vector[K * P] theta_vec = mu_vec + L_kron * to_vector(Z[i]); + for (k in 1:K) { + for (p in 1:P) { + Theta[i, k, p] = theta_vec[(k - 1) * P + p]; + } + } + } +} + +model { + // Priors on M (matches JAGS mu.par ~ dmnorm structure) + for (k in 1:K) { + for (p in 1:P) { + M[k, p] ~ normal(mu_hyp_mean[p], mu_hyp_sd[p]); + } + } + + L_Omega_B ~ lkj_corr_cholesky(lkj_B_eta); + L_Omega_P ~ lkj_corr_cholesky(lkj_P_eta); + L_Omega_eps ~ lkj_corr_cholesky(lkj_eps_eta); + + tau_B ~ cauchy(0, tau_B_scale); + tau_P ~ cauchy(0, tau_P_scale); + tau_eps ~ cauchy(0, tau_eps_scale); + + to_vector(Z) ~ std_normal(); + + matrix[K, K] L_Sigma_eps = diag_pre_multiply(tau_eps, L_Omega_eps); + + // Likelihood — compute log(y(t)) directly (matching JAGS) + for (i in 1:N) { + if (n_obs[i] > 0) { + for (o in 1:n_obs[i]) { + vector[K] mu_log_o; + vector[K] y_log_o; + for (k in 1:K) { + // Theta[i, k, *] holds: + // [1] = log(y0) + // [2] = log(y1 - y0) -> y1 = y0 + exp(par[2]) + // [3] = log(t1) + // [4] = log(alpha) + // [5] = log(shape - 1) -> shape = exp(par[5]) + 1 + real log_y0_o = Theta[i, k, 1]; + // y1 = y0 + exp(par2) = exp(log_y0) + exp(par2) + // log(y1) = log_sum_exp(log_y0, par2) -- numerically stable + real log_y1_o = log_sum_exp(log_y0_o, Theta[i, k, 2]); + real t1_o = exp(Theta[i, k, 3]); + real alpha_o = exp(Theta[i, k, 4]); + real shape_o = exp(Theta[i, k, 5]) + 1; + + mu_log_o[k] = log_two_phase_curve(time_obs[i, o], + log_y0_o, log_y1_o, + t1_o, alpha_o, shape_o); + y_log_o[k] = log_y[i, o, k]; + } + y_log_o ~ multi_normal_cholesky(mu_log_o, L_Sigma_eps); + } + } + } +} + +generated quantities { + // NOTE: Kinetics recomputed here intentionally — Stan scoping requires + // local variables to be redefined; this is not duplication that can be eliminated. + corr_matrix[K] Omega_B = multiply_lower_tri_self_transpose(L_Omega_B); + corr_matrix[P] Omega_P = multiply_lower_tri_self_transpose(L_Omega_P); + corr_matrix[K] Omega_eps = multiply_lower_tri_self_transpose(L_Omega_eps); + cov_matrix[K] Sigma_B = quad_form_diag(Omega_B, tau_B); + cov_matrix[P] Sigma_P = quad_form_diag(Omega_P, tau_P); + cov_matrix[K] Sigma_eps = quad_form_diag(Omega_eps, tau_eps); + + array[N, K] real y0; + array[N, K] real y1; + array[N, K] real t1; + array[N, K] real alpha; + array[N, K] real shape; + for (i in 1:N) { + for (k in 1:K) { + y0[i, k] = exp(Theta[i, k, 1]); + y1[i, k] = y0[i, k] + exp(Theta[i, k, 2]); + t1[i, k] = exp(Theta[i, k, 3]); + alpha[i, k] = exp(Theta[i, k, 4]); + shape[i, k] = exp(Theta[i, k, 5]) + 1; + } + } + + vector[N] log_lik; + { + // L_Sigma_eps rebuilt here per Stan scoping; intentional, not accidental. + matrix[K, K] L_Sigma_eps = diag_pre_multiply(tau_eps, L_Omega_eps); + for (i in 1:N) { + log_lik[i] = 0; + if (n_obs[i] > 0) { + for (o in 1:n_obs[i]) { + vector[K] mu_log_o; + vector[K] y_log_o; + for (k in 1:K) { + real log_y0_o = Theta[i, k, 1]; + real log_y1_o = log_sum_exp(log_y0_o, Theta[i, k, 2]); + real t1_o = exp(Theta[i, k, 3]); + real alpha_o = exp(Theta[i, k, 4]); + real shape_o = exp(Theta[i, k, 5]) + 1; + mu_log_o[k] = log_two_phase_curve(time_obs[i, o], + log_y0_o, log_y1_o, + t1_o, alpha_o, shape_o); + y_log_o[k] = log_y[i, o, k]; + } + log_lik[i] += multi_normal_cholesky_lpdf(y_log_o | mu_log_o, L_Sigma_eps); + } + } + } + } +} diff --git a/man/postprocess_stan_output.Rd b/man/postprocess_stan_output.Rd new file mode 100644 index 00000000..a85efe72 --- /dev/null +++ b/man/postprocess_stan_output.Rd @@ -0,0 +1,88 @@ +% 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{Post-process Stan output to sr_model format (cmdstanr version)} +\usage{ +postprocess_stan_output( + stan_fit, + ids, + antigens, + model = c("model_2", "model_1"), + stratification = "None", + param_names = c("y0", "y1", "t1", "alpha", "shape") +) +} +\arguments{ +\item{stan_fit}{CmdStanMCMC object from cmdstanr (Not rstan stanfit)} + +\item{ids}{subject IDs from attr(stan_data, "ids")} + +\item{antigens}{biomarker names from attr(stan_data, "antigens")} + +\item{model}{"model_1", "model_2"} + +\item{stratification}{label for this stratum} + +\item{param_names}{Character vector of parameter names to extract from the +Stan model's \code{generated quantities} block. Must match the variable +names exactly as declared in both \file{inst/stan/model_1.stan} and +\file{inst/stan/model_2.stan}. Defaults to +\code{c("y0", "y1", "t1", "alpha", "shape")}.} +} +\value{ +list with sr_tibble and cov_summaries +} +\description{ +Converts a CmdStanMCMC object (from cmdstanr's mod$sample()) into the +long-format tibble produced by run_mod(), so downstream plotting/summary +functions work without modification. +} +\examples{ +## Example: postprocess_stan_output() +## +## Convert a raw cmdstanr fit object into a list with sr_tibble +## (tidy parameter draws) and cov_summaries (covariance matrices). +## Requires a compiled cmdstan installation. + +if (interactive()) { +if (requireNamespace("cmdstanr", quietly = TRUE)) { + + set.seed(2026) + + sim_data <- sim_correlated_case_data( + n = 5, + omega_B = diag(2), + antigen_isos = c("IgG", "IgA"), + n_obs_per_subject = 5L + ) + + stan_data <- prep_data_stan(sim_data) + priors <- prep_priors_stan(model = "model_2") + + ## Compile and sample directly (without run_mod_stan wrapper) + mod <- cmdstanr::cmdstan_model( + system.file("stan", "model_2.stan", package = "shigella") + ) + raw_fit <- mod$sample( + data = c(stan_data, priors), + chains = 1, + iter_warmup = 200, + iter_sampling = 200, + refresh = 0, + show_messages = FALSE + ) + + ## Post-process the raw fit into tidy sr_model format + processed <- postprocess_stan_output( + stan_fit = raw_fit, + ids = attr(stan_data, "ids"), + antigens = attr(stan_data, "antigens"), + model = "model_2" + ) + + class(processed$sr_tibble) + names(processed$cov_summaries) +} +} +} diff --git a/man/prep_data_stan.Rd b/man/prep_data_stan.Rd new file mode 100644 index 00000000..c09f35fa --- /dev/null +++ b/man/prep_data_stan.Rd @@ -0,0 +1,62 @@ +% 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 the Stan backend} +\usage{ +prep_data_stan(data, drop_newperson = TRUE) +} +\arguments{ +\item{data}{either a \code{case_data} object (output of +\code{\link[=sim_correlated_case_data]{sim_correlated_case_data()}} or \code{\link[serodynamics:as_case_data]{serodynamics::as_case_data()}}) +or a \code{prepped_jags_data} object (output of +\code{\link[serodynamics:prep_data]{serodynamics::prep_data()}}).} + +\item{drop_newperson}{\link{logical} whether to drop the JAGS dummy +"newperson" row if it is present in a \code{prepped_jags_data} +input. Default \code{TRUE}. Has no effect when \code{data} is a +\code{case_data} object because the internal \code{prep_data()} call uses +\code{add_newperson = FALSE}.} +} +\value{ +a named \link{list} with attributes \code{ids} and \code{antigens}, +ready to pass to a compiled Stan model via +\code{cmdstanr::cmdstan_model()$sample()}. +} +\description{ +Converts case data into the structured list that the Stan models +(\code{model_1.stan}, \code{model_2.stan}) expect. Accepts either a raw +\code{case_data} object (the typical entry point) or a +\code{prepped_jags_data} object already produced by +\code{\link[serodynamics:prep_data]{serodynamics::prep_data()}} (the JAGS-side prep step). + +When given a \code{case_data} object, this function internally calls +\code{\link[serodynamics:prep_data]{serodynamics::prep_data()}} with \code{add_newperson = FALSE} (Stan +handles posterior prediction in the \verb{generated quantities} block +rather than via a dummy missing-data subject). + +The Stan models expect: +\itemize{ +\item \code{N}: number of subjects +\item \code{K}: number of antigen-isotype biomarkers +\item \code{P}: number of kinetic parameters (always 5) +\item \code{max_obs}: max number of observations per subject +\item \code{n_obs[N]}: actual number of observations per subject +\item \code{time_obs[N, max_obs]}: observation times (NA -> 0, ignored by +the likelihood via the \code{n_obs[i]} guard) +\item \code{log_y[N, max_obs, K]}: log-transformed antibody observations +} +} +\examples{ +## Example: prep_data_stan() +## +## Convert a case_data object directly into the list format that +## Stan models expect. Internally calls serodynamics::prep_data() +## with add_newperson = FALSE. + +sim <- sim_correlated_case_data(n = 5, seed = 2026) + +stan_data <- prep_data_stan(sim) + +str(stan_data) +} diff --git a/man/prep_priors_stan.Rd b/man/prep_priors_stan.Rd new file mode 100644 index 00000000..308ed2e3 --- /dev/null +++ b/man/prep_priors_stan.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prep_priors_stan.R +\name{prep_priors_stan} +\alias{prep_priors_stan} +\title{Prepare priors for Stan backend} +\usage{ +prep_priors_stan( + mu_hyp_mean = c(1, 7, 1, -4, -1), + mu_hyp_sd = c(5, 5, 5, 5, 5), + tau_P_scale = 1, + tau_B_scale = 1, + tau_eps_scale = 1, + lkj_P_eta = 2, + lkj_B_eta = 1, + lkj_eps_eta = 2, + model = c("model_2", "model_1") +) +} +\arguments{ +\item{mu_hyp_mean}{\link{numeric} length-5 prior mean for population params} + +\item{mu_hyp_sd}{\link{numeric} length-5 prior SD for population params. +Weakly informative, Stan-friendly. JAGS original was c(1, 316, 1, 32, 1). +5.0 on log-scale params covers ~5 orders of magnitude - plenty wide.} + +\item{tau_P_scale}{half-Cauchy scale for parameter SDs} + +\item{tau_B_scale}{half-Cauchy scale for biomarker SDs (Model 2 only)} + +\item{tau_eps_scale}{half-Cauchy scale for residual SDs} + +\item{lkj_P_eta}{LKJ shape for parameter correlation} + +\item{lkj_B_eta}{LKJ shape for biomarker correlation (Model 2 only)} + +\item{lkj_eps_eta}{LKJ shape for residual correlation} + +\item{model}{character: "model_1", "model_2"} +} +\value{ +named list with priors for the Stan data block +} +\description{ +Translates the JAGS prior specification into Stan's LKJ + half-Cauchy +decomposition. + +Defaults match the JAGS prior specification, with two adjustments for +Stan compatibility: +\itemize{ +\item mu_hyp_sd capped at ~5 (was 316 in JAGS) - Stan's HMC sampler +handles weakly-informative priors better with more reasonable scales. +JAGS Gibbs sampling tolerates wider priors, but Stan's gradient-based +HMC explores the tails too aggressively when sd is huge. +\item tau scales = 1.0 (was 2.5) - keeps initial steps reasonable +} +} +\examples{ +## Example: prep_priors_stan() +## +## Return the prior hyperparameters for the Kronecker correlated Stan model. + +priors <- prep_priors_stan(model = "model_2") + +str(priors) +} diff --git a/man/run_mod_stan.Rd b/man/run_mod_stan.Rd new file mode 100644 index 00000000..03270614 --- /dev/null +++ b/man/run_mod_stan.Rd @@ -0,0 +1,125 @@ +% 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 using the cmdstanr backend} +\usage{ +run_mod_stan( + data, + model = c("model_2", "model_1"), + chains = 4, + iter_sampling = 1000, + iter_warmup = 1000, + adapt_delta = 0.95, + max_treedepth = 12, + seed = sample.int(.Machine$integer.max, 1), + strat = NA, + parallel_chains = chains, + with_post = FALSE, + stan_dir = NULL, + compile_dir = NULL, + init = 0.1, + refresh = 200, + show_messages = TRUE, + ... +) +} +\arguments{ +\item{data}{case_data object (from sim_correlated_case_data() or +as_case_data())} + +\item{model}{character: "model_1", "model_2"} + +\item{chains}{Number of chains to run.} + +\item{iter_sampling}{Number of post-warmup iterations per chain.} + +\item{iter_warmup}{Number of warmup iterations per chain.} + +\item{adapt_delta}{Target average acceptance probability for Stan sampling.} + +\item{max_treedepth}{Maximum tree depth for Stan NUTS sampling.} + +\item{seed}{Random seed passed to Stan.} + +\item{strat}{optional stratification variable (default NA)} + +\item{parallel_chains}{Number of chains to run in parallel.} + +\item{with_post}{return raw CmdStanMCMC object as attribute (default FALSE)} + +\item{stan_dir}{Optional directory containing \code{model_*.stan} files. +If \code{NULL}, the function first looks for Stan files installed with the +package using \code{system.file("stan", ..., package = "shigella")}, then falls +back to \code{inst/stan} for interactive development.} + +\item{compile_dir}{directory where cmdstanr writes compiled binaries. +Default uses STAN_COMPILE_DIR env var, or +/tmp/\if{html}{\out{}}/cmdstan_bin.} + +\item{init}{initial value strategy. Numeric value scales down random init +(default 0.1 to avoid -inf in multi_normal_cholesky_lpdf)} + +\item{refresh}{Stan progress refresh interval.} + +\item{show_messages}{Logical; whether to show CmdStan messages.} + +\item{...}{additional priors passed to prep_priors_stan()} +} +\value{ +sr_model tibble +} +\description{ +Fits the two-phase antibody kinetics model using \strong{cmdstanr}. +The \code{compile_dir} argument allows compiled Stan binaries to be written to +a writable directory (default: \verb{/tmp}), which is useful on HPC systems +where the home directory is mounted non-executable. + +Output: an \code{sr_model} tibble with the same column schema as \code{run_mod()}, +so all existing plot / summary functions work unchanged. Stan-specific +attributes are also attached: +\itemize{ +\item \code{Omega_eps}, \code{Sigma_eps}: residual covariance (Model 2 only) +\item \code{Omega_B}, \code{Sigma_B}: biomarker covariance (Model 2 only) +\item \code{Omega_P}: parameter corr matrix (Model 2: single matrix; +Model 1: named list of K matrices, one per +biomarker) +\item \code{Sigma_P}: parameter covariance (Model 2 only) +\item \code{stan_fit}: raw CmdStanMCMC object (when with_post = TRUE) +} +} +\examples{ +## Example: run_mod_stan() +## +## Fit the Kronecker Stan model on a small synthetic dataset. +## Uses minimal MCMC settings so the example completes quickly. +## Requires a compiled cmdstan installation. + +if (interactive()) { + +if (requireNamespace("cmdstanr", quietly = TRUE)) { + + set.seed(2026) + + sim_data <- sim_correlated_case_data( + n = 5, + omega_B = matrix(c(1, 0.6, 0.6, 1), nrow = 2), + antigen_isos = c("IgG", "IgA"), + n_obs_per_subject = 5L + ) + + fit <- run_mod_stan( + data = sim_data, + model = "model_2", + chains = 1, + iter_warmup = 200, + iter_sampling = 200, + refresh = 0, + show_messages = FALSE + ) + + class(fit) + head(fit) +} +} +} diff --git a/man/run_phase0_diagnostic.Rd b/man/run_phase0_diagnostic.Rd new file mode 100644 index 00000000..fbd92cd9 --- /dev/null +++ b/man/run_phase0_diagnostic.Rd @@ -0,0 +1,79 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_phase0_diagnostic.R +\name{run_phase0_diagnostic} +\alias{run_phase0_diagnostic} +\title{Run Phase 0 interactive SLURM reproducibility diagnostic} +\usage{ +run_phase0_diagnostic( + n, + iter_warmup, + iter_sampling, + tag, + output_dir = "outputs/phase0", + true_rho_B = 0.6, + seed = 20260513L, + chains = 2L, + adapt_delta = 0.95, + max_treedepth = 12L, + compile_dir = NULL +) +} +\arguments{ +\item{n}{Number of subjects to simulate.} + +\item{iter_warmup}{Number of Stan warmup iterations per chain.} + +\item{iter_sampling}{Number of Stan sampling iterations per chain.} + +\item{tag}{File-naming tag, e.g. \code{"n5"} or \code{"n48"}.} + +\item{output_dir}{Directory for output files (default \code{"outputs/phase0"}).} + +\item{true_rho_B}{True Kronecker biomarker correlation (default \code{0.6}).} + +\item{seed}{Random seed for simulation and Stan (default \code{20260513}).} + +\item{chains}{Number of MCMC chains (default \code{2}).} + +\item{adapt_delta}{Stan \code{adapt_delta} (default \code{0.95}).} + +\item{max_treedepth}{Stan \code{max_treedepth} (default \code{12}).} + +\item{compile_dir}{Directory for compiled Stan binaries. If \code{NULL}, +defaults to \verb{/tmp//cmdstan_bin_phase0_}.} +} +\value{ +Invisibly returns the result bundle list, or \code{NULL} if the fit +crashed. +} +\description{ +Simulates correlated case data, fits model_2 under an interactive SLURM +allocation (\code{salloc}), extracts diagnostics, and saves a result bundle. +Use this to establish a Phase 0 baseline for comparing with Phase 1 +sbatch results and confirming determinism across allocation modes. +} +\examples{ +## Example: run_phase0_diagnostic() +## +## Runs the Phase 0 interactive SLURM reproducibility diagnostic. +## Requires a compiled cmdstan installation and is intended for use +## under an interactive SLURM allocation (salloc), not CI. + +if (interactive()) { + + if (requireNamespace("cmdstanr", quietly = TRUE)) { + + result <- run_phase0_diagnostic( + n = 5, + iter_warmup = 200, + iter_sampling = 200, + tag = "n5", + output_dir = file.path(tempdir(), "phase0"), + chains = 1L + ) + + names(result) # status, elapsed_min, diagnostic_summary, omega_B_summary, ... + } +} +} +\keyword{internal} diff --git a/man/run_phase1_diagnostic.Rd b/man/run_phase1_diagnostic.Rd new file mode 100644 index 00000000..9ef6eb9b --- /dev/null +++ b/man/run_phase1_diagnostic.Rd @@ -0,0 +1,86 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_phase1_diagnostic.R +\name{run_phase1_diagnostic} +\alias{run_phase1_diagnostic} +\title{Run Phase 1 SLURM single-job reproducibility diagnostic} +\usage{ +run_phase1_diagnostic( + n, + iter_warmup, + iter_sampling, + tag, + output_dir = "outputs/phase1", + phase0_dir = "outputs/phase0", + true_rho_B = 0.6, + seed = 20260513L, + chains = 2L, + adapt_delta = 0.95, + max_treedepth = 12L, + compile_dir = NULL +) +} +\arguments{ +\item{n}{Number of subjects to simulate (must match Phase 0).} + +\item{iter_warmup}{Number of Stan warmup iterations per chain.} + +\item{iter_sampling}{Number of Stan sampling iterations per chain.} + +\item{tag}{File-naming tag, e.g. \code{"n5"} or \code{"n48"}.} + +\item{output_dir}{Directory for output files (default \code{"outputs/phase1"}).} + +\item{phase0_dir}{Directory containing the Phase 0 result bundle +(default \code{"outputs/phase0"}). Used for the Phase 0 vs Phase 1 +comparison table.} + +\item{true_rho_B}{True Kronecker biomarker correlation (default \code{0.6}).} + +\item{seed}{Random seed - must match Phase 0 (default \code{20260513}).} + +\item{chains}{Number of MCMC chains (default \code{2}).} + +\item{adapt_delta}{Stan \code{adapt_delta} (default \code{0.95}).} + +\item{max_treedepth}{Stan \code{max_treedepth} (default \code{12}).} + +\item{compile_dir}{Directory for compiled Stan binaries. If \code{NULL}, +falls back to the \code{STAN_COMPILE_DIR} environment variable, then +\verb{/tmp//cmdstan_bin_phase1__}.} +} +\value{ +Invisibly returns the result bundle list, or \code{NULL} if the fit +crashed. +} +\description{ +Simulates correlated case data (identical seed and parameters to Phase 0), +fits model_2 inside a SLURM sbatch job, extracts diagnostics, and saves a +result bundle. Optionally compares output with a Phase 0 baseline to +isolate SLURM-vs-code attribution of any sampling issues. +} +\examples{ +## Example: run_phase1_diagnostic() +## +## Runs the Phase 1 SLURM single-job reproducibility diagnostic. +## Intended for use inside a SLURM sbatch job. Requires a compiled +## cmdstan installation and matching Phase 0 outputs for comparison. + +if (interactive()) { + + if (requireNamespace("cmdstanr", quietly = TRUE)) { + + result <- run_phase1_diagnostic( + n = 5, + iter_warmup = 200, + iter_sampling = 200, + tag = "n5", + output_dir = file.path(tempdir(), "phase1"), + phase0_dir = file.path(tempdir(), "phase0"), + chains = 1L + ) + + names(result) # status, elapsed_min, diagnostic_summary, omega_B_summary, ... + } +} +} +\keyword{internal} diff --git a/man/shigella-package.Rd b/man/shigella-package.Rd index 1f0e8287..ec9af0c6 100644 --- a/man/shigella-package.Rd +++ b/man/shigella-package.Rd @@ -4,9 +4,9 @@ \name{shigella-package} \alias{shigella} \alias{shigella-package} -\title{shigella: What the Package Does (One Line, Title Case)} +\title{shigella: Bayesian Modeling of Shigella Antibody Kinetics} \description{ -What the package does (one paragraph). +Tools for multivariate Bayesian hierarchical modeling of antibody response trajectories following confirmed Shigella infection, supporting kinetic parameter estimation and serosurveillance applications via Stan (cmdstanr) backends. } \seealso{ Useful links: diff --git a/man/sim_correlated_case_data.Rd b/man/sim_correlated_case_data.Rd new file mode 100644 index 00000000..26e734fb --- /dev/null +++ b/man/sim_correlated_case_data.Rd @@ -0,0 +1,142 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sim_correlated_case_data.R +\name{sim_correlated_case_data} +\alias{sim_correlated_case_data} +\title{Simulate correlated longitudinal case data} +\usage{ +sim_correlated_case_data( + n = 48, + mu = c(1, 7, 1, -4, -1), + tau_P = c(0.5, 0.7, 0.3, 1, 0.4), + tau_B = c(0.8, 0.8), + tau_eps = c(0.3, 0.3), + omega_P = diag(5), + omega_B = diag(2), + omega_eps = diag(2), + antigen_isos = c("biomarker_1", "biomarker_2"), + n_obs_per_subject = 5L, + time_grid = c(2, 7, 30, 90, 180), + seed = NULL +) +} +\arguments{ +\item{n}{\link{integer} number of individuals to simulate} + +\item{mu}{\link{numeric} length-P vector of population means on log scale +(defaults match JAGS prep_priors)} + +\item{tau_P}{\link{numeric} length-P vector of SDs across kinetic +parameters} + +\item{tau_B}{\link{numeric} length-K vector of SDs across biomarkers} + +\item{tau_eps}{\link{numeric} length-K vector of residual SDs} + +\item{omega_P}{\link{matrix} P x P parameter correlation matrix +(default: identity - no within-biomarker parameter correlation)} + +\item{omega_B}{\link{matrix} K x K biomarker correlation matrix +(default: identity - Scenario 2, residual correlation only)} + +\item{omega_eps}{\link{matrix} K x K residual correlation matrix +(default: identity - no residual correlation)} + +\item{antigen_isos}{\link{character} names for the K biomarkers} + +\item{n_obs_per_subject}{\link{integer} number of observations per subject +(default 5, matching the Shigella SOSAR cohort)} + +\item{time_grid}{\link{numeric} follow-up times in days +(default c(2, 7, 30, 90, 180))} + +\item{seed}{\link{integer} RNG seed} +} +\value{ +a \code{case_data} object plus attributes recording the truth: +\itemize{ +\item \code{"truth"} - list with mu, tau_P, tau_B, tau_eps, omega_P, +omega_B, omega_eps, sigma_P, sigma_B, sigma_eps +\item \code{"theta_true"} - N x P x K array of true subject parameters in +\strong{Stan's internal log-scale parameterisation}: +\code{log_y0}, \code{log_y1m0}, \code{log_t1}, \code{log_alpha}, \code{log_rm1}. +Fitted posterior summaries use natural-scale names +(\code{y0}, \code{y1}, \code{t1}, \code{alpha}, \code{shape}); apply \code{exp()} before comparing. +} +} +\description{ +Extends \code{\link[serodynamics:sim_case_data]{serodynamics::sim_case_data()}} to inject known correlation +structure at two levels: +\enumerate{ +\item \strong{Parameter-level (between-biomarker) correlation} via a Kronecker +covariance on the vectorised per-subject parameter matrix. +\item \strong{Residual-level correlation} via multivariate log-normal +observation noise. +} + +\strong{Data-generating process} + +\emph{Subject parameters.} Let \eqn{\Theta_i} be the \eqn{P \times K} +matrix of log-scale kinetic parameters for subject \eqn{i} +(rows = parameters, columns = biomarkers). Draw +\deqn{ + \mathrm{vec}(\Theta_i) \sim + \mathcal{N}\!\bigl(\mathrm{vec}(M),\; + \Sigma_B \otimes \Sigma_P\bigr), +} +where \eqn{M} is a \eqn{P \times K} population-mean matrix, +\eqn{\Sigma_P = \mathrm{diag}(\tau_P)\,\Omega_P\,\mathrm{diag}(\tau_P)} +is the \eqn{P \times P} within-biomarker parameter covariance, and +\eqn{\Sigma_B = \mathrm{diag}(\tau_B)\,\Omega_B\,\mathrm{diag}(\tau_B)} +is the \eqn{K \times K} between-biomarker covariance. + +\emph{Observation model.} For each subject \eqn{i}, time \eqn{t}, and +biomarker \eqn{k}, +\deqn{ + \log y_{i,t,k} = \log \mu_{i,t,k} + \varepsilon_{i,t,k}, + \quad + \boldsymbol{\varepsilon}_{i,t} \sim + \mathcal{N}(\mathbf{0},\, \Sigma_\varepsilon), +} +where +\eqn{\Sigma_\varepsilon = + \mathrm{diag}(\tau_\varepsilon)\,\Omega_\varepsilon\, + \mathrm{diag}(\tau_\varepsilon)}. + +\emph{Two-phase kinetics.} The deterministic log-mean follows +\deqn{ + \log \mu_{i,t,k} = + \begin{cases} + \log(y0_{ik}) + \beta_{ik}\,t, & t \le t1_{ik},\\[4pt] + \dfrac{1}{1-s_{ik}} + \log\!\bigl(y1_{ik}^{1-s_{ik}} + - (1-s_{ik})\,\alpha_{ik}(t - t1_{ik})\bigr), + & t > t1_{ik}, + \end{cases} +} +with growth rate +\eqn{\beta_{ik} = (\log y1_{ik} - \log y0_{ik})\,/\,t1_{ik}} +and shape \eqn{s_{ik} = \exp(\mathtt{log\_rm1}_{ik}) + 1 > 1}. + +This is the data-generating process for a Kronecker-correlated simulation +study. +} +\examples{ +## Example +## +## Generate synthetic Shigella antibody-kinetics data with a known +## IgG-IgA correlation rho_B. + +set.seed(2026) + +omega_B <- matrix(c(1.0, 0.6, + 0.6, 1.0), nrow = 2) + +sim_data <- sim_correlated_case_data( + n = 5, + omega_B = omega_B, + antigen_isos = c("IgG", "IgA"), + n_obs_per_subject = 5L +) + +head(sim_data) +} diff --git a/man/write_status.Rd b/man/write_status.Rd new file mode 100644 index 00000000..f40f40af --- /dev/null +++ b/man/write_status.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/write_status.R +\name{write_status} +\alias{write_status} +\title{Write a step-status line to a diagnostic log file} +\usage{ +write_status(status_file, step, msg = "") +} +\arguments{ +\item{status_file}{Path to the status log file (character).} + +\item{step}{Character label for the current step.} + +\item{msg}{Optional message string (default \code{""}).} +} +\description{ +Appends a timestamped status entry to a log file. Intended for use in +long-running HPC diagnostic scripts so crash location is visible even +when the script dies mid-way. +} +\examples{ +log_file <- tempfile(fileext = ".txt") +write_status(log_file, "INIT", "starting") +write_status(log_file, "DONE") +readLines(log_file) +} diff --git a/outputs/ci/phase0_n5_run26435798026/BIMODALITY_VERDICT.txt b/outputs/ci/phase0_n5_run26435798026/BIMODALITY_VERDICT.txt new file mode 100644 index 00000000..0fd62598 --- /dev/null +++ b/outputs/ci/phase0_n5_run26435798026/BIMODALITY_VERDICT.txt @@ -0,0 +1,19 @@ +PER-CHAIN MEDIANS: + Chain 1: median=+0.045 95% CrI [-0.599, +0.657] + Chain 2: median=+0.241 95% CrI [+0.200, +0.343] + (overall: Rhat=1.766 ESS_bulk=16) + +PAIRS PLOT INTERPRETATION: + Bundle stores only Omega_B[1,2] draws; M[2,k] / tau_B not available. + Hypothesis TRUE => opposite-sign chain medians; bimodal overall density + (two peaks straddling 0); per-chain density unimodal but at opposite modes. + At n=48: M[2,2] (log-boost, biomarker-2) anti-correlated with Omega_B[1,2] + in joint scatter is the defining feature of the sign-flip. + +VERDICT: WEAKLY BIMODAL + +CAVEAT: + n=5 is weakly informative; bimodality may show only partial chain separation + (Rhat 1.1-1.5). For n=48 look for: (i) bimodal Omega_B[1,2] marginal; + (ii) M[2,2] anti-correlated with Omega_B[1,2]; (iii) Rhat > 1.1 with + ESS_bulk < 200 despite adequate iteration count. diff --git a/outputs/ci/phase0_n5_run26435798026/PHASE0_STATUS.txt b/outputs/ci/phase0_n5_run26435798026/PHASE0_STATUS.txt new file mode 100644 index 00000000..52612165 --- /dev/null +++ b/outputs/ci/phase0_n5_run26435798026/PHASE0_STATUS.txt @@ -0,0 +1,14 @@ +[2026-05-26 06:21:12] STEP=INIT | Phase 0 started +[2026-05-26 06:21:12] STEP=LOAD_PACKAGES | logging +[2026-05-26 06:21:12] STEP=LOAD_PACKAGES | OK +[2026-05-26 06:21:12] STEP=COMPILE_DIR | setting up +[2026-05-26 06:21:12] STEP=COMPILE_DIR | OK +[2026-05-26 06:21:12] STEP=SIMULATE | running +[2026-05-26 06:21:13] STEP=SIMULATE | OK +[2026-05-26 06:21:13] STEP=FIT | running +[2026-05-26 06:24:09] STEP=FIT | OK +[2026-05-26 06:24:09] STEP=DIAG | extracting +[2026-05-26 06:24:09] STEP=DIAG | OK +[2026-05-26 06:24:09] STEP=SAVE | writing rds +[2026-05-26 06:24:09] STEP=SAVE | OK +[2026-05-26 06:24:09] STEP=DONE | Phase 0 completed successfully diff --git a/outputs/ci/phase0_n5_run26435798026/SUMMARY.txt b/outputs/ci/phase0_n5_run26435798026/SUMMARY.txt new file mode 100644 index 00000000..211a7140 --- /dev/null +++ b/outputs/ci/phase0_n5_run26435798026/SUMMARY.txt @@ -0,0 +1,16 @@ +RUN_ID: 26435798026 +N: 5 +ITER_WARMUP: 500 +ITER_SAMPLING: 500 +CHAINS: 2 +STATUS: OK +ELAPSED_MIN: 2.93 +TRUE_RHO_B: +0.600 +POST_MEDIAN_RHO_B: +0.230 +POST_LO_2.5: -0.516 +POST_HI_97.5: +0.603 +ESS_BULK: 16 +RHAT: 1.766 +DIVERGENT: 6 / 1000 +TREEDEPTH: 500 / 1000 +VERDICT: PATHOLOGICAL diff --git a/outputs/ci/phase0_n5_run26435798026/diagnostic_bimodality_pairs.png b/outputs/ci/phase0_n5_run26435798026/diagnostic_bimodality_pairs.png new file mode 100644 index 00000000..fb97fc5d Binary files /dev/null and b/outputs/ci/phase0_n5_run26435798026/diagnostic_bimodality_pairs.png differ diff --git a/outputs/ci/phase0_n5_run26435798026/env_versions.rds b/outputs/ci/phase0_n5_run26435798026/env_versions.rds new file mode 100644 index 00000000..56afc1e6 Binary files /dev/null and b/outputs/ci/phase0_n5_run26435798026/env_versions.rds differ diff --git a/outputs/ci/phase0_n5_run26435798026/one_fit_n5_ci.rds b/outputs/ci/phase0_n5_run26435798026/one_fit_n5_ci.rds new file mode 100644 index 00000000..d5d04cbd Binary files /dev/null and b/outputs/ci/phase0_n5_run26435798026/one_fit_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26435798026/one_fit_n5_ci_diag.rds b/outputs/ci/phase0_n5_run26435798026/one_fit_n5_ci_diag.rds new file mode 100644 index 00000000..e8240457 Binary files /dev/null and b/outputs/ci/phase0_n5_run26435798026/one_fit_n5_ci_diag.rds differ diff --git a/outputs/ci/phase0_n5_run26435798026/run.log b/outputs/ci/phase0_n5_run26435798026/run.log new file mode 100644 index 00000000..14a3e6fb --- /dev/null +++ b/outputs/ci/phase0_n5_run26435798026/run.log @@ -0,0 +1,189 @@ + +── PHASE 0: INTERACTIVE SLURM REPRODUCIBILITY TEST (n5_ci) ───────────────────── +Purpose: fit via salloc to compare determinism with Phase 1 sbatch +Started at: 2026-05-26 06:21:12 +Host: runnervmg397c + R R version 4.6.0 (2026-04-24) + cmdstanr 0.9.0.9000 + posterior 1.7.0 + shigella 0.0.0.9009 + serodynamics 0.0.0.9054 + cmdstan 2.38.0 + +compile_dir: /tmp/runner/cmdstan_bin_phase0_n5_ci +existing files: 0 +n_subjects: 5, rows: 50, true rho_B: 0.600 +ℹ Using Stan file: '/home/runner/work/_temp/Library/shigella/stan/model_2.stan' +ℹ Compile output directory: '/tmp/runner/cmdstan_bin_phase0_n5_ci' +ℹ Compiling model_2 (or using cache)... +ℹ Sampling model_2 with 2 chains... +Running MCMC with 2 parallel chains... + +Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 133, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 133, column 2 to column 43) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/Rtmpc5wkq7/model-1cf6798834be.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 1 finished in 51.4 seconds. +Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 2 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 2 finished in 153.1 seconds. + +Both chains finished successfully. +Mean chain execution time: 102.3 seconds. +Total execution time: 153.2 seconds. + +Warning: 6 of 1000 (1.0%) transitions ended with a divergence. +See https://mc-stan.org/misc/warnings for details. + +Warning: 500 of 1000 (50.0%) transitions hit the maximum treedepth limit of 12. +See https://mc-stan.org/misc/warnings for details. + +Fit elapsed: 2.93 min +Warning: 6 of 1000 (1.0%) transitions ended with a divergence. +See https://mc-stan.org/misc/warnings for details. + +Warning: 500 of 1000 (50.0%) transitions hit the maximum treedepth limit of 12. +See https://mc-stan.org/misc/warnings for details. + +Omega_B[1,2] posterior summary: +# A tibble: 1 × 8 + variable median mean sd `2.5%` `97.5%` ess_bulk rhat + +1 Omega_B[1,2] 0.230 0.155 0.264 -0.516 0.603 15.6 1.77 +saved -> outputs/ci/phase0_n5_run26435798026/one_fit_n5_ci.rds + +── PHASE 0 RESULT SUMMARY ────────────────────────────────────────────────────── +Status: OK +Elapsed: 2.93 min +True rho_B: +0.600 +Recovered median: +0.230 [-0.516, 0.603] +ESS_bulk: 16 +R-hat: 1.766 +Divergent: 6 / 1000 +Max-treedepth hits: 500 / 1000 + +── NEXT STEP ── + +1. Inspect outputs/ci/phase0_n5_run26435798026/one_fit_n5_ci.rds + +logs/phase0/*.log +2. If divergent rate <= 5% AND R-hat <= 1.01: +-> Proceed to Phase 1 (sbatch slurm/phase1_single.sbatch) +3. If divergent rate > 10% OR R-hat > 1.02: +-> Skip Phase 1-3, jump to Phase 4 diagnosis. +=== SUMMARY.txt === +RUN_ID: 26435798026 +N: 5 +ITER_WARMUP: 500 +ITER_SAMPLING: 500 +CHAINS: 2 +STATUS: OK +ELAPSED_MIN: 2.93 +TRUE_RHO_B: +0.600 +POST_MEDIAN_RHO_B: +0.230 +POST_LO_2.5: -0.516 +POST_HI_97.5: +0.603 +ESS_BULK: 16 +RHAT: 1.766 +DIVERGENT: 6 / 1000 +TREEDEPTH: 500 / 1000 +VERDICT: PATHOLOGICAL + +=== End SUMMARY.txt === diff --git a/outputs/ci/phase0_n5_run26435798026/sim_data_n5_ci.rds b/outputs/ci/phase0_n5_run26435798026/sim_data_n5_ci.rds new file mode 100644 index 00000000..88d35165 Binary files /dev/null and b/outputs/ci/phase0_n5_run26435798026/sim_data_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26470701379/BIMODALITY_VERDICT.txt b/outputs/ci/phase0_n5_run26470701379/BIMODALITY_VERDICT.txt new file mode 100644 index 00000000..986fdd50 --- /dev/null +++ b/outputs/ci/phase0_n5_run26470701379/BIMODALITY_VERDICT.txt @@ -0,0 +1,19 @@ +PER-CHAIN MEDIANS: + Chain 1: median=+0.025 95% CrI [-0.577, +0.729] + Chain 2: median=-0.078 95% CrI [-0.377, +0.128] + (overall: Rhat=1.161 ESS_bulk=12) + +PAIRS PLOT INTERPRETATION: + Bundle stores only Omega_B[1,2] draws; M[2,k] / tau_B not available. + Hypothesis TRUE => opposite-sign chain medians; bimodal overall density + (two peaks straddling 0); per-chain density unimodal but at opposite modes. + At n=48: M[2,2] (log-boost, biomarker-2) anti-correlated with Omega_B[1,2] + in joint scatter is the defining feature of the sign-flip. + +VERDICT: STRONGLY BIMODAL + +CAVEAT: + n=5 is weakly informative; bimodality may show only partial chain separation + (Rhat 1.1-1.5). For n=48 look for: (i) bimodal Omega_B[1,2] marginal; + (ii) M[2,2] anti-correlated with Omega_B[1,2]; (iii) Rhat > 1.1 with + ESS_bulk < 200 despite adequate iteration count. diff --git a/outputs/ci/phase0_n5_run26470701379/PHASE0_STATUS.txt b/outputs/ci/phase0_n5_run26470701379/PHASE0_STATUS.txt new file mode 100644 index 00000000..c73047a1 --- /dev/null +++ b/outputs/ci/phase0_n5_run26470701379/PHASE0_STATUS.txt @@ -0,0 +1,14 @@ +[2026-05-26 19:37:51] STEP=INIT | Phase 0 started +[2026-05-26 19:37:51] STEP=LOAD_PACKAGES | logging +[2026-05-26 19:37:52] STEP=LOAD_PACKAGES | OK +[2026-05-26 19:37:52] STEP=COMPILE_DIR | setting up +[2026-05-26 19:37:52] STEP=COMPILE_DIR | OK +[2026-05-26 19:37:52] STEP=SIMULATE | running +[2026-05-26 19:37:52] STEP=SIMULATE | OK +[2026-05-26 19:37:52] STEP=FIT | running +[2026-05-26 19:57:54] STEP=FIT | OK +[2026-05-26 19:57:54] STEP=DIAG | extracting +[2026-05-26 19:57:54] STEP=DIAG | OK +[2026-05-26 19:57:54] STEP=SAVE | writing rds +[2026-05-26 19:57:54] STEP=SAVE | OK +[2026-05-26 19:57:54] STEP=DONE | Phase 0 completed successfully diff --git a/outputs/ci/phase0_n5_run26470701379/SUMMARY.txt b/outputs/ci/phase0_n5_run26470701379/SUMMARY.txt new file mode 100644 index 00000000..2ff2fea4 --- /dev/null +++ b/outputs/ci/phase0_n5_run26470701379/SUMMARY.txt @@ -0,0 +1,16 @@ +RUN_ID: 26470701379 +N: 5 +ITER_WARMUP: 500 +ITER_SAMPLING: 500 +CHAINS: 2 +STATUS: OK +ELAPSED_MIN: 20.03 +TRUE_RHO_B: +0.600 +POST_MEDIAN_RHO_B: -0.043 +POST_LO_2.5: -0.523 +POST_HI_97.5: +0.677 +ESS_BULK: 12 +RHAT: 1.161 +DIVERGENT: 1 / 1000 +TREEDEPTH: 500 / 1000 +VERDICT: PATHOLOGICAL diff --git a/outputs/ci/phase0_n5_run26470701379/diagnostic_bimodality.log b/outputs/ci/phase0_n5_run26470701379/diagnostic_bimodality.log new file mode 100644 index 00000000..9f89a9b0 --- /dev/null +++ b/outputs/ci/phase0_n5_run26470701379/diagnostic_bimodality.log @@ -0,0 +1,5 @@ +null device + 1 +Pairs plot: outputs/ci/phase0_n5_run26470701379/diagnostic_bimodality_pairs.png +Verdict: outputs/ci/phase0_n5_run26470701379/BIMODALITY_VERDICT.txt +VERDICT: STRONGLY BIMODAL diff --git a/outputs/ci/phase0_n5_run26470701379/diagnostic_bimodality_pairs.png b/outputs/ci/phase0_n5_run26470701379/diagnostic_bimodality_pairs.png new file mode 100644 index 00000000..075a2879 Binary files /dev/null and b/outputs/ci/phase0_n5_run26470701379/diagnostic_bimodality_pairs.png differ diff --git a/outputs/ci/phase0_n5_run26470701379/env_versions.rds b/outputs/ci/phase0_n5_run26470701379/env_versions.rds new file mode 100644 index 00000000..56afc1e6 Binary files /dev/null and b/outputs/ci/phase0_n5_run26470701379/env_versions.rds differ diff --git a/outputs/ci/phase0_n5_run26470701379/one_fit_n5_ci.rds b/outputs/ci/phase0_n5_run26470701379/one_fit_n5_ci.rds new file mode 100644 index 00000000..6a1d58fe Binary files /dev/null and b/outputs/ci/phase0_n5_run26470701379/one_fit_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26470701379/one_fit_n5_ci_diag.rds b/outputs/ci/phase0_n5_run26470701379/one_fit_n5_ci_diag.rds new file mode 100644 index 00000000..871fc3d0 Binary files /dev/null and b/outputs/ci/phase0_n5_run26470701379/one_fit_n5_ci_diag.rds differ diff --git a/outputs/ci/phase0_n5_run26470701379/run.log b/outputs/ci/phase0_n5_run26470701379/run.log new file mode 100644 index 00000000..0b640ff1 --- /dev/null +++ b/outputs/ci/phase0_n5_run26470701379/run.log @@ -0,0 +1,259 @@ + +── PHASE 0: INTERACTIVE SLURM REPRODUCIBILITY TEST (n5_ci) ───────────────────── +Purpose: fit via salloc to compare determinism with Phase 1 sbatch +Started at: 2026-05-26 19:37:51 +Host: runnervmg397c + R R version 4.6.0 (2026-04-24) + cmdstanr 0.9.0.9000 + posterior 1.7.0 + shigella 0.0.0.9009 + serodynamics 0.0.0.9054 + cmdstan 2.38.0 + +compile_dir: /tmp/runner/cmdstan_bin_phase0_n5_ci +existing files: 0 +n_subjects: 5, rows: 50, true rho_B: 0.600 +ℹ Using Stan file: '/home/runner/work/_temp/Library/shigella/stan/model_2.stan' +ℹ Compile output directory: '/tmp/runner/cmdstan_bin_phase0_n5_ci' +ℹ Compiling model_2 (or using cache)... +ℹ Sampling model_2 with 2 chains... +Running MCMC with 2 parallel chains... + +Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 133, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 133, column 2 to column 43) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 1 finished in 67.1 seconds. +Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpxuPwIK/model-1a853798fed1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 2 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 2 finished in 1175.5 seconds. + +Both chains finished successfully. +Mean chain execution time: 621.3 seconds. +Total execution time: 1175.6 seconds. + +Warning: 1 of 1000 (0.0%) transitions ended with a divergence. +See https://mc-stan.org/misc/warnings for details. + +Warning: 500 of 1000 (50.0%) transitions hit the maximum treedepth limit of 15. +See https://mc-stan.org/misc/warnings for details. + +Fit elapsed: 20.03 min +Warning: 1 of 1000 (0.0%) transitions ended with a divergence. +See https://mc-stan.org/misc/warnings for details. + +Warning: 500 of 1000 (50.0%) transitions hit the maximum treedepth limit of 15. +See https://mc-stan.org/misc/warnings for details. + +Omega_B[1,2] posterior summary: +# A tibble: 1 × 8 + variable median mean sd `2.5%` `97.5%` ess_bulk rhat + +1 Omega_B[1,2] -0.0430 -0.0230 0.279 -0.523 0.677 11.5 1.16 +saved -> outputs/ci/phase0_n5_run26470701379/one_fit_n5_ci.rds + +── PHASE 0 RESULT SUMMARY ────────────────────────────────────────────────────── +Status: OK +Elapsed: 20.03 min +True rho_B: +0.600 +Recovered median: -0.043 [-0.523, 0.677] +ESS_bulk: 12 +R-hat: 1.161 +Divergent: 1 / 1000 +Max-treedepth hits: 500 / 1000 + +── NEXT STEP ── + +1. Inspect outputs/ci/phase0_n5_run26470701379/one_fit_n5_ci.rds + +logs/phase0/*.log +2. If divergent rate <= 5% AND R-hat <= 1.01: +-> Proceed to Phase 1 (sbatch slurm/phase1_single.sbatch) +3. If divergent rate > 10% OR R-hat > 1.02: +-> Skip Phase 1-3, jump to Phase 4 diagnosis. +=== SUMMARY.txt === +RUN_ID: 26470701379 +N: 5 +ITER_WARMUP: 500 +ITER_SAMPLING: 500 +CHAINS: 2 +STATUS: OK +ELAPSED_MIN: 20.03 +TRUE_RHO_B: +0.600 +POST_MEDIAN_RHO_B: -0.043 +POST_LO_2.5: -0.523 +POST_HI_97.5: +0.677 +ESS_BULK: 12 +RHAT: 1.161 +DIVERGENT: 1 / 1000 +TREEDEPTH: 500 / 1000 +VERDICT: PATHOLOGICAL + +=== End SUMMARY.txt === diff --git a/outputs/ci/phase0_n5_run26470701379/sim_data_n5_ci.rds b/outputs/ci/phase0_n5_run26470701379/sim_data_n5_ci.rds new file mode 100644 index 00000000..88d35165 Binary files /dev/null and b/outputs/ci/phase0_n5_run26470701379/sim_data_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26549772922/BIMODALITY_VERDICT.txt b/outputs/ci/phase0_n5_run26549772922/BIMODALITY_VERDICT.txt new file mode 100644 index 00000000..a7d6e9f7 --- /dev/null +++ b/outputs/ci/phase0_n5_run26549772922/BIMODALITY_VERDICT.txt @@ -0,0 +1,19 @@ +PER-CHAIN MEDIANS: + Chain 1: median=+0.389 95% CrI [-0.107, +0.695] + Chain 2: median=+0.296 95% CrI [+0.102, +0.461] + (overall: Rhat=1.562 ESS_bulk=7) + +PAIRS PLOT INTERPRETATION: + Bundle stores only Omega_B[1,2] draws; M[2,k] / tau_B not available. + Hypothesis TRUE => opposite-sign chain medians; bimodal overall density + (two peaks straddling 0); per-chain density unimodal but at opposite modes. + At n=48: M[2,2] (log-boost, biomarker-2) anti-correlated with Omega_B[1,2] + in joint scatter is the defining feature of the sign-flip. + +VERDICT: WEAKLY BIMODAL + +CAVEAT: + n=5 is weakly informative; bimodality may show only partial chain separation + (Rhat 1.1-1.5). For n=48 look for: (i) bimodal Omega_B[1,2] marginal; + (ii) M[2,2] anti-correlated with Omega_B[1,2]; (iii) Rhat > 1.1 with + ESS_bulk < 200 despite adequate iteration count. diff --git a/outputs/ci/phase0_n5_run26549772922/PHASE0_STATUS.txt b/outputs/ci/phase0_n5_run26549772922/PHASE0_STATUS.txt new file mode 100644 index 00000000..53288350 --- /dev/null +++ b/outputs/ci/phase0_n5_run26549772922/PHASE0_STATUS.txt @@ -0,0 +1,14 @@ +[2026-05-28 01:54:06] STEP=INIT | Phase 0 started +[2026-05-28 01:54:06] STEP=LOAD_PACKAGES | logging +[2026-05-28 01:54:06] STEP=LOAD_PACKAGES | OK +[2026-05-28 01:54:06] STEP=COMPILE_DIR | setting up +[2026-05-28 01:54:06] STEP=COMPILE_DIR | OK +[2026-05-28 01:54:06] STEP=SIMULATE | running +[2026-05-28 01:54:07] STEP=SIMULATE | OK +[2026-05-28 01:54:07] STEP=FIT | running +[2026-05-28 02:13:58] STEP=FIT | OK +[2026-05-28 02:13:58] STEP=DIAG | extracting +[2026-05-28 02:13:58] STEP=DIAG | OK +[2026-05-28 02:13:58] STEP=SAVE | writing rds +[2026-05-28 02:13:58] STEP=SAVE | OK +[2026-05-28 02:13:58] STEP=DONE | Phase 0 completed successfully diff --git a/outputs/ci/phase0_n5_run26549772922/SUMMARY.txt b/outputs/ci/phase0_n5_run26549772922/SUMMARY.txt new file mode 100644 index 00000000..50fdf167 --- /dev/null +++ b/outputs/ci/phase0_n5_run26549772922/SUMMARY.txt @@ -0,0 +1,16 @@ +RUN_ID: 26549772922 +N: 5 +ITER_WARMUP: 500 +ITER_SAMPLING: 500 +CHAINS: 2 +STATUS: OK +ELAPSED_MIN: 19.85 +TRUE_RHO_B: +0.600 +POST_MEDIAN_RHO_B: +0.305 +POST_LO_2.5: -0.056 +POST_HI_97.5: +0.675 +ESS_BULK: 7 +RHAT: 1.562 +DIVERGENT: 1 / 1000 +TREEDEPTH: 999 / 1000 +VERDICT: PATHOLOGICAL diff --git a/outputs/ci/phase0_n5_run26549772922/diagnostic_bimodality.log b/outputs/ci/phase0_n5_run26549772922/diagnostic_bimodality.log new file mode 100644 index 00000000..ac442403 --- /dev/null +++ b/outputs/ci/phase0_n5_run26549772922/diagnostic_bimodality.log @@ -0,0 +1,5 @@ +null device + 1 +Pairs plot: outputs/ci/phase0_n5_run26549772922/diagnostic_bimodality_pairs.png +Verdict: outputs/ci/phase0_n5_run26549772922/BIMODALITY_VERDICT.txt +VERDICT: WEAKLY BIMODAL diff --git a/outputs/ci/phase0_n5_run26549772922/diagnostic_bimodality_pairs.png b/outputs/ci/phase0_n5_run26549772922/diagnostic_bimodality_pairs.png new file mode 100644 index 00000000..24e988d2 Binary files /dev/null and b/outputs/ci/phase0_n5_run26549772922/diagnostic_bimodality_pairs.png differ diff --git a/outputs/ci/phase0_n5_run26549772922/env_versions.rds b/outputs/ci/phase0_n5_run26549772922/env_versions.rds new file mode 100644 index 00000000..cd486c87 Binary files /dev/null and b/outputs/ci/phase0_n5_run26549772922/env_versions.rds differ diff --git a/outputs/ci/phase0_n5_run26549772922/one_fit_n5_ci.rds b/outputs/ci/phase0_n5_run26549772922/one_fit_n5_ci.rds new file mode 100644 index 00000000..ba536a03 Binary files /dev/null and b/outputs/ci/phase0_n5_run26549772922/one_fit_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26549772922/one_fit_n5_ci_diag.rds b/outputs/ci/phase0_n5_run26549772922/one_fit_n5_ci_diag.rds new file mode 100644 index 00000000..c898277f Binary files /dev/null and b/outputs/ci/phase0_n5_run26549772922/one_fit_n5_ci_diag.rds differ diff --git a/outputs/ci/phase0_n5_run26549772922/run.log b/outputs/ci/phase0_n5_run26549772922/run.log new file mode 100644 index 00000000..af9941e8 --- /dev/null +++ b/outputs/ci/phase0_n5_run26549772922/run.log @@ -0,0 +1,194 @@ + +── PHASE 0: INTERACTIVE SLURM REPRODUCIBILITY TEST (n5_ci) ───────────────────── +Purpose: fit via salloc to compare determinism with Phase 1 sbatch +Started at: 2026-05-28 01:54:06 +Host: runnervm3jyl0 + R R version 4.6.0 (2026-04-24) + cmdstanr 0.9.0.9000 + posterior 1.7.0 + shigella 0.0.0.9009 + serodynamics 0.0.0.9055 + cmdstan 2.38.0 + +compile_dir: /tmp/runner/cmdstan_bin_phase0_n5_ci +existing files: 0 +n_subjects: 5, rows: 50, true rho_B: 0.600 +ℹ Using Stan file: '/home/runner/work/_temp/Library/shigella/stan/model_2.stan' +ℹ Compile output directory: '/tmp/runner/cmdstan_bin_phase0_n5_ci' +ℹ Compiling model_2 (or using cache)... +ℹ Sampling model_2 with 2 chains... +Running MCMC with 2 parallel chains... + +Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 133, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 133, column 2 to column 43) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -nan, but must be finite! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -nan, but must be finite! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -nan, but must be finite! (in '/tmp/Rtmp23ZyS6/model-1b6c2f09fce1.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 2 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 2 finished in 1160.3 seconds. +Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 1 finished in 1166.5 seconds. + +Both chains finished successfully. +Mean chain execution time: 1163.4 seconds. +Total execution time: 1166.6 seconds. + +Warning: 1 of 1000 (0.0%) transitions ended with a divergence. +See https://mc-stan.org/misc/warnings for details. + +Warning: 999 of 1000 (100.0%) transitions hit the maximum treedepth limit of 15. +See https://mc-stan.org/misc/warnings for details. + +Fit elapsed: 19.85 min +Warning: 1 of 1000 (0.0%) transitions ended with a divergence. +See https://mc-stan.org/misc/warnings for details. + +Warning: 999 of 1000 (100.0%) transitions hit the maximum treedepth limit of 15. +See https://mc-stan.org/misc/warnings for details. + +Omega_B[1,2] posterior summary: +# A tibble: 1 × 8 + variable median mean sd `2.5%` `97.5%` ess_bulk rhat + +1 Omega_B[1,2] 0.305 0.318 0.181 -0.0557 0.675 7.21 1.56 +saved -> outputs/ci/phase0_n5_run26549772922/one_fit_n5_ci.rds + +── PHASE 0 RESULT SUMMARY ────────────────────────────────────────────────────── +Status: OK +Elapsed: 19.85 min +True rho_B: +0.600 +Recovered median: +0.305 [-0.056, 0.675] +ESS_bulk: 7 +R-hat: 1.562 +Divergent: 1 / 1000 +Max-treedepth hits: 999 / 1000 + +── NEXT STEP ── + +1. Inspect outputs/ci/phase0_n5_run26549772922/one_fit_n5_ci.rds + +logs/phase0/*.log +2. If divergent rate <= 5% AND R-hat <= 1.01: +-> Proceed to Phase 1 (sbatch slurm/phase1_single.sbatch) +3. If divergent rate > 10% OR R-hat > 1.02: +-> Skip Phase 1-3, jump to Phase 4 diagnosis. +=== SUMMARY.txt === +RUN_ID: 26549772922 +N: 5 +ITER_WARMUP: 500 +ITER_SAMPLING: 500 +CHAINS: 2 +STATUS: OK +ELAPSED_MIN: 19.85 +TRUE_RHO_B: +0.600 +POST_MEDIAN_RHO_B: +0.305 +POST_LO_2.5: -0.056 +POST_HI_97.5: +0.675 +ESS_BULK: 7 +RHAT: 1.562 +DIVERGENT: 1 / 1000 +TREEDEPTH: 999 / 1000 +VERDICT: PATHOLOGICAL + +=== End SUMMARY.txt === diff --git a/outputs/ci/phase0_n5_run26549772922/sim_data_n5_ci.rds b/outputs/ci/phase0_n5_run26549772922/sim_data_n5_ci.rds new file mode 100644 index 00000000..88d35165 Binary files /dev/null and b/outputs/ci/phase0_n5_run26549772922/sim_data_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26551575231/PHASE0_STATUS.txt b/outputs/ci/phase0_n5_run26551575231/PHASE0_STATUS.txt new file mode 100644 index 00000000..1ffa5cad --- /dev/null +++ b/outputs/ci/phase0_n5_run26551575231/PHASE0_STATUS.txt @@ -0,0 +1,8 @@ +[2026-05-28 02:50:39] STEP=INIT | Phase 0 started +[2026-05-28 02:50:39] STEP=LOAD_PACKAGES | logging +[2026-05-28 02:50:40] STEP=LOAD_PACKAGES | OK +[2026-05-28 02:50:40] STEP=COMPILE_DIR | setting up +[2026-05-28 02:50:40] STEP=COMPILE_DIR | OK +[2026-05-28 02:50:40] STEP=SIMULATE | running +[2026-05-28 02:50:40] STEP=SIMULATE | OK +[2026-05-28 02:50:40] STEP=FIT | running diff --git a/outputs/ci/phase0_n5_run26551575231/diagnostic_bimodality.log b/outputs/ci/phase0_n5_run26551575231/diagnostic_bimodality.log new file mode 100644 index 00000000..a86df8e2 --- /dev/null +++ b/outputs/ci/phase0_n5_run26551575231/diagnostic_bimodality.log @@ -0,0 +1,2 @@ +Error: Missing SUMMARY.txt: outputs/ci/phase0_n5_run26551575231/SUMMARY.txt +Execution halted diff --git a/outputs/ci/phase0_n5_run26551575231/env_versions.rds b/outputs/ci/phase0_n5_run26551575231/env_versions.rds new file mode 100644 index 00000000..cd486c87 Binary files /dev/null and b/outputs/ci/phase0_n5_run26551575231/env_versions.rds differ diff --git a/outputs/ci/phase0_n5_run26551575231/one_fit_n5_ci.rds b/outputs/ci/phase0_n5_run26551575231/one_fit_n5_ci.rds new file mode 100644 index 00000000..153c32ba Binary files /dev/null and b/outputs/ci/phase0_n5_run26551575231/one_fit_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26551575231/run.log b/outputs/ci/phase0_n5_run26551575231/run.log new file mode 100644 index 00000000..18a89d9e --- /dev/null +++ b/outputs/ci/phase0_n5_run26551575231/run.log @@ -0,0 +1,116 @@ + +── PHASE 0: INTERACTIVE SLURM REPRODUCIBILITY TEST (n5_ci) ───────────────────── +Purpose: fit via salloc to compare determinism with Phase 1 sbatch +Started at: 2026-05-28 02:50:39 +Host: runnervm3jyl0 + R R version 4.6.0 (2026-04-24) + cmdstanr 0.9.0.9000 + posterior 1.7.0 + shigella 0.0.0.9009 + serodynamics 0.0.0.9055 + cmdstan 2.38.0 + +compile_dir: /tmp/runner/cmdstan_bin_phase0_n5_ci +existing files: 0 +n_subjects: 5, rows: 50, true rho_B: 0.600 +ℹ Using Stan file: '/home/runner/work/_temp/Library/shigella/stan/model_2.stan' +ℹ Compile output directory: '/tmp/runner/cmdstan_bin_phase0_n5_ci' +ℹ Compiling model_2 (or using cache)... +ℹ Sampling model_2 with 2 chains... +Running MCMC with 2 parallel chains... + +Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 133, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 133, column 2 to column 43) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[3] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 134, column 2 to column 43) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[3] is 0, but must be positive! (in '/tmp/RtmptqsCBe/model-1a91408b3c73.stan', line 134, column 2 to column 43) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 1 finished in 9472.1 seconds. +Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling) diff --git a/outputs/ci/phase0_n5_run26551575231/sim_data_n5_ci.rds b/outputs/ci/phase0_n5_run26551575231/sim_data_n5_ci.rds new file mode 100644 index 00000000..88d35165 Binary files /dev/null and b/outputs/ci/phase0_n5_run26551575231/sim_data_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26589037994/BIMODALITY_VERDICT.txt b/outputs/ci/phase0_n5_run26589037994/BIMODALITY_VERDICT.txt new file mode 100644 index 00000000..0716f08c --- /dev/null +++ b/outputs/ci/phase0_n5_run26589037994/BIMODALITY_VERDICT.txt @@ -0,0 +1,19 @@ +PER-CHAIN MEDIANS: + Chain 1: median=+0.183 95% CrI [+0.047, +0.310] + Chain 2: median=+0.038 95% CrI [-0.425, +0.570] + (overall: Rhat=1.355 ESS_bulk=22) + +PAIRS PLOT INTERPRETATION: + Bundle stores only Omega_B[1,2] draws; M[2,k] / tau_B not available. + Hypothesis TRUE => opposite-sign chain medians; bimodal overall density + (two peaks straddling 0); per-chain density unimodal but at opposite modes. + At n=48: M[2,2] (log-boost, biomarker-2) anti-correlated with Omega_B[1,2] + in joint scatter is the defining feature of the sign-flip. + +VERDICT: WEAKLY BIMODAL + +CAVEAT: + n=5 is weakly informative; bimodality may show only partial chain separation + (Rhat 1.1-1.5). For n=48 look for: (i) bimodal Omega_B[1,2] marginal; + (ii) M[2,2] anti-correlated with Omega_B[1,2]; (iii) Rhat > 1.1 with + ESS_bulk < 200 despite adequate iteration count. diff --git a/outputs/ci/phase0_n5_run26589037994/PHASE0_STATUS.txt b/outputs/ci/phase0_n5_run26589037994/PHASE0_STATUS.txt new file mode 100644 index 00000000..9f66b8d3 --- /dev/null +++ b/outputs/ci/phase0_n5_run26589037994/PHASE0_STATUS.txt @@ -0,0 +1,14 @@ +[2026-05-28 16:53:42] STEP=INIT | Phase 0 started +[2026-05-28 16:53:42] STEP=LOAD_PACKAGES | logging +[2026-05-28 16:53:42] STEP=LOAD_PACKAGES | OK +[2026-05-28 16:53:42] STEP=COMPILE_DIR | setting up +[2026-05-28 16:53:42] STEP=COMPILE_DIR | OK +[2026-05-28 16:53:42] STEP=SIMULATE | running +[2026-05-28 16:53:43] STEP=SIMULATE | OK +[2026-05-28 16:53:43] STEP=FIT | running +[2026-05-28 17:12:46] STEP=FIT | OK +[2026-05-28 17:12:46] STEP=DIAG | extracting +[2026-05-28 17:12:46] STEP=DIAG | OK +[2026-05-28 17:12:46] STEP=SAVE | writing rds +[2026-05-28 17:12:46] STEP=SAVE | OK +[2026-05-28 17:12:46] STEP=DONE | Phase 0 completed successfully diff --git a/outputs/ci/phase0_n5_run26589037994/SUMMARY.txt b/outputs/ci/phase0_n5_run26589037994/SUMMARY.txt new file mode 100644 index 00000000..e2a4c04d --- /dev/null +++ b/outputs/ci/phase0_n5_run26589037994/SUMMARY.txt @@ -0,0 +1,16 @@ +RUN_ID: 26589037994 +N: 5 +ITER_WARMUP: 500 +ITER_SAMPLING: 500 +CHAINS: 2 +STATUS: OK +ELAPSED_MIN: 19.04 +TRUE_RHO_B: +0.600 +POST_MEDIAN_RHO_B: +0.154 +POST_LO_2.5: -0.339 +POST_HI_97.5: +0.509 +ESS_BULK: 22 +RHAT: 1.355 +DIVERGENT: 12 / 1000 +TREEDEPTH: 500 / 1000 +VERDICT: PATHOLOGICAL diff --git a/outputs/ci/phase0_n5_run26589037994/diagnostic_bimodality.log b/outputs/ci/phase0_n5_run26589037994/diagnostic_bimodality.log new file mode 100644 index 00000000..07d4be1e --- /dev/null +++ b/outputs/ci/phase0_n5_run26589037994/diagnostic_bimodality.log @@ -0,0 +1,5 @@ +null device + 1 +Pairs plot: outputs/ci/phase0_n5_run26589037994/diagnostic_bimodality_pairs.png +Verdict: outputs/ci/phase0_n5_run26589037994/BIMODALITY_VERDICT.txt +VERDICT: WEAKLY BIMODAL diff --git a/outputs/ci/phase0_n5_run26589037994/diagnostic_bimodality_pairs.png b/outputs/ci/phase0_n5_run26589037994/diagnostic_bimodality_pairs.png new file mode 100644 index 00000000..486cfee9 Binary files /dev/null and b/outputs/ci/phase0_n5_run26589037994/diagnostic_bimodality_pairs.png differ diff --git a/outputs/ci/phase0_n5_run26589037994/env_versions.rds b/outputs/ci/phase0_n5_run26589037994/env_versions.rds new file mode 100644 index 00000000..cd486c87 Binary files /dev/null and b/outputs/ci/phase0_n5_run26589037994/env_versions.rds differ diff --git a/outputs/ci/phase0_n5_run26589037994/one_fit_n5_ci.rds b/outputs/ci/phase0_n5_run26589037994/one_fit_n5_ci.rds new file mode 100644 index 00000000..5fac6e06 Binary files /dev/null and b/outputs/ci/phase0_n5_run26589037994/one_fit_n5_ci.rds differ diff --git a/outputs/ci/phase0_n5_run26589037994/one_fit_n5_ci_diag.rds b/outputs/ci/phase0_n5_run26589037994/one_fit_n5_ci_diag.rds new file mode 100644 index 00000000..6ec65dc4 Binary files /dev/null and b/outputs/ci/phase0_n5_run26589037994/one_fit_n5_ci_diag.rds differ diff --git a/outputs/ci/phase0_n5_run26589037994/run.log b/outputs/ci/phase0_n5_run26589037994/run.log new file mode 100644 index 00000000..e338c04f --- /dev/null +++ b/outputs/ci/phase0_n5_run26589037994/run.log @@ -0,0 +1,179 @@ + +── PHASE 0: INTERACTIVE SLURM REPRODUCIBILITY TEST (n5_ci) ───────────────────── +Purpose: fit via salloc to compare determinism with Phase 1 sbatch +Started at: 2026-05-28 16:53:42 +Host: runnervm3jyl0 + R R version 4.6.0 (2026-04-24) + cmdstanr 0.9.0.9000 + posterior 1.7.0 + shigella 0.0.0.9009 + serodynamics 0.0.0.9055 + cmdstan 2.38.0 + +compile_dir: /tmp/runner/cmdstan_bin_phase0_n5_ci +existing files: 0 +n_subjects: 5, rows: 50, true rho_B: 0.600 +ℹ Using Stan file: '/home/runner/work/_temp/Library/shigella/stan/model_2.stan' +ℹ Compile output directory: '/tmp/runner/cmdstan_bin_phase0_n5_ci' +ℹ Compiling model_2 (or using cache)... +ℹ Sampling model_2 with 2 chains... +Running MCMC with 2 parallel chains... + +Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[5] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 134, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 133, column 2 to column 43) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 135, column 2 to column 47) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 1 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -nan, but must be finite! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 171, column 8 to column 63) +Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 1 +Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup) +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 133, column 2 to column 43) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 135, column 2 to column 47) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: +Chain 2 Exception: multi_normal_cholesky_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpxU8N3U/model-1aab6244a658.stan', line 171, column 8 to column 63) +Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, +Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified. +Chain 2 +Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 2 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 2 finished in 45.4 seconds. +Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup) +Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup) +Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) +Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) +Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) +Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) +Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) +Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) +Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) +Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) +Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) +Chain 1 finished in 1119.0 seconds. + +Both chains finished successfully. +Mean chain execution time: 582.2 seconds. +Total execution time: 1119.0 seconds. + +Warning: 12 of 1000 (1.0%) transitions ended with a divergence. +See https://mc-stan.org/misc/warnings for details. + +Warning: 500 of 1000 (50.0%) transitions hit the maximum treedepth limit of 15. +See https://mc-stan.org/misc/warnings for details. + +Fit elapsed: 19.04 min +Warning: 12 of 1000 (1.0%) transitions ended with a divergence. +See https://mc-stan.org/misc/warnings for details. + +Warning: 500 of 1000 (50.0%) transitions hit the maximum treedepth limit of 15. +See https://mc-stan.org/misc/warnings for details. + +Omega_B[1,2] posterior summary: +# A tibble: 1 × 8 + variable median mean sd `2.5%` `97.5%` ess_bulk rhat + +1 Omega_B[1,2] 0.154 0.123 0.204 -0.339 0.509 21.6 1.35 +saved -> outputs/ci/phase0_n5_run26589037994/one_fit_n5_ci.rds + +── PHASE 0 RESULT SUMMARY ────────────────────────────────────────────────────── +Status: OK +Elapsed: 19.04 min +True rho_B: +0.600 +Recovered median: +0.154 [-0.339, 0.509] +ESS_bulk: 22 +R-hat: 1.355 +Divergent: 12 / 1000 +Max-treedepth hits: 500 / 1000 + +── NEXT STEP ── + +1. Inspect outputs/ci/phase0_n5_run26589037994/one_fit_n5_ci.rds + +logs/phase0/*.log +2. If divergent rate <= 5% AND R-hat <= 1.01: +-> Proceed to Phase 1 (sbatch slurm/phase1_single.sbatch) +3. If divergent rate > 10% OR R-hat > 1.02: +-> Skip Phase 1-3, jump to Phase 4 diagnosis. +=== SUMMARY.txt === +RUN_ID: 26589037994 +N: 5 +ITER_WARMUP: 500 +ITER_SAMPLING: 500 +CHAINS: 2 +STATUS: OK +ELAPSED_MIN: 19.04 +TRUE_RHO_B: +0.600 +POST_MEDIAN_RHO_B: +0.154 +POST_LO_2.5: -0.339 +POST_HI_97.5: +0.509 +ESS_BULK: 22 +RHAT: 1.355 +DIVERGENT: 12 / 1000 +TREEDEPTH: 500 / 1000 +VERDICT: PATHOLOGICAL + +=== End SUMMARY.txt === diff --git a/outputs/ci/phase0_n5_run26589037994/sim_data_n5_ci.rds b/outputs/ci/phase0_n5_run26589037994/sim_data_n5_ci.rds new file mode 100644 index 00000000..88d35165 Binary files /dev/null and b/outputs/ci/phase0_n5_run26589037994/sim_data_n5_ci.rds differ diff --git a/scripts/diagnostic_bimodality.R b/scripts/diagnostic_bimodality.R new file mode 100644 index 00000000..eadd6ee3 --- /dev/null +++ b/scripts/diagnostic_bimodality.R @@ -0,0 +1,93 @@ +#!/usr/bin/env Rscript +# Bimodality / sign-label-ambiguity diagnostic for Omega_B[1,2] +# Usage: Rscript scripts/diagnostic_bimodality.R +# out_dir: run output directory containing SUMMARY.txt and one_fit_n_ci.rds +# If omitted, falls back to hardcoded n=5 path for local development. + +args <- commandArgs(trailingOnly = TRUE) +if (length(args) >= 1L) { + out_dir <- args[1L] +} else { + out_dir <- "outputs/ci/phase0_n5_run26435798026" +} + +summary_path <- file.path(out_dir, "SUMMARY.txt") +if (!file.exists(summary_path)) stop("Missing SUMMARY.txt: ", summary_path) + +lines <- readLines(summary_path) +n_line <- grep("^N:", lines, value = TRUE) +if (length(n_line) == 0L) stop("N: line not found in SUMMARY.txt") +N <- as.integer(sub("^N:\\s*", "", n_line[1L])) +if (is.na(N)) stop("Could not parse N from SUMMARY.txt") + +rds_path <- file.path(out_dir, sprintf("one_fit_n%d_ci.rds", N)) +if (!file.exists(rds_path)) stop("Missing RDS: ", rds_path) + +bundle <- readRDS(rds_path) +rho_all <- bundle$rho_B_posterior +n_chains <- bundle$fit_settings$chains +n_iter <- bundle$fit_settings$sampling + +if (is.null(rho_all)) stop("rho_B_posterior is NULL (fit may have crashed)") +if (length(rho_all) != n_chains * n_iter) + stop(sprintf("Expected %d draws, got %d", n_chains * n_iter, length(rho_all))) + +chain_draws <- lapply(seq_len(n_chains), + function(ch) rho_all[((ch-1)*n_iter+1):(ch*n_iter)]) +cq <- function(x, p) quantile(x, p, names=FALSE) +cstats <- lapply(chain_draws, function(x) c(med=median(x), lo=cq(x,.025), hi=cq(x,.975))) + +png_path <- file.path(out_dir, "diagnostic_bimodality_pairs.png") +png(png_path, width=900, height=420, res=120) +par(mfrow=c(1L, n_chains+1L), mar=c(4,4,3,1)) +plot(density(rho_all), main="Omega_B[1,2] all chains", + xlab=expression(rho[B]), lwd=2); abline(v=0, lty=2, col="grey60") +cols <- c("steelblue","tomato","forestgreen","goldenrod")[seq_len(n_chains)] +for (ch in seq_len(n_chains)) { + plot(density(chain_draws[[ch]]), + main=sprintf("Chain %d med=%+.3f", ch, cstats[[ch]]["med"]), + xlab=expression(rho[B]), col=cols[ch], lwd=2) + abline(v=0, lty=2, col="grey60") +} +dev.off() +cat("Pairs plot:", png_path, "\n") + +signs <- sapply(cstats, function(s) sign(s["med"])) +med_range <- diff(range(sapply(cstats, function(s) s["med"]))) +rhat <- bundle$omega_B_summary$rhat +ess <- bundle$omega_B_summary$ess_bulk +verdict <- if (length(unique(signs)) > 1L) "STRONGLY BIMODAL" else + if (med_range > .40 || (!is.null(rhat) && rhat > 1.20)) "WEAKLY BIMODAL" else + if (!is.null(rhat) && rhat < 1.10 && med_range < .20) "UNIMODAL" else + "INSUFFICIENT EVIDENCE" + +chain_lines <- sapply(seq_len(n_chains), function(ch) { + s <- cstats[[ch]] + sprintf(" Chain %d: median=%+.3f 95%% CrI [%+.3f, %+.3f]", ch, s["med"], s["lo"], s["hi"]) +}) + +txt <- c( + "PER-CHAIN MEDIANS:", + chain_lines, + sprintf(" (overall: Rhat=%.3f ESS_bulk=%.0f)", rhat, ess), + "", + "PAIRS PLOT INTERPRETATION:", + " Bundle stores only Omega_B[1,2] draws; M[2,k] / tau_B not available.", + " Hypothesis TRUE => opposite-sign chain medians; bimodal overall density", + " (two peaks straddling 0); per-chain density unimodal but at opposite modes.", + " At n=48: M[2,2] (log-boost, biomarker-2) anti-correlated with Omega_B[1,2]", + " in joint scatter is the defining feature of the sign-flip.", + "", + paste("VERDICT:", verdict), + "", + "CAVEAT:", + " n=5 is weakly informative; bimodality may show only partial chain separation", + " (Rhat 1.1-1.5). For n=48 look for: (i) bimodal Omega_B[1,2] marginal;", + " (ii) M[2,2] anti-correlated with Omega_B[1,2]; (iii) Rhat > 1.1 with", + " ESS_bulk < 200 despite adequate iteration count." +) + +verdict_path <- file.path(out_dir, "BIMODALITY_VERDICT.txt") +writeLines(txt, verdict_path) +cat("Verdict:", verdict_path, "\n") +cat("VERDICT:", verdict, "\n") diff --git a/scripts/phase0_interactive_reproducibility.R b/scripts/phase0_interactive_reproducibility.R new file mode 100644 index 00000000..c604f212 --- /dev/null +++ b/scripts/phase0_interactive_reproducibility.R @@ -0,0 +1,27 @@ +# ============================================================================ +# phase0_interactive_reproducibility.R — n = 5 pilot +# +# Execution: +# salloc --time=04:00:00 --cpus-per-task=2 --mem=10G +# Rscript scripts/phase0_interactive_reproducibility.R +# exit +# +# Purpose: reproduce a Phase 1 (sbatch) fit under interactive SLURM to confirm +# determinism across allocation modes. Run the n=5 version first as a fast +# smoke test before committing to the n=48 full-cohort run. +# ============================================================================ +local({ + flag <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE) + if (length(flag)) { + root <- dirname(dirname(normalizePath(sub("^--file=", "", flag)))) + if (file.exists(file.path(root, "DESCRIPTION"))) setwd(root) + } +}) +suppressPackageStartupMessages(library(shigella)) + +shigella:::run_phase0_diagnostic( + n = 5, + iter_warmup = 500, + iter_sampling = 500, + tag = "n5" +) diff --git a/scripts/phase0_interactive_reproducibility_n48.R b/scripts/phase0_interactive_reproducibility_n48.R new file mode 100644 index 00000000..50c89b62 --- /dev/null +++ b/scripts/phase0_interactive_reproducibility_n48.R @@ -0,0 +1,28 @@ +# ============================================================================ +# phase0_interactive_reproducibility_n48.R — n = 48 full cohort +# +# Execution: +# salloc --time=08:00:00 --cpus-per-task=2 --mem=20G +# Rscript scripts/phase0_interactive_reproducibility_n48.R +# exit +# +# Purpose: reproduce a Phase 1 (sbatch) fit under interactive SLURM to confirm +# determinism across allocation modes. n=48 version (full cohort size). +# Run phase0_interactive_reproducibility.R (n=5) first to confirm the +# pipeline works before committing to this longer run. +# ============================================================================ +local({ + flag <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE) + if (length(flag)) { + root <- dirname(dirname(normalizePath(sub("^--file=", "", flag)))) + if (file.exists(file.path(root, "DESCRIPTION"))) setwd(root) + } +}) +suppressPackageStartupMessages(library(shigella)) + +shigella:::run_phase0_diagnostic( + n = 48, + iter_warmup = 1000, + iter_sampling = 1000, + tag = "n48" +) diff --git a/scripts/phase1_single_diagnostic.R b/scripts/phase1_single_diagnostic.R new file mode 100644 index 00000000..cc2cc944 --- /dev/null +++ b/scripts/phase1_single_diagnostic.R @@ -0,0 +1,16 @@ +# ========================================================================== +# phase1_single_diagnostic.R — n = 5 pilot +# +# Launched by slurm/phase1_single.sbatch +# +# Purpose: fit model_2 inside a SLURM sbatch job and compare with the +# Phase 0 interactive baseline to isolate SLURM-vs-code attribution. +# ========================================================================== +suppressPackageStartupMessages(library(shigella)) + +shigella:::run_phase1_diagnostic( + n = 5, + iter_warmup = 500, + iter_sampling = 500, + tag = "n5" +) diff --git a/scripts/phase1_single_diagnostic_n48.R b/scripts/phase1_single_diagnostic_n48.R new file mode 100644 index 00000000..d1d303f4 --- /dev/null +++ b/scripts/phase1_single_diagnostic_n48.R @@ -0,0 +1,17 @@ +# ========================================================================== +# phase1_single_diagnostic_n48.R — n = 48 full cohort +# +# Launched by slurm/phase1_single_n48.sbatch +# +# Purpose: fit model_2 inside a SLURM sbatch job and compare with the +# Phase 0 interactive baseline to isolate SLURM-vs-code attribution. +# n=48 version (full cohort size). +# ========================================================================== +suppressPackageStartupMessages(library(shigella)) + +shigella:::run_phase1_diagnostic( + n = 48, + iter_warmup = 1000, + iter_sampling = 1000, + tag = "n48" +) diff --git a/slurm/phase1_single.sbatch b/slurm/phase1_single.sbatch new file mode 100644 index 00000000..9e8f52f7 --- /dev/null +++ b/slurm/phase1_single.sbatch @@ -0,0 +1,78 @@ +#!/bin/bash +# ========================================================================== +# phase1_single.sbatch — SLURM single-task diagnostic run +# +# Purpose: Run the SAME n=5 fit as phase0_interactive_reproducibility.R, but +# inside a SLURM compute job. Direct comparison isolates "is the issue +# SLURM or the code itself?" +# +# ========================================================================== + +set -euo pipefail + +# ========================================================================== +# 1. Environment activation +# ========================================================================== +source "$HOME/miniconda3/etc/profile.d/conda.sh" +conda activate r_chapter2 + +# ========================================================================== +# 2. Banner +# ========================================================================== +echo "========================================================================" +echo " PHASE 1 SLURM SINGLE JOB" +echo "========================================================================" +echo " Job ID: $SLURM_JOB_ID" +echo " Job name: $SLURM_JOB_NAME" +echo " Submit dir: $SLURM_SUBMIT_DIR" +echo " Compute node: $(hostname)" +echo " Started: $(date)" +echo " R binary: $(which R)" +echo " Rscript: $(which Rscript)" +echo " CPUs/task: $SLURM_CPUS_PER_TASK" +echo " Mem/node: ${SLURM_MEM_PER_NODE:-unset}" +echo "========================================================================" + +# ========================================================================== +# 3. Per-task compile dir — avoids /home noexec + array collision risk +# ========================================================================== +export STAN_COMPILE_DIR="/tmp/$USER/cmdstan_bin_phase1_${SLURM_JOB_ID}" +mkdir -p "$STAN_COMPILE_DIR" +echo "=== STAN_COMPILE_DIR=$STAN_COMPILE_DIR ===" + +# ========================================================================== +# 4. Disable thread oversubscription +# ========================================================================== +export OMP_NUM_THREADS=1 +export MKL_NUM_THREADS=1 +export OPENBLAS_NUM_THREADS=1 + +# ========================================================================== +# 5. Change to submit dir and run +# ========================================================================== +cd "$SLURM_SUBMIT_DIR" + +mkdir -p logs/phase1 outputs/phase1 + +Rscript scripts/phase1_single_diagnostic.R + +EXIT_CODE=$? + +# ========================================================================== +# 6. Footer with exit code +# ========================================================================== +echo "" +echo "========================================================================" +echo " PHASE 1 FINISHED" +echo "========================================================================" +echo " Exit code: $EXIT_CODE" +echo " Finished: $(date)" +echo " STAN_COMPILE_DIR final size: $(du -sh $STAN_COMPILE_DIR 2>/dev/null || echo 'gone')" +echo "========================================================================" + +# Cleanup compile dir only if success — keep on failure for inspection +if [ "$EXIT_CODE" -eq 0 ]; then + rm -rf "$STAN_COMPILE_DIR" +fi + +exit $EXIT_CODE diff --git a/slurm/phase1_single_n48.sbatch b/slurm/phase1_single_n48.sbatch new file mode 100644 index 00000000..7387bdb4 --- /dev/null +++ b/slurm/phase1_single_n48.sbatch @@ -0,0 +1,78 @@ +#!/bin/bash +# ========================================================================== +# phase1_single.sbatch — SLURM single-task diagnostic run +# +# Purpose: Run the SAME n=48 fit as phase0_interactive_reproducibility_n48.R, but +# inside a SLURM compute job. Direct comparison isolates "is the issue +# SLURM or the code itself?" +# +# ========================================================================== + +set -euo pipefail + +# ========================================================================== +# 1. Environment activation +# ========================================================================== +source "$HOME/miniconda3/etc/profile.d/conda.sh" +conda activate r_chapter2 + +# ========================================================================== +# 2. Banner +# ========================================================================== +echo "========================================================================" +echo " PHASE 1 SLURM SINGLE JOB" +echo "========================================================================" +echo " Job ID: $SLURM_JOB_ID" +echo " Job name: $SLURM_JOB_NAME" +echo " Submit dir: $SLURM_SUBMIT_DIR" +echo " Compute node: $(hostname)" +echo " Started: $(date)" +echo " R binary: $(which R)" +echo " Rscript: $(which Rscript)" +echo " CPUs/task: $SLURM_CPUS_PER_TASK" +echo " Mem/node: ${SLURM_MEM_PER_NODE:-unset}" +echo "========================================================================" + +# ========================================================================== +# 3. Per-task compile dir — avoids /home noexec + array collision risk +# ========================================================================== +export STAN_COMPILE_DIR="/tmp/$USER/cmdstan_bin_phase1_n48_${SLURM_JOB_ID}" +mkdir -p "$STAN_COMPILE_DIR" +echo "=== STAN_COMPILE_DIR=$STAN_COMPILE_DIR ===" + +# ========================================================================== +# 4. Disable thread oversubscription +# ========================================================================== +export OMP_NUM_THREADS=1 +export MKL_NUM_THREADS=1 +export OPENBLAS_NUM_THREADS=1 + +# ========================================================================== +# 5. Change to submit dir and run +# ========================================================================== +cd "$SLURM_SUBMIT_DIR" + +mkdir -p logs/phase1 outputs/phase1 + +Rscript scripts/phase1_single_diagnostic_n48.R + +EXIT_CODE=$? + +# ========================================================================== +# 6. Footer with exit code +# ========================================================================== +echo "" +echo "========================================================================" +echo " PHASE 1 FINISHED" +echo "========================================================================" +echo " Exit code: $EXIT_CODE" +echo " Finished: $(date)" +echo " STAN_COMPILE_DIR final size: $(du -sh $STAN_COMPILE_DIR 2>/dev/null || echo 'gone')" +echo "========================================================================" + +# Cleanup compile dir only if success — keep on failure for inspection +if [ "$EXIT_CODE" -eq 0 ]; then + rm -rf "$STAN_COMPILE_DIR" +fi + +exit $EXIT_CODE diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 00000000..80e674ac --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(shigella) + +test_check("shigella") diff --git a/tests/testthat/test-compute_kinetics_at_time.R b/tests/testthat/test-compute_kinetics_at_time.R new file mode 100644 index 00000000..aa37b48f --- /dev/null +++ b/tests/testthat/test-compute_kinetics_at_time.R @@ -0,0 +1,18 @@ +test_that("compute_kinetics_at_time aborts on degenerate shape == 1", { + expect_error( + shigella:::.compute_kinetics_at_time( + log_y0 = log(1), log_y1m0 = log(9), log_t1 = log(2), + log_alpha = log(0.5), log_rm1 = -Inf, tt = 5 + ), + regexp = "degenerate" + ) +}) + +test_that("compute_kinetics_at_time succeeds in growth phase even when log_rm1 = -Inf", { + expect_no_error( + shigella:::.compute_kinetics_at_time( + log_y0 = log(1), log_y1m0 = log(9), log_t1 = log(2), + log_alpha = log(0.5), log_rm1 = -Inf, tt = 1 # tt < t1_j = 2 + ) + ) +}) diff --git a/tests/testthat/test-postprocess_stan_output.R b/tests/testthat/test-postprocess_stan_output.R new file mode 100644 index 00000000..315e9e1e --- /dev/null +++ b/tests/testthat/test-postprocess_stan_output.R @@ -0,0 +1,215 @@ +test_that("postprocess_stan_output is callable", { + expect_true(is.function(postprocess_stan_output)) +}) + +test_that("summarize_matrix_array returns a list of matrices", { + # Build a minimal mock draws_array for array[2] corr_matrix[3] Omega_P + # Variable names: Omega_P[k,i,j] for k in 1:2, i in 1:3, j in 1:3 + var_names <- c() + for (k in 1:2) { + for (i in 1:3) { + for (j in 1:3) { + var_names <- c(var_names, sprintf("Omega_P[%d,%d,%d]", k, i, j)) + } + } + } + set.seed(1) + arr_data <- array( + runif(1 * 1 * length(var_names)), + dim = c(1L, 1L, length(var_names)), + dimnames = list( + iteration = "1", + chain = "1", + variable = var_names + ) + ) + result <- shigella:::.summarize_matrix_array(arr_data, "Omega_P", + n_arr = 2L, n_row = 3L, n_col = 3L) + + expect_type(result, "list") + expect_length(result, 2L) + expect_equal(dim(result[[1]]), c(3L, 3L)) + expect_equal(dim(result[[2]]), c(3L, 3L)) + # Each cell should equal the draw value (only one draw, so median = the value) + expect_equal(result[[1]][1, 1], + arr_data[1, 1, "Omega_P[1,1,1]"], + tolerance = 1e-10) +}) + +# ─── Mock Stan fit for fast (non-Stan) pipeline tests ─────────────────────── + +# Builds a minimal mock CmdStanMCMC-like list with just enough interface to +# satisfy postprocess_stan_output: $draws(variables) and $metadata(). +# Avoids running Stan in routine devtools::test() calls. +.make_mock_stan_fit <- function(N = 2L, K = 2L, + param_names = c("y0", "y1", "t1", + "alpha", "shape"), + n_iter = 4L, n_chains = 2L) { + P <- length(param_names) + + # Parameter variables: pname[subj, antigen] for all combinations + param_vars <- unlist(lapply(param_names, function(p) { + sprintf("%s[%d,%d]", p, + rep(seq_len(N), each = K), + rep(seq_len(K), N)) + })) + + # model_1 Omega_P: array[K] corr_matrix[P] -> Omega_P[k,i,j] + omega_P_vars <- unlist(lapply(seq_len(K), function(k) { + as.vector(outer(seq_len(P), seq_len(P), + function(i, j) sprintf("Omega_P[%d,%d,%d]", k, i, j))) + })) + + log_lik_vars <- sprintf("log_lik[%d]", seq_len(N)) + + all_vars <- c(param_vars, omega_P_vars, log_lik_vars) + + set.seed(42) + raw_arr <- array( + abs(rnorm(n_iter * n_chains * length(all_vars), mean = 1, sd = 0.2)), + dim = c(n_iter, n_chains, length(all_vars)) + ) + dimnames(raw_arr) <- list(NULL, NULL, all_vars) + draws_full <- posterior::as_draws_array(raw_arr) + + list( + draws = function(variables = NULL, ...) { + if (is.null(variables)) return(draws_full) + all_var_names <- posterior::variables(draws_full) + matched <- unlist(lapply(variables, function(v) { + grep(paste0("^", v, "(\\[|$)"), all_var_names, value = TRUE) + })) + if (length(matched) == 0) { + stop("No variables matched: ", paste(variables, collapse = ", ")) + } + posterior::subset_draws(draws_full, variable = matched) + }, + metadata = function() { + list(stan_variables = c(param_names, "Omega_P", "log_lik")) + } + ) +} + +test_that("postprocess_stan_output processes mock draws without Stan (model_1)", { + skip_if_not_installed("posterior") + + ids <- c("s1", "s2") + antigens <- c("IgG", "IgA") + param_names <- c("y0", "y1", "t1", "alpha", "shape") + n_iter <- 4L + n_chains <- 2L + + mock_fit <- .make_mock_stan_fit( + N = length(ids), + K = length(antigens), + param_names = param_names, + n_iter = n_iter, + n_chains = n_chains + ) + + result <- postprocess_stan_output( + stan_fit = mock_fit, + ids = ids, + antigens = antigens, + model = "model_1", + stratification = "test_stratum" + ) + + # Returns named list with sr_tibble and cov_summaries + expect_type(result, "list") + expect_named(result, c("sr_tibble", "cov_summaries")) + + # sr_tibble has expected columns + expect_s3_class(result$sr_tibble, "tbl_df") + expected_cols <- c("Iteration", "Chain", "Parameter", "Iso_type", + "Stratification", "Subject", "value") + expect_true(all(expected_cols %in% names(result$sr_tibble))) + + # All 5 parameters present for both subjects and antigens + expect_equal(length(unique(result$sr_tibble$Parameter)), length(param_names)) + expect_equal(sort(unique(result$sr_tibble$Subject)), sort(ids)) + expect_equal(sort(unique(result$sr_tibble$Iso_type)), sort(antigens)) + + # Row count: n_params × N × K × (n_iter × n_chains) + expect_equal( + nrow(result$sr_tibble), + length(param_names) * length(ids) * length(antigens) * n_iter * n_chains + ) + + # Numerical sanity: all values finite, ESS-free check + expect_true(all(is.finite(result$sr_tibble$value))) + + # Omega_P extracted (model_1 produces a named list of K matrices) + expect_true("Omega_P" %in% names(result$cov_summaries)) + omega_P <- result$cov_summaries$Omega_P + expect_type(omega_P, "list") + expect_equal(length(omega_P), length(antigens)) + for (mat in omega_P) { + expect_equal(dim(mat), c(length(param_names), length(param_names))) + expect_true(all(is.finite(mat))) + } +}) + +test_that("postprocess_stan_output produces sr_model output — model_2 (slow)", { + skip_if( + Sys.getenv("RUN_STAN_TESTS") != "true", + "Stan tests are skipped unless RUN_STAN_TESTS=true." + ) + skip_if_not_installed("cmdstanr") + + sim <- sim_correlated_case_data(n = 3, seed = 2026) + + fit <- run_mod_stan( + data = sim, + model = "model_2", + chains = 1, + iter_warmup = 100, + iter_sampling = 100, + refresh = 0, + show_messages = FALSE + ) + + expect_s3_class(fit, "sr_model") + expect_true("priors" %in% names(attributes(fit))) + expect_true("fitted_residuals" %in% names(attributes(fit))) +}) + +test_that("postprocess_stan_output extracts per-biomarker Omega_P — model_1 (slow)", { + skip_if( + Sys.getenv("RUN_STAN_TESTS") != "true", + "Stan tests are skipped unless RUN_STAN_TESTS=true." + ) + skip_if_not_installed("cmdstanr") + + sim <- sim_correlated_case_data(n = 3, seed = 2026) + + fit <- run_mod_stan( + data = sim, + model = "model_1", + chains = 1, + iter_warmup = 100, + iter_sampling = 100, + refresh = 0, + show_messages = FALSE + ) + + expect_s3_class(fit, "sr_model") + + # run_mod_stan flattens cov_summaries into individual top-level attributes + expect_true("Omega_P" %in% names(attributes(fit))) + omega_P <- attr(fit, "Omega_P") + + # model_1 Omega_P should be a named list (one 5x5 matrix per biomarker) + expect_type(omega_P, "list") + expect_equal(length(omega_P), 2L) # K=2 biomarkers (biomarker_1, biomarker_2) + expect_equal(names(omega_P), c("biomarker_1", "biomarker_2")) + for (mat in omega_P) { + expect_equal(dim(mat), c(5L, 5L)) + expect_equal(rownames(mat), c("y0", "y1", "t1", "alpha", "shape")) + expect_equal(colnames(mat), c("y0", "y1", "t1", "alpha", "shape")) + } + + # model_1 should NOT have Kronecker-only summaries + expect_false("Omega_B" %in% names(attributes(fit))) + expect_false("Omega_eps" %in% names(attributes(fit))) +}) diff --git a/tests/testthat/test-prep_data_stan.R b/tests/testthat/test-prep_data_stan.R new file mode 100644 index 00000000..04254e70 --- /dev/null +++ b/tests/testthat/test-prep_data_stan.R @@ -0,0 +1,29 @@ +test_that("prep_data_stan accepts a case_data object directly", { + sim <- sim_correlated_case_data(n = 5, seed = 2026) + + stan_data <- prep_data_stan(sim) + + expect_type(stan_data, "list") + expect_true(all(c("N", "K", "P", "max_obs", "n_obs", + "time_obs", "log_y") %in% names(stan_data))) + expect_equal(stan_data$N, 5) # newperson dropped via add_newperson=FALSE +}) + +test_that("prep_data_stan accepts a prepped_jags_data object", { + sim <- sim_correlated_case_data(n = 5, seed = 2026) + prepped <- serodynamics::prep_data(sim, add_newperson = FALSE) + + stan_data <- prep_data_stan(prepped) + + expect_type(stan_data, "list") + expect_equal(stan_data$N, 5) +}) + +test_that("prep_data_stan errors informatively on unsupported input class", { + bad_input <- list(not_a_real_jags_data = TRUE) + class(bad_input) <- "unrecognized_class" + expect_error( + prep_data_stan(bad_input), + regexp = "case_data" # message names the expected classes + ) +}) diff --git a/tests/testthat/test-prep_priors_stan.R b/tests/testthat/test-prep_priors_stan.R new file mode 100644 index 00000000..69b07d06 --- /dev/null +++ b/tests/testthat/test-prep_priors_stan.R @@ -0,0 +1,37 @@ +test_that("prep_priors_stan returns a list with expected names", { + priors <- prep_priors_stan(model = "model_2") + + expect_type(priors, "list") + expect_true(length(priors) > 0) +}) + +test_that("prep_priors_stan defaults are weakly informative", { + priors <- prep_priors_stan(model = "model_2") + + ## Prior SDs should not be absurdly diffuse (e.g., 316 from old JAGS + ## translation) — that caused Stan HMC to wander during warmup. + if (!is.null(priors$mu_hyp_sd)) { + expect_true(all(priors$mu_hyp_sd > 0)) + expect_true(all(priors$mu_hyp_sd <= 20)) + } +}) + +test_that("prep_priors_stan includes biomarker priors only for model_2", { + priors_model_2 <- prep_priors_stan(model = "model_2") + priors_model_1 <- prep_priors_stan(model = "model_1") + + expect_true("tau_B_scale" %in% names(priors_model_2)) + expect_true("lkj_B_eta" %in% names(priors_model_2)) + expect_false("tau_B_scale" %in% names(priors_model_1)) + expect_false("lkj_B_eta" %in% names(priors_model_1)) +}) + +test_that("prep_priors_stan model_2 structure is stable", { + priors <- prep_priors_stan(model = "model_2") + expect_equal( + names(priors), + c("mu_hyp_mean", "mu_hyp_sd", "tau_P_scale", "tau_eps_scale", + "lkj_P_eta", "tau_B_scale", "lkj_B_eta", "lkj_eps_eta") + ) + expect_true(all(vapply(priors, is.numeric, logical(1L)))) +}) diff --git a/tests/testthat/test-run_mod_stan.R b/tests/testthat/test-run_mod_stan.R new file mode 100644 index 00000000..e0e7b630 --- /dev/null +++ b/tests/testthat/test-run_mod_stan.R @@ -0,0 +1,40 @@ +test_that("run_mod_stan has required arguments", { + fn_args <- names(formals(run_mod_stan)) + expect_true("data" %in% fn_args) + expect_true("model" %in% fn_args) + expect_true("seed" %in% fn_args) +}) + +test_that("run_mod_stan completes a minimal fit (slow)", { + skip_if( + Sys.getenv("RUN_STAN_TESTS") != "true", + "Stan tests are skipped unless RUN_STAN_TESTS=true." + ) + skip_if_not_installed("cmdstanr") + + sim <- sim_correlated_case_data(n = 3, seed = 2026) + + warnings_seen <- character() + fit <- withCallingHandlers( + run_mod_stan( + data = sim, + model = "model_2", + chains = 1, + iter_warmup = 200, + iter_sampling = 100, + adapt_delta = 0.99, + max_treedepth = 15, + refresh = 0, + show_messages = FALSE + ), + warning = function(w) { + warnings_seen <<- c(warnings_seen, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + + # Low-iteration smoke fit may emit convergence warnings — that is expected. + # Assert that the function ran and returned a valid object; do not assert + # on the absence of warnings since they are diagnostic signals, not errors. + expect_s3_class(fit, "sr_model") +}) diff --git a/tests/testthat/test-run_phase0_diagnostic.R b/tests/testthat/test-run_phase0_diagnostic.R new file mode 100644 index 00000000..7541d4ef --- /dev/null +++ b/tests/testthat/test-run_phase0_diagnostic.R @@ -0,0 +1,30 @@ +test_that("run_phase0_diagnostic has required arguments", { + fn_args <- names(formals(shigella:::run_phase0_diagnostic)) + expect_true("n" %in% fn_args) + expect_true("tag" %in% fn_args) + expect_true("compile_dir" %in% fn_args) +}) + +test_that("run_phase0_diagnostic runs with minimal Stan settings", { + skip_if_not(Sys.getenv("RUN_STAN_TESTS") == "true", + "Skipping Stan test (set RUN_STAN_TESTS=true to enable)") + + out_dir <- file.path(tempdir(), "phase0_test") + on.exit(unlink(out_dir, recursive = TRUE)) + + result <- run_phase0_diagnostic( + n = 3, + iter_warmup = 50, + iter_sampling = 50, + tag = "test", + output_dir = out_dir, + chains = 1L + ) + + expect_true(is.list(result) || is.null(result)) + if (!is.null(result)) { + expect_equal(result$status, "OK") + expect_equal(result$n_subjects, 3L) + expect_true(file.exists(file.path(out_dir, "one_fit_test.rds"))) + } +}) diff --git a/tests/testthat/test-run_phase1_diagnostic.R b/tests/testthat/test-run_phase1_diagnostic.R new file mode 100644 index 00000000..7e9cfb81 --- /dev/null +++ b/tests/testthat/test-run_phase1_diagnostic.R @@ -0,0 +1,35 @@ +test_that("run_phase1_diagnostic has required arguments", { + fn_args <- names(formals(shigella:::run_phase1_diagnostic)) + expect_true("n" %in% fn_args) + expect_true("tag" %in% fn_args) + expect_true("compile_dir" %in% fn_args) + expect_true("phase0_dir" %in% fn_args) +}) + +test_that("run_phase1_diagnostic runs with minimal Stan settings", { + skip_if_not(Sys.getenv("RUN_STAN_TESTS") == "true", + "Skipping Stan test (set RUN_STAN_TESTS=true to enable)") + + out_dir <- file.path(tempdir(), "phase1_test") + ph0_dir <- file.path(tempdir(), "phase0_test") + on.exit(unlink(c(out_dir, ph0_dir), recursive = TRUE)) + + result <- run_phase1_diagnostic( + n = 3, + iter_warmup = 50, + iter_sampling = 50, + tag = "test", + output_dir = out_dir, + phase0_dir = ph0_dir, + chains = 1L + ) + + expect_true(is.list(result) || is.null(result)) + if (!is.null(result)) { + expect_equal(result$status, "OK") + expect_equal(result$n_subjects, 3L) + # run_phase1_diagnostic saves as one_fit_{tag}_jobid_{timestamp}.rds + rds_files <- list.files(out_dir, pattern = "^one_fit_test.*\\.rds$") + expect_true(length(rds_files) > 0) + } +}) diff --git a/tests/testthat/test-sim_correlated_case_data.R b/tests/testthat/test-sim_correlated_case_data.R new file mode 100644 index 00000000..4f4bb7ae --- /dev/null +++ b/tests/testthat/test-sim_correlated_case_data.R @@ -0,0 +1,82 @@ +test_that("sim_correlated_case_data returns a case_data object", { + sim <- sim_correlated_case_data(n = 5, seed = 2026) + + expect_s3_class(sim, "case_data") + expect_true(nrow(sim) > 0) + expect_true(all(c("id", "timeindays", "antigen_iso", "value") + %in% names(sim))) +}) + +test_that("sim_correlated_case_data attaches the truth attributes", { + omega_B <- matrix(c(1, 0.5, 0.5, 1), nrow = 2) + sim <- sim_correlated_case_data(n = 5, omega_B = omega_B, seed = 2026) + + truth <- attr(sim, "truth") + expect_type(truth, "list") + expect_equal(truth$omega_B, omega_B) + expect_true(!is.null(attr(sim, "theta_true"))) +}) + +test_that("sim_correlated_case_data supports n = 1", { + sim <- sim_correlated_case_data(n = 1, seed = 2026) + + expect_s3_class(sim, "case_data") + expect_equal(length(unique(sim$id)), 1) +}) + +test_that(".validate_corr_matrix rejects non-square omega", { + # 2x3 matrix fails the square check before .validate_corr_matrix even runs + # (caught by .validate_sim_inputs dimension check) + bad_omega <- matrix(0, nrow = 2, ncol = 3) + expect_error( + sim_correlated_case_data(n = 2, omega_B = bad_omega, seed = 1), + regexp = "omega_B" + ) +}) + +test_that(".validate_corr_matrix rejects non-symmetric omega", { + bad_omega <- matrix(c(1, 0.5, 0.3, 1), nrow = 2) # asymmetric + expect_error( + sim_correlated_case_data(n = 2, omega_B = bad_omega, seed = 1), + regexp = "symmetric" + ) +}) + +test_that(".validate_corr_matrix rejects non-unit-diagonal omega", { + bad_omega <- matrix(c(2, 0, 0, 2), nrow = 2) # diagonal != 1 + expect_error( + sim_correlated_case_data(n = 2, omega_B = bad_omega, seed = 1), + regexp = "unit diagonal" + ) +}) + +test_that(".validate_corr_matrix rejects non-PSD omega", { + # Symmetric, unit diagonal, but min eigenvalue = -1 + bad_omega <- matrix(c(1, 2, 2, 1), nrow = 2) + expect_error( + sim_correlated_case_data(n = 2, omega_B = bad_omega, seed = 1), + regexp = "positive semi-definite" + ) +}) + +test_that("sim_correlated_case_data theta_true structure is stable", { + sim <- sim_correlated_case_data(n = 5, seed = 2026) + theta <- attr(sim, "theta_true") + + expect_equal(dim(theta), c(5L, 5L, 2L)) + expect_equal( + dimnames(theta), + list( + subject = as.character(1:5), + param = c("log_y0", "log_y1m0", "log_t1", "log_alpha", "log_rm1"), + biomarker = c("biomarker_1", "biomarker_2") + ) + ) + + truth <- attr(sim, "truth") + expect_equal( + names(truth), + c("mu", "tau_P", "tau_B", "tau_eps", "omega_P", "omega_B", "omega_eps", + "sigma_P", "sigma_B", "sigma_eps") + ) +}) diff --git a/tests/testthat/test-write_status.R b/tests/testthat/test-write_status.R new file mode 100644 index 00000000..ff1d1ad1 --- /dev/null +++ b/tests/testthat/test-write_status.R @@ -0,0 +1,37 @@ +test_that("write_status appends a line to the log file", { + log_file <- tempfile(fileext = ".txt") + on.exit(unlink(log_file)) + + write_status(log_file, "TEST", "hello") + + lines <- readLines(log_file) + expect_length(lines, 1L) + expect_match(lines[1], "STEP=TEST") + expect_match(lines[1], "hello") +}) + +test_that("write_status defaults msg to empty string", { + log_file <- tempfile(fileext = ".txt") + on.exit(unlink(log_file)) + + write_status(log_file, "STEP_A") + + lines <- readLines(log_file) + expect_length(lines, 1L) + expect_match(lines[1], "STEP=STEP_A") +}) + +test_that("write_status appends multiple calls in order", { + log_file <- tempfile(fileext = ".txt") + on.exit(unlink(log_file)) + + write_status(log_file, "FIRST") + write_status(log_file, "SECOND") + write_status(log_file, "THIRD") + + lines <- readLines(log_file) + expect_length(lines, 3L) + expect_match(lines[1], "FIRST") + expect_match(lines[2], "SECOND") + expect_match(lines[3], "THIRD") +}) diff --git a/tmp_parse.py b/tmp_parse.py new file mode 100644 index 00000000..26ca905e --- /dev/null +++ b/tmp_parse.py @@ -0,0 +1,50 @@ +import json +from collections import Counter + +path = '/home/runner/.claude/projects/-home-runner-work-shigella-shigella/2caf7287-97e8-4e7d-ac4a-b168a221f2d7.jsonl' +bash_cmds = Counter() +raw_cmds = [] +mcp_tools = Counter() + +with open(path) as f: + for line in f: + try: + obj = json.loads(line) + except Exception: + continue + msg = obj.get('message', {}) + if msg.get('role') != 'assistant': + continue + for block in msg.get('content', []): + if not isinstance(block, dict) or block.get('type') != 'tool_use': + continue + name = block.get('name', '') + inp = block.get('input', {}) + if name == 'Bash': + cmd = inp.get('command', '').strip() + raw_cmds.append(cmd) + tokens = cmd.split() + if tokens: + lead = tokens[0] + while '=' in lead and len(tokens) > 1: + tokens = tokens[1:] + lead = tokens[0] + sub = tokens[1] if len(tokens) > 1 else '' + key = (lead + ' ' + sub).strip() + bash_cmds[key] += 1 + elif name.startswith('mcp__'): + mcp_tools[name] += 1 + +print('=== BASH (lead+sub) ===') +for k, v in bash_cmds.most_common(40): + print(f'{v:3d} {k}') + +print() +print('=== RAW COMMANDS ===') +for cmd in raw_cmds: + print(repr(cmd[:150])) + +print() +print('=== MCP ===') +for k, v in mcp_tools.most_common(20): + print(f'{v:3d} {k}') diff --git a/vignettes/articles/_metadata.yml b/vignettes/articles/_metadata.yml new file mode 100644 index 00000000..bc68c864 --- /dev/null +++ b/vignettes/articles/_metadata.yml @@ -0,0 +1,20 @@ +format: + html: + toc: true + toc-depth: 3 + number-sections: true + embed-resources: true + theme: cosmo + code-fold: true + fig-width: 7 + fig-height: 4.5 + docx: + toc: true + number-sections: true +execute: + echo: false + warning: false + message: false +editor: + markdown: + wrap: 72 diff --git a/vignettes/articles/chapter2.qmd b/vignettes/articles/chapter2.qmd new file mode 100644 index 00000000..aedd4597 --- /dev/null +++ b/vignettes/articles/chapter2.qmd @@ -0,0 +1,480 @@ +--- +title: "Chapter 2 — Correlated Multivariate Antibody Kinetics" +subtitle: "A Kronecker-structured Bayesian hierarchical model: methodology and Phase 0/1 simulation diagnostics" +author: "Kwan Ho Lee" +date: today +--- + +```{r setup} +#| include: false +options(knitr.kable.NA = "") # nolint: undesirable_function_linter +``` + +# Introduction {#sec-intro} + +## Motivation + +Chapter 1 established a univariate Bayesian hierarchical model for the +post-symptom-onset antibody trajectory of a single antigen-isotype +biomarker, fitting each biomarker independently. +This independence +assumption is biologically unrealistic: within the same subject, IgG and +IgA responses share the same infection event, the same immune system, +and the same antigenic stimulation, so their kinetic parameters are +expected to co-vary. +Ignoring this dependence discards information that +could improve parameter estimation precision and downstream +seroincidence inference. + +## Aim of Chapter 2 + +Chapter 2 extends the Chapter 1 framework to a multivariate hierarchical +model in which the random-effect covariance of the per-subject log-scale +kinetic parameters has a Kronecker product structure, +$\Sigma = \Sigma_B \otimes \Sigma_P$. Here $\Sigma_B$ captures +correlation between biomarkers and $\Sigma_P$ captures correlation +between kinetic parameters, reducing the unconstrained $5K \times 5K$ +joint covariance to a structured, interpretable factorization. + +## Scope of this document + +@sec-methods derives the model and maps each mathematical object to its +implementation in the `shigella` R package and the Stan model file. +@sec-results presents the Phase 0 (interactive SLURM) and Phase 1 (batch +SLURM) reproducibility check at two sample sizes relevant to the Chapter +1 cohort. @sec-discussion interprets the n = 48 fit pathology and +proposes next steps. + +# Methods {#sec-methods} + +## The Chapter 1 univariate kinetic model + +Following Teunis et al. (2016), the post-symptom-onset antibody +concentration $y(t)$ in a subject's blood at time $t$ is modeled as a +deterministic two-phase trajectory. +This formulation builds on the +within-host rise--decay seroresponse models developed by de Graaf et al. +(2014) and Teunis et al. (2016): + +$$ +y(t) \;=\; +\begin{cases} +y_0 \exp(\beta t), & 0 \le t \le t_1 \quad \text{(rise phase)} \\[6pt] +\left[ y_1^{1-r} - (1-r)\,\alpha\,(t - t_1) \right]^{1/(1-r)}, & t > t_1 \quad \text{(decay phase)} +\end{cases} +$$ + +where $\beta = \log(y_1 / y_0) / t_1$. The five per-subject kinetic +parameters are $(y_0, y_1, t_1, \alpha, r)$, denoting baseline level, +peak level, time-to-peak, decay rate, and decay shape. + +For hierarchical modeling, these are reparameterized into unconstrained +log-scale coordinates: + +$$ +\boldsymbol{\theta}_{i,k} +\;=\; +\bigl( \log y_{0}, \;\; \log(y_1 - y_0), \;\; \log t_1, \;\; \log\alpha, \;\; \log(r - 1) \bigr)_{i,k} +\;\in\; \mathbb{R}^{5}, +$$ + +for subject $i$ and biomarker $k$. + +## The multivariate extension {#sec-mv} + +Let $K$ index biomarkers (antigen-isotype pairs; $K=2$ for the IgG/IgA +pair) and $P = 5$ the number of kinetic parameters. +Stack the per-subject log-scale parameters across biomarkers: + +$$ +\boldsymbol{\theta}_i \;=\; +\bigl( \boldsymbol{\theta}_{i,1}^{\top},\,\boldsymbol{\theta}_{i,2}^{\top},\,\ldots,\,\boldsymbol{\theta}_{i,K}^{\top} \bigr)^{\top} +\;\in\; \mathbb{R}^{KP}. +$$ + +The Chapter 2 random-effect prior is + +$$ +\boldsymbol{\theta}_i \;\sim\; \mathcal{N}_{KP}\bigl( \boldsymbol{\mu},\; \Sigma_B \otimes \Sigma_P \bigr), +\qquad i = 1, \ldots, n, +$$ {#eq-kron} + +where $\otimes$ denotes the Kronecker product, +$\Sigma_B \in \mathbb{R}^{K \times K}$ is the between-biomarker +covariance, and $\Sigma_P \in \mathbb{R}^{P \times P}$ is the +between-parameter covariance. +This is the multivariate hierarchical +parameterization used to separate covariance across response dimensions +from covariance across kinetic parameters (Gelman et al., 2014). + +The Kronecker structure factorizes the cross-(biomarker, parameter) +covariance as + +$$ +\mathrm{Cov}(\theta_{i, k, p},\, \theta_{i, k', p'}) +\;=\; \Sigma_B[k, k'] \cdot \Sigma_P[p, p']. +$$ + +## Identifiability and the LKJ–scale parameterization + +The decomposition $\Sigma_B \otimes \Sigma_P$ is non-unique up to a +scalar: for any $c \neq 0$, the pair $(c\Sigma_B, c^{-1}\Sigma_P)$ +yields the same product. +To resolve this, $\Sigma_B$ and $\Sigma_P$ are +constrained to be correlation matrices and the marginal variances +are absorbed into separate positive scale vectors: + +$$ +\Sigma_B \;=\; \mathrm{diag}(\boldsymbol{\tau}_B)\,\Omega_B\,\mathrm{diag}(\boldsymbol{\tau}_B), +\quad +\Sigma_P \;=\; \mathrm{diag}(\boldsymbol{\tau}_P)\,\Omega_P\,\mathrm{diag}(\boldsymbol{\tau}_P), +$$ + +with $\Omega_B, \Omega_P$ correlation matrices (unit diagonal) and +$\boldsymbol{\tau}_B \in \mathbb{R}_+^K$, +$\boldsymbol{\tau}_P \in \mathbb{R}_+^P$. + +## Priors + +| Parameter | Prior | Rationale | +|---------------------------|------------------|---------------------------| +| $\boldsymbol{\mu} \in \mathbb{R}^{KP}$ | $\mathcal{N}(\mathbf{0}, 5\,\mathbf{I})$ componentwise | Weakly informative, on the log-scale parameter space | +| $\Omega_B$ | $\mathrm{LKJ}(\eta = 2)$ | Mild concentration around independence; permits either sign of correlation | +| $\Omega_P$ | $\mathrm{LKJ}(\eta = 2)$ | Same rationale | +| $\boldsymbol{\tau}_B$ | $\mathrm{half\text{-}}\mathcal{N}(0,\,1)$ componentwise | Positive scale; weakly informative | +| $\boldsymbol{\tau}_P$ | $\mathrm{half\text{-}}\mathcal{N}(0,\,1)$ componentwise | Same | + +The LKJ prior is used because it defines a proper prior over correlation +matrices and allows transparent control over concentration around the +identity matrix (Lewandowski et al., 2009). +I do not impose a sign +constraint on $\rho_B$: although IgG and IgA responses share the same +infection event, class switching and mucosal/systemic response timing +can differ, so early post-infection data may plausibly support either +positive or negative within-subject IgG--IgA correlation (Stavnezer et +al., 2008; Mattoo & Cherry, 2005). + +The target estimand of interest for this calibration study is the +off-diagonal entry $\rho_B \equiv \Omega_B[1,2]$, the between-biomarker +correlation. + +## Mapping mathematics to code {#sec-code-map} + +| Mathematical object | Implemented in | +|------------------------------------|------------------------------------| +| Two-phase trajectory $y(t)$ + log-scale parameterization | `inst/stan/model_2.stan` (transformed parameters block) | +| Synthetic data drawn from @eq-kron with known $\rho_B$ | `R/sim_correlated_case_data.R` | +| Construction of $\Sigma_B \otimes \Sigma_P$ from $(\Omega_B, \Omega_P, \boldsymbol{\tau}_B, \boldsymbol{\tau}_P)$ | `R/build_sigma_matrices.R` | +| Per-subject parameter draws $\boldsymbol{\theta}_i$ | `R/draw_subject_params.R` | +| Mean trajectory $\log \mu_k(t)$ | `R/compute_log_mu_k.R` | +| Stan data list assembly | `R/prep_data_stan.R` | +| Prior hyperparameter list | `R/prep_priors_stan.R` | +| MCMC fit invocation (cmdstanr) | `R/run_mod_stan.R` | +| Posterior summary + diagnostics | `R/postprocess_stan_output.R`, `R/extract_model1_omega_p_stan.R` | +| Phase 0 (interactive SLURM) workflow | `R/run_phase0_diagnostic.R` | +| Phase 1 (batch SLURM) workflow | `R/run_phase1_diagnostic.R` | + +The Stan implementation should use a non-centered hierarchical +parameterization, because non-centered forms often reduce difficult +funnel-like posterior geometry in hierarchical models fitted with HMC +(Betancourt & Girolami, 2015). +In code, this means sampling +standard-normal latent variables and transforming them through the +Cholesky factor of the structured covariance matrix rather than sampling +subject-level parameters directly. + +## Simulation design + +The package function `sim_correlated_case_data()` generates data from +@eq-kron with a user-specified ground-truth $\rho_B$. Two sample sizes +are tested here: + +| Setting | $n$ subjects | Per-subject obs. | $\rho_B$ truth | Warmup / Sampling per chain | Chains | Seed | +|-----------|----------:|----------:|----------:|----------:|----------:|----------:| +| Light | 5 | 10 | $+0.6$ | 500 / 500 | 2 | 20260513 | +| Realistic | 48 | 10 | $+0.6$ | 1000 / 1000 | 2 | 20260513 | + +The two sample sizes are chosen as a sanity check ($n=5$) and as the +Chapter-1 IpaB analytical cohort size ($n=48$). +Each fit is run twice +with identical seed and model file but in two different execution modes +(interactive `salloc` vs batch `sbatch`) to verify execution determinism +(see @sec-results). + +For the later multi-replicate simulation study, performance summaries +should be reported with Monte Carlo standard errors, rather than +interpreted from a small pilot alone. +Following Morris et al. (2019), +the planned full simulation should summarize bias, mean squared error, +coverage, and the Monte Carlo standard error of each performance +measure. + +# Results {#sec-results} + +Sampler diagnostics are interpreted using standard HMC/NUTS guidance: +divergences indicate numerical integration failures in difficult +posterior geometry, maximum-treedepth saturation indicates that NUTS +repeatedly required the maximum allowed trajectory length, and +$\widehat R$ plus effective sample size summarize cross-chain +convergence and Monte Carlo precision (Hoffman & Gelman, 2014; +Betancourt, 2018; Vehtari et al., 2021). + +## Reproducibility at $n = 5$ {#sec-results-n5} + +@tbl-n5 reports the side-by-side comparison of Phase 0 and Phase 1 at +$n = 5$. + +```{r} +#| label: tbl-n5 +#| tbl-cap: "Phase 0 (interactive `salloc`) vs Phase 1 (batch `sbatch`) at +#| $n = 5$, $\\rho_B^{\\text{truth}} = +0.6$. Both phases use identical seed +#| and model file." + +cmp_n5 <- data.frame( + Metric = c("Status", "Elapsed (min)", "Posterior median $\\rho_B$", + "95% CrI lower (2.5%)", "95% CrI upper (97.5%)", + "ESS_bulk", "$\\hat R$", "Divergent transitions (/1000)", + "Max-treedepth hits (/1000)"), + Phase0 = c("OK", "9.54", "$-0.034$", "$-0.463$", "$+0.611$", + "6", "1.371", "13", "493"), + Phase1 = c("OK", "9.40", "$-0.034$", "$-0.463$", "$+0.611$", + "6", "1.371", "13", "493") +) + +knitr::kable(cmp_n5, + col.names = c("Metric", "Phase 0 (interactive)", + "Phase 1 (batch)"), + align = "lrr") +``` + +All numerical summaries are bit-for-bit identical between the two +execution modes. +Stan's HMC sampler is deterministic given identical +seed, model, data, and runtime, so this match confirms that: + +1. The cluster environment (SLURM batch vs interactive allocation, + compute-node `/tmp`, library paths) does not affect the fit. +2. Any pathology observed at higher $n$ is therefore a property of the + model–data combination, not an environment artifact. + +The $n = 5$ fit itself is uninformative ($\hat R = 1.37$, ESS = 6): at +this sample size the LKJ(2) prior dominates the posterior, which is +expected. The $n = 5$ run serves only as a reproducibility sanity check. + +## Reproducibility and fit pathology at $n = 48$ {#sec-results-n48} + +@tbl-n48 reports the same comparison at $n = 48$, the Chapter 1 IpaB +cohort size. + +```{r} +#| label: tbl-n48 +#| tbl-cap: "Phase 0 (interactive `salloc`) vs Phase 1 (batch `sbatch`) at +#| $n = 48$, $\\rho_B^{\\text{truth}} = +0.6$. Identical seed and model. Cells +#| flagged violate conventional diagnostic thresholds." + +cmp_n48 <- data.frame( + Metric = c("Status", "Elapsed (min)", "Posterior median $\\rho_B$", + "95% CrI lower (2.5%)", "95% CrI upper (97.5%)", + "ESS_bulk", "$\\hat R$", "Divergent transitions (/2000)", + "Max-treedepth hits (/2000)"), + Phase0 = c("OK", "144.54", "$-0.624$", "$-0.884$", "$-0.378$ ", + "2", "2.824 ", "1", "1999 "), + Phase1 = c("OK", "131.83", "$-0.624$", "$-0.884$", "$-0.378$", + "2", "2.824 ", "1", "1999 "), + Acceptable = c("OK", "—", "near $+0.6$", "—", "should contain $+0.6$", + "$> 400$", "$< 1.01$", "$< 5\\%$", "$< 10\\%$") +) +knitr::kable(cmp_n48, + col.names = c("Metric", "Phase 0 (interactive)", + "Phase 1 (batch)", "Acceptable"), + align = "lrrr") +``` + +Three observations from @tbl-n48 are central to the interpretation in +@sec-discussion: + +1. Phase 0 ≡ Phase 1 at$n = 48$ as well: all summaries match + bit-for-bit. Environment is again ruled out. +2. Posterior has the wrong sign: the true value is $\rho_B = +0.6$ but + the recovered median is $-0.624$. The 95% credible interval + $[-0.884,\,-0.378]$ does not contain the truth. +3. Severe sampler geometry pathology, but few divergences: + $\hat R = 2.82$ and 99.95% of transitions hit the maximum tree + depth, yet the divergent count is only 1 / 2000 (0.05%). This is a + specific signature: the sampler is unable to traverse the typical + set within the step budget, rather than falling into a funnel. + +# Discussion {#sec-discussion} + +## Departure from the over-parameterization pattern + +A naive prior expectation was that $n = 48$ would show weak +identification — the canonical symptom being a posterior that is wide +and roughly centered on the LKJ prior, with $\hat R \approx 1$ (chains +agree they are uninformative) and credible-interval coverage at the +nominal level. + +The observed pattern is qualitatively different: + +| Feature | Standard weak identification | Observed at $n = 48$ | +|------------------------|------------------------|------------------------| +| Posterior shape | Wide, prior-like | Narrow (sd ≈ 0.24), confidently wrong | +| 95% CrI for truth | Wide, covers truth | Narrow, excludes truth | +| Between-chain agreement | $\hat R \approx 1$ | $\hat R = 2.82$ | +| Treedepth hits | Moderate | 99.95% | +| Divergent count | Nonzero | 0.05% | + +The combination near-zero divergent count + near-100% treedepth + +high$\hat R$ points away from "data uninformative for $\rho_B$" and +toward "sampler trapped in a low-information region of parameter space." + +## Three candidate explanations + +After ruling out environment and code-determinism effects via +@sec-results, three candidate mechanisms remain. + +### Structural identifiability of $\Sigma_B$ vs $\Sigma_P$ at finite $n$ + +The Kronecker decomposition is asymptotically identified up to the +scalar ambiguity resolved by the correlation-matrix constraint +(@sec-mv). At infinite $n$ the data separate $\Omega_B$ from $\Omega_P$. +At finite $n$, however, multiple $(\Omega_B, \Omega_P)$ pairs may yield +observationally near-equivalent products $\Omega_B \otimes \Omega_P$. +The boundary of empirical identifiability under our priors and +observation density at $K = 2, P = 5$ may exceed $n = 48$. + +### Likelihood multimodality + +The LKJ(2) prior on a $K = 2$ correlation matrix is unimodal in +$\rho_B$, but the likelihood for $\Omega_B$ given a finite dataset may +not be. +With only two chains, a mode-finding sampler could locate two +distinct modes that yield large $\hat R$ even after warmup, as appears +to be the case here. + +### Simulation–model mismatch + +The least glamorous and most consequential possibility: that +`sim_correlated_case_data()` does not actually produce data drawn from +the exact structure that `inst/stan/model_2.stan` assumes. +Candidate +sources of mismatch include Cholesky-factor ordering conventions, +biomarker-vs-parameter axis ordering in the Kronecker product, or +noise-model parameterization. +This must be checked first, because if +confirmed it would make the other two analyses moot. + +## Open questions + +Distinguishing among the three mechanisms above will require further +investigation, and the appropriate sequencing of those investigations is +itself a question on which advisor input is welcomed. +A natural next +step is a targeted simulation-based calibration workflow: simulate from +known parameter values, refit the model, and verify that posterior +summaries recover the known truth across repeated datasets (Talts et +al., 2018). + +## Scope notes + +- The findings here pertain to Model 2 at $n = 48$ under the specific + prior and simulation configuration described in @sec-methods. + Implications for the Chapter 1 univariate model, if any, are outside + the scope of this document. +- Only one true value of $\rho_B$ ($+0.6$) and only $K = 2$ biomarkers + were tested. Phase 0 and Phase 1 use the same seed. + +# Appendix A — Software environment {.appendix .unnumbered} + +| Component | Version | +|---------------------|---------------------------------------------| +| R | 4.6.0 (2026-04-24) | +| cmdstanr | 0.8.0 | +| cmdstan | 2.38.0 | +| posterior | 1.7.0 | +| shigella | 0.0.0.9006 | +| serodynamics | 0.0.0.9050 | +| Compute environment | UC Davis Shiva HPC (`r_chapter2` conda env) | + +# Appendix B — Reproducibility {.appendix .unnumbered} + +Branch `chapter2-stan-simulation` of +`https://github.com/UCD-SERG/shigella`. + +``` bash +# Phase 0 (interactive SLURM via salloc) +salloc --time=04:00:00 --cpus-per-task=2 --mem=20G +Rscript scripts/phase0_interactive_reproducibility.R # n=5 +Rscript scripts/phase0_interactive_reproducibility_n48.R # n=48 +exit + +# Phase 1 (batch SLURM) +sbatch --time=04:00:00 --cpus-per-task=2 --mem=20G \ + --output=logs/phase1/phase1_%j.out \ + slurm/phase1_single.sbatch # n=5 +sbatch --time=08:00:00 --cpus-per-task=2 --mem=20G \ + --output=logs/phase1/phase1_n48_%j.out \ + slurm/phase1_single_n48.sbatch # n=48 +``` + +Result files: `outputs/phase{0,1}/one_fit_n{5,48}*.rds` and comparison +tables `outputs/phase1/p0_vs_p1_comparison_.rds`. + +# References {.unnumbered} + +::: {#refs} +Betancourt, M. (2018). *A conceptual introduction to Hamiltonian Monte +Carlo*. arXiv:1701.02434. https://doi.org/10.48550/arXiv.1701.02434 + +Betancourt, M., & Girolami, M. (2015). Hamiltonian Monte Carlo for +hierarchical models. In *Current Trends in Bayesian Methodology with +Applications* (pp. 79--101). Chapman & Hall/CRC. + +de Graaf, W. F., Kretzschmar, M. E. E., Teunis, P. F. M., & Diekmann, O. +(2014). *A two-phase within-host model for immune response and its +application to serological profiles of pertussis*. Epidemics, 9, 1--7. +https://doi.org/10.1016/j.epidem.2014.08.002 + +Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & +Rubin, D. B. (2014). *Bayesian Data Analysis* (3rd ed.). Chapman & +Hall/CRC. + +Hoffman, M. D., & Gelman, A. (2014). *The No-U-Turn Sampler: Adaptively +setting path lengths in Hamiltonian Monte Carlo*. Journal of Machine +Learning Research, 15, 1593--1623. + +Lewandowski, D., Kurowicka, D., & Joe, H. (2009). *Generating random +correlation matrices based on vines and extended onion method*. Journal +of Multivariate Analysis, 100(9), 1989--2001. +https://doi.org/10.1016/j.jmva.2009.04.008 + +Mattoo, S., & Cherry, J. D. (2005). *Molecular pathogenesis, +epidemiology, and clinical manifestations of respiratory infections due +to Bordetella pertussis and other Bordetella subspecies*. Clinical +Microbiology Reviews, 18(2), 326--382. +https://doi.org/10.1128/CMR.18.2.326-382.2005 + +Morris, T. P., White, I. R., & Crowther, M. J. (2019). *Using simulation +studies to evaluate statistical methods*. Statistics in Medicine, +38(11), 2074--2102. https://doi.org/10.1002/sim.8086 + +Stavnezer, J., Guikema, J. E. J., & Schrader, C. E. (2008). *Mechanism +and regulation of class switch recombination*. Annual Review of +Immunology, 26, 261--292. +https://doi.org/10.1146/annurev.immunol.26.021607.090248 + +Talts, S., Betancourt, M., Simpson, D., Vehtari, A., & Gelman, A. +(2018). *Validating Bayesian inference algorithms with simulation-based +calibration*. arXiv:1804.06788. +https://doi.org/10.48550/arXiv.1804.06788 + +Teunis, P. F. M., van Eijkeren, J. C. H., de Graaf, W. F., Bonačić +Marinović, A., & Kretzschmar, M. E. E. (2016). *Linking the seroresponse +to infection to within-host heterogeneity in antibody production*. +Epidemics, 16, 33--39. https://doi.org/10.1016/j.epidem.2016.04.001 + +Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. +(2021). *Rank-normalization, folding, and localization: An improved* +$\widehat R$ for assessing convergence of MCMC. Bayesian Analysis, +16(2), 667--718. https://doi.org/10.1214/20-BA1221 +:::