Getting started with hpsreg

What this package does

hpsreg fits a covariate-adjusted Bradley-Terry-Davidson model to hierarchically ordered clinical outcomes. Instead of picking one endpoint, you list outcomes in order of clinical priority. Each pair of subjects is compared level by level; the first level that separates them decides the pair. Only pairs that tie at every level count as ties.

The workflow has three stages:

  1. Design a hierarchy of outcomes with hierarchy() and the lvl_*() constructors.
  2. Score every pair of subjects with score_pairs().
  3. Estimate covariate effects with btd().
library(hpsreg)

Example data

simulate_hps() returns a synthetic stem-cell transplant cohort with the same column layout as a typical registry extract. The workflow below is identical for real data with columns of the same type.

dat <- simulate_hps(N = 200, seed = 11)
str(dat, list.len = 6)
#> 'data.frame':    200 obs. of  24 variables:
#>  $ id                 : int  1 2 3 4 5 6 7 8 9 10 ...
#>  $ donor_type         : Factor w/ 3 levels "MRD","MUD","HID": 2 2 1 2 2 3 2 2 3 2 ...
#>  $ age                : num  46 57 62 51 22 62 60 38 72 70 ...
#>  $ sex                : Factor w/ 2 levels "F","M": 2 1 2 2 2 2 1 1 1 1 ...
#>  $ dri                : Factor w/ 4 levels "Low","Intermediate",..: 1 2 2 2 1 4 2 2 1 1 ...
#>  $ hct_ci             : Factor w/ 2 levels "0-2",">=3": 2 2 1 2 2 2 2 2 1 2 ...
#>   [list output truncated]

Design a hierarchy

Levels are listed from highest to lowest priority. lvl_tte() compares time-to-event columns, lvl_ordinal() and lvl_binary() handle categorical outcomes, and lvl_continuous() handles continuous scores. The when = guard limits a level to subjects that meet a condition (used here so the GVHD resolution comparison only applies to subjects who had GVHD).

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")
)
h
#> <hierarchy> 4 levels (first level that separates a pair decides it)
#>  1. survival: tte [death_day], later better (complete follow-up)
#>  2. relapse: tte [relapse_day], later better (complete follow-up)
#>  3. GVHD grade: ordinal [gvhd_grade], lower better
#>  4. GVHD resolution: binary [gvhd_resolved], favourable = 1 (when gvhd_grade > 0)

hps_hierarchy() returns the fuller transplant hierarchy used in the transplant examples, adding GVHD onset and infection levels.

Censoring in time-to-event levels

For lvl_tte() levels, a subject whose event has not yet been observed (right-censored) is compared with the other subject using the Gehan rule: if the event time is beyond the censoring time, the pair is uninformative at that level and passes to the next one. The package tracks these pairs, and summary() reports what fraction of tied pairs are tied because of censoring rather than genuine equivalence:

Ties: 19.7% of pairs, of which 82% pure and 18% censoring

Heavy censoring weakens the survival level and pushes more discrimination onto the lower levels. The Gehan rule is conservative under heavy censoring; IPCW (Peron) scoring is not implemented.

Score the pairs

score_pairs() enumerates every pair of subjects and scores it against the hierarchy. The compiled kernel exits at the first separating level, and the resolution() summary shows where the discrimination comes from.

pw <- score_pairs(h, dat)
pw
#> <hb_pairs> 19,900 pairs among 200 subjects
#>   wins/losses: 16292 resolved, 3608 tied on every level (18.1%)
resolution(pw)
#>             level    n       prop
#> 1        survival 3790 0.19045226
#> 2         relapse 3547 0.17824121
#> 3      GVHD grade 8201 0.41211055
#> 4 GVHD resolution  754 0.03788945
#> 5          (tied) 3608 0.18130653

In this cohort survival is the decisive level for most pairs, and only a small fraction reach the GVHD-resolution comparison. Those remaining pairs are the ones that determine the tie parameter.

The scored object pw can be reused across as many covariate models as you like — scoring is the only step that depends on the hierarchy. This is why btd() accepts pairs separately from hierarchy: score once, then fit ~ donor_type, ~ donor_type + age, ~ donor_type + age + dri, and compare without rescoring.

