diff --git a/R/prep_priors.R b/R/prep_priors.R index 952b2cea..ef994b98 100644 --- a/R/prep_priors.R +++ b/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,7 +32,7 @@ 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)) @@ -32,6 +40,16 @@ prep_priors <- function(max_antigens) { 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 diff --git a/inst/examples/run-mod-ka-example.R b/inst/examples/run-mod-ka-example.R new file mode 100644 index 00000000..20c00aa7 --- /dev/null +++ b/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)) +