Package 'TemporalHazard'

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

object	A hazard object.
...	Unused; for S3 compatibility.

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

formula	Optional formula of the form Surv(time, status) ~ predictors. When provided, overrides direct time/status/x arguments and extracts from data. Example: hazard(Surv(time, status) ~ x1 + x2, data = df, dist = "weibull", fit = TRUE).
data	Optional data frame containing variables referenced in formula.
time	Numeric follow-up time vector.
status	Numeric or logical event indicator vector.
time_lower	Optional numeric lower bound vector for censoring intervals. Used when status == 2 (interval-censored); defaults to time if NULL.
time_upper	Optional numeric upper bound vector for censoring intervals. Used when status %in% c(-1, 2); defaults to time if NULL.
x	Optional design matrix (or data frame coercible to matrix).
time_windows	Optional numeric vector of strictly positive cut points for piecewise time-varying coefficients. When provided, each predictor column in x is expanded into one column per time window so each window gets its own coefficient.
theta	Optional numeric coefficient vector (starting values for optimization).
dist	Character baseline distribution label (default "weibull"). One of "weibull", "exponential", "loglogistic", "lognormal", or "multiphase". The single-distribution families differ in the shape the hazard traces over time; "multiphase" builds an additive N-phase hazard and requires phases. See the Baseline distributions section for what each means and when to use it.
phases	Optional named list of hzr_phase() objects specifying the phases for a multiphase model (dist = "multiphase"). See Examples.
fit	Logical; if TRUE, fit the model via maximum likelihood (default FALSE).
weights	Optional numeric vector of observation weights (non-negative). Each observation's log-likelihood contribution is multiplied by its weight. Use for severity-weighted repeated events. Default NULL (unit weights). Implements the SAS WEIGHT statement.
control	Named list of control options (see Details).
...	Additional named arguments retained for parity with legacy calling conventions.

Details

Control parameters:

• maxit: Maximum iterations (default 1000)

• n_starts: Number of optimization starts for multiphase fits (default 5). Each start after the first adds random noise to the initial values, drawn from the ambient RNG stream; call set.seed() before fitting for reproducible results.

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

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 \boldsymbol{\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").

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.

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

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

include_planned

Logical; if FALSE, only rows marked as implemented are returned.

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
)

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

Arguments

object	A fitted hazard object (with fit = TRUE).
n_boot	Integer: number of bootstrap replicates (default 200).
fraction	Numeric in (0, 1]: fraction of data to sample per replicate (default 1.0 for full bootstrap; < 1 for bagging).
seed	Optional integer random seed for reproducibility. When supplied, set.seed(seed) is called at function entry, jumping the global RNG to the seeded state; it is not restored on exit. Pass NULL (the default) to skip the set.seed() call and start from the caller's current RNG state. Note that the bootstrap consumes random numbers either way, so the global RNG state will advance during the call – seed = NULL avoids the reset at entry, not the advance during resampling.
verbose	Logical; if TRUE, print progress every 50 replicates.
x	An hzr_bootstrap object.
digits	Number of decimal places for formatting.
...	Additional arguments (ignored).

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

x	Numeric vector: the continuous covariate to calibrate.
event	Numeric vector: event indicator (1 = event, 0 = no event).
groups	Integer: number of quantile bins (default 10).
by	Optional factor or character vector for stratified calibration (SAS logitgr.sas functionality). If provided, calibration is computed within each stratum. Default NULL (no stratification).
link	Character: transform to apply to event probabilities. One of "logit" (default), "gompertz" (complementary log-log), or "cox".
time	Optional numeric vector: follow-up time, required when link = "cox". The Cox link computes $\log(\text{events} / \sum \text{time})$ (constant hazard rate).

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

p	Numeric vector of probabilities.
eps	Small positive tolerance in $(0, 0.5)$; default 1e-12.

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

time	Numeric vector of follow-up times.
event	Integer vector of event type indicators: 0 = censored, 1 = event type 1, 2 = event type 2, etc.
x	An hzr_competing_risks object.
digits	Number of decimal places for formatting.
...	Additional arguments (ignored).

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

object	A fitted hazard object (with fit = TRUE).
time	Numeric scalar: the horizon at which predicted survival is used to rank subjects into risk groups (e.g. time = 12 ranks by 12-month predicted survival). It does not restrict the event/expected counts, which are accumulated over each subject's full follow-up.
groups	Integer: number of risk groups (default 10 for deciles).
status	Optional numeric vector of event indicators (1 = event, 0 = censored). If NULL (default), extracted from the fitted object's stored data.
event_time	Optional numeric vector of observed event/censoring times. If NULL (default), extracted from the fitted object.

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

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

expected

Expected event count: the sum of each subject's predicted cumulative hazard at its own follow-up 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 at the horizon in the group.

mean_cumhaz

Mean predicted cumulative hazard at follow-up 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

