diff --git a/DESCRIPTION b/DESCRIPTION
index ea5d92de..a46b1417 100644
--- a/DESCRIPTION
+++ b/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."),
diff --git a/NEWS.md b/NEWS.md
index 89199522..5f0fe04d 100644
--- a/NEWS.md
+++ b/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)
diff --git a/R/Run_Mod.R b/R/Run_Mod.R
index 25e75d60..48d202bf 100644
--- a/R/Run_Mod.R
+++ b/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
 }
diff --git a/R/as_case_data.R b/R/as_case_data.R
index eb55b140..a1d66259 100644
--- a/R/as_case_data.R
+++ b/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
+    })()
 }
diff --git a/R/ex_popparams.R b/R/ex_popparams.R
new file mode 100644
index 00000000..7a565776
--- /dev/null
+++ b/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)
+} 
diff --git a/R/nepal_sees_jags_output.R b/R/nepal_sees_jags_output.R
index fb75a06c..eac82eda 100644
--- a/R/nepal_sees_jags_output.R
+++ b/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
diff --git a/R/param_recode.R b/R/param_recode.R
new file mode 100644
index 00000000..de05454b
--- /dev/null
+++ b/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"
+  )
+}
diff --git a/R/plot_predicted_curve.R b/R/plot_predicted_curve.R
index 3a0d4654..ae782595 100644
--- a/R/plot_predicted_curve.R
+++ b/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`
diff --git a/R/post_summ.R b/R/post_summ.R
index 011adc35..c7f7b43f 100644
--- a/R/post_summ.R
+++ b/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.
diff --git a/R/prep_popparams.R b/R/prep_popparams.R
new file mode 100644
index 00000000..8ecb6751
--- /dev/null
+++ b/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)
+} 
diff --git a/R/unpack_jags.R b/R/unpack_jags.R
new file mode 100644
index 00000000..84bef9c7
--- /dev/null
+++ b/R/unpack_jags.R
@@ -0,0 +1,100 @@
+#' @title Unpacking MCMC Object
+#' @author Sam Schildhauer
+#' @description
+#'  `unpack_jags()` takes a long-format MCMC sample (typically created by
+#'  applying [ggmcmc::ggs()] to the `mcmc` component of [run_mod] output)
+#'  and unpacks it into separate rows for individual-level curve parameters
+#'  and population-level hyperparameters/precision terms.
+#' @param data A [tibble::tbl_df] in [ggmcmc::ggs()] / MCMC-long format,
+#'   usually `ggmcmc::ggs(jags_post[["mcmc"]])` where `jags_post` comes from
+#'   [run_mod]. Must contain at least `Iteration`, `Chain`, `Parameter`,
+#'   and `value` columns.
+#' @returns A [tibble::tbl_df] that
+#' contains MCMC samples from the joint posterior distribution of the model
+#' with unpacked individual-level parameters (e.g., `y0`, `y1`, `t1`,
+#' `alpha`, `shape`) and population-level parameters (e.g., `mu.par`,
+#' `prec.par`, `prec.logy`), along with isotypes and subjects.
+#' @keywords internal
+unpack_jags <- function(data) {
+  
+  # Convert Parameter column to character if it's a factor (from ggmcmc::ggs)
+  if (is.factor(data$Parameter)) {
+    data$Parameter <- as.character(data$Parameter)
+  }
+
+  unpack_with_pattern <- function(data, filter_pattern, regex_pattern,
+                                  subject_repl, subnum_repl, param_fun) {
+    data |>
+      dplyr::filter(startsWith(.data$Parameter, paste0(filter_pattern, "["))) |>
+      dplyr::mutate(
+        Subject = gsub(regex_pattern, subject_repl, .data$Parameter),
+        Subnum = gsub(regex_pattern, subnum_repl, .data$Parameter),
+        Param = param_fun(.data$Parameter, regex_pattern)
+      )
+  }
+
+  # Regular expressions for unpacking
+  regex_twoidx <- "([[:alnum:].]+)\\[([0-9]+),([0-9]+)\\]"       # e.g. x[1,2]
+  regex_threeidx <- "([[:alnum:].]+)\\[([0-9]+),([0-9]+),([0-9]+)\\]"
+  # e.g. x[1,2,3]
+  regex_oneidx <- "([[:alnum:].]+)\\[([0-9]+)\\]"                 # e.g. x[1]
+
+  # Unpacking mu.par
+  # Separating population parameters from the rest of the data
+  jags_mupar <- unpack_with_pattern(
+    data = data,
+    filter_pattern = "mu.par",
+    regex_pattern = regex_twoidx,
+    subject_repl = "\\1",
+    subnum_repl = "\\2",
+    param_fun = function(param, pattern) {
+      param_recode(gsub(pattern, "\\3", param))
+    }
+  )
+
+  # Unpacking prec.par
+  jags_precpar <- unpack_with_pattern(
+    data = data,
+    filter_pattern = "prec.par",
+    regex_pattern = regex_threeidx,
+    subject_repl = "\\1",
+    subnum_repl = "\\2",
+    param_fun = function(param, pattern) {
+      paste0(
+        param_recode(gsub(pattern, "\\3", param)), ", ",
+        param_recode(gsub(pattern, "\\4", param))
+      )
+    }
+  )
+
+  # Unpacking prec.logy
+  jags_preclogy <- unpack_with_pattern(
+    data = data,
+    filter_pattern = "prec.logy",
+    regex_pattern = regex_oneidx,
+    subject_repl = "\\1",
+    subnum_repl = "\\2",
+    param_fun = function(param, pattern) {
+      "prec.logy"
+    }
+  )
+
+  # Working with jags unpacked ggs outputs to clarify parameter and subject
+  jags_unpack_params <- data |>
+    dplyr::mutate(
+      Subject = gsub(regex_twoidx, "\\2", .data$Parameter),
+      Subnum = gsub(regex_twoidx, "\\3", .data$Parameter),
+      Param = gsub(regex_twoidx, "\\1", .data$Parameter)
+    ) |> 
+    dplyr::filter(.data$Param %in% c("y0", "y1", "t1", "alpha", "shape"))
+
+  # Putting data frame together
+  jags_unpack_bind <- dplyr::bind_rows(
+    jags_unpack_params,
+    jags_mupar,
+    jags_precpar,
+    jags_preclogy
+  )
+
+  return(jags_unpack_bind)
+}
diff --git a/R/use_att_names.R b/R/use_att_names.R
index b7212858..6f651542 100644
--- a/R/use_att_names.R
+++ b/R/use_att_names.R
@@ -2,7 +2,7 @@
 #' @description
 #'  `use_att_names` takes prepared longitudinal data for antibody kinetic
 #'  modeling and names columns using attribute values to allow merging
-#'  with a modeled [run_mod()] output [dplyr::tbl_df]. The column names include
+#'  with a modeled [run_mod] output [tibble::tbl_df]. The column names include
 #'  `Subject`, `Iso_type`, `t`, and `result`. 
 #' @param data A [data.frame] raw longitudinal data that has been
 #' prepared for antibody kinetic modeling using [as_case_data()].
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 8ea3768d..a0716aa2 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -43,6 +43,7 @@ onwards
 params
 pkgdown
 prec
+precisions
 preinstalling
 rhat
 seroconversion
diff --git a/man/ex_popparams.Rd b/man/ex_popparams.Rd
new file mode 100644
index 00000000..b389d835
--- /dev/null
+++ b/man/ex_popparams.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ex_popparams.R
+\name{ex_popparams}
+\alias{ex_popparams}
+\title{Excluding population parameters}
+\usage{
+ex_popparams(x)
+}
+\arguments{
+\item{x}{A \link{data.frame} with a \code{Subject} variable.}
+}
+\value{
+A filtered \link{data.frame} excluding population parameters.
+}
+\description{
+\code{ex_popparams} excludes estimated population parameters from final output
+\link{data.frame}.
+}
+\author{
+Sam Schildhauer
+}
+\keyword{internal}
diff --git a/man/nepal_sees_jags_output.Rd b/man/nepal_sees_jags_output.Rd
index 2bd13915..222a6e09 100644
--- a/man/nepal_sees_jags_output.Rd
+++ b/man/nepal_sees_jags_output.Rd
@@ -5,7 +5,8 @@
 \alias{nepal_sees_jags_output}
 \title{SEES Typhoid run_mod jags output}
 \format{
-An S3 object of class \code{sr_model}: A \link[dplyr:tbl_df]{dplyr::tbl_df} that contains the
+An S3 object of class \code{sr_model}: A \link[tibble:tbl_df-class]{tibble::tbl_df} that contains
+the
 posterior predictive distribution of the person-specific parameters for a
 "new person" with no observed data (\code{Subject = "newperson"}) and posterior
 distributions of the person-specific parameters for two arbitrarily-chosen
diff --git a/man/param_recode.Rd b/man/param_recode.Rd
new file mode 100644
index 00000000..60e648b0
--- /dev/null
+++ b/man/param_recode.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/param_recode.R
+\name{param_recode}
+\alias{param_recode}
+\title{Parameter recode}
+\usage{
+param_recode(x)
+}
+\arguments{
+\item{x}{A \link{vector} of character numbers that represent parameters.}
+}
+\value{
+A \link{vector} with recoded values.
+}
+\description{
+\code{param_recode} recodes character numbers as their corresponding parameter.
+}
+\author{
+Sam Schildhauer
+}
+\keyword{internal}
diff --git a/man/plot_predicted_curve.Rd b/man/plot_predicted_curve.Rd
index 7e2fabde..0bb77d83 100644
--- a/man/plot_predicted_curve.Rd
+++ b/man/plot_predicted_curve.Rd
@@ -23,7 +23,7 @@ plot_predicted_curve(
 )
 }
 \arguments{
-\item{model}{An \code{sr_model} object (returned by \code{\link[=run_mod]{run_mod()}}) containing
+\item{model}{An \code{sr_model} object (returned by \link{run_mod}) containing
 samples from the posterior distribution of the model parameters.}
 
 \item{ids}{The participant IDs to plot; for example, \code{"sees_npl_128"}.}
@@ -31,7 +31,7 @@ samples from the posterior distribution of the model parameters.}
 \item{antigen_iso}{The antigen isotype to plot; for example, "HlyE_IgA" or
 "HlyE_IgG".}
 
-\item{dataset}{(Optional) A \link[dplyr:tbl_df]{dplyr::tbl_df} with observed antibody response
+\item{dataset}{(Optional) A \link[tibble:tbl_df-class]{tibble::tbl_df} with observed antibody response
 data.
 Must contain:
 \itemize{
diff --git a/man/post_summ.Rd b/man/post_summ.Rd
index 54980e08..239faa1b 100644
--- a/man/post_summ.Rd
+++ b/man/post_summ.Rd
@@ -12,7 +12,7 @@ post_summ(
 )
 }
 \arguments{
-\item{data}{A \link{list} outputted from \code{\link[=run_mod]{run_mod()}}.}
+\item{data}{A \link{list} outputted from \link{run_mod}.}
 
 \item{iso}{Specify \link{character} string to produce tables of only a
 specific antigen/antibody combination, entered with quotes. Default outputs
@@ -36,7 +36,7 @@ A \link{data.frame} summarizing estimate mean, standard deviation (SD),
 median, and quantiles (2.5\%, 25.0\%, 50.0\%, 75.0\%, 97.5\%).
 }
 \description{
-\code{post_summ()} takes a \link{list} output from \code{\link[=run_mod]{run_mod()}}
+\code{post_summ()} takes a \link{list} output from \link{run_mod}
 to summary table for parameter, antigen/antibody, and stratification
 combination.
 Defaults will produce every combination of antigen/antibody, parameters,
diff --git a/man/prep_popparams.Rd b/man/prep_popparams.Rd
new file mode 100644
index 00000000..b4f4d3d2
--- /dev/null
+++ b/man/prep_popparams.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/prep_popparams.R
+\name{prep_popparams}
+\alias{prep_popparams}
+\title{Preparing population parameters}
+\usage{
+prep_popparams(x)
+}
+\arguments{
+\item{x}{A \link{data.frame} with a \code{Subject} variable.}
+}
+\value{
+A filtered \link{data.frame} with the \code{Subject} variable renamed to
+\code{Population_Parameter}.
+}
+\description{
+\code{prep_popparams} filters a \link{data.frame} to only include population
+parameters and renames the \code{Subject} variable as \code{Population_Parameter}.
+}
+\author{
+Sam Schildhauer
+}
+\keyword{internal}
diff --git a/man/run_mod.Rd b/man/run_mod.Rd
index 08f91427..78db7781 100644
--- a/man/run_mod.Rd
+++ b/man/run_mod.Rd
@@ -14,6 +14,7 @@ run_mod(
   niter = 100,
   strat = NA,
   with_post = FALSE,
+  with_pop_params = FALSE,
   ...
 )
 }
@@ -42,6 +43,11 @@ should be included as an element of the \link{list} object returned by \code{run
 (see \code{Value} section below for details).
 Note: These objects can be large.}
 
+\item{with_pop_params}{A \link{logical} value specifying whether population
+level parameters should be included as an attribute entitled
+\code{population_params}. Excluded as default.
+Note: These objects can be large.}
+
 \item{...}{
   Arguments passed on to \code{\link[=prep_priors]{prep_priors}}
   \describe{
@@ -93,7 +99,7 @@ If specified, must be 2 values long:
   }}
 }
 \value{
-An \code{sr_model} class object: a subclass of \link[dplyr:tbl_df]{dplyr::tbl_df} that
+An \code{sr_model} class object: a subclass of \link[tibble:tbl_df-class]{tibble::tbl_df} that
 contains MCMC samples from the joint posterior distribution of the model
 parameters, conditional on the provided input \code{data},
 including the following:
@@ -120,6 +126,19 @@ including the following:
 \item \code{nIterations}: Number of iteration specified.
 \item \code{nBurnin}: Number of burn ins.
 \item \code{nThin}: Thinning number (niter/nmc).
+\item \code{population_params}: Optionally included modeled population parameters
+indexed by
+\code{Iteration},
+\code{Chain}, \code{Parameter}, \code{Iso_type}, and \code{Stratification}. Excluded as
+default. Includes the following modeled population parameters:
+\itemize{
+\item \code{mu.par} = The population means of the host-specific model parameters
+(on logarithmic scales).
+\item \code{prec.par} = The population precision matrix of the hyperparameters
+(with diagonal elements equal to inverse variances).
+\item \code{prec.logy} = A vector of population precisions (inverse variances),
+one per antigen/isotype combination.
+}
 \item \code{priors}: A \link{list} that summarizes the input priors, including:
 \itemize{
 \item \code{mu_hyp_param}
diff --git a/man/unpack_jags.Rd b/man/unpack_jags.Rd
new file mode 100644
index 00000000..f4a41c48
--- /dev/null
+++ b/man/unpack_jags.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/unpack_jags.R
+\name{unpack_jags}
+\alias{unpack_jags}
+\title{Unpacking MCMC Object}
+\usage{
+unpack_jags(data)
+}
+\arguments{
+\item{data}{A \link[tibble:tbl_df-class]{tibble::tbl_df} in \code{\link[ggmcmc:ggs]{ggmcmc::ggs()}} / MCMC-long format,
+usually \code{ggmcmc::ggs(jags_post[["mcmc"]])} where \code{jags_post} comes from
+\link{run_mod}. Must contain at least \code{Iteration}, \code{Chain}, \code{Parameter},
+and \code{value} columns.}
+}
+\value{
+A \link[tibble:tbl_df-class]{tibble::tbl_df} that
+contains MCMC samples from the joint posterior distribution of the model
+with unpacked individual-level parameters (e.g., \code{y0}, \code{y1}, \code{t1},
+\code{alpha}, \code{shape}) and population-level parameters (e.g., \code{mu.par},
+\code{prec.par}, \code{prec.logy}), along with isotypes and subjects.
+}
+\description{
+\code{unpack_jags()} takes a long-format MCMC sample (typically created by
+applying \code{\link[ggmcmc:ggs]{ggmcmc::ggs()}} to the \code{mcmc} component of \link{run_mod} output)
+and unpacks it into separate rows for individual-level curve parameters
+and population-level hyperparameters/precision terms.
+}
+\author{
+Sam Schildhauer
+}
+\keyword{internal}
diff --git a/man/use_att_names.Rd b/man/use_att_names.Rd
index f43fcef0..5367bf7f 100644
--- a/man/use_att_names.Rd
+++ b/man/use_att_names.Rd
@@ -16,7 +16,7 @@ The input \link{data.frame} with columns named after attributes.
 \description{
 \code{use_att_names} takes prepared longitudinal data for antibody kinetic
 modeling and names columns using attribute values to allow merging
-with a modeled \code{\link[=run_mod]{run_mod()}} output \link[dplyr:tbl_df]{dplyr::tbl_df}. The column names include
+with a modeled \link{run_mod} output \link[tibble:tbl_df-class]{tibble::tbl_df}. The column names include
 \code{Subject}, \code{Iso_type}, \code{t}, and \code{result}.
 }
 \keyword{internal}
diff --git a/tests/testthat/_snaps/as_case_data.md b/tests/testthat/_snaps/as_case_data.md
index 44130b20..f14b4b05 100644
--- a/tests/testthat/_snaps/as_case_data.md
+++ b/tests/testthat/_snaps/as_case_data.md
@@ -9,7 +9,7 @@
 
     list(names = c("Country", "id", "sample_id", "bldculres", "antigen_iso", 
     "studyvisit", "dayssincefeveronset", "result", "visit_num"), 
-        class = c("case_data", "spec_tbl_df", "tbl_df", "tbl", "data.frame"
-        ), id_var = "id", biomarker_var = "antigen_iso", timeindays = "dayssincefeveronset", 
+        class = c("case_data", "tbl_df", "tbl", "data.frame"), id_var = "id", 
+        biomarker_var = "antigen_iso", timeindays = "dayssincefeveronset", 
         value_var = "result")
 
diff --git a/tests/testthat/_snaps/darwin/run_mod/popparam-summary-stats.csv b/tests/testthat/_snaps/darwin/run_mod/popparam-summary-stats.csv
new file mode 100644
index 00000000..8a7b29a5
--- /dev/null
+++ b/tests/testthat/_snaps/darwin/run_mod/popparam-summary-stats.csv
@@ -0,0 +1,31 @@
+Parameter,mean,sd
+alpha,-4.8276,0.200665
+"alpha, alpha",2.56363,0.368743
+"alpha, shape",-0.610002,0.872237
+"alpha, t1",0.253244,0.910382
+"alpha, y0",-0.269962,0.623738
+"alpha, y1",-0.052442,0.187775
+shape,-0.627015,0.0691708
+"shape, alpha",-0.610002,0.872237
+"shape, shape",9.52829,2.70326
+"shape, t1",4.04278,3.12409
+"shape, y0",-2.43068,3.82725
+"shape, y1",-1.03352,0.430277
+t1,0.463903,0.0405599
+"t1, alpha",0.253244,0.910382
+"t1, shape",4.04278,3.12409
+"t1, t1",18.8339,5.49919
+"t1, y0",0.778855,4.44261
+"t1, y1",0.0564977,0.659988
+y0,0.78573,0.110742
+"y0, alpha",-0.269962,0.623738
+"y0, shape",-2.43068,3.82725
+"y0, t1",0.778855,4.44261
+"y0, y0",11.1517,4.71125
+"y0, y1",0.230998,0.401854
+y1,5.87676,0.185882
+"y1, alpha",-0.052442,0.187775
+"y1, shape",-1.03352,0.430277
+"y1, t1",0.0564977,0.659988
+"y1, y0",0.230998,0.401854
+"y1, y1",0.706827,0.126534
diff --git a/tests/testthat/_snaps/run_mod.md b/tests/testthat/_snaps/run_mod.md
index 14cf7042..7a53344b 100644
--- a/tests/testthat/_snaps/run_mod.md
+++ b/tests/testthat/_snaps/run_mod.md
@@ -1,64 +1,20 @@
-# results are consistent with simulated data
-
-    list(names = c("Iteration", "Chain", "Parameter", "Iso_type", 
-    "Stratification", "Subject", "value"), class = c("sr_model", 
-    "tbl_df", "tbl", "data.frame"), nChains = 2L, nParameters = 505L, 
-        nIterations = 10L, nBurnin = 200, nThin = 1, priors = list(
-            mu_hyp_param = c(1, 7, 1, -4, -1), prec_hyp_param = c(1, 
-            1e-05, 1, 0.001, 1), omega_param = c(1, 50, 1, 10, 1), 
-            wishdf_param = 20, prec_logy_hyp_param = c(4, 1)))
-
 # results are consistent with SEES data
 
     list(names = c("Iteration", "Chain", "Parameter", "Iso_type", 
     "Stratification", "Subject", "value"), class = c("sr_model", 
     "tbl_df", "tbl", "data.frame"), nChains = 2L, nParameters = 440L, 
-        nIterations = 100L, nBurnin = 20, nThin = 1, priors = list(
-            mu_hyp_param = c(1, 7, 1, -4, -1), prec_hyp_param = c(1, 
-            1e-05, 1, 0.001, 1), omega_param = c(1, 50, 1, 10, 1), 
-            wishdf_param = 20, prec_logy_hyp_param = c(4, 1)))
+        nIterations = 100L, nBurnin = 20, nThin = 1, description = "jags_post[[\"mcmc\"]]", 
+        priors = list(mu_hyp_param = c(1, 7, 1, -4, -1), prec_hyp_param = c(1, 
+        1e-05, 1, 0.001, 1), omega_param = c(1, 50, 1, 10, 1), wishdf_param = 20, 
+            prec_logy_hyp_param = c(4, 1)))
 
 # results are consistent with unstratified SEES data
 
     list(names = c("Iteration", "Chain", "Parameter", "Iso_type", 
     "Stratification", "Subject", "value"), class = c("sr_model", 
     "tbl_df", "tbl", "data.frame"), nChains = 2L, nParameters = 1890L, 
-        nIterations = 100L, nBurnin = 20, nThin = 1, priors = list(
-            mu_hyp_param = c(1, 7, 1, -4, -1), prec_hyp_param = c(1, 
-            1e-05, 1, 0.001, 1), omega_param = c(1, 50, 1, 10, 1), 
-            wishdf_param = 20, prec_logy_hyp_param = c(4, 1)))
-
-# results are consistent with unstratified SEES data with jags.post  included
-
-    WAoAAAACAAQEAgACAwAAAAITAAAACAAAABAAAAAHAAQACQAAAAlJdGVyYXRpb24ABAAJAAAA
-    BUNoYWluAAQACQAAAAlQYXJhbWV0ZXIABAAJAAAACElzb190eXBlAAQACQAAAA5TdHJhdGlm
-    aWNhdGlvbgAEAAkAAAAHU3ViamVjdAAEAAkAAAAFdmFsdWUAAAAQAAAABAAEAAkAAAAIc3Jf
-    bW9kZWwABAAJAAAABnRibF9kZgAEAAkAAAADdGJsAAQACQAAAApkYXRhLmZyYW1lAAAADQAA
-    AAEAAAACAAAADQAAAAEAAAdiAAAADQAAAAEAAABkAAAADgAAAAFANAAAAAAAAAAAAA4AAAAB
-    P/AAAAAAAAAAAAITAAAABQAAAA4AAAAFP/AAAAAAAABAHAAAAAAAAD/wAAAAAAAAwBAAAAAA
-    AAC/8AAAAAAAAAAAAA4AAAAFP/AAAAAAAAA+5Pi1iONo8T/wAAAAAAAAP1BiTdLxqfw/8AAA
-    AAAAAAAAAA4AAAAFP/AAAAAAAABASQAAAAAAAD/wAAAAAAAAQCQAAAAAAAA/8AAAAAAAAAAA
-    AA4AAAABQDQAAAAAAAAAAAAOAAAAAkAQAAAAAAAAP/AAAAAAAAAAAAQCAAAAAQAEAAkAAAAF
-    bmFtZXMAAAAQAAAABQAEAAkAAAAMbXVfaHlwX3BhcmFtAAQACQAAAA5wcmVjX2h5cF9wYXJh
-    bQAEAAkAAAALb21lZ2FfcGFyYW0ABAAJAAAADHdpc2hkZl9wYXJhbQAEAAkAAAATcHJlY19s
-    b2d5X2h5cF9wYXJhbQAAAP4AAAQCAAAB/wAAABAAAAAIAAQACQAAAAVuYW1lcwAEAAkAAAAF
-    Y2xhc3MABAAJAAAAB25DaGFpbnMABAAJAAAAC25QYXJhbWV0ZXJzAAQACQAAAAtuSXRlcmF0
-    aW9ucwAEAAkAAAAHbkJ1cm5pbgAEAAkAAAAFblRoaW4ABAAJAAAABnByaW9ycwAAAP4=
-
-# results are consistent with unstratified SEES data with modified   priors
-
-    WAoAAAACAAQEAgACAwAAAAITAAAACAAAABAAAAAHAAQACQAAAAlJdGVyYXRpb24ABAAJAAAA
-    BUNoYWluAAQACQAAAAlQYXJhbWV0ZXIABAAJAAAACElzb190eXBlAAQACQAAAA5TdHJhdGlm
-    aWNhdGlvbgAEAAkAAAAHU3ViamVjdAAEAAkAAAAFdmFsdWUAAAAQAAAABAAEAAkAAAAIc3Jf
-    bW9kZWwABAAJAAAABnRibF9kZgAEAAkAAAADdGJsAAQACQAAAApkYXRhLmZyYW1lAAAADQAA
-    AAEAAAACAAAADQAAAAEAAAdiAAAADQAAAAEAAABkAAAADgAAAAFANAAAAAAAAAAAAA4AAAAB
-    P/AAAAAAAAAAAAITAAAABQAAAA4AAAAFP/AAAAAAAABAEAAAAAAAAD/wAAAAAAAAwAgAAAAA
-    AAC/8AAAAAAAAAAAAA4AAAAFP4R64UeuFHs/Gjbi6xxDLT+EeuFHrhR7P1BiTdLxqfw/hHrh
-    R64UewAAAA4AAAAFP/AAAAAAAABANAAAAAAAAD/wAAAAAAAAQCQAAAAAAAA/8AAAAAAAAAAA
-    AA4AAAABQCQAAAAAAAAAAAAOAAAAAkAIAAAAAAAAP/AAAAAAAAAAAAQCAAAAAQAEAAkAAAAF
-    bmFtZXMAAAAQAAAABQAEAAkAAAAMbXVfaHlwX3BhcmFtAAQACQAAAA5wcmVjX2h5cF9wYXJh
-    bQAEAAkAAAALb21lZ2FfcGFyYW0ABAAJAAAADHdpc2hkZl9wYXJhbQAEAAkAAAATcHJlY19s
-    b2d5X2h5cF9wYXJhbQAAAP4AAAQCAAAB/wAAABAAAAAIAAQACQAAAAVuYW1lcwAEAAkAAAAF
-    Y2xhc3MABAAJAAAAB25DaGFpbnMABAAJAAAAC25QYXJhbWV0ZXJzAAQACQAAAAtuSXRlcmF0
-    aW9ucwAEAAkAAAAHbkJ1cm5pbgAEAAkAAAAFblRoaW4ABAAJAAAABnByaW9ycwAAAP4=
+        nIterations = 100L, nBurnin = 20, nThin = 1, description = "jags_post[[\"mcmc\"]]", 
+        priors = list(mu_hyp_param = c(1, 7, 1, -4, -1), prec_hyp_param = c(1, 
+        1e-05, 1, 0.001, 1), omega_param = c(1, 50, 1, 10, 1), wishdf_param = 20, 
+            prec_logy_hyp_param = c(4, 1)))
 
diff --git a/tests/testthat/_snaps/run_mod/popparam-summary-stats.csv b/tests/testthat/_snaps/run_mod/popparam-summary-stats.csv
new file mode 100644
index 00000000..18976545
--- /dev/null
+++ b/tests/testthat/_snaps/run_mod/popparam-summary-stats.csv
@@ -0,0 +1,31 @@
+Parameter,mean,sd
+alpha,-4.73607,0.219781
+"alpha, alpha",2.76907,0.688682
+"alpha, shape",-0.383234,0.533813
+"alpha, t1",0.448772,1.57899
+"alpha, y0",0.700184,0.757452
+"alpha, y1",-0.0212652,0.222195
+shape,-0.681647,0.0690278
+"shape, alpha",-0.383234,0.533813
+"shape, shape",6.58223,2.3041
+"shape, t1",0.914755,2.42633
+"shape, y0",-2.45229,1.61395
+"shape, y1",-0.877117,0.319837
+t1,0.591297,0.100501
+"t1, alpha",0.448772,1.57899
+"t1, shape",0.914755,2.42633
+"t1, t1",20.3291,6.01563
+"t1, y0",3.25694,6.00627
+"t1, y1",0.0188795,0.473277
+y0,0.726466,0.0841371
+"y0, alpha",0.700184,0.757452
+"y0, shape",-2.45229,1.61395
+"y0, t1",3.25694,6.00627
+"y0, y0",13.2226,7.023
+"y0, y1",0.220216,0.552936
+y1,5.7489,0.348058
+"y1, alpha",-0.0212652,0.222195
+"y1, shape",-0.877117,0.319837
+"y1, t1",0.0188795,0.473277
+"y1, y0",0.220216,0.552936
+"y1, y1",0.681022,0.11074
diff --git a/tests/testthat/test-internal-popparams.R b/tests/testthat/test-internal-popparams.R
new file mode 100644
index 00000000..2cdf1031
--- /dev/null
+++ b/tests/testthat/test-internal-popparams.R
@@ -0,0 +1,173 @@
+test_that("param_recode correctly recodes parameter indices", {
+  # Test valid parameter indices
+  expect_equal(serodynamics:::param_recode("1"), "y0")
+  expect_equal(serodynamics:::param_recode("2"), "y1")
+  expect_equal(serodynamics:::param_recode("3"), "t1")
+  expect_equal(serodynamics:::param_recode("4"), "alpha")
+  expect_equal(serodynamics:::param_recode("5"), "shape")
+  
+  # Test vector of valid indices
+  expect_equal(
+    serodynamics:::param_recode(c("1", "2", "3", "4", "5")),
+    c("y0", "y1", "t1", "alpha", "shape")
+  )
+})
+
+test_that("unpack_jags handles factor Parameter column", {
+  # Create a simple ggs-like object with factor Parameter column
+  test_data <- tibble::tibble(
+    Iteration = rep(1:5, 3),
+    Chain = rep(1, 15),
+    Parameter = factor(
+      rep(c("mu.par[1,1]", "prec.par[1,1,2]", "y0[1,1]"), each = 5)
+    ),
+    value = rnorm(15)
+  )
+  
+  # Should not error with factor column
+  result <- serodynamics:::unpack_jags(test_data)
+  
+  expect_true(is.data.frame(result))
+  expect_true("Subject" %in% names(result))
+  expect_true("Param" %in% names(result))
+})
+
+test_that("unpack_jags correctly unpacks mu.par parameters", {
+  test_data <- tibble::tibble(
+    Iteration = rep(1:3, 2),
+    Chain = rep(1, 6),
+    Parameter = rep(c("mu.par[1,1]", "mu.par[1,2]"), each = 3),
+    value = rnorm(6)
+  )
+  
+  result <- serodynamics:::unpack_jags(test_data)
+  
+  # Check that mu.par rows are present
+  mu_rows <- result[result$Subject == "mu.par", ]
+  expect_gt(nrow(mu_rows), 0)
+  
+  # Check Subject is correctly set
+  expect_true("mu.par" %in% result$Subject)
+  
+  # Check Param is correctly decoded
+  expect_true(all(mu_rows$Param %in% c("y0", "y1", "t1", "alpha", "shape")))
+})
+
+test_that("unpack_jags correctly unpacks prec.par parameters", {
+  test_data <- tibble::tibble(
+    Iteration = rep(1:3, 2),
+    Chain = rep(1, 6),
+    Parameter = rep(c("prec.par[1,1,2]", "prec.par[1,2,3]"), each = 3),
+    value = rnorm(6)
+  )
+  
+  result <- serodynamics:::unpack_jags(test_data)
+  
+  # Check that prec.par rows are present
+  prec_rows <- result[result$Subject == "prec.par", ]
+  expect_gt(nrow(prec_rows), 0)
+  
+  # Check Subject is correctly set
+  expect_true("prec.par" %in% result$Subject)
+  
+  # Check Param contains pairs (e.g., "y0, y1")
+  expect_true(any(grepl(",", prec_rows$Param)))
+})
+
+test_that("unpack_jags correctly unpacks prec.logy parameters", {
+  test_data <- tibble::tibble(
+    Iteration = rep(1:3, 2),
+    Chain = rep(1, 6),
+    Parameter = rep(c("prec.logy[1]", "prec.logy[2]"), each = 3),
+    value = rnorm(6)
+  )
+  
+  result <- serodynamics:::unpack_jags(test_data)
+  
+  # Check that prec.logy rows are present
+  logy_rows <- result[result$Subject == "prec.logy", ]
+  expect_gt(nrow(logy_rows), 0)
+  
+  # Check Subject is correctly set
+  expect_true("prec.logy" %in% result$Subject)
+})
+
+test_that("unpack_jags correctly unpacks individual-level parameters", {
+  test_data <- tibble::tibble(
+    Iteration = rep(1:3, 2),
+    Chain = rep(1, 6),
+    Parameter = rep(c("y0[1,1]", "alpha[2,1]"), each = 3),
+    value = rnorm(6)
+  )
+  
+  result <- serodynamics:::unpack_jags(test_data)
+  
+  # Check individual-level parameters are unpacked
+  expect_true("1" %in% result$Subject | "2" %in% result$Subject)
+  expect_true("y0" %in% result$Param | "alpha" %in% result$Param)
+})
+
+test_that("prep_popparams filters to population parameters only", {
+  test_data <- tibble::tibble(
+    Iteration = 1:10,
+    Chain = rep(1, 10),
+    Parameter = rep(c("y0", "mu.par"), each = 5),
+    Subject = rep(c("1", "mu.par"), each = 5),
+    value = rnorm(10)
+  )
+  
+  result <- serodynamics:::prep_popparams(test_data)
+  
+  # Should only have population parameters
+  expect_true(all(result$Population_Parameter %in% 
+                    c("mu.par", "prec.par", "prec.logy")))
+  
+  # Should rename Subject to Population_Parameter
+  expect_true("Population_Parameter" %in% names(result))
+  expect_false("Subject" %in% names(result))
+  
+  # Should have 5 rows (only mu.par rows)
+  expect_equal(nrow(result), 5)
+})
+
+test_that("ex_popparams excludes population parameters", {
+  test_data <- tibble::tibble(
+    Iteration = 1:10,
+    Chain = rep(1, 10),
+    Parameter = rep(c("y0", "mu.par"), each = 5),
+    Subject = rep(c("1", "mu.par"), each = 5),
+    value = rnorm(10)
+  )
+  
+  result <- serodynamics:::ex_popparams(test_data)
+  
+  # Should not have any population parameters
+  expect_false(any(result$Subject %in% c("mu.par", "prec.par", "prec.logy")))
+  
+  # Should have 5 rows (only y0 rows)
+  expect_equal(nrow(result), 5)
+  
+  # Should keep original structure
+  expect_equal(names(result), names(test_data))
+})
+
+test_that("prep_popparams and ex_popparams are complementary", {
+  test_data <- tibble::tibble(
+    Iteration = rep(1:5, 4),
+    Chain = rep(1, 20),
+    Parameter = rep(c("y0", "mu.par", "prec.par", "prec.logy"), each = 5),
+    Subject = rep(c("1", "mu.par", "prec.par", "prec.logy"), each = 5),
+    value = rnorm(20)
+  )
+  
+  pop_params <- serodynamics:::prep_popparams(test_data)
+  individual_params <- serodynamics:::ex_popparams(test_data)
+  
+  # Together they should account for all rows
+  expect_equal(nrow(pop_params) + nrow(individual_params), nrow(test_data))
+  
+  # No overlap in subjects
+  # Note: pop_params has Population_Parameter column, not Subject
+  expect_false(any(individual_params$Subject %in% 
+                     c("mu.par", "prec.par", "prec.logy")))
+})
diff --git a/tests/testthat/test-run_mod.R b/tests/testthat/test-run_mod.R
index b4eef502..89186940 100644
--- a/tests/testthat/test-run_mod.R
+++ b/tests/testthat/test-run_mod.R
@@ -3,6 +3,7 @@ test_that(
   code = {
     testthat::announce_snapshot_file("sim-strat-curve-params.csv")
     testthat::announce_snapshot_file("sim-strat-fitted_residuals.csv")
+    testthat::announce_snapshot_file("popparam-summary-stats.csv")
     withr::local_seed(1)
     strat1 <- serocalculator::typhoid_curves_nostrat_100 |>
       sim_case_data(n = 100,
@@ -22,26 +23,49 @@ test_that(
       nmc = 10,
       niter = 10, # Number of iterations
       strat = "strat", # Variable to be stratified
+      with_pop_params = TRUE
     ) |>
       suppressWarnings()
     
-    
-    results |>
-      attributes() |>
-      rlist::list.remove(c("row.names", "fitted_residuals")) |>
-      expect_snapshot_value(style = "deparse")
-    
     results |>
       expect_snapshot_data(
         "sim-strat-curve-params",
         variant = darwin_variant()
       )
     
+    # Testing attributes
+    results |>
+      attributes() |>
+      names() |>
+      expect_setequal(c("names", "row.names", "class", "nChains", 
+                        "nParameters", "nIterations", "nBurnin", "nThin",
+                        "description", "population_params", "priors", 
+                        "fitted_residuals"))
+    
+    # Verify class appears immediately after names and row.names
+    expect_equal(
+      names(attributes(results))[1:3],
+      c("names", "row.names", "class")
+    )
+    
     attributes(results)$fitted_residuals |>
       expect_snapshot_data(
         "sim-strat-fitted_residuals",
         variant = darwin_variant()
       )
+
+    # Testing for population parameters
+    attributes(results)$population_params |>
+      dplyr::group_by(Parameter) |>
+      dplyr::summarise(
+        mean = mean(value),
+        sd = sd(value),
+        .groups = "drop"
+      ) |>
+      dplyr::arrange(Parameter) |>
+      expect_snapshot_data("popparam-summary-stats", 
+        variant = darwin_variant()    
+      )
     
   }
 )
@@ -62,7 +86,7 @@ test_that(
       nburn = 10, # Number of unrecorded samples before sampling begins
       nmc = 100,
       niter = 100, # Number of iterations
-      strat = "bldculres", # Variable to be stratified
+      strat = "bldculres" # Variable to be stratified
     ) |>
       suppressWarnings()
     
@@ -82,6 +106,7 @@ test_that(
         "strat-fitted_residuals",
         variant = darwin_variant()
       )
+    
   }
 )
 
@@ -101,7 +126,7 @@ test_that(
       nburn = 10, # Number of unrecorded samples before sampling begins
       nmc = 100,
       niter = 100, # Number of iterations
-      strat = NA, # Variable to be stratified
+      strat = NA # Variable to be stratified
     ) |>
       suppressWarnings()
     
@@ -144,12 +169,7 @@ test_that(
       with_post = TRUE
     ) |>
       suppressWarnings()
-    
-    results |>
-      attributes() |>
-      rlist::list.remove(c("row.names", "jags.post", "fitted_residuals")) |>
-      expect_snapshot_value(style = "serialize")
-    
+
     results |>
       expect_snapshot_data(
         "nostrat-curve-params-withpost",
@@ -183,15 +203,11 @@ test_that(
     ) |>
       suppressWarnings()
     
-    results |>
-      attributes() |>
-      rlist::list.remove(c("row.names", "fitted_residuals")) |>
-      expect_snapshot_value(style = "serialize")
-    
     results |>
       expect_snapshot_data(
         "nostrat-curve-params-specpriors",
         variant = darwin_variant()
       )
+    
   }
 )