Package 'anumaan'

Description

Summarises draw-level R_ALL, R_KNOWN_OUTCOME, and R_NF values from compute_event_profile_probabilities() into posterior mean and credible interval per hospital x pathogen x profile combination, and derives explicit eligible_for_YLL/eligible_for_YLD flags (with reasons) rather than asserting universal usability.

Usage

aggregate_profiles_for_daly(
  profile_output,
  hospital_col = "hospital",
  pathogen_col = "pathogen",
  estimand = "observed_stewardship_event_mix",
  ci_level = 0.95,
  min_n_events = 10L,
  min_n_draws = 100L,
  sampler_acceptable = TRUE
)

Arguments

Details

Cohort definitions: R_ALL is the mean over truly all events in the hospital-pathogen-panel cohort (no outcome filter). R_KNOWN_OUTCOME restricts to events with a known patient outcome (this is what earlier versions of this function called R_ALL – renamed because it is not actually all events). R_NF restricts to the non-fatal subset of the known-outcome cohort. All three are computed per posterior draw and then summarised here; n_draws_all/n_draws_known_outcome/n_draws_nf count only the non-NA draws feeding each summary (a cohort that is empty for a given hospital-pathogen pair yields NA draws, which must not be silently counted as valid support).

Eligibility: eligible_for_profile_inference is a hard boolean, TRUE only for rows generated by a conditional (observed-plus-imputed) method – "conditional_analytic_identity" or "conditional_gibbs_correlated"; see compute_event_profile_probabilities(). FALSE for any other profile_generation_method value, including the legacy "unconditional_simulation_correlated_not_daly_eligible" tag from an older profile_output – callers must not rely on the descriptive profile_generation_method string alone. eligible_for_YLL/eligible_for_YLD additionally require: sampler_acceptable (passed in by the caller, e.g. from fitted_model$diagnostics$converged_structural); the relevant event cohort has at least min_n_events events; and at least min_n_draws valid posterior draws contributed to the relevant summary. exclusion_reason_YLL/exclusion_reason_YLD record why, distinguishing a cohort that is completely empty ("no_known_outcome_events" / "no_nonfatal_events") from one that is merely small ("too_few_known_outcome_events" / "too_few_nonfatal_events").

Value

Tibble with one row per hospital x pathogen x profile: R_ALL, R_KNOWN_OUTCOME, and R_NF posterior mean and credible interval, panel composition (n_profile_classes, panel_eligibility_method, classes_excluded, classes_excluded_reason), event and draw counts, and profile/YLL/YLD eligibility flags with reasons.

Description

Creates a configuration object that controls the behavior of the preprocessing pipeline. Users can override default settings for their specific use case.

Usage

amr_config(
  column_mappings = NULL,
  fuzzy_match = TRUE,
  strict_validation = FALSE,
  date_columns = NULL,
  hai_cutoff = 2,
  infer_department = TRUE,
  event_gap_days = 14,
  mortality_window = 14,
  age_bins = "GBD_standard",
  contaminant_method = "auto",
  mdr_definition = "CDC",
  intermediate_as_resistant = TRUE,
  verbose = TRUE
)

Arguments

Value

An amr_config object (list with class)

Examples

## Not run: 
# Use defaults
config <- amr_config()

# Customize for specific hospital
config <- amr_config(
  hai_cutoff = 3,
  mdr_definition = 5,
  strict_validation = FALSE
)

## End(Not run)

Description

Repeats the existing conditional truncated-MVN Gibbs profile completion for fixed posterior draws and a fixed, reproducibly selected subset of incomplete events. It does not refit Stan or alter the fitted model. The comparison therefore isolates finite Gibbs-chain length. The conditional completion algorithm is the correct profile-completion algorithm; this diagnostic checks whether its finite Monte Carlo run length is adequate for a particular fit, especially when posterior residual correlations are high.

Usage

assess_gibbs_profile_stability(
  fitted_model,
  schedules = list(baseline = list(burnin = 30L, kept = 50L), medium = list(burnin = 60L,
    kept = 100L), long = list(burnin = 100L, kept = 250L)),
  n_posterior_draws = 200L,
  max_events = 250L,
  seed = 123L,
  tolerance = 0.01,
  event_selection = c("stratified", "all")
)

Arguments

Details

The default 0.01/0.02 thresholds are numerical-stability conventions for this diagnostic, not universal statistical thresholds.

Value

A structured list with summary, event/aggregate comparisons, selected events, schedules, posterior draw indices, and a separate stability status.

Description

Quantifies uncertainty in resistance-profile probability estimates by repeatedly resampling isolate counts from a binomial distribution implied by the observed n_tested and n_resistant values, refitting the convex optimisation QP for each replicate, and returning percentile confidence intervals across replicates.

Usage

bootstrap_profiles_convex(
  marginals,
  coresistance_output = NULL,
  B = 500L,
  seed = 123L,
  alpha = 0.05,
  n_cores = 1L,
  exclude_near_zero = TRUE,
  top_n_classes = NULL,
  sigma_sq = 1,
  ridge = 1e-08,
  pathogen_col = "organism_name",
  class_col = "antibiotic_class",
  n_tested_col = "n_tested",
  n_resistant_col = "n_resistant",
  org_group_col = "org_group"
)

Arguments

Details

Resampling mechanism: For each bootstrap replicate $b$ and each (pathogen, class) cell:

$n_{\text{resistant}}^{(b)} \sim \text{Binomial}(n_{\text{tested}},\; \hat{r}_{kd})$

The bootstrap marginal is then $\hat{r}_{kd}^{(b)} = n_{\text{resistant}}^{(b)} / n_{\text{tested}}$. Pairwise co-resistance rates (if supplied) are resampled analogously. The QP is re-solved on each bootstrap marginal using the existing compute_resistance_profiles() engine.

Performance: For $n \leq 12$ classes the QP solves in milliseconds; B = 500 replicates will complete in seconds. For $n \geq 14$ classes use n_cores > 1 or reduce B.

Value

A named list, one entry per pathogen. Each entry is a tibble with columns:

profile

Profile label (e.g. "RSS").

probability_mean

Mean profile probability across B replicates.

probability_median

Median across replicates.

lower

Lower percentile (alpha / 2).

upper

Upper percentile (1 - alpha / 2).

n_replicates_converged

Number of replicates for which the QP converged (non-uniform solution).

convergence_rate

Proportion of replicates that converged.

The point-estimate profiles tibble (from a single run on the original marginals) is stored as the attribute "point_estimate".

Examples

## Not run: 
marg <- compute_marginal_resistance(amr_clean)
co_res <- compute_pairwise_coresistance(marg)

boot <- bootstrap_profiles_convex(
  marginals = marg$marginal,
  coresistance_output = co_res,
  B = 500,
  seed = 42
)

# 95% intervals for K. pneumoniae
boot[["Klebsiella pneumoniae"]]

# Point estimate
attr(boot[["Klebsiella pneumoniae"]], "point_estimate")

## End(Not run)

Description

Constructs the constraint matrix M and target vector v used in the simplex-constrained weighted least-squares QP:

Usage

build_constraint_matrix(profiles_enum, r_marg, co_mat = NULL)

Arguments

Details

$\min_p \|Mp - v\|^2_W \quad \text{s.t.} \quad p \ge 0,\; \sum p = 1$

Marginal rows: $M_{d,\delta} = 1$ iff class $d$ is resistant in profile $\delta$. Pairwise rows: $M_{d_1 d_2,\delta} = 1$ iff both classes are resistant. When a pairwise value is unavailable (too few co-tested isolates), the product of marginals (independence assumption) is substituted. Pairwise values that exceed $\min(P(A), P(B))$ are capped.

Value

Named list:

M

Numeric matrix $(n + n(n-1)/2) \times 2^n$.

v

Numeric vector of constraint targets.

constraint_names

Character vector labelling each row of M.

fallback_pairs

Character vector of class pairs that used the independence fallback.

capped_pairs

Named numeric vector of pairs whose pairwise value was capped, showing original and capped values.

Description

Validates that the profile probability distributions returned by compute_resistance_profiles() satisfy all mathematical requirements: non-negativity, sum-to-one, and reconstruction of the input marginal and pairwise resistance rates within a specified tolerance.

Usage

check_profile_constraints(
  profiles_output,
  marginals = NULL,
  pairwise = NULL,
  tolerance = 1e-06,
  pathogen_col = "pathogen",
  class_col = "antibiotic_class",
  rate_col = "marginal_resistance",
  class1_col = "antibiotic_class_1",
  class2_col = "antibiotic_class_2",
  pairwise_rate_col = "pairwise_resistance_prevalence"
)

Arguments

Details

The function uses two complementary sources of evidence:

