Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 39 additions & 6 deletions .github/copilot-instructions.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,17 @@

## Repository Overview

**serodynamics** is an R package for modeling longitudinal antibody responses to infection. It implements Bayesian MCMC methods using JAGS (Just Another Gibbs Sampler) to estimate antibody dynamic curve parameters including baseline concentration, peak concentration, time to peak, shape parameter, and decay rate.
**serodynamics** is an R package for modeling longitudinal antibody responses to infection. It implements Bayesian MCMC methods to estimate antibody dynamic curve parameters including baseline concentration, peak concentration, time to peak, shape parameter, and decay rate.

The package supports two Bayesian modeling backends:
- **JAGS** (Just Another Gibbs Sampler) via `runjags` - the original implementation
- **Stan** via `cmdstanr` - a modern, efficient alternative (optional)

- **Type**: R package (statistical modeling)
- **Size**: ~121MB, ~209 files, ~76 R source files, ~1,664 lines of R code
- **Language**: R (>= 4.1.0)
- **Key Dependencies**: runjags, rjags, JAGS 4.3.1, serocalculator, ggmcmc, dplyr, ggplot2
- **Optional Dependencies**: cmdstanr (for Stan support)
- **Lifecycle**: Experimental (under active development)

## Critical Setup Requirements
Expand Down Expand Up @@ -213,9 +218,9 @@ devtools::check_man()

- **Windows**: Install Rtools from https://cran.r-project.org/bin/windows/Rtools/ (choose version matching your R version)

### JAGS Installation (REQUIRED)
### JAGS Installation (REQUIRED for JAGS models)

**ALWAYS install JAGS before attempting to build, test, or run this package.** The package will fail without it.
**Install JAGS if you plan to use `run_mod()` with JAGS models.** JAGS is required for the original JAGS-based modeling functions but is not needed if you only use Stan models via `run_mod_stan()`.

#### Installing JAGS in Docker (if using rocker/verse)

Expand Down Expand Up @@ -250,6 +255,25 @@ runjags::findJAGS()
runjags::testjags()
```

### Stan Installation (OPTIONAL for Stan models)

**Install cmdstanr and CmdStan if you plan to use `run_mod_stan()`.** Stan support is optional and provides a modern alternative to JAGS.

```r
# Install cmdstanr from r-universe
install.packages("cmdstanr",
repos = c("https://stan-dev.r-universe.dev",
getOption("repos")))

# Then install CmdStan
cmdstanr::install_cmdstan()