Unadjusted win statistics

win_stats() is the descriptive companion to the model: the net benefit, win ratio, and win odds for one arm against another, with a per-level breakdown of where the wins and losses come from. It takes either scored pairs or a hierarchy, and needs a grouping variable.

win_stats(pw, dat$year)
#> Win statistics: '2011-2015' (n = 110) vs '2005-2010' (n = 90)
#> Cross-group pairs: win 0.398 | loss 0.417 | tie 0.185
#> 
#>    statistic estimate   lower  upper       z      p
#>  net_benefit -0.01889 -0.1764 0.1386 -0.2351 0.8141
#>    win_ratio  0.95470  0.6485 1.4050 -0.2349 0.8143
#>     win_odds  0.96290  0.7027 1.3200 -0.2350 0.8142
#> 
#> 95% CI; ratios use log-scale intervals. Contribution by level:
#>            level treatment_wins treatment_losses   tied
#>         survival         0.1280           0.0684 0.0000
#>          relapse         0.0367           0.1315 0.0000
#>       GVHD grade         0.2125           0.2001 0.0000
#>  GVHD resolution         0.0210           0.0171 0.0000
#>           (tied)         0.0000           0.0000 0.1847

Passing the hierarchy directly is convenient when you have not scored yet:

win_stats(h, dat$year, data = dat)
#> Win statistics: '2011-2015' (n = 110) vs '2005-2010' (n = 90)
#> Cross-group pairs: win 0.398 | loss 0.417 | tie 0.185
#> 
#>    statistic estimate   lower  upper       z      p
#>  net_benefit -0.01889 -0.1764 0.1386 -0.2351 0.8141
#>    win_ratio  0.95470  0.6485 1.4050 -0.2349 0.8143
#>     win_odds  0.96290  0.7027 1.3200 -0.2350 0.8142
#> 
#> 95% CI; ratios use log-scale intervals. Contribution by level:
#>            level treatment_wins treatment_losses   tied
#>         survival         0.1280           0.0684 0.0000
#>          relapse         0.0367           0.1315 0.0000
#>       GVHD grade         0.2125           0.2001 0.0000
#>  GVHD resolution         0.0210           0.0171 0.0000
#>           (tied)         0.0000           0.0000 0.1847

The overall advantage is concentrated at the survival and relapse levels, not in the GVHD levels.

Fit the model

btd() fits the covariate-adjusted model. Each coefficient is the log of the multiplicative effect of a one-unit increase in the covariate on the subject’s worth; positive values mean a better position in the hierarchy. The intercept is dropped because it is not identifiable from pairwise differences.

fit <- btd(~ donor_type + age + dri, data = dat, pairs = pw)
summary(fit)
#> Hierarchical Bradley-Terry-Davidson (BTD) model
#> 
#> Call:
#> btd(formula = ~donor_type + age + dri, data = dat, pairs = pw)
#> 
#> 200 subjects, 19,900 pairs (18.1% tied); engine: cpp (pairs stored), converged in 4 iterations
#> 
#>                  Estimate Std. Error z value Pr(>|z|)    
#> donor_typeMUD   -0.535047   0.224303  -2.385 0.017062 *  
#> donor_typeHID   -0.134535   0.284414  -0.473 0.636196    
#> age              0.004181   0.007264   0.576 0.564909    
#> driIntermediate -0.433764   0.311263  -1.394 0.163450    
#> driHigh         -0.644495   0.355564  -1.813 0.069894 .  
#> driVery high    -1.596841   0.457029  -3.494 0.000476 ***
#> log_delta       -0.767523   0.127828  -6.004 1.92e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Standard errors: sandwich (U-statistic projection over subjects).
#> 
#> Tie parameter delta = 0.4642; log-likelihood = -19978.56
#>   Ties: 18.1% of pairs, of which 100% pure and 0% censoring

You can also let btd() score the pairs for you by passing hierarchy = h instead of pairs = pw. For small cohorts the two give identical estimates.

fit2 <- btd(~ donor_type + age + dri, data = dat, hierarchy = h)
all.equal(coef(fit), coef(fit2))
#> [1] TRUE