time	Numeric vector of times (must be > 0).
t_half	Half-life: time at which $G(t_{1/2}) = 0.5$. Must be > 0.
nu	Time exponent controlling rate dynamics. SAS early: NU. SAS late: relates to GAMMA/ETA.
m	Shape exponent controlling the distributional form. SAS early: M. SAS late: relates to GAMMA/ALPHA.

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:

Case	Sign	Behavior
1	m > 0, nu > 0	Standard sigmoidal
1L	m = 0, nu > 0	Exponential-like (Weibull CDF)
2	m < 0, nu > 0	Heavy-tailed
2L	m < 0, nu = 0	Exponential decay
3	m > 0, nu < 0	Bounded cumulative
3L	m = 0, nu < 0	Bounded exponential

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

time	Numeric vector of times (must be > 0).
tau	Positive scalar scale parameter.
gamma	Positive scalar time exponent.
alpha	Non-negative scalar shape parameter (0 selects limiting case).
eta	Positive scalar outer exponent.

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

SAS name	R argument	Role
TAU	tau	Scale (time at which $(t/\tau) = 1$)
GAMMA	gamma	Power exponent on $t/\tau$
ALPHA	alpha	Shape (0 = exponential limiting case)
ETA	eta	Outer power exponent

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

object	A fitted hazard object (with fit = TRUE).
time_grid	Optional numeric vector of time points at which to evaluate the parametric model. If NULL (default), uses the sorted unique event times from the fitted data.

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

time	Numeric vector of follow-up times.
status	Numeric event indicator (1 = event, 0 = censored).
conf_level	Confidence level for the interval (default 0.95). The SAS default of 0.68268948 corresponds to a 1-SD interval.
event_only	Logical; if TRUE (default), only return rows at event times (where n_event > 0). If FALSE, return rows at all times reported by survival::survfit() (events and censoring times).

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

x	Numeric vector with positive values.

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

x	Numeric vector.

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

time	Numeric vector of follow-up times.
event	Numeric event indicator (1 = event, 0 = censored).
weight	Optional numeric vector of event weights (default 1). Weights are applied only to events (censored observations contribute zero weight). Use for severity-weighted repeated events.
conf_level	Confidence level for the interval (default 0.95).
x	An hzr_nelson object.
digits	Number of decimal places for formatting.
...	Additional arguments (ignored).

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

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

Arguments

type	Character; the phase's temporal shape — one of "cdf" (early resolving risk), "hazard" (accumulating G1 aging risk), "g3" (late rising risk, the original C/SAS late phase), or "constant" (flat background rate). See the Phase types section for what each means and when to use it.
t_half	Positive scalar; initial half-life (time at which $G(t_{1/2}) = 0.5$). Used for "cdf" and "hazard" phases. SAS early: THALF/RHO.
nu	Numeric scalar; initial time exponent. Used for "cdf" and "hazard" phases. SAS early: NU.
m	Numeric scalar; initial shape exponent. Used for "cdf" and "hazard" phases. SAS early: M.
tau	Positive scalar; scale parameter for "g3" phases. SAS late: TAU.
gamma	Positive scalar; time exponent for "g3" phases. SAS late: GAMMA.
alpha	Non-negative scalar; shape parameter for "g3" phases. When alpha > 0, the generic G3 formula is used; alpha = 0 gives the exponential limiting case. SAS late: ALPHA.
eta	Positive scalar; outer exponent for "g3" phases. SAS late: ETA.
formula	Optional one-sided formula (e.g. ~ age + nyha) for phase-specific covariates. When NULL (default), the phase inherits the global formula from hazard().
fixed	Character vector naming shape parameters to hold fixed during optimization. Valid names for "cdf"/"hazard": "t_half", "nu", "m", or "shapes" (shorthand for all three). Valid names for "g3": "tau", "gamma", "alpha", "eta", or "shapes" (shorthand for all four). Fixed parameters are held at their starting values; only mu (and covariates) are estimated. Ignored for "constant" phases. This mirrors the SAS/C HAZARD workflow where shapes are typically fixed and only scale parameters are estimated.
x	An hzr_phase object (for print.hzr_phase()).
...	Additional arguments (ignored).

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

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 \boldsymbol{\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")

# 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

time	Numeric vector of times (> 0).
t_half	Half-life parameter (> 0).
nu	Time exponent.
m	Shape parameter.
type	Phase type: "cdf" (early – uses $G(t)$), "hazard" (late – uses cumulative hazard from $h(t)$), or "constant" (flat rate – $\Phi = t$).

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

time	Numeric vector of times (> 0).
t_half	Half-life parameter (> 0).
nu	Time exponent.
m	Shape parameter.
type	Phase type: "cdf" (early – uses $G(t)$), "hazard" (late – uses cumulative hazard from $h(t)$), or "constant" (flat rate – $\Phi = t$).

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

fit	A fitted hazard object built via the ⁠formula = Surv(...) ~ predictors, data = df⁠ interface.
scope	Candidate set. NULL (default) uses every data-frame column not already in the model for every phase. For single-distribution fits, pass a one-sided formula (~ age + nyha) or a character vector of names. For multiphase fits, pass a named list of one-sided formulas keyed by phase.
data	Data frame the base fit was built on. Required for refits.
direction	Search strategy — one of "both" (default), "forward", or "backward". Controls whether variables may only enter, only leave, or both. See the Selection direction and criterion section.
criterion	Entry / retention rule — one of "wald" (default) or "aic". "wald" applies SAS-style p-value thresholds (slentry / slstay); "aic" adds or drops whenever it lowers the AIC. See the Selection direction and criterion section.
slentry	Entry p-value threshold for the Wald criterion. Default 0.30 matches SAS SLENTRY.
slstay	Retention p-value threshold for the Wald criterion. Default 0.20 matches SAS SLSTAY.
max_steps	Hard cap on total accepted actions. Emits a warning() if hit. Default 50.
max_move	Per-variable oscillation cap. When a variable has entered + exited more than max_move times it is frozen for the remainder of the run. Default 4.
force_in	Character vector of variables that must remain in the model. Such variables are still scored and reported in the selection trace, but are never dropped.
force_out	Character vector of variables that may never be considered as candidates.
trace	Logical; print step-by-step progress to the console. Default TRUE.
...	Unused.
x	An hzr_stepwise object.
object	An hzr_stepwise object.

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.

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 full candidate model and only drop variables — the weakest term leaves each step until all survivors clear the retention rule.

direction = "both" (default)

Two-way stepwise: after each entry, already-selected variables are re-tested and may be dropped. This is the SAS SELECTION = STEPWISE strategy. max_move caps how often a single variable may oscillate before it is frozen.

criterion = "wald" (default)

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. Note this differs algorithmically from C/SAS HAZARD, which selects on a score (Q) statistic evaluated without refitting; the two can take different step paths 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.

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)

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

