Package {TemporalHazard}

TemporalHazard: Temporal Parametric Hazard Modeling

Description

Native R implementation of the multiphase parametric hazard model of Blackstone, Naftel, and Turner (1986), with a focus on behavioral parity, transparent numerics, and reproducible validation against the original ‘C’/‘SAS’ HAZARD program. The package fits time-varying hazards as an additive sum of parametric phases (early, constant, and late risk streams), each carrying its own covariate effects.

The multiphase model

The total cumulative hazard decomposes additively across J phases:

H(t \mid \mathbf{x}) = \sum_{j=1}^{J} \mu_j(\mathbf{x}) \, \Phi_j(t)

where \mu_j(\mathbf{x}) = \exp(\alpha_j + \mathbf{x}_j^\top \beta_j) is the phase-specific log-linear scale (intercept plus covariate effects) and \Phi_j(t) is the temporal shape contributed by phase j, with its own parameters depending on the phase type (see the phase vocabulary below). The instantaneous hazard and survival follow directly:

h(t \mid \mathbf{x}) = \sum_{j=1}^{J} \mu_j(\mathbf{x}) \, \varphi_j(t), \qquad S(t \mid \mathbf{x}) = \exp\!\bigl(-H(t \mid \mathbf{x})\bigr)

with \varphi_j = d\Phi_j/dt. Each phase is specified with hzr_phase() and the model is fit by maximum likelihood with hazard().

Phase vocabulary

Here G(t) is the generalized temporal decomposition CDF computed by hzr_decompos(), and G_3(t) is the unbounded late-phase intensity from hzr_decompos_g3(). See vignette("mf-mathematical-foundations") for the full derivation.

SAS/C HAZARD bridge

The classic three-phase HAZARD model maps directly onto hzr_phase() calls:

TemporalHazard generalizes the fixed three-phase structure to N phases of any type. hzr_argument_mapping() gives the full parameter translation table between the SAS/C parameterization and the R arguments.

Main entry points

Model fitting

hazard() for building and fitting single- or multiphase models; hzr_phase() for specifying one phase; hzr_stepwise() for forward, backward, or bidirectional covariate selection.

Prediction

predict() on a fitted hazard object for survival, cumulative hazard, and per-phase decomposed hazard; summary() for coefficient tables with Wald inference.

Parametric family

hzr_decompos() for the early-phase (G1) decomposition G(t), g(t), h(t); hzr_decompos_g3() for the late-phase (G3) intensity; hzr_phase_cumhaz() and hzr_phase_hazard() for per-phase \Phi(t) and \varphi(t).

Diagnostics

hzr_kaplan(), hzr_nelson() for nonparametric references; hzr_gof() for goodness of fit; hzr_calibrate(), hzr_deciles() for calibration; hzr_bootstrap() for resampling CIs; hzr_competing_risks() for cumulative incidence.

Vignettes

vignette("getting-started")

First fit, end to end.

vignette("mf-mathematical-foundations")

The decomposition family and multiphase model, with derivations.

vignette("fitting-hazard-models")

Single-phase through multiphase fitting.

vignette("prediction-visualization")

Prediction types and decomposed-hazard plots.

vignette("inference-diagnostics")

Bootstrap CIs and diagnostics.

vignette("sas-to-r-migration")

Translating SAS/C HAZARD code.

Author(s)

Maintainer: John Ehrlinger john.ehrlinger@gmail.com [copyright holder]

Authors:

• John Ehrlinger john.ehrlinger@gmail.com [copyright holder]

Other contributors:

• R Core Team (seeding and inversion in .hzr_mt_normal()) [copyright holder]

• Makoto Matsumoto (Mersenne Twister in .hzr_mt_normal()) [copyright holder]

• Takuji Nishimura (Mersenne Twister in .hzr_mt_normal()) [copyright holder]

References

Blackstone EH, Naftel DC, Turner ME Jr. The decomposition of time-varying hazard into phases, each incorporating a separate stream of concomitant information. J Am Stat Assoc. 1986;81(395):615–624. doi:10.1080/01621459.1986.10478314

Rajeswaran J, Blackstone EH, Ehrlinger J, Li L, Ishwaran H, Parides MK. Probability of atrial fibrillation after ablation: Using a parametric nonlinear temporal decomposition mixed effects model. Stat Methods Med Res. 2018;27(1):126–141. doi:10.1177/0962280215623583

See Also

Useful links:

• https://ehrlinger.github.io/TemporalHazard/

• https://github.com/ehrlinger/TemporalHazard

• Report bugs at https://github.com/ehrlinger/TemporalHazard/issues

Coerce x to a validated numeric design matrix

Description

Accepts data.frame or numeric matrix input; validates dimensions and finiteness. Called by hazard() to normalise the x argument before any downstream use.

Usage

.hzr_as_design_matrix(x, n = NULL)

Arguments

Value

A numeric matrix with column names preserved (or added if absent).

Warn when a phase has effectively left the model

Description

Warns rather than stops: the fit is arithmetically fine and the other phases' estimates are usable. It is the unidentified parameters that must not be read as estimates, and which ones those are differs by mode, so the message says which.

Usage

.hzr_check_phase_identifiability(
  theta,
  time,
  phases,
  covariate_counts,
  x_list,
  tol = 1e-08,
  other_times = NULL
)

Arguments

Details

