Add chapter2.qmd: fake-data simulation and correlated model setup

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: serodynamics
 Title: What the Package Does (One Line, Title Case)
-Version: 0.0.0.9043
+Version: 0.0.0.9044
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),
@@ -31,7 +31,8 @@ Imports:
     tibble,
     tidyr,
     tidyselect,
-    utils
+    utils,
+    mvtnorm
 Suggests: 
     Hmisc,
     knitr,

NAMESPACE

@@ -3,9 +3,11 @@
 S3method(autoplot,case_data)
 export(as_case_data)
 export(autoplot)
+export(clean_priors)
 export(expect_snapshot_data)
 export(get_biomarker_levels)
 export(get_biomarker_names_var)
+export(inits_kron)
 export(initsfunction)
 export(load_data)
 export(plot_jags_Rhat)
@@ -17,10 +19,13 @@ export(post_summ)
 export(postprocess_jags_output)
 export(prep_data)
 export(prep_priors)
+export(prep_priors_multi_b)
 export(run_mod)
 export(serodynamics_example)
 export(sim_case_data)
 export(sim_n_obs)
+export(simulate_multi_b_long)
+export(write_model_ch2_kron)
 importFrom(dplyr,.data)
 importFrom(dplyr,all_of)
 importFrom(dplyr,filter)

NEWS.md

@@ -4,6 +4,7 @@
 
 ## New features
 