x	Object to test.

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,
  conf.type = c("log-log", "logit"),
  ...
)

Arguments

object	A hazard object.
newdata	Optional matrix or data frame of predictors. For types requiring time (e.g., "survival", "cumulative_hazard"), newdata should include a time column, or time will be taken from the fitted object's data.
type	Prediction type: "linear_predictor": Linear predictor eta = x*beta (not available for multiphase) "hazard": Instantaneous hazard. Single-distribution models return the hazard scale exp(eta); multiphase models return the additive hazard h(t|x) = sum_j mu_j(x) phi_j'(t) and so require time values (like "survival"/"cumulative_hazard"). decompose is not supported for "hazard". "survival": Survival probability S(t|x) = exp(-H(t|x)) "cumulative_hazard": Cumulative hazard H(t|x) at event times
decompose	Logical; if TRUE and the model is multiphase, return a data frame with per-phase cumulative hazard contributions alongside the total. Ignored for single-distribution models. Default FALSE.
se.fit	Logical; if TRUE, compute delta-method standard errors and confidence limits for each prediction. The return value becomes a data frame with columns fit, se.fit, lower, upper. Default FALSE. CLs are computed on the log-hazard / log-cumhaz scale and on the log(-log(survival)) scale so lower/upper stay inside the valid range of each prediction type; linear_predictor uses symmetric natural-scale CLs. For multiphase models, se.fit = TRUE combines with decompose = TRUE when type = "cumulative_hazard": the result is a long data frame with one row per prediction time and component (component in "total" plus each phase name) and columns fit, se.fit, lower, upper. Per-phase CLs use only that phase's parameters, so they do not sum to the total CL. The combination is not available for type = "survival" (per-phase survival is not additive).
level	Numeric confidence level in ⁠(0, 1)⁠; default 0.95. Only used when se.fit = TRUE.
conf.type	Transform for type = "survival" confidence limits when se.fit = TRUE: "log-log" (default) builds them on ⁠log(-log S)⁠ (the survival::survfit standard); "logit" builds them on logit(1 - S), reproducing SAS HAZARD's HAZPRED survival limits. Other types are unaffected (hazard/cumulative-hazard use a log scale that already matches HAZPRED). Only used when se.fit = TRUE.
...	Unused; included for S3 compatibility.

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

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

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 hzr_calibrate

Description

Print method for hzr_calibrate

Usage

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

Arguments

x	An hzr_calibrate object.
digits	Number of decimal places for formatting.
...	Additional arguments (ignored).

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

x	An hzr_deciles object.
digits	Number of decimal places for formatting.
...	Additional arguments (ignored).

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

x	An hzr_gof object.
digits	Number of decimal places for formatting.
...	Additional arguments (ignored).

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

x	An hzr_kaplan object.
digits	Number of decimal places for formatting.
n	Maximum rows to print (default 20).
...	Additional arguments (ignored).

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

---

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

fit	An hzr_stepwise object.

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)

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

object	A hazard object.
...	Unused; for S3 compatibility.

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

object	A hazard object.
...	Unused; for S3 compatibility.

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 requires numDeriv and an invertible Hessian; if either is unavailable the fit emits a warning and the conserved log_mu stays NA (the rest of the matrix is unaffected). Returns a scalar NA only when the model has not been fitted or no covariance matrix is available.

Examples

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

---