1. Constraint residuals already stored in profiles_output[[pathogen]]$constraint_residuals – these are $M \hat{p} - v$ computed at solve time and are available without re-supplying the original rates.

2. Direct reconstruction from the binary class-indicator columns in the profiles data frame when marginals and/or pairwise are supplied – independently verifies the residuals and catches any post-solve modifications to the probability vector.

Value

A tibble with one row per constraint per pathogen:

pathogen

Pathogen name.

constraint_type

One of "nonneg", "sum_to_one", "marginal", "pairwise".

constraint_name

Specific label (e.g. "marg_3GC", "pair_3GC_Carbapenems").

target

The constraint target value (NA for nonneg).

reconstructed

The value implied by $\hat{p}$.

abs_residual

Absolute difference between target and reconstructed.

pass

Logical. TRUE if abs_residual < tolerance or the structural check passed.

source

Either "stored_residuals" (from solve-time values) or "recomputed" (from direct reconstruction).

Also prints a per-pathogen summary of passed and failed checks.

Examples

## Not run: 
marg <- compute_marginal_resistance(amr_clean)
co_res <- compute_pairwise_coresistance(marg)
rp <- compute_resistance_profiles(marg, co_res)

# Using stored constraint residuals only
checks <- check_profile_constraints(rp)

# Full re-verification with original rates
checks <- check_profile_constraints(
  rp,
  marginals = marg$marginal,
  pairwise  = NULL,
  tolerance = 1e-5
)

# Inspect failures
checks[!checks$pass, ]

## End(Not run)

Description

Looks up inst/extdata/class_display_labels.csv; falls back to replacing underscores with spaces for any class not in the table (e.g. a pathogen-specific class not yet added there), so this never errors on an unmapped class – it just degrades to a readable-if-generic label.

Usage

class_display_label(x)

Arguments

Value

Character vector of display labels, same length as x.

Description

Human-readable class-pair label

Usage

class_pair_label(class_1, class_2, short = TRUE)

Arguments

Value

Character vector, e.g. "AMG x CARB".

Description

Same lookup/fallback behaviour as class_display_label(), but returns the abbreviated form (e.g. "3GC"). Falls back to the first 4 characters of the underscore-cleaned name, upper-cased, for any class not in the reference table.

Usage

class_short_label(x)

Arguments

Value

Character vector of short labels, same length as x.

Description

Fits Weibull, Lognormal, and Gamma distributions to a numeric vector and returns their AIC values for comparison.

Usage

compare_distribution_aic(x, fits = NULL)

Arguments

Value

Data frame with columns Weibull_AIC, Lognormal_AIC, Gamma_AIC.

Examples

compare_distribution_aic(rlnorm(200, meanlog = 2, sdlog = 0.5))

Description

For each event, retains every observed (tested) AST cell exactly and computes a probability distribution only over the genuinely untested (NA) antibiotic-class cells in that event's profile panel. Any enumerated profile inconsistent with the event's observed cells receives probability 0 – this function never overwrites an observed R/S result.

Usage

compute_event_profile_probabilities(
  fitted_model,
  n_posterior_draws_for_profiles = 2000L,
  outcome_col = NULL,
  nonfatal_values = c("Discharged", "Survived", "Alive", "discharged", "survived",
    "alive"),
  seed = 123L,
  n_gibbs_burnin = 10L,
  n_gibbs_kept = 20L,
  posterior_draw_indices = NULL,
  event_indices = NULL
)

Arguments

Details

Identity residual structure (fitted_model$residual_structure == "identity"): classes are conditionally independent given the linear predictor, so the missing-cell probabilities are computed analytically as $P(Y_{ed}=1 \mid \theta) = \Phi(\mu_{ed})$ for each missing class $d$, and the profile probability over the missing dimensions is the exact product of per-class Bernoulli probabilities – no latent-variable simulation is used, so there is no added latent-profile-simulation noise. Rows carry profile_generation_method = "conditional_analytic_identity".

Correlated residual structure: the missing latent dimensions are sampled conditional on the observed resistance SIGNS of the tested classes via a Gibbs sampler on the truncated multivariate normal – see .gibbs_conditional_profile_probs(). This replaces the earlier unconditional latent-$Z$ simulation, which resampled tested cells along with untested ones and was therefore not observed-plus-imputed. Rows carry profile_generation_method = "conditional_gibbs_correlated" and ARE eligible for downstream DALY use (subject to the same panel-support and sampler-acceptability gates as identity-residual profiles – see aggregate_profiles_for_daly()). Because Gibbs is only run for n_gibbs_burnin + n_gibbs_kept iterations per draw rather than to full convergence, correlated-residual profiles carry additional Monte Carlo error beyond identity-residual profiles at the same n_posterior_draws_for_profiles; increase n_gibbs_kept (and/or n_gibbs_burnin) for lower-noise estimates at higher compute cost.

Both residual structures produce, per event per posterior draw, a profile probability vector over the panel's $2^{D_p}$ enumerated profiles that (a) sums to exactly 1 and (b) assigns exactly 0 to any profile inconsistent with that event's observed cells – enforced by construction in both cases (an exact product for identity; Gibbs never visiting an inconsistent pattern for correlated), and checked with a hard stop() for the identity path since it has no other source of Monte Carlo noise to blur a genuine bug.

Class panels: the antibiotic-class panel enumerated for each hospital x pathogen pair is drawn from the approved eligibility rules computed at fit time (fitted_model$eligibility_report: marginal n_tested/n_resistant/n_susceptible thresholds, plus pairwise co-testing sufficiency for correlated-residual fits) – not from "tested at least once". See .resolve_profile_class_panel().

Estimand: The posterior distribution over the observed event case-mix, conditional on each event's own observed AST results. This is labelled "observed_stewardship_event_mix".

Value

Named list: event_profiles (event-level posterior mean observed-plus-imputed profile probabilities, with n_classes_observed/n_classes_missing per event) and aggregate_draws (per-draw R_ALL [truly all events in the hp pair], R_KNOWN_OUTCOME [known-outcome cohort], and R_NF [nonfatal cohort], used by aggregate_profiles_for_daly() for credible intervals). Both tibbles contain all $2^D$ profiles per hospital-pathogen pair and carry profile_generation_method and panel_reason columns.

Description

Distributes pooled YLL/YLD totals to individual hospitals proportionally by deaths and discharges, and computes per-1 000 rates.

Usage

compute_hospital_daly(
  hospital_counts,
  total_deaths,
  total_discharged,
  yll_base,
  yll_associated,
  yll_attributable,
  yld_base,
  yld_associated,
  yld_attributable
)

Arguments

Value

Data frame with hospital-level YLL, YLD, DALY (absolute and per 1 000 cases) columns.

Description

Collapses drug-level susceptibility data to antibiotic-class level per isolate ($R_{e,k,c} = 1$ if resistant to any drug in class $c$), then computes marginal resistance for every pathogen x class combination found in the data.

Usage

compute_marginal_resistance(
  data,
  pathogen_col = "organism_name",
  org_group_col = "org_group",
  isolate_col = "isolate_id",
  antibiotic_class_col = "antibiotic_class",
  antibiotic_value_col = "antibiotic_value",
  zero_threshold = 0,
  min_n_tested = 30,
  facility_col = NULL,
  facility_name = NULL,
  outcome_col = NULL,
  outcome_value = NULL
)

Arguments

Details

Classes whose marginal resistance is at or below zero_threshold are listed in $near_zero as a flag for downstream use – they are not removed here.

Value

Named list:

marginal

Data frame with columns: pathogen_col, org_group_col, antibiotic_class_col, n_tested, n_resistant, marginal_resistance. Sorted descending by marginal_resistance within each pathogen.

near_zero

Subset of marginal where marginal_resistance <= zero_threshold. These classes are candidates for exclusion in downstream profiling.

class_long

Collapsed isolate x pathogen x class data frame (columns: isolate_col, pathogen_col, org_group_col, antibiotic_class_col, class_result). Pass this directly to compute_pairwise_coresistance().

Examples

## Not run: 
marg <- compute_marginal_resistance(
  data                 = amr_clean,
  pathogen_col         = "organism_name",
  org_group_col        = "org_group",
  isolate_col          = "isolate_id",
  antibiotic_class_col = "antibiotic_class",
  antibiotic_value_col = "antibiotic_value"
)

marg$marginal # full marginal resistance table
marg$near_zero # classes flagged as near-zero

## End(Not run)

Description

Takes the wide-format tibble produced by preprocess_for_profiles() and computes marginal resistance rates per pathogen x antibiotic class, optionally stratified by geography, year, and/or patient outcome. The counting unit is always isolate_id: a patient with multiple isolates contributes one count per isolate.

Usage

