Monitor population parameters

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: serodynamics
 Title: What the Package Does (One Line, Title Case)
-Version: 0.0.0.9047
+Version: 0.0.0.9048
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),

NEWS.md

@@ -15,6 +15,7 @@
 
 ## New features
 
+* Including optional population parameters as attributes in run_mod output. (#141)
 * Made "newperson" optional in `prep_data()` (#73)
 * 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)

R/Run_Mod.R

@@ -26,7 +26,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.
-#' @returns An `sr_model` class object: a subclass of [dplyr::tbl_df] that
+#' @param with_pop_params A [logical] value specifying whether population 
+#' level parameters should be included as an attribute entitled 
+#' `population_params`. Excluded as default.
+#' Note: These objects can be large.
+#' @returns An `sr_model` class object: a subclass of [tibble::tbl_df] that
 #' contains MCMC samples from the joint posterior distribution of the model
 #' parameters, conditional on the provided input `data`, 
 #' including the following:
@@ -49,6 +53,17 @@
 #'   - `nIterations`: Number of iteration specified.
 #'   - `nBurnin`: Number of burn ins.
 #'   - `nThin`: Thinning number (niter/nmc).
+#'   - `population_params`: Optionally included modeled population parameters
+#'   indexed by 
+#'   `Iteration`, 
+#'   `Chain`, `Parameter`, `Iso_type`, and `Stratification`. Excluded as 
+#'   default. Includes the following modeled population parameters:
+#'     - `mu.par` = The population means of the host-specific model parameters 
+#'     (on logarithmic scales).
+#'     - `prec.par` = The population precision matrix of the hyperparameters
+#'     (with diagonal elements equal to inverse variances).  
+#'     - `prec.logy` = A vector of population precisions (inverse variances), 
+#'     one per antigen/isotype combination.  
 #'   - `priors`: A [list] that summarizes the input priors, including:
 #'     - `mu_hyp_param`
 #'     - `prec_hyp_param`
@@ -71,6 +86,7 @@ run_mod <- function(data,
                     niter = 100,
                     strat = NA,
                     with_post = FALSE,
+                    with_pop_params = FALSE,
                     ...) {
   ## Conditionally creating a stratification list to loop through
   if (is.na(strat)) {
@@ -80,15 +96,14 @@ run_mod <- function(data,
   }
 
   ## Creating a shell to output results
-  jags_out <- data.frame(
-    "Iteration" = NA,
-    "Chain" = NA,
-    "Parameter" = NA,
-    "value" = NA,
-    "Parameter_sub" = NA,
-    "Subject" = NA,
-    "Iso_type" = NA,
-    "Stratification" = NA
+  jags_out <- tibble::tibble(
+    "Iteration" = integer(),
+    "Chain" = integer(),
+    "value" = numeric(),
+    "Subject" = character(),
+    "Parameter" = character(),
+    "Iso_type" = character(),
+    "Stratification" = character()
   )
 
   ## Creating output list for jags.post
@@ -118,6 +133,11 @@ run_mod <- function(data,
     nthin <- round(niter / nmc) # thinning needed to produce nmc from niter
 
     tomonitor <- c("y0", "y1", "t1", "alpha", "shape")
+    # Conditional statement for including population parameters
+    if (with_pop_params) {
+      tomonitor <- c(tomonitor, "mu.par", "prec.par",
+                     "prec.logy")
+    }
 
     jags_post <- runjags::run.jags(
       model = file_mod,
@@ -138,61 +158,106 @@ run_mod <- function(data,
     jags_post_final[[i]] <- jags_post
 
     # Unpacking and cleaning MCMC output.
-    jags_unpack <- ggmcmc::ggs(jags_post[["mcmc"]])
+    jags_packed <- ggmcmc::ggs(jags_post[["mcmc"]])
 
     # Adding attributes
-    mod_atts <- attributes(jags_unpack)
-    # Only keeping necessary attributes
-    mod_atts <- mod_atts[4:8]
-
+    mod_atts <- attributes(jags_packed)
+    # Select necessary attributes by name for robustness across platforms
+    mod_atts <- mod_atts[c("nChains", "nParameters", "nIterations",
+                           "nBurnin", "nThin", "description")]
+
     # extracting antigen-iso combinations to correctly number
-    # then by the order they are estimated by the program.
+    # them by the order they are estimated by the program.
     iso_dat <- data.frame(attributes(longdata)$antigens)
-    iso_dat <- iso_dat |> dplyr::mutate(Subnum = as.numeric(row.names(iso_dat)))
-    # Working with jags unpacked ggs outputs to clarify parameter and subject
-    jags_unpack <- jags_unpack |>
-      dplyr::mutate(
-        Subnum = sub(".*,", "", .data$Parameter),
-        Parameter_sub = sub("\\[.*", "", .data$Parameter),
-        Subject = sub("\\,.*", "", .data$Parameter)
-      ) |>
-      dplyr::mutate(
-        Subnum = as.numeric(sub("\\].*", "", .data$Subnum)),
-        Subject = sub(".*\\[", "", .data$Subject)
-      )
-    # Merging isodat in to ensure we are classifying antigen_iso
-    jags_unpack <- dplyr::left_join(jags_unpack, iso_dat, by = "Subnum")
+    iso_dat <- iso_dat |> dplyr::mutate(Subnum = row.names(iso_dat))
+
+    # Unpacking the mcmc object
+    jags_unpacked <- unpack_jags(jags_packed)
+
+    # Merging isodat in to ensure we are classifying antigen_iso. 
+    jags_unpacked <- dplyr::left_join(jags_unpacked, iso_dat, 
+                                      by = "Subnum")
+
+    # Adding in ID name
     ids <- data.frame(attr(longdata, "ids")) |>
-      mutate(Subject = as.character(dplyr::row_number()))
-    jags_unpack <- dplyr::left_join(jags_unpack, ids, by = "Subject")
-    jags_final <- jags_unpack |>
-      dplyr::select(!c("Subnum", "Subject")) |>
+      dplyr::mutate(Subject = as.character(dplyr::row_number()))
+    jags_final <- dplyr::left_join(jags_unpacked, ids, 
+                                   by = "Subject") |>
       dplyr::rename(c("Iso_type" = "attributes.longdata..antigens",
-                      "Subject" = "attr.longdata...ids.."))
+                      "Subject_mcmc" = "attr.longdata...ids..")) |>
+      # Subject handling:
+      # * For individual-level parameters, the left join above finds a row in
+      #   attr(longdata, "ids"), so Subject_mcmc contains the subject ID.
+      # * For population-level parameters (e.g., mu.par, prec.par, prec.logy),
+      #   there is no matching ID, so Subject_mcmc is NA. In that case we
+      #   intentionally keep the original Subject value from unpack_jags,
+      #   which holds the parameter name and serves as the identifier.
+      dplyr::mutate(Subject_mcmc = ifelse(is.na(.data$Subject_mcmc), 
+                                          .data$Subject, 
+                                          .data$Subject_mcmc)) |>
+      # At this point, jags_unpacked contains:
+      # * Parameter: original JAGS parameter names (e.g., "y0[1,2]")
+      # * Param: cleaned parameter names used elsewhere in the package.
+      # We drop the original JAGS-style Parameter column and keep the
+      # cleaned names, then rename Param back to Parameter so that all
+      # downstream code consistently uses a single "Parameter" column
+      # containing cleaned parameter names.
+      # Drop the temporary Subject (now only used as a fallback for population
+      # parameters) and rename Subject_mcmc back to Subject for downstream use.
+      dplyr::select(!c("Subnum", "Subject", "Parameter")) |>
+      dplyr::rename(c("Subject" = "Subject_mcmc",
+                      "Parameter" = "Param"))
+
     # Creating a label for the stratification, if there is one.
     # If not, will add in "None".
     jags_final$Stratification <- i
-    ## Creating output
-    jags_out <- data.frame(rbind(jags_out, jags_final))
+    # Creating output as a data frame with the
+    # jags output results for each stratification rbinded.
+    jags_out <- dplyr::bind_rows(jags_out, jags_final)
   }
-  # Ensuring output does not have any NAs
-  jags_out <- jags_out[complete.cases(jags_out), ]
-  # Outputting the finalized jags output as a data frame with the
-  # jags output results for each stratification rbinded.
-
-  # Making output a tibble and restructing.
-  jags_out <- dplyr::as_tibble(jags_out)  |>
-    select(!c("Parameter")) |>
-    rename("Parameter" = "Parameter_sub")
+
+  if (with_pop_params) {
+    # Preparing population parameters
+    population_params <- prep_popparams(jags_out)
+    population_params <- population_params[, c(
+      "Iteration", "Chain", "Parameter", "Iso_type",
+      "Stratification", "Population_Parameter", "value"
+    )]
+
+    # Taking out population parameters
+    jags_out <- ex_popparams(jags_out)
+  }
+
+  # Making output a tibble and restructuring.
   jags_out <- jags_out[, c("Iteration", "Chain", "Parameter", "Iso_type",
                            "Stratification", "Subject", "value")]
-  current_atts <- attributes(jags_out) 
-  current_atts <- c(current_atts, mod_atts)
-  attributes(jags_out) <- current_atts
+  current_atts <- attributes(jags_out)
+  # Explicitly build the attribute list in the correct order to ensure that
+  # `class` appears immediately after `names` and `row.names`.
+  # The dplyr operations above can carry ggmcmc attributes (nChains, etc.)
+  # from jags_packed into jags_out, placing `class` at the end. We use
+  # mod_atts (named selection from jags_packed) as the authoritative source
+  # for the ggmcmc-style metadata attributes.
+  new_atts <- list(
+    names = current_atts$names,
+    row.names = current_atts$row.names,
+    class = union("sr_model", current_atts$class),
+    nChains = mod_atts$nChains,
+    nParameters = mod_atts$nParameters,
+    nIterations = mod_atts$nIterations,
+    nBurnin = mod_atts$nBurnin,
+    nThin = mod_atts$nThin,
+    description = mod_atts$description
+  )
+  attributes(jags_out) <- new_atts
 
-  # Adding priors
+  # Adding population parameters optionally and priors in as attributes
+  if (with_pop_params) {
+    jags_out <- jags_out |>
+      structure(population_params = population_params)
+  }
   jags_out <- jags_out |>
-    structure("priors" = attributes(priorspec)$used_priors)
+    structure(priors = attributes(priorspec)$used_priors)
 
   # Calculating fitted and residuals
   # Renaming columns using attributes from as_case_data
@@ -206,7 +271,5 @@ run_mod <- function(data,
     jags_out <- jags_out |>
       structure(jags.post = jags_post_final)
   }
-  jags_out <- jags_out |>
-    structure(class = union("sr_model", class(jags_out)))
   jags_out
 }

R/as_case_data.R

@@ -51,10 +51,24 @@ as_case_data <- function(
       visit_num = dplyr::row_number()
     ) |>
     serocalculator::set_id_var(id_var) |>
-    structure(
-      class = union("case_data", class(data)),
-      biomarker_var = biomarker_var,
-      timeindays = time_in_days,
-      value_var = value_var
-    )
+    (\(x) {
+      current_atts <- attributes(x)
+      # Explicitly order attributes so that class comes immediately after
+      # names and row.names, which is the expected position for tibble
+      # subclasses. tibble::as_tibble() on a subclass can reorder
+      # standard attributes, so we rebuild them in the correct order.
+      # Use current_atts$id_var to preserve any standardization done by
+      # serocalculator::set_id_var (e.g., renaming the ID column to "id").
+      new_atts <- list(
+        names = current_atts$names,
+        row.names = current_atts$row.names,
+        class = union("case_data", current_atts$class),
+        id_var = current_atts$id_var,
+        biomarker_var = biomarker_var,
+        timeindays = time_in_days,
+        value_var = value_var
+      )
+      attributes(x) <- new_atts
+      x
+    })()
 }

R/ex_popparams.R

@@ -0,0 +1,13 @@
+#' @title Excluding population parameters
+#' @author Sam Schildhauer
+#' @description
+#'  `ex_popparams` excludes estimated population parameters from final output
+#'  [data.frame].
+#' @param x A [data.frame] with a `Subject` variable.
+#' @returns A filtered [data.frame] excluding population parameters.
+#' @keywords internal
+ex_popparams <- function(x) { 
+  x <- x |>
+    dplyr::filter(!(.data$Subject %in% c("mu.par", "prec.par", "prec.logy")))
+  return(x)
+} 

R/nepal_sees_jags_output.R

@@ -8,7 +8,8 @@
 #' which is the diagnosis type (typhoid or
 #' paratyphoid). Keeping only IDs `"newperson"`, `"sees_npl_1"`, `"sees_npl_2"`.
 #'
-#' @format An S3 object of class `sr_model`: A [dplyr::tbl_df] that contains the
+#' @format An S3 object of class `sr_model`: A [tibble::tbl_df] that contains 
+#'   the
 #'   posterior predictive distribution of the person-specific parameters for a
 #'   "new person" with no observed data (`Subject = "newperson"`) and posterior
 #'   distributions of the person-specific parameters for two arbitrarily-chosen

R/param_recode.R

@@ -0,0 +1,17 @@
+#' @title Parameter recode
+#' @author Sam Schildhauer
+#' @description
+#'  `param_recode` recodes character numbers as their corresponding parameter.
+#' @param x A [vector] of character numbers that represent parameters.
+#' @returns A [vector] with recoded values.
+#' @keywords internal
+param_recode <- function(x) {
+  dplyr::recode(
+    x,
+    "1" = "y0",
+    "2" = "y1",
+    "3" = "t1",
+    "4" = "alpha",
+    "5" = "shape"
+  )
+}

R/plot_predicted_curve.R

@@ -7,12 +7,12 @@
 #' and shows all individual 
 #' sampled curves.
 #'
-#' @param model An `sr_model` object (returned by [run_mod()]) containing 
+#' @param model An `sr_model` object (returned by [run_mod]) containing 
 #'   samples from the posterior distribution of the model parameters.
 #' @param ids The participant IDs to plot; for example, `"sees_npl_128"`.
 #' @param antigen_iso  The antigen isotype to plot; for example, "HlyE_IgA" or 
 #' "HlyE_IgG".
-#' @param dataset (Optional) A [dplyr::tbl_df] with observed antibody response 
+#' @param dataset (Optional) A [tibble::tbl_df] with observed antibody response
 #' data. 
 #' Must contain:
 #'   - `timeindays`

R/post_summ.R

@@ -2,7 +2,7 @@
 #' @title Summary Table of Jags Posterior Estimates
 #' @author Sam Schildhauer
 #' @description
-#'  `post_summ()` takes a [list] output from [serodynamics::run_mod()]
+#'  `post_summ()` takes a [list] output from [run_mod]
 #'  to summary table for parameter, antigen/antibody, and stratification
 #'  combination.
 #'  Defaults will produce every combination of antigen/antibody, parameters,
@@ -14,7 +14,7 @@
 #'  - t1 = time to peak
 #'  - r = shape parameter
 #'  - alpha = decay rate
-#' @param data A [list] outputted from [run_mod()].
+#' @param data A [list] outputted from [run_mod].
 #' @param iso Specify [character] string to produce tables of only a
 #' specific antigen/antibody combination, entered with quotes. Default outputs
 #' all antigen/antibody combinations.

R/prep_popparams.R

@@ -0,0 +1,15 @@
+#' @title Preparing population parameters
+#' @author Sam Schildhauer
+#' @description
+#'  `prep_popparams` filters a [data.frame] to only include population 
+#'  parameters and renames the `Subject` variable as `Population_Parameter`.
+#' @param x A [data.frame] with a `Subject` variable.
+#' @returns A filtered [data.frame] with the `Subject` variable renamed to
+#' `Population_Parameter`.
+#' @keywords internal
+prep_popparams <- function(x) { 
+  x <- x |>
+    dplyr::filter(.data$Subject %in% c("mu.par", "prec.par", "prec.logy")) |>
+    dplyr::rename(Population_Parameter = "Subject")
+  return(x)
+}