Re-fit a subset of strata, with log-likelihood graphs attached

Description

Companion to est_seroincidence_by_cached(). Caching a fit built with build_graph = TRUE is impractical over many strata, because a log-likelihood plot is attached to each one; the usual pattern is to cache the whole fit with build_graph = FALSE and then re-fit, with graphs, only the handful of strata you actually want to plot.

Usage

refit_strata(pop_data, strata_ids, strata = names(strata_ids), ...)

Arguments

pop_data	a data.frame() of cross-sectional serology data: the same data the full fit was built from.
strata_ids	a data.frame() whose columns are the stratifying variables and whose rows select the strata to re-fit. Matched against pop_data with dplyr::semi_join().
strata	a character() of the stratifying variable names. Defaults to the columns of strata_ids.
…	Arguments passed on to est_seroincidence_by curve_strata_varnames A subset of strata. Values must be variable names in curve_params. Default = "". noise_strata_varnames A subset of strata. Values must be variable names in noise_params. Default = "". 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. lambda_start starting guess for incidence rate, in events/year. antigen_isos Character vector with one or more antibody names. Must match pop_data 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 est_seroincidence() sr_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: a character() vector indicating antigen-isotype combinations iter: an integer() vector indicating MCMC sampling iterations y0: baseline antibody level at $t=0$ ($y(t=0)$) y1: antibody peak level (ELISA units) t1: duration of infection alpha: antibody decay rate (1/days for the current longitudinal parameter sets) r: shape factor of antibody decay cluster_var optional name(s) of the variable(s) in pop_data containing cluster identifiers for clustered sampling designs (e.g., households, schools). Can be a single variable name (character string) or a vector of variable names for multi-level clustering (e.g., c(“school”, “classroom”)). When provided, standard errors will be adjusted for within-cluster correlation using cluster-robust variance estimation. When fitting more than one antigen_isos at once, this argument also has a second use. log_likelihood() combines biomarkers by summing their marginal log-likelihoods, which is only valid if those contributions are independent. Two biomarker readings from the same person usually aren’t, since they share an infection history. Pass the id column returned by ids_varname() (e.g. cluster_var = ids_varname(pop_data)) to get a cluster-robust standard error that accounts for this within-person correlation. This is a distinct concern from a genuine sampling-cluster design, and multiple cluster_var values are grouped by their intersection rather than combined additively: if every subject belongs to exactly one sampling cluster, cluster_var = c(“cluster”, “id”) reduces to clustering by id alone, and does not add the sampling-cluster correction on top — see issue #543. See issue #645. stratum_var optional name of the variable in pop_data containing stratum identifiers. Used in combination with cluster_var for stratified cluster sampling designs. noise_params a data.frame() (or tibble::tibble()) containing the following variables, specifying noise parameters for each antigen isotype: antigen_iso: antigen isotype whose noise parameters are being specified on each row nu: biological noise eps: measurement noise y.low: lower limit of detection for the current antigen isotype y.high: upper limit of detection for the current antigen isotype verbose logical: if TRUE, print verbose log information to console

Details

Strata are selected by the columns that define them rather than by index or by label, because count_strata() names strata positionally (paste(“Stratum”, row_number())). A subset therefore renumbers them, so labels and positions from the full fit do not carry over.

Value

a “seroincidence.by” object covering only the selected strata, with graphs attached, or NULL if strata_ids has no rows.

See Also

est_seroincidence_by_cached()

Examples

library("serocalculator")

ests <- est_seroincidence_by(
  pop_data = sees_pop_data_pk_100,
  sr_params = typhoid_curves_nostrat_100,
  noise_params = example_noise_params_pk,
  strata = "catchment",
  antigen_isos = c("HlyE_IgG", "HlyE_IgA"),
  build_graph = FALSE,
  iterlim = 5
)

# Re-fit just the strata whose optimizer flagged a possible problem.
problem_ids <- ests |>
  summary() |>
  tibble::as_tibble() |>
  dplyr::filter(.data$nlm.convergence.code > 2) |>
  dplyr::select("catchment")

refit_strata(
  pop_data = sees_pop_data_pk_100,
  strata_ids = problem_ids,
  sr_params = typhoid_curves_nostrat_100,
  noise_params = example_noise_params_pk,
  antigen_isos = c("HlyE_IgG", "HlyE_IgA"),
  iterlim = 5
)