compute_marginals_from_data(
  data_wide,
  col_map,
  panel_map,
  stratify_by = NULL,
  outcome_col = NULL,
  min_n_tested = 30L,
  external_marginals = NULL,
  ext_col_map = list(pathogen_col = "pathogen", class_col = "antibiotic_class",
    geography_col = "geography", year_col = "year", rate_col = "resistance_prevalence")
)

Arguments

Details

Optionally, externally modelled marginals (e.g. from GBD ST-GPR) can be supplied to override the locally computed values while retaining the pairwise correlation structure from the local data.

Value

Tibble with columns: pathogen, antibiotic_class, [stratum columns], n_tested, n_resistant, marginal_resistance, marginal_source ("computed" or "external").

Description

For every pathogen and every pair of antibiotic classes $(c_i, c_j)$ that were both tested:

$T_{k,i,j} = \sum_e \mathbf{1}(c_i \text{ tested} \land c_j \text{ tested})$

$R_{k,i,j} = \sum_e \mathbf{1}(R_{e,k,c_i}=1 \land R_{e,k,c_j}=1)$

$\text{Prev}_{k,i,j} = R_{k,i,j} \;/\; T_{k,i,j}$

Usage

compute_pairwise_coresistance(
  marginal_output,
  pathogen_col = "organism_name",
  org_group_col = "org_group",
  isolate_col = "isolate_id",
  antibiotic_class_col = "antibiotic_class",
  min_co_tested = 10,
  facility_col = NULL,
  facility_name = NULL,
  outcome_col = NULL,
  outcome_value = NULL
)

Arguments

Details

Matrices are computed for all tested classes – no filtering by marginal resistance or GBD core list is applied here.

Value

Named list. One entry per pathogen (keyed by pathogen name), each containing:

prevalence

Square symmetric matrix of pairwise co-resistance rates. Diagonal and cells with n < min_co_tested are NA.

T_matrix

Integer matrix of co-tested isolate counts.

R_matrix

Integer matrix of co-resistant isolate counts.

classes

Character vector of class names (row/column order).

Examples

## Not run: 
marg <- compute_marginal_resistance(amr_clean, ...)
co_res <- compute_pairwise_coresistance(marg)

# Prevalence matrix for K. pneumoniae
co_res[["Klebsiella pneumoniae"]]$prevalence

# Co-tested counts
co_res[["Klebsiella pneumoniae"]]$T_matrix

## End(Not run)

Description

For each antibiotic-class pair (A, B) and each pathogen [x stratum], this function:

1. Computes the Pearson correlation $\rho_{AB}$ between binary resistance indicators across co-tested isolates.

2. Back-calculates pairwise prevalence using the GBD formula:

   $P(A \cap B) = P(A) P(B) + \rho_{AB} \sqrt{P(A)(1-P(A)) \cdot P(B)(1-P(B))}$

   where $P(A)$ and $P(B)$ come from marginals (which may include externally overridden values).

3. Caps the result to $\min(P(A), P(B))$ and floors it to 0. Cells below the floor use the independence product as fallback.

Usage

compute_pairwise_from_data(
  data_wide,
  marginals,
  col_map,
  panel_map,
  stratify_by = NULL,
  outcome_col = NULL,
  min_co_tested = 10L
)

Arguments

Value

Tibble with one row per pathogen x class-pair [x stratum]:

pathogen, antibiotic_class_1, antibiotic_class_2

Keys.

n_co_tested

Isolates tested for both classes.

rho

Pearson correlation between binary resistance indicators.

pairwise_prevalence

Back-calculated P(AnB), capped and floored.

method

"pearson_back_calc" or "independence_fallback".

Description

Does NOT classify from a single tail probability. Groups supported statistics into families (marginal, pairwise, profile [complete-profile and resistant-count statistics – both address joint multidrug-resistance burden], hospital_heterogeneity, cluster [admission and patient clustering]) and flags a family when either (a) any statistic in it is SEVERE (ppc_tail_probability < tail_severe or > 1 - tail_severe), or (b) the FRACTION of that family's supported statistics that are extreme exceeds max_fraction_core_extreme. Thresholds are configurable, not universal hardcoded truths, and are applied as discrepancy FLAGS, not frequentist rejection tests.

Usage

compute_posterior_predictive_status(ppc_statistics, thresholds = list())

Arguments

Details

Statistics with essentially ZERO posterior-predictive variance (replicated_sd at or near 0 – e.g. a discrete statistic that lands on the same value in every replicated state, such as max(resistant_count) pinned at its ceiling for a small class panel) are excluded from extreme/severe classification (though still returned in ppc_statistics with their true, unmodified ppc_tail_probability): mean(T_rep >= T_obs) is exactly 1 whenever every replicate exactly EQUALS the observed value, which is a well-known degeneracy of posterior-predictive tail probabilities for near-constant discrete statistics, not evidence of misfit – a discrepancy statistic that cannot vary under the posterior predictive distribution carries no discriminating information about model fit.

Value

List with status (one of "pass", "warning_marginal_ppc", "warning_pairwise_ppc", "warning_profile_ppc", "warning_hospital_heterogeneity_ppc", "warning_cluster_ppc", "fail_major_ppc_misfit", "insufficient_ppc_support"), reasons (character vector of every triggered family-level flag), thresholds_used, and family_status (per-family n/n_extreme/n_severe/fraction_extreme).

Description

Two experiments that share the same prior-generative structure (same residual structure, same beta/tau/LKJ priors, same declared random-effect block structure, same fixed-effect design-column definitions, same D) can safely REUSE a prior predictive check's results rather than rerunning simulate_probit_prior_predictive for every experiment – prior predictive checking depends only on the generative structure, never on the fitted posterior. If ANY component changes (a prior, a scaling, the declared block structure), the fingerprint changes and any cached result must be invalidated.

Usage

compute_prior_predictive_fingerprint(fitted_model)

Arguments

Value

List with fingerprint (a single character hash via rlang::hash(), already an anumaan dependency – no new package dependency introduced) and components (the exact named list that was hashed, for audit/debugging).

Description

Unlike posterior predictive checking, the goal here is NOT to make the simulated distribution match the observed data closely – it is to flag obviously implausible/extreme implications of the CHOSEN priors, before conditioning on any outcomes. Computes, per prior state and averaged across states: the fraction of event-class probabilities $\Phi(\mu_{ed})$ below 0.001 or above 0.999 (near-deterministic implied resistance), the fraction of generated events that are degenerate (all classes resistant or all susceptible), and the spread of per-"hospital" (the model's primary declared random-effect grouping) mean implied probability – a very large spread implies an implausibly extreme facility-to-facility prior.

Usage

compute_prior_predictive_status(prior_draws, thresholds = list())

Arguments

Value

List with status (one of "pass", "warning_extreme_prior_predictions", "warning_excessive_hospital_heterogeneity", "warning_degenerate_profiles", "fail_implausible_prior_predictive", "insufficient_prior_check"), reasons, thresholds_used, and summary (the averaged plausibility fractions).

Description

Drives ppc_draws$generate_state(s) once per posterior state (see simulate_probit_posterior_predictive), streaming each state's complete and observation-masked replicate through every requested discrepancy statistic and accumulating only the compact per-state scalar contributions – never an S x N_events x D array. Statistic families (Stan User's Guide-style posterior predictive checks, adapted to AMR resistance data):

marginal

Per-class resistance rate, pooled and per eligible hospital x pathogen x class cell (reusing fitted_model$eligibility_report$marginal).

resistant_count

Number of resistant classes per event, among whatever was tested for that event, summarised (mean/variance/median/ p90/max); plus, for fully-observed eligible panels, the full count distribution.

pairwise

RR/RS/SR/SS co-resistance proportions for every sufficiently co-tested class pair (reusing fitted_model$eligibility_report$pairwise).

complete_profile

All-susceptible/all-resistant/most-common-profile frequency, distinct-profile count, Shannon entropy, Simpson concentration, and the top-3 most frequent OBSERVED profiles (profile_top_observed_frequency – a purely data-dominant descriptive summary, deliberately NOT a claim of clinical importance; see the file header of probit_posterior_predictive.R for the distinction and the deferred clinical_priority_profiles extension point), restricted to adequately-supported fully-observed hospital x pathogen panels (enumerate_binary_profiles- style class panel via .resolve_profile_class_panel()).

hospital_heterogeneity

SD/IQR/MAD/range of per-hospital resistance proportions among hospitals meeting min_hospital_support, per class.

admission_clustering, patient_clustering

Within-/between-group same-class agreement and Hamming distance, only when the relevant role is declared eligible (see random_effect_roles and random_effect_eligibility_audit below); skipped cleanly (support_status != "supported"), never as a failure, when unsupported.

Usage

