| Title: | Fast Survival Analysis and Simulation for Clinical Trials |
|---|---|
| Description: | Provides fast alternatives to standard survival analysis functions in the 'survival' package, together with tools for time-to-event trial simulation and sequential analysis. The estimation and testing functions cover a single-time-point Kaplan-Meier estimator (survfit_fast()), log-rank tests including weighted and stratified variants (survdiff_fast()), a closed-form hazard ratio estimator based on the Pike-Halley Estimator method (coxph_fast()), restricted mean survival time (rmst_fast()), window mean survival time (wmst_fast()), milestone survival comparison (milestone_fast()), median survival time (medsurv_fast()), the max-combo test (maxcombo_fast()), the robust modestly-weighted log-rank test (rmw_fast()), the weighted Kaplan-Meier (Pepe-Fleming) test (wkm_fast()), the average hazard with survival weight (ahsw_fast()), and the Kalbfleisch-Prentice average hazard ratio (ahr_fast()). Besides vectors, these functions accept a model formula with a data frame, and their results can be returned as data frames. The simulation layer generates individual patient data (simdata_fast()), determines the calendar time of each analysis from combined event and calendar-time rules (cutoff_fast()), performs interim or sequential analyses (analysis_fast()), compares several experimental arms with a shared control (pairwise_fast()), applies treatment switching at progression or after an interim analysis (switch_fast()), and aggregates operating characteristics (simsummary_fast()). A visualization layer assembles design-stage scenarios (gen_scenario_fast()) and builds analysis-stage Kaplan-Meier curves (kmcurve_fast()), each with plot and print methods. All functions are designed for repeated evaluation inside large simulation loops, such as adaptive sample-size re-estimation, probability-of-success calculations, and regional consistency evaluation in multi-regional trials. Core computations are implemented in 'C++' via 'Rcpp' for maximum performance. Methodological background is described in Collett (2014, ISBN:9781439856789). |
| Authors: | Gosuke Homma [aut, cre] |
| Maintainer: | Gosuke Homma <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.2.0 |
| Built: | 2026-11-04 09:29:14 UTC |
| Source: | https://github.com/gosukehommaex/fastsurvival |
Estimates the average hazard ratio of Kalbfleisch and Prentice (1981) between
two groups over the time interval from 0 to tau, based on the
Kaplan-Meier estimator of each group's survival function. This is the
Kaplan-Meier (unweighted) special case of the estimator implemented in the
archived AHR package, restricted to two groups and recoded in C++ for
use inside simulation loops.
ahr_fast( time, event, group, control, side = 2, conf.level = 0.95, tau = NULL, null.ahr = 1, presorted = FALSE, data = NULL )ahr_fast( time, event, group, control, side = 2, conf.level = 0.95, tau = NULL, null.ahr = 1, presorted = FALSE, data = NULL )
time |
vector of right-censored event times
Alternatively, a formula |
event |
0/1 (or logical) event indicators, 1 for an event |
group |
vector with exactly two distinct values identifying the groups. |
control |
the value of |
side |
1 for a one-sided test in the direction of treatment benefit
(average hazard ratio below |
conf.level |
confidence level for the confidence interval (default 0.95) |
tau |
upper limit of the interval over which the average hazard ratio is
computed. If |
null.ahr |
value of the average hazard ratio under the null hypothesis used for the Z statistic and p-value (default 1) |
presorted |
if |
data |
A data frame in which the variables of a formula given as
|
The estimator works with the group shares of the total hazard. Writing
S1 and S2 for the two survival functions, the reference-group
share is theta1 = -integral(S2 dS1) / (1 - S1(tau) S2(tau)) over
[0, tau], the comparison-group share is theta2 = 1 - theta1,
and the average hazard ratio is ahr = theta2 / theta1. Under
proportional hazards with hazard ratio psi, ahr estimates
psi. A value above 1 indicates higher hazard (worse survival) in the
comparison group. The variance of theta1 is the direct
Greenwood-based estimator. The primary test is on the theta (group-share)
scale, as in Kalbfleisch and Prentice (1981) and Dormuth et al. (2024,
eq. 5): the comparison-group share is compared with its null value (0.5 when
null.ahr = 1). An equivalent test and a confidence interval on the
log(ahr) scale are also reported.
The integral uses the right-continuous Kaplan-Meier estimate S2 at
each jump of S1, exactly as the AHR package does, and
theta2 is obtained as 1 - theta1 rather than by a separate
integral. When both groups have events at the same time (tied event times,
as with data recorded in whole days or months), the two shares computed this
way are not symmetric in the groups, so exchanging control does not
give exactly the reciprocal average hazard ratio. With continuous times
there are no such ties and the estimate is symmetric. Specify
control as the actual reference group so that the result matches the
AHR package with that group as the reference.
This is distinct from ahsw_fast, which estimates the
Uno-Horiguchi average hazard with survival weight.
An object of class "ahr_fast", a list with elements
ahr (the average hazard ratio, comparison vs reference),
log.ahr, se.loghr, lower, upper,
conf.level, z and p.value (the primary test on the
theta / group-share scale, as in Dormuth et al. 2024 eq. 5),
z.loghr and p.value.loghr (the equivalent test on the
log(ahr) scale), se.theta (standard error of the tested
comparison-group share), null.share, null.ahr, theta
(the two group shares), var.theta1, var.theta2,
side, tau, n (the two group sizes) and groups.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Kalbfleisch, J. D., & Prentice, R. L. (1981). Estimation of the average hazard ratio. Biometrika, 68(1), 105-112.
Dormuth, I., Pauly, M., Rauch, G., & Herrmann, C. (2024). Sample size calculation under nonproportional hazards using average hazard ratios. Biometrical Journal, 66(6), e202300271.
set.seed(1) n <- 200 time1 <- rexp(n, 0.1) time2 <- rexp(n, 0.18) cens <- rexp(2 * n, 0.05) obs <- pmin(c(time1, time2), cens) event <- as.integer(c(time1, time2) <= cens) group <- rep(c(0, 1), each = n) ahr_fast(obs, event, group, control = 0, tau = 8) # Formula interface d <- data.frame(obs, event, group) ahr_fast(Surv(obs, event) ~ group, data = d, tau = 8)set.seed(1) n <- 200 time1 <- rexp(n, 0.1) time2 <- rexp(n, 0.18) cens <- rexp(2 * n, 0.05) obs <- pmin(c(time1, time2), cens) event <- as.integer(c(time1, time2) <= cens) group <- rep(c(0, 1), each = n) ahr_fast(obs, event, group, control = 0, tau = 8) # Formula interface d <- data.frame(obs, event, group) ahr_fast(Surv(obs, event) ~ group, data = d, tau = 8)
Computes the average hazard with survival weight (AHSW) of Uno and Horiguchi
for two groups and the between-group contrasts. The average hazard on the
window from 0 to tau is the ratio of the cumulative event probability
at tau to the restricted mean survival time at tau, both based
on the Kaplan-Meier estimate. The function returns the per-group average
hazard, the ratio of average hazards (RAH, treatment over control) on the log
scale, and the difference of average hazards (DAH, treatment minus control)
on the identity scale, each with a confidence interval and a test (one- or
two-sided, following side).
The C++ backend walks the pooled sorted data once per group, so the function
is suitable for simulation loops with presorted = TRUE.
ahsw_fast( time, event, group, control, side = 2, conf.level = 0.95, tau, presorted = FALSE, data = NULL )ahsw_fast( time, event, group, control, side = 2, conf.level = 0.95, tau, presorted = FALSE, data = NULL )
time |
A numeric vector of follow-up times for all subjects.
Alternatively, a formula |
event |
An integer or numeric vector of event indicators
(1 = event, 0 = censored), aligned with |
group |
A vector of group labels aligned with |
control |
A scalar value indicating which level of |
side |
1 for a one-sided test in the direction of treatment benefit
(lower average hazard in the treatment group, i.e. a ratio below 1 and a
negative difference) or 2 for a two-sided test (default 2). The reported
p-values follow this choice; the confidence intervals are always two-sided
at |
conf.level |
A single numeric value in (0, 1) specifying the confidence level. Defaults to 0.95. |
tau |
A single positive numeric value, the truncation time point for the average hazard. A warning is given when it exceeds the largest observed time of either group, because the Kaplan-Meier curve is not estimated beyond that time. |
presorted |
A logical value. If |
data |
A data frame in which the variables of a formula given as
|
The average hazard with survival weight is
AH(tau) = (1 - S(tau)) / integral over [0, tau] of S(u) du,
the ratio of the cumulative event probability at tau to the restricted
mean survival time at tau. It can be read as a general censoring-free
incidence rate on the window from 0 to tau and stays interpretable
under non-proportional hazards. This is a different quantity from the average
hazard ratio of Kalbfleisch and Prentice (see ahr_fast), which
averages the time-varying ratio of hazards rather than forming a single
average hazard per group and then contrasting.
Writing the treatment and control average hazards as a1 and a0, the ratio contrast is RAH = a1 / a0, formed on the log scale with variance v_Q1 / n1 + v_Q0 / n0, and the difference contrast is DAH = a1 - a0, with variance v_U1 / n1 + v_U0 / n0, using the independence of the two groups. The per-group variance terms v_Q (log scale) and v_U (identity scale) follow the asymptotic variance of Uno and Horiguchi, computed from the Nelson-Aalen increments, the running restricted mean survival time and the at-risk fraction. The confidence interval for RAH is exponentiated from the log scale, and the p-values are based on the normal approximation.
When presorted = TRUE, the inputs are assumed to be sorted in
ascending order of time, so the internal order() call is
skipped. Splitting into groups preserves the ascending order within each
group.
An object of class "ahsw_fast", a named numeric vector
containing the per-group average hazards (ah.ctrl, ah.trt),
the ratio contrast (rah, rah.lower, rah.upper,
p.rah), and the difference contrast (dah, dah.lower,
dah.upper, p.dah). The truncation time and confidence level
are stored as attributes tau and conf.level, and the
control label is also stored. Returns NA values (still with
class "ahsw_fast") when either group has no events up to tau
(a zero average hazard) or a non-finite variance.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Uno, H., & Horiguchi, M. (2023). Ratio and difference of average hazard with survival weight: new measures to quantify survival benefit of new therapy. Statistics in Medicine, 42(7), 936-952.
rmst_fast for the restricted mean survival time, and
survfit_fast for the Kaplan-Meier estimate at a time point.
library(survival) # Average hazard contrasts on the ovarian data ahsw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, tau = 600) # presorted = TRUE: sort once outside, reuse inside a loop ord <- order(ovarian$futime) ahsw_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord], control = 1, tau = 600, presorted = TRUE) # Formula interface ahsw_fast(Surv(futime, fustat) ~ rx, data = ovarian, tau = 600) # Validation against survAH if (requireNamespace("survAH", quietly = TRUE)) { arm <- as.numeric(ovarian$rx == 2) survAH::ah2(time = ovarian$futime, status = ovarian$fustat, arm = arm, tau = 600) }library(survival) # Average hazard contrasts on the ovarian data ahsw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, tau = 600) # presorted = TRUE: sort once outside, reuse inside a loop ord <- order(ovarian$futime) ahsw_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord], control = 1, tau = 600, presorted = TRUE) # Formula interface ahsw_fast(Surv(futime, fustat) ~ rx, data = ovarian, tau = 600) # Validation against survAH if (requireNamespace("survAH", quietly = TRUE)) { arm <- as.numeric(ovarian$rx == 2) survAH::ah2(time = ovarian$futime, status = ovarian$fustat, arm = arm, tau = 600) }
Performs interim or sequential analyses of simulated two-group
time-to-event data at one or more analysis times ("looks"). Each look is
defined by a target cumulative number of events (information-based timing),
by a calendar time (calendar-based timing), or by a calendar cutoff given
separately for every simulated trial (for example by
cutoff_fast).
At every look the data are administratively censored at the corresponding
calendar cutoff, and the requested statistics are computed for each
simulated trial by reusing survdiff_fast,
coxph_fast, rmst_fast, survfit_fast,
maxcombo_fast, ahsw_fast,
milestone_fast, rmw_fast,
ahr_fast, medsurv_fast,
wkm_fast, and wmst_fast. Optionally the same
statistics are also reported within each subgroup. The censoring and time
sorting are handled by a C++ backend, and the analysis cores are called with
presorted = TRUE, so each look avoids a redundant sort.
analysis_fast( data, control, event.looks = NULL, time.looks = NULL, stat = "logrank", tau = NULL, t.eval = NULL, conf.level = 0.95, side = 2, by.subgroup = FALSE, weight = c("logrank", "fh", "mwlrt", "gehan", "tarone-ware"), rho = 0, gamma = 0, t_star = NULL, strata = NULL, ms.method = c("wald", "loglog", "mover"), s_star = 0.5, mc.rho = c(0, 0, 1, 1), mc.gamma = c(0, 1, 0, 1), abseps = 1e-05, maxpts = 25000, medsurv.method = c("km", "nph"), medsurv.bw = NULL, wkm.weight = c("PF", "sqrtPF", "constant"), wmst.tau1 = 0, wmst.tau2 = NULL, cutoff.looks = NULL, mc.alpha = NULL )analysis_fast( data, control, event.looks = NULL, time.looks = NULL, stat = "logrank", tau = NULL, t.eval = NULL, conf.level = 0.95, side = 2, by.subgroup = FALSE, weight = c("logrank", "fh", "mwlrt", "gehan", "tarone-ware"), rho = 0, gamma = 0, t_star = NULL, strata = NULL, ms.method = c("wald", "loglog", "mover"), s_star = 0.5, mc.rho = c(0, 0, 1, 1), mc.gamma = c(0, 1, 0, 1), abseps = 1e-05, maxpts = 25000, medsurv.method = c("km", "nph"), medsurv.bw = NULL, wkm.weight = c("PF", "sqrtPF", "constant"), wmst.tau1 = 0, wmst.tau2 = NULL, cutoff.looks = NULL, mc.alpha = NULL )
data |
A data frame from |
control |
A scalar value indicating which level of |
event.looks |
A vector of positive whole numbers, the target cumulative
event counts, one per look. Mutually exclusive with |
time.looks |
A numeric vector of calendar times, one per look.
Mutually exclusive with |
stat |
A character vector naming the statistics to compute. Any subset
of |
tau |
A single positive numeric value, the restriction horizon for
|
t.eval |
A single positive numeric value, the landmark time for
|
conf.level |
A single numeric value in (0, 1), the confidence level of
all the confidence-interval columns (and of the |
side |
An integer, either 1 or 2. When |
by.subgroup |
A logical value. When |
weight |
A character string naming the weight scheme for the
|
rho |
A numeric Fleming-Harrington first parameter, used only when
|
gamma |
A numeric Fleming-Harrington second parameter, used only when
|
t_star |
A single non-negative numeric value, the timepoint of the
modestly-weighted log-rank test. Required only when |
strata |
An optional character vector naming one or more subgroup
columns of |
ms.method |
A character string naming the inference method for the
|
s_star |
A single numeric value in (0, 1], the survival-probability
threshold of the modestly-weighted component of the |
mc.rho |
A numeric vector of Fleming-Harrington first parameters for the
|
mc.gamma |
A numeric vector of Fleming-Harrington second parameters for
the |
abseps |
A single positive numeric value, the absolute error tolerance
passed to the multivariate normal integration of the |
maxpts |
A single positive integer, the maximum number of function
evaluations for the quasi-Monte-Carlo integration used by the
|
medsurv.method |
A character string naming the variance method for the
|
medsurv.bw |
An optional single positive numeric value, the kernel
bandwidth for the median hazard used by |
wkm.weight |
A character string naming the weight for the |
wmst.tau1 |
A single non-negative numeric value, the lower window limit
for the |
wmst.tau2 |
A single positive numeric value, the upper window limit for
the |
cutoff.looks |
A numeric matrix of per-simulation calendar cutoffs with
one row per simulated trial and one column per look, such as the output of
|
mc.alpha |
An optional numeric vector of nominal levels in (0, 1), one
per look or a single value for all looks, at which the |
The input data is the data frame returned by
simdata_fast for a two-group trial. The columns sim,
group, accrual_time, tte, and event are
required and must not contain missing values. The tte column must be
numeric and non-negative, event must be numeric or logical and coded
as 0 or 1, and group must have exactly two distinct values, one of
which is control. For a
multi-arm trial, subset the data to two arms first or use
pairwise_fast.
For a look at calendar time cutoff, each subject with accrual time
a contributes only if enrolled by then (a <= cutoff). The
observed time at the look is min(tte, cutoff - a) and the observed
event indicator is the original event when the natural event or
dropout occurred on or before cutoff (a + tte <= cutoff), and
zero otherwise (administrative censoring at the look).
When event.looks is supplied, the calendar cutoff for a target of
d events is the calendar time of the d-th event in that
simulated trial, counted over the whole trial population. If a simulation
contains fewer than d events, the target is never reached: the full
data are used, so the statistics are those of the final data, reached
is FALSE, and cutoff is NA. When several events occur
at the calendar time of the d-th event, all of them are included, so
the number of events analyzed can exceed the target. When time.looks is supplied, the cutoff is the specified
calendar time and reached is always TRUE. In both cases the
cutoff is determined once on the whole population and then used for the
overall analysis and for every subgroup analysis at that look.
When cutoff.looks is supplied, the calendar cutoff of each look is
given separately for every simulated trial, typically by
cutoff_fast. This analyzes an endpoint at cutoffs determined
by another endpoint (for example overall survival at the event-driven looks
of progression-free survival), by a subset of the subjects, or by a
combination of event and calendar-time rules. A missing or infinite entry
marks a look whose trigger was not met in that simulation: as for an
unreached event.looks target, the full data are used, reached
is FALSE, and cutoff is NA. In contrast,
pairwise_fast reports NA counts and statistics for an
unreached look, because its contrasts share one cutoff and have no common
final data.
Exactly one of event.looks, time.looks, and
cutoff.looks must be supplied.
The statistics are selected with stat, which may name one or more of
"logrank", "coxph", "rmst", "km",
"maxcombo", "ahsw", "milestone", "rmw",
"ahr", "medsurv", "wkm", and "wmst". The
references for each statistic are given on the help page of the
corresponding function listed under See Also.
The "logrank" statistic is configurable. By default it is the
ordinary unweighted, unstratified two-group log-rank test and reproduces the
behavior of earlier versions of this function exactly. A non-default
weight selects a weighted log-rank test (Fleming-Harrington,
modestly-weighted, Gehan-Breslow, or Tarone-Ware) for non-proportional
hazards, and a non-NULL strata selects the stratified test,
summing the per-stratum contributions. The two options combine to give the
stratified weighted log-rank test. Whatever configuration is chosen, the
result is written to the same logrank.z, logrank.chisq, and
logrank.p columns; the unweighted unstratified case is the K = 1
single-stratum degenerate form of the general statistic. The columns
rho, gamma, and t_star parametrize weight, and
strata names one or more subgroup columns of data used as the
stratification variable. The stratification is determined on the whole cut
data, independently of the population marginalization, so a stratified
overall analysis is the canonical primary test. When strata names
the subgroup column that defines a subgroup population, the stratum is
constant within that population and the stratified test there reduces to
the ordinary one. Several columns define their combinations as strata.
The same strata also stratify the
"coxph" statistic, which is then the stratified Pike-Halley estimate
of coxph_fast (a common hazard ratio with a separate baseline
hazard in each stratum).
The "maxcombo" statistic is the max-combo test of mc.rho and
mc.gamma Fleming-Harrington weights. Its maxcombo.stat is the
most extreme component (min of the component Z-scores when
side = 1, so that a negative value favors treatment, and the maximum
absolute component when side = 2), and maxcombo.p is the joint
multivariate-normal p-value, which already follows side. As in
maxcombo_fast, the p-value of a one-sided test with two or
three weights is computed by the deterministic TVPACK algorithm, and
otherwise by the randomized quasi-Monte-Carlo GenzBretz algorithm of
mvtnorm, which uses R's random-number generator; call set.seed()
beforehand for p-values that are reproducible to the last digit. The
correlation matrix of the default weights is singular, and the GenzBretz
p-values then carry an integration error of the order of 1e-4 with the
default maxpts (see maxcombo_fast).
The multivariate normal integration dominates the computing time of the
"maxcombo" statistic in a large simulation. When only the decision
at a nominal level matters, mc.alpha avoids most of these
integrals. With K weights and p_min the smallest
component p-value (pnorm(m) for side = 1 and
2 * pnorm(-m) for side = 2, where m is
maxcombo.stat), the max-combo p-value satisfies
p_min <= p <= min(1, K * p_min) (the Bonferroni inequality). When
p_min > mc.alpha, p_min is reported; when
K * p_min <= mc.alpha, K * p_min is reported; otherwise the
p-value is computed by integration. A bound is reported only when it is on
the same side of mc.alpha as the exact p-value, so the rejection
decisions maxcombo.p <= mc.alpha (as in simsummary_fast
with alpha = mc.alpha) that the bounds settle are those of the exact
p-values, and the reported value is a bound unless
maxcombo.p.exact is TRUE. The remaining p-values are
integrated as without mc.alpha. When the randomized GenzBretz
algorithm is used, the random numbers of each integral change with the
number of integrals, so near the level a decision can differ from that of
a run without mc.alpha, as it can between two runs with different
seeds.
Use exact p-values (the default mc.alpha = NULL) when the p-values
themselves, or decisions at other levels such as Bonferroni-adjusted levels,
are needed.
The "ahsw" statistic is the average hazard with survival weight of Uno
and Horiguchi on the window from 0 to tau. It reports the per-group
average hazards, the ratio (RAH) and difference (DAH) contrasts with their
confidence intervals, and p-values for both contrasts that follow
side, matching ahsw_fast. As there, a group without
an event up to tau has an average hazard of 0, and the contrasts
are then NA. The benefit direction is a
ratio below 1 (a negative log ratio) and a negative difference, so the
one-sided p-values are lower-tail probabilities.
The "milestone" statistic compares the Kaplan-Meier survival
probabilities of the two groups at the milestone timepoint tau. It
reports the per-group survival, the difference (treatment minus control) with
its confidence interval, the test statistic, and the p-value. The inference
method is selected with ms.method ("wald", "loglog", or
"mover"), matching milestone_fast. The benefit direction
is a positive difference (higher treatment survival). For "wald" and
"mover" a positive Z favors treatment; for "loglog" the
statistic is the difference of the complementary log-log transforms, so a
negative Z favors treatment, and the one-sided p-value is the lower tail as
in milestone_fast. The "mover" statistic is obtained
by inverting the MOVER interval at conf.level, so its p-value depends
on conf.level, consistently with the interval (see
milestone_fast). When the Kaplan-Meier estimate of a group
is 0 or 1 at tau, its one-sample interval degenerates to the
estimate, so the interval of the difference is still reported, while the
"loglog" statistic is NA.
The "rmw" statistic is the robust modestly-weighted log-rank test of
Magirr and Öhrn, the maximum of the standard log-rank component
and a single modestly-weighted component with survival-threshold
s_star. Its
rmw.stat is the most extreme standardized component (the minimum when
side = 1, so a negative value favors treatment, and the maximum
absolute component when side = 2), and rmw.p is the joint
bivariate-normal p-value, which already follows side, matching
rmw_fast.
The "ahr" statistic is the Kalbfleisch-Prentice average hazard ratio
over the window from 0 to tau. It reports the average hazard ratio
(treatment relative to control), the two group shares of the total hazard,
the test statistic on the group-share scale, and the p-value. The benefit
direction is an average hazard ratio below 1 (a negative Z), as in
ahr_fast.
The "medsurv" statistic compares the non-parametric Kaplan-Meier
median survival times of the two groups. It reports the per-group medians,
their difference (treatment minus control) with a plain-scale confidence
interval, the Wald statistic, and the p-value. The variance method is
selected with medsurv.method ("km" or "nph"), matching
medsurv_fast; medsurv.bw optionally overrides the kernel
bandwidth used by medsurv.method = "km". The benefit direction is a
positive difference (a longer median under treatment), so a positive Z favors
treatment. When the Kaplan-Meier estimate of a group stays at 0.5 up to an
infinite observed time, as for subjects with tte = Inf at a look
that is not reached, the median of that group is not defined and the
"medsurv" columns that depend on it are NA.
The "wkm" statistic is the weighted Kaplan-Meier (Pepe-Fleming) test.
It reports the weighted integrated survival difference (treatment minus
control) with its confidence interval, the Wald statistic, and the p-value.
The weight is selected with wkm.weight ("PF", "sqrtPF",
or "constant"), matching wkm_fast. The benefit direction
is a positive weighted difference (a positive Z favors treatment). The
integral runs over the observed range, so all the "wkm" columns are
NA when an observed time is infinite, as for a subject without a
finite event or dropout time (tte = Inf) at a look that is not
reached.
The "wmst" statistic compares the window mean survival time, the area
under the Kaplan-Meier curve between wmst.tau1 and wmst.tau2.
It reports the per-group WMST values, their difference (treatment minus
control) with its confidence interval, the Wald statistic, and the p-value,
matching wmst_fast. With wmst.tau1 = 0 the window mean
survival time is the restricted mean survival time at wmst.tau2. The
benefit direction is a positive difference (a positive Z favors treatment).
When by.subgroup = TRUE, the output is given in long form with a
population column. Each (sim, look) produces one row for the
whole trial (population = "overall") plus one row per subgroup level
of each subgroup factor in data. Subgroup factors are the columns
named subgroup or subgroup1, subgroup2, and so on. Each
factor is marginalized separately, so a factor with levels 1 and 2 yields
the populations subgroup_1 and subgroup_2 (or
subgroup1_1, subgroup1_2, and so on for numbered factors).
The look cutoff is always determined on the whole population, so subgroup
rows at a given look share the same cutoff, reached, and
look.value; their n.enrolled and n.event are the counts
within that subgroup. When by.subgroup = FALSE (default) the output
has no population column and one row per (sim, look), matching
the whole-population analysis.
A data frame. When by.subgroup = FALSE, it has
nsim * length(looks) rows. When by.subgroup = TRUE, it has
nsim * length(looks) * (1 + total subgroup levels) rows and an
extra population column placed after look.value. The common
columns are sim, look (1-based look index),
look.value (the requested event count or calendar time, or the
"look.value" attribute of cutoff.looks), optionally population, cutoff (the
calendar time used, NA when the look was not reached),
reached, n.enrolled, n.event, n.dropout (the
number of enrolled subjects whose dropout occurred on or before the cutoff)
and n.pipeline (n.enrolled - n.event - n.dropout, the
subjects still in follow-up at the cutoff), followed by the
columns of the requested statistics. A statistic that cannot be computed
for a row (no events, or an empty group) is NA. The statistic
columns are logrank.z, logrank.chisq, and logrank.p
for "logrank"; cox.coef, cox.hr, cox.se,
cox.z, cox.p, cox.lower, and cox.upper for
"coxph"; rmst.ctrl, rmst.trt, rmst.diff,
rmst.diff.lower, rmst.diff.upper, rmst.z, and
rmst.p for "rmst"; km.surv.ctrl and
km.surv.trt for "km"; maxcombo.stat and
maxcombo.p (and maxcombo.p.exact when mc.alpha is
supplied) for "maxcombo"; and ahsw.ah.ctrl,
ahsw.ah.trt, ahsw.rah, ahsw.rah.lower,
ahsw.rah.upper, ahsw.p.rah, ahsw.dah,
ahsw.dah.lower, ahsw.dah.upper, and ahsw.p.dah for
"ahsw"; milestone.surv.ctrl, milestone.surv.trt,
milestone.diff, milestone.diff.lower,
milestone.diff.upper, milestone.z, and milestone.p
for "milestone"; rmw.stat and rmw.p for "rmw";
and ahr.ahr, ahr.theta.ctrl, ahr.theta.trt,
ahr.z, and ahr.p for "ahr";
medsurv.ctrl, medsurv.trt, medsurv.diff,
medsurv.diff.lower, medsurv.diff.upper, medsurv.z, and
medsurv.p for "medsurv"; wkm.wdiff, wkm.lower,
wkm.upper, wkm.z, and wkm.p for "wkm"; and
wmst.ctrl, wmst.trt, wmst.diff,
wmst.diff.lower, wmst.diff.upper, wmst.z, and
wmst.p for "wmst". The Z columns
logrank.z, cox.z, and
rmst.z carry the natural sign of each test, and the p-value columns
follow side.
simdata_fast, survdiff_fast,
coxph_fast, rmst_fast,
survfit_fast, maxcombo_fast,
ahsw_fast, milestone_fast,
rmw_fast, ahr_fast,
medsurv_fast, wkm_fast,
wmst_fast, cutoff_fast,
pairwise_fast, switch_fast.
df <- simdata_fast( nsim = 50, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(list(0.10, 0.07), 0.05), prevalence = c(0.5, 0.5), seed = 1 ) # Whole-population analysis at two event-based looks res1 <- analysis_fast(df, control = 1, event.looks = c(80, 140)) head(res1) # Stratified log-rank and stratified Cox on the subgroup factor res2 <- analysis_fast(df, control = 1, time.looks = 24, stat = c("logrank", "coxph"), strata = "subgroup") head(res2) # Fleming-Harrington G(0, 1) weighted log-rank for delayed effects res3 <- analysis_fast(df, control = 1, time.looks = 24, stat = "logrank", weight = "fh", rho = 0, gamma = 1) head(res3) # Max-combo and AHSW together res4 <- analysis_fast(df, control = 1, time.looks = 24, stat = c("maxcombo", "ahsw"), tau = 18) head(res4) # Milestone survival, robust modestly-weighted, and average hazard ratio res5 <- analysis_fast(df, control = 1, time.looks = 24, stat = c("milestone", "rmw", "ahr"), tau = 18, ms.method = "loglog", s_star = 0.5, side = 1) head(res5)df <- simdata_fast( nsim = 50, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(list(0.10, 0.07), 0.05), prevalence = c(0.5, 0.5), seed = 1 ) # Whole-population analysis at two event-based looks res1 <- analysis_fast(df, control = 1, event.looks = c(80, 140)) head(res1) # Stratified log-rank and stratified Cox on the subgroup factor res2 <- analysis_fast(df, control = 1, time.looks = 24, stat = c("logrank", "coxph"), strata = "subgroup") head(res2) # Fleming-Harrington G(0, 1) weighted log-rank for delayed effects res3 <- analysis_fast(df, control = 1, time.looks = 24, stat = "logrank", weight = "fh", rho = 0, gamma = 1) head(res3) # Max-combo and AHSW together res4 <- analysis_fast(df, control = 1, time.looks = 24, stat = c("maxcombo", "ahsw"), tau = 18) head(res4) # Milestone survival, robust modestly-weighted, and average hazard ratio res5 <- analysis_fast(df, control = 1, time.looks = 24, stat = c("milestone", "rmw", "ahr"), tau = 18, ms.method = "loglog", s_star = 0.5, side = 1) head(res5)
The results of the estimation and testing functions have, besides their
own class, the class "fastsurv". Its as.data.frame() method
returns every estimate and test of a result as a row of a data frame with
the same columns for all the functions, so that results can be stacked with
rbind(). coef() returns the estimates and confint()
their confidence limits.
## S3 method for class 'fastsurv' as.data.frame(x, row.names = NULL, optional = FALSE, ...) ## S3 method for class 'fastsurv' coef(object, ...) ## S3 method for class 'fastsurv' confint(object, parm, level = NULL, ...)## S3 method for class 'fastsurv' as.data.frame(x, row.names = NULL, optional = FALSE, ...) ## S3 method for class 'fastsurv' coef(object, ...) ## S3 method for class 'fastsurv' confint(object, parm, level = NULL, ...)
x, object
|
A result of an estimation or testing function, of class
|
row.names |
|
optional |
Ignored; the column names are fixed. |
... |
Further arguments, currently ignored. |
parm |
The estimates for which limits are returned, as terms or as positions among the estimates. All estimates by default. |
level |
|
The data frame has the columns term, estimate,
std.error, conf.low, conf.high, statistic, and
p.value, with one row for each quantity in the table below. The
terms are the names of the elements of the result. The estimate, its
standard error, and its confidence limits are on the same scale. The
statistic and the p-value are those of the test that the function reports
for that quantity, following its side, as printed by the
print() method. Entries that the function does not compute are
NA.
survfit_fast |
surv |
survdiff_fast |
logrank, or the name of the
weight (test only) |
coxph_fast |
coef (the log hazard ratio) |
rmst_fast |
rmst; or rmst.ctrl,
rmst.trt, diff, log.ratio |
wmst_fast |
wmst; or wmst.control,
wmst.treatment, diff |
medsurv_fast |
median; or median.control,
median.treatment, diff |
milestone_fast |
surv.control,
surv.treatment, diff |
maxcombo_fast |
maxcombo and the component
statistics, such as FH(0,1) (test only) |
rmw_fast |
rmw, logrank, mwlrt
(test only) |
wkm_fast |
wdiff |
ahsw_fast |
ah.ctrl, ah.trt,
log.rah, dah |
ahr_fast |
log.ahr, theta
|
Ratios are reported on the log scale, on which their confidence intervals
are symmetric and their tests are computed: coef of
coxph_fast() is the log hazard ratio (as from coef() of
survival::coxph() with the group variable as a factor), and
log.ratio, log.rah, and log.ahr are the logarithms of
the ratio of restricted mean survival times, the ratio of average hazards,
and the average hazard ratio. Take
exp() of the estimate and of the confidence limits for the ratio.
rmst_fast() and ahsw_fast() do not store the standard error of
a log ratio or of the difference of average hazards; it is recovered from
the confidence limits, which are symmetric on these scales, as
(conf.high - conf.low) / (2 q), where q is the normal quantile
of the confidence level, and the statistic of log.rah and dah
is the estimate divided by this standard error.
The rows of a test only carry the statistic: the value of
survdiff_fast() (the Z-score when side = 1, the chi-square
when side = 2), the max-combo and rMW statistics with their
components, which have no p-value of their own. For milestone_fast(),
the statistic of diff is that of method, which for
"loglog" is on the complementary log-log scale; a standard error is
given only for method = "wald", whose interval is symmetric. For
ahr_fast(), theta is the treatment group's share of the total
hazard with its standard error and the primary test of
ahr_fast(), and log.ahr carries the equivalent test on the log
scale.
coef() returns the estimate column named by term, and
confint() the confidence limits as a matrix, for the rows that are
estimates. Both give an error for a test only. The limits are those computed
by the function at its conf.level; they cannot be recomputed at
another level, so level must be NULL or equal to that
conf.level.
as.data.frame() returns a data frame with the columns
described in Details. coef() returns a named numeric vector, and
confint() a matrix with one row per estimate and two columns
labeled by the percentages of the limits, as stats::confint().
library(survival) fit <- coxph_fast(Surv(futime, fustat) ~ rx, data = ovarian) as.data.frame(fit) coef(fit) exp(confint(fit)) as.data.frame(rmst_fast(Surv(futime, fustat) ~ rx, data = ovarian, tau = 600)) # Results of a loop (with the vector interface) stacked into a data frame res <- lapply(1:3, function(s) { set.seed(s) i <- sample(nrow(ovarian), replace = TRUE) fit_s <- coxph_fast(ovarian$futime[i], ovarian$fustat[i], ovarian$rx[i], control = 1) cbind(sim = s, as.data.frame(fit_s)) }) do.call(rbind, res)library(survival) fit <- coxph_fast(Surv(futime, fustat) ~ rx, data = ovarian) as.data.frame(fit) coef(fit) exp(confint(fit)) as.data.frame(rmst_fast(Surv(futime, fustat) ~ rx, data = ovarian, tau = 600)) # Results of a loop (with the vector interface) stacked into a data frame res <- lapply(1:3, function(s) { set.seed(s) i <- sample(nrow(ovarian), replace = TRUE) fit_s <- coxph_fast(ovarian$futime[i], ovarian$fustat[i], ovarian$rx[i], control = 1) cbind(sim = s, as.data.frame(fit_s)) }) do.call(rbind, res)
Returns the Kaplan-Meier curves of a kmcurve_fast object as a data
frame, one row per group and distinct observed time, for drawing the curves
with other graphics software or for tabulation.
## S3 method for class 'kmcurve_fast' as.data.frame( x, row.names = NULL, optional = FALSE, conf.level = 0.95, conf.type = "log", ... )## S3 method for class 'kmcurve_fast' as.data.frame( x, row.names = NULL, optional = FALSE, conf.level = 0.95, conf.type = "log", ... )
x |
A |
row.names |
|
optional |
Ignored; the column names are fixed. |
conf.level |
Confidence level of the pointwise limits. Defaults to 0.95. |
conf.type |
Type of the pointwise limits: |
... |
Further arguments, currently ignored. |
For each group the rows are the distinct observed times (event and
censoring times) in ascending order, preceded by a row at time 0 with
survival 1, so that a step function drawn through the rows starts at 1 and
ends at the largest observed time. At each time, n.risk is the
number of subjects with an observed time at or after it, n.event and
n.censor are the numbers of events and censored times at it, and
surv is the Kaplan-Meier estimate just after it. std.err is
the Greenwood standard error on the survival scale, the std.err
reported by summary() of a survival::survfit() object, and
lower and upper are the pointwise confidence limits of the type
conf.type, as drawn by plot.kmcurve_fast. At a time at
which every subject at risk fails, the estimate is 0 and std.err and
the limits are 0, as in survfit_fast, where summary()
reports a NaN standard error and NA limits. Apart from the
row at time 0, the times, counts, and survival estimates are those reported
by summary() of a survival::survfit() object with
censored = TRUE.
A data frame with the columns group (the group label, as
printed by print.kmcurve_fast), time,
n.risk, n.event, n.censor, surv,
std.err, lower, and upper. In the two-group case the
rows of the control group come first.
kmcurve_fast, plot.kmcurve_fast
set.seed(1) n <- 50 t0 <- rexp(n, log(2) / 12) t1 <- rexp(n, log(2) / 18) cens <- runif(2 * n, 0, 30) time <- pmin(c(t0, t1), cens) event <- as.integer(c(t0, t1) <= cens) group <- rep(c(0, 1), each = n) km <- as.data.frame(kmcurve_fast(time, event, group, control = 0)) head(km)set.seed(1) n <- 50 t0 <- rexp(n, log(2) / 12) t1 <- rexp(n, log(2) / 18) cens <- runif(2 * n, 0, 30) time <- pmin(c(t0, t1), cens) event <- as.integer(c(t0, t1) <= cens) group <- rep(c(0, 1), each = n) km <- as.data.frame(kmcurve_fast(time, event, group, control = 0)) head(km)
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.
coxph_fast( time, event, group, control, side = 2, conf.level = 0.95, presorted = FALSE, strata = NULL, data = NULL )coxph_fast( time, event, group, control, side = 2, conf.level = 0.95, presorted = FALSE, strata = NULL, data = NULL )
time |
A numeric vector of follow-up times for all subjects (pooled
over both groups).
Alternatively, a formula |
event |
An integer or numeric vector of event indicators
(1 = event, 0 = censored), aligned with |
group |
A vector of group labels aligned with |
control |
A scalar value indicating which level of |
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: |
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 |
strata |
An optional vector of stratum labels aligned with
|
data |
A data frame in which the variables of a formula given as
|
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(...)).
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:
coefLog 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 .95Lower 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 .95Upper 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.
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.
coxph for the standard iterative Cox estimator.
print.coxph_fast for the print method.
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 fit_cox <- summary(coxph(Surv(futime, fustat) ~ rx, data = ovarian)) cat("coxph_fast HR :", fit_fast["exp(coef)"], "\n") cat("coxph HR :", fit_cox$coefficients[, "exp(coef)"], "\n") # Stratified by residual disease, compared with coxph(... + strata()) coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, strata = ovarian$resid.ds) coef(coxph(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian, ties = "breslow")) # 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) # 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) exp(confint(fit_f)) # 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 ) }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 fit_cox <- summary(coxph(Surv(futime, fustat) ~ rx, data = ovarian)) cat("coxph_fast HR :", fit_fast["exp(coef)"], "\n") cat("coxph HR :", fit_cox$coefficients[, "exp(coef)"], "\n") # Stratified by residual disease, compared with coxph(... + strata()) coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, strata = ovarian$resid.ds) coef(coxph(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian, ties = "breslow")) # 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) # 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) exp(confint(fit_f)) # 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 ) }
Computes the calendar time of each analysis look in every simulated trial
from a combination of trigger rules: a target number of events, a planned
calendar time, a maximum calendar time, a minimum time after the previous
look, and a minimum follow-up after a given number of subjects are enrolled.
The cutoffs are computed for all simulated trials in a single C++ pass and
are returned as a matrix that can be passed to the cutoff.looks
argument of analysis_fast and pairwise_fast, or
to the cutoff argument of switch_fast. Separating the
determination of the analysis time from the analysis itself lets the same
cutoffs drive several endpoints (for example overall survival analyzed at
the progression-free survival event-driven looks), several contrasts, or a
modification of the data after an interim analysis.
cutoff_fast( data, event.looks = NULL, time.looks = NULL, max.time = NULL, min.gap = NULL, min.enrolled = NULL, min.followup = NULL, event.subset = NULL, tte.col = "tte", event.col = "event" )cutoff_fast( data, event.looks = NULL, time.looks = NULL, max.time = NULL, min.gap = NULL, min.enrolled = NULL, min.followup = NULL, event.subset = NULL, tte.col = "tte", event.col = "event" )
data |
A data frame of simulated trials, such as the output of
|
event.looks |
Target cumulative event counts. Either a vector with one
positive whole number per look ( |
time.looks |
Planned calendar times, one per look ( |
max.time |
Maximum calendar times, one per look ( |
min.gap |
Minimum time after the previous look, one per look
( |
min.enrolled |
Number of enrolled subjects after which the minimum
follow-up |
min.followup |
Minimum follow-up after the |
event.subset |
An optional logical vector with one element per row of
|
tte.col |
A single character string naming the observed-time column
(numeric and non-negative) used to count events. Defaults to
|
event.col |
A single character string naming the event-indicator column
(numeric or logical, 0 or 1) used to count events. Defaults to
|
For look l of a simulated trial the candidate calendar times are
t_cal, the planned calendar time time.looks[l];
t_event, the calendar time (accrual time plus observed time)
of the d-th counted event, where d is event.looks[l]
(or the entry of the event.looks matrix for that simulation and
look);
t_gap, the cutoff of the previous look plus min.gap[l]
(zero plus min.gap[1] at the first look);
t_enr, the accrual time of the min.enrolled[l]-th
enrolled subject plus min.followup[l];
and the cutoff is
min(max(t_cal, t_event, t_gap, t_enr), max.time[l]), where a
condition that is not supplied (or is NA at that look) is left out.
This is the rule of get_analysis_date() in the simtrial package
(with max.time playing the role of its
max_extension_for_target_event), applied to every simulated trial at
once. For example, event.looks = 300 with max.time = 48 is "300
events or month 48, whichever comes first", and event.looks = 300
with time.looks = 36 is "300 events but not before month 36".
An event or enrollment target that is not met in the simulated data has no
finite calendar time. Unless a max.time caps the look, its cutoff is
then NA, which analysis_fast reports as
reached = FALSE with the statistics of the full data, as it does for
an unreached event.looks target. In this respect the function differs
from simtrial, which uses the last observed event time for an unmet target.
Events are counted on the columns named by tte.col and
event.col, so a look can be triggered by an endpoint other than the
one analyzed (for example the e1_tte and e1_event columns of
an illness-death simulation). By default every subject's events are
counted; event.subset restricts the count to a subset of rows, such as
the control group (data$group == 1), one subgroup, or the two arms of
the primary contrast in a multi-arm trial. When the subjects tie at the
target event time, all events at that time are included in the analysis,
so the analyzed event count can exceed the target, as in
analysis_fast.
A matrix event.looks with one row per simulated trial gives
simulation-specific event targets, for example targets re-estimated at an
interim analysis.
Each look is determined by its own conditions; only min.gap refers
to the previous look. A look with max.time as its only condition is
placed at that time. When the previous look is not reached, t_gap is
infinite, so a look with min.gap is not reached either unless its
max.time caps it. The looks are not reordered: when the conditions
of a later look give an earlier cutoff than the previous look (for example
event.looks = c(150, 220) with max.time = c(NA, 30)), a
warning is given. Supplying min.gap (for example 0) at the later
looks keeps them in order unless their max.time is earlier.
A numeric matrix with one row per simulated trial and one column per
look, holding the calendar cutoffs (NA for a look that is not
reached). The row names are the simulation identifiers, and the attribute
"look.value" holds, for each look, the event target when a common
target is supplied, otherwise the planned calendar time, otherwise
NA. analysis_fast reports this attribute in its
look.value column.
analysis_fast, pairwise_fast,
switch_fast, simdata_fast.
df <- simdata_fast(nsim = 20, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.median = list(12, 18), d.hazard = 0.01, seed = 1) # 150 events or month 30, whichever comes first, then a final look at # 220 events but at least 6 months after the interim cut <- cutoff_fast(df, event.looks = c(150, 220), max.time = c(30, NA), min.gap = c(NA, 6)) head(cut) res <- analysis_fast(df, control = 1, cutoff.looks = cut) head(res) # Looks triggered by the events of the control group only cut_c <- cutoff_fast(df, event.looks = 100, event.subset = df$group == 1) head(cut_c)df <- simdata_fast(nsim = 20, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.median = list(12, 18), d.hazard = 0.01, seed = 1) # 150 events or month 30, whichever comes first, then a final look at # 220 events but at least 6 months after the interim cut <- cutoff_fast(df, event.looks = c(150, 220), max.time = c(30, NA), min.gap = c(NA, 6)) head(cut) res <- analysis_fast(df, control = 1, cutoff.looks = cut) head(res) # Looks triggered by the events of the control group only cut_c <- cutoff_fast(df, event.looks = 100, event.subset = df$group == 1) head(cut_c)
The estimation and testing functions and kmcurve_fast take
vectors of times, event indicators, and group labels, which suits repeated
calls inside a simulation loop. For the analysis of a single data set they
also accept a model formula with a Surv(time, event) response and a
data frame, as the functions of the survival package do:
coxph_fast(Surv(time, status) ~ arm, data = d) survdiff_fast(Surv(time, status) ~ arm + strata(region), data = d) rmst_fast(Surv(time, status) ~ 1, data = d, tau = 12)
The formula is given as the first argument (by position, or by the name of
that argument, such as time; there is no argument named
formula), and the data frame as data or as the second
argument. The other arguments (side,
tau, conf.level, and so on) are given by name as with the
vector interface, and the result is the same as that of the vector
interface with the vectors evaluated from the formula.
The response must be Surv(time, event) with right-censored data,
where event is 1 (or TRUE) for an event and 0 (or
FALSE) for a censored time; an expression such as
Surv(time, status == 2) may be used. A status coded 1 (censored)
and 2 (event), as in data sets of the survival package, gives an
error that suggests this expression. The call to Surv() is
read from the formula and not evaluated, so the survival package
need not be attached. Other forms of Surv(), such as counting-process
data, are not supported.
The right-hand side is 1 for a single group, or a single group
variable with exactly two distinct values for a two-group comparison. Which
of the two is allowed follows the vector interface of each function:
survfit_fast takes only ~ 1; rmst_fast,
wmst_fast, medsurv_fast, and
kmcurve_fast take ~ 1 or a group variable; the other
functions require a group variable. survdiff_fast and
coxph_fast also accept strata() terms, such as
arm + strata(region) or arm + strata(region, sex); several
stratification variables define a stratum for each combination of their
values. Interactions, further covariates, and other formula operators are
not supported, because the functions compare two groups without
adjustment.
With a group variable, control defaults to its first level: the
first level of a factor that occurs in the data, and otherwise the smallest
value. This is the reference level of a model formula with the group
variable as a factor. control can be given to choose the other group
as the control group. A numeric group variable is treated as two groups,
whereas survival::coxph() treats it as a continuous covariate, so
the coefficient of survival::coxph() equals the log hazard ratio of
the two groups only when the two values differ by 1.
The variables are evaluated in data, and variables not found there in
the environment of the formula. Missing values in the variables of the
formula are an error, as in the vector interface; remove the incomplete rows
first, for example with na.omit().
The formula interface adds the cost of evaluating the formula to every call,
so the vector interface, with presorted = TRUE where available,
remains the choice inside simulation loops.
as.data.frame.fastsurv for the results as data
frames.
library(survival) # The same hazard ratio from the formula and from the vectors coxph_fast(Surv(futime, fustat) ~ rx, data = ovarian) coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1) # Stratified log-rank test survdiff_fast(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian) # Single-group and two-group restricted mean survival time rmst_fast(Surv(futime, fustat) ~ 1, data = ovarian, tau = 600) rmst_fast(Surv(futime, fustat) ~ rx, data = ovarian, tau = 600, control = 2)library(survival) # The same hazard ratio from the formula and from the vectors coxph_fast(Surv(futime, fustat) ~ rx, data = ovarian) coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1) # Stratified log-rank test survdiff_fast(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian) # Single-group and two-group restricted mean survival time rmst_fast(Surv(futime, fustat) ~ 1, data = ovarian, tau = 600) rmst_fast(Surv(futime, fustat) ~ rx, data = ovarian, tau = 600, control = 2)
Assembles one or more two-group time-to-event scenarios into a single
scenario_fast object for design-stage exploration. Each scenario is a
complete set of simdata_fast arguments, so the same object can
be drawn with plot.scenario_fast, summarized with
print.scenario_fast, and passed to simdata_fast
to generate data. This is a design-helper that only collects and merges
arguments; it performs no simulation and no compiled computation.
gen_scenario_fast(scenarios, shared = list(), labels = NULL)gen_scenario_fast(scenarios, shared = list(), labels = NULL)
scenarios |
A non-empty list. Each element is a named list of
|
shared |
A named list of |
labels |
Optional character vector of scenario labels, one per scenario. |
Each element of scenarios is a named list of arguments that defines one
scenario. Those arguments override the shared arguments in shared, so
that parameters held constant across scenarios (sample size, accrual, dropout)
are written once in shared and only the varying parameters are written
per scenario. A hazard and a median are alternative specifications of the
same quantity, so a scenario that gives one of them also replaces the other
in shared: e.median in a scenario removes a shared
e.hazard, and likewise for d.*, h01.*, h02.*,
h12.*, and h12.switch.*. The survival specification follows simdata_fast
exactly: e.hazard or e.median given as a two-element list
(control first, treatment second) for the two groups, with e.time
supplying the breakpoints of a piecewise-exponential hazard (last element
Inf). A single hazard or median gives an exponential group.
Two optional fields may appear inside a scenario list and are interpreted by
this function rather than passed on: label sets the scenario label, and
null (logical) flags a null scenario for later filtering. Every other
field is treated as a simdata_fast argument. The scenario label
is taken from labels, then from the label field, then from the
names of scenarios, then from a default "Scenario k".
A factorial set of scenarios is built by constructing the scenarios
list with the base tools, for example mapping over the rows of an
expand.grid of the parameters that vary.
An object of class scenario_fast: a list with element
scenarios (a named list, one entry per scenario, each holding the
merged argument list args, the label, and the null
flag) and element shared (the shared arguments).
simdata_fast, plot.scenario_fast,
print.scenario_fast
# Three two-group scenarios sharing accrual and sample size: a proportional- # hazards case, a delayed-effect case, and a crossing-hazards case. scn <- gen_scenario_fast( scenarios = list( "Proportional" = list(e.median = list(12, 18)), "Delayed effect" = list( e.hazard = list(log(2) / 12, c(log(2) / 12, log(2) / 22)), e.time = c(0, 6, Inf) ), "Crossing" = list( e.hazard = list(log(2) / 12, c(log(2) / 7, log(2) / 24)), e.time = c(0, 5, Inf) ) ), shared = list(n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12) ) print(scn) plot(scn) # Generate data for every scenario with a short composition step. From here # analysis_fast() and simsummary_fast() are applied per scenario as usual. sim_list <- lapply(scn$scenarios, function(s) { do.call(simdata_fast, c(s$args, list(nsim = 100, seed = 1))) }) nrow(sim_list[[1]])# Three two-group scenarios sharing accrual and sample size: a proportional- # hazards case, a delayed-effect case, and a crossing-hazards case. scn <- gen_scenario_fast( scenarios = list( "Proportional" = list(e.median = list(12, 18)), "Delayed effect" = list( e.hazard = list(log(2) / 12, c(log(2) / 12, log(2) / 22)), e.time = c(0, 6, Inf) ), "Crossing" = list( e.hazard = list(log(2) / 12, c(log(2) / 7, log(2) / 24)), e.time = c(0, 5, Inf) ) ), shared = list(n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12) ) print(scn) plot(scn) # Generate data for every scenario with a short composition step. From here # analysis_fast() and simsummary_fast() are applied per scenario as usual. sim_list <- lapply(scn$scenarios, function(s) { do.call(simdata_fast, c(s$args, list(nsim = 100, seed = 1))) }) nrow(sim_list[[1]])
Builds the Kaplan-Meier step functions for a single realized data set (for
example one replicate from simdata_fast) and returns a
kmcurve_fast object that can be drawn with
plot.kmcurve_fast and summarized with
print.kmcurve_fast. Both the two-group case and the single-group
case are supported. This is an analysis-stage and plotting helper, written in
plain R for use on a single data set, and is not intended to be called inside
a simulation loop. The fast point estimator survfit_fast remains
the tool for repeated evaluation inside loops.
kmcurve_fast(time, event, group = NULL, control = NULL, data = NULL)kmcurve_fast(time, event, group = NULL, control = NULL, data = NULL)
time |
A numeric vector of follow-up times for all subjects.
Alternatively, a formula |
event |
An integer or numeric vector of event indicators (1 = event,
0 = censored), aligned with |
group |
A vector identifying the groups, aligned with |
control |
The value of |
data |
A data frame in which the variables of a formula given as
|
For each group the object stores the distinct event times, the number of
events and the number at risk at those times, the Kaplan-Meier survival
estimate, the Greenwood standard error, the Nelson-Aalen increments
d_i / Y_i (used by the plot method to estimate a smoothed time-varying
hazard ratio), and the censoring times. The pooled data are also retained so
that the plot and print methods can reuse coxph_fast,
rmst_fast, and ahr_fast for the constant Cox
hazard ratio, the restricted mean survival time, and the average hazard
ratio. When group is NULL or has a single distinct value, a
single-group object is built and the plot method draws one survival curve
with no hazard ratio.
An object of class kmcurve_fast: a list with a logical
two_group flag, the group values and their labels, the per-group
Kaplan-Meier summaries in element km, the pooled data, and the
per-group sample sizes n and event counts events. The
curves are returned as a data frame by
as.data.frame.kmcurve_fast.
plot.kmcurve_fast, print.kmcurve_fast,
as.data.frame.kmcurve_fast, survfit_fast,
coxph_fast,
rmst_fast, ahr_fast
set.seed(1) n <- 150 t0 <- rexp(n, log(2) / 12) t1 <- rexp(n, log(2) / 18) cens <- runif(2 * n, 0, 30) time <- pmin(c(t0, t1), cens) event <- as.integer(c(t0, t1) <= cens) group <- rep(c(0, 1), each = n) fit <- kmcurve_fast(time, event, group, control = 0) print(fit) plot(fit) plot(fit, hr = TRUE, rmst = TRUE, tau = 18) plot(fit, risk_table = TRUE, conf.style = "band") head(as.data.frame(fit)) # Formula interface d <- data.frame(time, event, group) kmcurve_fast(Surv(time, event) ~ group, data = d) # Single-group object: pass group = NULL (or a single-level group) fit1 <- kmcurve_fast(time[group == 0], event[group == 0]) plot(fit1, rmst = TRUE, tau = 18)set.seed(1) n <- 150 t0 <- rexp(n, log(2) / 12) t1 <- rexp(n, log(2) / 18) cens <- runif(2 * n, 0, 30) time <- pmin(c(t0, t1), cens) event <- as.integer(c(t0, t1) <= cens) group <- rep(c(0, 1), each = n) fit <- kmcurve_fast(time, event, group, control = 0) print(fit) plot(fit) plot(fit, hr = TRUE, rmst = TRUE, tau = 18) plot(fit, risk_table = TRUE, conf.style = "band") head(as.data.frame(fit)) # Formula interface d <- data.frame(time, event, group) kmcurve_fast(Surv(time, event) ~ group, data = d) # Single-group object: pass group = NULL (or a single-level group) fit1 <- kmcurve_fast(time[group == 0], event[group == 0]) plot(fit1, rmst = TRUE, tau = 18)
Computes the max-combo test, the maximum over a set of Fleming-Harrington weighted log-rank statistics, for comparing survival between two groups under non-proportional hazards. The C++ backend evaluates every weighted numerator and the full between-scheme covariance matrix in a single scan over the pooled sorted data, and the p-value is obtained from the multivariate normal distribution implied by the correlation of the component statistics. The test is robust to the shape of the hazard difference because the most extreme of several complementary weights is taken, with the multiplicity accounted for through the joint distribution.
maxcombo_fast( time, event, group, control, side = 2, rho = c(0, 0, 1, 1), gamma = c(0, 1, 0, 1), presorted = FALSE, abseps = 1e-05, maxpts = 25000, data = NULL )maxcombo_fast( time, event, group, control, side = 2, rho = c(0, 0, 1, 1), gamma = c(0, 1, 0, 1), presorted = FALSE, abseps = 1e-05, maxpts = 25000, data = NULL )
time |
A numeric vector of follow-up times for all subjects.
Alternatively, a formula |
event |
An integer or numeric vector of event indicators
(1 = event, 0 = censored), aligned with |
group |
A vector of group labels aligned with |
control |
A scalar value indicating which level of |
side |
An integer, either 1 or 2. If |
rho |
A numeric vector of Fleming-Harrington first parameters, one per
component weight. Defaults to |
gamma |
A numeric vector of Fleming-Harrington second parameters, one
per component weight, aligned with |
presorted |
A logical value. If |
abseps |
A single positive numeric value, the absolute error tolerance
passed to the multivariate normal integration. Defaults to 1e-5. The
quasi-Monte-Carlo integration (four or more weights, or a two-sided test)
stops at |
maxpts |
A single positive integer, the maximum number of function evaluations for the quasi-Monte-Carlo integration used with four or more weights or with a two-sided test (see Details). Defaults to 25000. Larger values reduce the integration error at a higher computing cost. |
data |
A data frame in which the variables of a formula given as
|
Each component is a Fleming-Harrington G(rho, gamma) weighted log-rank
statistic Z_k = U_k / sqrt(V_kk), where U_k = sum w_k (O_1 - E_1) and the
weight is w_k = S(t-)^rho_k (1 - S(t-))^gamma_k evaluated from the
left-continuous pooled Kaplan-Meier estimate S(t-). The between-scheme
covariance is V_ab = sum w_a w_b v, with v the hypergeometric variance
increment shared by all schemes, so the diagonal of V reproduces the
single-scheme weighted variances and the off-diagonal entries give the
correlation matrix R of the component Z-scores. The sign convention matches
survdiff_fast: a component Z is negative when the treatment
group is favored.
The max-combo statistic and its p-value depend on side. When
side = 1, the statistic is the most negative component, min_k Z_k, so
that a negative value favors the treatment group in the same way as
survdiff_fast with side = 1. The one-sided p-value is
1 - P(G_1 >= m, ..., G_K >= m) for G distributed as multivariate normal with
mean zero and correlation R, where m = min_k Z_k. When side = 2, the
statistic is max_k abs(Z_k) and the p-value is
1 - P(-m <= G_1 <= m, ..., -m <= G_K <= m).
The joint normal probability is evaluated by dimension. With a single weight
the univariate normal is used. With two or three weights and a one-sided
test, where the integration region is an orthant (all upper limits
infinite), the deterministic
TVPACK algorithm is used. For the two-sided test, whose
region is a bounded rectangle, and for four or more weights, the randomized
quasi-Monte-Carlo GenzBretz algorithm is used. It
uses R's random-number generator, so call set.seed() beforehand for
reproducible p-values, and it stops when its error estimate falls below
abseps or after maxpts function evaluations, whichever comes
first.
The weights of G(0,1) and G(1,0) add up to the weight of G(0,0) at every
time, because (1 - S) + S = 1, so the numerator of the G(0,0) component is
the sum of the other two. The correlation matrix of the default weights, or
of any set that contains these three, is therefore singular, and the
quasi-Monte-Carlo integral converges slowly. With the default maxpts
the absolute error of the p-value is typically of the order of 1e-4 rather
than abseps. In a check with the default weights and a one-sided
test, the p-values of 20,000 simulated analyses (two looks of 10,000 trials
with a delayed treatment effect) were computed twice with different seeds;
the two values differed by less than 1e-3 in every analysis and by at most
about 3e-4 in 99 percent of them. A larger maxpts reduces the error
at a higher computing cost. In a simulation study, the integration error affects only
the decisions for p-values within about that error of the nominal level.
When presorted = TRUE, the inputs are assumed to be sorted in ascending
order of time and the internal order() call is skipped, which is
useful inside simulation loops where the data are generated in sorted order.
An object of class "maxcombo_fast", a named numeric vector of
length two with elements statistic (the max-combo statistic;
min_k Z_k when side = 1, so a negative value favors treatment, and
max_k abs(Z_k) when side = 2) and p.value. The component
Z-scores are stored in the attribute z, their correlation matrix in
corr, the Fleming-Harrington parameters in rho and
gamma, the requested side, and the total sample size in
n. Returns NA values (still with class
"maxcombo_fast") when any component variance is zero or not finite.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Karrison, T. G. (2016). Versatile tests for comparing survival curves based on weighted log-rank statistics. The Stata Journal, 16(3), 678-690.
Lin, R. S., Lin, J., Roychoudhury, S., et al. (2020). Alternative analysis methods for time to event endpoints under nonproportional hazards: a comparative analysis. Statistics in Biopharmaceutical Research, 12(2), 187-198.
survdiff_fast for the single-scheme weighted log-rank test.
library(survival) # Standard four-weight max-combo, one-sided fit <- maxcombo_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1) fit["statistic"] fit["p.value"] # Two-sided test maxcombo_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 1, side = 2) # Custom weight set: proportional plus late-difference only maxcombo_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 1, rho = c(0, 0), gamma = c(0, 1)) # presorted = TRUE: sort once outside, reuse inside a loop ord <- order(ovarian$futime) maxcombo_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord], control = 1, presorted = TRUE) # Formula interface maxcombo_fast(Surv(futime, fustat) ~ rx, data = ovarian, side = 1) # Cross-check against simtrial::maxcombo (one-sided) if (requireNamespace("simtrial", quietly = TRUE)) { df <- data.frame( stratum = "All", treatment = ifelse(ovarian$rx == 2, "experimental", "control"), tte = ovarian$futime, event = ovarian$fustat ) simtrial::maxcombo(df, rho = c(0, 0, 1, 1), gamma = c(0, 1, 0, 1), return_corr = TRUE)$p_value }library(survival) # Standard four-weight max-combo, one-sided fit <- maxcombo_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1) fit["statistic"] fit["p.value"] # Two-sided test maxcombo_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 1, side = 2) # Custom weight set: proportional plus late-difference only maxcombo_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 1, rho = c(0, 0), gamma = c(0, 1)) # presorted = TRUE: sort once outside, reuse inside a loop ord <- order(ovarian$futime) maxcombo_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord], control = 1, presorted = TRUE) # Formula interface maxcombo_fast(Surv(futime, fustat) ~ rx, data = ovarian, side = 1) # Cross-check against simtrial::maxcombo (one-sided) if (requireNamespace("simtrial", quietly = TRUE)) { df <- data.frame( stratum = "All", treatment = ifelse(ovarian$rx == 2, "experimental", "control"), tte = ovarian$futime, event = ovarian$fustat ) simtrial::maxcombo(df, rho = c(0, 0, 1, 1), gamma = c(0, 1, 0, 1), return_corr = TRUE)$p_value }
Estimates the Kaplan-Meier median survival time for a single group, or the difference in median survival time between a treatment group and a control group, together with standard errors, confidence intervals and a Wald test for the difference. The estimator is designed for repeated evaluation inside simulation loops, with the single scan over the sorted data performed in C++.
medsurv_fast( time, event, group = NULL, control = NULL, side = 2, conf.level = 0.95, conf.type = "log", method = c("km", "nph"), bw = NULL, presorted = FALSE, data = NULL )medsurv_fast( time, event, group = NULL, control = NULL, side = 2, conf.level = 0.95, conf.type = "log", method = c("km", "nph"), bw = NULL, presorted = FALSE, data = NULL )
time |
Numeric vector of event or censoring times.
Alternatively, a formula |
event |
Integer vector, 1 for an event and 0 for censoring. |
group |
Optional grouping vector with exactly two distinct levels for a two-group comparison. If omitted, a single-group median is returned. |
control |
The level of |
side |
Either 2 for a two-sided test or 1 for a one-sided test of treatment superiority (difference greater than 0). |
conf.level |
Confidence level for the intervals. |
conf.type |
Confidence interval type for each group median, either
|
method |
Variance method, either |
bw |
Optional kernel bandwidth for the hazard at the median, used only
when |
presorted |
Logical; set to |
data |
A data frame in which the variables of a formula given as
|
The median in each group is the first event time at which the Kaplan-Meier
estimate drops to 0.5 or below. When the estimate equals 0.5 exactly over an
interval, the midpoint of that interval is returned (the next event time, or
the last observed time when no later event exists, closes the interval). This
is the convention used by survfit, including its tolerance of
sqrt(.Machine$double.eps) when comparing the estimate with 0.5. The
point estimate is the same for both variance methods. When the estimate
stays at 0.5 up to an infinite observed time (for example a subject without
a finite event or dropout time in simulated data), the midpoint is not
defined, and the median of that group and the statistics that depend on it
are NA.
Two variance methods are available through the method argument. With
method = "km" the variance of the estimated median follows the
counting process delta method, var(median) = greenwood_sum / hazard^2, where
the instantaneous hazard at the median is obtained by a Ramlau-Hansen kernel
smoother with an Epanechnikov kernel and bandwidth bw. With
method = "nph" the variance reproduces the computation used by
nph::nphparams with param_type = "Q" and
haz_method = "local": a local constant hazard at the median and a
sum 1 / (Y - k)^2 variance increment. The two methods give the same
median but generally different standard errors, since the variance of a
quantile depends on the local hazard estimate. When the Kaplan-Meier and
Nelson-Aalen medians coincide, which is the usual case, method = "nph"
reproduces the nph::nphparams standard error and p-value to numerical
precision.
The difference is computed as treatment minus control, so a positive difference indicates a longer median survival time under treatment.
A named numeric vector of class "medsurv_fast". For a single
group the elements are the median, its standard error and confidence
limits. For two groups the elements are the control and treatment medians,
their difference, the standard errors, the confidence limits and the Wald
statistics for the difference.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Fleming, T. R., & Harrington, D. P. (1991). Counting Processes and Survival Analysis. New York: John Wiley & Sons.
set.seed(1) n <- 200 g <- rep(0:1, each = n / 2) tt <- rexp(n, rate = ifelse(g == 0, 0.1, 0.07)) cc <- rexp(n, rate = 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) medsurv_fast(time, event, group = g, control = 0) medsurv_fast(time, event, group = g, control = 0, method = "nph") # Formula interface d <- data.frame(time, event, g) medsurv_fast(Surv(time, event) ~ g, data = d)set.seed(1) n <- 200 g <- rep(0:1, each = n / 2) tt <- rexp(n, rate = ifelse(g == 0, 0.1, 0.07)) cc <- rexp(n, rate = 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) medsurv_fast(time, event, group = g, control = 0) medsurv_fast(time, event, group = g, control = 0, method = "nph") # Formula interface d <- data.frame(time, event, g) medsurv_fast(Surv(time, event) ~ g, data = d)
Compares the Kaplan-Meier survival probabilities of two groups at a
prespecified milestone timepoint. The point estimate of interest is the
difference in milestone survival, treatment minus control. Three inference
methods are provided. The "wald" method uses the unpooled Greenwood
variance directly. The "loglog" and "mover" methods build the
confidence interval for the difference with the method of variance estimates
recovery (MOVER), recovering the variance from the one-sample complementary
log-log and log transformed confidence intervals respectively. When the
estimate of a group is 0 or 1 at tau, its one-sample interval
degenerates to the estimate, so the interval of the difference is still
given; the "loglog" test statistic is then NA. See Tang
(2021) for the MOVER difference interval and Tang (2022) for the use of
milestone survival in trial design.
milestone_fast( time, event, group, control, side = 2, conf.level = 0.95, tau, method = c("wald", "loglog", "mover"), presorted = FALSE, data = NULL )milestone_fast( time, event, group, control, side = 2, conf.level = 0.95, tau, method = c("wald", "loglog", "mover"), presorted = FALSE, data = NULL )
time |
A numeric vector of follow-up times.
Alternatively, a formula |
event |
An integer vector of event indicators, 1 for an event and 0 for a censored observation. |
group |
A vector with exactly two distinct values identifying the group. |
control |
The value of |
side |
1 for a one-sided test in the direction of treatment benefit
(treatment milestone survival larger than control) or 2 for a two-sided
test (default 2). The confidence interval is always reported as a
two-sided interval at |
conf.level |
The confidence level for the reported intervals. |
tau |
The milestone timepoint at which the survival probabilities are compared. A single positive number. A warning is given when it exceeds the largest observed time of either group, because the Kaplan-Meier curve is not estimated beyond that time. |
method |
The inference method for the difference in milestone survival,
one of |
presorted |
Logical. If |
data |
A data frame in which the variables of a formula given as
|
The "wald" statistic is the difference divided by its Greenwood
standard error. The "loglog" statistic is the difference of the
complementary log-log transforms of the two estimates divided by its
standard error, so it is negative when the treatment group has the higher
milestone survival, and its one-sided p-value is the lower tail. The
"mover" statistic is the upper (1 - conf.level) / 2 quantile
of the standard normal distribution (1.96 for conf.level = 0.95)
times the difference divided by the half-width of the MOVER interval on the
side of zero. It is obtained by inverting the interval at
conf.level, so the statistic and its p-value depend on
conf.level, and the two-sided p-value is below
1 - conf.level exactly when the interval excludes zero. For
"wald" and "mover" a positive statistic favors treatment and
the one-sided p-value is the upper tail.
An object of class "milestone_fast", a list with the
per-group milestone survival estimates and standard errors, the difference
estimate with its confidence interval, the test statistic, and the p-value.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Tang, Y. (2021). Some new confidence intervals for Kaplan-Meier based estimators from one and two sample survival data. Statistics in Medicine, 40(23), 4961-4976.
Tang, Y. (2022). Complex survival trial design by the product integration method. Statistics in Medicine, 41(4), 798-814.
set.seed(1) time <- c(rexp(50, 0.1), rexp(50, 0.07)) event <- rep(1, 100) group <- rep(c(0, 1), each = 50) milestone_fast(time, event, group, control = 0, tau = 10, method = "loglog") # Formula interface d <- data.frame(time, event, group) milestone_fast(Surv(time, event) ~ group, data = d, tau = 10)set.seed(1) time <- c(rexp(50, 0.1), rexp(50, 0.07)) event <- rep(1, 100) group <- rep(c(0, 1), each = 50) milestone_fast(time, event, group, control = 0, tau = 10, method = "loglog") # Formula interface d <- data.frame(time, event, group) milestone_fast(Surv(time, event) ~ group, data = d, tau = 10)
Runs the two-group analysis of analysis_fast once for each
experimental arm against a common control arm, on multi-arm simulated data
such as the output of simdata_fast with a length-greater-than-two
n. The contrasts are analyzed at a shared set of looks, and the results
are stacked into one long data frame with an arm column, optionally with
a Bonferroni-adjusted p-value across the contrasts.
pairwise_fast( data, control, event.looks = NULL, time.looks = NULL, primary = NULL, arms = NULL, stat = "logrank", adjust = c("none", "bonferroni"), p.col = NULL, cutoff.looks = NULL, ... )pairwise_fast( data, control, event.looks = NULL, time.looks = NULL, primary = NULL, arms = NULL, stat = "logrank", adjust = c("none", "bonferroni"), p.col = NULL, cutoff.looks = NULL, ... )
data |
A data frame of simulated trial data with columns |
control |
The group label of the control arm. |
event.looks |
A numeric vector of target cumulative event counts for the
primary contrast, one per look. Mutually exclusive with |
time.looks |
A numeric vector of calendar times, one per look. Mutually
exclusive with |
primary |
The group label of the experimental arm whose control-versus-arm
comparison defines the shared calendar cutoff. Required with
|
arms |
A vector of experimental arm labels to compare against the control.
Defaults to every group other than |
stat |
The test statistic passed to |
adjust |
Multiplicity adjustment across contrasts. Either |
p.col |
The name of the p-value column to adjust. By default it is the
column matching the chosen |
cutoff.looks |
A numeric matrix of per-simulation calendar cutoffs with
one row per simulated trial and one column per look, as returned by
|
... |
Further arguments passed to |
Three timing regimes are supported, exactly one of which must be requested.
With time.looks, every contrast is analyzed at the same fixed calendar
time or times, so all contrasts share the same data cutoff by construction.
With event.looks, the analysis is event-driven and primary names
the experimental arm whose control-versus-arm comparison defines the data
cutoff. The primary contrast is analyzed at the requested cumulative event
counts, its per-simulation calendar cutoffs are recorded, and then every
contrast, including the primary one, is analyzed at those same cutoffs. This
reproduces the standard design in which the primary event-driven analysis
fixes a single data cutoff at which all comparisons are performed. For a
simulation in which the primary event target is not reached at a look, that
look is marked with reached = FALSE, and its counts and statistics are
NA for every contrast. The shared cutoffs are computed by cutoff_fast on the
events of the control and primary arms.
With cutoff.looks, every contrast is analyzed at the supplied
per-simulation calendar cutoffs, for example those returned by
cutoff_fast for a combination of event and calendar-time rules.
A look that is not reached in a simulation (NA cutoff) is marked with
reached = FALSE and NA counts and statistics, as in the
event-driven regime. This differs from analysis_fast, which
analyzes the full data at an unreached look; in
simsummary_fast an unreached look of pairwise_fast
output therefore never crosses a boundary.
The Bonferroni option multiplies each p-value by the number of contrasts and
caps it at one, controlling the family-wise error rate across the
control-versus-arm comparisons at each look. Multiplicity across looks is a
separate matter handled by group-sequential boundaries in
simsummary_fast, not by this adjustment; given the output of
this function, simsummary_fast summarizes each arm separately (the
Bonferroni-adjusted p-values can be used there through
p.col = "p.adj").
This is a single-endpoint helper: it reads the tte and event
columns and does not support subgroups. Comparisons for a second endpoint are
obtained by calling pairwise_fast again on that endpoint's columns.
A data frame with one row per contrast, look, and simulation. The
first column is arm (the experimental arm compared against
control), followed by the columns returned by
analysis_fast: sim, look, look.value,
cutoff, reached, the enrollment and event counts, and the
statistic columns for the chosen stat. When adjust =
"bonferroni" a p.adj column is appended.
simdata_fast, analysis_fast,
cutoff_fast
# Three-arm trial: control (group 1) and two experimental arms. dfk <- simdata_fast( nsim = 100, n = c(120, 120, 120), a.time = c(0, 12), a.rate = 360 / 12, e.median = list(12, 16, 20), seed = 8 ) # Fixed calendar look at month 30, Bonferroni across the two contrasts. pw <- pairwise_fast(dfk, control = 1, time.looks = 30, stat = "logrank", side = 1, adjust = "bonferroni") head(pw) # Event-driven: arm 3 is the primary contrast and its 200th control-plus-arm-3 # event fixes the cutoff at which both contrasts are analyzed. pw2 <- pairwise_fast(dfk, control = 1, event.looks = 200, primary = 3, stat = "logrank", side = 1, adjust = "bonferroni") head(pw2)# Three-arm trial: control (group 1) and two experimental arms. dfk <- simdata_fast( nsim = 100, n = c(120, 120, 120), a.time = c(0, 12), a.rate = 360 / 12, e.median = list(12, 16, 20), seed = 8 ) # Fixed calendar look at month 30, Bonferroni across the two contrasts. pw <- pairwise_fast(dfk, control = 1, time.looks = 30, stat = "logrank", side = 1, adjust = "bonferroni") head(pw) # Event-driven: arm 3 is the primary contrast and its 200th control-plus-arm-3 # event fixes the cutoff at which both contrasts are analyzed. pw2 <- pairwise_fast(dfk, control = 1, event.looks = 200, primary = 3, stat = "logrank", side = 1, adjust = "bonferroni") head(pw2)
Draws the Kaplan-Meier curves of a kmcurve_fast object. With
rmst = TRUE the area under each curve up to tau is shaded to
visualize the restricted mean survival time. With hr = TRUE a second
panel shows the smoothed time-varying hazard ratio together with the constant
Cox proportional-hazards estimate and, optionally, the average hazard ratio,
in the style of stylized nonproportional-hazards figures. For a single-group
object only one survival curve is drawn; hr = TRUE is ignored because a
hazard ratio requires two groups. With risk_table = TRUE the numbers at
risk are printed below the survival panel.
## S3 method for class 'kmcurve_fast' plot( x, hr = FALSE, rmst = FALSE, tau = NULL, bw = NULL, n_grid = 200L, conf.int = TRUE, hr_max = NULL, cox_line = TRUE, ahr_line = FALSE, mark_time = TRUE, conf.level = 0.95, conf.type = "log", col_control = "#1b9e77", col_treat = "#d95f02", col_hr = "#3060c0", col_cox = "black", col_ahr = "grey40", lty_control = 1, lty_treat = 1, lty_hr = 1, lwd_control = 2, lwd_treat = 2, lwd_hr = 2, base_size = 1, legend_pos = "topright", xlab = "Time", ylab = "Survival probability", main_km = "Kaplan-Meier", main_hr = "Hazard ratio", group_labels = NULL, risk_table = FALSE, risk_times = NULL, conf.style = c("lines", "band"), ... )## S3 method for class 'kmcurve_fast' plot( x, hr = FALSE, rmst = FALSE, tau = NULL, bw = NULL, n_grid = 200L, conf.int = TRUE, hr_max = NULL, cox_line = TRUE, ahr_line = FALSE, mark_time = TRUE, conf.level = 0.95, conf.type = "log", col_control = "#1b9e77", col_treat = "#d95f02", col_hr = "#3060c0", col_cox = "black", col_ahr = "grey40", lty_control = 1, lty_treat = 1, lty_hr = 1, lwd_control = 2, lwd_treat = 2, lwd_hr = 2, base_size = 1, legend_pos = "topright", xlab = "Time", ylab = "Survival probability", main_km = "Kaplan-Meier", main_hr = "Hazard ratio", group_labels = NULL, risk_table = FALSE, risk_times = NULL, conf.style = c("lines", "band"), ... )
x |
A |
hr |
Logical. If |
rmst |
Logical. If |
tau |
Upper limit used for the restricted mean survival time and the average hazard ratio. Defaults to the smaller of the two groups' largest observed times. |
bw |
Kernel bandwidth for the smoothed hazard ratio. Defaults to a value derived from the spread of the event times. |
n_grid |
Number of grid points for the smoothed hazard ratio. |
conf.int |
Logical. If |
hr_max |
Upper limit of the hazard-ratio axis. Defaults to a value derived from the estimated hazard ratio. |
cox_line |
Logical. If |
ahr_line |
Logical. If |
mark_time |
Logical. If |
conf.level |
Confidence level for the survival confidence limits. |
conf.type |
Type of survival confidence limits: one of |
col_control, col_treat, col_hr, col_cox, col_ahr
|
Colors for the control curve, the treatment curve, the smoothed hazard ratio, the Cox line, and the average-hazard-ratio line. |
lty_control, lty_treat
|
Line types for the control and treatment curves. |
lty_hr |
Line type for the smoothed hazard ratio. |
lwd_control, lwd_treat
|
Line widths for the control and treatment curves. |
lwd_hr |
Line width for the lines in the hazard-ratio panel. |
base_size |
Relative size of all text. |
legend_pos |
Keyword giving the legend position in both panels, passed to
|
xlab, ylab
|
Axis labels for the survival panel. |
main_km, main_hr
|
Panel titles. |
group_labels |
Optional length-two character vector of legend labels for
the control and treatment groups. They also label the rows of the
number-at-risk table, which are otherwise labeled |
risk_table |
Logical. If |
risk_times |
Times at which the numbers at risk are given. Defaults to the tick marks of the time axis. |
conf.style |
How the pointwise confidence limits are drawn when
|
... |
Further arguments, currently ignored. |
The smoothed hazard ratio is obtained by kernel smoothing the Nelson-Aalen
increments of each group. For group j the kernel hazard estimate is
the sum over event times of (1 / bw) K((t - t_i) / bw) (d_i / Y_i),
with an Epanechnikov kernel K, and the hazard ratio is the treatment
estimate divided by the control estimate. The bandwidth bw controls the
width of the averaging window and materially affects the result, so it is
exposed as an argument. To limit boundary bias the hazard ratio is drawn over
an interior time window, and points where the control hazard estimate is
effectively zero are dropped.
The constant Cox hazard ratio is taken from coxph_fast, the
average hazard ratio from ahr_fast, and the restricted mean
survival time difference from rmst_fast.
The number at risk at a time is the number of subjects whose observed time is
at or after it, as reported by summary() of a
survival::survfit() object. The table has one row per group, in the
colors of the curves, at the tick marks of the time axis or at
risk_times; the bottom margin (and, for long group labels, the left
margin) is enlarged to hold it. The curves themselves are available as a data
frame from as.data.frame.kmcurve_fast.
The input x, invisibly.
kmcurve_fast, print.kmcurve_fast,
as.data.frame.kmcurve_fast
Draws, for each scenario in a scenario_fast object, the analytic
control and experimental survival curves together with the piecewise-constant
hazard ratio on a secondary axis, in the style of the stylized
nonproportional-hazards figures used in the design literature. The curves are
computed analytically from the piecewise-exponential parameters, so no
simulation is involved.
## S3 method for class 'scenario_fast' plot( x, which = NULL, tmax = NULL, hr_max = NULL, n_time = 481L, mfrow = NULL, col_control = "#1b9e77", col_treat = "#d95f02", col_hr = "#3060c0", lty_control = 1, lty_treat = 1, lty_hr = 2, lwd_control = 2, lwd_treat = 2, lwd_hr = 1.5, base_size = 1, xlab = "Time", ylab = "Survival probability", legend_panel = 1L, legend_pos = "topright", ... )## S3 method for class 'scenario_fast' plot( x, which = NULL, tmax = NULL, hr_max = NULL, n_time = 481L, mfrow = NULL, col_control = "#1b9e77", col_treat = "#d95f02", col_hr = "#3060c0", lty_control = 1, lty_treat = 1, lty_hr = 2, lwd_control = 2, lwd_treat = 2, lwd_hr = 1.5, base_size = 1, xlab = "Time", ylab = "Survival probability", legend_panel = 1L, legend_pos = "topright", ... )
x |
A |
which |
Optional integer or character vector selecting a subset of scenarios to draw. |
tmax |
Upper limit of the time axis. Defaults to a value derived from the scenarios. |
hr_max |
Upper limit of the hazard-ratio axis. Defaults to a value derived from the scenarios. |
n_time |
Number of points on the time grid. |
mfrow |
Optional length-two panel layout |
col_control, col_treat, col_hr
|
Colors for the control survival curve, the experimental survival curve, and the hazard ratio. |
lty_control, lty_treat, lty_hr
|
Line types for the control survival curve, the experimental survival curve, and the hazard ratio. |
lwd_control, lwd_treat, lwd_hr
|
Line widths for the control survival curve, the experimental survival curve, and the hazard ratio. |
base_size |
Relative size of all text (titles, axis labels, axis
annotation, and the legend). The default |
xlab, ylab
|
Axis labels for the survival panels. |
legend_panel |
Index (within the drawn subset) of the panel that carries the legend. |
legend_pos |
Keyword giving the legend position, passed to
|
... |
Further arguments, currently ignored. |
Each scenario is shown in its own panel arranged on a grid. The hazard ratio is rescaled into the survival plotting region and read on the right-hand axis; a reference line marks a hazard ratio of one. A single-arm scenario (no treatment specification) is drawn as a single survival curve with no hazard ratio.
The input x, invisibly.
gen_scenario_fast, print.scenario_fast
Draws the stopping probabilities of a simsummary_fast result
by look. With what = "cumulative" (default) one panel shows, for
every block (arm or population), the cumulative probability of stopping
for efficacy, which at the final look is the rejection rate, and, when
futility boundaries were applied, the cumulative probability of stopping
for futility. With what = "stage" each block has a panel of stacked
bars of the probabilities of stopping for efficacy and for futility at each
look.
## S3 method for class 'simsummary_fast' plot( x, what = c("cumulative", "stage"), xaxis = c("look", "events", "time"), col = NULL, ylim = NULL, legend_pos = "topleft", xlab = NULL, ylab = NULL, main = NULL, ... )## S3 method for class 'simsummary_fast' plot( x, what = c("cumulative", "stage"), xaxis = c("look", "events", "time"), col = NULL, ylim = NULL, legend_pos = "topleft", xlab = NULL, ylab = NULL, main = NULL, ... )
x |
An object of class |
what |
|
xaxis |
The horizontal axis: |
col |
Colors of the blocks, recycled. Defaults to a qualitative palette. |
ylim |
Limits of the vertical axis. Defaults to the range from 0 to slightly above the largest probability drawn. |
legend_pos |
Keyword giving the legend position, passed to
|
xlab, ylab, main
|
Axis labels and title. Defaults depend on |
... |
Further arguments, currently ignored. |
The plot uses the look rows of each block, which must keep all the looks of
the design in their original order and equal the rows of the stored result,
as in the print method (see print.simsummary_fast): a
selection of columns or of some of the looks, rows in a different order, or
rows of another result combined by rbind() give an error. The overall
rows are not used, so they are not compared. A result created by
FastSurvival 1.1.0, which does not store the looks and rows, gives an error
that asks for simsummary_fast() to be run again. The cumulative
probability of stopping for futility is the cumulative sum of
prob.stop.futility over the looks.
The looks are placed on the horizontal axis by their index
(xaxis = "look"), by the expected number of events at the look
("events", the column n.event.mean), or by the expected
calendar time of the look ("time", the column cutoff.mean).
The input x, invisibly.
simsummary_fast, print.simsummary_fast.
df <- simdata_fast( nsim = 200, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(0.05, 0.035), seed = 1 ) res <- analysis_fast(df, control = 1, event.looks = c(60, 105, 150), stat = c("logrank", "coxph"), side = 1) oc <- simsummary_fast(res, eff.col = "logrank.z", efficacy = c(NA, -2.96, -1.97), fut.col = "cox.coef", futility = c(log(1.2), NA, NA), direction = "lower") plot(oc) plot(oc, what = "stage", xaxis = "events")df <- simdata_fast( nsim = 200, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(0.05, 0.035), seed = 1 ) res <- analysis_fast(df, control = 1, event.looks = c(60, 105, 150), stat = c("logrank", "coxph"), side = 1) oc <- simsummary_fast(res, eff.col = "logrank.z", efficacy = c(NA, -2.96, -1.97), fut.col = "cox.coef", futility = c(log(1.2), NA, NA), direction = "lower") plot(oc) plot(oc, what = "stage", xaxis = "events")
The p-value column is labeled Pr(>|z|) for a two-sided test and
Pr(<z) for a one-sided test, whose p-value is the lower tail in the
direction of treatment benefit.
## S3 method for class 'ahr_fast' print(x, digits = max(3L, getOption("digits") - 3L), ...)## S3 method for class 'ahr_fast' print(x, digits = max(3L, getOption("digits") - 3L), ...)
x |
an object of class |
digits |
number of significant digits to print |
... |
further arguments (currently ignored) |
x, invisibly
Formats and prints an ahsw_fast object. The header shows the
truncation time and the control label. The body shows the per-group average
hazard with survival weight, followed by the between-group contrasts: the
ratio of average hazards (treatment over control) and the difference of
average hazards (treatment minus control), each with a confidence interval
and a p-value that follows the side recorded in the object.
## S3 method for class 'ahsw_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'ahsw_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
library(survival) fit <- ahsw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, tau = 600) print(fit)library(survival) fit <- ahsw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, tau = 600) print(fit)
Formats and prints a coxph_fast object similarly to
summary(survival::coxph(...)), showing the point estimate of the
log hazard ratio, the hazard ratio, the standard error on the log scale,
the Wald z-statistic, the p-value that follows the side recorded in
the object, and the Wald confidence interval for the hazard ratio.
The p-value column is labeled Pr(>|z|) for a two-sided test and
Pr(<z) for a one-sided test, whose p-value is the lower tail in the
direction of treatment benefit.
## S3 method for class 'coxph_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'coxph_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments passed to |
Invisibly returns x.
library(survival) fit <- coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1) print(fit)library(survival) fit <- coxph_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1) print(fit)
Prints a compact per-group summary of a kmcurve_fast object: the group
label and role, the sample size, the number of events, and the median
survival read from the Kaplan-Meier curve with the convention of
survival::survfit. In the two-group case the constant
Cox proportional-hazards estimate of the hazard ratio is reported below the
table.
## S3 method for class 'kmcurve_fast' print(x, ...)## S3 method for class 'kmcurve_fast' print(x, ...)
x |
A |
... |
Further arguments, currently ignored. |
The input x, invisibly.
kmcurve_fast, plot.kmcurve_fast
Formats and prints a maxcombo_fast object. It shows the total sample
size, a table of the Fleming-Harrington component Z-scores, and the
max-combo statistic with its p-value. For a one-sided test the statistic is
the most negative component, so a negative value favors the treatment
group, matching the sign convention of survdiff_fast. For a
two-sided test the statistic is the largest component in absolute value.
## S3 method for class 'maxcombo_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'maxcombo_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
library(survival) fit <- maxcombo_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1) print(fit)library(survival) fit <- maxcombo_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1) print(fit)
Formats and prints a medsurv_fast object in the same layout as the
other two-group estimation summaries in the package. The header shows the
control label and the inference settings. The body shows the per-group
median survival with its confidence interval, followed, for a two-group
object, by the difference contrast (treatment minus control) with a
confidence interval, the test statistic, and the p-value.
The p-value column is labeled Pr(>|z|) for a two-sided test and
Pr(>z) for a one-sided test, whose p-value is the upper tail in the
direction of treatment benefit.
## S3 method for class 'medsurv_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'medsurv_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
set.seed(1) time <- c(rexp(50, 0.1), rexp(50, 0.07)) status <- rep(1, 100) group <- rep(c(0, 1), each = 50) print(medsurv_fast(time, status, group, control = 0))set.seed(1) time <- c(rexp(50, 0.1), rexp(50, 0.07)) status <- rep(1, 100) group <- rep(c(0, 1), each = 50) print(medsurv_fast(time, status, group, control = 0))
Formats and prints a milestone_fast object in the same layout as the
other two-group estimation summaries in the package. The header shows the
milestone timepoint, the control label, and the inference settings. The body
shows the per-group milestone survival with its confidence interval, followed
by the difference contrast (treatment minus control) with a confidence
interval, the test statistic, and the p-value. The p-value follows the
alternative recorded in the object. Its column is labeled
Pr(>|z|) for a two-sided test, and for a one-sided test
Pr(>z) (upper tail) with "wald" and "mover" or
Pr(<z) (lower tail) with "loglog", the direction of treatment
benefit.
## S3 method for class 'milestone_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'milestone_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
set.seed(1) time <- c(rexp(50, 0.1), rexp(50, 0.07)) status <- rep(1, 100) group <- rep(c(0, 1), each = 50) print(milestone_fast(time, status, group, control = 0, tau = 10, method = "loglog"))set.seed(1) time <- c(rexp(50, 0.1), rexp(50, 0.07)) status <- rep(1, 100) group <- rep(c(0, 1), each = 50) print(milestone_fast(time, status, group, control = 0, tau = 10, method = "loglog"))
Formats and prints an rmst_fast object. In single-group mode it shows
the restricted mean survival time, its Greenwood standard error, and the
Wald confidence interval at the requested horizon. In two-group mode it
shows the per-group restricted mean survival times together with the
difference (treatment minus control) and ratio (treatment over control)
contrasts, each with a Wald z-statistic and a p-value that follows the
side recorded in the object.
The p-value column is labeled Pr(>|z|) for a two-sided test and
Pr(>z) for a one-sided test, whose p-value is the upper tail in the
direction of treatment benefit.
## S3 method for class 'rmst_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'rmst_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
set.seed(42) t_raw <- rexp(100, rate = 1 / 10) e_raw <- rbinom(100, 1, 0.7) # Single-group print(rmst_fast(t_raw, e_raw, tau = 10)) # Two-group comparison set.seed(7) n <- 200 time <- c(rexp(n, 0.10), rexp(n, 0.07)) event <- rbinom(2 * n, 1, 0.8) group <- rep(0:1, each = n) print(rmst_fast(time, event, group = group, control = 0, tau = 10))set.seed(42) t_raw <- rexp(100, rate = 1 / 10) e_raw <- rbinom(100, 1, 0.7) # Single-group print(rmst_fast(t_raw, e_raw, tau = 10)) # Two-group comparison set.seed(7) n <- 200 time <- c(rexp(n, 0.10), rexp(n, 0.07)) event <- rbinom(2 * n, 1, 0.8) group <- rep(0:1, each = n) print(rmst_fast(time, event, group = group, control = 0, tau = 10))
Formats and prints an rmw_fast object, showing the standardized
log-rank and modestly-weighted component Z-scores, their null correlation,
the survival-probability threshold s_star, the combined test
statistic, and the corresponding p-value.
## S3 method for class 'rmw_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'rmw_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
library(survival) fit <- rmw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1, s_star = 0.5) print(fit)library(survival) fit <- rmw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1, s_star = 0.5) print(fit)
Prints a compact one-row-per-scenario summary of a scenario_fast
object: the scenario label, the null flag, the total sample size, the control
and experimental medians, the hazard ratio at the start and end of the time
window, and whether the two survival curves cross. The medians and hazard
ratios are read from the analytic survival curves, so they are defined for
both exponential and piecewise-exponential scenarios.
## S3 method for class 'scenario_fast' print(x, tmax = NULL, ...)## S3 method for class 'scenario_fast' print(x, tmax = NULL, ...)
x |
A |
tmax |
Upper limit of the time window used to read medians, the end-of- window hazard ratio, and the crossing flag. Defaults to a value derived from the scenarios. |
... |
Further arguments, currently ignored. |
The input x, invisibly.
gen_scenario_fast, plot.scenario_fast
Formats an object returned by simsummary_fast in the style of a
group-sequential design report. After a header with the simulation count and
the boundary settings, two look-by-look tables are shown: a stopping-boundary
table (information fraction, events, sample size, the efficacy and futility
boundaries, and the cumulative efficacy crossing probability) and an
analysis-timing table (sample size, events, dropouts, pipeline, analysis time,
and the per-look efficacy and futility crossing probabilities). An overall
block reports the rejection rate and the expected counts and timing at the
stopping look.
## S3 method for class 'simsummary_fast' print(x, digits = 4, ...)## S3 method for class 'simsummary_fast' print(x, digits = 4, ...)
x |
An object of class |
digits |
A single positive integer, the number of decimal places used for the printed probabilities. Defaults to 4. |
... |
Further arguments, currently ignored. |
Column labels follow the convention of group-sequential design software:
Events (s), Sample (n), Dropouts (d), Pipeline
(the enrolled count minus events minus dropouts), Analysis Time (mean
calendar time), and Info. Frac. (the information fraction, computed as
the mean events at a look divided by the mean events at the final look, or
from look.value when no event count is available). Probabilities are
printed to digits decimal places and counts and times to fewer. Because
the summary is a Monte Carlo estimate under a single data-generating truth, it
does not carry the separate null and alternative columns or the alpha and beta
spending of an analytic design report. When the summary has several blocks
(arms of pairwise_fast output or subgroup populations), each
block is printed under its own heading. The object is a data frame with an
additional class, so the unrounded values remain available by subsetting
it. A selection of rows that keeps, in every block, all the looks in their
original order (for example one population or one arm) is printed as a
report. Any other selection of rows, such as some of the looks or rows in a
different order, is printed as an ordinary data frame, because the
boundaries and the information fractions refer to all the looks of the
design. A selection of columns keeps the class but not the boundary
settings, and is printed as an ordinary data frame, as is a selection
without any look row. Only a selection by x[i, ] keeps the boundary
settings; a selection by subset() or by x[i, j], even with all
the columns, is printed as an ordinary data frame. Results combined by
rbind() keep the attributes of the first result only, which include a
copy of its rows. The report is therefore printed only when the rows of
every block equal those of the first result; rows that differ, such as
those of a population summarized with other boundaries, make the whole
object print as an ordinary data frame, even when the block has the same
name in both results.
The object x, invisibly.
simsummary_fast, plot.simsummary_fast.
Formats and prints a survdiff_fast object similarly to
print(survival::survdiff(...)), showing the observed and expected
event counts for the control and treatment groups, the per-group
contributions (O-E)^2 / E and (O-E)^2 / V, the test
statistic, and the corresponding p-value. For a weighted log-rank test the
header names the weight scheme (and notes stratification), and only the
observed event counts are shown alongside the weighted statistic, since a
single unweighted expected count is not defined for a weighted test.
## S3 method for class 'survdiff_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'survdiff_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
library(survival) fit <- survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 2) print(fit)library(survival) fit <- survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 2) print(fit)
Formats and prints a survfit_fast object similarly to
print(summary(survival::survfit(...), times = t_eval)), showing the
Kaplan-Meier survival estimate, the Greenwood standard error on the
survival scale, and the confidence interval at the requested evaluation
time.
## S3 method for class 'survfit_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'survfit_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
set.seed(42) t_raw <- rexp(100, rate = 1 / 10) e_raw <- rbinom(100, 1, 0.7) ord <- order(t_raw) fit <- survfit_fast(t_raw[ord], e_raw[ord], t_eval = 10) print(fit)set.seed(42) t_raw <- rexp(100, rate = 1 / 10) e_raw <- rbinom(100, 1, 0.7) ord <- order(t_raw) fit <- survfit_fast(t_raw[ord], e_raw[ord], t_eval = 10) print(fit)
Formats and prints a wkm_fast object in the same layout as the other
two-group summaries in the package. The header shows the control label and
the inference settings. The body shows the weighted integrated survival
difference (treatment minus control) with a confidence interval, the test
statistic, and the p-value.
The p-value column is labeled Pr(>|z|) for a two-sided test and
Pr(>z) for a one-sided test, whose p-value is the upper tail in the
direction of treatment benefit.
## S3 method for class 'wkm_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'wkm_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
set.seed(1) n <- 200 g <- rep(0:1, each = n / 2) tt <- c(rexp(n / 2, 0.1), rexp(n / 2, 0.07)) cc <- rexp(n, 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) print(wkm_fast(time, event, group = g, control = 0))set.seed(1) n <- 200 g <- rep(0:1, each = n / 2) tt <- c(rexp(n / 2, 0.1), rexp(n / 2, 0.07)) cc <- rexp(n, 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) print(wkm_fast(time, event, group = g, control = 0))
Formats and prints a wmst_fast object in the same layout as the other
two-group summaries in the package. The header shows the window, the control
label and the inference settings. The body shows the per-group window mean
survival time with its confidence interval, followed, for a two-group object,
by the difference contrast (treatment minus control) with a confidence
interval, the test statistic, and the p-value.
The p-value column is labeled Pr(>|z|) for a two-sided test and
Pr(>z) for a one-sided test, whose p-value is the upper tail in the
direction of treatment benefit.
## S3 method for class 'wmst_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)## S3 method for class 'wmst_fast' print(x, digits = max(1L, getOption("digits") - 3L), ...)
x |
An object of class |
digits |
Number of significant digits to display. Defaults to the
global option |
... |
Additional arguments (currently unused). |
Invisibly returns x.
set.seed(1) n <- 200 g <- rep(0:1, each = n / 2) tt <- c(rexp(n / 2, 0.1), rexp(n / 2, 0.07)) cc <- rexp(n, 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) print(wmst_fast(time, event, group = g, control = 0, tau1 = 2, tau2 = 12))set.seed(1) n <- 200 g <- rep(0:1, each = n / 2) tt <- c(rexp(n / 2, 0.1), rexp(n / 2, 0.07)) cc <- rexp(n, 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) print(wmst_fast(time, event, group = g, control = 0, tau1 = 2, tau2 = 12))
Computes the restricted mean survival time (RMST) up to a horizon tau
from the Kaplan-Meier estimator. With a single group (the default), it
returns the RMST with its Greenwood-type standard error and a Wald
confidence interval. When a group is supplied, it additionally
returns the two-group contrasts: the RMST difference (treatment minus
control) and the RMST ratio (treatment over control), each with a standard
error, confidence interval, and test (one- or two-sided, following
side). The C++ backend integrates
the survival step function in a single scan and is reused once per group, so
the function is suitable for simulation loops with presorted = TRUE.
rmst_fast( time, event, group = NULL, control = NULL, side = 2, conf.level = 0.95, tau, presorted = FALSE, data = NULL )rmst_fast( time, event, group = NULL, control = NULL, side = 2, conf.level = 0.95, tau, presorted = FALSE, data = NULL )
time |
A numeric vector of event or censoring times.
Alternatively, a formula |
event |
An integer or numeric vector of event indicators
(1 = event, 0 = censored), aligned with |
group |
An optional vector of group labels aligned with |
control |
A scalar value indicating which level of |
side |
1 for a one-sided test in the direction of treatment benefit
(longer restricted mean survival in the treatment group) or 2 for a
two-sided test (default 2). Applies to the two-group contrasts; the
confidence intervals are always two-sided at |
conf.level |
A single numeric value in (0, 1) specifying the confidence level. Defaults to 0.95. |
tau |
A single positive numeric value specifying the restriction horizon. A warning is given when it exceeds the largest observed time (of either group in the two-group case), because the Kaplan-Meier curve is not estimated beyond that time. |
presorted |
A logical value. If |
data |
A data frame in which the variables of a formula given as
|
The RMST is the area under the Kaplan-Meier curve from 0 to tau:
RMST(tau) = integral over [0, tau] of S(u) du.
The variance follows the Greenwood-type estimator
Var[RMST] = sum_{t_i <= tau} A_i^2 d_i / (n_i (n_i - d_i)),
where A_i = integral over [t_i, tau] of S(u) du is the area to the right of
event time t_i. This matches the restricted mean reported by
survfit and by the survRM2 package.
When group is supplied, the treatment group is the level of
group that is not equal to control. Writing the treatment and
control RMST as r1 and r0 with variances v1 and v0, the difference contrast
is diff = r1 - r0 with Var[diff] = v1 + v0, using the independence of the
two groups. The ratio contrast is formed on the log scale by the delta
method, Var[log(r1 / r0)] = v1 / r1^2 + v0 / r0^2, with the confidence
interval exponentiated back to the ratio scale. These match the unadjusted
contrasts reported by survRM2.
When presorted = TRUE, the input vectors are assumed to be sorted in
ascending order of time; splitting into groups preserves the
ascending order within each group, so no re-sorting is performed. When
presorted = FALSE (default), sorting is handled internally. In
simulation loops where the data are generated in sorted order,
presorted = TRUE avoids one O(n log n) pass.
An object of class "rmst_fast", a named numeric vector. In
single-group mode it has length 4 with elements rmst,
std.err, lower, and upper. In two-group mode it
contains the per-group RMST (rmst.ctrl, rmst.trt), the
difference contrast (diff, se.diff, diff.lower,
diff.upper, z.diff, p.diff), and the ratio contrast
(ratio, ratio.lower, ratio.upper, z.ratio,
p.ratio). The restriction horizon, confidence level, and test side
are stored as attributes tau, conf.level, and side;
in two-group mode the control label is also stored. Returns
NA_real_ values (still with class "rmst_fast") when n
is zero in single-group mode.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Royston, P., & Parmar, M. K. B. (2013). Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome. BMC Medical Research Methodology, 13, 152.
Uno, H., Claggett, B., Tian, L., et al. (2014). Moving beyond the hazard ratio in quantifying the between-group difference in survival analysis. Journal of Clinical Oncology, 32(22), 2380-2385.
survfit_fast for the Kaplan-Meier estimate at a single time
point.
set.seed(42) t_raw <- rexp(100, rate = 1 / 10) e_raw <- rbinom(100, 1, 0.7) # Single-group RMST rmst_fast(t_raw, e_raw, tau = 10) # Single-group, pre-sorted (sort once, reuse in a loop) ord <- order(t_raw) rmst_fast(t_raw[ord], e_raw[ord], tau = 10, presorted = TRUE) # Two-group comparison (difference and ratio) set.seed(7) n <- 200 time <- c(rexp(n, 0.10), rexp(n, 0.07)) event <- rbinom(2 * n, 1, 0.8) group <- rep(0:1, each = n) rmst_fast(time, event, group = group, control = 0, tau = 10) # Formula interface, and the result as a data frame d <- data.frame(time, event, group) as.data.frame(rmst_fast(Surv(time, event) ~ group, data = d, tau = 10)) # Validation against survRM2 if (requireNamespace("survRM2", quietly = TRUE)) { survRM2::rmst2(time, event, group, tau = 10)$unadjusted.result }set.seed(42) t_raw <- rexp(100, rate = 1 / 10) e_raw <- rbinom(100, 1, 0.7) # Single-group RMST rmst_fast(t_raw, e_raw, tau = 10) # Single-group, pre-sorted (sort once, reuse in a loop) ord <- order(t_raw) rmst_fast(t_raw[ord], e_raw[ord], tau = 10, presorted = TRUE) # Two-group comparison (difference and ratio) set.seed(7) n <- 200 time <- c(rexp(n, 0.10), rexp(n, 0.07)) event <- rbinom(2 * n, 1, 0.8) group <- rep(0:1, each = n) rmst_fast(time, event, group = group, control = 0, tau = 10) # Formula interface, and the result as a data frame d <- data.frame(time, event, group) as.data.frame(rmst_fast(Surv(time, event) ~ group, data = d, tau = 10)) # Validation against survRM2 if (requireNamespace("survRM2", quietly = TRUE)) { survRM2::rmst2(time, event, group, tau = 10)$unadjusted.result }
Computes the robust modestly-weighted (rMW) log-rank test of Magirr and
Öhrn, which combines the standard log-rank test with a single
modestly-weighted log-rank test. The test statistic is the maximum of the
two standardized components, evaluated against their joint null distribution.
Because the standard log-rank statistic is included as one of the two
components and the components are strongly correlated under the null, the
multiplicity adjustment is small, so the test loses little power relative to
the log-rank test in worst-case scenarios while gaining substantial power
under delayed effects. The two component numerators, their variances, and
their null covariance are computed in a single pass by the C++ core
rmw_core.
rmw_fast( time, event, group, control, side = 2, s_star = 0.5, presorted = FALSE, data = NULL )rmw_fast( time, event, group, control, side = 2, s_star = 0.5, presorted = FALSE, data = NULL )
time |
A numeric vector of follow-up times for all subjects.
Alternatively, a formula |
event |
An integer or numeric vector of event indicators
(1 = event, 0 = censored), aligned with |
group |
A vector of group labels aligned with |
control |
A scalar value indicating which level of |
side |
An integer, either 1 or 2 (default 2). If |
s_star |
A single numeric value in |
presorted |
A logical value. If |
data |
A data frame in which the variables of a formula given as
|
The first component is the ordinary log-rank statistic, with weight one at
every event time. The second component is a modestly-weighted log-rank
statistic with weight min(1 / S(t-), 1 / s_star), where S(t-)
is the left-continuous pooled Kaplan-Meier estimate just prior to each event
time and s_star is a survival-probability threshold. This is the
survival-threshold parameterization of the modestly-weighted test of Magirr
and Burman, in which the weight is capped at 1 / s_star. It differs
from the timepoint parameterization exposed by survdiff_fast() with
weight = "mwlrt", where the cap is derived from a timepoint
t_star. The choice s_star = 0.5 caps the weight at 2 and is the
value used in the original rMW article.
Writing the two standardized components as Z_lr and Z_mw, under
the null hypothesis of equal survival the pair is asymptotically bivariate
normal with zero means, unit variances, and correlation
rho = C / sqrt(V_lr V_mw), where C is the covariance of the two
numerators returned by the C++ core. When side = 1 the statistic is
min(Z_lr, Z_mw), so a protective treatment effect (fewer events than
expected) yields a small, negative value, and the one-sided p-value is
P(min(Z_lr, Z_mw) <= observed) under the joint null. When
side = 2 the statistic is max(abs(Z_lr), abs(Z_mw)) and the
two-sided p-value is P(max(abs(Z_lr), abs(Z_mw)) >= observed). The
joint normal probability is evaluated with mvtnorm::pmvnorm, using the
exact TVPACK algorithm for the one-sided region, an orthant (all
upper limits infinite), and the
deterministic Miwa algorithm for the two-sided rectangle.
When presorted = TRUE, the input vectors are assumed to be sorted by
ascending time and the internal order() call is skipped, which
avoids one O(n log n) pass in simulation loops where the data are already
generated in time order.
An object of class "rmw_fast", a length-two numeric vector
c(statistic, p.value) with attributes z (the named component
Z-scores c(logrank, mwlrt)), corr (the 2 by 2 null
correlation matrix of the two components), s_star, O1 (the
observed number of events in the treatment group), side, and
n (the total sample size). Returns NA statistic and p-value
(still with class "rmw_fast") when either component variance is zero
(e.g., no event time at which both groups are at risk).
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Magirr, D., & Öhrn, F. (2026). Robust modestly weighted log-rank tests. Pharmaceutical Statistics, 25(1), e70066.
Magirr, D., & Burman, C.-F. (2019). Modestly weighted logrank tests. Statistics in Medicine, 38(20), 3782-3790.
survdiff_fast for the standard and weighted log-rank tests.
print.rmw_fast for the print method.
library(survival) # One-sided robust modestly-weighted test with s_star = 0.5 fit <- rmw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1, s_star = 0.5) fit # The log-rank component matches survdiff_fast with weight = "logrank" z_lr <- as.numeric(survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1)) cat("rmw_fast log-rank component:", attr(fit, "z")[["logrank"]], "\n") cat("survdiff_fast log-rank Z :", z_lr, "\n") # Two-sided test rmw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 2, s_star = 0.5) # presorted = TRUE: sort once outside, reuse inside a loop ord <- order(ovarian$futime) rmw_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord], control = 1, side = 1, s_star = 0.5, presorted = TRUE) # Formula interface rmw_fast(Surv(futime, fustat) ~ rx, data = ovarian, side = 1)library(survival) # One-sided robust modestly-weighted test with s_star = 0.5 fit <- rmw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1, s_star = 0.5) fit # The log-rank component matches survdiff_fast with weight = "logrank" z_lr <- as.numeric(survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 1)) cat("rmw_fast log-rank component:", attr(fit, "z")[["logrank"]], "\n") cat("survdiff_fast log-rank Z :", z_lr, "\n") # Two-sided test rmw_fast(ovarian$futime, ovarian$fustat, ovarian$rx, control = 1, side = 2, s_star = 0.5) # presorted = TRUE: sort once outside, reuse inside a loop ord <- order(ovarian$futime) rmw_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord], control = 1, side = 1, s_star = 0.5, presorted = TRUE) # Formula interface rmw_fast(Surv(futime, fustat) ~ rx, data = ovarian, side = 1)
Simulates time-to-event trial data for one or two groups across many
simulated trials, with piecewise accrual, piecewise-exponential survival and
dropout, and optional subgroups defined by a prevalence specification. The
entire generation pipeline (accrual, survival, dropout, derived columns, and
two-group interleaving) runs in a single C++ kernel that materializes the
output data frame once, avoiding intermediate R-level vector operations and
copies. The random-number stream is consumed in the same order as a per-group
reference implementation, so results are reproducible from seed.
simdata_fast( nsim = 1000, n, alloc = c(1, 1), a.time, a.rate = NULL, a.prop = NULL, e.hazard = NULL, e.median = NULL, e.time = NULL, d.hazard = NULL, d.median = NULL, d.time = NULL, seed = NULL, prevalence = NULL, fixed.alloc = FALSE, h01.hazard = NULL, h01.median = NULL, h01.time = NULL, h02.hazard = NULL, h02.median = NULL, h02.time = NULL, h12.hazard = NULL, h12.median = NULL, h12.time = NULL, switch.prop = NULL, h12.switch.hazard = NULL, h12.switch.median = NULL, h12.switch.time = NULL, switch.clock = "reset", stream = NULL )simdata_fast( nsim = 1000, n, alloc = c(1, 1), a.time, a.rate = NULL, a.prop = NULL, e.hazard = NULL, e.median = NULL, e.time = NULL, d.hazard = NULL, d.median = NULL, d.time = NULL, seed = NULL, prevalence = NULL, fixed.alloc = FALSE, h01.hazard = NULL, h01.median = NULL, h01.time = NULL, h02.hazard = NULL, h02.median = NULL, h02.time = NULL, h12.hazard = NULL, h12.median = NULL, h12.time = NULL, switch.prop = NULL, h12.switch.hazard = NULL, h12.switch.median = NULL, h12.switch.time = NULL, switch.clock = "reset", stream = NULL )
nsim |
Number of simulated trials. |
n |
Either a single total sample size (split by |
alloc |
A length-two allocation ratio (control, treatment), used when
|
a.time |
A numeric vector of accrual-interval breakpoints. |
a.rate |
Absolute accrual rates (subjects per unit time), interpreted in
one of two ways. With length |
a.prop |
Accrual proportions, one per accrual interval (length
|
e.hazard |
Survival hazard(s). A scalar or vector (piecewise, with
|
e.median |
Survival median(s); an alternative to |
e.time |
Survival breakpoints for piecewise hazards (last element
|
d.hazard |
Dropout hazard(s), same structure as |
d.median |
Dropout median(s); an alternative to |
d.time |
Dropout breakpoints for piecewise hazards. |
seed |
Optional integer seed for the |
prevalence |
Optional subgroup prevalence specification (numeric
vector, list of vectors, array, or a named |
fixed.alloc |
Logical; when |
h01.hazard |
Transition hazard(s) for the non-terminal (intermediate)
event (state 0 to state 1) in the illness-death model. A scalar or vector
for one group, or a two-element list for two groups. Supplying any of
|
h01.median |
Median(s) for the intermediate event; an alternative to
|
h01.time |
Breakpoints for a piecewise |
h02.hazard |
Transition hazard(s) for the terminal event without a prior
intermediate event (state 0 to state 2). Same group and piecewise
conventions as |
h02.median |
Median(s) for the direct terminal event; an alternative to
|
h02.time |
Breakpoints for a piecewise |
h12.hazard |
Transition hazard(s) for the terminal event after an
intermediate event (state 1 to state 2) for subjects who do not switch.
Defaults to |
h12.median |
Median(s) for the post-event terminal event; an alternative
to |
h12.time |
Breakpoints for a piecewise |
switch.prop |
Probability that a subject whose intermediate event is observed (on or before dropout) switches treatment at that event, a scalar or a two-element list (control, treatment). Defaults to zero (no switching); the treatment group is typically left at zero. |
h12.switch.hazard |
Transition hazard(s) from state 1 to state 2 for
subjects who switch. Required when any |
h12.switch.median |
Median(s) for the post-switch terminal event; an
alternative to |
h12.switch.time |
Breakpoints for a piecewise |
switch.clock |
Time origin for the post-event hazards. Currently only
|
stream |
Optional non-negative whole number selecting an independent
|
For each subject the observed time-to-event is
tte = pmin(surv_time, dropout_time) and event is 1 when the
survival time occurs first. The calendar time of the observed event is
accrual_time + tte.
A hazard of zero in the last piece (for example a cure fraction, as in
e.hazard = c(0.1, 0) with e.time = c(0, 24, Inf)) gives an
infinite latent time to some subjects. A subject with neither a finite
survival time nor a finite dropout time has tte = Inf and
event = 0, and an analysis at a finite cutoff censors the subject at
the cutoff.
The total enrolled is fixed at sum(n). With a.rate the rates are
absolute (subjects per unit time): when the accrual period is fully specified
the rates must accrue exactly sum(n), and when one extra rate is given
the end of the final interval is solved so the total is met. With a.prop
the values are relative proportions that distribute sum(n) across the
fully specified intervals. Each accrual interval receives a deterministic
number of subjects (the rate or proportion times the group total, rounded to
keep the per-group total exact), placed uniformly within the interval.
Survival and dropout are exponential when a single hazard (or median) is
supplied and piecewise-exponential when a vector is supplied together with
the corresponding e.time or d.time breakpoints, whose last
element must be Inf. Group-specific parameters are supplied as a
two-element list (control first, treatment second).
When prevalence is supplied the trial has subgroups. A numeric vector
defines a single factor; a list of numeric vectors defines several
independent factors; a multi-dimensional array defines the joint
distribution of correlated factors. Per-cell hazards may be supplied as a
list with one element per cell (cells in column-major order, the first
factor varying fastest). With fixed.alloc = TRUE the subgroup
sizes are deterministic and the fixed labels are assigned to subjects in a
random order within each simulated trial, so subgroup membership does not
depend on the accrual time; otherwise subgroup membership is drawn from the
prevalence distribution.
The number of groups is determined as follows. A length-two n always
gives two groups. A scalar n gives two groups, with the total split
by alloc, in any of these cases: alloc is supplied explicitly;
the prevalence is group-specific; without subgroups, any of e.hazard,
e.median, d.hazard, and d.median is a list; with
subgroups, one of them is a length-two list with a list element (per-cell
values within a group, as in list(list(0.10, 0.08, 0.06), 0.05)).
Otherwise a scalar n gives one group, and with subgroups a list is
read as one element per cell of that group. So with subgroups and a scalar
n, e.hazard = list(0.10, 0.05) means per-cell hazards of a
single group unless alloc is supplied, in which case it means
per-group hazards. In a two-group simulation each specification is either
shared by both groups (not a list) or a list of length two (control first),
and each group's element may be a per-cell list when there are subgroups.
Without subgroups, a length-two n requires at least one of the
survival and dropout specifications to be such a list.
The breakpoints e.time and d.time follow the same structure
(per group in a two-group simulation, per cell in a one-group simulation with
subgroups).
When n is a vector of length greater than two together with a per-arm
survival list, the simulation is a multi-arm trial. Each arm is generated in
turn with the validated single-group kernel over a common accrual window, and
the arms are stacked into one data frame with a group column labeled 1
to length(n) in the order of n. Per-arm survival is supplied as
an e.hazard or e.median list with one element per arm, and
optional dropout as a shared value or a per-arm list through d.hazard
or d.median. The arms share the master seed, so the result is
reproducible. A multi-arm design is analyzed as a set of pairwise contrasts by
subsetting the output to the control arm and one other arm and calling
analysis_fast once per contrast, which
pairwise_fast does for every experimental arm. Multi-arm mode does not support
subgroups or the illness-death model, which remain two-group.
The illness-death model (activated by h01.* or h02.*) also
accepts prevalence and fixed.alloc. The transition hazards
(h01.*, h02.*, h12.*, h12.switch.*),
switch.prop, and the dropout specification then follow the rules of
e.hazard with subgroups: in a two-group simulation each is shared or
a list with one element per group, and each group's element is shared or a
list with one element per subgroup cell, for example
h01.hazard = list(list(0.10, 0.06), 0.05); with a scalar n
and no alloc, a list without list elements holds one element per
cell of a single group. All subjects share one accrual process and the
subgroup cells are assigned after accrual, so the subgroups enroll over the
same calendar in proportion to their prevalence. This gives a mixture of
illness-death models without generating the subgroups separately and
adjusting their accrual rates. With a single cell (for example
prevalence = 1) the data equal those without prevalence, with
the subgroup column added.
A data.frame with nsim * sum(n) rows. The columns are
sim, group, any subgroup columns, accrual_time,
surv_time, dropout_time, tte, event, and
calendar_time.
In the illness-death model the columns are instead sim,
group, any subgroup columns, accrual_time,
e1_surv_time,
e2_surv_time, dropout_time, e1_tte, e1_event,
e2_tte, e2_event, e1_calendar_time,
e2_calendar_time, intermediate, switched, and
switch_time, where e1 is the first (state-0 exit) endpoint
and e2 is the terminal endpoint. In oncology e1 is
progression-free survival, e2 is overall survival, and
intermediate flags progression. The column intermediate
describes the latent process (an intermediate event before the terminal
event, e1_surv_time < e2_surv_time), also when it would occur
after dropout; the observed intermediate event is intermediate == 1
with e1_event == 1. A switch (switched == 1, at
switch_time, the intermediate-event time) occurs only after an
observed intermediate event.
The output has nsim * sum(n) rows, so memory rather than time limits
the number of simulated trials in one call. A large study can be run in
batches: each batch is generated with the same seed and its own
stream, analyzed with analysis_fast, and only the
analysis results are kept. Because each stream is an independent sequence of
the dqrng generator, the data of batch b are the same whether
the batches are run one after another or distributed over parallel workers
(for example with parallel::mclapply() or the future framework), and
whatever the number of workers. The simulation identifiers start at 1 in
every batch, so they are renumbered before the batch results are combined.
Functions that draw further random numbers, such as switch_fast,
continue the same stream within a batch. The max-combo p-values of
analysis_fast with four or more weights (or a two-sided test)
are computed with R's own random-number generator, so a batch should also
call set.seed() before the analysis for p-values that are identical
across runs.
Fleischer, F., Gaschler-Markefski, B., & Bluhmki, E. (2009). A statistical model for the dependence between progression-free survival and overall survival. Statistics in Medicine, 28(21), 2669-2686.
analysis_fast, cutoff_fast,
switch_fast, pairwise_fast
# Batches from independent random-number streams: the simulations of each # batch are renumbered before the analysis results are combined. run_batch <- function(b, nsim = 20) { d <- simdata_fast(nsim = nsim, n = c(100, 100), a.time = c(0, 12), a.rate = 200 / 12, e.median = list(12, 16), seed = 2026, stream = b) r <- analysis_fast(d, control = 1, event.looks = 120) r$sim <- r$sim + (b - 1) * nsim r } res <- do.call(rbind, lapply(1:3, run_batch)) range(res$sim) # One-group simulation, simple exponential, no dropout df1 <- simdata_fast( nsim = 100, n = 50, a.time = c(0, 12), a.rate = 50 / 12, e.median = 18, seed = 1 ) head(df1) # Accrual rate with the final interval computed from the total: 20 per unit # time for the first 12 units, then 30 per unit time until 500 are enrolled df1b <- simdata_fast( nsim = 100, n = 500, a.time = c(0, 12), a.rate = c(20, 30), e.median = 18, seed = 1 ) head(df1b) # Accrual by proportion: 30% enrolled in [0, 6], 70% in [6, 12] df1c <- simdata_fast( nsim = 100, n = 50, a.time = c(0, 6, 12), a.prop = c(0.3, 0.7), e.median = 18, seed = 1 ) head(df1c) # Two-group simulation, simple exponential, with dropout df2 <- simdata_fast( nsim = 100, n = c(60, 60), a.time = c(0, 6, 12), a.rate = c(8, 12), e.median = list(18, 24), d.hazard = list(0.01, 0.01), seed = 2 ) head(df2) # One factor with three levels: single subgroup column df3 <- simdata_fast( nsim = 100, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(list(0.10, 0.08, 0.06), 0.05), prevalence = c(0.5, 0.3, 0.2), seed = 3 ) head(df3) # The same trial specified by the total sample size and the allocation ratio df3b <- simdata_fast( nsim = 100, n = 300, alloc = c(1, 1), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(list(0.10, 0.08, 0.06), 0.05), prevalence = c(0.5, 0.3, 0.2), seed = 3 ) identical(df3, df3b) # Two independent factors (2 x 2): columns subgroup1 and subgroup2. # Four cells in column-major order: (1,1), (2,1), (1,2), (2,2). df4 <- simdata_fast( nsim = 100, n = 200, a.time = c(0, 12), a.rate = 200 / 12, e.hazard = list(0.10, 0.08, 0.07, 0.05), prevalence = list(c(0.5, 0.5), c(0.6, 0.4)), seed = 4 ) head(df4) # Two correlated factors via a joint-distribution array (2 x 2) df5 <- simdata_fast( nsim = 100, n = 200, a.time = c(0, 12), a.rate = 200 / 12, e.hazard = 0.08, prevalence = array(c(0.40, 0.10, 0.15, 0.35), dim = c(2, 2)), seed = 5 ) head(df5) # Two correlated endpoints, no switching (the maximal-independence model of # Fleischer et al., 2009, because h12 defaults to h02). In oncology e1 is PFS # and e2 is OS; here the # control has faster intermediate events and faster direct terminal events. dfid <- simdata_fast( nsim = 100, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, h01.median = list(8, 12), h02.median = list(24, 30), seed = 6 ) head(dfid) # With treatment switching: 40 percent of control subjects who reach the # intermediate event switch and then follow a more favorable post-event hazard dfsw <- simdata_fast( nsim = 100, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, h01.median = list(8, 12), h02.median = list(24, 30), switch.prop = list(0.4, 0), h12.switch.median = list(36, 36), seed = 7 ) head(dfsw) # Illness-death model with two subgroups (prevalence 0.4 and 0.6): the # treatment effect on progression is larger in subgroup 2. dfsg <- simdata_fast( nsim = 100, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, h01.median = list(8, list(10, 14)), h02.median = list(24, 30), prevalence = c(0.4, 0.6), seed = 8 ) table(dfsg$group, dfsg$subgroup) / 100 # Three-arm trial (one control and two treatment arms) analyzed as pairwise # contrasts against the shared control. dfk <- simdata_fast( nsim = 100, n = c(120, 120, 120), a.time = c(0, 12), a.rate = 360 / 12, e.median = list(12, 16, 20), seed = 8 ) # Control arm (group 1) versus treatment arm 2, one-sided log-rank at month 24. sub12 <- dfk[dfk$group %in% c(1, 2), ] res12 <- analysis_fast(sub12, control = 1, time.looks = 24, side = 1) head(res12)# Batches from independent random-number streams: the simulations of each # batch are renumbered before the analysis results are combined. run_batch <- function(b, nsim = 20) { d <- simdata_fast(nsim = nsim, n = c(100, 100), a.time = c(0, 12), a.rate = 200 / 12, e.median = list(12, 16), seed = 2026, stream = b) r <- analysis_fast(d, control = 1, event.looks = 120) r$sim <- r$sim + (b - 1) * nsim r } res <- do.call(rbind, lapply(1:3, run_batch)) range(res$sim) # One-group simulation, simple exponential, no dropout df1 <- simdata_fast( nsim = 100, n = 50, a.time = c(0, 12), a.rate = 50 / 12, e.median = 18, seed = 1 ) head(df1) # Accrual rate with the final interval computed from the total: 20 per unit # time for the first 12 units, then 30 per unit time until 500 are enrolled df1b <- simdata_fast( nsim = 100, n = 500, a.time = c(0, 12), a.rate = c(20, 30), e.median = 18, seed = 1 ) head(df1b) # Accrual by proportion: 30% enrolled in [0, 6], 70% in [6, 12] df1c <- simdata_fast( nsim = 100, n = 50, a.time = c(0, 6, 12), a.prop = c(0.3, 0.7), e.median = 18, seed = 1 ) head(df1c) # Two-group simulation, simple exponential, with dropout df2 <- simdata_fast( nsim = 100, n = c(60, 60), a.time = c(0, 6, 12), a.rate = c(8, 12), e.median = list(18, 24), d.hazard = list(0.01, 0.01), seed = 2 ) head(df2) # One factor with three levels: single subgroup column df3 <- simdata_fast( nsim = 100, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(list(0.10, 0.08, 0.06), 0.05), prevalence = c(0.5, 0.3, 0.2), seed = 3 ) head(df3) # The same trial specified by the total sample size and the allocation ratio df3b <- simdata_fast( nsim = 100, n = 300, alloc = c(1, 1), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(list(0.10, 0.08, 0.06), 0.05), prevalence = c(0.5, 0.3, 0.2), seed = 3 ) identical(df3, df3b) # Two independent factors (2 x 2): columns subgroup1 and subgroup2. # Four cells in column-major order: (1,1), (2,1), (1,2), (2,2). df4 <- simdata_fast( nsim = 100, n = 200, a.time = c(0, 12), a.rate = 200 / 12, e.hazard = list(0.10, 0.08, 0.07, 0.05), prevalence = list(c(0.5, 0.5), c(0.6, 0.4)), seed = 4 ) head(df4) # Two correlated factors via a joint-distribution array (2 x 2) df5 <- simdata_fast( nsim = 100, n = 200, a.time = c(0, 12), a.rate = 200 / 12, e.hazard = 0.08, prevalence = array(c(0.40, 0.10, 0.15, 0.35), dim = c(2, 2)), seed = 5 ) head(df5) # Two correlated endpoints, no switching (the maximal-independence model of # Fleischer et al., 2009, because h12 defaults to h02). In oncology e1 is PFS # and e2 is OS; here the # control has faster intermediate events and faster direct terminal events. dfid <- simdata_fast( nsim = 100, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, h01.median = list(8, 12), h02.median = list(24, 30), seed = 6 ) head(dfid) # With treatment switching: 40 percent of control subjects who reach the # intermediate event switch and then follow a more favorable post-event hazard dfsw <- simdata_fast( nsim = 100, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, h01.median = list(8, 12), h02.median = list(24, 30), switch.prop = list(0.4, 0), h12.switch.median = list(36, 36), seed = 7 ) head(dfsw) # Illness-death model with two subgroups (prevalence 0.4 and 0.6): the # treatment effect on progression is larger in subgroup 2. dfsg <- simdata_fast( nsim = 100, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, h01.median = list(8, list(10, 14)), h02.median = list(24, 30), prevalence = c(0.4, 0.6), seed = 8 ) table(dfsg$group, dfsg$subgroup) / 100 # Three-arm trial (one control and two treatment arms) analyzed as pairwise # contrasts against the shared control. dfk <- simdata_fast( nsim = 100, n = c(120, 120, 120), a.time = c(0, 12), a.rate = 360 / 12, e.median = list(12, 16, 20), seed = 8 ) # Control arm (group 1) versus treatment arm 2, one-sided log-rank at month 24. sub12 <- dfk[dfk$group %in% c(1, 2), ] res12 <- analysis_fast(sub12, control = 1, time.looks = 24, side = 1) head(res12)
Aggregates the per-simulation, per-look output of analysis_fast
into operating characteristics: the rejection rate, the futility-stopping
rate, the distribution of the stopping look, and the expected analysis timing
(events and calendar time at stopping). A group-sequential design is summarized
by applying the supplied per-look boundaries in sequence, and a fixed design is
the single-look (K = 1) degenerate case of the same logic. This function
consumes boundaries computed elsewhere (for example by gsDesign or rpact) and
does not compute or spend alpha itself; it estimates the stopping probabilities
by Monte Carlo over the simulated trials rather than by the analytic numerical
integration used by those packages, so the two agree only up to Monte Carlo
error and converge as the number of simulations grows.
simsummary_fast( data, eff.col = NULL, efficacy = NULL, fut.col = eff.col, futility = NULL, direction = c("lower", "upper"), p.col = NULL, alpha = NULL )simsummary_fast( data, eff.col = NULL, efficacy = NULL, fut.col = eff.col, futility = NULL, direction = c("lower", "upper"), p.col = NULL, alpha = NULL )
data |
A data frame from |
eff.col |
A single character naming the statistic column compared with
|
efficacy |
A numeric vector of efficacy boundaries, one per look, on the
scale of |
fut.col |
A single character naming the statistic column compared with
|
futility |
A numeric vector of futility boundaries, one per look, on the
scale of |
direction |
A single string, either |
p.col |
A single character naming the p-value column compared with
|
alpha |
A numeric vector of per-look nominal significance levels, one per
look. Entries may be |
The boundaries are supplied per look, one value for each distinct value of the
look column in data. Two boundary modes are available and exactly
one must be used.
In the Z mode the efficacy boundary is compared with the column named by
eff.col and the futility boundary with the column named by
fut.col. The two columns may differ, which lets the efficacy and
futility rules live on different scales. For example a beta-spending efficacy
boundary can be applied to the standardized statistic logrank.z while a
futility boundary expressed as a log hazard ratio is applied to
cox.coef. The crossing direction is set by direction: with
"lower" (the natural sign of logrank.z, cox.z, and
cox.coef, where treatment benefit is negative) the efficacy boundary is
crossed when the efficacy statistic is at or below efficacy and the
futility boundary is crossed when the futility statistic is at or above
futility; with "upper" the inequalities are reversed. Either
boundary vector may contain NA at some looks to omit that rule there, so
efficacy-only and futility-only looks are expressed by placing NA in the
other vector. A look with NA on both rules can never stop the trial.
In the p mode the p-value in the column named by p.col is compared with
the per-look nominal level alpha, rejecting when the p-value is at or
below the level. Futility is not used in the p mode, and alpha may
contain NA to omit the efficacy test at a look.
For each simulated trial the looks are examined in order. The trial stops at
the first look whose statistic crosses a boundary. Crossing the efficacy
boundary is a rejection of the null hypothesis; crossing the futility boundary
(Z mode only) is a stop without rejection. When both are crossed at the same
look the efficacy stop takes precedence. A look whose relevant statistic is
NA, or whose boundary is NA, triggers no crossing of that rule
and the trial continues. A trial that reaches the final look without crossing
the efficacy boundary does not reject.
Each look's prob.stop.efficacy is the marginal probability of stopping
for efficacy for the first time at that look, that is the probability of not
stopping at any earlier look and crossing the efficacy boundary at this one.
This is the stage-wise rejection contribution of a group-sequential design; the
cumulative sum cum.reject is the cumulative power up to and including
that look, and the total over all looks is rejection.rate. These are the
same quantities that gsDesign and rpact report, estimated here by simulation.
A simulation whose event target was not reached at a look (reached =
FALSE in analysis_fast output) carries the statistics of its
final data and an NA cutoff at that look. Its statistics are compared
with the boundary like any other, and it is left out of cutoff.mean
(the mean over the simulations with a finite cutoff) but not out of the
count means. In pairwise_fast output the statistics and
counts of an unreached look are NA, so that look never crosses a
boundary and is left out of the count means.
The rejection rate is the type I error under a null data-generating truth and
the power under an alternative truth, but because the function does not know the
truth used to generate data it is reported neutrally as the rejection
rate and its interpretation is left to the user.
When data carries a population column (the long form produced by
analysis_fast with by.subgroup = TRUE), the same
boundaries are applied within each population and the output has one block of
rows per population. Likewise, when data carries an arm column
(the output of pairwise_fast), each experimental arm's
contrast against the control is summarized separately, and the output gains
an arm column with one block of rows per arm (and population). Within
a block, data must have at most one row per simulation and look.
An object of class "simsummary_fast": a data frame with one row
per population and look plus an overall summary row appended after each
population's looks. The columns are arm (only when data has
an arm column), population, look (the look
index, or "overall" on the summary row), optionally look.value,
n.enrolled.mean and n.event.mean (the mean enrolled and event
counts at that look, or at the stopping look on the summary row),
n.dropout.mean and n.pipeline.mean (the mean dropout count and
pipeline count n.enrolled - n.event - n.dropout, when those columns
are present in data), cutoff.mean (the mean calendar time,
likewise), prob.stop.efficacy,
prob.stop.futility, prob.stop.any, and cum.reject. On the
overall row prob.stop.efficacy is the total rejection rate,
prob.stop.futility the total futility rate, prob.stop.any their
sum, and cum.reject again the total rejection rate; its timing columns
are the expected counts and calendar time at the stopping look. The number of
simulations is stored in the attribute nsim (one value per block when
the blocks differ) and the boundary settings, with the look labels and a
copy of the rows of the result, in the attribute boundary.
analysis_fast, simdata_fast,
print.simsummary_fast, plot.simsummary_fast.
df <- simdata_fast( nsim = 200, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(0.05, 0.035), seed = 1 ) res <- analysis_fast(df, control = 1, event.looks = c(60, 105, 150), stat = c("logrank", "coxph"), side = 1) # Efficacy on the standardized log-rank Z, futility on the log hazard ratio, # with a futility-only first look and efficacy-only later looks simsummary_fast(res, eff.col = "logrank.z", efficacy = c(NA, -2.96, -1.97), fut.col = "cox.coef", futility = c(log(1.2), NA, NA), direction = "lower") # p-value boundaries instead simsummary_fast(res, p.col = "logrank.p", alpha = c(0.0006, 0.0151, 0.0245))df <- simdata_fast( nsim = 200, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.hazard = list(0.05, 0.035), seed = 1 ) res <- analysis_fast(df, control = 1, event.looks = c(60, 105, 150), stat = c("logrank", "coxph"), side = 1) # Efficacy on the standardized log-rank Z, futility on the log hazard ratio, # with a futility-only first look and efficacy-only later looks simsummary_fast(res, eff.col = "logrank.z", efficacy = c(NA, -2.96, -1.97), fut.col = "cox.coef", futility = c(log(1.2), NA, NA), direction = "lower") # p-value boundaries instead simsummary_fast(res, p.col = "logrank.p", alpha = c(0.0006, 0.0151, 0.0245))
Computes the log-rank test statistic for comparing survival curves between
two groups. Returns either a one-sided Z-score or a two-sided chi-square
statistic. The C++ backend uses a two-pointer merge scan over pooled sorted
vectors, eliminating the rank() call that dominates the pure-R
implementation. When a strata argument is supplied, the stratified
log-rank test is computed instead, matching
survdiff with a strata() term. A non-default
weight argument selects a weighted log-rank test (Fleming-Harrington,
modestly-weighted, Gehan-Breslow, or Tarone-Ware) for non-proportional
hazards.
survdiff_fast( time, event, group, control, side = 2, presorted = FALSE, strata = NULL, weight = c("logrank", "fh", "mwlrt", "gehan", "tarone-ware"), rho = 0, gamma = 0, t_star = NULL, data = NULL )survdiff_fast( time, event, group, control, side = 2, presorted = FALSE, strata = NULL, weight = c("logrank", "fh", "mwlrt", "gehan", "tarone-ware"), rho = 0, gamma = 0, t_star = NULL, data = NULL )
time |
A numeric vector of follow-up times for all subjects.
Alternatively, a formula |
event |
An integer or numeric vector of event indicators
(1 = event, 0 = censored), aligned with |
group |
A vector of group labels aligned with |
control |
A scalar value indicating which level of |
side |
An integer, either 1 or 2 (default 2). If |
presorted |
A logical value. If |
strata |
An optional vector of stratum labels aligned with |
weight |
A character string naming the weight scheme. |
rho |
A numeric Fleming-Harrington first parameter, used only when
|
gamma |
A numeric Fleming-Harrington second parameter, used only when
|
t_star |
A single non-negative numeric value, the timepoint of the
modestly-weighted log-rank test. Required only when |
data |
A data frame in which the variables of a formula given as
|
The log-rank statistic is computed as:
Z = (O_1 - E_1) / sqrt(V_1)
where O_1 is the observed number of events in the treatment group, E_1 is the expected number under the null hypothesis of equal survival, and V_1 is the hypergeometric variance. Tied event times are handled correctly: all subjects sharing the same event time form a tied block, and the block is processed atomically in the two-pointer merge.
When strata is NULL (default), the ordinary two-group log-rank
test is computed by the C++ core logrank_core, which walks the pooled
sorted data with a single two-pointer scan, maintaining running at-risk
counts per group. No rank vector is constructed, so the dominant O(n log n)
cost of rank() in the pure-R version is removed.
When strata is supplied, the stratified log-rank test is computed by
the C++ core stratified_logrank_core. The contributions O_1, E_1, and
V_1 are accumulated within each stratum and then summed across strata, so
the overall statistic is Z = (sum O_1 - sum E_1) / sqrt(sum V_1). A stratum
that contains only one group contributes zero to all three totals, the same
convention used by survdiff. This is the standard
stratified log-rank test, equivalent to a log-rank test that conditions on
the stratum at each event time.
When a non-default weight is requested, the weighted log-rank test is
computed by the C++ core weighted_logrank_core. The statistic is
Z = U / sqrt(V) with U = sum w_j (O_1j - E_1j) and
V = sum w_j^2 n_0j n_1j O_j (n_j - O_j) / (n_j^2 (n_j - 1)), where the weight
w_j is one of the schemes named by weight. The Fleming-Harrington and
modestly-weighted schemes use the left-continuous pooled Kaplan-Meier
estimate S(t-), initialized at 1 and updated after each event time, so the
first event always has S(t-) = 1. The modestly-weighted scheme determines
its weight cap in a first pass over the event times before accumulating U
and V in a second pass; the other schemes accumulate in a single pass. When
both weight and strata are supplied, the stratified weighted
log-rank test is computed by the C++ core
stratified_weighted_logrank_core: each stratum is an independent
weighted log-rank test whose weights come from that stratum's own pooled
Kaplan-Meier estimate, and the per-stratum U and V are summed before
standardizing once as Z = sum U / sqrt(sum V).
When presorted = TRUE, the input vectors are assumed to be sorted and
the internal order() call is skipped. In the unstratified case the
assumed order is ascending time. In the stratified case the assumed
order is by stratum first and by ascending time within each stratum,
so that rows of the same stratum are contiguous. When presorted =
FALSE (default), sorting is handled internally. In simulation loops where
the data are generated in the required order, setting presorted =
TRUE avoids one O(n log n) pass.
The returned object has class "survdiff_fast" and is a length-one
numeric value (Z-score or chi-square) with the underlying counts O_0, E_0,
O_1, E_1, V_1, the requested side, and the total sample size stored
as attributes. When strata is supplied, the number of strata is also
stored as the attribute strata. A print() method formats the
result similarly to print(survival::survdiff(...)), displaying
observed and expected event counts for both the control and treatment
groups.
An object of class "survdiff_fast", which is a length-one
numeric value with attributes O0, E0, O1, E1,
V1, side, and n, plus strata (the number of
strata) when strata is supplied, or weight (the scheme name)
when a non-default weight is used. The numeric value is the Z-score
when side = 1, or the chi-square statistic when side = 2.
For a weighted test the value is U / sqrt(V) (or its square), V1
holds the weighted variance V, O0 and O1 hold the raw
observed event counts, and E0 and E1 are NA because a
single unweighted expected count is not defined for a weighted test.
Returns NA_real_ (still with class "survdiff_fast") when the
variance is zero (e.g., no event time at which both groups are at risk).
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Gehan, E. A. (1965). A generalized Wilcoxon test for comparing arbitrarily singly-censored samples. Biometrika, 52(1-2), 203-224.
Mantel, N. (1966). Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemotherapy Reports, 50(3), 163-170.
Tarone, R. E., & Ware, J. (1977). On distribution-free tests for equality of survival distributions. Biometrika, 64(1), 156-160.
Fleming, T. R., & Harrington, D. P. (1991). Counting Processes and Survival Analysis. New York: John Wiley & Sons.
Magirr, D., & Burman, C.-F. (2019). Modestly weighted logrank tests. Statistics in Medicine, 38(20), 3782-3790.
survdiff for the standard implementation.
print.survdiff_fast for the print method.
library(survival) # Two-sided test: compare with survdiff fit <- survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 2) fit chisq_ref <- survdiff(Surv(futime, fustat) ~ rx, data = ovarian)$chisq cat("survdiff_fast chi-square:", as.numeric(fit), "\n") cat("survdiff chi-square:", chisq_ref, "\n") # One-sided test (Z-score) survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 1) # presorted = TRUE: sort once outside, reuse inside a loop ord <- order(ovarian$futime) survdiff_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord], control = 2, side = 2, presorted = TRUE) # Stratified log-rank test: compare with survdiff + strata() fit_str <- survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 2, strata = ovarian$resid.ds) chisq_str <- survdiff(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian)$chisq cat("stratified survdiff_fast:", as.numeric(fit_str), "\n") cat("stratified survdiff :", chisq_str, "\n") # Weighted log-rank tests for non-proportional hazards # Fleming-Harrington G(0, 1), emphasizing late differences survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 1, weight = "fh", rho = 0, gamma = 1) # Modestly-weighted log-rank test with t_star = 365 survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 1, weight = "mwlrt", t_star = 365) # Stratified weighted log-rank test: Fleming-Harrington G(0,1) within strata survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 1, weight = "fh", rho = 0, gamma = 1, strata = ovarian$resid.ds) # Formula interface (the first level of rx, 1, is the control group) survdiff_fast(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian) # Speed comparison against survdiff() if (requireNamespace("microbenchmark", quietly = TRUE)) { microbenchmark::microbenchmark( survdiff_fast = survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 2), survdiff = survdiff(Surv(futime, fustat) ~ rx, data = ovarian), times = 1000 ) }library(survival) # Two-sided test: compare with survdiff fit <- survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 2) fit chisq_ref <- survdiff(Surv(futime, fustat) ~ rx, data = ovarian)$chisq cat("survdiff_fast chi-square:", as.numeric(fit), "\n") cat("survdiff chi-square:", chisq_ref, "\n") # One-sided test (Z-score) survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 1) # presorted = TRUE: sort once outside, reuse inside a loop ord <- order(ovarian$futime) survdiff_fast(ovarian$futime[ord], ovarian$fustat[ord], ovarian$rx[ord], control = 2, side = 2, presorted = TRUE) # Stratified log-rank test: compare with survdiff + strata() fit_str <- survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 2, strata = ovarian$resid.ds) chisq_str <- survdiff(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian)$chisq cat("stratified survdiff_fast:", as.numeric(fit_str), "\n") cat("stratified survdiff :", chisq_str, "\n") # Weighted log-rank tests for non-proportional hazards # Fleming-Harrington G(0, 1), emphasizing late differences survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 1, weight = "fh", rho = 0, gamma = 1) # Modestly-weighted log-rank test with t_star = 365 survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 1, weight = "mwlrt", t_star = 365) # Stratified weighted log-rank test: Fleming-Harrington G(0,1) within strata survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 1, weight = "fh", rho = 0, gamma = 1, strata = ovarian$resid.ds) # Formula interface (the first level of rx, 1, is the control group) survdiff_fast(Surv(futime, fustat) ~ rx + strata(resid.ds), data = ovarian) # Speed comparison against survdiff() if (requireNamespace("microbenchmark", quietly = TRUE)) { microbenchmark::microbenchmark( survdiff_fast = survdiff_fast(ovarian$futime, ovarian$fustat, ovarian$rx, 2, side = 2), survdiff = survdiff(Surv(futime, fustat) ~ rx, data = ovarian), times = 1000 ) }
Computes the Kaplan-Meier survival probability at a specified time point, together with a standard error and confidence interval based on Greenwood's variance formula. The C++ backend performs binary search for the evaluation cutoff and accumulates the Kaplan-Meier product and Greenwood sum in a single scan over event positions only, without constructing intermediate vectors.
survfit_fast( t_sorted, e_sorted, t_eval, conf.level = 0.95, conf.type = "log", presorted = TRUE, data = NULL )survfit_fast( t_sorted, e_sorted, t_eval, conf.level = 0.95, conf.type = "log", presorted = TRUE, data = NULL )
t_sorted |
A numeric vector of event or censoring times. Must be
sorted in ascending order when |
e_sorted |
An integer or numeric vector of event indicators
(1 = event, 0 = censored), aligned with |
t_eval |
A single numeric value specifying the time point at which the survival probability is evaluated. |
conf.level |
A single numeric value in (0, 1) specifying the confidence level. Defaults to 0.95. |
conf.type |
A character string specifying the confidence interval type.
Must be one of |
presorted |
A logical value. If |
data |
A data frame in which the variables of a formula given as
|
The Kaplan-Meier estimate at time t_eval is defined as the
product-limit estimator evaluated at the largest observed event time less
than or equal to t_eval. If t_eval is smaller than the first
observed event time, S(t) = 1 and the standard error is zero.
The standard error is estimated by Greenwood's formula:
SE[S(t)] = S(t) * sqrt(sum_{t_i <= t, d_i > 0} d_i / (n_i * (n_i - d_i)))
where d_i is the number of events and n_i is the number at risk at time t_i.
The output field std.err is SE[S(t)], the standard error on the
survival scale, which is the std.err reported by summary()
of a survfit object. The std.err component
of the survfit object itself is on the cumulative-hazard scale,
SE[S(t)] / S(t), whatever conf.type is.
When S(t_eval) = 0 (all subjects have experienced the event by
t_eval), the standard error is zero and the confidence interval
collapses to [0, 0], consistent with survfit.
When presorted = TRUE (default), t_sorted and e_sorted
are expected to be sorted in ascending order of time; the order is checked
(a linear-time test) and an error is raised for unsorted input. When
presorted = FALSE, the vectors are sorted internally before
computation. When the standard error is zero (for example S(t) = 1
before the first event), the confidence interval collapses to the estimate.
Three confidence interval types are supported via conf.type:
"plain": Linear interval on the survival scale,
S(t) +/- z * SE. The bounds are clipped to [0, 1].
"log": Interval on the log scale (default in
survfit),
S(t) * exp(+/- z * SE / S(t)). The upper bound is capped at 1, as in
survfit.
"log-log": Interval on the complementary log-log scale,
S(t)^exp(+/- z * SE / (S(t) * log(S(t)))).
The returned object has class "survfit_fast" and is a named numeric
vector of length 4 with the evaluation time t_eval, the confidence
level conf.level, and the confidence interval type conf.type
stored as attributes. A print() method formats the result similarly
to print(summary(survival::survfit(...))).
An object of class "survfit_fast", which is a named numeric
vector of length 4 with elements surv, std.err,
lower, and upper, representing the Kaplan-Meier survival
estimate, the Greenwood standard error SE[S(t)], and the lower and upper
confidence limits at t_eval. The evaluation time, confidence
level, and confidence interval type are stored as attributes
t_eval, conf.level, and conf.type.
Returns a vector of NA_real_ values (still with class
"survfit_fast") when n is zero.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457-481.
survfit for the standard Kaplan-Meier estimator.
print.survfit_fast for the print method.
set.seed(42) t_raw <- rexp(100, rate = 1 / 10) e_raw <- rbinom(100, 1, 0.7) # presorted = TRUE (default): sort once outside, reuse inside a loop ord <- order(t_raw) t_s <- t_raw[ord] e_s <- e_raw[ord] survfit_fast(t_s, e_s, t_eval = 10, conf.type = "plain") survfit_fast(t_s, e_s, t_eval = 10, conf.type = "log") survfit_fast(t_s, e_s, t_eval = 10, conf.type = "log-log") # presorted = FALSE: sort internally, convenient for one-off calls survfit_fast(t_raw, e_raw, t_eval = 10, presorted = FALSE) # Validation against survival::survfit library(survival) fit <- survfit(Surv(t_raw, e_raw) ~ 1, conf.type = "plain") summary(fit, times = 10) # Formula interface (the data are sorted internally) d <- data.frame(time = t_raw, status = e_raw) survfit_fast(Surv(time, status) ~ 1, data = d, t_eval = 10)set.seed(42) t_raw <- rexp(100, rate = 1 / 10) e_raw <- rbinom(100, 1, 0.7) # presorted = TRUE (default): sort once outside, reuse inside a loop ord <- order(t_raw) t_s <- t_raw[ord] e_s <- e_raw[ord] survfit_fast(t_s, e_s, t_eval = 10, conf.type = "plain") survfit_fast(t_s, e_s, t_eval = 10, conf.type = "log") survfit_fast(t_s, e_s, t_eval = 10, conf.type = "log-log") # presorted = FALSE: sort internally, convenient for one-off calls survfit_fast(t_raw, e_raw, t_eval = 10, presorted = FALSE) # Validation against survival::survfit library(survival) fit <- survfit(Surv(t_raw, e_raw) ~ 1, conf.type = "plain") summary(fit, times = 10) # Formula interface (the data are sorted internally) d <- data.frame(time = t_raw, status = e_raw) survfit_fast(Surv(time, status) ~ 1, data = d, t_eval = 10)
Applies a treatment-switching rule to simulated trial data, such as the
output of simdata_fast, by changing only the outcomes that
occur after each subject's switch. Subjects of the selected groups switch
either at their intermediate event (for example disease progression), at
an opening time after an interim analysis (crossover at a milestone), or at
the later of the two, and the remaining time to the terminal event after the
switch is either multiplied by an acceleration factor or redrawn from a new
(piecewise) exponential hazard. All simulated trials are processed at once
with vectorized operations. Because the outcomes observed before the switch
are kept exactly, an analysis at or before the opening time gives the same
result before and after switching, so a crossover that depends on an interim
result can be simulated by analyzing the interim, deciding per simulated
trial, and switching afterward.
switch_fast( data, group, prob = 1, when = c("intermediate", "cutoff", "later"), cutoff = NULL, delay = 0, sims = NULL, aft.factor = NULL, hazard = NULL, median = NULL, time = NULL, seed = NULL, stream = NULL )switch_fast( data, group, prob = 1, when = c("intermediate", "cutoff", "later"), cutoff = NULL, delay = 0, sims = NULL, aft.factor = NULL, hazard = NULL, median = NULL, time = NULL, seed = NULL, stream = NULL )
data |
A data frame from |
group |
A vector of group labels whose subjects may switch, for example
the control group |
prob |
A single probability in [0, 1] that an eligible subject switches. Defaults to 1. |
when |
A character string, one of |
cutoff |
Per-simulation calendar times at which switching opens, used
with |
delay |
A single non-negative number added to |
sims |
An optional logical vector with one element per simulated trial
(in the order of the sorted distinct values of |
aft.factor |
A single positive acceleration factor applied to the
remaining time after the switch. Exactly one of |
hazard |
Post-switch hazard, a single value or a vector for a piecewise
exponential distribution with breakpoints |
median |
A single post-switch median; an alternative to a single
|
time |
Breakpoints for a piecewise |
seed |
Optional integer seed for the |
stream |
Optional non-negative whole number selecting a |
The function works on the latent columns of simdata_fast. For
a single-endpoint simulation these are surv_time and
dropout_time; for an illness-death simulation they are
e1_surv_time (time to the first event, such as progression-free
survival), e2_surv_time (time to the terminal event, such as overall
survival), dropout_time, and intermediate. The observed
columns (tte, event, and calendar_time, or their
e1_ and e2_ counterparts) are recomputed from the modified
latent times, and the columns switched and switch_time (time
from accrual to the switch) record the switches. The columns sim,
group, accrual_time, and the subgroup columns are unchanged.
The switch time s, measured from accrual, is
when = "intermediate": the intermediate event time
e1_surv_time of a subject whose intermediate event occurred
(intermediate == 1); illness-death data only;
when = "cutoff": the opening time cutoff + delay of the
subject's simulated trial minus the accrual time (zero for a subject
accrued after the opening);
when = "later": the later of the intermediate event time and
the opening time, for a subject whose intermediate event occurred;
illness-death data only.
A subject of the selected group switches with probability
prob when the switch occurs before the terminal event and before
dropout (s < e2_surv_time, or s < surv_time for
single-endpoint data, and s < dropout_time) and the subject has not
already switched (switched == 1, for example at progression through
the switch.prop argument of simdata_fast). With when = "cutoff" or "later", the switch
never occurs before the opening time, so every outcome observed by calendar
time cutoff + delay is left unchanged; subjects accrued after the
opening are also eligible, as in a protocol that opens crossover from that
time on.
The effect of switching on the terminal time T > s is
aft.factor = f: T_new = s + f * (T - s), the causal
accelerated failure time model of the rank-preserving structural failure
time (RPSFT) method, under which f > 1 prolongs the remaining
survival;
hazard (or median): T_new = s + R, where the
remaining time R is drawn from the (piecewise) exponential
distribution with hazard hazard and breakpoints time
measured from the switch.
In illness-death data with when = "cutoff", a switch can occur before
the intermediate event. The acceleration factor is then applied in the same
way to the first-event time when it lies after the switch, which keeps the
first event no later than the terminal event; a redrawn hazard is not
available in this case because the latent time to the intermediate event
after the terminal event is not stored.
sims restricts the switching to some simulated trials, for example
those in which the interim analysis crossed an efficacy boundary. It is
available with when = "cutoff" or "later" only, so that the
decision cannot change outcomes observed before it is made.
The random numbers (one uniform per row when prob < 1, and one
exponential per row when hazard or median is used) are drawn
with dqrng, in row order, so the result is reproducible from
seed and stream, or from the state left by a seeded
simdata_fast call.
The data frame data with the latent and observed columns of
the switched subjects updated, and the columns switched (0 or 1)
and switch_time (time from accrual to the switch, NA for
subjects who did not switch) added or updated.
Robins, J. M., & Tsiatis, A. A. (1991). Correcting for non-compliance in randomized trials using rank preserving structural failure time models. Communications in Statistics - Theory and Methods, 20(8), 2609-2631.
simdata_fast, cutoff_fast,
analysis_fast.
# Crossover at an interim analysis: control subjects still on study at the # 150th event switch to the experimental treatment, which multiplies their # remaining survival time by 1.5. df <- simdata_fast(nsim = 50, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.median = list(12, 18), d.hazard = 0.01, seed = 1) ia <- cutoff_fast(df, event.looks = 150) dfs <- switch_fast(df, group = 1, when = "cutoff", cutoff = ia, aft.factor = 1.5) mean(dfs$switched[dfs$group == 1]) # The interim analysis is unchanged; the final analysis is diluted. r0 <- analysis_fast(df, control = 1, cutoff.looks = ia) r1 <- analysis_fast(dfs, control = 1, cutoff.looks = ia) all.equal(r0$logrank.z, r1$logrank.z) fa0 <- analysis_fast(df, control = 1, event.looks = 220) fa1 <- analysis_fast(dfs, control = 1, event.looks = 220) c(mean(fa0$logrank.z), mean(fa1$logrank.z)) # Switching at progression in an illness-death simulation: 50 percent of # the control subjects who progress switch, with a post-progression median # of 15 instead of 10. id <- simdata_fast(nsim = 50, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, h01.median = list(6, 9), h02.median = list(20, 24), h12.median = list(10, 10), seed = 2) ids <- switch_fast(id, group = 1, prob = 0.5, when = "intermediate", median = 15) table(ids$group, ids$switched)# Crossover at an interim analysis: control subjects still on study at the # 150th event switch to the experimental treatment, which multiplies their # remaining survival time by 1.5. df <- simdata_fast(nsim = 50, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, e.median = list(12, 18), d.hazard = 0.01, seed = 1) ia <- cutoff_fast(df, event.looks = 150) dfs <- switch_fast(df, group = 1, when = "cutoff", cutoff = ia, aft.factor = 1.5) mean(dfs$switched[dfs$group == 1]) # The interim analysis is unchanged; the final analysis is diluted. r0 <- analysis_fast(df, control = 1, cutoff.looks = ia) r1 <- analysis_fast(dfs, control = 1, cutoff.looks = ia) all.equal(r0$logrank.z, r1$logrank.z) fa0 <- analysis_fast(df, control = 1, event.looks = 220) fa1 <- analysis_fast(dfs, control = 1, event.looks = 220) c(mean(fa0$logrank.z), mean(fa1$logrank.z)) # Switching at progression in an illness-death simulation: 50 percent of # the control subjects who progress switch, with a post-progression median # of 15 instead of 10. id <- simdata_fast(nsim = 50, n = c(150, 150), a.time = c(0, 12), a.rate = 300 / 12, h01.median = list(6, 9), h02.median = list(20, 24), h12.median = list(10, 10), seed = 2) ids <- switch_fast(id, group = 1, prob = 0.5, when = "intermediate", median = 15) table(ids$group, ids$switched)
Computes the weighted Kaplan-Meier (WKM) statistic of Pepe and Fleming for comparing two survival curves. The statistic is the weighted integral of the difference between the Kaplan-Meier estimates of the treatment and control groups over the observed range, standardized to a Wald z statistic. Unlike weighted log-rank tests, this test targets the integrated difference in survival and is sensitive to differences even when the hazard functions cross. The single scan over the sorted data is performed in C++ for use inside simulation loops.
wkm_fast( time, event, group, control = NULL, side = 2, conf.level = 0.95, weight = c("PF", "sqrtPF", "constant"), presorted = FALSE, data = NULL )wkm_fast( time, event, group, control = NULL, side = 2, conf.level = 0.95, weight = c("PF", "sqrtPF", "constant"), presorted = FALSE, data = NULL )
time |
Numeric vector of event or censoring times.
Alternatively, a formula |
event |
Integer vector, 1 for an event and 0 for censoring. |
group |
Grouping vector with exactly two distinct levels. |
control |
The level of |
side |
Either 2 for a two-sided test or 1 for a one-sided test of treatment superiority (weighted difference greater than 0). |
conf.level |
Confidence level for the interval of the weighted difference. |
weight |
Weight function, one of |
presorted |
Logical; set to |
data |
A data frame in which the variables of a formula given as
|
The default weight is the Pepe-Fleming combined censoring weight
w(t) = n G1(t) G2(t) / (n1 G1(t) + n2 G2(t)), where G1 and G2 are the
Kaplan-Meier estimates of the censoring survival functions of the two
groups. This weight gives less weight to regions with heavy censoring and
stabilizes the variance in the
tail. The choice weight = "sqrtPF" uses its square root, and
weight = "constant" uses a weight of 1, in which case the numerator
reduces to the difference in restricted mean survival time over the observed
range. With weight = "PF" the result reproduces
nphsim::wkm.Stat for data without tied times.
Because the integral extends to the largest observed time of the pooled data, the curve of a group whose follow-up ends earlier with a censored time is carried forward flat beyond that time. The Pepe-Fleming weight is zero there, but the constant and square-root weights are not, so with these weights the statistic can be biased, even under the null hypothesis, when follow-up differs between the groups. The Pepe-Fleming weight is recommended in that case.
The integral runs up to the largest observed time, so it is not defined
when a time is infinite (for example a subject without a finite event or
dropout time in simulated data), and all the statistics are then NA.
The weighted difference is computed as treatment minus control, so a positive value and a positive z indicate longer survival under treatment.
A named numeric vector of class "wkm_fast" with the weighted
integrated difference, its standard error and confidence limits, and the
Wald statistics for the test.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Pepe, M. S., & Fleming, T. R. (1989). Weighted Kaplan-Meier statistics: a class of distance tests for censored survival data. Biometrics, 45(2), 497-507.
Pepe, M. S., & Fleming, T. R. (1991). Weighted Kaplan-Meier statistics: large sample and optimality considerations. Journal of the Royal Statistical Society. Series B (Methodological), 53(2), 341-352.
set.seed(1) n <- 300 g <- rep(0:1, each = n / 2) tt <- c(rexp(n / 2, log(2) / 12), rexp(n / 2, log(2) / 16)) cc <- rexp(n, rate = 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) wkm_fast(time, event, group = g, control = 0) # Formula interface d <- data.frame(time, event, g) wkm_fast(Surv(time, event) ~ g, data = d)set.seed(1) n <- 300 g <- rep(0:1, each = n / 2) tt <- c(rexp(n / 2, log(2) / 12), rexp(n / 2, log(2) / 16)) cc <- rexp(n, rate = 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) wkm_fast(time, event, group = g, control = 0) # Formula interface d <- data.frame(time, event, g) wkm_fast(Surv(time, event) ~ g, data = d)
Estimates the window mean survival time (WMST) for a single group, or the
difference in WMST between a treatment group and a control group, together
with standard errors, confidence intervals and a Wald test for the
difference. The WMST is the area under the Kaplan-Meier curve between the
window limits tau1 and tau2. With tau1 = 0 it is the
restricted mean survival time at tau2, and the result then matches
survRM2::rmst2. The single scan over the sorted data is performed in
C++ for use inside simulation loops.
wmst_fast( time, event, group = NULL, control = NULL, tau1 = 0, tau2 = NULL, side = 2, conf.level = 0.95, presorted = FALSE, data = NULL )wmst_fast( time, event, group = NULL, control = NULL, tau1 = 0, tau2 = NULL, side = 2, conf.level = 0.95, presorted = FALSE, data = NULL )
time |
Numeric vector of event or censoring times.
Alternatively, a formula |
event |
Integer vector, 1 for an event and 0 for censoring. |
group |
Optional grouping vector with exactly two distinct levels for a two-group comparison. If omitted, a single-group WMST is returned. |
control |
The level of |
tau1 |
Lower limit of the window. Defaults to 0. |
tau2 |
Upper limit of the window. If omitted, the largest time common
to both groups is used (the minimum over groups of the maximum observed
time), or the maximum observed time for a single group. A supplied value
larger than that gives a warning, because the Kaplan-Meier curve is not
estimated beyond it. The window must be finite: when that default is
infinite (for example when each group has a subject without a finite
event or dropout time in simulated data), |
side |
Either 2 for a two-sided test or 1 for a one-sided test of treatment superiority (difference greater than 0). |
conf.level |
Confidence level for the intervals. |
presorted |
Logical; set to |
data |
A data frame in which the variables of a formula given as
|
The variance of the WMST estimate uses the counting process (Greenwood) form
in which each event time at or before tau2 contributes the squared
remaining window area. Event times at or before tau1 contribute the
full window area, since the hazard before tau1 scales the survival
curve throughout the window.
The difference is computed as treatment minus control, so a positive difference indicates a longer mean survival time within the window under treatment.
A named numeric vector of class "wmst_fast". For a single
group the elements are the WMST, its standard error and confidence limits.
For two groups the elements are the control and treatment WMST values,
their difference, the standard errors, the confidence limits and the Wald
statistics for the difference.
The object also has the class "fastsurv", whose
as.data.frame(), coef(), and confint() methods
are described in as.data.frame.fastsurv.
Paukner, M., & Chappell, R. (2021). Window mean survival time. Statistics in Medicine, 40(25), 5521-5533.
set.seed(1) n <- 200 g <- rep(0:1, each = n / 2) tt <- c(rexp(n / 2, 0.1), rexp(n / 2, 0.07)) cc <- rexp(n, rate = 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) wmst_fast(time, event, group = g, control = 0, tau1 = 2, tau2 = 12) # Formula interface d <- data.frame(time, event, g) wmst_fast(Surv(time, event) ~ g, data = d, tau1 = 2, tau2 = 12)set.seed(1) n <- 200 g <- rep(0:1, each = n / 2) tt <- c(rexp(n / 2, 0.1), rexp(n / 2, 0.07)) cc <- rexp(n, rate = 0.02) time <- pmin(tt, cc) event <- as.integer(tt <= cc) wmst_fast(time, event, group = g, control = 0, tau1 = 2, tau2 = 12) # Formula interface d <- data.frame(time, event, g) wmst_fast(Surv(time, event) ~ g, data = d, tau1 = 2, tau2 = 12)