Three conditions, not two. absent and saturated are per-phase. The third is a property of the times: when their own relative range falls below tol they separate no phase from any other, and the per-phase messages would name a cause that did not occur (#211). It is asked of the times rather than inferred from the phases, because a saturated or absent phase is flat across times that are perfectly well spread.

absent is measured on share, which does not depend on how the times are spread, so it is always reported. The flatness verdict is withheld when other_times vary, since the measures here cannot see those.

Value

The share data.frame, invisibly; called for the warning.

Check the SAS objective's data preconditions at entry

Description

Both conditions objective = "sas" imposes (no left-censored rows, and a positive width on every interval row) are pure functions of the data, so they hold or fail identically at every start. Evaluated inside the objective they reach the user through .hzr_optim_multiphase()'s per-start tryCatch, which frames them as "produced no usable fit from N starts" and invites raising n_starts, a remedy that cannot work. hazard() calls this once, before any optimization, so a data defect is reported as one.

Usage

.hzr_check_sas_data(
  status,
  time,
  time_lower,
  time_upper,
  objective,
  row_ids = NULL
)

Arguments

Details

This does not replace the guards inside the objective and the gradient. Those remain because the gradient is reachable without hazard() (the score test calls it directly), so entry validation is not guaranteed to have run. See .hzr_check_sas_status().

Bounds are normalised here exactly as .hzr_logl_multiphase() normalises them, so the check cannot disagree with the objective about which rows are offenders. Indices are reported against the data, not against the interval subset the inner guard sees.

Value

NULL, invisibly; called for its side effect.

Reject row types the SAS objective has no counterpart for

Description

⁠PROC HAZARD⁠ has no left-censoring statement, so no SAS run corresponds to a fit containing left-censored rows. That is a data defect rather than parameter infeasibility, so it stops rather than returning -Inf.

Usage

.hzr_check_sas_status(status, objective)

Arguments

Details

Called from BOTH .hzr_logl_multiphase() and .hzr_gradient_multiphase(). Guarding only the objective would leave the gradient computing happily for data the objective refuses. The gradient is reachable on its own, for instance from .hzr_score_test(), so the objective's refusal is not guaranteed to come first.

Value

NULL, invisibly; called for its side effect.

Apply the Conservation of Events adjustment to one phase's log_mu

Description

Given the current theta vector, analytically solve the fixmu phase's log_mu so that total predicted events = total observed events.

Usage

.hzr_conserve_events(
  theta,
  fixmu_phase,
  fixmu_pos,
  time,
  status,
  phases,
  covariate_counts,
  x_list,
  total_events,
  weights = NULL,
  time_lower = NULL,
  details = FALSE
)

Arguments

Details

This is called BEFORE each likelihood evaluation inside the optimizer, matching the C HAZARD CONSRV entry point.

Value

Updated theta vector with fixmu phase's log_mu adjusted, or the list described under details.

Log-scale terms shared by the m > 0 branches of hzr_decompos()

Description

Cases 1 (nu > 0) and 3 (nu < 0) evaluate the same intermediate quantities; only the sign of rho and the final assembly of G differ. Forming them directly overflows well inside the range an optimizer can reach: 2^m is Inf from m = 1024, and bt^(-1/nu) overflows sooner still (from m ~ 750 at t/t_half = 0.5, m*nu = 3). btnu then becomes Inf and btnu^(-1/m) collapses silently to 0.

Usage

.hzr_decompos_g1_logs(time, t_half, nu, m)

Arguments

Details

The m in rho's (2^m - 1)/m factor cancels the explicit m multiplier, so exactly

m \, b(t)^{-1/\nu} = (t_{1/2}/t)^{1/\nu} (2^m - 1),

and log(btnu) follows from a softplus of the log of that product. Every intermediate below stays finite for any m > 0 a double can hold.

Value

List with log_S (= -\log(\mathrm{btnu})/m), log_g, a, the softplus argument, so that log(btnu) = log1pexp(a), and the terms log_g is built from (e, log_bt, log_rho), and log_btnu, log1pexp(a) itself, so a caller need not form it again.

Vcov submatrix for the delta-method sandwich

Description

Restricts the sandwich to the parameters the prediction uses. A parameter held fixed (fixed, from the fit's fixed_mask; e.g. fixed = "shapes") with an NA row is dropped as known. (The CoE-conserved log_mu is one in an object saved before #652 whose full-information recompute failed; the caller warns. A fit made since has no vcov in that case.) A fixed one with a variance, such as a shape a g3 constraint derives, stays in. An NA variance on a parameter that is NOT fixed was masked by .hzr_safe_solve() as non-positive: it withholds the standard error if the prediction uses it (#586). Returns NULL (with a warning) when CLs cannot be computed. Shared by the aggregate and decomposed se.fit paths.

Usage

.hzr_free_vcov(
  vcov_mat,
  p,
  unused = integer(0),
  fixed = NULL,
  param_names = NULL
)

Arguments

Value

list(vcov_use, free_idx), or NULL if unusable.

Finite-difference derivatives of G3 phase Phi and phi w.r.t. shape params

Description

Computes the G3 cumulative intensity, its time derivative, and their partial derivatives with respect to log_tau, gamma, alpha, and eta using central finite differences.

Usage

.hzr_g3_phase_derivatives(time, tau, gamma, alpha, eta, h = 1e-05)

Arguments

Value

Named list with Phi, phi, and 8 derivative vectors.

Gradient of the multiphase log-likelihood

Description

Computes the score vector d\ell / d\theta using analytic chain-rule formulas for log_mu and beta parameters, and finite-difference derivatives (via .hzr_phase_derivatives()) for shape parameters (log_t_half, nu, m): central, except that the m stencil stays on the side of m near m = 0, where the phase family has a cusp.

Usage

.hzr_gradient_multiphase(
  theta,
  time,
  status,
  time_lower = NULL,
  time_upper = NULL,
  x = NULL,
  weights = NULL,
  phases,
  covariate_counts,
  x_list,
  objective = c("likelihood", "sas"),
  sanitize = TRUE,
  phase_cache = NULL,
  ...
)

Arguments

Value

Numeric vector of length length(theta): the gradient. With sanitize = TRUE (the default) a component that cannot be evaluated is 0, and so is the whole vector at an infeasible point (guards the optimizer); with sanitize = FALSE those are NA.

Find the position of each phase's log_mu in the full theta vector

Description

Find the position of each phase's log_mu in the full theta vector

Usage

.hzr_log_mu_positions(phases, covariate_counts)

Arguments

Value

Named integer vector: position of log_mu for each phase.

Interval-censored log-likelihood contribution

Description

The single place the interval-censored contribution is written. Both .hzr_logl_multiphase() and the finite-difference closure inside .hzr_gradient_multiphase() delegate here, so the optimizer cannot step by the gradient of a different objective than the one it evaluates.

Usage

.hzr_logl_interval(
  cumhaz_lower,
  cumhaz_upper,
  lower,
  upper,
  weights,
  objective = c("likelihood", "sas")
)

Arguments

Details

Callers pass only the interval rows, already subset by status == 2, so this helper never sees the status mask and cannot disagree with a caller about which rows are intervals.

Value

Scalar summed contribution; -Inf for infeasible parameters.

Log-likelihood for multiphase additive hazard model

Description

Log-likelihood for multiphase additive hazard model

Usage

.hzr_logl_multiphase(
  theta,
  time,
  status,
  time_lower = NULL,
  time_upper = NULL,
  x = NULL,
  weights = NULL,
  phases,
  covariate_counts,
  x_list,
  objective = c("likelihood", "sas"),
  return_gradient = FALSE,
  return_hessian = FALSE,
  ...
)

Arguments

Value

Scalar log-likelihood. Returns -Inf for infeasible parameters.

Normal deviates from a private Mersenne Twister

Description

Reproduces ⁠set.seed(seed); rnorm(n)⁠ under R's default RNGkind(), bit for bit, without touching .Random.seed. Written out rather than called because calling it means writing to the global environment, and changing the numbers would change which optimum existing multiphase fits reach. The 32-bit words are held in doubles: every product stays below 2^53, so the arithmetic is exact.

Usage

.hzr_mt_normal(seed, n)

Arguments

Details

The generator is that of Matsumoto and Nishimura (1998). The seeding, the zero-word substitution and the inversion are derived from R's own src/main/RNG.c and src/nmath/snorm.c, and a test pins them against set.seed() and rnorm(). See system.file("COPYRIGHTS", package = "TemporalHazard").

Value

Numeric vector of length n.

References

Matsumoto M, Nishimura T (1998). Mersenne Twister: a 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Transactions on Modeling and Computer Simulation, 8(1), 3-30. doi:10.1145/272991.272995

Compute total cumulative hazard H(t|x) for multiphase model

Description

Compute total cumulative hazard H(t|x) for multiphase model

Usage

.hzr_multiphase_cumhaz(
  time,
  theta,
  phases,
  covariate_counts,
  x_list,
  per_phase = FALSE
)

Arguments

Value

If per_phase = FALSE: numeric vector of length n (total H(t|x)). If per_phase = TRUE: named list with ⁠$total⁠ and one element per phase.

Compute total instantaneous hazard h(t|x) for multiphase model

Description

Compute total instantaneous hazard h(t|x) for multiphase model

Usage

.hzr_multiphase_hazard(time, theta, phases, covariate_counts, x_list)

Arguments

Value

Numeric vector of length n.

Transform multiphase theta from internal to user-facing scale

Description

Converts log_mu -> mu, log_t_half -> t_half; nu, m, beta pass through.

Usage

.hzr_multiphase_theta_user(theta, phases, covariate_counts)

Arguments

Value

Named numeric vector on the reporting scale.

Fit a multiphase additive hazard model via maximum likelihood

Description

Assembles starting values from phase specifications, resolves per-phase design matrices, and delegates to .hzr_optim_generic().

Usage

.hzr_optim_multiphase(
  time,
  status,
  time_lower = NULL,
  time_upper = NULL,
  x = NULL,
  theta_start = NULL,
  weights = NULL,
  control = list(),
  phases,
  objective = c("likelihood", "sas"),
  formula_global = NULL,
  data = NULL
)

Arguments

Value

List with par (internal scale), value, convergence, vcov, etc.

Parse Surv() formula for hazard modeling

Description

Extracts time, status, time_lower, time_upper, and predictors from a formula of the form Surv(time, status) ~ x1 + x2 + .... Supports right-censored, left-censored, interval-censored, and counting-process (start-stop) data.

Usage

.hzr_parse_formula(formula, data)

Arguments

Details

For counting-process (start-stop) data, use Surv(start, stop, event). The start times are returned as time_lower and stop times as time, enabling the likelihood to compute H(stop) - H(start) per epoch.

Surv() and this package code censoring status differently, so the returned status is translated, not passed through:

Under type = "interval", Surv() reuses the time2 column to hold the status of any non-interval row, so an upper bound is read only where the row is genuinely interval-censored.

Value

A list with elements: time, status, time_lower, time_upper, x

Parse hazard binary output

Description

Parses tabular output from the C hazard binary into a structured data frame. Expects columns: parameter, estimate, StdErr (or SE), z-value (or z), p-value (or pval). The parser is flexible and case-insensitive for column matching.

Usage

.hzr_parse_hazard_output(output_text, theta = NULL)

Arguments

Value

Data frame with columns: parameter, estimate, std_err, z_stat, p_value

Derivatives of phase cumulative and instantaneous hazard w.r.t. shape params

Description

Computes \Phi_j(t), \phi_j(t), and their derivatives with respect to log_t_half, nu, and m using finite differences on hzr_decompos(): central in log_t_half and nu, and in m central except near m = 0, where the stencil keeps the sign of m. hzr_decompos() changes formula at m = 0 and the two sides meet in a cusp, so for m >= 0 a stencil that would reach 0 becomes one-sided forward (second order), and for m < 0 the step is capped at 1% of ⁠|m|⁠ (floored at 1e-10) and becomes one-sided backward if it would still reach 0. The log_t_half derivative is a central difference in log_t_half itself, so its step is proportional to t_half at every scale. It is NaN where the two points of the difference cannot both be used: t_half too small or too large to step, or a point at which hzr_decompos() fails. A "hazard" phase's value is -log_surv from hzr_decompos(), which keeps its accuracy far past saturation, so the derivative does too.

Usage

.hzr_phase_derivatives(
  time,
  t_half,
  nu,
  m,
  type = c("cdf", "hazard", "constant")
)

Arguments

Value

Named list:

Phi

Cumulative hazard contribution \Phi(t).

phi

Instantaneous hazard contribution \phi(t) = d\Phi/dt.

dPhi_dlog_thalf

d\Phi / d(\log t_{1/2}).

dPhi_dnu

d\Phi / d\nu.

dPhi_dm

d\Phi / dm.

dphi_dlog_thalf

d\phi / d(\log t_{1/2}).

dphi_dnu

d\phi / d\nu.

dphi_dm

d\phi / dm.

Build a logical mask of free (non-fixed) parameters in the full theta vector

Description

Returns a logical vector the same length as the full theta where TRUE means the parameter is free (to be optimized) and FALSE means it is held fixed at its starting value.

Usage

.hzr_phase_free_mask(phases, covariate_counts)

Arguments

Value

Logical vector of length sum(.hzr_phase_n_params(...)).

Does each phase have shape parameters at all?

Description

A constant phase is mu and nothing else. The saturated identifiability message says mu survives while the shape parameters go flat, which is vacuous for a phase that has none (the wording defect in #211).

Usage

.hzr_phase_has_shape(phases)

Arguments

Details

This asks whether the parameters exist, not whether they are free. A phase whose shapes are pinned with hzr_phase(fixed = ) still has them, and the message is true and useful there: it is usually why they were pinned.

Value

Logical vector, one element per phase, in phases order.

Total number of parameters for a phase

Description

Returns the count of free parameters in the theta sub-vector for one phase: 1 (log_mu) + n_shape + n_covariates.

Usage

.hzr_phase_n_params(phase, n_covariates = 0L)

Arguments

Value

Integer.

Number of shape parameters for a phase

Description

Returns 3 (t_half, nu, m) for "cdf" and "hazard" phases, 4 (tau, gamma, alpha, eta) for "g3", and 0 for "constant".

Usage

.hzr_phase_n_shape(phase)

Arguments

Value

Integer: 3, 4 or 0.

Diagnose phases that have effectively left the model

Description

Two distinct failure modes, with different consequences, both silent:

Usage

.hzr_phase_shares(theta, time, phases, covariate_counts, x_list)

Arguments

Details

"absent"

The phase contributes essentially none of \Lambda at any observed time; it has not started by the end of follow-up. Its mu and its shape are unidentified.

"saturated"

The phase's \Phi is effectively constant across the observed times: a cdf phase whose half-life is far shorter than the first observation has already finished. It then contributes \mu \cdot \Phi \approx \mu, a constant offset, so mu stays well identified while the shape parameters (t_half, nu, m) go exactly flat: the likelihood is unchanged whether they are pinned or fitted.

Share is taken of \Lambda, not of h, because every row type's contribution runs through \Lambda(t). A phase can supply almost none of the instantaneous hazard late in follow-up and still be perfectly well identified through the offset it already contributed, which is why the hazard is the wrong basis for this test.

The maximum over times is the right summary rather than the mean: a phase mattering only in the first week is identified, and a mean over long follow-up would hide it.

Value

data.frame with one row per phase: share (largest share of \Lambda at any observed time) and variation (relative range of the phase's contribution across observed times, NA when the phase carries covariates; mu then varies by row and the two sources of variation cannot be separated from the contribution alone).

The shape is the same whatever happens: if no observed time carries a usable total, both columns are NA rather than the frame being NULL. fit$phase_share is then one type for a caller to handle rather than two, and "could not be measured" stays distinct from "measured as zero".

Extract starting values from a phase specification

Description

Returns initial theta sub-vector on the estimation (internal) scale: log(mu), then log(t_half), nu, m for "cdf"/"hazard" phases or log(tau), gamma, alpha, eta for "g3" (nothing for "constant"), followed by zeros for covariate coefficients.

Usage

.hzr_phase_start(phase, n_covariates = 0L, mu_start = 0.1)

Arguments

Value

Unnamed numeric vector of starting values.

Generate parameter names for a phase's sub-vector

Description

Builds the named labels for the parameter block that belongs to a single phase in the full theta vector. The block layout is:

Usage

.hzr_phase_theta_names(phase, phase_name, covariate_names = character(0))

Arguments

Details

• For "cdf"/"hazard": ⁠[log_mu, log_t_half, nu, m, beta_1, ..., beta_p]⁠

• For "g3": ⁠[log_mu, log_tau, gamma, alpha, eta, beta_1, ..., beta_p]⁠

• For "constant": ⁠[log_mu, beta_1, ..., beta_p]⁠

Value

Character vector of parameter names.

Apply a scale transform + back-transform for CLs

Description

Apply a scale transform + back-transform for CLs

Usage

.hzr_predict_cl_from_se(fit, se_nat, level, scale)

Arguments

Value

data.frame with columns fit, se.fit, lower, upper.

Jacobian of multiphase cumulative-hazard predictions

Description

Only cumulative_hazard and survival are delegated here; hazard and linear_predictor are rejected upstream for multiphase models.

Usage

.hzr_predict_jacobian_multiphase(
  theta,
  time,
  phases,
  covariate_counts,
  x_list,
  p,
  per_phase = FALSE
)

Arguments

Value

Numeric ⁠n x p⁠ Jacobian of H(t|x) (default), or a named list of per-phase ⁠n x p⁠ Jacobians when per_phase = TRUE.

Numeric Jacobian of a single-distribution prediction w.r.t. theta

Description

Numeric Jacobian of a single-distribution prediction w.r.t. theta

Usage

.hzr_predict_jacobian_numeric(predict_fn, theta)

Arguments

Value

Numeric n x p Jacobian.

Jacobian of Weibull predictions with respect to theta

Description

Jacobian of Weibull predictions with respect to theta

Usage

.hzr_predict_jacobian_weibull(type, theta, time, x, p)

Arguments

Value

Numeric n x p Jacobian, with respect to log(mu), nu and beta.

Per-row standard errors via the delta method

Description

Per-row standard errors via the delta method

Usage

.hzr_predict_se_from_jacobian(J, vcov)

Arguments

Value

Numeric vector of length n.

Compute point predictions + delta-method CLs

Description

Dispatched from predict.hazard() when se.fit = TRUE.

Usage

.hzr_predict_with_se(
  object,
  type,
  time = NULL,
  x = NULL,
  x_list = NULL,
  cov_counts = NULL,
  phases = NULL,
  level = 0.95,
  diff_fn,
  conf_type = c("log-log", "logit")
)

Arguments

Details

DELTA-METHOD TARGET

The quantity we actually differentiate depends on the prediction type and the CL scale the type uses:

type = "cumulative_hazard" -> target = H, log-scale CL type = "survival" -> target = H, log-log-survival CL (final fit is exp(-H); reported se.fit is S * se(H), the SE of S) type = "hazard" -> target = exp(eta) (single-dist only), log-scale CL type = "linear_predictor" -> target = eta, natural-scale CL

The caller supplies diff_fn(theta) that returns the target vector of length n. For Weibull and multiphase we build J analytically; for exp / loglogistic / lognormal we fall through to a numeric jacobian of diff_fn.

Value

data.frame with columns fit, se.fit, lower, upper.

Per-phase + total cumulative-hazard predictions with delta-method CLs

Description

Long-format companion to .hzr_predict_with_se() for multiphase models under decompose = TRUE. Computes log-scale CLs for the total cumulative hazard and for each phase's additive contribution ⁠H_j(t) = mu_j(x) Phi_j(t)⁠. Per-phase CLs use only that phase's parameter block, so they do NOT sum to the total CL (cross-phase covariance contributes only to the total).

Usage

.hzr_predict_with_se_decomposed(
  object,
  time,
  x_list,
  cov_counts,
  phases,
  level = 0.95
)

Arguments

Value

Long data.frame(time, component, fit, se.fit, lower, upper); component is an ordered factor with levels c("total", names(phases)).

Select the phase whose log_mu will be solved by conservation

Description

Chooses the phase contributing the largest share of total cumulative hazard at the current theta, matching the C HAZARD SETCOE strategy.

Usage

.hzr_select_fixmu_phase(
  theta,
  time,
  status,
  phases,
  covariate_counts,
  x_list,
  weights = NULL,
  time_lower = NULL
)

Arguments

Details

When one phase dominates by an extreme margin (>10x the median contribution) it is typically due to poorly-scaled shape parameters at the starting theta rather than a genuine signal that the phase should absorb all events. Selecting such a phase as the fixmu phase traps the optimizer: CoE keeps rescaling that phase's log_mu upward while the optimizer tries to drive it to near-zero. We therefore exclude extreme outliers and select among the remaining phases. If all phases are outliers (or only one phase exists) the largest is used as a fallback.

Value

Character: name of the phase to fix.

Split the full theta vector into per-phase sub-vectors

Description

Split the full theta vector into per-phase sub-vectors

Usage

.hzr_split_theta(theta, phases, covariate_counts)

Arguments

Value

Named list of numeric vectors, one per phase.

Deterministic multi-start perturbations

Description

Draws the offsets that shift the starting values for every multiphase start after the first. They are the numbers ⁠set.seed(seed); rnorm()⁠ gives under R's default generators, computed by .hzr_mt_normal() without R's stream, so a fit is reproducible without set.seed(), does not depend on RNGkind(), and neither reads nor writes .Random.seed.

Usage

.hzr_start_perturbations(n_starts, n_par, seed, sd = 0.5)

Arguments

Value

A list of max(0, n_starts - 1) numeric vectors of length n_par; empty for a single start, which is never perturbed.

Extract parameters from a phase sub-vector (internal scale)

Description

Extract parameters from a phase sub-vector (internal scale)

Usage

.hzr_unpack_phase_theta(theta_j, phase)

Arguments

Value

Named list with log_mu, log_t_half, nu, m, beta.

Validate a list of phase specifications

Description

Checks that phases is a non-empty named list of hzr_phase objects. Auto-names unnamed phases as phase_1, phase_2, etc.

Usage

.hzr_validate_phases(phases)

Arguments

Value

The validated (and possibly auto-named) list, invisibly.

Compute log(1 + exp(x)) with overflow/underflow protection

Description

Compute log(1 + exp(x)) with overflow/underflow protection

Usage

.log1pexp(x)

Arguments

Value

Numeric vector: log(1 + exp(x)).

Compute log(exp(x) - 1) with numerical stability

Description

Compute log(exp(x) - 1) with numerical stability

Usage

.log_expm1(x)

Arguments

Value

Numeric vector: log(exp(x) - 1).

AVC: Atrioventricular Canal Repair

Description

Survival data for 310 patients who underwent repair of atrioventricular septal defects (congenital heart disease) at the Cleveland Clinic between 1977 and 1993. Exhibits two identifiable hazard phases: an early post-operative risk and a constant late phase.

Usage

avc

Format

A data frame with 310 rows and 11 variables:

study

Patient identifier

status

NYHA functional class (1–4). A covariate, not a censoring status: the event indicator is dead.

inc_surg

Surgical grade of AV valve incompetence

opmos

Date of operation (months since January 1967)

age

Age at repair (months)

mal

Malalignment indicator (0/1)

com_iv

Interventricular communication indicator (0/1)

orifice

Associated cardiac anomaly indicator (0/1)

dead

Death indicator (1 = dead, 0 = censored)

int_dead

Follow-up interval to death or last contact (months)

op_age

Interaction term: opmos x age

Source

Blackstone, Naftel, and Turner (1986) doi:10.1080/01621459.1986.10478314. Cleveland Clinic Foundation.

See Also

vignette("fitting-hazard-models"), vignette("prediction-visualization")

Other datasets: cabgkul, omc, tga, uslife2023, valves

Examples

data(avc)
avc <- na.omit(avc)

# Kaplan-Meier survival
km <- survival::survfit(survival::Surv(int_dead, dead) ~ 1, data = avc)
plot(km, xlab = "Months after AVC repair", ylab = "Survival",
     main = "AVC: Kaplan-Meier survival estimate")


# Two-phase hazard fit (early CDF + constant -- what AVC supports)
fit <- hazard(
  survival::Surv(int_dead, dead) ~ 1, data = avc,
  dist = "multiphase",
  phases = list(
    early    = hzr_phase("cdf", t_half = 0.5, nu = 1, m = 1),
    constant = hzr_phase("constant")
  ),
  fit = TRUE, control = list(n_starts = 5, maxit = 1000)
)
summary(fit)

CABGKUL: Primary Isolated Coronary Artery Bypass Grafting (KU Leuven)

Description

Survival data for 5,880 patients who underwent primary isolated CABG at KU Leuven, Belgium, between 1971 and July 1987. The simplest dataset structure (intercept-only, right-censored) with large sample size exercising all three temporal hazard phases.

Usage

cabgkul

Format

A data frame with 5880 rows and 2 variables:

int_dead

Follow-up interval to death or last contact (months)

dead

Death indicator (1 = dead, 0 = censored)

Source

KU Leuven cardiac surgery registry. Primary benchmark dataset for C binary parity testing.

See Also

vignette("fitting-hazard-models")

Other datasets: avc, omc, tga, uslife2023, valves

Examples

data(cabgkul)

# Kaplan-Meier survival
km <- survival::survfit(survival::Surv(int_dead, dead) ~ 1, data = cabgkul)
plot(km, xlab = "Months after CABG", ylab = "Survival",
     main = "CABGKUL: Kaplan-Meier survival (n = 5,880)")


# Single-phase Weibull fit with parametric overlay
fit <- hazard(survival::Surv(int_dead, dead) ~ 1, data = cabgkul,
              dist = "weibull", theta = c(mu = 0.10, nu = 1.0), fit = TRUE)
t_grid <- seq(0.01, max(cabgkul$int_dead) * 0.9, length.out = 200)
surv   <- predict(fit, newdata = data.frame(time = t_grid),
                  type = "survival")
plot(km, xlab = "Months after CABG", ylab = "Survival",
     main = "CABGKUL: Weibull vs. Kaplan-Meier")
lines(t_grid, surv, col = "blue", lwd = 2)
legend("bottomleft", c("KM", "Weibull"), col = c("black", "blue"),
       lty = 1, lwd = c(1, 2))

Extract coefficients from hazard model

Description

Extract coefficients from hazard model

Usage

## S3 method for class 'hazard'
coef(object, ...)

Arguments

Value

A named numeric vector of fitted parameter estimates, or NULL if the model has not been fitted (fit = FALSE).

Examples

set.seed(1)
fit <- hazard(time = rexp(30, 0.5), status = rep(1L, 30),
              theta = c(0.3, 1.0), dist = "weibull", fit = TRUE)
coef(fit)

Build and optionally fit a hazard model

Description

Creates a hazard object and optionally fits it via maximum likelihood. This mirrors the argument-oriented workflow of the legacy HAZARD C/SAS implementation: supply starting values in theta and the function will optimize to produce fitted estimates.

Usage

hazard(
  formula = NULL,
  data = NULL,
  time = NULL,
  status = NULL,
  time_lower = NULL,
  time_upper = NULL,
  x = NULL,
  time_windows = NULL,
  theta = NULL,
  dist = "weibull",
  phases = NULL,
  fit = FALSE,
  weights = NULL,
  control = list(),
  objective = c("likelihood", "sas"),
  ...
)

Arguments

Details

Control parameters:

• maxit: Maximum iterations of the quasi-Newton (BFGS) optimizer (default 1000), applied to each start. A fit whose optimizer reports convergence but fails SAS's gradient test is continued with stats::nlm() under the same limit (see "Convergence"), so raising maxit lets that continuation run further too. The Nelder-Mead warm-up that a multiphase fit with fixed parameters may run first has its own limit, which maxit does not change. It must be a single finite number of at least 1 (a fraction is truncated, as PROC HAZARD truncates MAXITER): anything else stops hazard() at once, fitted or not (#541). A fit with no iterations is its starting values, not an estimate; to evaluate a model at parameters you supply, use hzr_evaluate().

• n_starts: Number of optimization starts for multiphase fits (default 5). Each start after the first offsets the initial values. The offsets are computed by the package without using R's random number stream, so a multiphase fit is reproducible without set.seed() and neither reads nor changes that stream.

• start_seed: Seed selecting the ensemble of multi-start offsets (default 3). Any whole number within integer range, negative included. Fits are reproducible at any value; vary it to probe a different set of starts when a fit is suspected of sitting in a local optimum, and compare the resulting objective values. A fractional value is rejected rather than truncated, since 3.9 and 3 would otherwise select the same ensemble without saying so.

• phase_share_tol: Threshold for the multiphase identifiability warning (default 1e-8). A phase is reported as having left the model when it supplies less than this share of the cumulative hazard at every observed time (it never started; neither its mu nor its shape is identified), or when its contribution varies by less than this relative amount across them (it finished before the first observation and acts as a constant offset; mu stays identified, the shape parameters do not). A third condition is a property of the observed times rather than of any phase: when their own relative range falls below this threshold and no phase's contribution varies across them, they separate no phase from any other and none of the shapes is identified. That verdict is withheld when the fit supplies counting-process entry times or interval bounds that add enough evaluation points the measure did not see. "Enough" is counted, not assumed: with the event times tied, each added point supplies one more functional of the parameters, and one more than the number of added points must reach the number of free parameters. A zero entry time, or a bound equal to an event time, adds nothing and does not count. The count is deliberately conservative: it does not credit the extra functional that a zero entry time makes separately observable, so a fit can be warned about while being identified. For the same reason the per-phase measures can overstate what an interval-censored or left-truncated fit loses, since a phase flat across the event times may still be identified through the bounds; the saturated warning says so where such points exist, rather than claiming the likelihood is unchanged (#228). The measured shares are kept on the fit as fit$fit$phase_share. Raise it to catch marginal phases, set it to 0 to silence the check.

• reltol: Relative convergence tolerance on the objective the optimizer minimizes, the negative log-likelihood (default 1e-5). The stored fit$fit$objective has the opposite sign (see Value). BFGS stops when an iteration reduces it by less than ⁠reltol * (|objective| + reltol)⁠, so the stopping gap grows with the size of the log-likelihood: about 0.0024 at a log-likelihood of -240. On a flat surface plain BFGS can stop that far short of the optimum and still report convergence, so hazard() then applies SAS/C HAZARD's relative-gradient test and, when the stop fails it, continues with stats::nlm(); see the "Convergence" section.

• conserve: Apply Conservation of Events (dist = "multiphase" only; default TRUE). CoE counts exact events, so it is automatically disabled whenever any status falls outside {0, 1} (which interval or left censoring guarantees) and whenever the model has fewer than two phases. On a fitted multiphase object the outcome is recorded next to the request, both fields living under fit$spec$control:

  • fit$spec$control$conserve_applied: logical, whether CoE was actually applied;

  • fit$spec$control$conserve_disabled_reason: one of "not_requested", "unsupported_censoring", "single_phase", "no_events", "setup_failed", or NA when CoE was applied.

  Read fit$spec$control$conserve_applied, not fit$spec$control$conserve: the latter says only what you asked for.

• shape_param_count: The number of baseline parameters at the front of theta: the scale and any shape parameters, so 2 for "weibull" and 1 for "exponential". For a single-distribution model only. The fit does not use it: theta is checked against the distribution's own count, and so are hzr_stepwise(), its refits and its score test. Only an internal helper that writes a SAS input deck for parity checks reads it. A multiphase fit derives its own layout, so nothing reads it there.

The elements above are accepted without a warning: maxit and reltol for every model, shape_param_count for a single-distribution model, and n_starts, start_seed, phase_share_tol and conserve for dist = "multiphase" (#376). Any other element draws one warning that names it and says why it has no effect, and the fit proceeds unchanged, as stats::optim() does for unknown control names. The element is dropped before the fit, so a name such as n_starts_extra cannot be read as n_starts, and fit$spec$control keeps none of the ignored elements. That covers:

• a name no fit reads, such as the misspelling n_startz, and an unnamed element;

• abstol (read only by a bounded optimizer no fit uses), method, condition, nocov and nocor, which earlier versions documented as accepted without reading them;

• fix and quasi, which no fit has ever read. A fit given fix was never constrained, so results obtained with it may be affected; hold a parameter with hzr_phase(fixed = ) on a multiphase phase. A single-distribution model has no mechanism for fixing a parameter.

• a multiphase element such as n_starts given to a single-distribution fit, and shape_param_count given to a multiphase one.

No name in control is an error. A bad maxit stops hazard() at once (see above); a bad value for another element the fit reads still stops a fit (fit = TRUE) where it is read. hzr_stepwise() and hzr_bootstrap() pass control to every candidate refit, and an error there would count as a failed candidate, so a screen would report success having tested nothing.

SAS ⁠PROC HAZARD⁠ options with no control equivalent: NOCOV and NOCOR only suppress printed output, and hazard() prints nothing until asked. ⁠CONDITION=⁠ stops SAS's optimizer on a condition-number test that hazard() does not have. QUASI is hazard()'s optimizer already: the fits it runs use BFGS (a multiphase fit may run a Nelder-Mead warm-up first, and a stop that fails SAS's gradient test continues with stats::nlm()). SAS jobs often write ⁠STEEPEST QUASI⁠ together, steepest descent first; there is no steepest-descent option and no two-stage strategy. The multiphase likelihood is multimodal, so a different descent path can land on a different optimum: a fit translated from a job using STEEPEST may not reproduce SAS's estimates, and hzr_translate_sas() records the keyword as untranslated rather than dropping it.

Censoring status coding:

• 1: Exact event at time

• 0: Right-censored at time

• -1: Left-censored with upper bound at time_upper \(or time\)

• 2: Interval-censored in the interval \(time_lower, time_upper\)

Time-varying coefficients:

• If time_windows is supplied, predictors are expanded to piecewise window interactions so each window has its own coefficient vector.

• This is implemented as design-matrix expansion, so the existing likelihood engines remain unchanged.

Value

An object of class hazard, a named list with components: call (the matched call), spec (model specification: dist, control, time_windows, phases), data (input data: time, status, x, weights, etc.), fit (optimisation results: theta, objective, the maximized log-likelihood, so larger is better, or PROC HAZARD's objective under objective = "sas" (summary() reports which, as its objective); converged, se, vcov, counts, message, and rel_gradient, rel_gradient_reason and polish_code, the SAS/C acceptance test described under "Convergence"; all NULL when fit = FALSE; multiphase fits add starts, one row per optimisation start with its status ("ok", "nonconverged", "infeasible", "nonfinite" or "error"), objective (on the same scale as the fit's, so the best start has the largest; NA unless the start reached a point where the likelihood is defined, and recorded as the optimizer returned it, before any Conservation of Events adjustment to the conserved phase's scale), convergence (the optim code, 0 for success), whether it was the best and so the reported fit, and the message of any error. A start that stops at maxit has a finite objective and can win the selection, so status distinguishes it from one that converged. Every fit carries weak, which takes one of three values: a list, when the likelihood is near-flat along a parameter combination, giving the params spanning that direction, their squared loadings (weights), the strongest pairwise correlation among them, the Hessian rcond and n_directions, the number of near-flat directions found (when a single parameter carries the direction on its own, the list has one params entry, single = TRUE, its estimate, and correlation = NA, and se_metric, its standard error on the log scale; only a g3 shape, gamma, alpha or eta, or a multiphase phase's log_mu, log_t_half or log_tau, is named this way, the last three when their standard error exceeds 10 log units, and with weights = NA); NULL when the fit was examined and is well identified; and NA when the check could not run because no usable Hessian was available, which includes an unfitted object, an install without the suggested numDeriv, a Hessian that is not positive-definite (unless a ridge was found, which is reported), and an ill-conditioned fit in which an estimated parameter has no finite variance and no ridge was found among the others, since that parameter was not examined. Test with is.list(fit$fit$weak), not !is.null(): the NA case has not been examined and must not be read as a clean result), engine (implementation tag, "native-r-m2"), and two fields recording what the fit did not do: degraded, a character vector of the steps not performed, in the fixed order "fitting", "standard_errors", "conserved_phase_variance", "weak_direction_check", "boundary_check", "conservation_of_events", and empty when nothing was lost; and degraded_causes, a character vector with the same names giving the reason for each. print() and summary() always show them as a "Not done in this run" block, which reads "none" when nothing was lost. fit$fit$weak is NA exactly when "weak_direction_check" is listed.

fit$fit$boundary is its sibling and takes the same three states, for phases whose fit sits where their parameterisation breaks down: NULL when the check ran and found nothing, a list of records when it found something, and NA when it did not run — and it is NA exactly when "boundary_check" is listed in degraded, with the reason in degraded_causes. Each record carries mechanism, phase, parameter and a printable detail. The mechanisms are "unbounded_phase", a "hazard" phase whose t_half is below the first observed time; "phase_discontinuity", a "cdf" or "hazard" phase whose shape has collapsed to a step the observed times cannot resolve (it can lie inside the data); and two made at setup, before the fit: "g3_alpha_one" (a g3 phase with alpha fixed at 1, re-expressed as PROC HAZARD does, see hzr_phase()) and "g3_fixge2_alpha_start" (a free alpha start moved to 2/3 under constraint = "eta_gamma"); and two made after it. "g3_corner_supremum" is a g3 phase under constraint = "eta_gamma" whose estimated gamma converged below 1000 although the log-likelihood is higher at a much larger gamma, with any fixed tau or alpha held. The estimate is then not the maximum-likelihood one, and the supremum may lie at gamma = Inf. A fixed gamma, and a fit with left- or interval-censored rows, are not examined. "coe_no_events_left" is a fit under Conservation of Events whose log-likelihood is higher with the conserved phase's scale sent to zero, where no events remain for that phase. Both records also carry gain and the higher point as certificate_theta with its certificate_loglik (the corner record adds gamma_hat); the estimates are not changed. "coe_phase_vanished" is a fit under Conservation of Events that has itself run to that boundary, or beside it: the conserved phase's largest share of the cumulative hazard is below control$phase_share_tol. Its record carries share and tol, and warned_by = "phase_share": the identifiability warning has already reported the phase, so this record raises no warning of its own. A fit can carry several records, in no guaranteed order, so select them by mechanism. Only rows the likelihood reads count as observed times. A fit with any record that has no warned_by raises one warning whose classes are "hzr_" plus each such record's mechanism, all inheriting "hzr_boundary", so one handler catches the whole family.

The fitted model

For dist = "multiphase" the cumulative hazard, instantaneous hazard, and survival are

H(t \mid \mathbf{x}) = \sum_{j=1}^{J} \mu_j(\mathbf{x}) \, \Phi_j(t), \qquad h(t \mid \mathbf{x}) = \sum_{j=1}^{J} \mu_j(\mathbf{x}) \, \varphi_j(t), \qquad S(t \mid \mathbf{x}) = \exp\!\bigl(-H(t \mid \mathbf{x})\bigr)

where \mu_j(\mathbf{x}) = \exp(\alpha_j + \mathbf{x}_j^\top \beta_j) and the temporal shapes \Phi_j, \varphi_j are set by each phase's type (see hzr_phase()). The proportional-hazards single-phase families ("weibull", "exponential") are the special case J = 1, with covariates acting multiplicatively on one temporal shape. The "loglogistic" (proportional-odds) and "lognormal" (accelerated-failure-time) families place covariates differently (on the odds of failure and the log-time location, respectively), so they are separate parameterizations, not special cases of this additive form. Parameters are estimated on an unconstrained internal scale (e.g. \log\mu, \log t_{1/2}) and transformed back for reporting; see vignette("mf-mathematical-foundations").

Convergence

The optimizer stops when an iteration improves the log-likelihood by less than control$reltol relative to its size, and on a flat ridge that can happen well short of the maximum. SAS/C HAZARD accepts an optimum only on a different test, the relative gradient \max_i |g_i| \max(|x_i|, 1) / \max(|\ell|, 1) \le \epsilon^{1/3}, about 6e-6, and hazard() applies it too: when the optimizer reports convergence and the test fails, the fit is continued with stats::nlm() at SAS's tolerances, and the continued point is kept if it improves the log-likelihood.

Every fit records the result in fit$fit$rel_gradient and, when the continuation improved the fit, its termination code in fit$fit$polish_code. print() and summary() show both. rel_gradient is NA when the test was not applied (the optimizer did not report convergence) or the gradient cannot be evaluated at the estimates; NA is never reported as a pass. Neither is it always a failure: some routes to it, such as a non-converged stop, do say the estimates are unreliable, while others, such as a finite-difference score that needed a point where the log-likelihood is not usable, say nothing against them. Which route it took is recorded in fit$fit$rel_gradient_reason, NA_character_ when the test did run, and print() and summary() show it. Under Conservation of Events the conserved scale is re-solved from the other parameters at every step, as SAS/C does, and the score used by the optimizer, the continuation and the test includes how that scale moves with them. A fit can still end with the test not met, under Conservation of Events or without it. A warning is raised only for code 4, the iteration limit (raise control$maxit), and code 5, where the log-likelihood kept rising along some direction and the model may have no maximum. Codes 2 and 3, where SAS/C prints a caution, are recorded without one, except for the single-distribution fits: there a converged fit with a relative gradient above 1e-3 that did not stop on code 4 or 5 (in practice code 2 or 3, or no code, which is recorded only when the continuation's point was kept: none is when it found no better point, stats::nlm() raised an error, or its minimum was not finite) warns, with class "hzr_possible_false_maximum", that the fit may not be a maximum, and suggests other starting values or centring or rescaling the covariates. Badly scaled covariates can put a good fit above 1e-3, where a restart does not help and rescaling does. hzr_bootstrap() counts the replicates in which any fit (the base refit, a stepwise refit or the final fit) meets the same rule, and warns once. Such a stop can be a false maximum far below the best one from an ordinary start, and the gradient test alone cannot always tell it from a good fit, so converged is left as it is. The warning is not a guarantee: a fit can stop short of its maximum with a relative gradient below 1e-3 and raise nothing. The test is relative to the size of the log-likelihood, so a fit that meets it is within SAS's tolerance of the maximum, not exactly at it.

The optimizer treats a score it cannot use as zero, so it can stop on those zeros far from the maximum. For the single-distribution fits (exponential, Weibull, lognormal and loglogistic), a stop where the score at the estimates is not finite, or could not be computed, is therefore reported as converged = FALSE, with a warning of class "hzr_unverified_convergence"; rel_gradient_reason names which.

Baseline distributions

The dist argument selects the parametric form of the baseline hazard. The four single-distribution families differ in the shape the hazard traces over follow-up; choose by what the risk is expected to do over time. "multiphase" is the general additive model that lets several such shapes coexist.

"weibull": monotone rising or falling hazard (default)

The workhorse parametric model: H(t \mid \mathbf{x}) = (\mu t)^\nu \exp(\eta), with hazard h \propto t^{\nu - 1}. The single shape \nu makes risk increase over time (\nu > 1), decrease (\nu < 1), or stay flat (\nu = 1). Use it as the default when a single monotone trend describes the hazard. The scale \mu is the baseline at \mathbf{x} = 0, so a covariate far from zero (a calendar year, an age in days) can push it beyond what a number can hold: it is then reported as Inf, as 0, or as a value too small to keep its digits, with a warning of class "hzr_unrepresentable_scale". The fit keeps \log \mu and its standard error, which summary() shows as a log(mu) row and predict() reads; coef() shows \mu as it is, and vcov() gives NA for it. Centering or rescaling the covariate reports \mu itself, without changing the model.

"exponential": constant hazard

The memoryless special case \nu = 1: a time-invariant baseline rate, H(t \mid \mathbf{x}) = \mu t \exp(\eta). Use it when the event rate does not change with follow-up time (the constant background risk also appears as the "constant" phase in a multiphase model).

"loglogistic": unimodal (rise-then-fall) hazard

A log-logistic proportional-odds form (covariates act multiplicatively on the odds of failure, \exp(\eta), not as an AFT time shift) whose hazard rises to a single peak and then declines when the shape exceeds 1 (and is monotone decreasing otherwise), with heavier tails than the log-normal. Use it when risk climbs to an early peak and then eases off.

"lognormal": early-peaking, resolving hazard

An accelerated-failure-time form in which \log time is Gaussian; the hazard rises to an early peak and then decays toward zero. Use it for risk that is concentrated early and resolves over time.

"multiphase": additive N-phase hazard

Sums several phase shapes into one model, H = \sum_j \mu_j(\mathbf{x}) \Phi_j(t), so the overall hazard can fall, level off, and rise again within one fit. Requires phases; see hzr_phase() for the available phase shapes. This is the form that reproduces the classic C/SAS HAZARD models.

Note

objective = "sas" exists to reproduce legacy ⁠PROC HAZARD⁠ runs and must not be used for new analyses. It is a density, not a probability: it is inconsistent for wide intervals, where the two forms differ materially (22 log-likelihood units on the esophagectomy reference fit). The default is the statistically correct interval likelihood. See inst/dev/SAS-INTERVAL-OBJECTIVE-DESIGN.md for the derivation and the four-reference evidence.

References

Blackstone EH, Naftel DC, Turner ME Jr. The decomposition of time-varying hazard into phases, each incorporating a separate stream of concomitant information. J Am Stat Assoc. 1986;81(395):615–624. doi:10.1080/01621459.1986.10478314

Rajeswaran J, Blackstone EH, Ehrlinger J, Li L, Ishwaran H, Parides MK. Probability of atrial fibrillation after ablation: Using a parametric nonlinear temporal decomposition mixed effects model. Stat Methods Med Res. 2018;27(1):126–141. doi:10.1177/0962280215623583

See Also

predict.hazard() for survival/cumulative-hazard predictions, summary.hazard() for model summaries, hzr_phase() for specifying multiphase temporal shapes.

Vignettes with worked examples: vignette("fitting-hazard-models"): single-phase through multiphase fitting, vignette("prediction-visualization"): prediction types and decomposed hazard plots, vignette("inference-diagnostics"): bootstrap CIs and model diagnostics.

Examples

# -- Univariable Weibull ----------------------------------------------
set.seed(1)
time   <- rexp(50, rate = 0.3)
status <- sample(0:1, 50, replace = TRUE, prob = c(0.3, 0.7))
fit <- hazard(time = time, status = status,
              theta = c(0.3, 1.0), dist = "weibull", fit = TRUE)
summary(fit)

# -- Formula interface with covariates --------------------------------
set.seed(1001)
n   <- 180
dat <- data.frame(
  time   = rexp(n, rate = 0.35) + 0.05,
  status = rbinom(n, size = 1, prob = 0.6),
  age    = rnorm(n, mean = 62, sd = 11),
  nyha   = sample(1:4, n, replace = TRUE),
  shock  = rbinom(n, size = 1, prob = 0.18)
)

fit2 <- hazard(
  survival::Surv(time, status) ~ age + nyha + shock,
  data    = dat,
  theta   = c(mu = 0.25, nu = 1.10, beta1 = 0, beta2 = 0, beta3 = 0),
  dist    = "weibull",
  fit     = TRUE,
  control = list(maxit = 300)
)
summary(fit2)


# -- Parametric survival with Kaplan-Meier overlay -----------------
if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(ggplot2)

  # Parametric curve on a fine grid at median covariate profile
  t_grid   <- seq(0.05, max(dat$time), length.out = 80)
  curve_df <- data.frame(
    time = t_grid, age = median(dat$age), nyha = 2, shock = 0
  )
  curve_df$survival <- predict(fit2, newdata = curve_df,
                               type = "survival") * 100

  # Kaplan-Meier empirical overlay
  km    <- survival::survfit(survival::Surv(time, status) ~ 1, data = dat)
  km_df <- data.frame(time = km$time, survival = km$surv * 100)

  ggplot() +
    geom_step(data = km_df, aes(time, survival, colour = "Kaplan-Meier")) +
    geom_line(data = curve_df, aes(time, survival,
                                   colour = "Parametric (Weibull)")) +
    scale_colour_manual(
      values = c("Parametric (Weibull)" = "#0072B2",
                 "Kaplan-Meier"         = "#D55E00")
    ) +
    scale_y_continuous(limits = c(0, 100)) +
    labs(x = "Months after surgery", y = "Freedom from death (%)",
         colour = NULL) +
    theme_minimal() +
    theme(legend.position = "bottom")
}



# -- Multiphase model (two phases) ---------------------------------
fit_mp <- hazard(
  survival::Surv(time, status) ~ 1,
  data   = dat,
  dist   = "multiphase",
  phases = list(
    early = hzr_phase("cdf", t_half = 0.5, nu = 2, m = 0,
                       fixed = "shapes"),
    late  = hzr_phase("cdf", t_half = 5,   nu = 1, m = 0,
                       fixed = "shapes")
  ),
  fit     = TRUE,
  control = list(n_starts = 5, maxit = 1000)
)
summary(fit_mp)

# -- Per-phase decomposed cumulative hazard ------------------------
if (requireNamespace("ggplot2", quietly = TRUE)) {
  t_grid <- seq(0.01, max(dat$time), length.out = 100)
  decomp <- predict(fit_mp, newdata = data.frame(time = t_grid),
                    type = "cumulative_hazard", decompose = TRUE)

  df_long <- data.frame(
    time = rep(decomp$time, 3),
    cumhaz = c(decomp$total, decomp$early, decomp$late),
    component = rep(c("Total", "Early (cdf)", "Late (cdf)"),
                    each = nrow(decomp))
  )
  df_long$component <- factor(df_long$component,
    levels = c("Total", "Early (cdf)", "Late (cdf)"))

  ggplot2::ggplot(df_long,
    ggplot2::aes(x = time, y = cumhaz, colour = component,
                 linewidth = component)) +
    ggplot2::geom_line() +
    ggplot2::scale_colour_manual(values = c(
      "Total" = "black", "Early (cdf)" = "#0072B2",
      "Late (cdf)" = "#D55E00"
    )) +
    ggplot2::scale_linewidth_manual(values = c(
      "Total" = 1.2, "Early (cdf)" = 0.6, "Late (cdf)" = 0.6
    )) +
    ggplot2::labs(
      x = "Time", y = "Cumulative hazard H(t)",
      colour = NULL, linewidth = NULL,
      title = "Multiphase decomposition: early + late"
    ) +
    ggplot2::theme_minimal() +
    ggplot2::theme(legend.position = "bottom")
}

Legacy HAZARD to TemporalHazard argument mapping

Description

Returns a formal mapping table that defines how legacy SAS HAZARD/C-style inputs map to hazard(...) arguments in this package.

Usage

hzr_argument_mapping(include_planned = TRUE)

Arguments

Value

A data frame with one row per mapping rule.

See Also

hazard() for the target arguments, hzr_phase() for phase construction, and hzr_decompos()/hzr_decompos_g3() for the parameter-level early/late translations.

Examples

hzr_argument_mapping()
hzr_argument_mapping(include_planned = FALSE)

Bootstrap resampling for hazard model coefficients

Description

Resample data with replacement, refit the hazard model on each replicate, and accumulate coefficient distributions. Returns a tidy data frame of per-replicate estimates with summary statistics. This is the R equivalent of the SAS bootstrap.hazard.sas macro.

Usage

hzr_bootstrap(
  object,
  n_boot = 200L,
  fraction = 1,
  seed = NULL,
  verbose = FALSE,
  scope = NULL,
  direction = c("both", "forward", "backward"),
  criterion = c("score", "wald", "aic"),
  slentry = 0.3,
  slstay = 0.2,
  max_steps = 50L,
  max_move = 4L,
  force_in = character(),
  force_out = character(),
  id = NULL,
  ...
)

## S3 method for class 'hzr_bootstrap'
print(x, digits = 4, ...)

Arguments

Details

When scope is supplied, each replicate instead runs a fresh hzr_stepwise() selection on the resampled data (starting from a fixed-shape refit of object) instead of refitting object's exact formula. This is the R equivalent of the SAS ⁠%HAZBOOT⁠ macro: fit the hazard shape with no covariates (fixing it via hzr_phase(..., fixed = "shapes")), then bootstrap-screen candidate covariates for how often they enter the model. summary$pct then reports the selection frequency across replicates, and summary$mean/sd/⁠ci_*⁠ describe the coefficient distribution conditional on selection.

A replicate in which any fit (the base refit, a stepwise refit or the final fit) may not be a maximum, by the rule hazard() warns on with class "hzr_possible_false_maximum", is counted, and hzr_bootstrap() warns once with that count. Such replicates are kept in the pooled results. Two more kinds of replicate are counted, warned about once and kept in the same way: one whose base refit or final fit stopped on nlm() code 4 or 5 and failed the relative-gradient test, and, in selection mode, one whose screen had a candidate refit fail, so that candidate was never tested.

Value

A list with class "hzr_bootstrap" containing:

replicates

Data frame with columns replicate, parameter, and estimate, one row per parameter per successful replicate.

summary

Data frame with columns parameter, n, pct, mean, sd, min, max, ci_lower, ci_upper, one row per parameter. In mode = "select", pct is the selection frequency and the other statistics are conditional on selection. A free parameter whose sd is 0, to within rounding, across two or more replicates draws a warning naming it: the replicates did not re-estimate it. Parameters the fit holds fixed are exempt: those held by hzr_phase(fixed = ) and shapes a constraint derives. A conserved log_mu is not exempt: it is re-solved in every replicate.

n_success

Number of successfully converged replicates.

n_failed

Number of replicates that failed: the refit stopped with an error, returned something other than a fit, returned a non-finite objective or one at the optimizer's -1e10 sentinel (which stands in for a likelihood that could not be evaluated), did not converge, returned a finite objective but no parameter estimates, or reported convergence on a likelihood that is flat or still rising past its estimate. A single-distribution refit that ends where the likelihood is not defined reports no objective, so it fails as a non-finite objective, not at the sentinel (#486).

failure_reasons

Named integer vector counting why replicates failed, most common first: the refit's error message (or "error with an empty message"), "refit returned a <class>, not a fit object" ("an <class>" when the class begins with a vowel; or "... with no \code{fit}, not a fit object"), "refit returned no parameter estimates", or "non-finite objective (did not converge)" (the reason for a single-distribution refit that ends where the likelihood is not defined, #486), or "objective at the optimizer's -1e10 sentinel (no log-likelihood)" (a refit whose reported objective is still the sentinel), or "refit did not converge (converged = FALSE)" (a refit that reports converged = FALSE, including one stopped on a score that is not finite, #518), or "likelihood flat or still rising (no unique finite maximum)" (a refit that reports convergence where stepping its estimates further out, along a coefficient or the direction its standard errors or Hessian say is flattest, does not lower the log-likelihood. Under separation in the resample an estimate is at infinity; with collinear columns the estimates are finite but not determined. Such replicates also draw a warning naming their count, #654), or, in select mode, "base refit did not converge" (a base refit with a non-finite objective) or "base refit did not converge (converged = FALSE)" (a base refit that reports converged = FALSE: the screen measures its candidates against the base, so the replicate is not screened). It sums to n_failed, and is an empty named integer vector, never NULL, when none failed. When every replicate fails, hzr_bootstrap() also warns, naming the most common reason.

n_uncomputable_replicates

Select mode only: number of otherwise successful replicates whose screen stopped because no remaining candidate could be tested (its score statistic, or for a removal its Wald statistic, could not be computed), rather than because no candidate met slentry or slstay. A non-zero count means every reported selection frequency is biased: a candidate such a replicate could not test for entry counts as not selected, and one it could not test for removal as selected. A replicate that went on after deciding a variable without a Wald test (wald_no_variance) is not counted here unless it also stopped, and gets a warning of its own either way. Always 0 in refit mode.

uncomputable_reasons

Select mode only: named integer vector counting why candidate scores were unavailable, summed over every replicate. information_indefinite is the one to read first: it marks candidates whose effect is too large for the score test's approximation at zero, typically strong variables. Those are refit and Wald-tested automatically, so a candidate reaching this count is one whose refit also failed and which therefore went untested, understating its selection frequency. rows_differ (under criterion = "aic") and loglik_below_base (under any criterion) mark entries a replicate declined without comparing them. Such an entry counts as not selected unless a later step of the same replicate tested and entered it, so these can understate a selection frequency; a warning gives how many replicates completed after declining one. The tally counts attempts, not distinct candidates. Empty in refit mode.

n_nonmonotone_replicates

Select mode only: number of otherwise successful replicates in which a forward step lowered the log-likelihood. Entered models are nested, so this cannot occur at the optimum; such a replicate continued from a refit that did not converge, and its later selections are pooled on the same footing as any other. Always 0 in refit mode.

n_wald_fallback_replicates, n_wald_fallbacks

Select mode only: the number of otherwise successful replicates that entered at least one variable on a Wald test instead of the score statistic, and the total number of such entries across all replicates. The score criterion declines a candidate whose observed information is indefinite at beta = 0 (which happens when the effect is large), so those candidates are refit and Wald-tested rather than dropped. A high count means much of the selection was decided by a different criterion from the one requested, which matters most here: these entries drive the pooled selection frequencies on the same footing as every other. Always 0 in refit mode.

mode

"refit" (fixed-formula bootstrap) or "select" (embedded stepwise selection).

scope

Only present when mode == "select": the candidate scope used.

unresolved

Only present when mode == "select": the ⁠$scope$unresolved⁠ record of the up-front hzr_stepwise() screen on the original data, the names in force_in, force_out and a character scope that matched nothing and were ignored. print() shows a line when any is non-empty.

Selection mode is experimental

Everything reached through scope (the selection arguments, and the summary$pct selection frequencies they produce) is new and should be treated as unstable. The fixed-formula bootstrap (scope = NULL) is not affected and has been stable since 0.9.3.

Two reasons to expect change. The first is that the design is still being read off real runs rather than settled in advance: production screens have already moved the defaults once and turned up several ways a screen could report success while selecting nothing.

The second is scale, and it is the one to plan around. A screen over a large candidate pool runs for hours, and this function writes nothing until its final replicate, so a run that dies late loses everything. There is no built-in way to split one screen across processes and combine the parts. If you are running at that scale, drive hzr_bootstrap() in chunks from your own script and pool the replicates yourself: deriving each chunk's seed from its chunk number, offsetting replicate ids so a variable selected in two chunks is not counted once, and recomputing frequencies from the pooled replicates rather than averaging across chunks. Whatever eventually covers that inside the package may well change this function's interface.

See Also

hazard() for model fitting, vcov.hazard() for Hessian-based standard errors, hzr_stepwise() for the selection procedure used when scope is supplied.

Examples

data(avc)
avc <- na.omit(avc)
fit <- hazard(
  survival::Surv(int_dead, dead) ~ age + mal,
  data  = avc,
  dist  = "weibull",
  theta = c(mu = 0.01, nu = 0.5, 0, 0),
  fit   = TRUE
)
bs <- hzr_bootstrap(fit, n_boot = 50, seed = 123)
print(bs)

# Embedded stepwise selection: screen candidate covariates for how
# often they enter the model across resamples (R equivalent of SAS
# %HAZBOOT).
base <- hazard(
  survival::Surv(int_dead, dead) ~ 1,
  data  = avc,
  dist  = "weibull",
  theta = c(mu = 0.01, nu = 0.5),
  fit   = TRUE
)
bs_sel <- hzr_bootstrap(base, n_boot = 20, seed = 123,
                         scope = ~ age + mal,
                         slentry = 0.3, slstay = 0.2)
print(bs_sel)

Calibrate a continuous variable against an outcome

Description

Group a continuous covariate into quantile bins, compute the event probability (or hazard rate) per bin, and apply a link transform (logit, Gompertz, or Cox). This is the R equivalent of the SAS logit.sas and logitgr.sas macros.

Usage

hzr_calibrate(
  x,
  event,
  groups = 10L,
  by = NULL,
  link = c("logit", "gompertz", "cox"),
  time = NULL
)

Arguments

Details

Use this function before model entry to assess whether the covariate relationship with the outcome is approximately linear on the link scale. If the transformed probabilities are roughly linear against the group means, the covariate can enter the model untransformed. Curvature suggests a transformation (log, quadratic) may improve fit.

Value

A data frame with one row per group (or per group-by-stratum combination) and columns:

group

Integer group label.

by

Stratum level (only present when by is provided).

n

Number of observations in the group.

events

Number of events.

mean

Mean of x within the group.

min

Minimum of x within the group.

max

Maximum of x within the group.

prob

Event probability (events / n), or for Cox link: events / sum(time).

link_value

Transformed probability on the chosen link scale.

See Also

hzr_deciles() for model-based calibration after fitting.

Examples

data(avc)
avc <- na.omit(avc)

# Logit calibration of age
cal <- hzr_calibrate(avc$age, avc$dead, groups = 10)
print(cal)


if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(ggplot2)
  ggplot(cal, aes(mean, link_value)) +
    geom_point(size = 3) +
    geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
                linetype = "dashed") +
    labs(x = "Age at repair (months)", y = "Logit(P(death))") +
    theme_minimal()
}

Clamp probabilities away from 0 and 1

Description

Bounds a probability vector to [\epsilon,\, 1 - \epsilon] so that downstream \log p or \log(1 - p) evaluations stay finite. Used throughout the likelihood and prediction code to guard survival and CDF values against exact 0 or 1.

Usage

hzr_clamp_prob(p, eps = 1e-12)

Arguments

Value

Numeric vector bounded to ⁠[eps, 1 - eps]⁠.

See Also

hzr_log1pexp() and hzr_log1mexp() for the companion stable-logarithm primitives.

Examples

hzr_clamp_prob(c(0, 0.5, 1))

Competing risks cumulative incidence

Description

Compute cumulative incidence functions for multiple competing event types using the Aalen-Johansen estimator with Greenwood variance. This is the R equivalent of the SAS markov.sas macro.

Usage

hzr_competing_risks(time, event)

## S3 method for class 'hzr_competing_risks'
print(x, digits = 4, ...)

Arguments

Details

Unlike the naive 1 - KM estimator (which overestimates incidence when competing risks exist), this provides the correct marginal cumulative incidence for each event type.

Value

A data frame with one row per unique event time and columns:

time

Event time.

n_risk

Number at risk.

n_event_1, n_event_2, ...

Events of each type at this time.

n_censor

Number censored at this time.

surv

Overall event-free survival (freedom from all events).

incid_1, incid_2, ...

Cumulative incidence for each event type.

se_surv

Standard error of overall survival.

se_1, se_2, ...

Standard error of each cumulative incidence.

Note

This estimator is unweighted: every observation contributes a unit count to the at-risk and event tallies. There is no weights argument, so case or inverse-probability weights are not yet supported for competing-risks incidence (unlike hazard(), which accepts weights). Pre-expand weighted rows to individual records if an approximate weighted estimate is needed.

References

Aalen O, Johansen S (1978). An empirical transition matrix for non-homogeneous Markov chains based on censored observations. Scand J Statist 5(3):141–150.

Kalbfleisch JD, Prentice RL (1980). The Statistical Analysis of Failure Time Data. Wiley, New York.

See Also

hzr_kaplan() for single-event survival estimation.

Examples

data(valves)
valves_cc <- na.omit(valves)
# Combine death and PVE into a competing risks event variable
# 0 = censored, 1 = death, 2 = PVE
event_cr <- ifelse(valves_cc$dead == 1, 1L,
                   ifelse(valves_cc$pve == 1, 2L, 0L))
time_cr <- pmin(valves_cc$int_dead, valves_cc$int_pve)
cr <- hzr_competing_risks(time_cr, event_cr)
head(cr)

Decile-of-risk calibration

Description

Partition observations into groups (default 10) by predicted risk and compare observed vs. expected event counts in each group. Good calibration means the two track each other across the risk spectrum.

Usage

hzr_deciles(object, time, groups = 10L, status = NULL, event_time = NULL)

Arguments

Details

This reproduces the SAS deciles.hazard.sas macro. All subjects are ranked by predicted survival at the horizon time and split into equal-sized risk groups. Within each group the expected event count is the sum of each subject's predicted cumulative hazard at its own follow-up time, and the observed count is its number of events; under conservation of events the group totals sum to the total observed events. A subject with an entry time (time_lower on a left-truncated fit) is at risk only after it, so it contributes its cumulative hazard at exit minus its value at entry. For a weighted fit both counts carry the case weights, since such a fit conserves weighted events. Both follow hzr_gof(). The horizon therefore only stratifies subjects into risk groups; it does not restrict or exclude any subject, and the expected/observed totals are independent of it.

A multiphase fit that dropped the rows where a phase covariate was missing is refused with an error, as hzr_gof() refuses it: its design matrices hold fewer rows than there are subjects, so no subject's predicted cumulative hazard can be formed. Refit it on complete cases.

Value

A data frame with one row per risk group and columns:

group

Integer group label (1 = lowest risk, ranked by predicted survival at time).

n

Number of observations in the group.

events

Observed event count in the group (all events over follow-up), weighted by the case weights for a weighted fit.

expected

Expected event count: the sum of each subject's predicted cumulative hazard at its own follow-up time, minus its value at its entry time when it has one, weighted by the case weights for a weighted fit.

observed_rate

Observed event rate (events / n; for a weighted fit, events per unit of case weight).

expected_rate

Expected event rate (expected / n; for a weighted fit, per unit of case weight).

chi_sq

Chi-square contribution: (events - expected)^2 / expected. For a weighted fit the denominator is the Poisson variance of the weighted count, the sum of each subject's squared weight times its cumulative hazard, so the statistic does not change when every weight is multiplied by the same constant.

p_value

Upper-tail p-value from the chi-square test for this group (1 df).

mean_survival

Mean predicted survival probability at the horizon in the group.

mean_cumhaz

Mean predicted cumulative hazard at follow-up in the group: unweighted, and without the entry-time correction, so for a left-truncated or weighted fit it is not expected / n.

An attribute "overall" is attached with the overall chi-square statistic, degrees of freedom, and p-value.

See Also

predict.hazard() for the prediction types used internally.

Examples

data(avc)
avc <- na.omit(avc)
fit <- hazard(
  survival::Surv(int_dead, dead) ~ age + mal,
  data  = avc,
  dist  = "weibull",
  theta = c(mu = 0.01, nu = 0.5, beta_age = 0, beta_mal = 0),
  fit   = TRUE
)
cal <- hzr_deciles(fit, time = 120)
print(cal)

Generalized temporal decomposition

Description

Computes the cumulative distribution G(t), density g(t), and hazard h(t) = g(t)/(1 - G(t)) for the parametric family defined by half-life, time exponent, and shape. It supplies the "cdf" and "hazard" phase shapes of a multiphase hazard model; the "g3" late phase comes from hzr_decompos_g3(), and the "constant" phase is linear in time.

Usage

hzr_decompos(time, t_half, nu, m)

Arguments

Value

A named list with four numeric vectors, each the same length as time:

G

Cumulative distribution G(t) \in [0, 1].

g

Density g(t) = dG/dt \ge 0. The "early" phase temporal pattern.

h

Hazard h(t) = g(t)/(1 - G(t)) \ge 0. The "late" phase temporal pattern.

log_surv

\log(1 - G(t)), computed from each case's own log-scale terms rather than from G, so it keeps its accuracy where G rounds to 1. -log_surv is the cumulative hazard of a "hazard" phase, and h is computed as \exp(\log g - \log(1 - G)).

Parameter mapping from SAS/C HAZARD

The original C code used separate parameterizations for early (DELTA, RHO/THALF, NU, M) and late (TAU, GAMMA, ALPHA, ETA) phases. The early phase maps onto the three parameters here: DELTA must be 0, and RHO is fixed by THALF, NU and M. The late phase does not: it is a separate four-parameter shape computed by hzr_decompos_g3(). See hzr_argument_mapping() for the full translation table.

Valid parameter combinations

Six cases are defined by the signs of nu and m:

The combination m < 0 and nu < 0 is undefined and raises an error. nu = 0 is supported only with m < 0 (Case 2L, the exponential-decay limit); nu = 0 with m >= 0 has no usable limiting form and raises an error.

Mathematical form

The construction fixes a rate \rho so that G(t_{1/2}) = 0.5 exactly. For the base case (m > 0,\ \nu > 0):

\rho = \nu \, t_{1/2} \left(\frac{2^m - 1}{m}\right)^{\!\nu}, \qquad b(t) = \frac{\nu t}{\rho}

The CDF and density are then

G(t) = \bigl(1 + m \, b(t)^{-1/\nu}\bigr)^{-1/m}, \qquad g(t) = \bigl(1 + m \, b(t)^{-1/\nu}\bigr)^{-1/m - 1} \, b(t)^{-1/\nu - 1} / \rho

and the hazard is h(t) = g(t) / (1 - G(t)). The remaining five cases in the table above arise as limits (m \to 0, \nu \to 0) or sign reflections of this base form; the implementation dispatches to the appropriate branch after inspecting the signs of nu and m. See vignette("mf-mathematical-foundations") for the full derivation of every case.

References

Blackstone EH, Naftel DC, Turner ME Jr. The decomposition of time-varying hazard into phases, each incorporating a separate stream of concomitant information. J Am Stat Assoc. 1986;81(395):615–624. doi:10.1080/01621459.1986.10478314

Rajeswaran J, Blackstone EH, Ehrlinger J, Li L, Ishwaran H, Parides MK. Probability of atrial fibrillation after ablation: Using a parametric nonlinear temporal decomposition mixed effects model. Stat Methods Med Res. 2018;27(1):126–141. doi:10.1177/0962280215623583

See Also

hzr_phase_cumhaz() for the phase-level cumulative hazard contribution, hzr_argument_mapping() for SAS/C parameter mapping, hzr_phase() for specifying phases in hazard() models.

vignette("mf-mathematical-foundations") for the full derivation.

Examples

t_grid <- seq(0.1, 10, by = 0.1)

# Case 1: standard sigmoidal (m > 0, nu > 0)
d1 <- hzr_decompos(t_grid, t_half = 3, nu = 2, m = 1)
plot(t_grid, d1$G, type = "l", main = "CDF (m=1, nu=2)")

# Case 1L: Weibull-like (m = 0, nu > 0)
d1L <- hzr_decompos(t_grid, t_half = 3, nu = 2, m = 0)

# Case 2: heavy-tailed (m < 0, nu > 0)
d2 <- hzr_decompos(t_grid, t_half = 3, nu = 2, m = -1)

Late-phase (G3) temporal decomposition

Description

Computes the cumulative intensity G_3(t) and its derivative g_3(t) = dG_3/dt for the late-phase parametric family used in the original Blackstone C/SAS HAZARD code. Unlike hzr_decompos() (which computes the early-phase G1, a bounded CDF), this function can produce unbounded values, making it suitable for modelling increasing late risk.

Usage

hzr_decompos_g3(time, tau, gamma, alpha, eta)

Arguments

Value

A named list with two numeric vectors, each the same length as time:

G3

Cumulative intensity G_3(t) \ge 0 (may exceed 1).

g3

Derivative g_3(t) = dG_3/dt \ge 0.

Mathematical form

When \alpha > 0:

G_3(t) = \bigl(\bigl((t/\tau)^\gamma + 1\bigr)^{1/\alpha} - 1\bigr)^\eta

When \alpha = 0 (limiting exponential case):

G_3(t) = \bigl(\exp\bigl((t/\tau)^\gamma\bigr) - 1\bigr)^\eta

Parameter mapping from SAS/C HAZARD

References

Blackstone EH, Naftel DC, Turner ME Jr. The decomposition of time-varying hazard into phases, each incorporating a separate stream of concomitant information. J Am Stat Assoc. 1986;81(395):615–624. doi:10.1080/01621459.1986.10478314

See Also

hzr_decompos() for the early-phase (G1) decomposition, hzr_phase_cumhaz() for phase-level cumulative hazard helpers.

Examples

t_grid <- seq(0.1, 10, by = 0.1)

# Weibull-like: alpha = 1 gives G3(t) = (t/tau)^(gamma*eta)
d <- hzr_decompos_g3(t_grid, tau = 1, gamma = 3, alpha = 1, eta = 1)
plot(t_grid, d$G3, type = "l", main = "G3: power law (gamma=3)")

# General case with alpha > 0
d2 <- hzr_decompos_g3(t_grid, tau = 2, gamma = 2, alpha = 0.5, eta = 1)

Evaluate a hazard model at parameters you supply

Description

Computes a model's objective, and optionally its hazard and cumulative hazard, at parameters you supply rather than at parameters fitted from the data. The objective is the log-likelihood, except for a model built with objective = "sas" on data with interval-censored rows, where it is PROC HAZARD's interval-mean-hazard objective (see objective under Value). This is what a parity check needs: the objective at another program's converged estimates, evaluated by this package's own code.

Usage

hzr_evaluate(object, theta, times = NULL)

Arguments

Details

The result is not a fit and does not pretend to be one. It carries no standard errors, no convergence status and no covariance matrix, because nothing was estimated: print() says so on its first line, and the parameters are labelled as supplied. Where a fitted model is what you want, use hazard(..., fit = TRUE).

A phase built with hzr_phase(constraint = ) derives one of its shapes from the others, and that rule is applied to theta here as the fit applies it, so the derived entry you pass is replaced rather than used as given. At a fitted model's own estimates this returns that fit's objective, under Conservation of Events too: since #362 the fit's objective is recomputed at the estimates it returns, except where the fit warns that it could not, and then the two can differ.

A theta that passes the input checks but that the likelihood cannot evaluate – an overflowing rate, a shape outside the family, or a value past a guard that stops short of where the log-likelihood itself overflows – gives -Inf, with a warning of class "hzr_evaluate_not_finite", for every distribution. A Weibull mu or nu at or below 0 is still refused outright by those input checks.

Value

An object of class hzr_evaluation: a list with theta, the parameters the likelihood was evaluated at – the vector supplied, with any constrained entry replaced by the value its phase derives, which is warned about as hazard() warns; logLik, the log-likelihood there, or, for a model built with objective = "sas", that objective, which is not a log-likelihood wherever the data have interval-censored rows; objective, "likelihood" or "sas", saying which of the two logLik holds; dist; n_obs and n_events, the rows the likelihood scored and the exact events among them (a left- or interval-censored row counts in n_obs, not in n_events); and, when times was given, curve, a data frame of time, hazard and cumulative_hazard, which is NA throughout where logLik is -Inf.

See Also

hazard() to fit a model, hzr_theta_names() for the parameter order.

Examples

data(avc, package = "TemporalHazard")
spec <- hazard(survival::Surv(int_dead, dead) ~ 1, data = avc,
               dist = "weibull", theta = c(1, 1), fit = FALSE)
# The likelihood at parameters from somewhere else, e.g. a SAS run.
hzr_evaluate(spec, theta = c(0.05, 0.9))

Goodness-of-fit: observed vs. predicted events

Description

Compare a fitted hazard model with the data two ways: its survival curve against the nonparametric Kaplan-Meier estimate, and the number of events it expects against the number observed, tallied over follow-up. This is the R equivalent of the SAS hazplot.sas macro and implements the conservation-of-events diagnostic.

Usage

hzr_gof(object, time_grid = NULL)

Arguments

Details

The diagnostic is for right-censored data: every stored status must be 0 (censored) or 1 (event). A fit with any left-censored (status -1) or interval-censored (status 2) row is refused with an error. So is a multiphase fit that dropped the rows where a covariate was missing: its design matrices then hold fewer rows than there are patients, so no patient-by-patient tally can be formed. Refit it on complete cases.

At each time point the function computes:

• The Kaplan-Meier survival and cumulative hazard.

• The parametric survival and cumulative hazard from the fitted model at the covariate means (and per-phase components for multiphase models). This is the curve to plot against the Kaplan-Meier estimate.

• Cumulative observed events vs. cumulative expected events. Each patient's expected count is their own cumulative hazard, from their own covariates, at the end of their follow-up, less their cumulative hazard at entry when the fit is left truncated (time_lower on a status 0 or 1 row). These are summed over the patients leaving follow-up at each time.

• The running residual (expected minus observed).

The conservation-of-events principle says a model fit by maximum likelihood predicts as many events as were observed: add up every patient's cumulative hazard over their follow-up and you get the event count back. The final residual is then zero and the printed "Conservation ratio (E/O)" is 1. A multiphase fit with Conservation of Events applied (control = list(conserve = TRUE), the default) meets the identity by construction. Weibull and exponential fits meet it at the exact maximum, so at a converged fit E/O sits close to 1, off only by how far short of the maximum the optimizer stopped. The log-logistic and log-normal models carry no such identity, and for them E/O is a check of calibration in total.

The parametric curve is a different quantity. For a model with covariates it belongs to one "mean patient" with average covariates, and the mean patient's cumulative hazard is not the average of the patients' cumulative hazards, so par_cumhaz does not enter the expected count. For an intercept-only model every patient shares that curve. Without entry times each patient's expected count is the curve at their exit time; with entry times it is the curve's rise from entry to exit, so the two differ. The means are those of the design-matrix columns, taken phase by phase for a multiphase fit, whether its covariates enter globally, through the phase formulas or both, so a factor enters as the proportion of patients in each level. With time_windows, a phase built on the global covariates carries their means in the window that contains each time; a phase formula's own columns keep their plain means.

For a weighted fit both tallies carry the case weights: observed events are \sum_i w_i d_i and expected events \sum_i w_i H_i, the form in which a weighted fit conserves events. The n_risk, n_event, n_censor and Kaplan-Meier columns are unweighted.

Each patient is placed at the Kaplan-Meier time survival::survfit() gives them, which merges exit times closer together than its tolerance. On the default grid every patient lands on a grid point, and hzr_gof() warns if one cannot be placed. With a custom time_grid, a patient is counted in both tallies only if that time, or failing it their own follow-up time, falls on a grid point. A grid such as seq() over the follow-up holds few exit times, so hzr_gof() warns with the number left out, and the totals, residual and E/O then describe only the patients counted. Include the exit times in time_grid to count everyone.

Value

A data frame with one row per time point and columns:

time

Evaluation time.

n_risk

Number at risk at this time: subjects still in follow-up. A subject with an entry time is at risk only after it.

n_event

Number of events at this time.

n_censor

Number censored at this time.

km_surv

Kaplan-Meier survival estimate, using the counting-process risk set when the fit has entry times.

km_cumhaz

Kaplan-Meier cumulative hazard (-\log(\text{km\_surv})).

par_surv

Parametric survival from the fitted model, at the covariate means for a model with covariates. For plotting against km_surv; not used for cum_expected.

par_cumhaz

Parametric cumulative hazard, at the covariate means for a model with covariates. For plotting; not used for cum_expected.

cum_observed

Cumulative observed events to this time, weighted by the case weights for a weighted fit. With a custom time_grid, only the patients whose exit is a grid point up to this time are counted (see Details).

cum_expected

Cumulative expected events: over the patients leaving follow-up by this time (with a custom time_grid, the same patients as cum_observed), the sum of each patient's own cumulative hazard at exit minus that at entry, weighted by the case weights for a weighted fit. With time_windows, both cumulative hazards use the patient's covariate window at exit, as the likelihood does.

residual

Expected minus observed (cum_expected - cum_observed).

For multiphase models, additional columns are appended for each phase: par_cumhaz_<phase>, also at the covariate means.

An attribute "summary" is attached with scalar diagnostics: total observed events, total expected events, the final residual, and n_tallied, the number of patients the two tallies cover, out of n.

See Also

hzr_deciles() for decile-of-risk calibration, predict.hazard() for prediction types.

Examples

data(avc)
avc <- na.omit(avc)
fit <- hazard(
  survival::Surv(int_dead, dead) ~ age + mal,
  data  = avc,
  dist  = "weibull",
  theta = c(mu = 0.01, nu = 0.5, beta_age = 0, beta_mal = 0),
  fit   = TRUE
)
gof <- hzr_gof(fit)
print(gof)

# Expected events are summed per patient.  This fit has no entry times,
# so the total is the sum of each patient's own cumulative hazard at
# their follow-up time:
nd <- avc[, c("age", "mal")]
nd$time <- avc$int_dead
c(hzr_gof = attr(gof, "summary")$total_expected,
  predict = sum(predict(fit, newdata = nd, type = "cumulative_hazard")))

# Plot observed vs expected events
if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(ggplot2)
  ggplot(gof, aes(x = time)) +
    geom_line(aes(y = cum_observed), colour = "#D55E00") +
    geom_line(aes(y = cum_expected), colour = "#0072B2") +
    labs(x = "Time", y = "Cumulative events") +
    theme_minimal()
}

Kaplan-Meier survival with exact logit confidence limits

Description

Compute the product-limit (Kaplan-Meier) survival estimate with logit-transformed confidence limits that respect the [0, 1] boundary. This is the R equivalent of the SAS kaplan.sas macro.

Usage

hzr_kaplan(time, status, conf_level = 0.95, event_only = TRUE)

Arguments

Details

The standard Greenwood confidence interval can exceed [0, 1] in the tails. The logit-transformed interval avoids this by working on the log-odds scale:

\text{CL}_{\text{lower}} = S / \bigl(S + (1-S)\, \exp(z_\alpha\,\text{SI})\bigr)

\text{CL}_{\text{upper}} = S / \bigl(S + (1-S)\, \exp(-z_\alpha\,\text{SI})\bigr)

where \text{SI} = \sqrt{V_P - 1} / (1 - S), V_P is the cumulative Greenwood variance product, and z_\alpha is the normal quantile for the requested confidence level.

Value

A data frame with one row per time point and columns:

time

Event/censoring time.

n_risk

Number at risk at start of interval.

n_event

Number of events at this time.

n_censor

Number censored at this time.

survival

Kaplan-Meier survival estimate.

std_err

Standard error of survival (Greenwood).

cl_lower

Lower confidence limit (logit-transformed).

cl_upper

Upper confidence limit (logit-transformed).

cumhaz

Cumulative hazard = -\log(S).

hazard

Interval hazard rate = \log(S_{t-1} / S_t) / \Delta t.

density

Probability density estimate = (S_{t-1} - S_t) / \Delta t.

life

Restricted mean survival time (area under curve to this time).

References

Kaplan EL, Meier P (1958). Nonparametric estimation from incomplete observations. J Am Stat Assoc 53(282):457–481. doi:10.1080/01621459.1958.10501452

Greenwood M (1926). The natural duration of cancer. Reports on Public Health and Medical Subjects 33:1–26.

See Also

hzr_gof() for parametric vs. nonparametric comparison.

Examples

data(cabgkul)
km <- hzr_kaplan(cabgkul$int_dead, cabgkul$dead)
head(km)


if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(ggplot2)
  ggplot(km, aes(time)) +
    geom_step(aes(y = survival * 100)) +
    geom_ribbon(aes(ymin = cl_lower * 100, ymax = cl_upper * 100),
                stat = "identity", alpha = 0.2) +
    labs(x = "Months", y = "Survival (%)") +
    theme_minimal()
}

Numerically stable log(1 - exp(-x)) for x > 0

Description

Computes

\mathrm{log1mexp}(x) = \log\!\bigl(1 - e^{-x}\bigr), \qquad x > 0

accurately across the full range. Following Mächler (2012), the evaluation switches at x = \log 2: for small x it uses \log(-\mathrm{expm1}(-x)) (avoiding cancellation as x \to 0), and for larger x it uses \mathrm{log1p}(-e^{-x}) (avoiding loss of precision as x \to \infty). Non-positive or non-finite inputs return NA.

Usage

hzr_log1mexp(x)

Arguments

Value

Numeric vector with element-wise \log(1 - e^{-x}).

References

Mächler M (2012). Accurately Computing \log(1 - e^{-|a|}). R package Rmpfr vignette. https://CRAN.R-project.org/package=Rmpfr/vignettes/log1mexp-note.pdf

See Also

hzr_log1pexp() for the softplus \log(1 + e^{x}), hzr_clamp_prob() for boundary-safe probabilities.

Examples

hzr_log1mexp(c(0.01, 0.5, 5))

Numerically stable log(1 + exp(x))

Description

Computes the "softplus" function

\mathrm{log1pexp}(x) = \log\!\bigl(1 + e^{x}\bigr)

without overflow. A naive log(1 + exp(x)) overflows to Inf for x \gtrsim 710 even though the true value is \approx x. The identity \log(1 + e^{x}) = x + \log(1 + e^{-x}) for x > 0 keeps every intermediate finite.

Usage

hzr_log1pexp(x)

Arguments

Value

Numeric vector with element-wise \log(1 + e^{x}).

References

Mächler M (2012). Accurately Computing \log(1 - e^{-|a|}). R package Rmpfr vignette. https://CRAN.R-project.org/package=Rmpfr/vignettes/log1mexp-note.pdf

See Also

hzr_log1mexp() for the complementary \log(1 - e^{-x}), hzr_clamp_prob() for boundary-safe probabilities.

Examples

hzr_log1pexp(c(-50, 0, 50))

Wayne Nelson cumulative hazard estimator with lognormal confidence limits

Description

Compute the Nelson-Aalen cumulative hazard estimate with lognormal confidence limits. Supports weighted events for severity-adjusted analyses of repeated/recurrent events. This is the R equivalent of the SAS nelsonl.sas macro.

Usage

hzr_nelson(time, event, weight = NULL, conf_level = 0.95)

## S3 method for class 'hzr_nelson'
print(x, digits = 4, ...)

Arguments

Details

Unlike survival::survfit() which uses the Breslow estimator with Greenwood variance, this function uses the Wayne Nelson estimator with lognormal confidence limits that are always non-negative.

Value

A data frame with one row per unique event time and columns:

time

Event time.

n_risk

Number at risk.

n_event

Number of events at this time.

weight_sum

Sum of event weights at this time.

cumhaz

Nelson cumulative hazard estimate.

std_err

Standard error.

cl_lower

Lower lognormal confidence limit.

cl_upper

Upper lognormal confidence limit.

hazard

Interval hazard rate.

cum_events

Cumulative (weighted) event count.

References

Nelson W (1972). Theory and applications of hazard plotting for censored failure data. Technometrics 14(4):945–966. doi:10.1080/00401706.1972.10488991

Aalen O (1978). Nonparametric inference for a family of counting processes. Ann Statist 6(4):701–726. doi:10.1214/aos/1176344247

See Also

hzr_kaplan() for survival estimation.

Examples

data(cabgkul)
nel <- hzr_nelson(cabgkul$int_dead, cabgkul$dead)
head(nel)

Specify a single hazard phase

Description

Creates an hzr_phase object describing one term in a multiphase additive cumulative hazard model. Pass a list of these to the phases argument of hazard() when dist = "multiphase".

Usage

hzr_phase(
  type = c("cdf", "hazard", "constant", "g3"),
  t_half = 1,
  nu = 1,
  m = 0,
  tau = 1,
  gamma = 1,
  alpha = 1,
  eta = 1,
  formula = NULL,
  fixed = character(0),
  constraint = c("none", "alpha_gamma_eta", "eta_gamma")
)

## S3 method for class 'hzr_phase'
print(x, ...)

Arguments

Value

An S3 object of class "hzr_phase" with elements:

type

Phase type string.

t_half

Initial half-life (cdf/hazard phases).

nu

Initial time exponent (cdf/hazard phases).

m

Initial shape exponent (cdf/hazard phases).

tau

Scale parameter (g3 phases).

gamma

Time exponent (g3 phases).

alpha

Shape parameter (g3 phases).

eta

Outer exponent (g3 phases).

formula

Phase-specific formula or NULL.

fixed

Character vector of fixed parameter names (may be empty).

constraint

The shape constraint (g3 phases).

Role in the multiphase model

Each phase is one term j in the additive cumulative hazard

H(t \mid \mathbf{x}) = \sum_{j=1}^{J} \mu_j(\mathbf{x}) \, \Phi_j(t)

where \mu_j(\mathbf{x}) = \exp(\alpha_j + \mathbf{x}_j^\top \beta_j) is the phase-specific log-linear scale and \Phi_j(t) is the temporal shape selected by type (below). The t_half/nu/m (or g3 tau/gamma/alpha/eta) arguments set the starting values for that shape; formula attaches the covariates \mathbf{x}_j that enter \mu_j.

Phase types

The type argument chooses the temporal shape \Phi_j(t) for the phase. Each captures a qualitatively different pattern of risk over time; a typical clinical model combines an early, a constant, and a late phase so that the total hazard can fall, level off, and rise again.

"cdf": early, resolving risk

Named for the cumulative distribution function: the phase contributes \Phi(t) = G(t), the bounded CDF of the temporal decomposition (0 at t = 0, rising to a ceiling of 1). Because it saturates, the hazard it adds, \mu\,g(t), peaks early and then decays toward zero (the signature of a one-time insult that patients either succumb to or survive past, e.g. peri-operative mortality). Shape set by t_half, nu, m. SAS/C equivalent: the Early (G1) phase.

"hazard": accumulating aging risk (G1 family)

Named because the phase contributes a cumulative hazard built from the same G1 family: \Phi(t) = -\log(1 - G(t)), which is unbounded and monotone increasing. Its hazard \mu\,h(t) rises without leveling off, so it models risk that grows as subjects age. This is an alternative late-risk form derived from G1; for the original SAS/C late phase prefer "g3". Shape set by t_half, nu, m.

"g3": late, rising risk (original C/SAS late phase)

Named for the G3 (third) decomposition family used by the original HAZARD program for the late phase. It contributes \Phi(t) = G_3(t) from hzr_decompos_g3(), an unbounded intensity with its own four-parameter shape (tau, gamma, alpha, eta) that is more flexible than the G1-derived "hazard" form for capturing accelerating late mortality (e.g. structural valve deterioration years after surgery). Use this when reproducing classic three-phase HAZARD models. SAS/C equivalent: the Late (G3) phase.

"constant": flat background rate

A time-invariant hazard: \Phi(t) = t, so the added hazard \mu is constant (the exponential model). It represents the steady, ongoing risk present at all follow-up times, independent of how long ago the time origin was. Takes no shape parameters; only its scale \mu (and any covariates) is estimated. SAS/C equivalent: the Constant (G2) phase.

The shape derivative \varphi_j = d\Phi_j/dt (which forms the instantaneous hazard contribution \mu_j\,\varphi_j(t)) is g(t) for "cdf", h(t) for "hazard", g_3(t) for "g3", and 1 for "constant".

See Also

hazard() for fitting multiphase models, hzr_decompos() for the underlying parametric family, hzr_phase_cumhaz() and hzr_phase_hazard() for computing \Phi(t) and \phi(t) from these specifications.

vignette("fitting-hazard-models") for multiphase fitting examples, vignette("mf-mathematical-foundations") for the mathematical framework.

Examples

# Classic 3-phase Blackstone pattern
early <- hzr_phase("cdf",      t_half = 0.5, nu = 2, m = 0)
const <- hzr_phase("constant")
late  <- hzr_phase("g3", tau = 1, gamma = 3, alpha = 1, eta = 1)

# Fix all shapes (C/SAS-style: only estimate mu)
early_fixed <- hzr_phase("cdf", t_half = 0.5, nu = 2, m = 0,
                          fixed = "shapes")
late_fixed  <- hzr_phase("g3", tau = 1, gamma = 3, alpha = 1, eta = 1,
                          fixed = "shapes")

# Derive alpha from gamma and eta (SAS/C FIXGAE2)
late_gae2 <- hzr_phase("g3", tau = 14, gamma = 22, eta = 0.18,
                        constraint = "alpha_gamma_eta")

# Fix only some parameters
early_partial <- hzr_phase("cdf", t_half = 0.5, nu = 2, m = 0,
                            fixed = c("nu", "m"))

# Phase with specific covariates
early_cov <- hzr_phase("cdf", t_half = 0.5, nu = 2, m = 0,
                        formula = ~ age + shock)

# Use in hazard():
# hazard(Surv(time, status) ~ age, data = dat,
#        dist = "multiphase",
#        phases = list(early = early, constant = const, late = late))

Cumulative hazard contribution from a single phase

Description

Computes \Phi_j(t) for one phase in the additive model H(t|x) = \sum_j \mu_j(x) \Phi_j(t).

Usage

hzr_phase_cumhaz(
  time,
  t_half = 1,
  nu = 1,
  m = 0,
  type = c("cdf", "hazard", "constant")
)

Arguments

Details

• "cdf": \Phi(t) = G(t). Bounded [0, 1]. Models early risk that resolves over time.

• "hazard": \Phi(t) = -\log(1 - G(t)). Monotone increasing. Models late or aging risk. This is the cumulative hazard derived from the hazard function h(t), since \int_0^t h(s)\,ds = -\log(1 - G(t)).

• "constant": \Phi(t) = t. Ignores t_half, nu, m. Equivalent to exponential (constant hazard rate).

Value

Numeric vector of cumulative hazard contributions \Phi(t), same length as time.

See Also

hzr_decompos() for the underlying parametric family, hzr_phase_hazard() for the instantaneous hazard contribution.

Examples

t_grid <- seq(0.1, 10, by = 0.1)
phi_early <- hzr_phase_cumhaz(t_grid, t_half = 2, nu = 2, m = 0,
                               type = "cdf")
phi_late  <- hzr_phase_cumhaz(t_grid, t_half = 5, nu = 1, m = 0,
                               type = "hazard")
phi_const <- hzr_phase_cumhaz(t_grid, type = "constant")

Instantaneous hazard contribution from a single phase

Description

Computes \phi_j(t) = d\Phi_j/dt for one phase, the derivative of the cumulative hazard contribution returned by hzr_phase_cumhaz().

Usage

hzr_phase_hazard(
  time,
  t_half = 1,
  nu = 1,
  m = 0,
  type = c("cdf", "hazard", "constant")
)

Arguments

Details

• "cdf": \phi(t) = g(t) (density).

• "hazard": \phi(t) = h(t) = g(t)/(1-G(t)).

• "constant": \phi(t) = 1.

Value

Numeric vector of instantaneous hazard contributions \phi(t), same length as time.

See Also

hzr_decompos() for the underlying parametric family, hzr_phase_cumhaz() for the cumulative version.

Examples

t_grid <- seq(0.1, 10, by = 0.1)
phi_early <- hzr_phase_hazard(t_grid, t_half = 2, nu = 2, m = 0,
                               type = "cdf")
phi_late  <- hzr_phase_hazard(t_grid, t_half = 5, nu = 1, m = 0,
                               type = "hazard")

Read a SAS outhaz estimate dataset

Description

⁠PROC HAZARD⁠'s ⁠outhaz=⁠ dataset stores the converged estimates and the asymptotic variance-covariance matrix at full double precision, where the printed .lst carries about seven significant figures. For any quantity the dataset holds, it is the better parity reference: print precision stops being the binding constraint and optimizer convergence takes over.

Usage

hzr_read_outhaz(path)

Arguments

Details

The log-likelihood is not stored here; take it from the .lst.

Value

An hzr_outhaz object: a list with estimates (named numeric), status (named integer, 1 free / 0 fixed), vcov (matrix over free parameters, dimnames set) and flags (named numeric of model-structure flags). When no parameter is free, or the dataset was written with NOCOV, vcov is a 0x0 matrix rather than NULL, so check its dimensions rather than is.null().

Experimental

This function is experimental and its return shape is expected to change. The result carries the fitted model, so predict.hzr_outhaz() predicts from it (that is what the hzr_translate_sas(librefs = ) path emits), but it is not a hazard object and none of the other hazard methods apply to it. The ⁠_STATUS_⁠ coding is asserted against a synthetic fixture, so a real ⁠OUTHAZ=⁠ file using a different convention would yield an empty vcov alongside a fully populated estimates.

Examples

f <- system.file("extdata", "outhaz-fixture.rds", package = "TemporalHazard")
if (nzchar(f)) str(hzr_read_outhaz(f))

Build repeated-event segments, reproducing the SAS ⁠%repeat⁠ macro

Description

Converts a long data set of candidate event times into one row per inter-event segment, ready to fit as a repeated-events model. This is a native R implementation of the SAS macro ⁠%repeat⁠ used by the repeated-events HAZARD jobs, whose fit input was seldom saved.

Usage

hzr_repeated_events(data, id, time, followup, indicator)

Arguments

Details

The input holds one row per candidate event per subject, with an indicator naming the event of interest: a value of 1 marks an event, and 0 or NA marks its absence. For each subject the function keeps the event rows, guarantees at least one row, appends a censored row at the end of follow-up when the last event happened earlier, and then lags the event time within subject so that each row describes the interval since the previous event.

Two details are inherited from the macro and are not R conventions. An indicator code that is neither 0, 1 nor missing is treated as neither an event nor an absence, and rows are ordered with missing times first, as SAS sorts them.

The result can contain segments of zero length, where iv_start is not less than time: for example, an event that shares its time with the subject's previous row. They are kept, as the macro keeps them. A model that treats iv_start as the time a subject enters the risk set needs each segment to end after it starts, so adjust or remove these rows before such a fit. The SAS jobs this function reproduces moved time forward by about one hour on each.

data may not already contain a column named rcensor, first, last, event, event_no, iv_start, iv_seg or renewal: this function creates all eight and would silently overwrite an existing one of the same name.

Value

A data frame with one row per inter-event segment. Every input column is retained, with time and indicator altered where the macro alters them, and the following columns added:

event

1 when the row is an event of interest, otherwise 0.

event_no

Running count of event repeats within the subject.

rcensor

1 when the row is right censored at the end of follow-up, otherwise 0. This is set both for a genuine censoring row and for an event row whose event time equals the end of follow-up. rcensor and event are therefore not mutually exclusive; see the Note below.

iv_start

Interval from time zero to the start of the segment.

iv_seg

Duration of the segment, time minus iv_start.

renewal

Segment number under the modulated renewal formulation.

first, last

1 when the row was the subject's first or last row at stage 4, otherwise 0. The row stage 4 appends to pad a subject's censored tail inherits its source row's flags, so a subject can have two rows with last == 1 together, or two with first == 1 together.

Note

rcensor is not a plain censoring indicator: a row can have event == 1 and rcensor == 1 together, when a subject's last event happens to fall exactly at the end of follow-up. This is faithful to the SAS macro, whose own header warns that this combination is not compatible with a HAZARD-style fit. Do not pass rcensor directly as a censoring indicator without first accounting for the overlap with event.

Examples

# One row per candidate event per subject: s1 has two events, s2 one,
# and s3 none (its only row is a non-event, so it contributes no rows
# of its own before the function pads it with a censored one below).
events <- data.frame(
  id = c("s1", "s1", "s2", "s3"),
  t = c(1, 3, 2, NA),
  fu = c(10, 10, 10, 10),
  ev = c(1, 1, 1, 0)
)
# The result appends a censored row at end of follow-up for any subject
# whose last event happened before `fu`, so every subject gets at least
# one row even when, like s3, they never had an event.
hzr_repeated_events(events, id = "id", time = "t", followup = "fu", indicator = "ev")

Stepwise covariate selection for a parametric hazard model

Description

Run forward, backward, or two-way stepwise selection on an existing hazard fit using score (Q) statistics, Wald p-values, or AIC deltas as the entry / retention criterion. Phase-specific entry is supported for multiphase models: a covariate can enter one phase and not another.

Usage

hzr_stepwise(
  fit,
  scope = NULL,
  data,
  direction = c("both", "forward", "backward"),
  criterion = c("score", "wald", "aic"),
  slentry = 0.3,
  slstay = 0.2,
  max_steps = 50L,
  max_move = 4L,
  force_in = character(),
  force_out = character(),
  trace = TRUE,
  ...
)

## S3 method for class 'hzr_stepwise'
print(x, ...)

## S3 method for class 'hzr_stepwise'
summary(object, ...)

## S3 method for class 'summary.hzr_stepwise'
print(x, ...)

## S3 method for class 'hzr_stepwise'
as.data.frame(x, ...)

Arguments

Details

The steps data frame has columns:

step_num

Integer sequence starting at 1.

action

"enter", "drop", or "frozen".

variable

The term entered, dropped or frozen, as the model's terms() label: a non-syntactic column ⁠_X1⁠ is recorded as `_X1` on every row, whether the scope named it "_X1", "`_X1`" or in a formula.

phase

Phase name (multiphase) or NA_character_.

criterion

The criterion actually applied to this step: "score", "wald", or "aic". Under criterion = "score" the drop rows read "wald", because score is entry-only.

score

Winning score used for the decision.

stat, df

Test statistic and degrees of freedom.

stat_type

What stat is on this row, and so which reference distribution recomputes its p-value: "score_q" (chi-square on df), "wald_z" (standard normal) or "wald_chisq" (chi-square on df). df alone does not distinguish them; a scalar Wald is reported as a z, not as its square, so it and a score Q are both recorded at df = 1 while calling for different distributions. It also identifies the rows the Wald fallback rescued, but only among entry rows: under criterion = "score" those are the rows with action == "enter" & stat_type == "wald_z". Drop rows are always Wald-tested under that criterion (removal follows SAS and is tested on the current model's Wald p-value), so they read "wald_z" whether or not the fallback ever fired. See $criteria$n_wald_fallbacks for the count.

p_value, delta_aic

Always populated when computable, regardless of the active criterion.

logLik, aic, n_coef

Goodness-of-fit diagnostics of the model after this step. Under objective = "sas", on a step whose model reads an interval-censored row (one of positive weight, not dropped by a phase design), logLik holds PROC HAZARD's objective, not a log-likelihood, and aic is computed from it. A step can change which rows are read (see n_rows), so on a step that starts or stops reading such rows delta_logLik is the difference of two different quantities. The trace's final line labels the final model's value, and criterion = "aic" warns once (class hzr_stepwise_sas_objective) as soon as a step's model reads such a row.

n_rows

Number of rows in the fit's data after this step (rows given weight 0 are counted). A multiphase fit drops every row where a variable in the model is missing, so entering a variable with missing values shrinks the sample every later step is tested on, and dropping it grows the sample back. A step that changes the rows of positive weight raises a warning of class hzr_stepwise_sample_changed naming the step, the variable and the row counts before and after. Which variables enter is not changed by this: compare n_rows across steps to see how many rows each test rested on.

Value

An object of class c("hzr_stepwise", "hazard"), the final fit augmented with:

steps

Data frame with one row per accepted / frozen action; see Details.

scope

Record of the candidate scope, plus force_in, force_out, and the frozen set. Four fields answer four different questions about the pins, and none substitutes for another: force_in and force_out are the names as given, including any that pinned nothing; force_in_resolved and force_out_resolved are what those names resolved to, as a column or term label, so a bare "_X1" appears there as "`_X1`"; unresolved is which of them resolved to nothing; and frozen is what the loop held in that the caller never named. The resolved fields carry the resolved names only: the frozen set is not merged into them, because frozen already records it and merging would list variables nobody asked for (#451). Resolving is not applying: a name that resolves to a variable the model does not contain is recorded here and still pins nothing, because force_in only keeps a variable that is already in. The resolved fields are not aligned element for element with the as-given ones, which stay longer by every name that resolved to nothing. A variable in frozen is in the final model if it was in the model when it was frozen, and out of it otherwise; see the The frozen set section. unresolved is a list with elements force_in, force_out and scope, each the names that could not be used and were therefore ignored (character() when none were). That is usually a name matching neither a column of data nor a term label; for force_in and force_out it also covers a name that IS a column but which no model formula can hold, such as "." or "", since such a column can never become a model term and so can never be pinned. The trace, and so print() and summary(), carries a line for each non-empty one, and a screen whose character scope was emptied this way says so where it stops. candidates, force_in and force_out here are the arguments as given, so a character candidates still lists the names that were ignored, as force_in and force_out do; read unresolved for which those were (#451).

criteria

Named list of the threshold / direction settings actually applied, plus n_uncomputable_scores (how many candidate scores were NA, counted once per step: an entry the score test could not score under criterion = "score"; an entry under criterion = "wald", or a removal under any criterion, whose Wald statistic could not be computed for want of a variance, wald_no_variance; or an entry under criterion = "aic" whose fit had no finite objective, nonfinite, or was fitted on different rows from the current model because the candidate is missing on some, rows_differ; or an entry under criterion = "aic" or "wald", or rescued by the Wald fallback under "score", whose refit ended below the current model's log-likelihood, which it contains, so that the refit cannot have converged, loglik_below_base), uncomputable_reasons (a named integer vector of why), wald_untested_removals and wald_untested_entries (the "var" / "var@phase" tokens of variables kept in, or left out, on a step whose Wald test for them could not be computed; a variable tested, or refused as loglik_below_base, at a later step is not listed. Entries are listed under criterion = "wald" only: under "score" an entry no test could reach is reported by its reason, such as fallback_no_variance) and stopped_uncomputable (TRUE when the last iteration had candidates for entry or for removal and could test none of them). Read uncomputable_reasons before treating an unscored candidate as a bad one: information_indefinite marks candidates whose effect is too large for the score test's approximation at zero, which are typically the strongest variables on offer rather than degenerate ones. Candidates with that cause, or with coefficient_diverging, information_nonpositive or nuisance_singular, are refit and tested by Wald automatically, counted in n_wald_fallbacks. nuisance_singular is a fault of the score test at the current model, whose information matrix could not be formed or inverted, so it applies to every candidate at that step and each of them is refitted. A candidate still reaches uncomputable_reasons when that refit fails, or when its cause is any other, for which no refit is attempted. Read uncomputable_reasons for which one it was in any given run. For every criterion it also carries refit_failures (the "var" / "var@phase" tokens of candidate moves whose refit errored, failed to converge, or was refused because the move would not change the model – a drop that removes no design column, #320), refit_failure_reasons (why each one failed or was refused: the refit's error message, that it did not converge, or that the move changes nothing; named by the same tokens), n_refit_failures, and stopped_refit_failed (TRUE when the run ended on an iteration in which a refit failed or a move was refused. A refit failure is a screen that could not test its candidates, rather than one that tested them and liked none; a refusal is determinate – the move was tested and would have left the model as it was. Read refit_failure_reasons for which it was). Check it before reading a zero-row steps as an honest null result.

trace_msg

Character vector of the trace lines, captured regardless of the trace flag.

elapsed

difftime from start to finish.

final_call

The call that produced this result.

print.hzr_stepwise returns x invisibly.

summary.hzr_stepwise returns a summary.hzr_stepwise object (extends summary.hazard) with ⁠$stepwise_steps⁠ and ⁠$stepwise_trace⁠ appended.

print.summary.hzr_stepwise returns x invisibly.

as.data.frame.hzr_stepwise returns the ⁠$steps⁠ data frame.

Selection direction and criterion

Two arguments shape the search. direction decides which moves are allowed at each step; criterion decides how a candidate move is scored and whether it is accepted.

direction = "forward"

Start from the base model and only add variables; the best eligible candidate enters each step until none clears the entry rule. Variables never leave once in.

direction = "backward"

Start from the base model, which must already hold every candidate, and only drop variables; the weakest term leaves each step until all survivors clear the retention rule. scope is not read, so a non-empty one is an error: protect terms with force_in.

direction = "both" (default)

Two-way stepwise: on every iteration, whether or not a variable entered, every term in the model, the base model's included, is re-tested and may be dropped unless it is in force_in or was frozen by max_move at an earlier iteration (see the The frozen set section). scope limits what may enter, not what may leave. This is the SAS SELECTION = STEPWISE strategy. max_move caps how often a single variable may oscillate before it is frozen.

criterion = "score" (default)

Accept moves on SAS-style significance thresholds, using the score (Q) statistic of the candidate coefficient; this reproduces C/SAS HAZARD's SELECTION statistic. Q is evaluated at the current model's MLE with the candidate's coefficient pinned at zero, so no candidate refit is needed: the reduced-model information is inverted once per step and reused across every candidate. Only the winner is refit. A candidate enters if its p-value is below slentry.

Score is an entry criterion; the drop path never refit per candidate in the first place, so removals are tested on the current model's Wald p-value against slstay, as SAS does.

For single-distribution fits, the score criterion computes the observed information numerically via the suggested numDeriv package and errors with a clear message if it is not installed; a multiphase fit uses the analytic Hessian instead and does not need it.

Following SAS, the variance used during selection is approximate: shaping-parameter covariances are ignored. This affects selection only; final-model standard errors are unchanged and still come from the full Hessian. Candidates must be single-column numeric main-effect terms; a factor is rejected with an error rather than skipped.

criterion = "wald"

Accept moves on SAS-style significance thresholds, using the Wald \chi^2 of the affected coefficient(s): a candidate enters if its p-value is below slentry, and a term is dropped if its p-value rises above slstay. Entry candidates are scored from a refit that adds the candidate (so its new coefficient can be tested); drop candidates are scored from the current model's Wald p-values without a per-candidate refit, and a single refit is run only after a drop is chosen. This was the default before version 1.2.0. It differs algorithmically from C/SAS HAZARD, so the two criteria can take different step paths (and select different variable sets) even when they converge to a similar final model.

criterion = "aic"

Accept any move with \Delta\mathrm{AIC} < 0 (a strictly better penalised fit), ignoring slentry / slstay. Entry candidates use the actual \Delta\mathrm{AIC} from the candidate refit; drop candidates use a Wald-to-likelihood-ratio approximation, \Delta\mathrm{AIC} \approx W - 2\,\mathrm{df}, computed from the current model without a per-candidate refit (the chosen drop is refit afterwards). Use this for a non-significance-based, information-criterion search.

The frozen set

A variable that has moved (entered or left) more than max_move times is frozen: it is held in the state it is in, in the model or out of it, for the rest of the run, and listed in ⁠$scope$frozen⁠.

The freeze takes effect at the end of the iteration in which the variable reached the cap. A two-way iteration (direction = "both") makes a forward step and then a backward step, so a variable that reaches the cap by entering and is dropped by the same iteration's backward step is frozen out. Its "frozen" row in ⁠$steps⁠ follows the "drop". A variable is therefore in the final model exactly when it was in the model when it was frozen, and ⁠$scope$frozen⁠ agrees with the final model.

Earlier versions froze such a variable as it entered and then dropped it anyway, so ⁠$scope$frozen⁠ could name a variable the final model did not contain (#580).

This is not PROC HAZARD's MOVE rule. PROC HAZARD makes one move per step, a removal before any entry; it counts a variable's exits (and, under NOSTEPWISE, its entries too), separately for each phase; and a variable at its limit can neither leave nor enter. Outside NOSTEPWISE it therefore freezes a variable only as it leaves. max_move counts entries and exits alike, by name across phases, which is why the SAS translator records ⁠MOVE=⁠ rather than mapping it. Adopting PROC HAZARD's rule would change which variables are selected, and is deferred; see https://github.com/ehrlinger/TemporalHazard/issues/378 and https://github.com/ehrlinger/TemporalHazard/issues/379.

See Also

hazard() for the base model and hzr_phase() for multiphase scopes; stepwise_trace() to retrieve the captured selection log.

Examples

data(avc)
avc <- na.omit(avc)
base <- hazard(survival::Surv(int_dead, dead) ~ age,
               data = avc, dist = "weibull", fit = TRUE,
               theta = c(mu = 0.01, nu = 0.5, 0))

sw <- hzr_stepwise(base, scope = ~ age + mal,
                   data = avc, direction = "forward")
print(sw)

Parameter names for a phase specification, in theta order

Description

Returns the names of the theta vector a multiphase hazard() fit builds from phases, in the order theta requires. Use it to check a hand-written starting vector against the specification it is meant to go with, before any fit runs.

Usage

hzr_theta_names(phases, covariates = NULL)

Arguments

Details

theta is positional and its entries are not on a common scale. For an early + late pair the layout is

early.log_mu, early.log_t_half, early.nu, early.m,
late.log_mu,  late.log_tau,     late.gamma, late.alpha, late.eta

so the late phase logs mu and tau while carrying gamma, alpha and eta on the natural scale. Wrapping the wrong element in log() produces a fit, not an error, which is why a comment describing the order is not enough and this returns the real thing. The order is a property of the specification, so it changes the moment a phase is added, removed or retyped.

Phase names come from names(phases); unnamed phases are labelled phase_1, phase_2, ... by the same validation hazard() applies, so the labels here are the labels a fit will use.

Phase names must be unique, and "total" and "time" are reserved; this applies the same validation hazard() does, so both reject them identically.

Value

A character vector whose order is the required theta order. Its length is the number of parameters the specification implies, so length(hzr_theta_names(phases)) is the length theta must have.

See Also

hzr_phase() for building a specification, hazard() for the fit whose theta this names.

Examples

phases <- list(early = hzr_phase("cdf"), late = hzr_phase("g3"))
hzr_theta_names(phases)

# Check a hand-written starting vector before fitting.
theta0 <- c(log(0.05), log(0.2), 0, -0.4, log(0.03), log(1), 1, 1, 1)
stopifnot(length(theta0) == length(hzr_theta_names(phases)))
setNames(theta0, hzr_theta_names(phases))

# With covariates on one phase.
hzr_theta_names(phases, covariates = list(early = c("age", "sex")))

Translate a SAS HAZARD job into a Quarto document

Description

Reads a SAS program containing ⁠PROC HAZARD⁠ and/or ⁠PROC HAZPRED⁠ blocks and emits a Quarto document of the equivalent hazard() and predict.hazard() calls.

Usage

hzr_translate_sas(path, out_dir = NULL, librefs = NULL)

Arguments

Details

Experimental: a job that translates does render (the emitted hazard() chunk binds its fit and asks for an actual fit), but this is a translation aid, not a turnkey reproduction, and some SAS constructs are refused rather than translated. See the Experimental section below.

Constructs the translator does not cover are recorded on the returned object and rendered as visible callouts, never dropped. A ⁠PROC HAZPRED⁠ job whose ⁠INHAZ=⁠ fitted model cannot be located emits a stop(), so the document fails to render rather than reporting over a model it did not load.

A phase-statement variable that ⁠PROC HAZARD⁠ cannot read as a name, such as AGE*SEX, LOG(AGE) or ⁠B SEX⁠, is left out of the model. (LOG() is read as LOG, as ⁠PROC HAZARD⁠ reads it: its lexer treats ⁠)⁠ as whitespace and ( as a change of state, hazard_l.l:32 and ⁠:56⁠.) A phase variable must be a NAME, ⁠[_A-Z][_A-Z0-9]*⁠ (hazard_l.l:39; phasevar : NAME, hazard_y.y:213), so ⁠PROC HAZARD⁠ rejects such a job at parse. The document raises a warning() above the fit, and the operand is recorded in ⁠$untranslated⁠.

Two more syntax errors take the same route. A value on a ⁠PROC HAZARD⁠ option that takes none, such as NOCOV=1 (hazard_y.y:65-76), is ignored in the fit. A TIME, EVENT, RCENSOR, LCENSOR or WEIGHT statement with more than one operand, which takes exactly one name (hazard_y.y:106-127), is fitted on the first. Both are warned about above the fit and recorded in ⁠$untranslated⁠. TIME or EVENT with no operand, and no other statement supplying the variable, leaves nothing to fit, so that job's fit chunk is a stop().

A job may contain more than one ⁠PROC HAZARD⁠ and/or ⁠PROC HAZPRED⁠ block. Every block is preserved: the first of a kind keeps the bare chunk name (fit, pred, pred_haz), later ones get fit_2, fit_3, pred_2, and so on, so the emitted document contains every model, not just the last one seen. job$outhaz is therefore a character vector, one element per ⁠PROC HAZARD⁠ block that set ⁠OUTHAZ=⁠ (in the order the blocks appear). Each ⁠PROC HAZPRED⁠ block's ⁠INHAZ=⁠ is resolved independently against that vector, matching the most recently written ⁠OUTHAZ=⁠ at that point in the file, mirroring SAS itself, where a later ⁠OUTHAZ=⁠ write overwrites the dataset an earlier one wrote under the same name. If a job's own ⁠OUTHAZ=⁠ values don't cover it, librefs is tried next; distinct external ⁠INHAZ=⁠ values each get their own loaded-fit chunk. When neither resolves and the job holds more than one local fit, which fit a predict() call belongs to is genuinely unknown; the emitted call falls back to referencing fit and the ambiguity is recorded in untranslated, never guessed at silently.

Value

An hzr_sas_job object, invisibly. ⁠$calls⁠ holds the emitted calls keyed by chunk label, ⁠$grid⁠ the last prediction grid seen and ⁠$inhaz⁠ the first unresolved ⁠INHAZ=⁠ (not all of each, when a job has several), ⁠$untranslated⁠ the recorded gaps, ⁠$coverage⁠ the token counts, ⁠$notes⁠ any callout attached to a chunk by label, emitted immediately above that chunk in the rendered document, and ⁠$not_run⁠ the chunks SAS never ran, named by label, each holding the ⁠%MACRO⁠ the job defines but never calls around it. Those chunks are rendered with eval: false; a caller evaluating ⁠$calls⁠ itself should skip them.

What the emitted document needs to run

Base R at the version this package declares in DESCRIPTION (Depends), plus TemporalHazard itself. The emitted chunks never reach into this package's internals, so the document renders in your session, and in a colleague's, without anything further installed, with one exception. A translated SELECTION screen on a multiphase job with interval-censored rows (an ICENSOR job; LCENSOR is left truncation and does not need it) needs the suggested package numDeriv. Without it a screen that tests an entry stops and says so. A screen that completes but cannot test a removal (a BACKWARD screen, whose base has no usable variance) is caught by the emitted check, which warns and names the variables.

Experimental

The emitted document renders: the hazard() chunk binds its fit to a name and passes fit = TRUE, so the predict() chunks have something to predict from. A SELECTION statement is translated into an hzr_stepwise() call carrying the job's own candidates, per-variable flags and thresholds (#160); the screen is real, and may select a different model than ⁠PROC HAZARD⁠ did, which the emitted document says in a callout above the chunk. SELECTION options with no faithful translation are refused outright rather than mistranslated, emitting a stop() in place of the fit: FAST, MAXVARS, RESTRICT, a per-variable ⁠MOVE=⁠ or ⁠ORDER=⁠, and a variable held by ⁠/I⁠ in one phase but movable in another. So is LCENSOR combined with ICENSOR, which one time_lower argument cannot express (#155). A job with ICENSOR is fitted with objective = "sas", the interval term ⁠PROC HAZARD⁠ accumulates, so it reproduces ⁠PROC HAZARD⁠'s estimates; the value it reports is then ⁠PROC HAZARD⁠'s objective, not a log-likelihood, and a callout above the chunk says so (#543). Degenerate intervals are resolved as ⁠PROC HAZARD⁠ does before its fit: a lower bound equal to the time makes the row an exact event, and one that is missing, negative or after the time drops the row from the fit (not from your data frame); the status chunk warns with both counts. Prediction grids the parser cannot resolve are refused whole, and the predict() chunks that would have read such a grid become a stop() naming it, rather than a predict(newdata = ) over a name no chunk builds. An unresolved ⁠INHAZ=⁠ stops the render on purpose.

On a fit loaded from an external ⁠INHAZ=⁠ dataset, point predictions work but se.fit = TRUE is refused when ⁠PROC HAZARD⁠ estimated a late shape parameter on a composite scale (the generic unconstrained three-phase case, not an exotic one). A translated ⁠PROC HAZPRED⁠ block asks for confidence limits unless the SAS job says NOCL, so such a job stops at its predict() chunks.

Repeated events

A job that builds its fit input with the SAS macro ⁠%repeat⁠ has that call translated to hzr_repeated_events(). The emitted chunk renames the function's output columns to the names the job gives them, upper-cased like every name the translator emits, so the fit reads them. An input column named like one of those outputs is dropped first, with a warning, because the macro overwrites it. The macro's input is built by the job's own DATA steps, which are not translated, so the document stops until that input is assigned. Any step between the macro and the fit that names the macro's output, or uses a macro variable that might, is not translated either; a plain ⁠PROC SORT⁠ that only reorders rows is the one step let through. Its chunk stops and quotes the step, for the reader to replace with R code or to delete if the step leaves the output unchanged. The same holds between two ⁠%repeat⁠ calls when the second reads the first's output. A step that changes the output without naming it, such as a macro that writes it internally, is not detected.

Comparing a translated fit against a SAS listing

The emitted hzr_phase("g3", ...) call carries the TAU, GAMMA, ALPHA and ETA that PARMS specified. ⁠PROC HAZARD⁠ does not always report the values it was given: when ALPHA is fixed at 1 it rewrites the late phase before fitting, pinning TAU at 1 and folding the two shape parameters into one exponent, because at ⁠alpha = 1, tau = 1⁠ the G3 form collapses to t^{\gamma\eta} and only that product is identified. With ETA fixed it reports GAMMA as \gamma\eta and ETA as 1; otherwise it reports ETA as \gamma\eta and GAMMA as 1.

So a listing can print a GAMMA/ETA pair that differs from the job's own PARMS text, and from the emitted call, while describing the same fitted model. Compare the product \gamma\eta, not the two parameters separately, and expect TAU to read 1 and fixed. The log-likelihood, MUL and every prediction agree regardless. Most jobs are unaffected in practice: ⁠GAMMA = 1 FIXGAMMA⁠ is the usual companion to ⁠ALPHA = 1 FIXALPHA⁠, and it makes the rewrite an identity.

The one case that is a genuine model difference (ALPHA fixed at 1 with GAMMA and ETA both estimated, where ⁠PROC HAZARD⁠ fixes ETA and fits one parameter fewer than hazard() would) is recorded in ⁠$untranslated⁠ rather than left to be discovered in the comparison.

⁠$coverage⁠ counts tokens the parser recognised; it is not evidence that the emitted calls execute, and a job can report full coverage with an empty ⁠$untranslated⁠ while its document still errors on render. See the 1.2.2 NEWS.md entry. The function's API, the hzr_sas_job field layout and the emitted document format are all expected to change.

Examples

job <- hzr_translate_sas(
  system.file("extdata", "hz-example.sas", package = "TemporalHazard")
)

Test if an object is an hzr_phase

Description

Test if an object is an hzr_phase

Usage

is_hzr_phase(x)

Arguments

Value

Logical scalar.

Examples

is_hzr_phase(hzr_phase("cdf"))
is_hzr_phase("not a phase")

OMC: Open Mitral Commissurotomy

Description

Data for 339 patients who underwent open mitral commissurotomy at the University of Alabama Birmingham. Records repeated thromboembolic events (up to 3 per patient) as indicators; the event dates are not included.

Usage

omc

Format

A data frame with 339 rows and 7 variables:

study

Patient identifier

te1

Indicator for first thromboembolic event

te2

Indicator for second thromboembolic event

te3

Indicator for third thromboembolic event

int_dead

Follow-up interval to death or last contact (months)

dead

Death indicator (1 = dead, 0 = censored)

opdjul

Operation date (Julian)

Source

University of Alabama Birmingham cardiac surgery registry.

See Also

Other datasets: avc, cabgkul, tga, uslife2023, valves

Predict from a hazard model object

Description

Produces prediction outputs from a hazard object. Supports multiple prediction types including linear predictor, hazard, survival probability, and cumulative hazard.

Usage

## S3 method for class 'hazard'
predict(
  object,
  newdata = NULL,
  type = c("hazard", "linear_predictor", "survival", "cumulative_hazard"),
  decompose = FALSE,
  se.fit = FALSE,
  level = 0.95,
  conf.type = c("log-log", "logit"),
  ...
)

Arguments

predict(fit, newdata, type = "survival", se.fit = TRUE,
        level = 2 * stats::pnorm(1) - 1, conf.type = "logit")

Details

For Weibull models with survival or cumulative_hazard predictions:

• Cumulative hazard: H(t|x) = (mu*t)^nu * exp(eta)

• Survival: S(t|x) = exp(-H(t|x))

Time values must be positive and finite. If newdata contains a time column, it will be used; otherwise, the time vector from the fitted object is used. For models fit with time_windows, predictions for type = "linear_predictor" or "hazard" also require time values (via newdata$time or fitted-time fallback) so window-specific coefficients can be selected.

See the section "How newdata is evaluated" for what is recomputed from newdata and when predict() warns.

Value

When se.fit = FALSE (default), a numeric vector of predictions. When se.fit = TRUE, a data frame with columns fit, se.fit, lower, upper (delta-method point estimate, standard error, and confidence limits at level). se.fit is always the standard error of fit on its own scale: for type = "survival" that is S * se(H), the delta method applied to S = exp(-H). The survival limits are built from se(H) on the conf.type scale, so they are not ⁠fit +/- z * se.fit⁠. ⁠PROC HAZPRED⁠ prints no standard error, only the limits. For multiphase type = "cumulative_hazard" with decompose = TRUE, a long data frame (time, component, fit, se.fit, lower, upper); with decompose = TRUE and se.fit = FALSE, a wide data frame of per-phase contributions.

How newdata is evaluated

predict() evaluates the model's formulas on newdata as given, as stats::predict.lm() does. It uses the fit's factor levels and contrasts, and the centering, basis and knots that a top-level scale(), poly(), ns() or bs() term recorded. Everything else is recomputed from newdata, so a prediction can differ from the fit without any error. predict() warns, naming the cause, in three such cases. The predicted values are the same with or without the warning.

• A term that computes a statistic over the rows. In I(age - mean(age)), I(scale(age)^2) or I(as.integer(factor(grp))), the mean, the scaling or the factor coding comes from newdata's rows, so a row's prediction depends on which other rows are given. Compute such a variable in the data before fitting, and supply it in newdata. The check knows a fixed list of R's functions, among them mean(), median(), min(), max(), quantile(), sd(), IQR(), ave(), rank(), length(), scale(), factor() and cut() with a count of breaks. A function not on it, including one you write, is recomputed from newdata's rows just the same, with no warning: the list is a floor, not a boundary.

• A column of another type than the fit saw. A numeric column given as character compares as text ("154.6" > 50 is FALSE), and a difftime in other units is used in those units. The check compares against the fitting data the fit kept, which fits saved before TemporalHazard 1.1.0 do not have, so those fits are not checked. A factor given as its level labels, or an integer for a double, is not a mismatch.

• A design rebuilt under this session's contrasts. A formula fit saved by version 1.2.10 or earlier kept no record of its contrasts, and its design is rebuilt under options(contrasts =). predict() warns when that option names a function other than contr.treatment or contr.poly and the rebuilt design is used. A redefined contr.treatment is not detected.

Two cases are not detected:

• A constant the formula reads from its environment, such as cutoff in I(age > cutoff), is read when you predict, so a value changed since the fit is used.

• A comparison of strings, such as I(grp > "b"), follows the session's collation (LC_COLLATE), which can order strings differently from the session that fitted the model.

See Also

hazard() for model fitting, summary.hazard() for model summaries, hzr_phase() for multiphase temporal shapes.

vignette("prediction-visualization") for detailed prediction workflows including decomposed hazard plots and patient-specific curves.

Examples

# -- Basic predictions ------------------------------------------------
set.seed(1)
fit <- hazard(time = rexp(50, 0.3), status = rep(1L, 50),
              theta = c(0.3, 1.0), dist = "weibull", fit = TRUE)
predict(fit, type = "survival")
predict(fit, newdata = data.frame(time = c(1, 2, 5)),
        type = "cumulative_hazard")

# -- Patient-specific survival curves ---------------------------------
set.seed(1001)
n   <- 180
dat <- data.frame(
  time   = rexp(n, rate = 0.35) + 0.05,
  status = rbinom(n, size = 1, prob = 0.6),
  age    = rnorm(n, mean = 62, sd = 11),
  nyha   = sample(1:4, n, replace = TRUE),
  shock  = rbinom(n, size = 1, prob = 0.18)
)
fit2 <- hazard(
  survival::Surv(time, status) ~ age + nyha + shock,
  data  = dat,
  theta = c(mu = 0.25, nu = 1.10, beta1 = 0, beta2 = 0, beta3 = 0),
  dist  = "weibull", fit = TRUE
)

new_patients <- data.frame(
  time = c(0.5, 1.5, 3.0),
  age  = c(50, 65, 75),
  nyha = c(1, 3, 4),
  shock = c(0, 0, 1)
)
# Compute predictions from the clean covariate frame before adding columns
surv   <- predict(fit2, newdata = new_patients, type = "survival")
cumhaz <- predict(fit2, newdata = new_patients, type = "cumulative_hazard")
new_patients$survival          <- surv
new_patients$cumulative_hazard <- cumhaz
new_patients


# -- Grouped survival curves ---------------------------------------
if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(ggplot2)

  t_grid <- seq(0.05, max(dat$time), length.out = 80)
  profiles <- data.frame(
    label = c("Low risk (age 50, NYHA I)",
              "High risk (age 75, NYHA IV)"),
    age   = c(50, 75),
    nyha  = c(1, 4),
    shock = c(0, 1)
  )

  curve_list <- lapply(seq_len(nrow(profiles)), function(i) {
    nd <- data.frame(
      time  = t_grid,
      age   = profiles$age[i],
      nyha  = profiles$nyha[i],
      shock = profiles$shock[i]
    )
    nd$survival <- predict(fit2, newdata = nd, type = "survival") * 100
    nd$profile  <- profiles$label[i]
    nd
  })
  curve_df <- do.call(rbind, curve_list)

  ggplot(curve_df, aes(time, survival, colour = profile)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 100)) +
    labs(x = "Months after surgery",
         y = "Freedom from death (%)",
         title = "Predicted survival by risk profile",
         colour = NULL) +
    theme_minimal()
}



# -- Multiphase predictions with decomposition --------------------
set.seed(42)
n   <- 200
dat <- data.frame(
  time   = rexp(n, rate = 0.25) + 0.01,
  status = rbinom(n, size = 1, prob = 0.65)
)
fit_mp <- hazard(
  survival::Surv(time, status) ~ 1,
  data   = dat,
  dist   = "multiphase",
  phases = list(
    early = hzr_phase("cdf", t_half = 0.5, nu = 2, m = 0,
                       fixed = "shapes"),
    late  = hzr_phase("cdf", t_half = 5,   nu = 1, m = 0,
                       fixed = "shapes")
  ),
  fit     = TRUE,
  control = list(n_starts = 5, maxit = 1000)
)

t_grid <- seq(0.01, max(dat$time) * 0.9, length.out = 100)
nd     <- data.frame(time = t_grid)

# Overall survival
predict(fit_mp, newdata = nd, type = "survival")

# Per-phase decomposed cumulative hazard
decomp <- predict(fit_mp, newdata = nd,
                  type = "cumulative_hazard", decompose = TRUE)
head(decomp)

Predictions from a fit loaded out of a SAS ⁠OUTHAZ=⁠ dataset

Description

Rebuilds the multiphase model the ⁠OUTHAZ=⁠ dataset describes (which phases are in it, their shapes, and the fitted parameter vector) and then predicts exactly as predict.hazard() does.

Usage

## S3 method for class 'hzr_outhaz'
predict(
  object,
  newdata,
  type = c("hazard", "survival", "cumulative_hazard"),
  decompose = FALSE,
  se.fit = FALSE,
  level = 0.95,
  conf.type = c("log-log", "logit"),
  ...
)

Arguments

predict(fit, newdata, type = "survival", se.fit = TRUE,
        level = 2 * stats::pnorm(1) - 1, conf.type = "logit")

Details

newdata is required. An ⁠OUTHAZ=⁠ dataset holds a converged model and no data, so there is no fitted time vector to fall back on; without newdata the prediction would be over nothing.

Value

What predict.hazard() returns for a fit with no covariates: a numeric vector of predictions when se.fit = FALSE, and when se.fit = TRUE a data frame with columns fit, se.fit, lower and upper (the confidence limits at level). The decomposed long-format return of predict.hazard() is not reachable here, because decompose = TRUE is refused.

What is reconstructed, and what is refused

A phase is in the model when its intercept row (E0, C0, L0) carries ⁠_STATUS_ = 1⁠, which is the test ⁠PROC HAZPRED⁠ itself applies. The early phase becomes hzr_phase("cdf") and the late phase hzr_phase("g3"), with the shape estimates read off the DELTA/THALF/NU/M and TAU/GAMMA/ALPHA/ETA rows.

These cases error rather than return a number that looks like a prediction: a dataset carrying covariates (their coefficients cannot be matched to newdata columns from the file alone); a non-zero DELTA (the early-phase time transformation is not implemented); and a G1FLAG that disagrees with the signs of the M and NU estimates.

With se.fit = TRUE there are three more, because the stored covariance is on SAS's estimation scale and has to be mapped onto this package's:

• a fit constrained by FIXMNU1, which ties M to 1/NU;

• a fit estimating GAMMA, ALPHA or ETA on one of PROC HAZARD's composite late-phase scales, log(GAMMA*ETA - 2) or log(GAMMA*ETA/ALPHA - 2) rather than log() of the parameter. This is the ordinary unconstrained late phase, not an exotic case: with G3FLAG = 1 and no FIXGE2/FIXGAE2, ETA and ALPHA are always on a composite scale and GAMMA is whenever ETA is fixed. So an early + constant + late fit with any free late shape gets point predictions but no standard errors;

• a fit under FIXGE2 or FIXGAE2 where a late parameter is derived from an estimated one (ETA = 2/GAMMA, GAMMA = 2/ETA, ALPHA = GAMMA*ETA/2). The derived parameter has no covariance row, so its contribution to the variance would be dropped without trace.

See Also

hzr_read_outhaz(), predict.hazard().

Examples

f <- system.file("extdata", "outhaz-fixture.rds", package = "TemporalHazard")
if (nzchar(f)) {
  fit <- hzr_read_outhaz(f)
  predict(fit, newdata = data.frame(time = c(1, 6, 12)), type = "survival")
}

Print method for fitted hazard models

Description

Compact one-block summary of a fitted hazard object: sample size, number of predictors, distribution, theta vector, and log-likelihood, followed by the "Not done in this run" block described in hazard(). S3 dispatch only: users call print(fit) rather than invoking this directly.

Usage

## S3 method for class 'hazard'
print(x, ...)

Arguments

Value

x, invisibly.

Print method for hzr_calibrate

Description

Print method for hzr_calibrate

Usage

## S3 method for class 'hzr_calibrate'
print(x, digits = 3, ...)

Arguments

Value

The object x of class c("hzr_calibrate", "data.frame"), invisibly. The data frame has one row per quantile group and columns: group (group index), n (group size), events (event count), mean, min, max (covariate summary within group), prob (observed event probability), link_value (transformed probability on the link scale). When stratified via the by argument, a by column is also present. Attributes: "link" (the transform applied, e.g. "logit") and "groups" (number of quantile bins).

Print method for hzr_deciles

Description

Print method for hzr_deciles

Usage

## S3 method for class 'hzr_deciles'
print(x, digits = 3, ...)

Arguments

Value

The object x of class c("hzr_deciles", "data.frame"), invisibly. The data frame has one row per risk group and columns: group (integer group index, 1 = lowest risk), n (group size), events (observed event count; weighted for a weighted fit), expected (expected event count from model predictions, net of any entry time, and weighted as events is), observed_rate, expected_rate (events / n; per unit of case weight for a weighted fit, and NA for a group whose weights sum to 0), chi_sq (per-group (O-E)^2/V contribution, where V is the Poisson variance of the weighted count, equal to E for an unweighted fit), p_value (1-df chi-square upper-tail p), mean_survival, mean_cumhaz (mean predicted values in group). An "overall" attribute contains the omnibus chi-square test (fields: chi_sq, df, p_value, time, groups, total_events, total_expected, n_included, n_excluded).

Print method for hzr_gof

Description

Print method for hzr_gof

Usage

## S3 method for class 'hzr_gof'
print(x, digits = 3, ...)

Arguments

Value

The object x of class c("hzr_gof", "data.frame"), invisibly. The data frame has one row per time point and columns: time, n_risk, n_event, n_censor, km_surv (Kaplan-Meier survival), km_cumhaz, par_surv (parametric survival), par_cumhaz, cum_observed (cumulative observed events), cum_expected (cumulative expected events from model), residual (cum_expected - cum_observed). Multiphase models additionally include par_cumhaz_<phase> columns for per-phase cumulative hazard contributions. A "summary" attribute contains scalar diagnostics: total_observed, total_expected, final_residual, dist, n, n_tallied. When n_tallied is below n, the printed totals cover only those subjects, and a note says so.

Print method for hzr_kaplan

Description

Print method for hzr_kaplan

Usage

## S3 method for class 'hzr_kaplan'
print(x, digits = 4, n = 20, ...)

Arguments

Value

The object x of class c("hzr_kaplan", "data.frame"), invisibly. The data frame has one row per event time (or all times when event_only = FALSE) and columns: time (follow-up time), n_risk (number at risk), n_event (events in interval), n_censor (censored observations in interval), survival (Kaplan-Meier survival estimate), std_err (Greenwood standard error on log-hazard scale), cl_lower, cl_upper (logit-transformed confidence limits on the survival scale), cumhaz (Nelson-Aalen cumulative hazard), hazard (interval hazard estimate), density (estimated event density), life (life-table life expectancy contribution).

Print method for hzr_sas_job

Description

Shows what hzr_translate_sas() made of one SAS job: the chunks it emitted, how many of the job's tokens it mapped, and one line per construct it could not translate. Read the untranslated lines before trusting the document – each one names something the emitted R does differently from PROC HAZARD, or not at all.

Usage

## S3 method for class 'hzr_sas_job'
print(x, ...)

Arguments

Value

The object x, invisibly.

See Also

hzr_translate_sas(), which builds the object.

Print method for hazard summary objects

Description

Formatted console display of summary.hazard() output: distribution, phase list (for multiphase), coefficient table with standard errors, and log-likelihood. When the post-fit Hessian is ill-conditioned or not positive-definite, a note warns that the standard errors may be unreliable, and a further note names the parameters spanning a weakly identified direction when one was found. A "Not done in this run" block is always printed: it lists each step this fit did not perform, with the reason, and reads "none" when nothing was lost. S3 dispatch only: users call print(summary(fit)) rather than invoking this directly.

Usage

## S3 method for class 'summary.hazard'
print(x, ...)

Arguments

Value

x, invisibly.

Extract the captured console trace from an hzr_stepwise fit

Description

Every run of hzr_stepwise() records the header, per-step lines, and final summary regardless of the trace flag. This accessor returns the full character vector for display or logging.

Usage

stepwise_trace(fit)

Arguments

Value

Character vector, one element per console line.

See Also

hzr_stepwise(), which produces the object this accessor reads.

Examples

data(avc)
avc <- na.omit(avc)
base <- hazard(survival::Surv(int_dead, dead) ~ age,
               data = avc, dist = "weibull", fit = TRUE,
               theta = c(mu = 0.01, nu = 0.5, 0))

sw <- hzr_stepwise(base, scope = ~ age + mal,
                   data = avc, direction = "forward")
cat(stepwise_trace(sw), sep = "\n")

Summarize a hazard model

Description

Returns a compact summary of a hazard object, including model metadata, fit diagnostics, and coefficient-level statistics when available.

Usage

## S3 method for class 'hazard'
summary(object, ...)

Arguments

Value

An object of class summary.hazard. Its log_lik is the log-likelihood at the estimates, and NA for a fit with objective = "sas" that read an interval-censored row (one of positive weight, not dropped by a phase design): there the fitted objective is PROC HAZARD's interval-mean-hazard objective, not a log-likelihood. That value is always in objective_value, and objective says which of the two it is ("likelihood" or "sas"). For a Weibull fit whose scale mu, or its variance, cannot be represented, the coefficient table has a log(mu) row under mu, with the standard error of log(mu); it is not tested against 0.

See Also

hazard() for model fitting, predict.hazard() for predictions.

vignette("fitting-hazard-models") for fitting workflows, vignette("inference-diagnostics") for bootstrap CIs and diagnostics.

Examples

# -- Single-phase Weibull summary ------------------------------------
set.seed(1)
fit <- hazard(time = rexp(30, 0.5), status = rep(1L, 30),
              theta = c(0.3, 1.0), dist = "weibull", fit = TRUE)
summary(fit)


# -- Multiphase model summary ----------------------------------------
set.seed(42)
n   <- 200
dat <- data.frame(
  time   = rexp(n, rate = 0.25) + 0.01,
  status = rbinom(n, size = 1, prob = 0.65)
)
fit_mp <- hazard(
  survival::Surv(time, status) ~ 1,
  data   = dat,
  dist   = "multiphase",
  phases = list(
    early = hzr_phase("cdf", t_half = 0.5, nu = 2, m = 0,
                       fixed = "shapes"),
    late  = hzr_phase("cdf", t_half = 5,   nu = 1, m = 0,
                       fixed = "shapes")
  ),
  fit     = TRUE,
  control = list(n_starts = 5, maxit = 1000)
)
summary(fit_mp)

TGA: Transposition of the Great Arteries

Description

Survival data for 470 patients who underwent the arterial switch operation for transposition of the great arteries at Boston Children's Hospital and Children's Hospital of Philadelphia. Used for sensitivity analysis and internal validation examples.

Usage

tga

Format

A data frame with 470 rows and 14 variables:

study

Patient identifier

simple

Simple TGA indicator (0/1)

dextroin

D-looped transposition indicator (0/1)

ca_1rl2c

Coronary artery pattern indicator

hyaaproc

Hybrid approach procedure indicator (0/1)

no_tca

No total circulatory arrest indicator (0/1)

tca_time

Total circulatory arrest time (minutes)

age_days

Age at operation (days)

arciopyr

Aortic cross-clamp time per year

dead

Death indicator (1 = dead, 0 = censored)

int_dead

Follow-up interval to death or last contact (months)

source

Source institution (BCH or CHOP)

ca1_2_l

Coronary artery configuration (1/2/L)

opyear

Year of operation

Source

Boston Children's Hospital and Children's Hospital of Philadelphia.

See Also

Other datasets: avc, cabgkul, omc, uslife2023, valves

US Life Table 2023: All-Interval-Censored SAS Parity Anchor

Description

The published NCHS United States life table for 2023, expressed on a synthetic radix of 100,000, as a ⁠PROC HAZARD⁠ job consumes it: one interval-censored row per year of age, weighted by the number of deaths falling in that year. Aggregate published counts only: no patient-level data and no PHI.

Usage

uslife2023

Format

A data frame with 124 rows and 3 variables:

age_l

Lower bound of the age interval, in years (integer)

age_u

Upper bound of the age interval, in years (integer); age_u - age_l is exactly 1 on every row

d_all

Deaths in the interval on a 100,000 radix. This is the ICENSOR variable, and it is a weight, not an indicator (a fact the weighted life table forces and 335 occurrences of ⁠icensor icens_wt=il_dead⁠ across the SAS corpus corroborate). Sums to 100000.0125; ranges from 0.2352 to 3620.335

Details

This is the reference fixture for objective = "sas". It is the cleanest available anchor for that form for two reasons. Every row is interval-censored, so nothing else dilutes the signal; and every interval is exactly one year wide, so \log(u - l) = 0 and the width term switches off, isolating the S(u)\,\Delta\Lambda core. It also settles by itself the rival hypothesis that SAS bakes in a constant width divisor: that reading is off by 593,146 log-likelihood units here.

Two rows of the source table (ages 119–120 and 124–125) carry d_all == 0 and are dropped by the SAS job's own ⁠IF D_ALL=0 THEN DELETE⁠, leaving 124 of 126. The age grid is therefore not contiguous, which is why the fixture carries explicit age_l/age_u columns rather than an implied one-year step.

Source

National Center for Health Statistics, United States Life Tables, 2023. Built by data-raw/make_data.R from inst/extdata/uslife2023.csv, which needs no SAS license to reproduce.

References

The SAS reference fit is distributions/hz.icall.lst in that study, whose printed log-likelihood is -410414 with ⁠Number of events conserved⁠ = 100000.

See Also

hazard for the objective argument this fixture anchors.

Other datasets: avc, cabgkul, omc, tga, valves

Examples

data(uslife2023)

# Every interval is one year wide -- the property that makes this the anchor.
stopifnot(all(uslife2023$age_u - uslife2023$age_l == 1))

# Deaths by age, on the published 100,000 radix.
plot(uslife2023$age_l, uslife2023$d_all, type = "h",
     xlab = "Age (years)", ylab = "Deaths per 100,000",
     main = "US Life Table 2023: deaths by year of age")

Valves: Primary Heart Valve Replacement

Description

Data for 1,533 patients who underwent primary heart valve replacement. The largest multivariable example dataset with multiple endpoints including death, prosthetic valve endocarditis (PVE), bioprosthesis degeneration, and reoperation.

Usage

valves

Format

A data frame with 1533 rows and 19 variables:

age_cop

Age at operation (years)

nyha

NYHA functional class (1–4)

mitral

Mitral valve position indicator (0/1)

double_

Double valve replacement indicator (0/1)

ao_pinc

Aortic position, incompetence (0/1)

black

Black race indicator (0/1)

i_path

Ischemic pathology indicator (0/1)

nve

Native valve endocarditis indicator (0/1)

mechvalv

Mechanical valve indicator (0/1)

male

Male sex indicator (0/1)

int_dead

Follow-up interval to death or last contact (months)

dead

Death indicator (1 = dead, 0 = censored)

int_pve

Follow-up interval to PVE or last contact (months)

pve

PVE indicator (1 = PVE, 0 = censored)

bio

Bioprosthesis indicator (0/1)

int_rdg

Follow-up interval to degeneration or last contact (months)

reop_dg

Reoperation for degeneration indicator (0/1)

int_reop

Follow-up interval to reoperation or last contact (months)

reop

Reoperation indicator (0/1)

Source

Cleveland Clinic Foundation heart valve replacement registry.

See Also

vignette("fitting-hazard-models"), vignette("prediction-visualization")

Other datasets: avc, cabgkul, omc, tga, uslife2023

Examples

data(valves)
valves_cc <- na.omit(valves)

# Kaplan-Meier for two endpoints
km_death <- survival::survfit(
  survival::Surv(int_dead, dead) ~ 1, data = valves_cc)
km_pve <- survival::survfit(
  survival::Surv(int_pve, pve) ~ 1, data = valves_cc)

plot(km_death, xlab = "Months after valve replacement", ylab = "Survival",
     main = "Valves: Death and PVE endpoints")
lines(km_pve, col = "red")
legend("bottomleft", c("Death", "PVE"), col = c("black", "red"), lty = 1)

Extract variance-covariance matrix from hazard model

Description

Returns the estimated variance-covariance matrix of the fitted coefficients.

Usage

## S3 method for class 'hazard'
vcov(object, ...)

Arguments

Value

A numeric matrix containing the estimated variance-covariance matrix of the fitted coefficients, with rows and columns named by the coefficient labels (phase-prefixed for multiphase models, e.g. early.x). Rows and columns for parameters held fixed (e.g. fixed shape parameters) are NA because they carry no variance; the finite free-parameter block is still usable. For Conservation-of-Events fits the conserved phase log_mu normally carries a variance: it is removed from the optimizer search but the vcov is the full-information matrix at the optimum (the CoE solution is the unconstrained MLE). That recomputation needs a Hessian (analytic, or from numDeriv for left- or interval-censored rows) that can be inverted; if it fails the fit warns, records why in degraded_causes, and has no covariance matrix, so vcov() returns NA. For a Weibull fit whose scale mu, or its variance, cannot be represented, mu's row and column are NA, the reason is recorded in the fit's degraded_causes, and summary() shows the variance of log(mu) instead. Returns a scalar NA only when the model has not been fitted or no covariance matrix is available.

Examples

set.seed(1)
fit <- hazard(time = rexp(30, 0.5), status = rep(1L, 30),
              theta = c(0.3, 1.0), dist = "weibull", fit = TRUE)
vcov(fit)