Package {TemporalHazard}

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).

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
)

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.

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 central-difference derivatives (via .hzr_phase_derivatives()) for shape parameters (log_t_half, nu, m).

Usage

.hzr_gradient_multiphase(
  theta,
  time,
  status,
  time_lower = NULL,
  time_upper = NULL,
  x = NULL,
  weights = NULL,
  phases,
  covariate_counts,
  x_list,
  ...
)

Arguments

Value

Numeric vector of length length(theta) – the gradient. Returns a zero vector if any component is non-finite (guards optimizer).

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.

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,
  return_gradient = FALSE,
  return_hessian = FALSE,
  ...
)

Arguments

Value

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

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,
  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.

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 t_half, nu, and m using central finite differences on hzr_decompos(). The log_t_half derivative is obtained via the chain rule: d\Phi/d(\log t_{1/2}) = t_{1/2} \cdot d\Phi/dt_{1/2}.

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(...)).

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, 0 for "constant".

Usage

.hzr_phase_n_shape(phase)

Arguments

Value

Integer: 3 or 0.

Extract starting values from a phase specification

Description

Returns initial theta sub-vector on the estimation (internal) scale: log(mu), log(t_half), nu, m, followed by zeros for covariate coefficients.

Usage

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

Arguments

Value

Named 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 "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
)

Arguments

Value

Numeric n x p Jacobian of H(t|x).

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.

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
)

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)) 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.

Select the phase whose log_mu will be solved by conservation

Description

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

Usage

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

Arguments

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.

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)

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, 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, 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

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(),
  ...
)

Arguments

Details

Control parameters:

• maxit: Maximum iterations (default 1000)

• reltol: Relative parameter change tolerance (default 1e-5)

• abstol: Absolute gradient norm tolerance (default 1e-6)

• method: Optimization method: "bfgs" or "nm" (default "bfgs")

• condition: Condition number control (default 14)

• nocov, nocor: Suppress covariance/correlation output (legacy; no-op in M2)

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, converged, se, vcov, counts, message; all NULL when fit = FALSE), and engine (implementation tag, "native-r-m2").

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.

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
)

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

Arguments

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.

n_success

Number of successfully converged replicates.

n_failed

Number of replicates that failed to converge.

See Also

hazard() for model fitting, vcov.hazard() for Hessian-based standard errors.

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)

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

Clamp probabilities away from 0 and 1

Usage

hzr_clamp_prob(p, eps = 1e-12)

Arguments

Value

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

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.

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 implements the workflow of the SAS deciles.hazard.sas macro. Patients are ranked by predicted cumulative hazard at a specified time point, grouped into quantile bins, and each bin is tested with a chi-square goodness-of-fit statistic. Subjects censored before the requested horizon are excluded from the observed-vs-expected comparison.

Value

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

group

Integer group label (1 = lowest risk).

n

Number of observations in the group.

events

Observed event count by the requested horizon (event indicator = 1 and event time <= time).

expected

Expected event count by the requested horizon, computed as the sum of individual event probabilities (1 - S(time)).

observed_rate

Observed event rate (events / n).

expected_rate

Expected event rate (expected / n).

chi_sq

Chi-square contribution: (events - expected)^2 / expected.

p_value

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

mean_survival

Mean predicted survival probability in the group.

mean_cumhaz

Mean predicted cumulative hazard in the group.

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. This single function generates all temporal phase shapes used in multiphase hazard models.

Usage

hzr_decompos(time, t_half, nu, m)

Arguments

Value

A named list with three 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.

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. Both collapse onto the three parameters here. 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.

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)

Goodness-of-fit: observed vs. predicted events

Description

Compare a fitted hazard model against the nonparametric Kaplan-Meier estimate by computing observed and expected (parametric) event counts at each distinct event time. 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

At each observed event time the function computes:

• The Kaplan-Meier survival and cumulative hazard.

• The parametric survival and cumulative hazard from the fitted model (and per-phase components for multiphase models).

• Cumulative observed events vs. cumulative expected events (sum of individual cumulative hazards for those exiting the risk set at each time).

• The running residual (expected minus observed).

Perfect model fit implies the expected and observed event counts track each other (residual near zero). This is the conservation-of-events principle.

Value

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

time

Evaluation time.

n_risk

Number at risk (Kaplan-Meier).

n_event

Number of events at this time.

n_censor

Number censored at this time.

km_surv

Kaplan-Meier survival estimate.

km_cumhaz

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

par_surv

Parametric survival from the fitted model.

par_cumhaz

Parametric cumulative hazard.

cum_observed

Cumulative observed events to this time.

cum_expected

Cumulative expected events (sum of individual cumulative hazards for observations exiting the risk set).

residual

Expected minus observed (cum_expected - cum_observed).

For multiphase models, additional columns are appended for each phase: par_cumhaz_<phase>.

An attribute "summary" is attached with scalar diagnostics: total observed events, total expected events, and the final residual.

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)

# 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).

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

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

Usage

hzr_log1mexp(x)

Arguments

Value

Numeric vector with element-wise log(1 - exp(-x)).

Examples

hzr_log1mexp(c(0.01, 0.5, 5))

Numerically stable log(1 + exp(x))

Description

Numerically stable log(1 + exp(x))

Usage

hzr_log1pexp(x)

Arguments

Value

Numeric vector with element-wise log(1 + exp(x)).

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.

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)
)

## 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).

Phase types

"cdf"

Early risk that resolves over time. \Phi(t) = G(t), bounded [0, 1]. SAS equivalent: Early / G1 phase.

"hazard"

Late or aging risk that accumulates. \Phi(t) = -\log(1 - G(t)), monotone increasing. SAS equivalent: Late / G3 phase.

"constant"

Flat background hazard rate. \Phi(t) = t. No shape parameters are estimated. SAS equivalent: Constant / G2 phase.

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")

# 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")

Stepwise covariate selection for a parametric hazard model

Description

Run forward, backward, or two-way stepwise selection on an existing hazard fit using 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("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

Variable affected.

phase

Phase name (multiphase) or NA_character_.

criterion

"wald" or "aic".

score

Winning score used for the decision.

stat, df

Wald statistic and degrees of freedom.

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.

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.

criteria

Named list of the threshold / direction settings actually applied.

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.

Examples

data(avc)
avc <- na.omit(avc)
base <- hazard(survival::Surv(int_dead, dead) ~ age,
               data = avc, dist = "weibull", fit = TRUE)

sw <- hzr_stepwise(base, scope = ~ age + nyha,
                   data = avc, direction = "forward",
                   control = list(n_starts = 1))
print(sw)

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. Contains repeated thromboembolic events (up to 3 per patient) with left censoring, exercising the interval censoring likelihood.

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, 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,
  ...
)

Arguments

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.

Value

Numeric vector of predictions.

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)

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. 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), expected (expected event count from model predictions), observed_rate, expected_rate (events / n), chi_sq (per-group (O-E)^2/E contribution), 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.

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 hazard summary objects

Description

Formatted console display of summary.hazard() output: distribution, phase list (for multiphase), coefficient table with standard errors, and log-likelihood. 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.

Examples

data(avc)
avc <- na.omit(avc)
base <- hazard(survival::Surv(int_dead, dead) ~ age,
               data = avc, dist = "weibull", fit = TRUE)

sw <- hzr_stepwise(base, scope = ~ age + nyha,
                   data = avc, direction = "forward",
                   control = list(n_starts = 1))
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.

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 ------------------------------------
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, valves

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

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, or NA if the model has not been fitted or the covariance matrix is unavailable.

Examples

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