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.2.0.9022 |
Built: | 2024-11-16 03:06:48 UTC |
Source: | https://github.com/ucd-serg/serocalculator |
seroincidence.by
objectExtract or replace parts of a seroincidence.by
object
## S3 method for class 'seroincidence.by' x[i, ...]
## S3 method for class 'seroincidence.by' x[i, ...]
x |
the object to subset/replace elements of |
i |
the indices to subset/replace |
... |
passed to |
the subset specified
kinetics of the antibody (ab) response (power function decay)
ab(t, par, ...)
ab(t, par, ...)
t |
age at infection? |
par |
parameters |
... |
arguments passed to |
a matrix()
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 <- "https://osf.io/download/rtw5k/" %>% readr::read_rds() %>% as_curve_params() print(curve_data)
library(magrittr) curve_data <- "https://osf.io/download/rtw5k/" %>% readr::read_rds() %>% 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 <- "https://osf.io/download//hqy4v/" %>% readr::read_rds() %>% as_noise_params() print(noise_data)
library(magrittr) noise_data <- "https://osf.io/download//hqy4v/" %>% readr::read_rds() %>% 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 <- "https://osf.io/download//n6cp3/" %>% readr::read_rds() %>% as_pop_data() print(xs_data)
library(magrittr) xs_data <- "https://osf.io/download//n6cp3/" %>% readr::read_rds() %>% 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
Note that 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 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 = load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) %>% # Reduce dataset for this example autoplot() curve
library(dplyr) library(ggplot2) library(magrittr) curve = load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) %>% # Reduce dataset for this example 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) xs_data <- load_pop_data( file_path = "https://osf.io/download//n6cp3/", age = "Age", id = "index_id", value = "result", standardize = TRUE ) xs_data %>% autoplot(strata = "Country", type = "density") xs_data %>% autoplot(strata = "Country", type = "age-scatter")
library(dplyr) library(ggplot2) xs_data <- load_pop_data( file_path = "https://osf.io/download//n6cp3/", age = "Age", id = "index_id", value = "result", standardize = TRUE ) xs_data %>% autoplot(strata = "Country", type = "density") xs_data %>% autoplot(strata = "Country", 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 <- load_pop_data("https://osf.io/download//n6cp3/") curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est1 <- est.incidence( pop_data = xs_data %>% filter(Country == "Pakistan"), curve_param = curve, noise_param = noise %>% filter(Country == "Pakistan"), antigen_isos = c("HlyE_IgG", "HlyE_IgA"), build_graph = TRUE ) # Plot the log-likelihood curve autoplot(est1)
library(dplyr) library(ggplot2) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est1 <- est.incidence( pop_data = xs_data %>% filter(Country == "Pakistan"), curve_param = curve, noise_param = noise %>% filter(Country == "Pakistan"), 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 <- "https://osf.io/download//n6cp3/" %>% load_pop_data() curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curve, noise_params = noise %>% filter(Country == "Pakistan"), 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 <- "https://osf.io/download//n6cp3/" %>% load_pop_data() curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curve, noise_params = noise %>% filter(Country == "Pakistan"), 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, width = 0.001, ...)
## S3 method for class 'summary.seroincidence.by' autoplot(object, xvar, alpha = 0.7, shape = 1, width = 0.001, ...)
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 |
width |
width for jitter |
... |
unused |
a ggplot2::ggplot()
object
library(dplyr) library(ggplot2) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curve, noise_params = noise %>% filter(Country == "Pakistan"), antigen_isos = c("HlyE_IgG", "HlyE_IgA"), #num_cores = 8 #Allow for parallel processing to decrease run time ) est2sum <- summary(est2) autoplot(est2sum, "catchment")
library(dplyr) library(ggplot2) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curve, noise_params = noise %>% filter(Country == "Pakistan"), antigen_isos = c("HlyE_IgG", "HlyE_IgA"), #num_cores = 8 #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(dplyr) xs_data <- readr::read_rds("https://osf.io/download//n6cp3/") %>% as_pop_data() check_pop_data(xs_data, verbose = TRUE)
library(dplyr) xs_data <- readr::read_rds("https://osf.io/download//n6cp3/") %>% 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 = pop_data$antigen_iso %>% unique(), 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 = pop_data$antigen_iso %>% unique(), 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 <- load_pop_data("https://osf.io/download//n6cp3/") curves <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est1 <- est.incidence( pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curves, noise_params = noise %>% filter(Country == "Pakistan"), antigen_isos = c("HlyE_IgG", "HlyE_IgA"), iterlim = 5 # limit iterations for the purpose of this example ) summary(est1)
library(dplyr) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") curves <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est1 <- est.incidence( pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curves, noise_params = noise %>% filter(Country == "Pakistan"), antigen_isos = c("HlyE_IgG", "HlyE_IgA"), iterlim = 5 # limit iterations for the purpose of this example ) 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 <- load_pop_data("https://osf.io/download//n6cp3/") curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% # Reduce dataset for the purposes of this example: slice(1:100, .by = antigen_iso) noise <- load_noise_params("https://osf.io/download//hqy4v/") est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curve, noise_params = noise %>% filter(Country == "Pakistan"), 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 ) summary(est2)
library(dplyr) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% # Reduce dataset for the purposes of this example: slice(1:100, .by = antigen_iso) noise <- load_noise_params("https://osf.io/download//hqy4v/") est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curve, noise_params = noise %>% filter(Country == "Pakistan"), 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 ) summary(est2)
A subset of noise parameter estimates from the SEES study, for examples and testing.
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.
Retrieves additional data from internet. The data format must be .RDS or a zipped .RDS. The purpose of this function is to download data such as longitudinal response parameters from an online repository or cross-sectional population data.
Data for this package is available at: https://osf.io/ne8pc/files/osfstorage
You can save the data into your chosen directory using the optional savePath argument. Specify the file path and the file name.
Large datasets may timeout. If so, you can increase the download time by updating the maximum timeout time in the code below. (Ex: increase from 300 to 1000)
options(timeout = max(300, getOption("timeout")))
get_additional_data(fileURL, savePath = NULL)
get_additional_data(fileURL, savePath = NULL)
fileURL |
URL of the file to be downloaded. |
savePath |
Folder directory and filename to save the downloaded and unzipped (if needed) file. File is saved only
if this argument is not |
the R object stored in the file indicated by the fileURL
input
Data object
## Not run: curve_param_samples <- get_additional_data( fileURL = "https://osf.io/download/bhfvx" ) # optionally, save the data to disk curve_param_samples <- get_additional_data( fileURL = "https://osf.io/download/bhfvx", savePath = "~/Downloads/curv_params.rds" ) ## End(Not run)
## Not run: curve_param_samples <- get_additional_data( fileURL = "https://osf.io/download/bhfvx" ) # optionally, save the data to disk curve_param_samples <- get_additional_data( fileURL = "https://osf.io/download/bhfvx", savePath = "~/Downloads/curv_params.rds" ) ## End(Not run)
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 <- load_pop_data("https://osf.io/download//n6cp3/") # Load curve parameters and subset for the purposes of this example dmcmc <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # 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 = dmcmc, 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 <- load_pop_data("https://osf.io/download//n6cp3/") # Load curve parameters and subset for the purposes of this example dmcmc <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # 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 = dmcmc, 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 )
graph.curve.params( curve_params, antigen_isos = unique(curve_params$antigen_iso), verbose = FALSE )
curve_params |
a |
antigen_isos |
antigen isotypes |
verbose |
verbose output |
a ggplot2::ggplot()
object
## Not run: plot1 <- graph.curve.params(curve_params) print(plot1) ## End(Not run)
## Not run: plot1 <- graph.curve.params(curve_params) print(plot1) ## End(Not run)
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("https://osf.io/download/rtw5k/") print(curve)
curve <- load_curve_params("https://osf.io/download/rtw5k/") 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("https://osf.io/download//hqy4v/") print(noise)
noise <- load_noise_params("https://osf.io/download//hqy4v/") 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("https://osf.io/download//n6cp3/") print(xs_data)
xs_data <- load_pop_data("https://osf.io/download//n6cp3/") 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 in longitudinal parameters dmcmc <- load_curve_params("https://osf.io/download/rtw5k") xs_data <- "https://osf.io/download//n6cp3/" %>% load_pop_data() # 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 = dmcmc, noise_params = cond, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), lambda = 0.1 ) %>% print()
library(dplyr) library(tibble) # load in longitudinal parameters dmcmc <- load_curve_params("https://osf.io/download/rtw5k") xs_data <- "https://osf.io/download//n6cp3/" %>% load_pop_data() # 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 = dmcmc, noise_params = cond, antigen_isos = c("HlyE_IgG", "HlyE_IgA"), lambda = 0.1 ) %>% print()
generate random sample from baseline distribution
mk_baseline(kab, n = 1, blims, ...)
mk_baseline(kab, n = 1, blims, ...)
kab |
index for which row of antibody baseline limits to read from |
n |
number of observations |
blims |
range of possible baseline antibody levels |
... |
not currently used |
a numeric()
vector
Graph an antibody decay curve model
plot_curve_params_one_ab( object, verbose = FALSE, alpha = 0.4, n_curves = 100, n_points = 1000, log_x = FALSE, log_y = TRUE, rows_to_graph = seq_len(min(n_curves, nrow(object))), xlim = c(10^-1, 10^3.1), ... )
plot_curve_params_one_ab( object, verbose = FALSE, alpha = 0.4, n_curves = 100, n_points = 1000, log_x = FALSE, log_y = TRUE, rows_to_graph = seq_len(min(n_curves, nrow(object))), xlim = c(10^-1, 10^3.1), ... )
object |
a |
verbose |
verbose output |
alpha |
(passed to
|
n_curves |
how many curves to plot (see details). |
n_points |
Number of points to interpolate along the x axis
(passed to |
log_x |
should the x-axis be on a logarithmic scale ( |
log_y |
should the Y-axis be on a logarithmic scale
(default, |
rows_to_graph |
which rows of |
xlim |
range of x values to graph |
... |
Arguments passed on to
|
n_curves
and rows_to_graph
In most cases, curve_params
will contain too many rows of MCMC
samples for all of these samples to be plotted at once.
Setting the n_curves
argument to a value smaller than the
number of rows in curve_params
will cause this function to select
the first n_curves
rows to graph.
Setting n_curves
larger than the number of rows in ' will
result all curves being plotted.
If the user directly specifies the rows_to_graph
argument,
then n_curves
has no effect.
a ggplot2::ggplot()
object
library(dplyr) # loads the `%>%` operator and `dplyr::filter()` load_curve_params("https://osf.io/download/rtw5k/") %>% dplyr::filter(antigen_iso == "HlyE_IgG") %>% serocalculator:::plot_curve_params_one_ab()
library(dplyr) # loads the `%>%` operator and `dplyr::filter()` load_curve_params("https://osf.io/download/rtw5k/") %>% dplyr::filter(antigen_iso == "HlyE_IgG") %>% serocalculator:::plot_curve_params_one_ab()
seroincidence
ObjectCustom print()
function to show output of the seroincidence calculator est.incidence()
.
## S3 method for class 'seroincidence' print(x, ...)
## S3 method for class 'seroincidence' print(x, ...)
x |
A list containing output of function |
... |
Additional arguments affecting the summary produced. |
## Not run: # Estimate seroincidence seroincidence <- est.incidence.by(...) # Print the seroincidence object to the console print(seroincidence) # Or simply type (appropriate print method will be invoked automatically) seroincidence ## End(Not run)
## Not run: # Estimate seroincidence seroincidence <- est.incidence.by(...) # Print the seroincidence object to the console print(seroincidence) # Or simply type (appropriate print method will be invoked automatically) seroincidence ## End(Not run)
seroincidence.by
ObjectCustom print()
function to show output of the seroincidence calculator est.incidence.by()
.
## S3 method for class 'seroincidence.by' print(x, ...)
## S3 method for class 'seroincidence.by' print(x, ...)
x |
A list containing output of function |
... |
Additional arguments affecting the summary produced. |
## Not run: # Estimate seroincidence seroincidence <- est.incidence.by(...) # Print the seroincidence object to the console print(seroincidence) # Or simply type (appropriate print method will be invoked automatically) seroincidence ## End(Not run)
## Not run: # Estimate seroincidence seroincidence <- est.incidence.by(...) # Print the seroincidence object to the console print(seroincidence) # Or simply type (appropriate print method will be invoked automatically) seroincidence ## End(Not run)
Custom print()
function for "summary.seroincidence.by" objects (constructed by summary.seroincidence.by()
)
## S3 method for class 'summary.seroincidence.by' print(x, ...)
## S3 method for class 'summary.seroincidence.by' print(x, ...)
x |
A "summary.seroincidence.by" object (constructed by |
... |
Additional arguments affecting the summary produced. |
## Not run: # Estimate seroincidence seroincidence <- est.incidence.by(...) # Calculate summary statistics for the seroincidence object seroincidenceSummary <- summary(seroincidence) # Print the summary of seroincidence object to the console print(seroincidenceSummary) # Or simply type (appropriate print method will be invoked automatically) seroincidenceSummary ## End(Not run)
## Not run: # Estimate seroincidence seroincidence <- est.incidence.by(...) # Calculate summary statistics for the seroincidence object seroincidenceSummary <- summary(seroincidence) # Print the summary of seroincidence object to the console print(seroincidenceSummary) # Or simply type (appropriate print method will be invoked automatically) seroincidenceSummary ## End(Not run)
take a random sample from longitudinal parameter set given age at infection, for a list of antibodies
row_longitudinal_parameter(age, antigen_isos, nmc, npar, ...)
row_longitudinal_parameter(age, antigen_isos, nmc, npar, ...)
age |
age at infection |
antigen_isos |
antigen isotypes |
nmc |
mcmc sample to use |
npar |
number of parameters |
... |
passed to |
an array of parameters: c(y0,b0,mu0,mu1,c1,alpha,shape)
A subset of data from the SEES data, for examples and testing.
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
A subset of data from the SEES data, for examples and testing.
sees_pop_data_pk_100_old_names
sees_pop_data_pk_100_old_names
sees_pop_data_pk_100_old_names
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
This package translates antibody levels measured in a (cross-sectional) population sample into an estimate of the frequency with which seroconversions (infections) occur in the sampled population.
_PACKAGE
Peter Teunis [email protected]
Doug Ezra Morrison [email protected]
Kristen Aiemjoy [email protected]
Kristina Lai [email protected]
Methods for estimating seroincidence
Teunis, P. F. M., and J. C. H. van Eijkeren. "Estimation of seroconversion rates for infectious diseases: Effects of age and noise." Statistics in Medicine 39, no. 21 (2020): 2799-2814.
Teunis, P. F. M., J. C. H. van Eijkeren, W. F. de Graaf, A. Bonačić Marinović, and M. E. E. Kretzschmar. "Linking the seroresponse to infection to within-host heterogeneity in antibody production." Epidemics 16 (2016): 33-39.
Applications
Aiemjoy, Kristen, Jessica C. Seidman, Senjuti Saha, Sira Jam Munira, Mohammad Saiful Islam Sajib, Syed Muktadir Al Sium, Anik Sarkar et al. "Estimating typhoid incidence from community-based serosurveys: a multicohort study." The Lancet Microbe 3, no. 8 (2022): e578-e587.
Aiemjoy, Kristen, John Rumunu, Juma John Hassen, Kirsten E. Wiens, Denise Garrett, Polina Kamenskaya, Jason B. Harris et al. "Seroincidence of enteric fever, Juba, South Sudan." Emerging infectious diseases 28, no. 11 (2022): 2316.
Monge, S., Teunis, P. F., Friesema, I., Franz, E., Ang, W., van Pelt, W., Mughini-Gras, L. "Immune response-eliciting exposure to Campylobacter vastly exceeds the incidence of clinically overt campylobacteriosis but is associated with similar risk factors: A nationwide serosurvey in the Netherlands" Journal of Infection, 2018, 1–7, doi:10.1016/j.jinf.2018.04.016
Kretzschmar, M., Teunis, P. F., Pebody, R. G. "Incidence and reproduction numbers of pertussis: estimates from serological and social contact data in five European countries" PLoS Medicine 7, no. 6 (June 1, 2010):e1000291. doi:10.1371/journal.pmed.1000291.
Simonsen, J., Strid, M. A., Molbak, K., Krogfelt, K. A., Linneberg, A., Teunis, P. "Sero-epidemiology as a tool to study the incidence of Salmonella infections in humans" Epidemiology and Infection 136, no. 7 (July 1, 2008): 895–902. doi:10.1017/S0950268807009314
Simonsen, J., Teunis, P. F., van Pelt, W., van Duynhoven, Y., Krogfelt, K. A., Sadkowska-Todys, M., Molbak K. "Usefulness of seroconversion rates for comparing infection pressures between countries" Epidemiology and Infection, April 12, 2010, 1–8. doi:10.1017/S0950268810000750.
Falkenhorst, G., Simonsen, J., Ceper, T. H., van Pelt, W., de Valk, H., Sadkowska-Todys, M., Zota, L., Kuusi, M., Jernberg, C., Rota, M. C., van Duynhoven, Y. T., Teunis, P. F., Krogfelt, K. A., Molbak, K. "Serological cross-sectional studies on salmonella incidence in eight European countries: no correlation with incidence of reported cases" BMC Public Health 12, no. 1 (July 15, 2012): 523–23. doi:10.1186/1471-2458-12-523.
Teunis, P. F., Falkenhorst, G., Ang, C. W., Strid, M. A., De Valk, H., Sadkowska-Todys, M., Zota, L., Kuusi, M., Rota, M. C., Simonsen, J. B., Molbak, K., Van Duynhoven, Y. T., van Pelt, W. "Campylobacter seroconversion rates in selected countries in the European Union" Epidemiology and Infection 141, no. 10 (2013): 2051–57. doi:10.1017/S0950268812002774.
de Melker, H. E., Versteegh, F. G., Schellekens, J. F., Teunis, P. F., Kretzschmar, M. "The incidence of Bordetella pertussis infections estimated in the population from a combination of serological surveys" The Journal of Infection 53, no. 2 (August 1, 2006): 106–13. doi:10.1016/j.jinf.2005.10.020
Quantification of seroresponse
de Graaf, W. F., Kretzschmar, M. E., Teunis, P. F., Diekmann, O. "A two-phase within-host model for immune response and its application to serological profiles of pertussis" Epidemics 9 (2014):1–7. doi:10.1016/j.epidem.2014.08.002.
Berbers, G. A., van de Wetering, M. S., van Gageldonk, P. G., Schellekens, J. F., Versteegh, F. G., Teunis, P.F. "A novel method for evaluating natural and vaccine induced serological responses to Bordetella pertussis antigens" Vaccine 31, no. 36 (August 12, 2013): 3732–38. doi:10.1016/j.vaccine.2013.05.073.
Versteegh, F. G., Mertens, P. L., de Melker, H. E., Roord, J. J., Schellekens, J. F., Teunis, P. F. "Age-specific long-term course of IgG antibodies to pertussis toxin after symptomatic infection with Bordetella pertussis" Epidemiology and Infection 133, no. 4 (August 1, 2005): 737–48.
Teunis, P. F., van der Heijden, O. G., de Melker, H. E., Schellekens, J. F., Versteegh, F. G., Kretzschmar, M. E. "Kinetics of the IgG antibody response to pertussis toxin after infection with B. pertussis" Epidemiology and Infection 129, no. 3 (December 10, 2002):479. doi:10.1017/S0950268802007896.
Makes a cross-sectional data set (age, y(t) set) and adds noise, if desired.
sim.cs( lambda = 0.1, n.smpl = 100, age.rng = c(0, 20), age.fx = NA, antigen_isos, n.mc = 0, renew.params = FALSE, add.noise = FALSE, curve_params, noise_limits, format = "wide", verbose = FALSE, ... )
sim.cs( lambda = 0.1, n.smpl = 100, age.rng = c(0, 20), age.fx = NA, antigen_isos, n.mc = 0, renew.params = FALSE, add.noise = FALSE, curve_params, noise_limits, format = "wide", verbose = FALSE, ... )
lambda |
a |
n.smpl |
number of samples to simulate |
age.rng |
age range of sampled individuals, in years |
age.fx |
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.mc |
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 <- load_curve_params("https://osf.io/download/rtw5k") # 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.cs( curve_params = dmcmc, lambda = lambda, n.smpl = nrep, age.rng = lifespan, antigen_isos = antibodies, n.mc = 0, renew.params = TRUE, add.noise = TRUE, noise_limits = dlims, format = "long" )
# Load curve parameters dmcmc <- load_curve_params("https://osf.io/download/rtw5k") # 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.cs( curve_params = dmcmc, lambda = lambda, n.smpl = nrep, age.rng = lifespan, antigen_isos = antibodies, n.mc = 0, renew.params = TRUE, add.noise = TRUE, noise_limits = dlims, format = "long" )
Simulate multiple data sets
sim.cs.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, renew.params = TRUE, add.noise = TRUE, verbose = FALSE, ... )
sim.cs.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, renew.params = TRUE, add.noise = TRUE, 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 |
renew.params |
whether to generate a new parameter set for each infection
|
add.noise |
a |
verbose |
whether to report verbose information |
... |
Arguments passed on to
|
output: (age, y(t) set)
simcs.tinf( lambda, n.smpl, age.rng, age.fx = NA, antigen_isos, n.mc = 0, renew.params = FALSE, ... )
simcs.tinf( lambda, n.smpl, age.rng, age.fx = NA, antigen_isos, n.mc = 0, renew.params = FALSE, ... )
lambda |
seroconversion rate (in events/person-day) |
n.smpl |
number of samples n.smpl (= nr of simulated records) |
age.rng |
age range to use for simulating data, in days |
age.fx |
age.fx for parameter sample (age.fx = NA for age at infection) |
antigen_isos |
Character vector with one or more antibody names. Values must match |
n.mc |
|
renew.params |
|
... |
arguments passed to |
an array()
simulate antibody kinetics of y over a time interval
simresp.tinf( lambda, t.end, age.fx, antigen_isos, n.mc = 0, renew.params, predpar, ... )
simresp.tinf( lambda, t.end, age.fx, antigen_isos, n.mc = 0, renew.params, predpar, ... )
lambda |
seroconversion rate (1/days), |
t.end |
end of time interval (beginning is time 0) in days(?) |
age.fx |
parameter estimates for fixed age (age.fx in years) or not. when age.fx = NA then age at infection is used. |
antigen_isos |
antigen isotypes |
n.mc |
a posterior sample may be selected (1:4000), or not when n.mc = 0 a posterior sample is chosen at random. |
renew.params |
At infection, a new parameter sample may be generated (when |
predpar |
an
|
... |
Arguments passed on to
|
This function returns a list()
with:
t = times (in days, birth at day 0),
b = bacteria level, for each antibody signal (not used; probably meaningless),
y = antibody level, for each antibody signal
smp = whether an infection involves a big jump or a small jump
t.inf = times when infections have occurred.
Generic method for extracting strata from objects. See strata.seroincidence.by()
strata(x)
strata(x)
x |
an object |
the strata of x
Strata
attribute from an object, if presentExtract the Strata
attribute from an object, if present
## S3 method for class 'seroincidence.by' strata(x)
## S3 method for class 'seroincidence.by' strata(x)
x |
any R object |
a tibble::tibble()
with strata in rows, or
NULL
if x
does not have a "strata"
attribute
summary()
method for pop_data
objects
## S3 method for class 'pop_data' summary(object, strata = NULL, ...) ## S3 method for class 'summary.pop_data' print(x, ...)
## S3 method for class 'pop_data' summary(object, strata = NULL, ...) ## S3 method for class 'summary.pop_data' print(x, ...)
object |
a |
strata |
a |
... |
unused |
x |
an object of class |
a summary.pop_data
object, which is a list containing two summary tables:
age_summary
summarizing age
ab_summary
summarizing value
, stratified by antigen_iso
library(dplyr) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") summary(xs_data, strata = "Country")
library(dplyr) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") summary(xs_data, strata = "Country")
This function is a summary()
method for seroincidence
objects.
## S3 method for class 'seroincidence' summary(object, coverage = 0.95, ...)
## S3 method for class 'seroincidence' summary(object, coverage = 0.95, ...)
object |
a |
coverage |
desired confidence interval coverage probability |
... |
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 <- load_pop_data("https://osf.io/download//n6cp3/") curves <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est1 <- est.incidence( pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curves, noise_params = noise %>% filter(Country == "Pakistan"), antigen_isos = c("HlyE_IgG", "HlyE_IgA") ) summary(est1)
library(dplyr) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") curves <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est1 <- est.incidence( pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curves, noise_params = noise %>% filter(Country == "Pakistan"), 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, showDeviance = TRUE, showConvergence = TRUE, ... )
## S3 method for class 'seroincidence.by' summary( object, confidence_level = 0.95, showDeviance = TRUE, showConvergence = TRUE, ... )
object |
A dataframe containing output of function |
confidence_level |
desired confidence interval coverage probability |
showDeviance |
Logical flag ( |
showConvergence |
Logical flag ( |
... |
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 showDeviance = TRUE
) Negative log likelihood (NLL) at estimated (maximum likelihood)
lambda
)
nlm.convergence.code
(included if showConvergence = 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 <- load_pop_data("https://osf.io/download//n6cp3/") curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curve, noise_params = noise %>% filter(Country == "Pakistan"), antigen_isos = c("HlyE_IgG", "HlyE_IgA"), #num_cores = 8 # Allow for parallel processing to decrease run time ) summary(est2) ## Not run: # estimate seroincidence seroincidence <- est.incidence.by(...) # calculate summary statistics for the seroincidence object seroincidenceSummary <- summary(seroincidence) ## End(Not run)
library(dplyr) xs_data <- load_pop_data("https://osf.io/download//n6cp3/") curve <- load_curve_params("https://osf.io/download/rtw5k/") %>% filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>% slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example noise <- load_noise_params("https://osf.io/download//hqy4v/") est2 <- est.incidence.by( strata = c("catchment"), pop_data = xs_data %>% filter(Country == "Pakistan"), curve_params = curve, noise_params = noise %>% filter(Country == "Pakistan"), antigen_isos = c("HlyE_IgG", "HlyE_IgA"), #num_cores = 8 # Allow for parallel processing to decrease run time ) summary(est2) ## Not run: # estimate seroincidence seroincidence <- est.incidence.by(...) # calculate summary statistics for the seroincidence object seroincidenceSummary <- summary(seroincidence) ## End(Not run)
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
Warn about missing stratifying variables in a dataset
warn.missing.strata(data, strata, dataname)
warn.missing.strata(data, strata, dataname)
data |
the dataset that should contain the strata |
strata |
a |
dataname |
the name of the dataset, for use in warning messages if some strata are missing. |
a character()
vector of the subset of stratifying variables that are present in pop_data
## Not run: expected_strata <- data.frame(Species = "banana", type = "orchid") warn.missing.strata(iris, expected_strata, dataname = "iris") ## End(Not run)
## Not run: expected_strata <- data.frame(Species = "banana", type = "orchid") warn.missing.strata(iris, expected_strata, dataname = "iris") ## End(Not run)