Package {hpsreg}

hpsreg: Regression for Hierarchically Ordered Clinical Outcomes (Bradley-Terry-Davidson Model)

Description

A three-stage workflow:

1. Design the hierarchy with hierarchy() and the ⁠lvl_*()⁠ constructors (lvl_binary(), lvl_ordinal(), lvl_continuous(), lvl_tte(), lvl_custom()): binary, ordinal, continuous and time-to-event outcomes, tie margins, eligibility guards, groups.

2. Score pairs with score_pairs() (compiled kernel, early exit at the first level that separates a pair), summarised by resolution().

3. Estimate: covariate-adjusted Davidson model with sandwich standard errors, btd(); unadjusted win_stats(); checks with btd_calibration().

simulate_hps() and hps_hierarchy() provide a worked stem cell transplant example (worked example: system.file("examples", "hps-example.Rmd")).

Author(s)

Maintainer: Ravi Varadhan ravi.varadhan@jhu.edu

Covariate-adjusted analysis of hierarchical pairwise outcomes

Description

Fits the Bradley-Terry-Davidson model to the win / loss / tie outcomes produced by score_pairs(). For a pair of subjects with covariate difference x_i - x_j and d = (x_i - x_j)^\top \beta / 2,

P(i\ \mathrm{wins}) = e^{d}/q,\quad P(j\ \mathrm{wins}) = e^{-d}/q, \quad P(\mathrm{tie}) = \delta/q,

with q = e^{d} + e^{-d} + \delta. A coefficient is the log of the multiplicative effect of a one-unit increase of the covariate on a subject's "worth"; positive values mean a better position in the hierarchy. If no pair is tied the model reduces to Bradley-Terry (\delta = 0).

Usage

btd(
  formula,
  data,
  pairs = NULL,
  hierarchy = NULL,
  scale = TRUE,
  engine = c("cpp", "r"),
  storage = c("auto", "stored", "compact", "stream"),
  control = btd_control()
)

## S3 method for class 'btd'
summary(object, exponentiate = FALSE, ...)

Arguments

Details

Estimation is by Levenberg-Marquardt on standardised covariates (C++ engine, or a pure R reference implementation with engine = "r"). Standard errors are sandwich (robust) estimates based on the projection of the pairwise estimating equation onto subjects, which accounts for the dependence between pairs that share a subject (Hoeffding / U-statistic variance). They assume the complete set of N(N-1)/2 pairs; for a sampled or user-supplied subset the normalisation uses the average number of pairs per subject and should be treated as approximate. The model-based inverse information, which ignores that dependence, is returned as vcov_model for reference only.

Value

An object of class btd with coefficients, log_delta, delta, vcov (sandwich, on the original covariate scale, including log_delta), vcov_model, loglik, convergence information, the latent scores eta and the pairs used.

The object supports the usual generics: summary() (coefficient table with z-tests), coef(), vcov() (model = TRUE for the model-based information), confint() (Wald intervals), logLik(), nobs() (number of pairs; the subject count is fit$n_subjects), print(), and predict() (latent score for the training rows, or for a newdata frame).

Pair-level predicted probabilities are provided by btd_pair_probs(); calibration tables by btd_calibration().

Two workflows

btd() needs the scored pairs of a hierarchy. You can supply them or let btd() compute them:

1. Pre-computed pairs (pairs = pw). Score once with score_pairs(), inspect with resolution() or win_stats(), then fit several covariate models on the same pairs without rescoring.

2. Inline scoring (hierarchy = h). btd() calls score_pairs() for you; the pairs are kept in the fitted object (fit$pairs) for later use.

In both cases data must be the same data frame (same rows, same order) that is scored.

References

Varadhan R. "Hierarchical Pairwise Scoring: A Regression Framework for Composite Outcomes, with Application to Hematopoietic Stem Cell Transplantation." Manuscript in preparation, 2026.

Examples

## From scratch: define a hierarchy from data-frame columns, score, fit.
dat <- simulate_hps(N = 150, seed = 11)

h <- hierarchy(
  lvl_tte("death_day", better = "later", name = "survival"),
  lvl_tte("relapse_day", better = "later", name = "relapse"),
  lvl_ordinal("gvhd_grade", better = "lower", name = "GVHD grade"),
  lvl_binary("gvhd_resolved", favorable = 1, when = "gvhd_grade > 0",
             name = "GVHD resolution")
)

## Workflow 1: pre-computed pairs (reusable across models)
pw <- score_pairs(h, dat)
resolution(pw)
fit1 <- btd(~ donor_type + age + dri + hct_ci, data = dat, pairs = pw)
summary(fit1)
fit2 <- btd(~ donor_type + age, data = dat, pairs = pw)   # no rescoring

