Fast Closed-Form Hazard Ratio Estimation via the Pike-Halley Estimator

Estimates the hazard ratio for a two-group parallel trial using the Pike-Halley Estimator, a pure closed-form approximation to the Cox partial likelihood maximizer. The function returns the point estimate, its standard error on the log scale, and a Wald-type confidence interval, using output names consistent with summary(survival::coxph(...)). The C++ backend accepts pooled sorted vectors directly, performing group splitting and all accumulation in a single C++ pass without intermediate R-level vector copies.

Usage

coxph_fast(
  time,
  event,
  group,
  control,
  side = 2,
  conf.level = 0.95,
  presorted = FALSE,
  strata = NULL,
  data = NULL
)

Arguments

time: A numeric vector of follow-up times for all subjects (pooled over both groups). Alternatively, a formula Surv(time, event) ~ group, optionally with strata() terms, with the variables in data; see formula_interface.
event: An integer or numeric vector of event indicators (1 = event, 0 = censored), aligned with time.
group: A vector of group labels aligned with time. Any type that supports equality comparison is accepted.
control: A scalar value indicating which level of group represents the control group. Subjects with group != control are treated as the treatment group. With the formula interface it defaults to the first level of the group variable.
side: 1 for a one-sided test in the direction of treatment benefit (hazard ratio below 1, i.e. a negative coefficient) or 2 for a two-sided test (default 2). The returned vector has no p-value: side is stored as an attribute, and the print() method reports the p-value that follows it. The confidence interval is always two-sided at conf.level.
conf.level: A single numeric value in (0, 1) specifying the confidence level for the Wald interval. Defaults to 0.95.
presorted: A logical value. If TRUE, time, event, and group are assumed to be already sorted in ascending order of time (with strata, sorted by stratum and by time within stratum, so that the rows of each stratum are contiguous), and the internal order() call is skipped. If FALSE (default), sorting is handled internally. The order is checked, and an error is given when it does not hold.
strata: An optional vector of stratum labels aligned with time, without missing values. When supplied, the stratified estimator described in Details is computed. Several stratification factors can be combined with interaction(). NULL (default) gives the unstratified estimator.
data: A data frame in which the variables of a formula given as time are evaluated. Used only with the formula interface; see formula_interface.

Value

An object of class "coxph_fast", which is a named numeric vector of length 5 with elements matching the column names of summary(coxph(...))$coefficients and summary(coxph(...))$conf.int:

coef: Log hazard ratio log(theta_hat).
exp(coef): Hazard ratio theta_hat (point estimate).
se(coef): Standard error of coef on the log scale, equal to 1 / sqrt(I_0).
lower .95: Lower bound of the Wald confidence interval for the hazard ratio. The label reflects conf.level (e.g., "lower .90" when conf.level = 0.90).
upper .95: Upper bound of the Wald confidence interval.

With strata, the number of strata is stored in the attribute strata. Returns a vector of NA_real_ values (still with class "coxph_fast") when the estimate cannot be computed (e.g., no events, all events in one group, or I_0 = 0). The object also has the class "fastsurv", whose as.data.frame(), coef(), and confint() methods are described in as.data.frame.fastsurv.

Details

Let t_k (k = 1, ..., K) denote the distinct observed event times in the pooled sample. At each t_k, let n_T_k and n_C_k be the numbers at risk in the treatment and control groups just before t_k, and let O_T_k and O_C_k be the numbers of events in each group, with n_k = n_T_k + n_C_k and O_k = O_T_k + O_C_k. Define E_T = sum n_T_k O_k / n_k and E_C = sum n_C_k O_k / n_k as the log-rank expected event totals.

The Pike-Halley Estimator is obtained in three steps. First, the Pike anchor is computed as theta_0 = (O_T E_C) / (O_C E_T). Second, the score U_0, the observed information I_0, and the third-order curvature term J_0 of the Breslow partial likelihood are evaluated at eta_0 = log(theta_0):

p_k = n_T_k theta_0 / (n_C_k + n_T_k theta_0) U_0 = sum (O_T_k - O_k p_k) I_0 = sum O_k p_k (1 - p_k) J_0 = sum O_k p_k (1 - p_k) (1 - 2 p_k)

Third, a single Halley correction to the score is applied, expanded to second order so that it is in closed form (a third-order correction):

delta_hat = U_0 / I_0 - J_0 U_0^2 / (2 I_0^3) theta_hat = theta_0 exp(delta_hat)

This correction converges cubically, so the residual error |log(theta_hat) - log(theta_Cox)| is of the order of the cube of the error of the Pike anchor. Near the null hypothesis the anchor is already close to the Cox estimate and the residual error is negligible. At a fixed hazard ratio away from 1 the error of the Pike anchor relative to the Cox estimate does not vanish as the sample size grows, so the residual error levels off at a small value that depends on the hazard ratio and the censoring pattern instead of decreasing with the sample size.

The Wald standard error on the log scale is SE = 1 / sqrt(I_0), where I_0 is the observed information evaluated at the Pike anchor. The Wald confidence interval reported by summary(coxph(...)) uses the observed information at the maximum likelihood estimate instead. The two differ by an amount of the order of the error of the Pike anchor, which is small when the hazard ratio is not far from 1; Berry, Kitchin, and Mock (1991) found the bias of the Pike estimator minimal for hazard ratios below 3.

