20  Hazard and nonparametric curves

20.1 When to use it

The survival chapter starts with patient-level follow-up and asks what the observed cohort experienced. This chapter starts one step later: you already hold a grid of fitted values and want to show what a model says between the observed event times. That distinction matters. A smooth line is a model summary, not a more precise Kaplan-Meier curve.

Name the quantity before you choose its column. Survival probability asks what proportion remains event-free beyond a time. Cumulative hazard asks how much event intensity has accumulated through that time. Hazard rate asks when events are occurring fastest among patients still at risk. They can be derived from one fitted model, but their scales, units, and clinical readings differ.

Reach for hv_hazard() when you have a parametric prediction grid and want to show the smooth curve, often with the empirical Kaplan-Meier points overlaid so a reviewer can judge how well the model fits. Reach for the nonparametric temporal-trend curve (hv_nonparametric()) when your outcome is a prevalence or a continuous measurement tracked over follow-up with a formula-based confidence band rather than a survival function. Reach for the ordinal curve (hv_ordinal()) when the outcome has graded levels (regurgitation grade, severity class) and you want one probability curve per grade. All three take pre-computed data and hand back an S3 object. Inspect that object first, then call plot() to get the bare ggplot you decorate. If the empirical hazard shape suggests distinct early, constant, or late phases, continue to fit an additive hazard model; these displays come before that modeling decision.

20.2 The data it needs

hv_hazard() takes a parametric prediction grid: one row per time point, with an estimate column (survival, hazard, or cumhaz) and explicit lower/upper CI columns that you name in the call. Two helpers generate the demonstration data. sample_hazard_data() builds the prediction grid; sample_hazard_empirical() builds the binned Kaplan-Meier overlay (the discrete empirical estimate the smooth curve is meant to track). Using the same cohort size for both keeps the overlay aligned with the curve.

A third helper, sample_hazard_cohort(), returns the subject-level records the overlay was computed from: one row per patient with a time and a status (1 = event, 0 = censored, the same coding as the survival chapter), rather than a curve. hv_atrisk() reads that indicator to decide who leaves the risk set, so a cohort coded the other way round silently counts censoring as events. Called with the same arguments it describes the same cohort, which is what lets the next section put an honest numbers-at-risk table under the figure.

dat_hp <- sample_hazard_data(n = 500, time_max = 10)
emp_hp <- sample_hazard_empirical(n = 500, time_max = 10, n_bins = 6)
coh_hp <- sample_hazard_cohort(n = 500, time_max = 10)

haz_surv <- hv_hazard(
  dat_hp,
  estimate_col  = "survival",
  lower_col     = "surv_lower",
  upper_col     = "surv_upper",
  empirical     = emp_hp,
  emp_lower_col = "lower",
  emp_upper_col = "upper"
)
haz_surv$meta
$x_col
[1] "time"

$estimate_col
[1] "survival"

$lower_col
[1] "surv_lower"

$upper_col
[1] "surv_upper"

$group_col
NULL

$has_ci
[1] TRUE

$n_obs
[1] 500

$emp_x_col
[1] "time"

$emp_estimate_col
[1] "estimate"

$emp_lower_col
[1] "lower"

$emp_upper_col
[1] "upper"

$emp_group_col
NULL

$emp_geom
[1] "point"

$ref_x_col
[1] "time"

$ref_estimate_col
[1] "survival"

$ref_group_col
NULL
head(haz_surv$data)
        time survival surv_lower surv_upper    hazard haz_lower haz_upper
1 0.01000000 99.99558   99.93731        100 0.6629126 0.3996474  1.099602
2 0.03002004 99.97702   99.84413        100 1.1485818 0.6924408  1.905203
3 0.05004008 99.95054   99.75561        100 1.4829116 0.8939968  2.459770
4 0.07006012 99.91808   99.66720        100 1.7546549 1.0578216  2.910523
5 0.09008016 99.88059   99.57770        100 1.9896233 1.1994760  3.300275
6 0.11010020 99.83868   99.48662        100 2.1996335 1.3260840  3.648628
       cumhaz cumhaz_lower cumhaz_upper
1 0.004419417            0    0.5871202
2 0.022986980            0    1.3519229
3 0.049470012            0    1.9990187
4 0.081954223            0    2.5912321
5 0.119483723            0    3.1493081
6 0.161453396            0    3.6834317

