higher_criticism_test: Higher Criticism test for non-proportional hazard deviations

In survHC: Higher Criticism Tests for Non-Proportional Hazard Deviations

Higher Criticism test for non-proportional hazard deviations

Description

Compares two right-censored survival samples using Higher Criticism applied to exact per-interval hypergeometric p-values. The test is powerful against sparse departures, where the hazard difference is concentrated in a small and a priori unknown subset of the follow-up, a regime in which the log-rank test has little power.

Usage

higher_criticism_test(durations_A, ...)

## Default S3 method:
higher_criticism_test(
  durations_A,
  durations_B,
  event_observed_A = NULL,
  event_observed_B = NULL,
  alternative = c("both", "greater", "less"),
  gamma = 0.2,
  hc_version = c("dj2008", "dj2004", "beta", "star"),
  t_0 = -1,
  n_intervals_to_pool = NULL,
  n_permutations = 0,
  seed = NULL,
  ...
)

## S3 method for class 'formula'
higher_criticism_test(formula, data, subset, na.action, ...)

Arguments

durations_A, durations_B	Numeric vectors of event or censoring times. durations_A may instead be a formula, which selects the formula method.
event_observed_A, event_observed_B	Optional 0/1 vectors, 1 indicating an observed event and 0 a censored observation. Default: all events.
alternative	One of "both" (default), "greater" (excess hazard in group B) or "less" (excess hazard in group A).
gamma	Fraction of the smallest ordered p-values scanned, in (0, 1]. Default 0.2.
hc_version	Which Higher Criticism variant to use, matching the named methods of the Python package multitest (>= 0.2.2) and the hc_version argument of lifelines-hc. "dj2008" (default) uses u_i = i/N with denominator \sqrt{u_i(1-u_i)/N}; "dj2004" uses u_i = i/N with denominator \sqrt{p_{(i)}(1-p_{(i)})/N}, which is equation (2) of Kipnis, Galili and Yakhini (2026); "beta" uses u_i = i/(N+1) with denominator \sqrt{u_i(1-u_i)/(N+2)}; "star" uses the "beta" denominator but starts the scan at the first ordered p-value exceeding (1-\epsilon)/N. Ignored by berk_jones_test, fisher_combination_test and min_p_test.
t_0	Restrict the analysis to events at or before t_0. A negative value (default) applies no restriction.
n_intervals_to_pool	Number of equal-width intervals to pool event times into. Strongly recommended for continuous survival data, where each distinct time carries at most one event and the per-time hypergeometric p-values are too coarse to be informative. Values between 50 and 200 usually work well. NULL (default) applies no pooling.
n_permutations	Number of label permutations used to calibrate the p-value. 0 (default) skips calibration and returns NA.
seed	Optional integer seed for the permutation test. When supplied, the caller's random seed is restored on exit.
formula	A two-sided formula Surv(time, status) ~ group; see Details.
data	An optional data frame, list or environment holding the variables named in formula.
subset	An optional vector specifying a subset of observations to use.
na.action	A function describing what should happen when the data contain NAs. Defaults to getOption("na.action"), normally na.omit.
...	Further arguments. The formula method passes these on to the default method, which accepts none and signals an error if any are supplied, so that a misspelled argument name cannot be silently ignored.

Details

Two interfaces are available and return the same numbers. The default method takes the two samples as vectors, mirroring the Python package lifelines-hc. The formula method takes the usual survival formula Surv(time, status) ~ group, which needs the survival package for Surv itself; without it, cbind(time, status) ~ group and time ~ group also work, the latter treating every observation as an event.

The grouping variable must take exactly two values. Its first level is group A and its second is group B, so alternative = "greater" tests for excess hazard in the second level. Levels are taken in factor order for a factor, and otherwise in the order sort(unique(.)) gives: numerically for numbers, alphabetically for strings.

The two interfaces differ only in their treatment of missing values. The default method rejects NA durations; the formula method applies na.action, which by default drops incomplete rows.

Value

An object of class "htest", with the Higher Criticism score in statistic and the permutation p-value in p.value. The permutation statistics are returned in permutation_statistics. The formula method additionally records the two group levels, in the order (A, B), in group_levels.

References

Kipnis, A., Galili, B. and Yakhini, Z. (2026). Higher criticism for rare and weak non-proportional hazard deviations in survival analysis. Biometrika 113(1), asaf075. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1093/biomet/asaf075")}

Donoho, D. and Jin, J. (2004). Higher criticism for detecting sparse heterogeneous mixtures. The Annals of Statistics 32(3), 962–994. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1214/009053604000000265")}

Kipnis, A., Galili, B. and Yakhini, Z. (2026). lifelines-hc: Higher Criticism testing for sparse non-proportional hazard departures in Python. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.48550/arXiv.2610.07268")}

See Also

suspected_deviations to recover the intervals driving a detection; berk_jones_test, fisher_combination_test and min_p_test for alternative aggregations of the same per-interval p-values.

Examples

set.seed(1)
T_A <- rexp(200, rate = 1 / 10)
T_B <- rexp(200, rate = 1 / 10)
higher_criticism_test(T_A, T_B, n_intervals_to_pool = 50)

# with permutation calibration
higher_criticism_test(T_A, T_B, n_intervals_to_pool = 50,
                      n_permutations = 200, seed = 42)

# the same data through the formula interface
trial <- data.frame(time = c(T_A, T_B), status = 1,
                    arm = rep(c("A", "B"), c(200, 200)))
higher_criticism_test(cbind(time, status) ~ arm, data = trial,
                      n_intervals_to_pool = 50)

# Surv() needs the survival package, which is only suggested
if (requireNamespace("survival", quietly = TRUE)) {
  higher_criticism_test(survival::Surv(time, status) ~ arm, data = trial,
                        n_intervals_to_pool = 50)
}

survHC documentation built on Oct. 22, 2026, 1:06 a.m.