Fix mu_par indexing, add empty data validation, add multi-antigen test

- Fixed mu_par indexing in sample_predictive_stan: use [antigen, param] not [param, antigen]
- Corrected comment to reflect Stan's actual array declaration
- Added validation in prep_data_stan for empty input data
- Added multi-antigen test case for sample_predictive_stan
- Improved test clarity with explicit variable names and comments
- All 63 tests pass (0 errors, 0 warnings)
- R CMD check passes (0 errors, 0 warnings, 0 notes)
- Parallel validation passes (Code Review + CodeQL)

Agent-Logs-Url: https://github.com/ucdavis/serodynamics/sessions/b268f550-e4e5-4afd-b2a0-37c2a0a195d1

Co-authored-by: d-morrison <2474437+d-morrison@users.noreply.github.com>

Author: Copilot

R/prep_data_stan.R

@@ -55,6 +55,17 @@ prep_data_stan <- function(
 
   # Convert to Stan format
   # Stan requires explicit max dimensions
+  # Validate that we have at least one subject with observations
+  if (length(jags_data$nsmpl) == 0 || all(jags_data$nsmpl == 0)) {
+    cli::cli_abort(
+      c(
+        "No observations found in input data.",
+        "i" = "Stan models require at least one subject with observations.",
+        "i" = "Check that your input data is not empty."
+      )
+    )
+  }
+  
   max_nsmpl <- max(jags_data$nsmpl)
 
   # Create padded arrays (Stan doesn't handle ragged arrays like JAGS)

R/sample_predictive_stan.R

@@ -150,13 +150,13 @@ sample_predictive_stan <- function(
     # CmdStan names: y0[subj,k], y1[subj,k], etc.
     # We need to find columns matching pattern for antigen k
     # Extract population-level parameters for this antigen
-    # mu_par has dimensions [param, antigen] where param = 1:5
+    # mu_par has dimensions [antigen, param] where param = 1:5
     # (y0, y1, t1, alpha, shape)
-    y0_pop <- combined_draws[, paste0("mu_par[1,", k, "]")]
-    y1_pop <- combined_draws[, paste0("mu_par[2,", k, "]")]
-    t1_pop <- combined_draws[, paste0("mu_par[3,", k, "]")]
-    alpha_pop <- combined_draws[, paste0("mu_par[4,", k, "]")]
-    shape_pop <- combined_draws[, paste0("mu_par[5,", k, "]")]
+    y0_pop <- combined_draws[, paste0("mu_par[", k, ",1]")]
+    y1_pop <- combined_draws[, paste0("mu_par[", k, ",2]")]
+    t1_pop <- combined_draws[, paste0("mu_par[", k, ",3]")]
+    alpha_pop <- combined_draws[, paste0("mu_par[", k, ",4]")]
+    shape_pop <- combined_draws[, paste0("mu_par[", k, ",5]")]
 
     # Extract measurement error precision for this antigen
     prec_logy_k <- combined_draws[, paste0("prec_logy[", k, "]")]

tests/testthat/test-stan-functions.R

@@ -208,6 +208,74 @@ test_that("sample_predictive_stan works with fitted Stan model", {
   )
 })
 
+test_that("sample_predictive_stan works with multiple antigens", {
+  skip_if_not_installed("cmdstanr")
+  skip_if(
+    is.null(tryCatch(cmdstanr::cmdstan_version(), error = function(e) NULL)),
+    "CmdStan not installed"
+  )
+  
+  withr::local_seed(2)
+  # Use two antigens to test multi-antigen indexing
+  # n=15 provides sufficient data while keeping test runtime reasonable
+  # Two antigens (HlyE_IgA, HlyE_IgG) test parameter extraction for each antigen
+  case_data <- serocalculator::typhoid_curves_nostrat_100 |>
+    sim_case_data(n = 15, antigen_isos = c("HlyE_IgA", "HlyE_IgG"))
+  
+  # Fit model with posterior samples
+  model_output <- run_mod_stan(
+    data = case_data,
+    file_mod = serodynamics_example("model.stan"),
+    nchain = 2,
+    nadapt = 100,
+    niter = 10,
+    with_post = TRUE
+  ) |>
+    suppressWarnings()
+  
+  # Generate predictions
+  n_samples <- 8
+  time_points <- c(10, 50)
+  antigen_isos <- c("HlyE_IgA", "HlyE_IgG")
+  
+  predictions <- sample_predictive_stan(
+    model_output,
+    time_points = time_points,
+    n_samples = n_samples
+  )
+  
+  # Check structure
+  expect_s3_class(predictions, "posterior_predictive_stan")
+  
+  # Check samples array dimensions
+  # Dimensions: [n_samples, n_timepoints, n_antigens]
+  expect_equal(
+    dim(predictions$samples),
+    c(n_samples, length(time_points), length(antigen_isos))
+  )
+  
+  # Check time points
+  expect_equal(predictions$time_points, time_points)
+  
+  # Check summary structure - should have 2 antigens
+  expect_type(predictions$summary, "list")
+  expect_equal(length(predictions$summary), length(antigen_isos))
+  expect_true(all(antigen_isos %in% names(predictions$summary)))
+  
+  # Check that each antigen has summary stats for each time point
+  for (antigen in antigen_isos) {
+    summary_df <- predictions$summary[[antigen]]
+    expect_s3_class(summary_df, "data.frame")
+    expect_equal(nrow(summary_df), length(time_points))
+    expect_true(
+      all(
+        c("time_point", "mean", "median", "lower_95", "upper_95") %in%
+          names(summary_df)
+      )
+    )
+  }
+})
+
 test_that("sample_predictive_stan errors without posterior samples", {
   skip_if_not_installed("cmdstanr")
   skip_if(