An R6 class for analyzing correlates of protection using Bayesian methods. Fits a four-parameter logistic model to antibody titre and infection outcome data.
This class implements a Bayesian logistic model for analyzing correlates of protection. It provides methods for model fitting, prediction, performance evaluation, and visualization.
SeroCOP R6 Class for Correlates of Protection Analysis
This class provides a complete workflow for correlates of protection analysis:
Data input and validation
Bayesian model fitting using brms (which uses Stan backend)
Model diagnostics and validation
Performance metrics (ROC AUC, Brier score, LOO-CV)
Visualization
titreNumeric vector of antibody titres (log scale recommended)
infectedBinary vector of infection status (0/1)
groupOptional factor vector for hierarchical modeling
fitStan fit object (after fitting)
looLOO-CV object (after fitting)
loo_by_groupList of per-group LOO contributions (hierarchical models only)
priorsList of prior distributions for model parameters
weightsOptional numeric vector of observation weights
new()Create a new SeroCOP object
SeroCOP$new(titre, infected, group = NULL, weights = NULL)titreNumeric vector of antibody titres
infectedBinary vector (0/1) of infection outcomes
groupOptional factor or character vector for grouping (e.g., age groups)
weightsOptional numeric vector of non-negative observation weights.
When supplied, a weighted Bernoulli likelihood is used. If NULL
(the default) all observations are given equal weight.
definePrior()Define custom prior distributions for model parameters
SeroCOP$definePrior(
floor_alpha = NULL,
floor_beta = NULL,
ceiling_alpha = NULL,
ceiling_beta = NULL,
ec50_mean = NULL,
ec50_sd = NULL,
slope_mean = NULL,
slope_sd = NULL
)floor_alphaAlpha parameter for floor beta prior (default: 1)
floor_betaBeta parameter for floor beta prior (default: 9)
ceiling_alphaAlpha parameter for ceiling beta prior (default: 9)
ceiling_betaBeta parameter for ceiling beta prior (default: 1)
ec50_meanMean for ec50 normal prior (default: midpoint of titre range)
ec50_sdSD for ec50 normal prior (default: titre range / 4)
slope_meanMean for slope normal prior (default: 0)
slope_sdSD for slope normal prior (default: 2)
\dontrun{
model <- SeroCOP$new(titre = titre, infected = infected)
# Use default priors centered on data
model$definePrior()
# Custom priors
model$definePrior(
floor_alpha = 2, floor_beta = 18, # Stronger prior for low floor
ceiling_alpha = 18, ceiling_beta = 2, # Stronger prior for high ceiling
ec50_mean = 2.0, ec50_sd = 1.0, # Custom ec50 prior
slope_mean = 0, slope_sd = 1 # More conservative slope
)
}
fit_model()Fit the Bayesian logistic model
SeroCOP$fit_model(
chains = 4,
iter = 2000,
warmup = 1000,
thin = 1,
ceiling_group = FALSE,
...
)chainsNumber of MCMC chains (default: 4)
iterNumber of iterations per chain (default: 2000)
warmupNumber of warmup iterations (default: iter/2)
thinThinning interval (default: 1)
ceiling_groupLogical; when TRUE and a group variable is present,
the ceiling parameter also receives group-level random effects (default: FALSE).
...Additional arguments passed to brms::brm
predict()Get posterior predictions for infection probability
newdataOptional new titre values for prediction
groupOptional group level (character) for group-specific predictions.
If NULL (default), predictions are marginalised over all groups
(population-average curve). Supply a valid group level string (e.g.
"5-11 years") to obtain the group-specific curve. Ignored when
the model was not fitted with a group variable.
predict_protection()Extract probability of protection from the fitted model
newdataOptional vector of new titre values for prediction
groupOptional group level (character) for group-specific predictions.
If NULL (default), predictions are marginalised over all groups.
See predict() for details.
summary()Get summary statistics for model parameters
plot_group_curves()Plot group-specific correlates of risk curves with a population-weighted marginalised curve overlaid on each panel.
plot_cop_curves()Plot group-specific and marginalised correlates of protection (COP) curves. The COP at titre \(t\) for group \(g\) is defined as \(1 - P(\text{infection} \mid t, g) / \text{ceiling}_g\), i.e.\ the probability of protection conditional on exposure. The marginalised curve is the weighted average of group-specific COP curves.
extract_cop()Extract the correlate of protection conditional on exposure.
if (FALSE) { # \dontrun{
# Create synthetic data
set.seed(123)
n <- 200
titre <- rnorm(n, mean = 2, sd = 1.5)
prob <- 0.05 + 0.9 / (1 + exp(2 * (titre - 1.5)))
infected <- rbinom(n, 1, prob)
# Initialize and fit
model <- SeroCOP$new(titre = titre, infected = infected)
model$fit(chains = 4, iter = 2000)
# Get metrics
model$get_metrics()
# Plot results
model$plot_curve()
model$plot_roc()
} # }
## ------------------------------------------------
## Method `SeroCOP$definePrior`
## ------------------------------------------------
if (FALSE) { # \dontrun{
model <- SeroCOP$new(titre = titre, infected = infected)
# Use default priors centered on data
model$definePrior()
# Custom priors
model$definePrior(
floor_alpha = 2, floor_beta = 18, # Stronger prior for low floor
ceiling_alpha = 18, ceiling_beta = 2, # Stronger prior for high ceiling
ec50_mean = 2.0, ec50_sd = 1.0, # Custom ec50 prior
slope_mean = 0, slope_sd = 1 # More conservative slope
)
} # }