compute_probit_ppc_statistics(
  fitted_model_or_ppc_draws,
  n_states = 500L,
  seed = 123L,
  preserve_observation_mask = TRUE,
  statistics = c("marginal", "resistant_count", "pairwise", "complete_profile",
    "hospital_heterogeneity", "admission_clustering", "patient_clustering"),
  random_effect_roles = NULL,
  random_effect_eligibility_audit = NULL,
  min_hospital_support = 30L,
  min_complete_profile_events = 30L,
  ci_level = 0.95
)

Arguments

Value

A tibble with columns statistic_name, stratum, observed_value, replicated_mean, replicated_sd, replicated_q025, replicated_q50, replicated_q975, ppc_tail_probability, ppc_two_sided, n_replications, support_status.

References

Stan Development Team. "Posterior and Prior Predictive Checks." Stan User's Guide. https://mc-stan.org/docs/stan-users-guide/posterior-predictive-checks.html

Description

Calculates performance and processing statistics.

Usage

compute_processing_statistics(raw_data, processed_data, preprocessing_log)

Arguments

Value

List with statistics

Description

Combines the outputs of validate_marginal_calibration(), validate_pairwise_calibration(), validate_complete_profile_calibration(), and mask_and_validate_ast() into a single profile_validation_status field. This is deliberately separate from diagnostic_status (see fit_bayesian_multivariate_probit()): diagnostic_status describes HMC sampler/model-fitting health (R-hat, ESS, divergences, E-BFMI, tree-depth); profile_validation_status describes predictive agreement with the observed data and data support for that agreement. A model can converge cleanly (diagnostic_status = "pass") while its profile predictions are poorly calibrated, or vice versa, and the two should never be merged into one flag.

Usage

compute_profile_validation_status(
  marginal_tbl = NULL,
  pairwise_tbl = NULL,
  complete_profile_tbl = NULL,
  masked_summary = NULL,
  thresholds = list()
)

Arguments

Details

All thresholds are explicit parameters with documented defaults – none are silently hard-coded elsewhere. Pass thresholds to override any subset; unspecified entries keep their default.

In-sample masked validation can never produce a bare "pass" on its own. When masked_summary$validation_mode == "in_sample_no_refit" (see mask_and_validate_ast()), the masked cells were already seen during fitting, so predicting them well is not genuine validation. If every threshold is otherwise cleared, the status becomes "in_sample_validation_only" rather than "pass" – a real "pass" requires either validation_mode == "refit_masked" (a genuine holdout) or no masked-AST check at all (the other observed-versus-model checks being used purely descriptively).

Value

Named list: status (one of "pass", "in_sample_validation_only", "warning_marginal_miscalibration", "warning_pairwise_miscalibration", "warning_sparse_complete_profiles", "fail_masked_ast_validation", or "not_evaluated" when no inputs were supplied), reasons (character vector of every triggered condition – including "masked_validation_in_sample_only" as an informational reason even when a more severe status wins – the most severe status is returned but all are listed), and thresholds_used (the resolved threshold list, for audit trail).

Description

For each pathogen $k$ with $n_k$ antibiotic classes, enumerates all $2^{n_k}$ binary resistance profiles $\delta \in \{0,1\}^{n_k}$ and estimates their probabilities by solving a simplex-constrained weighted least-squares Quadratic Programme (GBD equation 7.5.1.3):

Usage

compute_resistance_profiles(
  marginal_output,
  coresistance_output,
  pathogens = NULL,
  top_n_pathogens = NULL,
  exclude_near_zero = TRUE,
  top_n_classes = NULL,
  sigma_sq = 1,
  ridge = 1e-08,
  pathogen_col = "organism_name",
  antibiotic_class_col = "antibiotic_class",
  facility_col = NULL,
  facility_name = NULL,
  outcome_col = NULL,
  outcome_value = NULL,
  n_cores = 1L
)

Arguments

Details

$\hat{p} = \arg\min_{p \in \Delta_{2^n}} \sum_{i=1}^{m} \frac{(m_i^\top p - v_i)^2}{\sigma_i^2}$

where $m = n(n+1)/2$ data-derived linear constraints encode n marginal resistance rates and n(n-1)/2 pairwise co-resistance rates, and $\Delta$ is the standard probability simplex.

Constraint rows in M

Marginal (rows 1 to n)

Row $d$: $M_{d,\delta} = 1$ iff $\delta_d = 1$ (i.e., class $d$ is resistant in profile $\delta$). Constraint: $\sum_\delta M_{d,\delta}\,p_\delta = \hat{r}_{kd}$.

Pairwise (rows n+1 to m)

Row $(d_1,d_2)$: $M_{d_1 d_2,\delta} = 1$ iff $\delta_{d_1} = 1 \land \delta_{d_2} = 1$. Constraint: $\sum_\delta M_{d_1 d_2,\delta}\,p_\delta = \hat{r}_{k,d_1 d_2}$. When a pairwise estimate is unavailable (too few co-tested isolates), the product of marginals (independence assumption) is used as fallback.

The QP is solved via quadprog::solve.QP. A small ridge term (ridge) is added to the Hessian to guarantee strict positive-definiteness. On solver failure the pathogen gets a uniform distribution over all profiles.

Performance Notes

This function has been optimized for speed with vectorized profile generation and label creation (10-100x faster than previous versions). However, computational complexity is still exponential in the number of classes:

• n <= 14: Fast (seconds to minutes)

• n = 15: Moderate (minutes)

• n >= 16: Slow and memory-intensive (use top_n_classes)

For large datasets with many pathogens, consider using top_n_classes to limit each pathogen to its most-tested classes (e.g., top_n_classes = 12).

Value

Named list, one entry per pathogen, each a list with:

profiles

Data frame: profile (character label, e.g.\ "RSR"), probability ($\hat{p}_\delta$), and one binary (0/1) integer column per antibiotic class indicating whether that class is resistant (1) or susceptible (0) in each profile.

classes

Character vector of antibiotic class names used (alphabetical; bit 0 = classes[1]).

n_classes

Integer. Number of classes used.

constraint_residuals

Named numeric vector of $m_i^\top \hat{p} - v_i$ for each constraint. Small absolute values indicate good constraint satisfaction.

Examples

## Not run: 
marg <- compute_marginal_resistance(amr_clean)
co_res <- compute_pairwise_coresistance(marg)
rp <- compute_resistance_profiles(marg, co_res)

# All profiles and probabilities for K. pneumoniae
rp[["Klebsiella pneumoniae"]]$profiles

# Constraint residuals (quality check)
rp[["Klebsiella pneumoniae"]]$constraint_residuals

# Single pathogen
rp_kp <- compute_resistance_profiles(
  marg, co_res,
  pathogens = "Klebsiella pneumoniae"
)

## End(Not run)

Description

Maps organism names to GBD RR pathogen categories and antibiotic class names to RR drug categories in a single pass. Either mapping is skipped silently when the corresponding column is absent from data.

Usage

daly_add_rr_mappings(
  data,
  organism_col = "organism_normalized",
  class_col = "antibiotic_class"
)

Arguments

Details

Pathogen mapping uses the organism taxonomy (get_organism_taxonomy()); drug class mapping uses inst/extdata/WHO_aware_class.csv.

Value

Data frame with rr_pathogen and/or rr_drug columns appended.

Description

For each resistance profile delta (from compute_resistance_profiles()), determines the profile-level RR_kd_LOS using the GBD max rule: RR_kd_LOS = max over c in C_R(d) of RR_kc_LOS [if C_R(d) non-empty] = 1 [if d = all-susceptible] where C_R(d) = {c : d_c = 1}. The CI reported for each profile is that of its dominant (max-RR) class.

Usage

daly_assign_rr_to_profiles(
  profiles_output,
  rr_table,
  pathogen_col = "pathogen",
  class_col = "antibiotic_class",
  rr_col = "RR_LOS",
  fallback_rr = 1,
  los_r_col = NULL,
  los_s_col = NULL,
  los_unit = c("years", "days")
)

Arguments

Details

If rr_table was fit with a syndrome_name filter, its RR values are syndrome-specific but get applied here to all profiles of a pathogen – i.e. it assumes syndrome-invariant LOS prolongation. Refit with syndrome_name = NULL for a RR pooled across syndromes.

Optionally, when los_r_col and los_s_col are supplied, also propagates the absolute mean LOS for the profile's dominant class (resistant and susceptible) as LOSR_years / LOSS_years. This does not change the RR assignment above; it is an additive lookup against the dominant class already selected by the max rule, using whichever absolute-LOS columns are present in rr_table (e.g. mean_LOS_R / mean_LOS_S from daly_fit_los_rr_distribution(), or mean_los_resistant / mean_los_susceptible from daly_fit_los_rr()). LOSR_years/LOSS_years are set to NA for the all-susceptible profile, which has no dominant class.