With strata, the estimator targets the stratified Cox model, in which each stratum has its own baseline hazard and the hazard ratio is common to all strata. The at-risk sets are formed within each stratum, the observed and expected totals O_T, O_C, E_T, and E_C are summed over strata to give the stratified Pike anchor, and U_0, I_0, and J_0 are those of the stratified Breslow partial likelihood (sums of the per-stratum terms). The Halley correction and the Wald interval are then applied unchanged, so the result approximates coxph(Surv(time, event) ~ group + strata(s), ties = "breslow"). With a fixed number of strata the approximation error shrinks at the same rate as in the unstratified case; it is larger when the strata are small and the hazard ratio is far from 1, because the stratified anchor is then further from the maximum likelihood estimate.

The C++ core (pihe_core) accepts the pooled sorted data together with an integer group indicator and performs group splitting, at-risk counting, and per-distinct-event-time accumulation in a single left-to-right scan. This eliminates the rev(cumsum(rev(...))), tapply(), which(), diff(), and group-split vector copies present in the pure-R version.

The returned object has class "coxph_fast" and is a named numeric vector of length 5. A print() method formats the result similarly to summary(coxph(...)).

References

Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society. Series B (Methodological), 34(2), 187-202.

Berry, G., Kitchin, R. M., & Mock, P. A. (1991). A comparison of two simple hazard ratio estimators based on the logrank test. Statistics in Medicine, 10(5), 749-755.

Examples

library(survival)

# Compare coxph_fast with coxph on the ovarian dataset.
# coxph() treats rx as numeric with rx=1 as the reference (control),
# so set control = 1 for a consistent comparison.
fit_fast <- coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1)
fit_fast
#> Pike-Halley estimator for the hazard ratio (two-group)
#> 
#>   control = 1
#>   alternative = two.sided
#> 
#> Coefficients:
#>          coef exp(coef) se(coef)      z Pr(>|z|)
#> group -0.5964    0.5508   0.5868 -1.016     0.31
#> 
#> Hazard ratio and 95% Wald confidence interval:
#>       exp(coef) exp(-coef) lower .95 upper .95
#> group    0.5508     1.8155    0.1744    1.7399

fit_cox <- summary(coxph(Surv(futime, fustat) ~ rx, data = ovarian))
cat("coxph_fast HR :", fit_fast["exp(coef)"], "\n")
#> coxph_fast HR : 0.5508019 
cat("coxph      HR :", fit_cox$coefficients[, "exp(coef)"], "\n")
#> coxph      HR : 0.5508019 

# Stratified by residual disease, compared with coxph(... + strata())
coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1,
           strata = ovarian$resid.ds)
#> Stratified Pike-Halley estimator for the hazard ratio (two-group, 2 strata)
#> 
#>   control = 1
#>   alternative = two.sided
#> 
#> Coefficients:
#>          coef exp(coef) se(coef)      z Pr(>|z|)
#> group -0.6626    0.5155   0.5945 -1.115    0.265
#> 
#> Hazard ratio and 95% Wald confidence interval:
#>       exp(coef) exp(-coef) lower .95 upper .95
#> group    0.5155     1.9399    0.1608    1.6529
coef(coxph(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian,
           ties = "breslow"))
#>        rx 
#> -0.662632 

# presorted = TRUE: sort once outside, reuse inside a loop
ord <- order(ovarian$futime)
coxph_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord],
           control = 1, presorted = TRUE)
#> Pike-Halley estimator for the hazard ratio (two-group)
#> 
#>   control = 1
#>   alternative = two.sided
#> 
#> Coefficients:
#>          coef exp(coef) se(coef)      z Pr(>|z|)
#> group -0.5964    0.5508   0.5868 -1.016     0.31
#> 
#> Hazard ratio and 95% Wald confidence interval:
#>       exp(coef) exp(-coef) lower .95 upper .95
#> group    0.5508     1.8155    0.1744    1.7399

# Formula interface, and the result as a data frame (log hazard ratio)
fit_f <- coxph_fast(Surv(futime, fustat) ~ rx, data = ovarian)
as.data.frame(fit_f)
#>   term estimate std.error  conf.low conf.high statistic   p.value
#> 1 coef -0.59638 0.5868442 -1.746574 0.5538135 -1.016249 0.3095107
exp(confint(fit_f))
#>          2.5 %   97.5 %
#> coef 0.1743704 1.739875

# \donttest{
# Speed comparison against coxph()
if (requireNamespace("microbenchmark", quietly = TRUE)) {
  microbenchmark::microbenchmark(
    coxph_fast = coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2),
    coxph      = coxph(Surv(futime, fustat) ~ rx, data = ovarian),
    times = 1000
  )
}
#> Unit: microseconds
#>        expr     min       lq      mean  median      uq      max neval cld
#>  coxph_fast  27.191  31.6875  38.75332  40.871  44.942   94.301  1000  a 
#>       coxph 686.744 715.5070 752.97335 726.068 739.308 4688.872  1000   b
# }