Explaining priors

R/prep_priors.R

@@ -1,15 +1,23 @@
 #' @title Prepare priors
+#'The function is preparing a set of hyperparameters for a Bayesian model that characterizes antibody (or biomarker) responses. 
+#'It assumes that each antigen-isotype is modeled by 5 parameters, which are transformed by taking their logarithm 
+#'(with the fifth parameter adjusted as shape-1). These 5 parameters are named in the code as:
+#'1. y0
+#'2. y1,
+#'3. t1,
+#'4. alpha, and
+#'5. shape 
 #'
 #' @param max_antigens [integer]: how many antigen-isotypes will be modeled
 #'
 #' @returns a [list] with elements:
-#' "n_params": how many parameters??
-#' - "mu.hyp": ??
-#' - "prec.hyp": ??
-#' - "omega" : ??
-#' - "wishdf": Wishart distribution degrees of freedom
-#' - "prec.logy.hyp": array of hyper-parameters for the precision
-#' (inverse variance) of the biomarkers, on log-scale
+#' n_params: Number of model parameters (log-transformed y0, y1, t1, alpha, shape-1).
+#' mu.hyp: Hyperparameter means for the log-transformed model parameters.
+#' prec.hyp: Diagonal precision matrix for the normal priors, setting the inverse variance for each parameter.
+#' omega: Diagonal scale matrix used with the Wishart prior to shape the covariance structure.
+#' wishdf: Degrees of freedom for the Wishart prior (set to 20).
+#' prec.logy.hyp: Hyperparameters (e.g., shape and rate) for the precision (inverse variance) of the log-scale biomarkers.
+
 #' @export
 #'
 #' @examples
@@ -24,14 +32,24 @@ prep_priors <- function(max_antigens) {
   prec_logy_hyp <- array(NA, dim = c(max_antigens, 2))
 
   # Fill parameter arrays
-  # the parameters are log(c(y0,  y1,    t1,  alpha, shape-1))
+ # the parameters are log(c(y0,  y1,    t1,  alpha, shape-1))
   for (k.test in 1:max_antigens) {
     mu_hyp[k.test, ] <- c(1.0, 7.0, 1.0, -4.0, -1.0)
     prec_hyp[k.test, , ] <- diag(c(1.0, 0.00001, 1.0, 0.001, 1.0))
     omega[k.test, , ] <- diag(c(1.0, 50.0, 1.0, 10.0, 1.0))
     wishdf[k.test] <- 20
     prec_logy_hyp[k.test, ] <- c(4.0, 1.0)
   }
+
+
+  # uninformative testing priors
+  # for (k.test in 1:max_antigens) {
+  #   mu_hyp[k.test, ] <- c(1.0, 1.0, 1.0, -4.0, -1.0)
+  #   prec_hyp[k.test, , ] <- diag(c(0.00001, 1e-8, 0.00001, 0.00001, 0.00001))
+  #   omega[k.test, , ] <- diag(c(1.0, 50.0, 1.0, 10.0, 1.0))
+  #   wishdf[k.test] <- 20
+  #   prec_logy_hyp[k.test, ] <- c(4.0, 1.0)
+  # }
 
   # test for change
 

inst/examples/run-mod-ka-example.R

@@ -0,0 +1,49 @@
+
+options(scipen = 999)
+library(runjags)
+
+    dataset <- serodynamics_example(
+      "SEES_Case_Nepal_ForSeroKinetics_02-13-2025.csv"
+    ) |>
+      readr::read_csv() |>
+      dplyr::mutate(
+        .by = person_id,
+        visit_num = dplyr::row_number()
+      ) |>
+      as_case_data(
+        id_var = "person_id",
+        biomarker_var = "antigen_iso",
+        value_var = "result",
+        time_in_days = "dayssincefeveronset"
+      )
+
+    results <- run_mod(
+      data = dataset, # The data set input
+      file_mod = serodynamics_example("model.jags"),
+      nchain = 2, # Number of mcmc chains to run
+      nadapt = 100, # Number of adaptations to run
+      nburn = 100, # Number of unrecorded samples before sampling begins
+      nmc = 100,
+      niter = 1000, # Number of iterations
+      strat = NA # Variable to be stratified
+    ) |>
+      suppressWarnings() |>
+      magrittr::use_series("curve_params")
+
+
+test.1 <- results %>% group_by(Iso_type, Parameter_sub) %>% summarise(mean = mean(value), sd = sd(value))
+
+wide_predpar_df <- results %>%
+  select(-Parameter) %>%
+  rename(parameter = Parameter_sub,
+         antigen_iso = Iso_type) %>%
+  pivot_wider(names_from = "parameter", values_from="value") %>%
+  rowwise() %>%
+  #mutate(y1 = y0+y1) %>%
+  droplevels() %>%
+  ungroup() %>%
+  rename(r = shape)
+
+
+autoplot(serocalculator::as_curve_params(wide_predpar_df))
+