+* Added Chapter 2 fake-data simulation and correlated Kronecker model (#133)
 * Including fitted and residual values as data frame in run_mod output. (#101)
 * Added  `plot_predicted_curve()` with support for faceting by multiple IDs (#68)
 * Replacing old data object with new run_mod output (#102)

R/Run_Mod.R

@@ -11,6 +11,22 @@
 #'  - t1 = time to peak
 #'  - shape = shape parameter
 #'  - alpha = decay rate
+#' @details
+#'  When `correlated = TRUE`, `run_mod()` fits a Chapter-2 Kronecker prior
+#'  across biomarkers: \eqn{\mathrm{Cov}(\mathrm{vec}(\Theta_i)) = 
+#'  \Sigma_P \otimes \Sigma_B}.
+#'  The likelihood for observed antibody data is unchanged; only the prior
+#'  covariance differs. Internally this mode:
+#'  \itemize{
+#'    \item calls \code{clean_priors()} on the base priors from 
+#'    \code{prep_priors()},
+#'    \item adds Kronecker hyperpriors via \code{prep_priors_multi_b()} 
+#'    (OmegaP, nuP, OmegaB, nuB)
+#'          and \code{n_blocks = <# biomarkers>},
+#'    \item uses \code{inits_kron()} to avoid conflicting legacy inits,
+#'    \item and monitors \code{TauB} and \code{TauP} in addition to the 
+#'    core parameters.
+#'  }
 #' @param data A [base::data.frame()] with the following columns.
 #' @param file_mod The name of the file that contains model structure.
 #' @param nchain An [integer] between 1 and 4 that specifies
@@ -26,6 +42,11 @@
 #' should be included as an element of the [list] object returned by `run_mod()`
 #' (see `Value` section below for details).
 #' Note: These objects can be large.
+#' @param correlated Logical; use Chapter-2 Kronecker prior across biomarkers.
+#'   Default FALSE (independence).
+#' @param file_mod_kron Path to a JAGS file for the Kronecker model.If
+#'   `correlated = TRUE` and this path does not exist, a temporary
+#'   `model_ch2_kron.jags` is written and used automatically.
 #' @returns An `sr_model` class object: a subclass of [dplyr::tbl_df] that
 #' contains MCMC samples from the joint posterior distribution of the model
 #' parameters, conditional on the provided input `data`, 
@@ -60,18 +81,22 @@
 #'   - An optional `"jags.post"` attribute, included when argument
 #'   `with_post` = TRUE.
 #' @inheritDotParams prep_priors
+#' @seealso clean_priors, prep_priors_multi_b, inits_kron, write_model_ch2_kron
 #' @export
 #' @example inst/examples/run_mod-examples.R
 run_mod <- function(data,
-                    file_mod = serodynamics_example("model.jags"),
-                    nchain = 4,
-                    nadapt = 0,
-                    nburn = 0,
-                    nmc = 100,
-                    niter = 100,
-                    strat = NA,
-                    with_post = FALSE,
-                    ...) {
+  file_mod = serodynamics_example("model.jags"),
+  nchain = 4,
+  nadapt = 0,
+  nburn = 0,
+  nmc = 100,
+  niter = 100,
+  strat = NA,
+  with_post = FALSE,
+  ...,
+  correlated = FALSE, # (used when correlated = TRUE)
+  file_mod_kron = "model_ch2_kron.jags"
+) {
   ## Conditionally creating a stratification list to loop through
   if (is.na(strat)) {
     strat_list <- "None"
@@ -106,8 +131,41 @@ run_mod <- function(data,
 
     # prepare data for modeline
     longdata <- prep_data(dl_sub)
-    priorspec <- prep_priors(max_antigens = longdata$n_antigen_isos,
-                             ...)
+
+    # ---------- CHOOSE MODEL/PRIORS DEPENDING ON `correlated` ----------
+    if (!correlated) {
+      # original (independence) behavior
+      priorspec <- prep_priors(
+        max_antigens = longdata$n_antigen_isos, ...
+      )
+      model_path <- file_mod                      # UNCHANGED for independence
+      init_fun   <- initsfunction
+      to_monitor <- c("y0", "y1", "t1", "alpha", "shape")
+    } else {
+      # CH2 behavior: Kronecker prior across biomarkers
+      base_priors <- prep_priors(
+        max_antigens = longdata$n_antigen_isos, ...
+      )
+      base_priors <- serodynamics::clean_priors(base_priors)       
+
+      kron_priors <- serodynamics::prep_priors_multi_b(            
+        n_blocks = longdata$n_antigen_isos
+      )
+      B_scalar   <- list(n_blocks = longdata$n_antigen_isos)       
+
+      priorspec  <- c(base_priors, kron_priors, B_scalar)          
+
+      # Changed: use file_mod_kron when correlated = TRUE
+      model_path <- file_mod_kron                   
+      if (!file.exists(model_path)) {               
+        model_path <- serodynamics::write_model_ch2_kron(
+          file.path(tempdir(), "model_ch2_kron.jags")
+        )
+      }
+
+      init_fun   <- function(chain) serodynamics::inits_kron(chain)    
+      to_monitor <- c("y0", "y1", "t1", "alpha", "shape", "TauB", "TauP") 
+    }
 
     # inputs for jags model
     nchains <- nchain # nr of MC chains to run simultaneously
@@ -117,19 +175,17 @@ run_mod <- function(data,
     niter <- niter # nr of iterations for posterior sample
     nthin <- round(niter / nmc) # thinning needed to produce nmc from niter
 
-    tomonitor <- c("y0", "y1", "t1", "alpha", "shape")
-
     jags_post <- runjags::run.jags(
-      model = file_mod,
+      model = model_path,
       data = c(longdata, priorspec),
-      inits = initsfunction,
+      inits = init_fun,
       method = "parallel",
       adapt = nadapt,
       burnin = nburnin,
       thin = nthin,
       sample = nmc,
       n.chains = nchains,
-      monitor = tomonitor,
+      monitor = to_monitor,
       summarise = FALSE
     )
     # Assigning the raw jags output to a list.
@@ -190,9 +246,11 @@ run_mod <- function(data,
   current_atts <- c(current_atts, mod_atts)
   attributes(jags_out) <- current_atts
 
-  # Adding priors
+  # Changed: Attach priors robustly
+  used_priors <- attr(priorspec, "used_priors", exact = TRUE)
+  if (is.null(used_priors)) used_priors <- names(priorspec)
   jags_out <- jags_out |>
-    structure("priors" = attributes(priorspec)$used_priors)
+    structure("priors" = used_priors)
 
   # Calculating fitted and residuals
   # Renaming columns using attributes from as_case_data

R/clean_priors.R

@@ -0,0 +1,26 @@
+#' @title Drop legacy/unused prior fields
+#' @author Kwan Ho Lee
+#' @description
+#'  `clean_priors()` removes elements from a priors list that are not
+#'  used by the Kronecker model. This helps avoid passing unused or
+#'  conflicting parameters (like `omega`, `wishdf`, or `prec.par`) into JAGS.
+#'
+#' @param x A [base::list()] of priors (e.g., from [prep_priors()]).
+#'
+#' @return A [base::list()] with the legacy fields removed.
+#'
+#' @export
+#' @example inst/examples/examples-clean_priors.R
+clean_priors <- function(x) {
+  drop <- intersect(names(x), c("omega", "wishdf", "Omega", "WishDF", 
+                                "prec.par"))
+  # keep a copy of attributes we care about
+  used <- attr(x, "used_priors", exact = TRUE)
+
+  y <- if (length(drop)) x[setdiff(names(x), drop)] else x
+
+  # restore attribute if it existed
+  if (!is.null(used)) attr(y, "used_priors") <- used
+
+  y
+}

R/inits_kron.R

@@ -0,0 +1,32 @@
+#' @title Safe inits for the Kronecker model
+#' @author Kwan Ho Lee
+#' @description
+#'  `inits_kron()` wraps a base initializer (if provided) and makes sure
+#'  that legacy pieces (`prec.par`) and Kronecker precision terms
+#'  (`TauB`, `TauP`) are not preset. This avoids conflicts when running
+#'  the multi-biomarker Kronecker model.
+#'
+#' @param chain Integer chain index passed through to the base inits function.
+#' @param base_inits A function with signature `function(chain)` that returns a
+#'  named list of initial values. Defaults to a simple RNG seed initializer.
+#'
+#' @return A [base::list()] of inits suitable for \pkg{runjags}.
+#'
+#' @export
+#' @example inst/examples/examples-inits_kron.R
+inits_kron <- function(chain,
+                       base_inits = NULL) {
+  if (is.null(base_inits)) {
+    base_inits <- function(ch) {
+      list(
+        .RNG.name = "base::Mersenne-Twister",
+        .RNG.seed = 123 + as.integer(ch)
+      )
+    }
+  }
+  z <- base_inits(chain)
+  z$TauB <- NULL
+  z$TauP <- NULL
+  z$prec.par <- NULL
+  z
+}

R/prep_priors_multi_b.R

@@ -0,0 +1,53 @@
+#' @title Priors for the Kronecker (multi-biomarker) model
+#' @author Kwan Ho Lee
+#' @description
+#' `prep_priors_multi_b()` builds Wishart hyperparameters for the
+#' within-biomarker precision matrix `T_P` and the across-biomarker precision
+#' matrix `T_B` used in the Kronecker prior `T = T_B %x% T_P`.
+#'
+#' @param n_blocks Integer scalar (B): number of biomarkers.
+#' @param omega_p_scale Numeric length-5 vector for the diagonal of Omega_P
+#'   (parameter scale).
+#' @param nu_p Numeric scalar: degrees of freedom for
+#'   `T_P ~ Wishart(Omega_P, nu_p)`.
+#' @param omega_b_scale Numeric length-`n_blocks` vector for the diagonal of
+#'   Omega_B (biomarker scale).
+#' @param nu_b Numeric scalar: degrees of freedom for
+#'   `T_B ~ Wishart(Omega_B, nu_b)`.
+#'
+#' @return A list with elements `OmegaP`, `nuP`, `OmegaB`, `nuB`.
+#' @export
+#' @example inst/examples/examples-prep_priors_multi_b.R
+prep_priors_multi_b <- function(
+  n_blocks,
+  omega_p_scale = rep(0.1, 5),
+  nu_p = 6,
+  omega_b_scale = rep(1.0, n_blocks),
+  nu_b = n_blocks + 1
+) {
+  # validations (cli::cli_abort is linter-approved)
+  if (!is.numeric(n_blocks) || length(n_blocks) != 1L || n_blocks < 1 ||
+        !isTRUE(all.equal(n_blocks, as.integer(n_blocks)))) {
+    cli::cli_abort("`n_blocks` must be a positive integer of length 1.")
+  }
+  if (!is.numeric(omega_p_scale) || length(omega_p_scale) != 5L) {
+    cli::cli_abort("`omega_p_scale` must be a numeric vector of length 5.")
+  }
+  if (!is.numeric(omega_b_scale) || length(omega_b_scale) != n_blocks) {
+    cli::cli_abort("`omega_b_scale` must be a numeric vector of length 
+                   `n_blocks`.")
+  }
+  if (!is.numeric(nu_p) || length(nu_p) != 1L) {
+    cli::cli_abort("`nu_p` must be a numeric scalar.")
+  }
+  if (!is.numeric(nu_b) || length(nu_b) != 1L) {
+    cli::cli_abort("`nu_b` must be a numeric scalar.")
+  }
+
+  list(
+    OmegaP = diag(omega_p_scale, nrow = 5),
+    nuP    = nu_p,
+    OmegaB = diag(omega_b_scale, nrow = n_blocks),
+    nuB    = nu_b
+  )
+}