20.3 Build it

The most common figure is a smooth parametric survival curve with its confidence band plus discrete Kaplan-Meier points and error bars. We named the estimate and CI columns when we built haz_surv; now plot() uses those stored roles. The curve is on the percentage scale, so 80 means 80% surviving, not 0.80.

plot(haz_surv) +
  scale_colour_manual(values = c("steelblue"), guide = "none") +
  scale_fill_manual(values = c("steelblue"), guide = "none") +
  scale_x_continuous(limits = c(0, 10), breaks = 0:10) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20),
                     labels = function(x) paste0(x, "%")) +
  labs(x = "Years", y = "Survival (%)") +
  theme_hv_manuscript()
Figure 20.1: Smooth parametric survival curve with its confidence band and overlaid binned Kaplan-Meier points with error bars

20.4 Read it

The point of the overlay is the comparison between the smooth curve and the discrete points. Look for:

  • Points sitting on the curve. If the empirical Kaplan-Meier estimates fall on or very near the parametric line, the model is describing the data well. If the points wander systematically off the curve (above it early, below it late), the parametric form is the wrong shape and you should reconsider it.
  • The error bars against the band. The empirical error bars and the smooth confidence ribbon should tell a consistent story. Where the cohort thins, both widen; a point with a huge error bar is a region the figure cannot really support.
  • The curve shape itself. A survival curve falls and flattens; a hazard curve often peaks early then settles; a cumulative hazard only rises. A shape that contradicts the biology is a signal to check the model before the figure.

20.5 Numbers at risk under the curve

A survival figure without a risk table asks the reader to take the right-hand tail on trust. hv_atrisk() counts the risk set at times you choose and hv_atrisk_compose() stacks it under the curve on a matched axis.

The table needs subject-level records, which is why coh_hp exists: a curve cannot tell you how many patients stood behind each part of it.

curve_hp <- plot(haz_surv) +
  scale_colour_manual(values = c("steelblue"), guide = "none") +
  scale_fill_manual(values = c("steelblue"), guide = "none") +
  scale_x_continuous(limits = c(0, 10), breaks = seq(0, 10, 2)) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20),
                     labels = function(x) paste0(x, "%")) +
  labs(x = "Years", y = "Survival (%)") +
  theme_hv_manuscript()

risk_hp <- hv_atrisk(coh_hp, time = "time", status = "status",
                     report_times = seq(0, 10, 2))

hv_atrisk_compose(curve_hp, risk_hp)
Figure 20.2: The parametric survival curve and its Kaplan-Meier overlay with an aligned numbers-at-risk table, counted from the subject-level cohort the overlay was computed from

Be precise about what the counts describe, because two things are drawn here and they are not the same object. The counts and the Kaplan-Meier points both come from coh_hp: refit a Kaplan-Meier to that cohort and you recover the overlay exactly. The smooth curve is the parametric model, and it sits up to about two percentage points away from those points across this range. So the table describes the observed cohort, which is the relationship you want. A risk table always counts patients, never model.

That matters more than it sounds. The counts have to come from the data the figure is about; a table generated from some other cohort of the same size would look entirely plausible and be meaningless. If your own curve comes from a fitted model, take the risk table from the records you fitted it to.

20.6 Variations

20.6.1 Hazard rate

Switch estimate_col to "hazard" and the matching CI columns for the instantaneous event rate. The sample stores \(100h(t)\), so its units are events per 100 patient-years. This is the curve to show when the question is when events occur fastest rather than how many survive.

haz_rate <- hv_hazard(
  dat_hp,
  estimate_col = "hazard",
  lower_col    = "haz_lower",
  upper_col    = "haz_upper"
)

plot(haz_rate) +
  scale_colour_manual(values = c("firebrick"), guide = "none") +
  scale_fill_manual(values = c("firebrick"), guide = "none") +
  scale_x_continuous(limits = c(0, 10), breaks = 0:10) +
  scale_y_continuous(limits = c(0, 30)) +
  labs(x = "Years", y = "Events per 100 patient-years") +
  theme_hv_manuscript()
Figure 20.3: Instantaneous event rate per 100 patient-years, showing when post-operative events occur fastest among patients still at risk

20.6.2 Cumulative hazard