Value

Named list (one entry per pathogen): original profiles data frame augmented with RR_LOS_profile, dominant_class, (if available) CI_lower_profile / CI_upper_profile, and (when los_r_col / los_s_col are supplied) LOSR_years / LOSS_years.

Description

Computes the case fatality ratio (CFR) for each pathogen (k) within a specified infectious syndrome (L) using facility-level microbiology data.

Usage

daly_calc_case_fatality(
  data,
  syndrome_col,
  syndrome_name,
  specimen_col,
  specimen_name,
  polymicrobial_col,
  patient_col,
  pathogen_col,
  outcome_col,
  death_value = "Died",
  glass_ref = NULL,
  facility_col = NULL,
  facility_name = NULL,
  pathogen_name = NULL
)

Arguments

Details

The unit of analysis is the patient. Each patient contributes total weight 1, distributed equally across their valid pathogens via a fractional weight $1/m_r$ (where $m_r$ is the number of distinct valid pathogens for patient r). This ensures polymicrobial patients are not double-counted.

For polymicrobial patients (polymicrobial_col == 1), only pathogens listed in the GLASS reference (glass_ref) for the given specimen type are retained before weighting. Monomicrobial patients (polymicrobial_col == 0) are never filtered.

When facility_col is supplied and facility_name is NULL, both per-facility and pooled CFR are returned. The pooled CFR sums weighted deaths and totals across facilities before dividing.

Value

A list:

cfr_table

Pooled CFR: pathogen_col, weighted_deaths, weighted_total, CFR_LK.

facility_level

Per-facility CFR (only when facility_col is supplied and facility_name is NULL).

Description

Computes CR_L, the factor that adjusts a hospital-derived CFR to account for infection cases managed outside the inpatient setting. The adjustment type for each syndrome is looked up from the adjustment_ref table (loaded from ‘inst/extdata/adjustment_for_CFR’).

Usage

daly_calc_cr_l(
  data,
  syndrome_col,
  syndrome_name,
  patient_col,
  visit_type_col,
  inpatient_value = "Inpatient",
  outpatient_value = "Outpatient",
  adjustment_ref,
  facility_col = NULL,
  facility_name = NULL
)

Arguments

Details

Three adjustment types are supported:

None

CR_L = 1. The hospital CFR applies directly (e.g., BSI, Meningitis, hospital-acquired infections).

Inpatient ratio

CR_L = (patients with $\ge$ 1 inpatient visit) / (all patients). Used when community cases are captured partly in outpatient data (e.g., community-acquired LRI, UTI).

Outpatient to inpatient ratio

CR_L = (patients with $\ge$ 1 outpatient AND $\ge$ 1 inpatient visit) / (patients with $\ge$ 1 outpatient visit). Used for syndromes where OP-to-IP transition captures disease severity (e.g., STI, Skin, Eye, Oral, Bone/joint infections).

Value

Data frame with columns syndrome (= syndrome_name), adjustment_type, CR_L, and (when facility_col is supplied) facility_col. Additional columns (n_inpatient, n_total, etc.) give the raw counts used.

Description

Computes D_J, the number of deaths for each underlying cause J from population-level vital registration or mortality data. Optionally stratified by grouping variables such as age group, sex, year, or location.

Usage

daly_calc_deaths_by_cause(
  pop_data,
  cause_col = "cause_of_death",
  deaths_col = NULL,
  groupby_cols = NULL
)

Arguments

Details

This function requires population-level data. If only facility-level data are available, use daly_calc_deaths_by_syndrome() with facility_data instead, which directly counts deaths by syndrome.

Value

Data frame with columns: cause_col, any groupby_cols, D_J (death count), D_J_method ("population"), D_J_confidence ("high").

References

Antimicrobial Resistance Collaborators. Global burden of bacterial antimicrobial resistance in 2019. Lancet. 2022.

Examples

## Not run: 
# One row per death record
d_j <- daly_calc_deaths_by_cause(
  pop_data  = vital_reg,
  cause_col = "icd10_cause"
)

# Pre-aggregated counts, stratified by age and sex
d_j <- daly_calc_deaths_by_cause(
  pop_data     = vital_reg,
  cause_col    = "icd10_cause",
  deaths_col   = "n_deaths",
  groupby_cols = c("Age_bin", "gender")
)

## End(Not run)

Description

Computes D_L, the number of deaths due to each infectious syndrome L.

Usage

daly_calc_deaths_by_syndrome(
  d_j = NULL,
  s_j = NULL,
  m_lj = NULL,
  pop_data = NULL,
  facility_data = NULL,
  cause_col = "cause_of_death",
  syndrome_col = "syndrome",
  syndrome = NULL,
  deaths_col = NULL,
  infection_flag_col = "is_infection_death",
  outcome_col = "final_outcome",
  death_value = "Died",
  patient_col = "patient_id",
  facility_col = NULL,
  facility_name = NULL,
  groupby_cols = NULL
)

Arguments

Details

Three operating modes, evaluated in priority order:

1. Pre-computed components (HIGH confidence): pass the outputs of daly_calc_deaths_by_cause(), daly_calc_infection_fraction(), and daly_calc_syndrome_fraction() via d_j, s_j, and m_lj. Computes: D_L = sum_J(D_J * S_J * M_LJ).

2. Raw population data (HIGH confidence): pass pop_data and all three components are computed internally before combining. Computes: D_L = sum_J(D_J * S_J * M_LJ).

3. Facility fallback (LOW confidence): when no population inputs are available, counts the number of unique patients who died for each syndrome directly from facility_data.

Use syndrome to restrict results to a single syndrome of interest. Use facility_col and facility_name to restrict the facility fallback to a specific site.

Value

Data frame with columns: syndrome_col, any groupby_cols, D_L (deaths by syndrome), D_L_method, D_L_confidence.

References

Antimicrobial Resistance Collaborators. Global burden of bacterial antimicrobial resistance in 2019. Lancet. 2022.

Examples

## Not run: 
# Mode 1: pass pre-computed components, all syndromes
d_j <- daly_calc_deaths_by_cause(pop_data = vr, cause_col = "icd10")
s_j <- daly_calc_infection_fraction(pop_data = vr, cause_col = "icd10")
m_lj <- daly_calc_syndrome_fraction(
  pop_data = vr, cause_col = "icd10", syndrome_col = "syndrome"
)
d_l <- daly_calc_deaths_by_syndrome(
  d_j = d_j, s_j = s_j, m_lj = m_lj,
  cause_col = "icd10", syndrome_col = "syndrome"
)

# Mode 2: raw population data, filter to one syndrome
d_l <- daly_calc_deaths_by_syndrome(
  pop_data     = vital_reg,
  cause_col    = "icd10_cause",
  syndrome_col = "infectious_syndrome",
  syndrome     = "Bloodstream infection"
)

# Mode 3: facility fallback, one facility, all syndromes
d_l <- daly_calc_deaths_by_syndrome(
  facility_data = amr_data,
  syndrome_col  = "specimen_normalized",
  outcome_col   = "final_outcome",
  death_value   = "Died",
  patient_col   = "patient_id",
  facility_col  = "location",
  facility_name = "Mumbai"
)

# Mode 3: facility fallback, all facilities, one syndrome, stratified by age
d_l <- daly_calc_deaths_by_syndrome(
  facility_data = amr_data,
  syndrome_col  = "specimen_normalized",
  syndrome      = "Blood",
  death_value   = "Died",
  groupby_cols  = c("Age_bin")
)

## End(Not run)

Description

Estimates the number of incident cases of syndrome L using:

Usage

daly_calc_incidence_from_cfr(
  deaths_L,
  cfr_lk_tbl,
  P_Lk_prime_tbl,
  CR_L = 1,
  pathogen_col = "pathogen",
  cfr_col = "CFR_LK",
  plk_col = "P_Lk_prime",
  facility_col = NULL,
  deaths_col = "deaths"
)

Arguments

Details

$I_L = \frac{D_L}{\text{CFR}_L \times \text{CR}_L}$

where the syndrome-level CFR is the pathogen-weighted average:

$\text{CFR}_L = \sum_k P'_{Lk} \times \text{CFR}_{Lk}$

Use this when you have population- or facility-level death counts and want to back-calculate incidence. For a direct patient count from facility data, use count_incident_cases() instead.

Value

Data frame with columns deaths, CFR_L, CR_L, I_L (incident cases), and facility_col when applicable.

Description

Computes S_J, the fraction of deaths for underlying cause J that were related to infection. Values range from 0 (no infection involvement) to 1 (all deaths for cause J involved infection).

Usage

