Postoperative risk rarely stays constant. It may be high immediately after an operation, settle to a low background rate after recovery, then rise again as the patient and the repair age. A Kaplan-Meier curve tells us how much survival remains. An additive hazard model asks the companion clinical question: which early, constant, or late process is contributing risk at each time?
Begin with the empirical curves. They are the evidence for choosing a parametric shape, not decoration added after the model is fit. If the empirical hazard falls sharply and then flattens, one monotone Weibull may be a useful summary. If it shows distinct regions, fit those phases explicitly and check the total curve against the empirical estimate.
TemporalHazard (Ehrlinger 2026) is the maintained R implementation of the Blackstone, Naftel, and Turner additive hazard model. It replaces the historical SAS/C HAZARD workflow in this recipe. We fit and inspect the R object here; the prediction chapter then turns it into curves.
Both examples use public data bundled with TemporalHazard. cabgkul contains 5,880 primary isolated CABG patients, with follow-up in months and a death indicator. avc contains 310 patients after atrioventricular canal repair and supports the two clinically identifiable phases used below. No patient extract is needed to run either recipe.
data(cabgkul, package = "TemporalHazard")
data(avc, package = "TemporalHazard")
avc <- na.omit(avc)
data.frame(
cohort = c("CABG", "AVC repair"),
patients = c(nrow(cabgkul), nrow(avc)),
deaths = c(sum(cabgkul$dead), sum(avc$dead)),
maximum_follow_up_months = c(max(cabgkul$int_dead), max(avc$int_dead))
) cohort patients deaths maximum_follow_up_months
1 CABG 5880 545 201.8276
2 AVC repair 305 68 170.5826hzr_kaplan() estimates survival and hzr_nelson() estimates cumulative hazard from the same subject-level records. Build both objects before choosing a distribution, then inspect the columns that will be plotted.
km <- hzr_kaplan(cabgkul$int_dead, cabgkul$dead)
na <- hzr_nelson(cabgkul$int_dead, cabgkul$dead)
head(as.data.frame(km)[c("time", "n_risk", "n_event", "survival",
"cl_lower", "cl_upper")]) time n_risk n_event survival cl_lower cl_upper
1 0.03285 5880 39 0.9933673 0.9909348 0.9951503
2 0.06571 5841 9 0.9918367 0.9891840 0.9938429
3 0.09856 5832 3 0.9913265 0.9886051 0.9934024
4 0.13142 5829 7 0.9901361 0.9872617 0.9923668
5 0.16427 5822 9 0.9886054 0.9855476 0.9910222
6 0.19713 5813 2 0.9882653 0.9851684 0.9907217head(as.data.frame(na)[c("time", "n_risk", "n_event", "cumhaz",
"cl_lower", "cl_upper")]) time n_risk n_event cumhaz cl_lower cl_upper
1 0.03285 5880 39 0.006632653 0.006305543 0.006972149
2 0.06571 5841 9 0.008173485 0.007711148 0.008656013
3 0.09856 5832 3 0.008687888 0.008122982 0.009281341
4 0.13142 5829 7 0.009888780 0.009236153 0.010574906
5 0.16427 5822 9 0.011434641 0.010703714 0.012201847
6 0.19713 5813 2 0.011778697 0.010978880 0.012620817The Kaplan-Meier panel shows the surviving fraction. The Nelson-Aalen panel turns the same follow-up into accumulated event intensity. Its local slope is the empirical clue to the hazard shape: changing slope means the instantaneous risk is changing with time.
ggplot(km, aes(time, survival * 100)) +
geom_ribbon(aes(ymin = cl_lower * 100, ymax = cl_upper * 100),
fill = "steelblue", alpha = 0.18) +
geom_step(colour = "steelblue", linewidth = 0.7) +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20),
labels = function(x) paste0(x, "%")) +
labs(x = "Months after CABG", y = "Freedom from death") +
theme_hv_manuscript()
The curve alone does not say how much cohort is left behind each part of it, and by the late follow-up here that is most of what a reader needs. hv_atrisk() counts the risk set at times you choose and hv_atrisk_compose() stacks it under the curve on a matched axis. It reads the same subject-level records hzr_kaplan() did, which is what makes the counts belong to this curve rather than to a cohort that happens to resemble it.
curve <- ggplot(km, aes(time, survival * 100)) +
geom_ribbon(aes(ymin = cl_lower * 100, ymax = cl_upper * 100),
fill = "steelblue", alpha = 0.18) +
geom_step(colour = "steelblue", linewidth = 0.7) +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20),
labels = function(x) paste0(x, "%")) +
scale_x_continuous(breaks = seq(0, 144, 48)) +
coord_cartesian(xlim = c(0, 168)) +
labs(x = "Months after CABG", y = "Freedom from death") +
theme_hv_manuscript()
risk_panel <- hv_atrisk(cabgkul, time = "int_dead", status = "dead",
report_times = seq(0, 144, 48))
hv_atrisk_compose(curve, risk_panel)
Read the two together. Freedom from death is around 60% at the right-hand edge, but 168 patients remain at twelve years against 5,880 at the start, so that part of the curve rests on a small and increasingly selected group. The widening confidence band says the same thing; the counts say it in a number a reviewer can check. Choose report_times that a clinical reader thinks in, and stop them where the cohort still supports an estimate rather than running them to the last observed event.
ggplot(na, aes(time, cumhaz)) +
geom_ribbon(aes(ymin = cl_lower, ymax = cl_upper),
fill = "firebrick", alpha = 0.18) +
geom_step(colour = "firebrick", linewidth = 0.7) +
labs(x = "Months after CABG", y = "Cumulative hazard, H(t)") +
theme_hv_manuscript()
Now fit the simple parametric description. The reviewed starting vector has a scale and shape value; it is an optimizer starting point, not a result to report. The estimates belong to the fitted object returned by hazard().
fit <- hazard(
time = cabgkul$int_dead,
status = cabgkul$dead,
theta = c(0.01, 1),
dist = "weibull",
fit = TRUE
)
summary(fit)hazard model summary
observations: 5880
predictors: 0
dist: weibull
engine: native-r-m2
converged: TRUE
log-lik: -3935.72
evaluations: fn=37, gr=9
Coefficients:
estimate std_error z_stat p_value
mu 0.0003613699 6.016407e-05 6.006407 1.896793e-09
nu 0.5944429760 2.346410e-02 25.334149 1.343733e-141data.frame(
converged = fit$fit$converged,
log_likelihood = fit$fit$objective,
function_evaluations = unname(fit$fit$counts["function"]),
gradient_evaluations = unname(fit$fit$counts["gradient"])
) converged log_likelihood function_evaluations gradient_evaluations
1 TRUE -3935.72 37 9The fitted Weibull shape is below 1, so its hazard declines with time. That is consistent with an early postoperative component, but the single curve cannot separate an early decline from a continuing background rate. Thus we treat it as a compact population summary and keep the empirical panels beside it as the fit check.
The AVC data have a clear early risk period and a much lower continuing rate. We represent those with an early CDF phase plus a constant phase. The early shape values below are the reviewed package example; fixed = "shapes" holds them fixed while the model estimates the scale of each phase. We do not invent new starting values.
fit_multiphase <- 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,
fixed = "shapes"
),
constant = hzr_phase("constant")
),
fit = TRUE,
control = list(n_starts = 3, maxit = 500)
)
summary(fit_multiphase)Multiphase hazard model (2 phases)
observations: 305
predictors: 0
dist: multiphase
phase 1: early - cdf (early risk)
phase 2: constant - constant (flat rate)
engine: native-r-m2
converged: TRUE
log-lik: -228.029
evaluations: fn=32, gr=10
Coefficients (internal scale):
Phase: early (cdf)
estimate std_error z_stat p_value
log_mu -1.4132735 0.1290435 -10.95192 6.50568e-28
log_t_half -0.6931472 NA NA NA
nu 1.0000000 NA NA NA
m 1.0000000 NA NA NA
Phase: constant (constant)
estimate std_error z_stat p_value
log_mu -7.609476 0.4495827 -16.92564 2.911483e-64Inspect convergence separately from the coefficient table. A populated summary is not enough: the optimizer must converge, and the covariance calculation must be positive definite for the free parameters.
data.frame(
converged = fit_multiphase$fit$converged,
positive_definite = fit_multiphase$fit$pd,
reciprocal_condition = fit_multiphase$fit$rcond,
log_likelihood = fit_multiphase$fit$objective
) converged positive_definite reciprocal_condition log_likelihood
1 TRUE TRUE 0.0784654 -228.0287theta <- coef(fit_multiphase)
check_times <- c(0.5, 12, max(avc$int_dead))
early_scale <- exp(theta["early.log_mu"])
constant_scale <- exp(theta["constant.log_mu"])
phase_contributions <- data.frame(
month = check_times,
early_hazard_per_month = early_scale * hzr_phase_hazard(
check_times, t_half = 0.5, nu = 1, m = 1, type = "cdf"
),
constant_hazard_per_month = constant_scale * hzr_phase_hazard(
check_times, type = "constant"
),
early_cumulative_hazard = early_scale * hzr_phase_cumhaz(
check_times, t_half = 0.5, nu = 1, m = 1, type = "cdf"
),
constant_cumulative_hazard = constant_scale * hzr_phase_cumhaz(
check_times, type = "constant"
)
)
phase_contributions month early_hazard_per_month constant_hazard_per_month
1 0.5000 1.216727e-01 0.0004957317
2 12.0000 7.787053e-04 0.0004957317
3 170.5826 4.157014e-06 0.0004957317
early_cumulative_hazard constant_cumulative_hazard
1 0.1216727 0.0002478658
2 0.2336116 0.0059487803
3 0.2426342 0.0845632012At 0.5 months, the evaluated early hazard contribution is much higher than the constant contribution on the same per-month scale. It then declines, while the constant contribution remains flat and continues to accumulate through late follow-up. The two phases overlap and add; the model does not switch from one to the other at 0.5 months. For this bounded early CDF phase, t_half = 0.5 means that half of the early phase’s eventual cumulative-hazard contribution has accumulated by 0.5 months.
A multiphase model can return numbers even when the data do not identify the requested phases. Stop and simplify the model when any of these occur:
fit$fit$converged is FALSE, especially after the iteration limit.fit$fit$pd is FALSE.Those are model warnings, not requests for a larger maxit. First ask whether the empirical curve contains the phase at all. Fixing reviewed shape parameters reduces the number of quantities the data must identify; adding an unsupported late phase does the opposite.