The "cumhaz" column equals \(100\{-\log[S(t)]\}\). Cumulative hazard \(H(t)\) is dimensionless; the helper multiplies it by 100 for display, but that scaling does not turn it into a percentage or a count per 100 patients. The local slope of \(100H(t)\) is the event rate per 100 patient-years. Reach for cumulative hazard when the question concerns accumulated event intensity rather than the surviving fraction.

haz_cumhaz <- hv_hazard(
  dat_hp,
  estimate_col  = "cumhaz",
  lower_col     = "cumhaz_lower",
  upper_col     = "cumhaz_upper"
)

plot(haz_cumhaz) +
  scale_colour_manual(values = c("darkorange"), guide = "none") +
  scale_fill_manual(values = c("darkorange"), guide = "none") +
  scale_x_continuous(limits = c(0, 10), breaks = 0:10) +
  labs(x = "Years", y = "Scaled cumulative hazard, 100 × H(t)") +
  theme_hv_manuscript()
Figure 20.4: Scaled cumulative hazard, 100 × H(t), a dimensionless quantity whose local slope is the event rate per 100 patient-years

20.6.3 Stratified by group

Pass group_col to compare two or more groups in one panel, each with its own curve, band, and empirical overlay. Look for curves whose bands stop overlapping in the windows where the treatment difference matters.

dat_strat <- sample_hazard_data(
  n = 400, time_max = 10,
  groups = c("No Takedown" = 1.0, "Takedown" = 0.65)
)
emp_strat <- sample_hazard_empirical(
  n = 400, time_max = 10, n_bins = 6,
  groups = c("No Takedown" = 1.0, "Takedown" = 0.65)
)

haz_strat <- hv_hazard(
  dat_strat,
  estimate_col  = "survival",
  lower_col     = "surv_lower",
  upper_col     = "surv_upper",
  group_col     = "group",
  empirical     = emp_strat,
  emp_lower_col = "lower",
  emp_upper_col = "upper"
)

plot(haz_strat) +
  scale_colour_manual(
    values = c("No Takedown" = "steelblue", "Takedown" = "firebrick"),
    name   = NULL
  ) +
  scale_fill_manual(
    values = c("No Takedown" = "steelblue", "Takedown" = "firebrick"),
    guide  = "none"
  ) +
  scale_x_continuous(limits = c(0, 10), breaks = 0:10) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20),
                     labels = function(x) paste0(x, "%")) +
  labs(x = "Years after Surgery", y = "Survival (%)") +
  theme_hv_manuscript() +
  theme(legend.position = "top")
Figure 20.5: Survival curves, bands, and empirical overlays for two groups in a single panel

20.6.4 Population life-table overlay

For age-stratified survival you often want to set the study curves against what the general population would do. Pass a life-table data frame to reference and set ref_group_col to draw the population survival as dashed lines per age group. The gap between a study curve and its dashed reference is the excess mortality attributable to the condition rather than to ageing.

dat_age <- sample_hazard_data(
  n = 600, time_max = 12,
  groups = c("<65" = 0.5, "65–80" = 1.0, "≥80" = 1.8)
)
emp_age <- sample_hazard_empirical(
  n = 600, time_max = 12, n_bins = 6,
  groups = c("<65" = 0.5, "65–80" = 1.0, "≥80" = 1.8)
)
lt <- sample_life_table(
  age_groups = c("<65", "65–80", "≥80"),
  age_mids   = c(55, 72, 85),
  time_max   = 12
)

haz_age <- hv_hazard(
  dat_age,
  estimate_col     = "survival",
  lower_col        = "surv_lower",
  upper_col        = "surv_upper",
  group_col        = "group",
  empirical        = emp_age,
  emp_lower_col    = "lower",
  emp_upper_col    = "upper",
  reference        = lt,
  ref_estimate_col = "survival",
  ref_group_col    = "group"
)