daly_calc_infection_fraction(
  pop_data,
  cause_col = "cause_of_death",
  infection_flag_col = "is_infection_death",
  groupby_cols = NULL
)

Arguments

Details

This function requires population-level data where each death record carries a binary flag indicating whether infection was involved. If only facility data are available, use daly_calc_deaths_by_syndrome() with facility_data instead.

Value

Data frame with columns: cause_col, any groupby_cols, D_J (total deaths), infection_deaths (infection-related death count), S_J (infection fraction 0-1), S_J_method, S_J_confidence.

References

Antimicrobial Resistance Collaborators. Global burden of bacterial antimicrobial resistance in 2019. Lancet. 2022.

Examples

## Not run: 
s_j <- daly_calc_infection_fraction(
  pop_data           = vital_reg,
  cause_col          = "icd10_cause",
  infection_flag_col = "is_infectious",
  groupby_cols       = c("Age_bin", "gender")
)

## End(Not run)

Description

Computes the population attributable fraction (PAF) for length of stay from a named list of profile data frames, each containing resistance-profile probabilities and profile-level relative risk of LOS.

Usage

daly_calc_paf_los(
  profiles_with_rr,
  probability_col = "probability",
  rr_profile_col = "RR_LOS_profile",
  profile_col = "profile"
)

Arguments

Details

For each pathogen the formula is: PAF = sum_d R'_kd * (RR_kd - 1) / (1 + sum_d R'_kd * (RR_kd - 1))

where R'_kd is the probability of resistance profile d and RR_kd is the corresponding relative LOS multiplier.

Value

Named list (one entry per pathogen) with columns profile, probability, rr_profile_col, numerator, PAF_LOS, and denominator added to the input profile data frame.

Description

Computes PAF_kd_mortality for each pathogen k and resistance profile delta using the GBD multi-exposure Levin formula, substituting the mortality odds ratio (OR_death) from fit_mortality_rr_logistic() in place of the LOS relative risk used by compute_paf_los():

Usage

daly_calc_paf_mortality(
  profiles_with_rr,
  probability_col = "probability",
  rr_profile_col = "RR_LOS_profile",
  profile_col = "profile"
)

Arguments

Details

$\text{PAF}_{kd,\text{mort}} = \frac{R'_{K\delta}\,(\text{OR}_{K\delta} - 1)} {1 + \sum_\delta R'_{K\delta}\,(\text{OR}_{K\delta} - 1)}$

The all-susceptible profile carries OR = 1 and contributes 0. The denominator equals $E[\text{OR}_k] = \sum_\delta R'_{K\delta} \cdot \text{OR}_{K\delta}$, numerically identical to the denominator produced by compute_paf_los().

Overall mortality PAF for pathogen k:

$\text{PAF}_{k,\text{mort}} = \sum_\delta \text{PAF}_{kd,\text{mort}} = \frac{\sum_\delta R'_{K\delta}(\text{OR}_{K\delta}-1)} {1 + \sum_\delta R'_{K\delta}(\text{OR}_{K\delta}-1)}$

Usage pipeline:

  # 1. Fit mortality OR per class
  mort_or <- fit_mortality_rr_logistic(data, ...)

  # 2. Assign OR to profiles via max rule (rr_col = "OR_death")
  profiles_with_or <- assign_rr_to_profiles(
      profiles_output,
      rr_table = mort_or,
      rr_col   = "OR_death"
  )

  # 3. Compute per-profile and overall mortality PAF
  paf_mort <- daly_calc_paf_mortality(profiles_with_or)

Value

Named list (one entry per pathogen) containing:

• per_profile: profile data frame augmented with numerator_mort (= $R'_{K\delta}(\text{OR}_{K\delta}-1)$), PAF_mortality (= numerator / denominator), and denominator_mort (= $1 + \sum_\delta$ numerator).

• PAF_k_mort: overall mortality PAF for pathogen k.

• denominator_mort: shared denominator $E[\text{OR}_k]$.

Description

Computes the fatal pathogen distribution for a given infectious syndrome (L) using facility-level microbiology data. This quantity represents the fractional contribution of each pathogen (k) to fatal infection cases, and is used in YLL / mortality burden estimation.

Usage

daly_calc_pathogen_fraction_fatal(
  data,
  syndrome_col,
  syndrome_name,
  specimen_col = NULL,
  specimen_name = NULL,
  polymicrobial_col,
  patient_col,
  pathogen_col,
  outcome_col,
  death_value = "Death",
  glass_ref = NULL,
  facility_col = NULL,
  facility_name = NULL,
  pathogen_name = NULL
)

Arguments

Details

The unit of analysis is the patient. Each patient contributes total weight 1, distributed equally across their valid pathogens. For a patient with $m_r$ valid pathogens, each pathogen receives weight $1/m_r$.

For polymicrobial patients (polymicrobial_col == 1), only pathogens listed in the GLASS reference (glass_ref) for the given specimen type are retained before weighting. Monomicrobial patients (polymicrobial_col == 0) are never filtered.

The pooled formula across facilities is:

$P_{LK}^{\text{pooled}} = \frac{\sum_f N^{F}_{f,L,K}}{\sum_f N^{F}_{f,L}}$

Value

A list:

P_Lk_fatal

Pooled P_LK data frame: pathogen_col, N_F_LK, N_F_L, P_Lk_fatal.

facility_level

Per-facility P_LK (only when facility_col is supplied and facility_name is NULL).

Description

Computes the profile-weighted expected proportion of deaths attributable to resistance (R_k) and the expected odds ratio of death (E[OR_death]) for each pathogen.

Usage

daly_calc_resistance_prevalence_fatal(
  profiles_with_rr,
  probability_col = "probability",
  rr_profile_col = "RR_LOS_profile"
)

Arguments

Value

Named list (one entry per pathogen) each containing per_profile data frame, R_k, and E_OR_k.

Description

Computes M_LJ, the fraction of infection-related deaths for underlying cause J that are attributed to infectious syndrome L. This distributes infection deaths across clinical syndrome categories (e.g., bloodstream infection, pneumonia, urinary tract infection).

Usage

daly_calc_syndrome_fraction(
  pop_data,
  cause_col = "cause_of_death",
  syndrome_col = "syndrome",
  infection_flag_col = "is_infection_death",
  groupby_cols = NULL
)

Arguments

Details

This function requires population-level data. If only facility data are available, use daly_calc_deaths_by_syndrome() with facility_data instead.

Value

Data frame with columns: cause_col, syndrome_col, any groupby_cols, infection_deaths_LJ (deaths for syndrome L given cause J), infection_deaths_J (total infection deaths for cause J), M_LJ (syndrome fraction 0-1), M_LJ_method, M_LJ_confidence.

References

Antimicrobial Resistance Collaborators. Global burden of bacterial antimicrobial resistance in 2019. Lancet. 2022.

Examples

## Not run: 
m_lj <- daly_calc_syndrome_fraction(
  pop_data           = vital_reg,
  cause_col          = "icd10_cause",
  syndrome_col       = "infectious_syndrome",
  infection_flag_col = "is_infectious",
  groupby_cols       = c("Age_bin")
)

## End(Not run)

Description

Primary, direct implementation of the profile-specific associated-YLD equation:

Usage

daly_calc_yld_associated(
  profiles_with_los,
  P_Lk_prime_tbl,
  yld_ref = NULL,
  DW_sepsis = NULL,
  pathogen_col = "pathogen",
  plk_col = "P_Lk_prime",
  incidence_col = "N_NF_L",
  facility_col = NULL,
  facility_name = NULL,
  facility_state_map = NULL,
  state_col = "state",
  state_name = NULL,
  probability_col = "probability",
  dominant_class_col = "dominant_class",
  losr_col = "LOSR_years"
)

Arguments

Details

YLD_assoc_k = I_L * P'_Lk * DW_L * sum_{delta in D_k} [R'_k_delta * LOSR_k,d*(delta)]

where D_k is the set of RESISTANT resistance profiles for pathogen k (the all-susceptible profile has no dominant class and is excluded – see "Details"). d*(delta) is the dominant drug class of profile delta (GBD max rule, assigned by daly_assign_rr_to_profiles()).

This does not require computing a baseline/pooled YLD first: unlike the pre-existing multiplicative pathway (a pooled baseline YLD multiplied by an associated-burden fraction), this is the direct equation, evaluated per pathogen (and per facility, if P_Lk_prime_tbl is facility-level) from its own inputs.

Component sourcing:

• I_L and P'_Lk: read directly from P_Lk_prime_tbl – the P_Lk_prime (pooled) or facility_level element already returned by daly_calc_pathogen_fraction_nonfatal(), which computes both the non-fatal incidence count and the pathogen fraction from the same observed cohort. Incidence is never back-calculated from deaths, CFR, or LOS.