## Workflow 2: let btd() score the pairs
fit3 <- btd(~ donor_type + age, data = dat, hierarchy = h)
all.equal(coef(fit2), coef(fit3))

## The packaged transplant hierarchy is a shortcut for the one above
fit4 <- btd(~ donor_type + age, data = dat, hierarchy = hps_hierarchy())

## How many threads can this build use? 0 means no OpenMP, in which
## case the engine runs serially regardless of the `threads` argument.
openmp_threads()


## Multiple threads: results are bitwise identical to a single-threaded
## fit, so this only changes the runtime, not the estimates.
n <- max(1L, min(4L, openmp_threads()))
fit_mt <- btd(~ donor_type + age, data = dat, hierarchy = h,
              control = btd_control(threads = n))
all.equal(coef(fit_mt), coef(fit3))

## The same setting can be supplied through an option, which is
## convenient when many calls are made in a script.
op <- options(hpsreg.threads = n)
fit_mt2 <- btd(~ donor_type + age, data = dat, hierarchy = h)
options(op)
all.equal(coef(fit_mt), coef(fit_mt2))

Calibration of the pairwise model

Description

Compares observed outcomes of the scored pairs with the model's predicted probabilities.

Usage

btd_calibration(object, by = c("level", "bin"), n_bins = 10L)

Arguments

Details

• by = "level": among pairs decided at each hierarchy level, the observed share won by the first-listed subject versus the model's P(i)/(P(i)+P(j)). Systematic differences between levels suggest that a single latent scale does not describe all levels.

• by = "bin": pairs are grouped by quantiles of the predicted latent difference; observed and predicted win / tie / loss fractions are compared. For fits made with streaming (no stored pairs) the bins are approximate quantile bins built from a fine histogram of the latent difference.

The returned object is plotted by plot.btd_calibration, which draws a dot plot for by = "level" and a three-panel observed-vs-predicted line plot for by = "bin".

Value

A btd_calibration object (a data frame with additional attributes) with columns

• by = "level": level, n, obs_win_i, pred_win_i.

• by = "bin": bin, n, obs_win_i, pred_win_i, obs_tie, pred_tie, obs_win_j, pred_win_j.

See Also

plot.btd_calibration for the graphical counterpart; btd() to fit the model; score_pairs() to score the pairs.

Examples

dat <- simulate_hps(N = 120, seed = 5)
h   <- hps_hierarchy()
fit <- btd(~ donor_type + dri, dat, hierarchy = h)

## Calibration by hierarchy level: for each level, the observed and
## predicted share of pairs won by the first subject.
cal_level <- btd_calibration(fit, by = "level")
cal_level
plot(cal_level)

## Calibration by quantile bins of the latent difference d:
## observed vs predicted win / tie / loss fractions.
cal_bin <- btd_calibration(fit, by = "bin", n_bins = 8)
cal_bin
plot(cal_bin)

Control parameters for the Bradley-Terry-Davidson fit

Description

Control parameters for the Bradley-Terry-Davidson fit

Usage

btd_control(
  max_iter = 200L,
  tol = 1e-06,
  lambda0 = 1e-06,
  verbose = FALSE,
  threads = NULL
)

Arguments

Value

A list.

Outcome levels for a hierarchy

Description

Constructors for the building blocks of a hierarchy(). Each level compares two subjects on one outcome and returns "first wins", "second wins" or "neutral" (the pair is passed to the next level).

Usage

lvl_binary(var, favorable = 1, when = NULL, name = var, na = "tie")

lvl_ordinal(
  var,
  better,
  order = NULL,
  threshold = 0,
  when = NULL,
  name = var,
  na = "tie"
)

lvl_continuous(var, better, threshold = 0, when = NULL, name = var, na = "tie")

lvl_tte(
  time,
  status = NULL,
  better,
  horizon = Inf,
  threshold = 0,
  when = NULL,
  name = time
)

lvl_custom(fun, name = "custom", when = NULL)

Arguments

Value

An object of class hb_level.

Eligibility guards

when is an expression evaluated in the data (e.g. when = gvhd_grade > 0). Inside package code, where an unquoted column name would trigger an "undefined global variable" note in ⁠R CMD check⁠, pass it as a string (when = "gvhd_grade > 0"); both forms are equivalent. If either subject of a pair is ineligible (guard FALSE or NA) the level is neutral for that pair. This expresses conditional comparisons such as "compare resolution of GVHD only for patients who had GVHD".

Missing values

For lvl_binary(), lvl_ordinal() and lvl_continuous() the argument na controls pairs in which exactly one subject has a missing value: "tie" (default) makes the level neutral, "best" / "worst" treat a missing value as the best / worst possible outcome. Pairs in which both are missing are always neutral.