plot(haz_age) +
  scale_colour_manual(
    values = c("<65" = "steelblue", "65–80" = "forestgreen",
               "≥80" = "firebrick"),
    name   = "Age group"
  ) +
  scale_fill_manual(
    values = c("<65" = "steelblue", "65–80" = "forestgreen",
               "≥80" = "firebrick"),
    guide  = "none"
  ) +
  scale_linetype_manual(
    values = c("<65" = "dashed", "65–80" = "dashed", "≥80" = "dashed"),
    name = "Population reference"
  ) +
  scale_x_continuous(limits = c(0, 12), breaks = seq(0, 12, 2)) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20),
                     labels = function(x) paste0(x, "%")) +
  labs(x = "Years", y = "Survival (%)",
       caption = "Dashed lines: US population life table") +
  theme_hv_manuscript() +
  theme(legend.position = "top")
Figure 20.6: Age-stratified study survival curves set against dashed population life-table references

20.6.5 Nonparametric temporal-trend curve

When the outcome is a prevalence or a continuous measurement tracked over follow-up rather than a survival function, hv_nonparametric() prepares the average curve with a formula-based CI ribbon. For probability outcomes, sample_nonparametric_curve_data() computes an effective-sample-size standard error and applies a normal quantile on the logit scale. For continuous outcomes, it applies the normal quantile to a fixed residual standard error divided by the square root of the effective sample size. sample_nonparametric_curve_points() generates the binned summary points that overlay the curve. Build the S3 object once with the matching lower_col/upper_col and optional data_points, then plot().

curve_dat <- sample_nonparametric_curve_data(
  n            = 500,
  time_max     = 12,
  outcome_type = "probability",
  ci_level     = 0.68
)
pts_dat <- sample_nonparametric_curve_points(n = 500, time_max = 12)

np <- hv_nonparametric(
  curve_data  = curve_dat,
  lower_col   = "lower",
  upper_col   = "upper",
  data_points = pts_dat
)

The bare panel shows the average curve, its CI ribbon, and the binned summary points in default colours. Read it the same way you read the empirical overlay in Figure 20.1: the curve should pass through the binned points. The ribbon reflects the helper’s formula-based standard error, not a changing at-risk count or a bootstrap distribution. A flat curve where you expected a trend is a hint that the outcome_type or CI columns are misspecified.

plot(np)

The shaded ribbon here is a 68% normal-quantile interval, approximately one standard error on the scale used by the helper. Pass ci_level = 0.95 to sample_nonparametric_curve_data() for a 95% interval.

plot(np) +
  scale_colour_manual(values = c("steelblue"), guide = "none") +
  scale_fill_manual(values   = c("steelblue"), guide = "none") +
  scale_x_continuous(limits = c(0, 12), breaks = seq(0, 12, 2),
                     labels = function(x) paste(x, "yr")) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.2),
                     labels = scales::percent) +
  labs(x = "Follow-up (years)", y = "Prevalence (%)") +
  theme_hv_manuscript()
Figure 20.7: Nonparametric temporal-trend curve with its 68% formula-based normal-quantile interval and binned summary points

Pass group_col to the constructor to compare two average curves in a single panel; each group gets its own CI ribbon.

curve_grp <- sample_nonparametric_curve_data(
  n            = 400,
  time_max     = 7,
  groups       = c("Ozaki" = 0.7, "CE-Pericardial" = 1.3),
  outcome_type = "continuous",
  ci_level     = 0.68
)
pts_grp <- sample_nonparametric_curve_points(
  n = 400, time_max = 7,
  groups = c("Ozaki" = 0.7, "CE-Pericardial" = 1.3),
  outcome_type = "continuous"
)
np_grp <- hv_nonparametric(
  curve_data  = curve_grp,
  group_col   = "group",
  lower_col   = "lower",
  upper_col   = "upper",
  data_points = pts_grp
)

p_haz_grp <- plot(np_grp) +
  scale_colour_manual(
    values = c("Ozaki" = "steelblue", "CE-Pericardial" = "firebrick"),
    name   = NULL
  ) +
  scale_fill_manual(
    values = c("Ozaki" = "steelblue", "CE-Pericardial" = "firebrick"),
    guide  = "none"
  ) +
  scale_x_continuous(limits = c(0, 7), breaks = 0:7) +
  labs(x = "Follow-up (years)", y = "AV Peak Gradient (mmHg)") +
  theme_hv_manuscript()

# Rising curves leave the upper-left empty; pin the key there.
hv_legend_inside(p_haz_grp, prefer = "topleft")
Figure 20.8: Two nonparametric average curves compared in one panel, each with its own formula-based normal-quantile interval