• DW_L: resolved exactly as in the package's pre-existing YLD weight mechanism – either the DW_sepsis scalar, or a lookup in yld_ref (loaded from inst/extdata/Proxy_YLD_per_case.csv) by state/facility. See "Disability weight caveat" below.

• R'_k_delta, d*(delta), and LOSR: read from profiles_with_los, the list returned by daly_assign_rr_to_profiles() when called with los_r_col/los_s_col set.

Why D_k excludes the all-susceptible profile: this matches the pre-existing convention elsewhere in this package that "associated" burden is a partition over resistant profiles only (the removed daly_calc_fraction_associated_yld()'s Fraction_k summed over resistant profiles only; daly_calc_paf_los()'s all-susceptible term already contributes exactly zero by construction, since RR=1 there). YLD_associated therefore represents the burden occurring specifically among pathogen-k patients with a resistant profile, not the total burden across all pathogen-k patients.

Disability weight caveat: Proxy_YLD_per_case.csv's DW_sepsis column is empirically a duration-inclusive "YLD per case" proxy (proxy_yld_per_case == yld_days_per_case / 365), not a pure disability weight. Because this function already multiplies by profile-specific LOSR_years separately, sourcing DW_L from yld_ref risks double-counting duration; a warning is emitted when DW_sepsis is not supplied directly. Prefer supplying a pure disability weight via the DW_sepsis scalar for this pathway.

Value

P_Lk_prime_tbl augmented with I_L, DW_L, profile_term (= sum_delta R'_k_delta * LOSR_k_delta), and YLD_associated.

Description

Primary, direct implementation of the profile-specific attributable-YLD equation:

Usage

daly_calc_yld_attributable(
  profiles_with_los,
  P_Lk_prime_tbl,
  yld_ref = NULL,
  DW_sepsis = NULL,
  pathogen_col = "pathogen",
  plk_col = "P_Lk_prime",
  incidence_col = "N_NF_L",
  facility_col = NULL,
  facility_name = NULL,
  facility_state_map = NULL,
  state_col = "state",
  state_name = NULL,
  probability_col = "probability",
  dominant_class_col = "dominant_class",
  losr_col = "LOSR_years",
  loss_col = "LOSS_years"
)

Arguments

Details

YLD_attrib_k = I_L * P'_Lk * DW_L * sum_{delta in D_k} [R'_k_delta * (LOSR_k,d*(delta) - LOSS_k,d*(delta))]

The excess (LOSR - LOSS) term makes this a genuine counterfactual: how much MORE disability burden exists because these infections were resistant, versus if they had been susceptible to the same (dominant) class. D_k excludes the all-susceptible profile from the sum entirely (rather than including a zero LOSR - LOSS term for it) – see daly_calc_yld_associated() for why.

Obtains I_L, P'_Lk, DW_L, R'_k_delta, d*(delta), LOSR, and LOSS the same way as daly_calc_yld_associated() – see that function's documentation for full component-by-component sourcing and the disability-weight caveat. Does not require a baseline YLD or a PAF_LOS computed first; unlike the pre-existing multiplicative pathway (daly_calc_paf_los() multiplied against a pooled baseline YLD), this is the direct equation. daly_calc_paf_los() remains available standalone for diagnostics / comparison against this direct calculation.

Value