Time-to-event levels

With status = NULL the time variable is read as "time of the event, NA if the event did not occur" and follow-up is assumed complete (no censoring); non-events are placed at horizon (or infinity). With a status column (1 = event, 0 = censored) pairs are compared with the Gehan rule: a pair is decided only if the ordering is certain under right censoring, otherwise it is uninformative and passed to the next level (and flagged, see score_pairs()). Follow-up beyond horizon is truncated and treated as complete at horizon. This rule is simple and biased toward the null under heavy censoring; inverse-probability-of-censoring scoring is not yet implemented.

Examples

lvl_tte("death_day", better = "later", name = "survival")
lvl_ordinal("gvhd_grade", better = "lower")
lvl_binary("resolved", favorable = 1, when = gvhd_grade > 0)
lvl_continuous("pain_score", better = "lower", threshold = 1)

Define a hierarchy of outcomes

Description

Levels are listed from highest to lowest priority. When two subjects are compared, the first level on which they differ decides who wins; if no level separates them the pair is a tie.

Usage

hierarchy(...)

Arguments

Value

An object of class hb_hierarchy.

Examples

h <- hierarchy(
  lvl_tte("death_day", better = "later", name = "survival"),
  lvl_group("GVHD",
    lvl_ordinal("gvhd_grade", better = "lower"),
    lvl_binary("gvhd_resolved", favorable = 1, when = gvhd_grade > 0)
  )
)
h

Data-generating coefficients for the transplant example

Description

Starting values for simulate_hps(): the negative log hazard ratios of published m-GRFS effects (uncalibrated), with fixed age and comorbidity effects. Use these as a convenient default; effects should be re-calibrated to a target if exact hazard ratios matter.

Usage

hps_beta()

Value

A named numeric vector.

Examples

hps_beta()

Composite m-GRFS Cox model (comparator for the example)

Description

Fits a Cox model for modified GVHD- and relapse-free survival: the first of death, relapse, grade >= 2 GVHD or grade >= 2 infection within the horizon. This is the conventional composite analysis that the hierarchical analysis is compared with in the example. Requires the survival package.

Usage

hps_cox_grfs(dat, covars, horizon = 180)

Arguments

Value

A list with fit, coef and se.

Examples

dat <- simulate_hps(N = 200, seed = 4)
hps_cox_grfs(dat, c("donor_type", "age"))$coef

Hierarchy for the stem cell transplant example

Description

Four levels, evaluated in this order: survival (later death better, non-death best), relapse (later relapse better, non-relapse best), GVHD (lower grade; then resolution; then later onset) and infection (lower grade; then resolution). Resolution and onset are only compared between patients who had the event.

Usage

hps_hierarchy()

Value

A hierarchy().

Examples

hps_hierarchy()

Group several levels into one reporting level

Description

A group is evaluated as consecutive levels (in the order given) but is reported as a single level, e.g. GVHD grade, then resolution, then onset time are three comparisons forming one "GVHD" level of the hierarchy.

Usage

lvl_group(name, ..., when = NULL)

Arguments

Value

An object of class hb_group.

OpenMP availability

Description

Number of threads OpenMP can use in this build of the package, or 0 if the package was compiled without OpenMP. The streaming engine uses one thread by default; set options(hpsreg.threads = n) (or btd_control(threads = n)) to use more. Results are bitwise identical for any number of threads.

Usage

openmp_threads()

Value

An integer.

Examples

openmp_threads()

Plot a calibration table

Description

Two graphical summaries are produced depending on the by argument of btd_calibration():

Usage

## S3 method for class 'btd_calibration'
plot(x, ...)

Arguments

Details

• by = "level": a dot plot with one row per hierarchy level, showing the observed and predicted win share for pairs resolved at that level. Lines of agreement between the two are what a well-calibrated hierarchy should look like.

• by = "bin": three line plots (win, tie, loss) of observed against predicted probabilities across quantile bins of the latent difference d.

The bin plot shows observed and predicted probabilities only; no binomial error bars are drawn. Pairs share subjects, so treating them as independent Bernoulli draws would understate the sampling variability. For a variability estimate use the subject-clustered standard errors from btd() rather than a pair-level binomial interval.

Value

x, invisibly.

See Also

btd_calibration()

Examples

dat <- simulate_hps(N = 120, seed = 5)
fit <- btd(~ donor_type + dri, dat, hierarchy = hps_hierarchy())

plot(btd_calibration(fit, by = "level"))
plot(btd_calibration(fit, by = "bin", n_bins = 8))

Predicted latent score and pair probabilities

Description

