Title: | Estimating Infection Rates from Serological Data |
---|---|
Description: | Translates antibody levels measured in cross-sectional population samples into estimates of the frequency with which seroconversions (infections) occur in the sampled populations. Replaces the previous `seroincidence` package. |
Authors: | Peter Teunis [aut, cph] (Author of the method and original code.), Kristina Lai [aut, cre], Chris Orwa [aut], Kristen Aiemjoy [aut], Douglas Ezra Morrison [aut] |
Maintainer: | Kristina Lai <[email protected]> |
License: | GPL-3 |
Version: | 1.3.0.9029 |
Built: | 2025-02-18 04:27:47 UTC |
Source: | https://github.com/ucd-serg/serocalculator |
Load antibody decay curve parameter
as_curve_params(data, antigen_isos = NULL)
as_curve_params(data, antigen_isos = NULL)
data |
|
antigen_isos |
a |
a curve_data
object
(a tibble::tbl_df with extra attribute antigen_isos
)
library(magrittr) curve_data <- serocalculator_example("example_curve_params.csv") %>% read.csv() %>% as_curve_params() print(curve_data)
library(magrittr) curve_data <- serocalculator_example("example_curve_params.csv") %>% read.csv() %>% as_curve_params() print(curve_data)
Load noise parameters
as_noise_params(data, antigen_isos = NULL)
as_noise_params(data, antigen_isos = NULL)
data |
|
antigen_isos |
|
a noise_params
object (a tibble::tbl_df with
extra attribute antigen_isos
)
library(magrittr) noise_data <- serocalculator_example("example_noise_params.csv") %>% read.csv() %>% as_noise_params() print(noise_data)
library(magrittr) noise_data <- serocalculator_example("example_noise_params.csv") %>% read.csv() %>% as_noise_params() print(noise_data)
Load a cross-sectional antibody survey data set
as_pop_data( data, antigen_isos = NULL, age = "Age", value = "result", id = "index_id", standardize = TRUE )
as_pop_data( data, antigen_isos = NULL, age = "Age", value = "result", id = "index_id", standardize = TRUE )
data |
|
antigen_isos |
|
age |
a |
value |
a |
id |
a |
standardize |
a |
a pop_data
object (a tibble::tbl_df
with extra attribute antigen_isos
)
library(magrittr) xs_data <- serocalculator_example("example_pop_data.csv") |> read.csv() |> as_pop_data() print(xs_data)
library(magrittr) xs_data <- serocalculator_example("example_pop_data.csv") |> read.csv() |> as_pop_data() print(xs_data)
graph antibody decay curves by antigen isotype
## S3 method for class 'curve_params' autoplot( object, antigen_isos = unique(object$antigen_iso), ncol = min(3, length(antigen_isos)), ... )
## S3 method for class 'curve_params' autoplot( object, antigen_isos = unique(object$antigen_iso), ncol = min(3, length(antigen_isos)), ... )
object |
a |
antigen_isos |
antigen isotypes to analyze (can subset |
ncol |
how many columns of subfigures to use in panel plot |
... |
Arguments passed on to
|
rows_to_graph
If you directly specify rows_to_graph
when calling this function,
the row numbers are enumerated separately for each antigen isotype;
in other words, for the purposes of this argument,
row numbers start over at 1 for each antigen isotype.
There is currently no way to specify different row numbers for different antigen isotypes;
if you want to do that,
you will could call plot_curve_params_one_ab()
directly for each antigen isotype
and combine the resulting panels yourself.
Or you could subset curve_params
manually,
before passing it to this function,
and set the n_curves
argument to Inf
.
a ggplot2::ggplot()
object
library(dplyr) library(ggplot2) library(magrittr) curve <- serocalculator_example("example_curve_params.csv") %>% read.csv() %>% as_curve_params() %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% autoplot() curve
library(dplyr) library(ggplot2) library(magrittr) curve <- serocalculator_example("example_curve_params.csv") %>% read.csv() %>% as_curve_params() %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% autoplot() curve
autoplot()
method for pop_data
objects
## S3 method for class 'pop_data' autoplot(object, log = FALSE, type = "density", strata = NULL, ...)
## S3 method for class 'pop_data' autoplot(object, log = FALSE, type = "density", strata = NULL, ...)
object |
A |
log |
whether to show antibody responses on logarithmic scale |
type |
an option to choose type of chart:
the current options are |
strata |
the name of a variable in |
... |
unused |
a ggplot2::ggplot object
library(dplyr) library(ggplot2) library(magrittr) xs_data <- serocalculator_example("example_pop_data.csv") %>% read.csv() %>% as_pop_data() xs_data %>% autoplot(strata = "catchment", type = "density") xs_data %>% autoplot(strata = "catchment", type = "age-scatter")
library(dplyr) library(ggplot2) library(magrittr) xs_data <- serocalculator_example("example_pop_data.csv") %>% read.csv() %>% as_pop_data() xs_data %>% autoplot(strata = "catchment", type = "density") xs_data %>% autoplot(strata = "catchment", type = "age-scatter")
Plot the log-likelihood curve for the incidence rate estimate
## S3 method for class 'seroincidence' autoplot(object, log_x = FALSE, ...)
## S3 method for class 'seroincidence' autoplot(object, log_x = FALSE, ...)
object |
a |
log_x |
should the x-axis be on a logarithmic scale ( |
... |
unused |
library(dplyr) library(ggplot2) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est1 <- est.incidence( pop_data = xs_data, curve_param = curve, noise_param = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), build_graph = TRUE ) # Plot the log-likelihood curve autoplot(est1)
library(dplyr) library(ggplot2) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est1 <- est.incidence( pop_data = xs_data, curve_param = curve, noise_param = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), build_graph = TRUE ) # Plot the log-likelihood curve autoplot(est1)
seroincidence.by
log-likelihoodsPlots log-likelihood curves by stratum, for seroincidence.by
objects
## S3 method for class 'seroincidence.by' autoplot(object, ncol = min(3, length(object)), ...)
## S3 method for class 'seroincidence.by' autoplot(object, ncol = min(3, length(object)), ...)
object |
a '"seroincidence.by"' object (from |
ncol |
number of columns to use for panel of plots |
... |
Arguments passed on to
|
an object of class "ggarrange"
, which is a ggplot2::ggplot()
or a list()
of ggplot2::ggplot()
s.
library(dplyr) library(ggplot2) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data, curve_params = curve, curve_strata_varnames= NULL, noise_strata_varnames = NULL, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), #num_cores = 8, #Allow for parallel processing to decrease run time build_graph = TRUE ) # Plot the log-likelihood curve autoplot(est2)
library(dplyr) library(ggplot2) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data, curve_params = curve, curve_strata_varnames= NULL, noise_strata_varnames = NULL, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), #num_cores = 8, #Allow for parallel processing to decrease run time build_graph = TRUE ) # Plot the log-likelihood curve autoplot(est2)
summary.seroincidence.by
objectsPlot method for summary.seroincidence.by
objects
## S3 method for class 'summary.seroincidence.by' autoplot( object, xvar, alpha = 0.7, shape = 1, dodge_width = 0.001, CIs = FALSE, ... )
## S3 method for class 'summary.seroincidence.by' autoplot( object, xvar, alpha = 0.7, shape = 1, dodge_width = 0.001, CIs = FALSE, ... )
object |
a |
xvar |
the name of a stratifying variable in |
alpha |
transparency for the points in the graph (1 = no transparency, 0 = fully transparent) |
shape |
shape argument for |
dodge_width |
width for jitter |
CIs |
logical, if |
... |
unused |
a ggplot2::ggplot()
object
library(dplyr) library(ggplot2) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data, curve_params = curve, noise_params = noise, curve_strata_varnames= NULL, noise_strata_varnames = NULL, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), num_cores = 2 # Allow for parallel processing to decrease run time ) est2sum <- summary(est2) autoplot(est2sum, "catchment")
library(dplyr) library(ggplot2) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data, curve_params = curve, noise_params = noise, curve_strata_varnames= NULL, noise_strata_varnames = NULL, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), num_cores = 2 # Allow for parallel processing to decrease run time ) est2sum <- summary(est2) autoplot(est2sum, "catchment")
Check the formatting of a cross-sectional antibody survey dataset.
check_pop_data(pop_data, verbose = FALSE)
check_pop_data(pop_data, verbose = FALSE)
pop_data |
dataset to check |
verbose |
whether to print an "OK" message when all checks pass |
NULL (invisibly)
library(magrittr) xs_data <- serocalculator_example("example_pop_data.csv") %>% read.csv() %>% as_pop_data() check_pop_data(xs_data, verbose = TRUE)
library(magrittr) xs_data <- serocalculator_example("example_pop_data.csv") %>% read.csv() %>% as_pop_data() check_pop_data(xs_data, verbose = TRUE)
This function models seroincidence using maximum likelihood estimation; that is, it finds the value of the seroincidence parameter which maximizes the likelihood (i.e., joint probability) of the data.
est.incidence( pop_data, curve_params, noise_params, antigen_isos = get_biomarker_names(pop_data), lambda_start = 0.1, stepmin = 1e-08, stepmax = 3, verbose = FALSE, build_graph = FALSE, print_graph = build_graph & verbose, ... )
est.incidence( pop_data, curve_params, noise_params, antigen_isos = get_biomarker_names(pop_data), lambda_start = 0.1, stepmin = 1e-08, stepmax = 3, verbose = FALSE, build_graph = FALSE, print_graph = build_graph & verbose, ... )
pop_data |
a data.frame with cross-sectional serology data per antibody and age, and additional columns |
curve_params |
a
|
noise_params |
a
|
antigen_isos |
Character vector with one or more antibody names. Values must match |
lambda_start |
starting guess for incidence rate, in years/event. |
stepmin |
A positive scalar providing the minimum allowable relative step length. |
stepmax |
a positive scalar which gives the maximum allowable
scaled step length. |
verbose |
logical: if TRUE, print verbose log information to console |
build_graph |
whether to graph the log-likelihood function across a range of incidence rates (lambda values) |
print_graph |
whether to display the log-likelihood curve graph in the course of running |
... |
Arguments passed on to
|
a "seroincidence"
object, which is a stats::nlm()
fit object with extra meta-data attributes lambda_start
, antigen_isos
, and ll_graph
library(dplyr) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est1 <- est.incidence( pop_data = xs_data, curve_params = curve, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), ) summary(est1)
library(dplyr) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est1 <- est.incidence( pop_data = xs_data, curve_params = curve, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), ) summary(est1)
Function to estimate seroincidences based on cross-sectional serology data and longitudinal response model.
est.incidence.by( pop_data, curve_params, noise_params, strata, curve_strata_varnames = strata, noise_strata_varnames = strata, antigen_isos = pop_data %>% pull("antigen_iso") %>% unique(), lambda_start = 0.1, build_graph = FALSE, num_cores = 1L, verbose = FALSE, print_graph = FALSE, ... )
est.incidence.by( pop_data, curve_params, noise_params, strata, curve_strata_varnames = strata, noise_strata_varnames = strata, antigen_isos = pop_data %>% pull("antigen_iso") %>% unique(), lambda_start = 0.1, build_graph = FALSE, num_cores = 1L, verbose = FALSE, print_graph = FALSE, ... )
pop_data |
a data.frame with cross-sectional serology data per
antibody and age, and additional columns corresponding to
each element of the |
curve_params |
a
|
noise_params |
a
|
strata |
a character vector of stratum-defining variables.
Values must be variable names in |
curve_strata_varnames |
A subset of |
noise_strata_varnames |
A subset of |
antigen_isos |
Character vector with one or more antibody names. Values must match |
lambda_start |
starting guess for incidence rate, in years/event. |
build_graph |
whether to graph the log-likelihood function across a range of incidence rates (lambda values) |
num_cores |
Number of processor cores to use for calculations when computing by strata. If set to more than 1 and package parallel is available, then the computations are executed in parallel. Default = 1L. |
verbose |
logical: if TRUE, print verbose log information to console |
print_graph |
whether to display the log-likelihood curve graph in the course of running |
... |
Arguments passed on to
|
If strata
is left empty, a warning will be produced,
recommending that you use est.incidence()
for unstratified analyses,
and then the data will be passed to est.incidence()
.
If for some reason you want to use est.incidence.by()
with no strata instead of calling est.incidence()
,
you may use NA
, NULL
, or ""
as the strata
argument to avoid that warning.
if strata
has meaningful inputs:
An object of class "seroincidence.by"
; i.e., a list of
"seroincidence"
objects from est.incidence()
, one for each stratum,
with some meta-data attributes.
if strata
is missing, NULL
, NA
, or ""
:
An object of class "seroincidence"
.
library(dplyr) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est2 <- est.incidence.by( strata = "catchment", pop_data = xs_data, curve_params = curve, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), # num_cores = 8 # Allow for parallel processing to decrease run time iterlim = 5 # limit iterations for the purpose of this example ) print(est2) summary(est2)
library(dplyr) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est2 <- est.incidence.by( strata = "catchment", pop_data = xs_data, curve_params = curve, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), # num_cores = 8 # Allow for parallel processing to decrease run time iterlim = 5 # limit iterations for the purpose of this example ) print(est2) summary(est2)
A subset of noise parameter estimates from the SEES study, for examples and testing, for Pakistan
example_noise_params_pk
example_noise_params_pk
example_noise_params_pk
A curve_params
object (from as_curve_params()
) with 4 rows and 7 columns:
which antigen and isotype are being measured (data is in long format)
Location for which the noise parameters were estimated
Lower limit of detection
Measurement noise, defined by a CV (coefficient of variation) as the ratio of the standard deviation to the mean for replicates. Note that the CV should ideally be measured across plates rather than within the same plate.
Biological noise: error from cross-reactivity to other antibodies. It is defined as the 95th percentile of the distribution of antibody responses to the antigen-isotype in a population with no exposure.
Upper limit of detection
Lab for which noise was estimated.
A subset of noise parameter estimates from the SEES study, for examples and testing.
example_noise_params_sees
example_noise_params_sees
example_noise_params_pk
A curve_params
object (from as_curve_params()
) with 4 rows and 7 columns:
which antigen and isotype are being measured (data is in long format)
Location for which the noise parameters were estimated
Lower limit of detection
Measurement noise, defined by a CV (coefficient of variation) as the ratio of the standard deviation to the mean for replicates. Note that the CV should ideally be measured across plates rather than within the same plate.
Biological noise: error from cross-reactivity to other antibodies. It is defined as the 95th percentile of the distribution of antibody responses to the antigen-isotype in a population with no exposure.
Upper limit of detection
Lab for which noise was estimated.
Extract biomarker levels
get_biomarker_levels(object, ...)
get_biomarker_levels(object, ...)
object |
a |
... |
unused |
the biomarker levels in object
sees_pop_data_100 |> get_biomarker_levels()
sees_pop_data_100 |> get_biomarker_levels()
Get biomarker variable name
get_biomarker_names_var(object, ...)
get_biomarker_names_var(object, ...)
object |
a |
... |
unused |
a character string identifying the biomarker names column in object
sees_pop_data_100 |> get_biomarker_names_var()
sees_pop_data_100 |> get_biomarker_names_var()
Get antibody measurement values
get_values(object, ...)
get_values(object, ...)
object |
a |
... |
unused |
a numeric vector of antibody measurement values
sees_pop_data_100 |> get_values()
sees_pop_data_100 |> get_values()
Extract antibody measurement values
get_values_var(object, ...)
get_values_var(object, ...)
object |
a |
... |
unused |
the name of the column in object
specified as containing
antibody abundance measurements
sees_pop_data_100 |> get_values_var()
sees_pop_data_100 |> get_values_var()
Graph log-likelihood of data
graph_loglik( pop_data, curve_params, noise_params, antigen_isos = pop_data %>% get_biomarker_levels(), x = 10^seq(-3, 0, by = 0.1), highlight_points = NULL, highlight_point_names = "highlight_points", log_x = FALSE, previous_plot = NULL, curve_label = paste(antigen_isos, collapse = " + "), ... )
graph_loglik( pop_data, curve_params, noise_params, antigen_isos = pop_data %>% get_biomarker_levels(), x = 10^seq(-3, 0, by = 0.1), highlight_points = NULL, highlight_point_names = "highlight_points", log_x = FALSE, previous_plot = NULL, curve_label = paste(antigen_isos, collapse = " + "), ... )
pop_data |
a |
curve_params |
a
|
noise_params |
a
|
antigen_isos |
Character vector listing one or more antigen isotypes.
Values must match |
x |
sequence of lambda values to graph |
highlight_points |
a possible highlighted value |
highlight_point_names |
labels for highlighted points |
log_x |
should the x-axis be on a logarithmic scale ( |
previous_plot |
if not NULL, the current data is added to the existing graph |
curve_label |
if not NULL, add a label for the curve |
... |
Arguments passed on to
|
library(dplyr) library(tibble) # Load cross-sectional data xs_data <- sees_pop_data_pk_100 # Load curve parameters and subset for the purposes of this example curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) # Load noise parameters cond <- tibble( antigen_iso = c("HlyE_IgG", "HlyE_IgA"), nu = c(0.5, 0.5), # Biologic noise (nu) eps = c(0, 0), # M noise (eps) y.low = c(1, 1), # Low cutoff (llod) y.high = c(5e6, 5e6)) # High cutoff (y.high) # Graph the log likelihood lik_HlyE_IgA <- # nolint: object_name_linter graph_loglik( pop_data = xs_data, curve_params = curve, noise_params = cond, antigen_isos = "HlyE_IgA", log_x = TRUE ) lik_HlyE_IgA # nolint: object_name_linter
library(dplyr) library(tibble) # Load cross-sectional data xs_data <- sees_pop_data_pk_100 # Load curve parameters and subset for the purposes of this example curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) # Load noise parameters cond <- tibble( antigen_iso = c("HlyE_IgG", "HlyE_IgA"), nu = c(0.5, 0.5), # Biologic noise (nu) eps = c(0, 0), # M noise (eps) y.low = c(1, 1), # Low cutoff (llod) y.high = c(5e6, 5e6)) # High cutoff (y.high) # Graph the log likelihood lik_HlyE_IgA <- # nolint: object_name_linter graph_loglik( pop_data = xs_data, curve_params = curve, noise_params = cond, antigen_isos = "HlyE_IgA", log_x = TRUE ) lik_HlyE_IgA # nolint: object_name_linter
Graph estimated antibody decay curve
graph.curve.params( curve_params, antigen_isos = unique(curve_params$antigen_iso), verbose = FALSE, show_quantiles = TRUE, show_all_curves = FALSE, alpha_samples = 0.3 )
graph.curve.params( curve_params, antigen_isos = unique(curve_params$antigen_iso), verbose = FALSE, show_quantiles = TRUE, show_all_curves = FALSE, alpha_samples = 0.3 )
curve_params |
a |
antigen_isos |
antigen isotypes |
verbose |
verbose output |
show_quantiles |
whether to show point-wise (over time) quantiles |
show_all_curves |
whether to show individual curves under quantiles |
alpha_samples |
|
a ggplot2::ggplot()
object
curve <- typhoid_curves_nostrat_100 |> dplyr::filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) plot1 <- graph.curve.params(curve) print(plot1) plot2 <- graph.curve.params(curve, show_all_curves = TRUE) show(plot2)
curve <- typhoid_curves_nostrat_100 |> dplyr::filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) plot1 <- graph.curve.params(curve) print(plot1) plot2 <- graph.curve.params(curve, show_all_curves = TRUE) show(plot2)
Load antibody decay curve parameter samples
load_curve_params(file_path, antigen_isos = NULL)
load_curve_params(file_path, antigen_isos = NULL)
file_path |
path to an RDS file containing MCMC samples of antibody decay curve parameters |
antigen_isos |
|
a curve_params
object (a tibble::tbl_df with extra attribute antigen_isos
)
curve <- load_curve_params(serocalculator_example("example_curve_params.rds")) print(curve)
curve <- load_curve_params(serocalculator_example("example_curve_params.rds")) print(curve)
Load noise parameters
load_noise_params(file_path, antigen_isos = NULL)
load_noise_params(file_path, antigen_isos = NULL)
file_path |
path to an RDS file containing biologic
and measurement noise of antibody decay curve parameters
|
antigen_isos |
|
a noise
object (a tibble::tbl_df
with extra attribute antigen_isos
)
noise <- load_noise_params(serocalculator_example("example_noise_params.rds")) print(noise)
noise <- load_noise_params(serocalculator_example("example_noise_params.rds")) print(noise)
Load a cross-sectional antibody survey data set
load_pop_data(file_path, ...)
load_pop_data(file_path, ...)
file_path |
path to an RDS file containing a cross-sectional antibody
survey data set, stored as a |
... |
Arguments passed on to
|
a pop_data
object (a tibble::tbl_df with extra attributes)
xs_data <- load_pop_data(serocalculator_example("example_pop_data.rds")) print(xs_data)
xs_data <- load_pop_data(serocalculator_example("example_pop_data.rds")) print(xs_data)
Calculates the log-likelihood of a set of cross-sectional antibody response
data, for a given incidence rate (lambda
) value.
log_likelihood( lambda, pop_data, curve_params, noise_params, antigen_isos = get_biomarker_levels(pop_data), verbose = FALSE, ... )
log_likelihood( lambda, pop_data, curve_params, noise_params, antigen_isos = get_biomarker_levels(pop_data), verbose = FALSE, ... )
lambda |
a numeric vector of incidence parameters, in events per person-year |
pop_data |
a |
curve_params |
a
|
noise_params |
a
|
antigen_isos |
Character vector listing one or more antigen isotypes.
Values must match |
verbose |
logical: if TRUE, print verbose log information to console |
... |
additional arguments passed to other functions (not currently used). |
the log-likelihood of the data with the current parameter values
library(dplyr) library(tibble) # Load cross-sectional data xs_data <- sees_pop_data_pk_100 # Load curve parameters and subset for the purposes of this example curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) # Load noise params cond <- tibble( antigen_iso = c("HlyE_IgG", "HlyE_IgA"), nu = c(0.5, 0.5), # Biologic noise (nu) eps = c(0, 0), # M noise (eps) y.low = c(1, 1), # low cutoff (llod) y.high = c(5e6, 5e6) ) # high cutoff (y.high) # Calculate log-likelihood ll_AG <- log_likelihood( pop_data = xs_data, curve_params = curve, noise_params = cond, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), lambda = 0.1 ) %>% print()
library(dplyr) library(tibble) # Load cross-sectional data xs_data <- sees_pop_data_pk_100 # Load curve parameters and subset for the purposes of this example curve <- typhoid_curves_nostrat_100 %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) # Load noise params cond <- tibble( antigen_iso = c("HlyE_IgG", "HlyE_IgA"), nu = c(0.5, 0.5), # Biologic noise (nu) eps = c(0, 0), # M noise (eps) y.low = c(1, 1), # low cutoff (llod) y.high = c(5e6, 5e6) ) # high cutoff (y.high) # Calculate log-likelihood ll_AG <- log_likelihood( pop_data = xs_data, curve_params = curve, noise_params = cond, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), lambda = 0.1 ) %>% print()
A subset of data from the SEES data, for examples and testing.
sees_pop_data_100
sees_pop_data_100
sees_pop_data_pk_100
A pop_data
object (from as_pop_data()
) with 200 rows and 8 columns:
Observation ID
Country where the participant was living
survey sampling cluster
survey catchment area
participant's age when sampled, in years
which antigen and isotype are being measured (data is in long format)
concentration of antigen isotype, in ELISA units
A subset of data from the SEES data, for examples and testing, data from Pakistan only.
sees_pop_data_pk_100
sees_pop_data_pk_100
sees_pop_data_pk_100
A pop_data
object (from as_pop_data()
) with 200 rows and 8 columns:
Observation ID
Country where the participant was living
survey sampling cluster
survey catchment area
participant's age when sampled, in years
which antigen and isotype are being measured (data is in long format)
concentration of antigen isotype, in ELISA units
The serocalculator package comes bundled with a number of sample files
in its inst/extdata
directory.
This serocalculator_example()
function make those sample files
easy to access.
serocalculator_example(file = NULL)
serocalculator_example(file = NULL)
file |
Name of file. If |
Adapted from readr::readr_example()
following the guidance in
https://r-pkgs.org/data.html#sec-data-example-path-helper.
a character string providing
the path to the file specified by file
,
or a vector or available files if file = NULL
.
serocalculator_example() serocalculator_example("example_pop_data.csv")
serocalculator_example() serocalculator_example("example_pop_data.csv")
Makes a cross-sectional data set (age, y(t) set) and adds noise, if desired.
sim_pop_data( lambda = 0.1, n_samples = 100, age_range = c(0, 20), age_fixed = NA, antigen_isos = intersect(get_biomarker_levels(curve_params), rownames(noise_limits)), n_mcmc_samples = 0, renew_params = FALSE, add_noise = FALSE, curve_params, noise_limits, format = "wide", verbose = FALSE, ... )
sim_pop_data( lambda = 0.1, n_samples = 100, age_range = c(0, 20), age_fixed = NA, antigen_isos = intersect(get_biomarker_levels(curve_params), rownames(noise_limits)), n_mcmc_samples = 0, renew_params = FALSE, add_noise = FALSE, curve_params, noise_limits, format = "wide", verbose = FALSE, ... )
lambda |
a |
n_samples |
number of samples to simulate |
age_range |
age range of sampled individuals, in years |
age_fixed |
specify the curve parameters to use by age (does nothing at present?) |
antigen_isos |
Character vector with one or more antibody names.
Values must match |
n_mcmc_samples |
how many MCMC samples to use:
|
renew_params |
whether to generate a new parameter set for each infection
|
add_noise |
a |
curve_params |
a
|
noise_limits |
biologic noise distribution parameters |
format |
a
|
verbose |
logical: if TRUE, print verbose log information to console |
... |
Arguments passed on to |
a tibble::tbl_df containing simulated cross-sectional serosurvey data, with columns:
age
: age (in days)
one column for each element in the antigen_iso
input argument
# Load curve parameters dmcmc <- typhoid_curves_nostrat_100 # Specify the antibody-isotype responses to include in analyses antibodies <- c("HlyE_IgA", "HlyE_IgG") # Set seed to reproduce results set.seed(54321) # Simulated incidence rate per person-year lambda <- 0.2 # Range covered in simulations lifespan <- c(0, 10) # Cross-sectional sample size nrep <- 100 # Biologic noise distribution dlims <- rbind( "HlyE_IgA" = c(min = 0, max = 0.5), "HlyE_IgG" = c(min = 0, max = 0.5) ) # Generate cross-sectional data csdata <- sim_pop_data( curve_params = dmcmc, lambda = lambda, n_samples = nrep, age_range = lifespan, antigen_isos = antibodies, n_mcmc_samples = 0, renew_params = TRUE, add_noise = TRUE, noise_limits = dlims, format = "long" )
# Load curve parameters dmcmc <- typhoid_curves_nostrat_100 # Specify the antibody-isotype responses to include in analyses antibodies <- c("HlyE_IgA", "HlyE_IgG") # Set seed to reproduce results set.seed(54321) # Simulated incidence rate per person-year lambda <- 0.2 # Range covered in simulations lifespan <- c(0, 10) # Cross-sectional sample size nrep <- 100 # Biologic noise distribution dlims <- rbind( "HlyE_IgA" = c(min = 0, max = 0.5), "HlyE_IgG" = c(min = 0, max = 0.5) ) # Generate cross-sectional data csdata <- sim_pop_data( curve_params = dmcmc, lambda = lambda, n_samples = nrep, age_range = lifespan, antigen_isos = antibodies, n_mcmc_samples = 0, renew_params = TRUE, add_noise = TRUE, noise_limits = dlims, format = "long" )
Simulate multiple data sets
sim_pop_data_multi( nclus = 10, lambdas = c(0.05, 0.1, 0.15, 0.2, 0.3), num_cores = max(1, parallel::detectCores() - 1), rng_seed = 1234, verbose = FALSE, ... )
sim_pop_data_multi( nclus = 10, lambdas = c(0.05, 0.1, 0.15, 0.2, 0.3), num_cores = max(1, parallel::detectCores() - 1), rng_seed = 1234, verbose = FALSE, ... )
nclus |
number of clusters |
lambdas |
#incidence rate, in events/person*year |
num_cores |
number of cores to use for parallel computations |
rng_seed |
starting seed for random number generator,
passed to |
verbose |
whether to report verbose information |
... |
Arguments passed on to
|
# Load curve parameters dmcmc <- typhoid_curves_nostrat_100 # Specify the antibody-isotype responses to include in analyses antibodies <- c("HlyE_IgA", "HlyE_IgG") # Set seed to reproduce results set.seed(54321) # Simulated incidence rate per person-year lambdas = c(.05, .1, .15, .2, .3) # Range covered in simulations lifespan <- c(0, 10); # Cross-sectional sample size nrep <- 100 # Biologic noise distribution dlims <- rbind( "HlyE_IgA" = c(min = 0, max = 0.5), "HlyE_IgG" = c(min = 0, max = 0.5) ) sim_pop_data_multi( curve_params = dmcmc, lambdas = lambdas, n_samples = nrep, age_range = lifespan, antigen_isos = antibodies, n_mcmc_samples = 0, renew_params = TRUE, add_noise = TRUE, noise_limits = dlims, format = "long", nclus = 10)
# Load curve parameters dmcmc <- typhoid_curves_nostrat_100 # Specify the antibody-isotype responses to include in analyses antibodies <- c("HlyE_IgA", "HlyE_IgG") # Set seed to reproduce results set.seed(54321) # Simulated incidence rate per person-year lambdas = c(.05, .1, .15, .2, .3) # Range covered in simulations lifespan <- c(0, 10); # Cross-sectional sample size nrep <- 100 # Biologic noise distribution dlims <- rbind( "HlyE_IgA" = c(min = 0, max = 0.5), "HlyE_IgG" = c(min = 0, max = 0.5) ) sim_pop_data_multi( curve_params = dmcmc, lambdas = lambdas, n_samples = nrep, age_range = lifespan, antigen_isos = antibodies, n_mcmc_samples = 0, renew_params = TRUE, add_noise = TRUE, noise_limits = dlims, format = "long", nclus = 10)
Generic method for extracting strata from objects.
See strata.seroincidence.by()
strata(x)
strata(x)
x |
an object |
the strata of x
This function is a summary()
method for seroincidence
objects.
## S3 method for class 'seroincidence' summary(object, coverage = 0.95, verbose = TRUE, ...)
## S3 method for class 'seroincidence' summary(object, coverage = 0.95, verbose = TRUE, ...)
object |
a |
coverage |
desired confidence interval coverage probability |
verbose |
whether to produce verbose messaging |
... |
unused |
a tibble::tibble()
containing the following:
est.start
: the starting guess for incidence rate
ageCat
: the age category we are analyzing
incidence.rate
: the estimated incidence rate, per person year
CI.lwr
: lower limit of confidence interval for incidence rate
CI.upr
: upper limit of confidence interval for incidence rate
coverage
: coverage probability
log.lik
:
log-likelihood of the data used in the call to est.incidence()
,
evaluated at the maximum-likelihood estimate of lambda
(i.e., at incidence.rate
)
iterations
: the number of iterations used
antigen_isos
: a list of antigen isotypes used in the analysis
nlm.convergence.code
:
information about convergence of the likelihood maximization procedure
performed by nlm()
(see "Value" section of stats::nlm()
, component code
);
codes 3-5 indicate issues:
1: relative gradient is close to zero, current iterate is probably solution.
2: successive iterates within tolerance, current iterate is probably solution.
3: Last global step failed to locate a point lower than x.
Either x is an approximate local minimum of the function,
the function is too non-linear for this algorithm,
or stepmin
in est.incidence()
(a.k.a., steptol
in stats::nlm()
) is too large.
4: iteration limit exceeded; increase iterlim
.
5: maximum step size stepmax
exceeded five consecutive times.
Either the function is unbounded below,
becomes asymptotic to a finite value from above in some direction,
or stepmax
is too small.
library(dplyr) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 |> filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est1 <- est.incidence( pop_data = xs_data, curve_params = curve, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA") ) summary(est1)
library(dplyr) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 |> filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk est1 <- est.incidence( pop_data = xs_data, curve_params = curve, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA") ) summary(est1)
"seroincidence.by"
ObjectsCalculate seroincidence from output of the seroincidence calculator
est.incidence.by()
.
## S3 method for class 'seroincidence.by' summary( object, confidence_level = 0.95, show_deviance = TRUE, show_convergence = TRUE, verbose = FALSE, ... )
## S3 method for class 'seroincidence.by' summary( object, confidence_level = 0.95, show_deviance = TRUE, show_convergence = TRUE, verbose = FALSE, ... )
object |
A dataframe containing output of function |
confidence_level |
desired confidence interval coverage probability |
show_deviance |
Logical flag ( |
show_convergence |
Logical flag ( |
verbose |
a logical scalar indicating whether to print verbose messages to the console |
... |
Additional arguments affecting the summary produced. |
A summary.seroincidence.by
object, which is a tibble::tibble,
with the following columns:
incidence.rate
maximum likelihood estimate of lambda
(seroincidence)
CI.lwr
lower confidence bound for lambda
CI.upr
upper confidence bound for lambda
Deviance
(included if show_deviance = TRUE
)
Negative log likelihood (NLL) at estimated (maximum likelihood) lambda
)
nlm.convergence.code
(included if show_convergence = TRUE
)
Convergence information returned by stats::nlm()
The object also has the following metadata
(accessible through base::attr()
):
antigen_isos
Character vector with names of input antigen isotypes
used in est.incidence.by()
Strata
Character with names of strata used in est.incidence.by()
library(dplyr) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 |> filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk # estimate seroincidence est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data, curve_params = curve, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), # num_cores = 8 # Allow for parallel processing to decrease run time ) # calculate summary statistics for the seroincidence object summary(est2)
library(dplyr) xs_data <- sees_pop_data_pk_100 curve <- typhoid_curves_nostrat_100 |> filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) noise <- example_noise_params_pk # estimate seroincidence est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data, curve_params = curve, noise_params = noise, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), # num_cores = 8 # Allow for parallel processing to decrease run time ) # calculate summary statistics for the seroincidence object summary(est2)
A subset of data from the SEES study, for examples and testing.
typhoid_curves_nostrat_100
typhoid_curves_nostrat_100
typhoid_curves_nostrat_100
A curve_params
object (from as_curve_params()
) with 500 rows and 7
columns:
which antigen and isotype are being measured (data is in long format)
MCMC iteration
Antibody concentration at t = 0 (start of active infection)
Antibody concentration at t = t1
(end of active infection)
Duration of active infection
Antibody decay rate coefficient
Antibody decay rate exponent parameter