data must be the same data frame (same rows, same order) that was scored. btd() checks the row count against pairs$N and errors if they differ, but it cannot detect a reordered frame — the row indices in pw would point at the wrong subjects. If you subset or reorder, rescore.

Interpreting delta

Every summary(fit) prints a delta value and its log. These describe the model’s tie propensity — a weight on the probability that a pair ties, independent of the covariate effects. It is not a probability: delta = 0.49 means the tie term contributes roughly half as much as an even-matched win, and delta = 0 reduces the model to ordinary Bradley-Terry (used automatically when there are no tied pairs). Use btd_calibration() to check whether a single delta fits all levels.

Standard errors

The default standard errors are sandwich (robust) estimates that account for the dependence between pairs sharing a subject. The model-based values, which assume independent pairs and are typically too small, are available for reference:

cbind(sandwich = sqrt(diag(vcov(fit))),
      model    = sqrt(diag(vcov(fit, model = TRUE))))
#>                    sandwich        model
#> donor_typeMUD   0.224302540 0.0264546845
#> donor_typeHID   0.284413658 0.0333297133
#> age             0.007263554 0.0008346674
#> driIntermediate 0.311263413 0.0368945227
#> driHigh         0.355563947 0.0399803355
#> driVery high    0.457029363 0.0544482623

confint(fit) returns Wald intervals for the coefficients (but not for log_delta). nobs(fit) is the number of pairs — N(N-1)/2 for a complete set — not the number of subjects; the subject count is fit$n_subjects.

Prediction and calibration

predict() returns the latent score (log-worth) for each subject, and btd_pair_probs() returns response-scale probabilities for pairs.

predict(fit, dat[1:3, ])
#>          1          2          3 
#>  0.2477370 -0.1400401  0.4159097
btd_pair_probs(fit, dat[1:3, ], dat[4:6, ])
#>       win_a       tie     win_b
#> 1 0.4903544 0.1851515 0.3244941
#> 2 0.3485638 0.1867963 0.4646399
#> 3 0.6597236 0.1600509 0.1802256

btd_calibration() compares observed and predicted win shares, by hierarchy level or by bin of the latent difference. In a well-calibrated hierarchy the two sets of points lie on top of each other. Observed consistently above predicted at the lower levels would suggest the hierarchy discriminates those pairs more sharply than a single latent scale allows.

btd_calibration(fit, by = "level")
#>             level    n obs_win_i pred_win_i
#> 1        survival 3790 0.5287599  0.5225149
#> 2         relapse 3547 0.5181844  0.5040202
#> 3      GVHD grade 8201 0.4334837  0.5002090
#> 4 GVHD resolution  754 0.4801061  0.5025660
plot(btd_calibration(fit, by = "level"))

Large cohorts

Every pair is scored, so N = 10000 means about 50 million pairs. btd() offers three storage modes, chosen automatically based on N:

All three give the same estimates up to numerical precision. The mode actually used is reported in the fit summary.

fit <- btd(~ donor_type + age, data = dat, hierarchy = h, storage = "compact")

score_pairs() refuses to hold more than max_pairs (default 5e7, about N = 10,000) pairs at once. For larger N, either raise max_pairs, sample with sample_pairs(), or pass the hierarchy directly to btd() — the compact and streaming engines score pairs on the fly and never materialise the full set.

Threads

The engine is parallelised across pairs with OpenMP when the build supports it. Results are bitwise identical regardless of thread count — threads only change runtime. openmp_threads() reports how many threads the build can use; 0 means OpenMP is not available (stock macOS toolchains, for example), in which case everything runs serially.

openmp_threads()
#> [1] 16

Set the thread count per call or globally:

fit_mt <- btd(~ donor_type + age, data = dat, hierarchy = h,
              control = btd_control(threads = 4))

options(hpsreg.threads = 4)          # applies to every subsequent call
fit_mt2 <- btd(~ donor_type + age, data = dat, hierarchy = h)
all.equal(coef(fit_mt), coef(fit_mt2))

Citation

If you use hpsreg in published work, please cite:

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

The same reference is available from R with citation("hpsreg").

Where to go next