predict() returns the latent score (x - \bar x)^\top \beta (centred at the training covariate means; only differences matter). btd_pair_probs() returns the model probabilities that subject A wins, that the pair is tied, and that subject B wins, for rows of a compared with the same rows of b.

Usage

## S3 method for class 'btd'
predict(object, newdata, ...)

btd_pair_probs(object, a, b)

Arguments

Value

A numeric vector, or a matrix with columns win_a, tie, win_b.

Examples

dat <- simulate_hps(N = 100, seed = 3)
fit <- btd(~ donor_type + age, dat, pairs = score_pairs(hps_hierarchy(), dat))
btd_pair_probs(fit, dat[1:3, ], dat[4:6, ])

Where are pairs decided?

Description

Number and share of pairs resolved at each level of the hierarchy, and the share left tied. Pass either scored pairs, or a hierarchy together with the data: in that case the pairs are generated and counted on the fly and are never stored, so this works for very large N.

Usage

resolution(pairs, by = c("level", "leaf"), data = NULL)

Arguments

Value

A data frame with columns level, n (pairs decided there) and prop (share of all pairs). The last row, (tied), counts pairs that no level separates.

Examples

dat <- simulate_hps(N = 80, seed = 2)
resolution(score_pairs(hps_hierarchy(), dat))
resolution(hps_hierarchy(), data = dat)      # streaming, nothing stored

Sample a random set of pairs

Description

Draws n distinct pairs uniformly from the N(N-1)/2 possible pairs without enumerating them. Useful when N is too large for all pairs.

Usage

sample_pairs(N, n)

Arguments

Value

A two-column integer matrix ⁠(i, j)⁠ with i < j.

Examples

set.seed(1); head(sample_pairs(1000, 5))

Score all pairs of subjects under a hierarchy

Description

Applies a hierarchy() to every pair of subjects (or to a supplied set of pairs). Built-in levels are evaluated by a compiled kernel that stops at the first level that separates a pair; lvl_custom() levels are evaluated in R on the pairs still unresolved.

Usage

score_pairs(hierarchy, data, pairs = NULL, max_pairs = 5e+07)

Arguments

Value

An object of class hb_pairs with components id_i, id_j (row indices), W (1 if i wins, -1 if j wins, 0 if tied on every level), level (index of the resolving top-level level, K + 1 if unresolved), leaf (finer index when the hierarchy contains groups, else NULL), uninf (1 if an earlier level was uninformative because of censoring), N, complete and the hierarchy.

Examples

dat <- simulate_hps(N = 100, seed = 1)
pw <- score_pairs(hps_hierarchy(), dat)
pw
resolution(pw)

Simulate a transplant cohort

Description

Simulates covariates and 180-day outcomes (death, relapse, graft-versus-host disease, infection) for the worked example. Outcomes share a latent frailty-like score driven by the covariates in beta.

Usage

simulate_hps(beta = hps_beta(), N = 613, seed = NULL)

Arguments

Value

A data frame with covariates (donor_type, age, sex, dri, hct_ci, cell_source, conditioning, year, sex_mismatch, diagnosis), the latent score, and outcome columns (death_day, relapse_day, gvhd_grade, gvhd_resolved, gvhd_onset_day, infection_grade, infection_resolved, infection_onset_day). NA in a ⁠*_day⁠ column means the event did not occur by day 180.

Examples

dat <- simulate_hps(N = 200, seed = 1)
head(dat[, c("donor_type", "age", "death_day", "gvhd_grade")])

Win statistics for two groups

Description

Summarises the cross-group pairs of a scored hierarchy: the proportion of pairs won by each group, and the derived net benefit, win ratio and win odds of treatment versus the reference group, with U-statistic (Hoeffding) variances that account for each subject appearing in many pairs.

Usage

win_stats(pairs, group, treatment = NULL, conf.level = 0.95, data = NULL)

Arguments

Details

With p_w and p_l the proportions of cross-group pairs won by the treatment and reference groups and p_t = 1 - p_w - p_l: net benefit p_w - p_l; win ratio p_w/p_l; win odds (p_w + p_t/2)/(p_l + p_t/2). Confidence intervals for the two ratios are computed on the log scale (delta method).

Requires the complete set of pairs. Pass a hierarchy together with data instead of scored pairs to compute everything on the fly without storing pairs (memory O(N)); this is the way to analyse very large N.

Value

An object of class hb_win with the summary table, the by_level contribution of each hierarchy level to wins and losses, and the counts.

Examples

dat <- simulate_hps(N = 150, seed = 7)
pw <- score_pairs(hps_hierarchy(), dat)
win_stats(pw, dat$year)
win_stats(hps_hierarchy(), dat$year, data = dat)   # streaming, same result