# Verify installation
cmdstanr::cmdstan_version()
```

See the [cmdstanr documentation](https://mc-stan.org/cmdstanr/) for more details and troubleshooting.

## Build and Development Workflow

### Initial Setup
Expand Down Expand Up @@ -374,10 +398,14 @@ Team members can trigger actions by commenting on PRs:

### Key Directories

- **R/**: Package source code (30 R files)
- **R/**: Package source code (30+ R files)
- `Run_Mod.R`: Main function to run JAGS Bayesian models
- `run_mod_stan.R`: Main function to run Stan Bayesian models (new)
- `as_case_data.R`: Convert data to case_data class
- `prep_data.r`, `prep_priors.R`: Data preparation for JAGS
- `prep_data.r`: Data preparation for JAGS
- `prep_data_stan.R`: Data preparation for Stan (new)
- `prep_priors.R`: Prior preparation for JAGS
- `prep_priors_stan.R`: Prior preparation for Stan (new)
- `sim_case_data.R`: Simulate case data for testing
- `post_summ.R`, `postprocess_jags_output.R`: Post-processing JAGS results
- `plot_*.R`: Diagnostic plotting functions (trace, density, Rhat, effective sample size)
Expand All @@ -397,8 +425,13 @@ Team members can trigger actions by commenting on PRs:
- **data-raw/**: Raw data processing scripts (not included in package build)

- **inst/**: Installed files
- `inst/extdata/`: JAGS model files (`model.jags`, `model.dobson.jags`), example CSV data
- `inst/extdata/`: Model files
- `model.jags`, `model.dobson.jags`: JAGS model files
- `model.stan`, `model.dobson.stan`: Stan model files (new)
- Example CSV data files
- `inst/examples/`: Example R scripts for documentation
- JAGS examples: `run_mod-examples.R`, `examples-prep_priors.R`
- Stan examples: `run_mod_stan-examples.R`, `examples-prep_priors_stan.R`, `examples-prep_data_stan.R` (new)
- `inst/WORDLIST`: Custom spelling dictionary

- **vignettes/**: Package vignettes
Expand Down
15 changes: 14 additions & 1 deletion .github/workflows/copilot-setup-steps.yml
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,19 @@ jobs:
runjags::testjags()
shell: Rscript {0}

# Install cmdstanr for Stan modeling support (optional)
# See .github/copilot-instructions.md "Stan Installation" section
- name: Install cmdstanr
run: |
install.packages("cmdstanr",
repos = c("https://stan-dev.r-universe.dev",
"https://cloud.r-project.org"),
verbose = TRUE)
cmdstanr::install_cmdstan(cores = 2)
cat("cmdstanr version:", as.character(packageVersion("cmdstanr")), "\n")
cat("CmdStan version:", cmdstanr::cmdstan_version(), "\n")
shell: Rscript {0}

# Verify JAGS installation
# See .github/copilot-instructions.md "JAGS Installation" section for details
- name: Verify JAGS installation
Expand Down Expand Up @@ -144,7 +157,7 @@ jobs:

# Display key installed packages
cat("Key installed packages:\n")
key_packages <- c("devtools", "rjags", "runjags", "rcmdcheck", "lintr", "spelling", "testthat")
key_packages <- c("devtools", "rjags", "runjags", "cmdstanr", "rcmdcheck", "lintr", "spelling", "testthat")
for (pkg in key_packages) {
if (requireNamespace(pkg, quietly = TRUE)) {
cat(" -", pkg, ":", as.character(packageVersion(pkg)), "\n")
Expand Down
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -28,3 +28,4 @@ docs/

**/*.quarto_ipynb
README.html
*.knit.*
4 changes: 3 additions & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -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."),
Expand Down Expand Up @@ -33,6 +33,7 @@ Imports:
tidyselect,
utils
Suggests:
cmdstanr,
Hmisc,
knitr,
readr,
Expand Down Expand Up @@ -63,6 +64,7 @@ Config/needs/test:
rlist
Config/needs/check: commonmark, xml2
Remotes:
stan-dev/cmdstanr,
ucd-serg/serocalculator
Depends:
R (>= 4.1.0)
Expand Down
3 changes: 3 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,11 @@ export(plot_predicted_curve)
export(post_summ)
export(postprocess_jags_output)
export(prep_data)
export(prep_data_stan)
export(prep_priors)
export(prep_priors_stan)
export(run_mod)
export(run_mod_stan)
export(serodynamics_example)
export(sim_case_data)
export(sim_n_obs)
Expand Down
9 changes: 9 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,14 @@
# serodynamics (development version)

## New features

* Added Stan support as an alternative to JAGS for Bayesian modeling
* New `run_mod_stan()` function for fitting models with Stan/cmdstanr
* New `prep_data_stan()` function to prepare data in Stan format
* New `prep_priors_stan()` function to prepare priors for Stan models
* New Stan model files: `inst/extdata/model.stan` and `inst/extdata/model.dobson.stan`
* Stan support is optional (cmdstanr in Suggests) and can be used alongside JAGS

* Added dev container configuration for persistent, cached development environment
that includes R, JAGS, and all dependencies preinstalled, making Copilot
Workspace sessions much faster.
Expand Down
76 changes: 76 additions & 0 deletions R/prep_data_stan.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
#' Prepare data for Stan
#'
#' @param dataframe a [data.frame] containing case data
#' @param biomarker_column [character] string indicating
#' which column contains antigen-isotype names
#' @param verbose whether to produce verbose messaging
#' @param add_newperson whether to add an extra record with missing data
#'
#' @returns a `prepped_stan_data` object (a [list] with Stan-formatted data)
#' @export
#'
#' @examples
#' set.seed(1)
#' raw_data <-
#' serocalculator::typhoid_curves_nostrat_100 |>
#' sim_case_data(n = 5)
#' prepped_data <- prep_data_stan(raw_data)
prep_data_stan <- function(
dataframe,
biomarker_column = get_biomarker_names_var(dataframe),
verbose = FALSE,
add_newperson = TRUE) {

# First use existing prep_data function to get the base structure
jags_data <- prep_data(
dataframe = dataframe,
biomarker_column = biomarker_column,
verbose = verbose,
add_newperson = add_newperson
)

# Convert to Stan format
# Stan requires explicit max dimensions
max_nsmpl <- max(jags_data$nsmpl)

# Create padded arrays (Stan doesn't handle ragged arrays like JAGS)
# We need to pad smpl.t and logy to max_nsmpl
nsubj <- jags_data$nsubj
n_antigen_isos <- jags_data$n_antigen_isos

# Initialize with zeros (will be ignored in model for obs > nsmpl[subj])
smpl_t_padded <- array(0, dim = c(nsubj, max_nsmpl))
logy_padded <- array(0, dim = c(nsubj, max_nsmpl, n_antigen_isos))

# Fill in actual data
for (subj in 1:nsubj) {
n_obs <- jags_data$nsmpl[subj]
if (n_obs > 0) {
smpl_t_padded[subj, 1:n_obs] <- jags_data$smpl.t[subj, 1:n_obs]
for (k in 1:n_antigen_isos) {
logy_padded[subj, 1:n_obs, k] <- jags_data$logy[subj, 1:n_obs, k]
}
}
}

stan_data <- list(
nsubj = nsubj,
n_antigen_isos = n_antigen_isos,
n_params = 5, # y0, y1, t1, alpha, shape
nsmpl = as.integer(jags_data$nsmpl),
max_nsmpl = as.integer(max_nsmpl),
smpl_t = smpl_t_padded,
logy = logy_padded
)

# Add attributes from JAGS data
stan_data <- stan_data |>
structure(
class = c("prepped_stan_data", "list"),
antigens = attributes(jags_data)$antigens,
n_antigens = attributes(jags_data)$n_antigens,
ids = attributes(jags_data)$ids
)

return(stan_data)
}
80 changes: 80 additions & 0 deletions R/prep_priors_stan.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
#' @title Prepare priors for Stan
#' @description
#' Takes multiple [vector] inputs to allow for modifiable priors for Stan
#' models. Converts JAGS precision-based priors to Stan covariance-based
#' priors.
#'
#' @inheritParams prep_priors
#'
#' @returns A "curve_params_priors_stan" object
#' (a subclass of [list] with the inputs to `prep_priors_stan()` attached
#' as [attributes] entry named `"used_priors"`), containing Stan-formatted
#' priors.
#' @export
#' @example inst/examples/examples-prep_priors_stan.R

prep_priors_stan <- function(
max_antigens,
mu_hyp_param = c(1.0, 7.0, 1.0, -4.0, -1.0),
prec_hyp_param = c(1.0, 0.00001, 1.0, 0.001, 1.0),
omega_param = c(1.0, 50.0, 1.0, 10.0, 1.0),
wishdf_param = 20,
prec_logy_hyp_param = c(4.0, 1.0)) {

# Input validation (same as prep_priors)
if (length(mu_hyp_param) != 5) {
cli::cli_abort("Need to specify 5 priors for {.arg mu_hyp_param}")
}
if (length(prec_hyp_param) != 5) {
cli::cli_abort("Need to specify 5 priors for {.arg prec_hyp_param}")
}
if (length(omega_param) != 5) {
cli::cli_abort("Need to specify 5 priors for {.arg omega_param}")
}
if (length(wishdf_param) != 1) {
cli::cli_abort("Need to specify 1 prior for {.arg wishdf_param}")
}
if (length(prec_logy_hyp_param) != 2) {
cli::cli_abort("Need to specify 2 priors for {.arg prec_logy_hyp_param}")
}

# Model parameters
n_params <- 5
mu_hyp <- array(NA, dim = c(max_antigens, n_params))
prec_hyp <- array(NA, dim = c(max_antigens, n_params, n_params))
omega <- array(NA, dim = c(max_antigens, n_params, n_params))
wishdf <- rep(NA, max_antigens)
prec_logy_hyp <- array(NA, dim = c(max_antigens, 2))

# Fill parameter arrays (same structure as JAGS)
for (k in 1:max_antigens) {
mu_hyp[k, ] <- mu_hyp_param
prec_hyp[k, , ] <- diag(prec_hyp_param)
omega[k, , ] <- diag(omega_param)
wishdf[k] <- wishdf_param
prec_logy_hyp[k, ] <- prec_logy_hyp_param
}

# Return results as a list (Stan model will handle conversion)
prepped_priors <- list(
"n_params" = n_params,
"mu_hyp" = mu_hyp,
"prec_hyp" = prec_hyp,
"omega" = omega,
"wishdf" = wishdf,
"prec_logy_hyp" = prec_logy_hyp
) |>
structure(class = c("curve_params_priors_stan", "list"))

# Add used priors as attributes
prepped_priors <- prepped_priors |>
structure("used_priors" = list(
mu_hyp_param = mu_hyp_param,
prec_hyp_param = prec_hyp_param,
omega_param = omega_param,
wishdf_param = wishdf_param,
prec_logy_hyp_param = prec_logy_hyp_param
))

return(prepped_priors)
}
Loading