20.6.6 Nonparametric ordinal-outcome curve

When the outcome has graded levels, hv_ordinal() prepares one probability curve per grade from a cumulative proportional-odds model (the prevalence of each regurgitation grade over time, say). The curve data is long format, one row per time by grade. sample_nonparametric_ordinal_data() generates the curves and sample_nonparametric_ordinal_points() the binned summary points.

ord_dat <- sample_nonparametric_ordinal_data(
  n = 800, time_max = 5,
  grade_labels = c("None", "Mild", "Moderate", "Severe")
)
ord_pts <- sample_nonparametric_ordinal_points(
  n = 800, time_max = 5,
  grade_labels = c("None", "Mild", "Moderate", "Severe")
)
head(ord_dat)
        time  estimate grade
1 0.01000000 0.6276258  None
2 0.01012532 0.6276898  None
3 0.01025221 0.6277545  None
4 0.01038069 0.6278200  None
5 0.01051078 0.6278864  None
6 0.01064250 0.6279535  None

Build the object with both curve data and summary points. Each line is one grade level; because every patient has exactly one grade at each time, the lines sum to roughly 1.0 at every time point. Read the figure as a stacked story: the "None" line starts high and declines as patients drift into worse grades, and the severe line creeps up. Colours run from grey for None through graduated severity to firebrick for Severe so the eye reads worsening as warming.

np_ord <- hv_ordinal(curve_data = ord_dat, data_points = ord_pts)
np_ord$meta
$x_col
[1] "time"

$estimate_col
[1] "estimate"

$grade_col
[1] "grade"

$n_obs
[1] 2000

$n_grades
[1] 4

$has_data_points
[1] TRUE
plot(np_ord) +
  scale_colour_manual(
    values = c("None" = "grey40", "Mild" = "steelblue",
               "Moderate" = "darkorange", "Severe" = "firebrick"),
    name = "AR Grade"
  ) +
  scale_x_continuous(limits = c(0, 5), breaks = 0:5) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.2),
                     labels = scales::percent) +
  labs(x = "Follow-up (years)", y = "Grade prevalence (%)") +
  theme_hv_manuscript() +
  theme(legend.position = c(0.75, 0.6))
Figure 20.9: One probability curve per regurgitation grade over time, coloured from grey to firebrick as severity worsens

To collapse adjacent grades, add their probabilities at each time before passing the data to the constructor. Do not simply drop the middle rows; then the displayed levels no longer account for every patient and will not sum to one. Here we combine None with Mild and Moderate with Severe.

ord_two <- transform(
  ord_dat,
  grade = ifelse(grade %in% c("None", "Mild"),
                 "None or Mild", "Moderate or Severe")
)
ord_two <- aggregate(estimate ~ time + grade, data = ord_two, FUN = sum)

np_ord_two <- hv_ordinal(curve_data = ord_two)

plot(np_ord_two) +
  scale_colour_manual(
    values = c("None or Mild" = "steelblue",
               "Moderate or Severe" = "firebrick"),
    name   = NULL
  ) +
  scale_x_continuous(limits = c(0, 5), breaks = 0:5) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.2),
                     labels = scales::percent) +
  labs(x = "Follow-up (years)", y = "Grade prevalence (%)") +
  theme_hv_manuscript() +
  theme(legend.position = "top")
Figure 20.10: Two-level ordinal version after summing adjacent grade probabilities; the two curves account for the full cohort at each time

20.7 Pitfalls

  • Mismatched CI columns. Every estimate column has its own pair of CI columns. Plot hazard with the survival bounds and the band will be nonsense. Name the matching lower_col/upper_col for the column you are plotting.
  • Forgetting the empirical overlay. A smooth parametric curve always looks convincing. Without the Kaplan-Meier points a reviewer cannot tell whether the model fits, so include empirical whenever you have it.
  • Ordinal lines that do not sum to one. If your grade curves drift above 1.0 or cross oddly, the wide-format grade columns were probably reshaped to long format incorrectly. Each time point’s grades should add to about 1.0.
  • Over-reading the CI level. The nonparametric ribbon defaults to a 68% normal-quantile interval, approximately one standard error on the helper’s scale. It is narrower than the 95% interval readers expect and is not a bootstrap interval. State the level and construction in the caption, or switch to ci_level = 0.95.