Function to estimate seroincidences based on cross-sectional serology data and longitudinal response model.
Arguments
- pop_data
a data.frame with cross-sectional serology data per antibody and age, and additional columns corresponding to each element of the
strata
input- curve_params
a
data.frame()
containing MCMC samples of parameters from the Bayesian posterior distribution of a longitudinal decay curve model. The parameter columns must be named:antigen_iso
: acharacter()
vector indicating antigen-isotype combinationsiter
: aninteger()
vector indicating MCMC sampling iterationsy0
: baseline antibody level at $t=0$ ($y(t=0)$)y1
: antibody peak level (ELISA units)t1
: duration of infectionalpha
: antibody decay rate (1/days for the current longitudinal parameter sets)r
: shape factor of antibody decay
- noise_params
a
data.frame()
(ortibble::tibble()
) containing the following variables, specifying noise parameters for each antigen isotype:antigen_iso
: antigen isotype whose noise parameters are being specified on each rownu
: biological noiseeps
: measurement noisey.low
: lower limit of detection for the current antigen isotypey.high
: upper limit of detection for the current antigen isotype
- strata
a character vector of stratum-defining variables. Values must be variable names in
pop_data
.- curve_strata_varnames
A subset of
strata
. Values must be variable names incurve_params
. Default = "".- noise_strata_varnames
A subset of
strata
. Values must be variable names innoise_params
. Default = "".- antigen_isos
Character vector with one or more antibody names. Values must match
pop_data
- 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
est.incidence()
- ...
Arguments passed on to
est.incidence
,stats::nlm
stepmin
A positive scalar providing the minimum allowable relative step length.
stepmax
a positive scalar which gives the maximum allowable scaled step length.
stepmax
is used to prevent steps which would cause the optimization function to overflow, to prevent the algorithm from leaving the area of interest in parameter space, or to detect divergence in the algorithm.stepmax
would be chosen small enough to prevent the first two of these occurrences, but should be larger than any anticipated reasonable step.typsize
an estimate of the size of each parameter at the minimum.
fscale
an estimate of the size of
f
at the minimum.ndigit
the number of significant digits in the function
f
.gradtol
a positive scalar giving the tolerance at which the scaled gradient is considered close enough to zero to terminate the algorithm. The scaled gradient is a measure of the relative change in
f
in each directionp[i]
divided by the relative change inp[i]
.iterlim
a positive integer specifying the maximum number of iterations to be performed before the program is terminated.
check.analyticals
a logical scalar specifying whether the analytic gradients and Hessians, if they are supplied, should be checked against numerical derivatives at the initial parameter values. This can help detect incorrectly formulated gradients or Hessians.
Value
if
strata
has meaningful inputs: An object of class"seroincidence.by"
; i.e., a list of"seroincidence"
objects fromest.incidence()
, one for each stratum, with some meta-data attributes.if
strata
is missing,NULL
,NA
, or""
: An object of class"seroincidence"
.
Details
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.
Examples
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
)
#> Warning: curve_params is missing all strata variables, and will be used unstratified.
#>
#> To avoid this warning, specify the desired set of stratifying variables in the `curve_strata_varnames` and `noise_strata_varnames` arguments to `est.incidence.by()`.
#> Warning: noise_params is missing all strata variables, and will be used unstratified.
#>
#> To avoid this warning, specify the desired set of stratifying variables in the `curve_strata_varnames` and `noise_strata_varnames` arguments to `est.incidence.by()`.
summary(est2)
#> Seroincidence estimated given the following setup:
#> a) Antigen isotypes : HlyE_IgG, HlyE_IgA
#> b) Strata : catchment
#>
#> Seroincidence estimates:
#> # A tibble: 2 × 13
#> Stratum catchment n est.start incidence.rate SE CI.lwr CI.upr
#> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Stratum 1 aku 294 0.1 0.118 0.00825 0.103 0.136
#> 2 Stratum 2 kgh 200 0.1 0.183 0.0139 0.157 0.212
#> # ℹ 5 more variables: coverage <dbl>, log.lik <dbl>, iterations <int>,
#> # antigen.isos <chr>, nlm.convergence.code <ord>