P_Lk_prime_tbl augmented with I_L, DW_L, profile_term (= sum_delta R'_k_delta * (LOSR_k_delta - LOSS_k_delta)), and YLD_attributable.

Description

Computes years of life lost (YLL) associated with AMR directly from patient-level facility records, without requiring population-level P_LK (syndrome fractions) or R_kd (resistance profile scalars).

Usage

daly_calc_yll_associated(
  data,
  outcome_col,
  death_value = "Death",
  pathogen_col,
  patient_col,
  age_bin_col,
  sex_col,
  facility_col = NULL,
  syndrome_col = NULL,
  syndrome_name = NULL,
  date_culture_col = NULL,
  specimen_col = NULL,
  poly_weight_method = "monomicrobial_proportion",
  min_mono_per_facility = 30L,
  gap_days = 14,
  le_path = system.file("extdata", "life_expectancy_all.xlsx", package = "anumaan"),
  male_value = "Male",
  female_value = "Female",
  age_bin_map = c(`<1` = "0-1"),
  stratify_by = NULL
)

Arguments

Details

For every fatal patient with pathogen k (and optionally syndrome L) the individual YLL contribution is:

$\text{YLL}_{r,k} = \text{LE}(\text{age\_bin}_r, \text{sex}_r) \times w_{r,k}$

where $w_{r,k}$ is the polymicrobial death weight for patient r and pathogen k (= 1 for monomicrobial, 0-1 for polymicrobial episodes). Total YLL associated:

$\text{YLL}_{\text{associated}} = \sum_{r,k} \text{YLL}_{r,k}$

Polymicrobial weights are computed via prep_flag_polymicrobial() + prep_compute_poly_weights(). When facility_col is provided the weights are derived per facility (reflecting local organism distributions), with automatic fallback to globally-pooled proportions for facilities whose monomicrobial reference pool is smaller than min_mono_per_facility. If date_culture_col is NULL polymicrobial flagging is skipped and all weights default to 1.

Value

A named list:

total

Scalar: total YLL associated across all pathogens and facilities.

per_pathogen

Data frame: YLL summed per pathogen k, pooled across facilities. Columns: pathogen_col, n_patients, YLL_associated_k.

by_age_sex

Data frame: YLL by age_bin_col x sex.

by_pathogen_age_sex

Data frame: YLL by pathogen x age_bin_col x sex.

by_facility

Data frame (only when facility_col is supplied): per-facility YLL, one row per facility x pathogen.

by_syndrome

Data frame (only when syndrome_col is supplied and syndrome_name is NULL): YLL by syndrome.

by_syndrome_pathogen

Data frame (only when syndrome_col supplied): YLL by syndrome x pathogen.

stratified

Data frame (only when stratify_by is supplied): YLL aggregated by the requested columns.

patient_data

The death-cohort data frame used for computation, with polymicrobial_weight, life_expectancy, and yll_contribution columns attached.

Examples

## Not run: 
yll <- daly_calc_yll_associated(
  data             = cohort_df,
  outcome_col      = "final_outcome",
  death_value      = "Death",
  pathogen_col     = "organism_name",
  patient_col      = "PatientInformation_id",
  age_bin_col      = "Age_bin",
  sex_col          = "gender",
  facility_col     = "center_name",
  syndrome_col     = "infectious_syndrome",
  date_culture_col = "culture_date",
  specimen_col     = "sample_type",
  stratify_by      = c("location", "infectious_syndrome")
)
yll$total
yll$per_pathogen
yll$by_age_sex
yll$stratified

## End(Not run)

Description

Takes the per-patient YLL data produced by daly_calc_yll_associated() (which already contains life_expectancy, death_weight, and yll_contribution) and multiplies by the mortality PAF from daly_calc_paf_mortality() to produce AMR-attributable YLL.

Usage

daly_calc_yll_attributable(
  yll_patient_data,
  paf_mort,
  pathogen_col,
  patient_col,
  age_bin_col,
  sex_col = ".sex_norm",
  resistance_profile_col = NULL,
  profile_col = "profile",
  facility_col = NULL,
  syndrome_col = NULL,
  stratify_by = NULL
)

Arguments

Details

Two PAF modes:

PAF_k scalar mode (default)

When resistance_profile_col is NULL, the overall mortality PAF per pathogen k (PAF_k_mort) from paf_mort is used as a scalar multiplier.

Per-profile mode

When resistance_profile_col is supplied, each patient row is matched to its resistance profile delta and the profile-specific PAF_mortality is applied. Patients whose profile does not appear in paf_mort receive NA and a warning is issued.

$\text{YLL}^{\text{attr}}_{i,k} = \text{yll\_contribution}_{i,k} \times \text{PAF}_{k(,\delta)}$

Value

A list:

total

Scalar: total AMR-attributable YLL.

per_pathogen

One row per pathogen: pathogen_col, n_patients, YLL_associated_k, PAF_k_mort, YLL_attributable_k.

by_age_sex

YLL attributable stratified by age_bin x sex.

by_pathogen_age_sex

YLL attributable by pathogen x age_bin x sex.

by_facility

Per-facility breakdown (when facility_col supplied).

by_syndrome

Per-syndrome breakdown (when syndrome_col supplied).

by_syndrome_pathogen

Per-syndrome x pathogen (when syndrome_col supplied).

stratified

User-defined stratification (when stratify_by supplied).

patient_data

Patient-level data augmented with PAF_kd and YLL_attributable_contribution.

Examples

## Not run: 
yll_assoc <- daly_calc_yll_associated(data = cohort, ...)
mort_or <- fit_mortality_rr_logistic(data = rr_data, ...)
profiles_or <- assign_rr_to_profiles(profiles_out,
  rr_table = mort_or,
  rr_col = "OR_death"
)
paf_mort <- daly_calc_paf_mortality(profiles_or)

yll_attr <- daly_calc_yll_attributable(
  yll_patient_data = yll_assoc$patient_data,
  paf_mort         = paf_mort,
  pathogen_col     = "organism_name",
  patient_col      = "PatientInformation_id",
  age_bin_col      = "Age_bin",
  sex_col          = ".sex_norm",
  facility_col     = "center_name",
  syndrome_col     = "infectious_syndrome",
  stratify_by      = c("location", "infectious_syndrome")
)

## End(Not run)

Description

Computes LOS with infection-type-specific clock start:

• CAI: LOS = date_discharge - date_admission

• HAI: LOS = date_discharge - date_culture

• Unknown / Not known: LOS = date_discharge - date_culture (if culture date present), else date_discharge - date_admission. If neither reference date is available the patient is excluded.

Only discharged patients are retained. A checkpoint reports how many are missing a discharge date; those patients fall back to los_col when provided, otherwise they are excluded. Rows with LOS <= 0 or LOS > max_los are dropped. Returns one row per patient.

Usage

daly_compute_patient_los(
  data,
  patient_id_col = "PatientInformation_id",
  facility_col = "center_name",
  facility_name = NULL,
  organism_col = "organism_name",
  syndrome_col = "syndrome",
  syndrome_name = NULL,
  date_admission_col = "date_of_admission",
  date_discharge_col = "final_outcome_date",
  date_culture_col = "date_of_first_positive_culture",
  final_outcome_col = "final_outcome",
  final_outcome_value = "Discharged",
  infection_type_derived_col = "infection_type_derived",
  los_col = NULL,
  max_los = 200
)

Arguments

Value

Data frame: one row per patient with patient_id_col, facility_col, infection_type_derived_col, LOS_days.

Description

Counts the number of unique patients per infectious syndrome directly from facility-level data. This is the direct-count approach to incidence – no CFR, no CR_L adjustment, no pathogen weighting. Use this when you want raw facility-reported case counts rather than the formula-derived estimate from calculate_incidence_L().

Usage

daly_count_incident_cases(
  data,
  syndrome_col,
  syndrome_name,
  patient_col,
  facility_col = NULL,
  facility_name = NULL,
  pathogen_col = NULL,
  pathogen_name = NULL
)

Arguments

Details

Results are returned facility-wise when facility_col is supplied and facility_name is NULL. When facility_name is specified, only that facility is returned. When no facility information is provided, a single pooled count across all records is returned.

Value

Data frame with columns: syndrome_col, n_cases (unique patient count), and facility_col if supplied.

Description

Derives the healthcare-associated infection (HAI) vs community-acquired infection (CAI) classification for each patient record in the mortality cohort, using admission-to-culture date differences. Also performs a data-quality check for deceased patients missing an outcome date.

Usage

daly_derive_hai_cai_for_mortality(
  data,
  infection_type_col = "type_of_infection",
  date_admission_col = "date_of_admission",
  date_culture_col = "date_of_first_positive_culture",
  final_outcome_col = "final_outcome",
  final_outcome_date_col = "final_outcome_date",
  death_value = "Death",
  hai_threshold_hours = 48,
  patient_id_col = "PatientInformation_id"
)

Arguments

Value

The input data frame with a derived infection_type column (values "HAI" or "CAI").

Description

Drops profiles whose resistant classes have no RR estimate in rr_table and re-normalises the remaining probabilities. The all-susceptible reference profile is always kept.

Usage

daly_filter_profiles_to_rr_classes(
  profiles_with_rr,
  rr_table,
  pathogen_col = "pathogen",
  class_col = "antibiotic_class",
  probability_col = "probability",
  fallback_rr = 1
)

Arguments

Value

Named list (one entry per pathogen) with filtered and re-normalised profiles.

Description

Estimates class-specific relative length of stay (RR_LOS) for each pathogen using a Gamma generalized linear model with log link:

Usage

daly_fit_los_rr(
  data,
  patient_id_col = "PatientInformation_id",
  facility_col = "center_name",
  organism_col = "organism_name",
  syndrome_col = "syndrome",
  infection_type_col = "infection_type",
  antibiotic_class_col = "antibiotic_class",
  antibiotic_name_col = "antibiotic_name",
  antibiotic_value_col = "antibiotic_value",
  date_admission_col = "date_of_admission",
  date_discharge_col = "final_outcome_date",
  date_culture_col = "date_of_first_positive_culture",
  final_outcome_col = "final_outcome",
  final_outcome_value = "Discharged",
  age_col = "Age",
  sex_col = "Gender",
  unit_type_col = NULL,
  comorbidity_col = NULL,
  syndrome_name = NULL,
  organism_name = NULL,
  hai_threshold_hours = 48,
  max_los = 365,
  min_n = 20L,
  min_resistant = 5L,
  min_susceptible = 5L
)

Arguments

Details

E(LOS_i) = exp(beta0 + beta1*Resistance_ci + beta2*HAI_i + beta3*Age_i + beta4*Sex_i + beta5*ICU_i + beta6*Comorbidity_i + centre fixed effects)

RR_LOS for a resistance class is exp(beta1).

Value

Data frame with pathogen, antibiotic_class, hospital (NA when facility_col is NULL), RR_LOS, CI_lower, CI_upper, beta_log_rr, mean_los_resistant, mean_los_susceptible, n_patients, n_resistant, n_susceptible, model_family, link_function, and syndrome_scope.

Description

For each pathogen x antibiotic class, splits patients into resistant (R) and susceptible (S) groups, fits a parametric distribution (default: gamma) to each group's LOS, and returns: RR_LOS = mean_LOS(R) / mean_LOS(S)

Usage

daly_fit_los_rr_distribution(
  data,
  patient_id_col = "PatientInformation_id",
  facility_col = "center_name",
  organism_col = "organism_name",
  syndrome_col = "syndrome",
  infection_type_col = "infection_type",
  antibiotic_class_col = "antibiotic_class",
  antibiotic_name_col = "antibiotic_name",
  antibiotic_value_col = "antibiotic_value",
  date_admission_col = "date_of_admission",
  date_discharge_col = "final_outcome_date",
  date_culture_col = "date_of_first_positive_culture",
  final_outcome_col = "final_outcome",
  final_outcome_value = "Discharged",
  syndrome_name = NULL,
  organism_name = NULL,
  facility_name = NULL,
  hai_threshold_hours = 48,
  distributions = "gamma",
  max_los = 365,
  min_n = 10L
)

Arguments

Details

This is an alternative to daly_fit_los_rr() that avoids regression and models the full LOS distribution per group. The output flat data frame is fully compatible with assign_rr_to_profiles() and filter_profiles_to_rr_classes().

Value

A data frame with one row per pathogen x class:

pathogen

Pathogen name.

antibiotic_class

Antibiotic class name.

RR_LOS

Fitted mean ratio: mean_LOS(R) / mean_LOS(S).

mean_LOS_R

Fitted mean LOS for resistant patients.

mean_LOS_S

Fitted mean LOS for susceptible patients.

best_dist_R

Distribution fitted to the R group.

best_dist_S

Distribution fitted to the S group.

n_R

Number of resistant patients.

n_S

Number of susceptible patients.

syndrome_scope

Syndrome filter applied, or "all".

See Also

daly_fit_los_rr, daly_assign_rr_to_profiles

Description

Fits a pathogen-class specific mixed-effects logistic regression model:

Usage

daly_fit_mortality_rr(
  data,
  patient_id_col = "PatientInformation_id",
  facility_col = "center_name",
  organism_col = "organism_name",
  syndrome_col = "syndrome",
  infection_type_col = "type_of_infection",
  antibiotic_class_col = "antibiotic_class",
  antibiotic_name_col = "antibiotic_name",
  antibiotic_value_col = "antibiotic_value",
  unit_type_col = "unit_type",
  date_admission_col = "date_of_admission",
  date_culture_col = "date_of_first_positive_culture",
  final_outcome_col = "final_outcome",
  final_outcome_date_col = "final_outcome_date",
  age_col = "Age",
  sex_col = "Gender",
  comorbidity_col = NULL,
  death_value = "Death",
  syndrome_name = NULL,
  organism_name = NULL,
  hai_threshold_hours = 48,
  icu_values = c("ICU", "Intensive Care", "Critical Care", "PICU", "NICU"),
  phi_threshold = 0.7,
  min_n = 20L,
  min_deaths = 10L
)

Arguments