Package 'prsladder'

Assess the Calibration of a Ladder

Description

The function compares the risks that each model of a ladder predicts with the risks observed in a cohort, at each horizon. The function makes the comparison for all persons together, along a smooth calibration curve, and within groups of predicted risk. Use it to find if a model predicts risks that are too high, too low, or too extreme for new persons. The PRS reporting standards ask for calibration in all analyses, and say that it is most relevant for persons outside the development sample (Wand et al. 2021).

The function has three modes. Without newdata and validation, it evaluates the development persons, and the estimates are apparent. With validation, it corrects the O/E ratio, the calibration slope, the smooth curve, and the summaries of the curve for optimism with the bootstrap. With newdata, it evaluates other persons. validation explains when to use each mode.

Usage

assess_calibration(
  ladder,
  horizons,
  groups = 10L,
  newdata = NULL,
  validation = NULL,
  min_events = 10L,
  min_at_risk = 50L,
  min_censoring_survival = 0.02
)

Arguments

ladder	(prs_ladder or prs_export) A ladder from fit_ladder() or its export from export_ladder(). An export holds no persons, and it needs newdata.
horizons	(numeric(); range ⁠(0, Inf)⁠) Distinct finite times since entry, in the units of the age columns.
groups	(integer(1); default 10L; range ⁠[1, Inf)⁠) Number of quantile groups of predicted risk. A value of 1L puts all persons in one group. Tied risks can give a model fewer groups (see Details).
newdata	(prs_cohort or NULL; default NULL) A different cohort, such as another biobank, a later period, or another ancestry group. Build it with prs_cohort(data, reference = ladder$cohort), or with prs_cohort(data, reference = ladder) for an export. Its scores, covariates, and index must be those of the development cohort. If NULL, the function evaluates the development persons. The persons must be separate from the development persons and from the GWAS data that gave the scores (Choi et al. 2020). For a ladder, the function gives a warning with the number of persons of newdata whose ids are also in ladder$cohort. Give newdata or validation, not both.
validation	(prs_bootstrap_plan or NULL; default NULL) A plan from bootstrap_plan() for the optimism correction of the estimates on the development persons (see the section "Optimism correction"). If NULL, the function gives the apparent or the external estimates.
min_events	(integer(1); default 10L; range ⁠[1, Inf)⁠) Minimum number of T2 events at or before a horizon.
min_at_risk	(integer(1); default 50L; range ⁠[1, Inf)⁠) Minimum number of persons whose time since entry is at least the horizon.
min_censoring_survival	(numeric(1); default 0.02; range ⁠(0, 1]⁠) Minimum censoring survival $G(h)$ at a horizon $h$. A value of 1 requires no censoring at or before the horizon.

Details

Persons and windows

At horizon $h$, the function evaluates only the persons with a prediction (see predict_risk()). The column n_removed counts the other persons. The horizon guard of horizon_support() applies to the evaluated persons, not to the full cohort. The function leaves out each horizon that fails the guard and gives one warning for all of them.

The window of each person starts at entry and stops at min(stop, start + h). A second event is an observed event only when it occurs in the window.

Calibration-in-the-large

The function compares the observed and expected numbers of events in the windows. This comparison adapts the observed-to-expected comparison of Crowson et al. (2016), which uses the full follow-up of each person, to windows that stop at $h$:

• The observed count $O$ is the number of events in the windows. It is the same for each model.

• The expected count $E$ is the sum over persons of the cumulative hazard of the model over the window. The cumulative hazard uses the Breslow baseline hazard of the fit. With this baseline, the sum of $E$ over the development persons equals $O$ also when event times are tied. The Efron baseline of predict_risk() gives a larger sum at tied event times. That baseline estimates the hazard of the continuous time that the ties hide. expected_risk keeps that baseline, because on the probability scale it matches the Kaplan-Meier observed_risk. Without tied event times, the two baselines are equal.

• The ratio oe_ratio is $O/E$. A ratio above one means that the model predicts too few events.

• The 95% confidence interval treats $O$ as Poisson, $O/E \times \exp(\pm 1.96 / \sqrt{O})$. The interval is on the log scale. oe_ratio_se is the delta-method standard error of the ratio on the ratio scale, $(O/E) / \sqrt{O}$. Because the guard requires at least min_events events, $O$ is positive and the interval is finite.

The expected count starts at the entry of each person. Because of this, the ratio needs no assumption about censoring other than the assumption of the Cox model. On the development persons, the ratio is one under the diagnosis index, because every window stops at $h$ on the time scale of the model. Under the recruitment index, the windows stop at different times, and the ratio is close to one. The apparent interval is still useful. It covers the O/E ratio of the fitted model in the population of the development persons. In a simulation of 200 development samples, it covered that ratio in 93% to 98% of the samples. Use a plan in validation, or a different cohort in newdata, to get an informative ratio.

With competing events (see prs_cohort()), the window of a person with a competing event stops at that event. The expected count then uses the cause-specific hazard of T2 only. Because of this, the ratio measures the calibration of the model of T2, and it does not use the competing model. The calibration slope below also measures the model of T2.

Risk scale

The risk-scale columns describe the same persons:

• expected_risk is the mean predicted risk from predict_risk().

• observed_risk is one minus the Kaplan-Meier survival at $h$ on time since entry.

• observed_risk_se is the Greenwood standard error of that estimate.

With competing events, expected_risk is the mean predicted cumulative incidence, which uses the models of both events. observed_risk is then the Aalen-Johansen cumulative incidence of T2 at $h$ (Aalen and Johansen 1978), and observed_risk_se is its infinitesimal jackknife standard error. The comparison of these two risks within groups of predicted risk is the calibration check of Wolbers et al. (2009).

The observed risk ignores the predictors. It assumes that censoring on time since entry is independent of the scores and the covariates.

Under the recruitment index, the incident and the prevalent persons can have different censoring on time since entry (see assess_discrimination()). When the evaluated persons hold both values of t1_status, the observed risk mixes the risks of the two values. Let $\hat F_s(h)$ be the Kaplan-Meier or Aalen-Johansen risk of the $n_s$ persons with value $s$ of t1_status, and $\mathrm{SE}_s$ its standard error. The observed risk and its standard error are

$\begin{aligned} \hat F &= \sum_s \frac{n_s}{n} \hat F_s(h), \\ \mathrm{SE}(\hat F) &= \sqrt{\sum_s \frac{n_s^2}{n^2} \mathrm{SE}_s^2 + \frac{1}{n} \sum_s \frac{n_s}{n} \left(\hat F_s(h) - \hat F\right)^2}. \end{aligned}$

The first term of the variance treats the sizes $n_s$ as constants. The second term is the variance of the shares $n_s / n$. When the mixture is not NA, it equals the mean of the censoring-weighted events, with one censoring distribution for each value of t1_status. It assumes that censoring is independent of the entry within each value.

Groups of predicted risk

For each model, the function puts the persons into groups quantile groups of predicted risk. Group 1 holds the lowest risks. In each group, observed_risk is the observed risk of the persons in the group, from the same estimator as for all persons. Without competing events, this comparison is the grouped Kaplan-Meier calibration check of Crowson et al. (2016). It is zero for a group with no events in the windows. It is NA when nobody in the group is at risk at $h$ and a person is censored at the last time, because no estimate then exists at $h$. A group whose persons all have T2 or a competing event before $h$ keeps its estimate. Without competing events, this estimate is one. Under the recruitment index, the observed risk of a group mixes the risks of its incident and its prevalent persons. The weights of the mixture are their shares in the group. It is then also NA when this rule gives NA for the incident or the prevalent persons of the group. The same rules apply to the observed risk of all persons.

Tied risks can make two or more quantiles equal, for example in a model with only a few distinct predicted risks. The groups between equal quantiles are then empty, and the function removes them. Persons with the same risk are always in the same group. That model has fewer than groups groups at that horizon. The function gives one warning that names each such model and horizon. The other models keep groups groups.

Smooth calibration curve

The smooth curve gives the observed risk as a function of the predicted risk (Austin et al. 2020). Let $\hat p_i$ be the predicted risk of person $i$ at $h$, and $x_i = \log(-\log(1 - \hat p_i))$. The function fits a Cox model of the windows on time since entry, with a natural spline of $x$. The spline has an inner knot at the 50th percentile of $x$ and outer knots at the 10th and 90th percentiles. This spline is the restricted cubic spline with 3 knots of Austin et al. (2020). The observed risk at $\hat p$ is

$1 - \exp\{-\hat H(h)\, e^{\hat g(x)}\}.$

Here $\hat H$ is the Efron baseline cumulative hazard and $\hat g$ is the spline of this model. Austin et al. (2020) fit the model to the full follow-up and take the curve at $h$. This function fits it to the windows that stop at $h$, as the O/E ratio and the slope do. The proportional-hazards assumption of the model then must hold only up to $h$.

Under the recruitment index with both values of t1_status, the Cox model has one baseline $\hat H_s$ for each value $s$. A logistic model of t1_status on the same spline of $x$ gives the probability $\hat P(s \mid x)$ of each value. The observed risk at $\hat p$ is then the mixture

$\sum_s \hat P(s \mid x) \left[1 - \exp\{-\hat H_s(h)\, e^{\hat g(x)}\}\right].$

This mixture is the observed risk of all persons with the predicted risk $\hat p$. The predicted risk can separate the incident and the prevalent persons, for example when the hazard rises with time since T1. The logistic model then gives $\hat P(s \mid x)$ close to 0 or 1. At such an $x$, the curve uses the baseline of one value of t1_status. The ranges of the predicted risks of the two values can then have a gap. Between the two ranges, the curve has no persons, and its shape there comes from the limit of the logistic fit.

Three summaries describe the distance between the curve and the diagonal. Let $d_i$ be the absolute difference between the observed risk at $\hat p_i$ and $\hat p_i$. ici is the mean of $d_i$ over the evaluated persons, the integrated calibration index. e50 and e90 are the median and the 90th percentile of $d_i$. The summaries use the observed risk of each person, not a grid. The element curve gives the curve at 200 risks evenly spaced from the 1st to the 99th percentile of the predicted risks.

The curve and its summaries are NA in these cases:

• The cohort declares competing events. The decile groups stay available.

• A predicted risk is 0 or 1, or the three knots are not distinct, for example when a model has few distinct risks.

• The Cox fit of the curve gives a warning. The function then gives one warning for all such models and horizons.

• A coefficient of the Cox fit or of the logistic fit is NA. The function gives no warning.

• The logistic fit of t1_status stops with an error, or gives a warning other than the warnings of separation. The function then gives one other warning for all such models and horizons.

The element curve has no rows for such a model and horizon. The summaries have no standard error or interval outside the corrected mode, and their ⁠_se⁠, ⁠_lower⁠, and ⁠_upper⁠ columns are then NA.

Calibration slope

For an external cohort, the calibration slope is the coefficient of the development linear predictor $\eta$ in a Cox model fitted to the windows of newdata:

• A slope of one means that the spread of the predictions is correct.

• A slope below one means that the predictions are too extreme.

• A slope above one means that the predictions are too moderate.

This slope is the Cox model form of the calibration slope of Crowson et al. (2016). The Cox model has the time scale of the fit, and with design = "index" in fit_ladder() also its strata. Its interval is the Wald interval, the slope plus or minus 1.96 standard errors. The slope is NA when the linear predictor is constant over the persons. When its Cox fit gives a warning, the slope is NA, and the function gives one warning for all such models and horizons. The slope is precise only when the linear predictor varies across persons. Its standard error is about $1 / (\mathrm{SD}(\eta) \sqrt{O})$. When the terms of a model have no effect, as can happen for base, the development coefficients are noise. The true slope is then zero, and the interval can contain both zero and one. On the development persons, the slope does not show overfitting. Over all follow-up, the slope is exactly one, because the model comes from these persons. When the hazards are proportional, a window that stops at $h$ moves the slope away from one. But its mean stays near one, also for an overfit model. In a simulation of an overfit model, the mean slope at 5 years was 0.97 on the development persons and 0.37 on new persons. Because of this, the slope columns are NA for apparent evaluation.

Optimism correction

With a plan from bootstrap_plan() in validation, the function corrects the measures at each horizon for optimism with the bootstrap (Harrell Jr. et al. 1996). These measures are the O/E ratio, the slope, ICI, E50, E90, and the curve. bootstrap_plan() describes the method, the random number streams, and the failed replicates. The column evaluation of models is then "optimism-corrected":

• The function corrects the O/E ratio on the log scale. oe_ratio is $\exp(a - o)$, where $a$ is the apparent log ratio and $o$ the mean optimism of the log ratio. The rows of optimism have the measure "log_oe_ratio".

• slope is the apparent slope at the horizon minus the mean optimism. Because the window moves the apparent slope away from one, the corrected slope at a horizon holds this noise.

• ici, e50, and e90 are each the apparent value minus its mean optimism, truncated at 0. Because the absolute value is not linear, these values differ from the summaries of the corrected curve.

• The corrected curve is the apparent curve minus the mean optimism at each of its 200 points, as in rms::calibrate(). Each replicate gives two curves of the refit at each point. The optimism at the point is the observed risk on the bootstrap sample minus the observed risk on the development persons. The element curve holds the apparent curve with the label "apparent" and the corrected curve with the label "optimism-corrected", on the same points. The optimism of the curve is not in optimism. At a point with fewer than half of the replicates, the corrected curve is NA.

• The element shrinkage gives the corrected slope over all follow-up. Its apparent value is 1, because each model comes from the development persons. This slope is the shrinkage factor of the linear predictor that new persons need (van Houwelingen and le Cessie 1990; Harrell Jr. et al. 1996). Report it as the measure of overfitting. The column heuristic_shrinkage of fit_ladder()$table gives an estimate without the bootstrap.

The function does not correct the other columns. observed_events, observed_risk, and observed_risk_se do not depend on a model, and expected_events and expected_risk keep their apparent values. The element groups holds the apparent groups, with the label "apparent".

A replicate in which the Cox fit of a curve fails gives no optimism for that curve and its summaries. The other measures of the replicate stay. The column n_boot_used of optimism counts the replicates used for each row.

The standard errors and the intervals of oe_ratio, slope, ici, e50, e90, and shrinkage depend on interval of the plan:

• "analytic" and "none": NA. The refit standard error of the slope is too large for the corrected slope. A Poisson interval of the corrected O/E ratio is not yet tested.

• "two-stage": the standard error is the standard deviation of the outer corrected estimates, and the limits are their 2.5% and 97.5% quantiles. For the O/E ratio, the quantiles are on the log scale, and the standard error is on the ratio scale. The limits of ici, e50, and e90 are truncated at 0. The curve has no interval.

Value

A list with the six elements below. Because each element has the same columns in every mode, the results of two modes stack with base::rbind(). The rows are in the order of horizon, then model as in ladder$models, then group or point of the curve. When no horizon passes the guard, models, groups, and curve have zero rows.

models

One row per horizon and model, with these columns:

model

Model name.

horizon

Horizon $h$.

evaluation

"apparent", "optimism-corrected" with validation, or "external" with newdata.

n

Number of persons evaluated.

n_removed

Number of persons without a prediction.

observed_events

Observed count $O$.

expected_events

Expected count $E$.

oe_ratio, oe_ratio_se, oe_ratio_lower, oe_ratio_upper

Ratio $O/E$, its standard error, and its 95% confidence interval.

expected_risk

Mean predicted risk.

observed_risk, observed_risk_se

Kaplan-Meier risk, or Aalen-Johansen cumulative incidence with competing events, and its standard error. Under the recruitment index, the mixture of these risks over the values of t1_status (see Details). NA when nobody is at risk at $h$ and a person is censored at the last time, for all persons or, under the recruitment index, for the persons of one value of t1_status.

slope, slope_se, slope_lower, slope_upper

Calibration slope, its standard error, and its 95% confidence interval, NA for apparent evaluation.

ici, e50, e90

Mean, median, and 90th percentile of the absolute difference between the observed risk of the curve and the predicted risk. Each has the columns ⁠_se⁠, ⁠_lower⁠, and ⁠_upper⁠, which are NA except with a two-stage plan.

groups

One row per horizon, model, and group, with these columns:

model, horizon, evaluation

As in models.

group

Group number, from 1 for the lowest risks. A model with tied risks can have fewer than groups groups.

n

Number of persons in the group.

n_events

Number of T2 events in the windows of the group.

expected_risk

Mean predicted risk of the group.

observed_risk, observed_risk_se

Observed risk of the group and its standard error.

curve

The smooth calibration curve at 200 points of each model and horizon, with the columns model, horizon, evaluation, predicted, and observed. With validation, each model and horizon has two curves on the same points, first the apparent curve and then the corrected curve. It has zero rows with competing events.

shrinkage

One row per model with the corrected slope over all follow-up, and the columns model, slope, slope_se, slope_lower, and slope_upper. It has zero rows without validation.

optimism

One row per corrected estimate, with the columns of the element optimism of assess_discrimination(). The rows are in the order of horizon, then the measures "log_oe_ratio", "slope", "ici", "e50", and "e90", then model. The rows of shrinkage follow, with the measure "shrinkage", the horizon NA, and the apparent value 1. It has zero rows without validation.

settings

A list with the elements of the element settings of assess_discrimination().

References

Aalen OO, Johansen S (1978). “An empirical transition matrix for non-homogeneous Markov chains based on censored observations.” Scandinavian Journal of Statistics, 5(3), 141–150.

Austin PC, Harrell Jr. FE, van Klaveren D (2020). “Graphical calibration curves and the integrated calibration index (ICI) for survival models.” Statistics in Medicine, 39(21), 2714–2742. doi:10.1002/sim.8570.

Choi SW, Mak TS, O'Reilly PF (2020). “Tutorial: a guide to performing polygenic risk score analyses.” Nature Protocols, 15(9), 2759–2772. doi:10.1038/s41596-020-0353-1.

Crowson CS, Atkinson EJ, Therneau TM (2016). “Assessing calibration of prognostic risk scores.” Statistical Methods in Medical Research, 25(4), 1692–1706. doi:10.1177/0962280213497434.

Harrell Jr. FE, Lee KL, Mark DB (1996). “Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.” Statistics in Medicine, 15(4), 361–387. doi:10.1002/(SICI)1097-0258(19960229)15:4<361::AID-SIM168>3.0.CO;2-4.

van Houwelingen JC, le Cessie S (1990). “Predictive value of statistical models.” Statistics in Medicine, 9(11), 1303–1325. doi:10.1002/sim.4780091109.

Wand H, Lambert SA, Tamburro C, Iacocca MA, O'Sullivan JW, Sillari C, Kullo IJ, Rowley R, Dron JS, Brockman D, Venner E, McCarthy MI, Antoniou AC, Easton DF, Hegele RA, Khera AV, Chatterjee N, Kooperberg C, Edwards K, Vlessis K, Kinnear K, Danesh JN, Parkinson H, Ramos EM, Roberts MC, Ormond KE, Khoury MJ, Janssens ACJW, Goddard KAB, Kraft P, MacArthur JAL, Inouye M, Wojcik GL (2021). “Improving reporting standards for polygenic scores in risk prediction studies.” Nature, 591(7849), 211–219. doi:10.1038/s41586-021-03243-6.

Wolbers M, Koller MT, Witteman JCM, Steyerberg EW (2009). “Prognostic models with competing risks: methods and application to coronary risk prediction.” Epidemiology, 20(4), 555–561. doi:10.1097/EDE.0b013e3181a39056.

See Also

horizon_support() to check horizons before the evaluation, validation for the three modes, bootstrap_plan() for the optimism correction, assess_discrimination(), and assess_thresholds().

Examples

#----------------------------------------------------------------------------
# assess_calibration() example
#----------------------------------------------------------------------------
# Fit a ladder on a development cohort.
dat <- simulate_prs_data(n = 2000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(onset = "onset_score", progression = "progression_score")
)
ladder <- fit_ladder(cohort)

# Build an external cohort with the rules of the development cohort.
external <- prs_cohort(
  simulate_prs_data(n = 2000L, seed = 2L, id_prefix = "E"),
  reference = cohort
)

# Assess calibration at 5 and 10 years in 5 groups of predicted risk.
calibration <- assess_calibration(
  ladder,
  horizons = c(5, 10),
  groups = 5L,
  newdata = external
)
calibration$models[c("model", "horizon", "oe_ratio", "slope", "ici")]

# Observed and expected risk by group for the model with both scores.
calibration$groups[calibration$groups$model == "onset+progression", ]

#----------------------------------------------------------------------------
# Optimism correction on the development persons
#----------------------------------------------------------------------------
# A small number of replicates keeps the example fast. An analysis needs
# 200 or more.
corrected <- assess_calibration(
  ladder,
  horizons = 10,
  validation = bootstrap_plan(n_boot = 20L, seed = 1L)
)
corrected$models[c("model", "oe_ratio", "slope", "ici", "e90")]
corrected$shrinkage

---

Assess the Discrimination of a Ladder

Description

The function reports two measures of predictive accuracy for each model of a ladder at each horizon, with standard errors and 95% confidence intervals. The time-dependent AUC measures discrimination, which is how well the predicted risks rank persons with an event above persons without one. The Brier score measures overall accuracy, which is the mean squared error of the predicted risks. The function also reports the difference in each measure between pairs of models, to show the gain from added scores.

The function has three modes. Without newdata and validation, it evaluates the development persons, and the estimates are apparent. With validation, it corrects the apparent estimates for optimism with the bootstrap. With newdata, it evaluates other persons. validation explains when to use each mode.

Usage

assess_discrimination(
  ladder,
  horizons,
  pairs = NULL,
  newdata = NULL,
  validation = NULL,
  min_events = 10L,
  min_at_risk = 50L,
  min_censoring_survival = 0.02
)

Arguments

ladder	(prs_ladder or prs_export) A ladder from fit_ladder() or its export from export_ladder(). An export holds no persons, and it needs newdata.
horizons	(numeric(); range ⁠(0, Inf)⁠) Distinct finite times since entry, in the units of the age columns.
pairs	(list of character(2) or NULL; default NULL) Pairs of models to compare, each as c(full, reduced) model names. The names must be in names(ladder$models) and must differ. Pairs need not be nested, and a pair must not occur twice. If NULL, the function compares every nested pair in the order of compare_ladder().
newdata	(prs_cohort or NULL; default NULL) A different cohort, such as another biobank, a later period, or another ancestry group. Build it with prs_cohort(data, reference = ladder$cohort), or with prs_cohort(data, reference = ladder) for an export. Its scores, covariates, and index must be those of the development cohort. If NULL, the function evaluates the development persons. The persons must be separate from the development persons and from the GWAS data that gave the scores (Choi et al. 2020). For a ladder, the function gives a warning with the number of persons of newdata whose ids are also in ladder$cohort. Give newdata or validation, not both.
validation	(prs_bootstrap_plan or NULL; default NULL) A plan from bootstrap_plan() for the optimism correction of the estimates on the development persons (see the section "Optimism correction"). If NULL, the function gives the apparent or the external estimates.
min_events	(integer(1); default 10L; range ⁠[1, Inf)⁠) Minimum number of T2 events at or before a horizon.
min_at_risk	(integer(1); default 50L; range ⁠[1, Inf)⁠) Minimum number of persons whose time since entry is at least the horizon.
min_censoring_survival	(numeric(1); default 0.02; range ⁠(0, 1]⁠) Minimum censoring survival $G(h)$ at a horizon $h$. A value of 1 requires no censoring at or before the horizon.

Details

Persons and outcome

At horizon $h$, the function evaluates only the persons with a prediction (see predict_risk()). The column n_removed counts the other persons. The horizon guard of horizon_support() applies to the evaluated persons. The function leaves out each horizon that fails the guard and gives one warning for all of them.

The outcome is on time since entry, stop - start. At horizon $h$, each person is in one of three classes:

• A case has an event at or before $h$.

• A control is event-free at $h$.

• A person censored at or before $h$ has weight zero.

Cases and controls get inverse probability of censoring weights from the Kaplan-Meier estimate of the censoring distribution, as in assess_thresholds(). Under the recruitment index, incident persons enter at T1 and prevalent persons enter at recruitment (see prs_cohort()). Because of this, the two groups of t1_status can have different censoring on time since entry. When the evaluated persons hold both groups, the function estimates one censoring distribution for each group. The element censoring_strata of settings records whether the evaluated cohort holds both groups. With prevalent = "drop" and under the diagnosis index, the function estimates one censoring distribution for all persons.

The weights assume that censoring is independent of the scores and the covariates. Under the recruitment index, they also assume that censoring is independent of the entry within each group of t1_status.

With competing events (see prs_cohort()), a case has T2 at or before $h$. The controls are the persons free of both events at $h$ and the persons with a competing event at or before $h$ (Blanche et al. 2013). A person with a competing event gets the weight of an event at that time, and the outcome $Y_i$ is 0.

The function computes all measures from the risks of predict_risk().

Measures

• auc is the weighted probability that a case has a higher predicted risk than a control (Blanche et al. 2013). Tied risks count one half. A value of 0.5 means no discrimination and a value of 1 means perfect discrimination.

• brier is the mean of $w_i (Y_i - r_i)^2$ over all evaluated persons (Gerds and Schumacher 2006). Here $w_i$ is the censoring weight, $Y_i$ is 1 for a case and 0 for a control, and $r_i$ is the predicted risk. The mean divides by the number of evaluated persons, which includes the censored persons with weight zero. A lower value is better. The Brier score changes with the proportion of cases, which increases with $h$. Compare Brier scores only at the same horizon.

Differences between models

Each difference is the measure of the full model minus the measure of the reduced model:

• A positive delta_auc favors the full model.

• A negative delta_brier favors the full model.

The function gives no p-value for a difference. For a binary outcome, a test of a difference in AUC between nested models is conservative on the development persons (Demler et al. 2012; Pepe et al. 2013). The likelihood-ratio test of compare_ladder() is the test of added information.

Standard errors and intervals

Each standard error is the standard deviation of an influence function over the persons, divided by $\sqrt{n}$. The influence function of an estimate is $n$ times the derivative of the estimate with respect to the weight of one person. This derivative includes the change in the Kaplan-Meier censoring weights. With one censoring distribution for each group of t1_status, the weight of a person changes only the censoring weights of its own group. For the censoring distribution, the variance of this influence function is the Greenwood variance. The Kaplan-Meier standard errors of cumulative_risk() and assess_calibration() also use the Greenwood variance. The influence function of a difference is the difference of the influence functions of the two models on the same persons. Because of this, its standard error includes the correlation between the two models. Each 95% confidence interval is the estimate plus or minus 1.96 standard errors. The function truncates the intervals of auc and brier to $[0, 1]$. A standard error of 0 gives NA limits, because no person of the sample can then change the estimate. The intervals of delta_auc and delta_brier are not truncated.

The standard errors treat the predicted risks as constant. This assumption holds for persons outside the development sample. On the development persons, the estimates are apparent and tend to be optimistic, and the intervals ignore the uncertainty of the fitted models.

Optimism correction

With a plan from bootstrap_plan() in validation, the function corrects each estimate for optimism with the bootstrap (Harrell Jr. et al. 1996). bootstrap_plan() describes the method, the random number streams, and the failed replicates. The column evaluation is then "optimism-corrected", and each estimate column holds the corrected estimate, the apparent estimate minus the mean optimism. The function truncates the corrected auc and brier to $[0, 1]$. Because the true value is in this range, the truncation never moves an estimate away from it. The element optimism holds the apparent estimate, the optimism, and its Monte Carlo standard error of each estimate. The columns n and n_removed keep their values.

The standard error and the interval depend on interval of the plan:

• "analytic": auc_se and brier_se are the influence-function standard errors of the apparent estimates. The interval is the corrected estimate plus or minus 1.96 of them, truncated to $[0, 1]$. The differences get NA standard errors and limits.

• "two-stage": the standard error is the standard deviation of the outer corrected estimates, and the limits are their 2.5% and 97.5% quantiles.

• "none": the standard errors and limits are NA.

For a model with an AUC near 0.5, the corrected AUC can be too high. Its analytic interval then covers less often than 95% (see the section "Status of the intervals" of validation).

The analytic interval of a difference ignores the variance of the optimism and of the fitted models. In a simulation with an added score of no effect, it covered the true difference in only 73% of the samples. Because of this, the function gives no analytic interval for a difference. The likelihood-ratio test of compare_ladder() tests whether an added score improves a model. No improvement in performance is the same null hypothesis as a coefficient of zero for the added score (Pepe et al. 2013).

Value

A list with the four elements below. Because each element has the same columns in every mode, the results of two modes stack with base::rbind(). The rows are in the order of horizon, then model as in ladder$models or pair as in pairs. When no horizon passes the guard, the data frames have zero rows.

models

One row per horizon and model, with these columns:

model

Model name.

horizon

Horizon $h$.

evaluation

"apparent", "optimism-corrected" with validation, or "external" with newdata.

n

Number of persons evaluated.

n_removed

Number of persons without a prediction.

auc, auc_se, auc_lower, auc_upper

Time-dependent AUC, its standard error, and its 95% confidence interval.

brier, brier_se, brier_lower, brier_upper

Brier score, its standard error, and its 95% confidence interval.

pairs

One row per horizon and pair, with these columns:

full, reduced

Model names of the pair.

horizon, evaluation, n, n_removed

As in models.

delta_auc, delta_auc_se, delta_auc_lower, delta_auc_upper

Difference in AUC, full minus reduced, its standard error, and its 95% confidence interval.

delta_brier, delta_brier_se, delta_brier_lower, delta_brier_upper

Difference in Brier score, full minus reduced, its standard error, and its 95% confidence interval.

optimism

One row per corrected estimate, in the order of horizon, then the measures auc, brier, delta_auc, and delta_brier, then model or pair. It has zero rows without validation.

model

Model name, NA for a pair.

full, reduced

Model names of a pair, NA for a model.

horizon

Horizon $h$.

measure

Name of the measure.

apparent

Apparent estimate.

optimism

Mean optimism over the replicates used. The corrected estimate is apparent - optimism, truncated to the range of the measure. NA when fewer than half of the replicates remain.

optimism_mcse

Monte Carlo standard error of optimism, the standard deviation of the optimism of the replicates divided by the square root of their number.

n_boot_used

Integer number of replicates used.

settings

A list with these elements:

evaluation

The label of the mode.

design

The argument design of fit_ladder().

horizons

The horizons that pass the guard.

min_events, min_at_risk, min_censoring_survival

The guard arguments.

censoring_strata

"t1_status" when the evaluated cohort holds both values of t1_status, else NULL. A horizon or a bootstrap resample whose persons hold one value uses one censoring distribution.

validation

The plan, or NULL without validation.

replicates

A data frame with one row per horizon and the columns horizon, n_boot_used, failed_refit, failed_spline, failed_guard, and failed_risk. The integer columns count the replicates used and the replicates that fail from each cause (see bootstrap_plan()). It has zero rows without validation.

References

Blanche P, Dartigues J, Jacqmin-Gadda H (2013). “Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks.” Statistics in Medicine, 32(30), 5381–5397. doi:10.1002/sim.5958.

Choi SW, Mak TS, O'Reilly PF (2020). “Tutorial: a guide to performing polygenic risk score analyses.” Nature Protocols, 15(9), 2759–2772. doi:10.1038/s41596-020-0353-1.

Demler OV, Pencina MJ, D'Agostino Sr. RB (2012). “Misuse of DeLong test to compare AUCs for nested models.” Statistics in Medicine, 31(23), 2577–2587. doi:10.1002/sim.5328.

Gerds TA, Schumacher M (2006). “Consistent estimation of the expected Brier score in general survival models with right-censored event times.” Biometrical Journal, 48(6), 1029–1040. doi:10.1002/bimj.200610301.

Harrell Jr. FE, Lee KL, Mark DB (1996). “Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.” Statistics in Medicine, 15(4), 361–387. doi:10.1002/(SICI)1097-0258(19960229)15:4<361::AID-SIM168>3.0.CO;2-4.

Pepe MS, Kerr KF, Longton G, Wang Z (2013). “Testing for improvement in prediction model performance.” Statistics in Medicine, 32(9), 1467–1482. doi:10.1002/sim.5727.

See Also

horizon_support() to check horizons before the evaluation, validation for the three modes, bootstrap_plan() for the optimism correction, compare_ladder() for the likelihood-ratio test of added scores, assess_calibration(), and assess_thresholds().

Examples

#----------------------------------------------------------------------------
# assess_discrimination() example
#----------------------------------------------------------------------------
# Fit a ladder on a development cohort.
dat <- simulate_prs_data(n = 2000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(onset = "onset_score", progression = "progression_score")
)
ladder <- fit_ladder(cohort)

# Build an external cohort with the rules of the development cohort.
external <- prs_cohort(
  simulate_prs_data(n = 2000L, seed = 2L, id_prefix = "E"),
  reference = cohort
)

# Assess discrimination at 5 and 10 years, and the gain from adding the
# progression score to the onset score.
discrimination <- assess_discrimination(
  ladder,
  horizons = c(5, 10),
  pairs = list(c("onset+progression", "onset")),
  newdata = external
)
discrimination$models
discrimination$pairs

#----------------------------------------------------------------------------
# Optimism correction on the development persons
#----------------------------------------------------------------------------
# A small number of replicates keeps the example fast. An analysis needs
# 200 or more.
corrected <- assess_discrimination(
  ladder,
  horizons = 10,
  validation = bootstrap_plan(n_boot = 20L, seed = 1L)
)
corrected$models
corrected$optimism

---

Assess Net Benefit for Decision Curves

Description

assess_net_benefit() gives the net benefit of each model of a ladder at risk thresholds, at each horizon. The net benefit weighs the true positives of a rule against its false positives. The weight of a false positive is the odds of the threshold, $p_t / (1 - p_t)$. This weight is the harm of one false positive in units of the benefit of one true positive. Its reciprocal is the largest number of false positives that a clinician accepts for one true positive. A plot of net benefit against the threshold is a decision curve (Vickers and Elkin 2006).

The function also gives the net benefit of treating all persons, the standardized net benefit, and the change in net benefit between pairs of models. plot_net_benefit() draws the decision curves.

The function has three modes. Without newdata and validation, it evaluates the development persons, and the estimates are apparent. With validation, it corrects the apparent estimates for optimism with the bootstrap. With newdata, it evaluates other persons. validation explains when to use each mode.

Usage

assess_net_benefit(
  ladder,
  horizons,
  thresholds,
  pairs = NULL,
  newdata = NULL,
  validation = NULL,
  min_events = 10L,
  min_at_risk = 50L,
  min_censoring_survival = 0.02
)

Arguments

ladder	(prs_ladder or prs_export) A ladder from fit_ladder() or its export from export_ladder(). An export holds no persons, and it needs newdata.
horizons	(numeric(); range ⁠(0, Inf)⁠) Distinct finite times since entry, in the units of the age columns.
thresholds	(numeric(); range ⁠(0, 1)⁠) Strictly increasing risk thresholds $p_t$. The function applies the same thresholds at every horizon. The argument has no default. See the section "Thresholds" in Details.
pairs	(list of character(2) or NULL; default NULL) Pairs of models to compare, each as c(full, reduced) model names. The names must be in names(ladder$models) and must differ. Pairs need not be nested, and a pair must not occur twice. If NULL, the function compares every nested pair in the order of compare_ladder().
newdata	(prs_cohort or NULL; default NULL) A different cohort, such as another biobank, a later period, or another ancestry group. Build it with prs_cohort(data, reference = ladder$cohort), or with prs_cohort(data, reference = ladder) for an export. Its scores, covariates, and index must be those of the development cohort. If NULL, the function evaluates the development persons. The persons must be separate from the development persons and from the GWAS data that gave the scores (Choi et al. 2020). For a ladder, the function gives a warning with the number of persons of newdata whose ids are also in ladder$cohort. Give newdata or validation, not both.
validation	(prs_bootstrap_plan or NULL; default NULL) A plan from bootstrap_plan() for the optimism correction of the estimates on the development persons (see the section "Optimism correction"). If NULL, the function gives the apparent or the external estimates.
min_events	(integer(1); default 10L; range ⁠[1, Inf)⁠) Minimum number of T2 events at or before a horizon.
min_at_risk	(integer(1); default 50L; range ⁠[1, Inf)⁠) Minimum number of persons whose time since entry is at least the horizon.
min_censoring_survival	(numeric(1); default 0.02; range ⁠(0, 1]⁠) Minimum censoring survival $G(h)$ at a horizon $h$. A value of 1 requires no censoring at or before the horizon.

Details

Persons and weights

The function evaluates the persons of assess_thresholds(), with the same horizon guard and the same censoring weights. See the section "Persons and weights" of assess_thresholds(). Under the recruitment index, the weights use one censoring distribution for each group of t1_status, as in assess_discrimination(). With competing events, a person with a competing event by the horizon is a non-event. The prevalence is then the cumulative incidence of T2 by the horizon (Vickers et al. 2008; van Geloven et al. 2022).

The horizon is time since entry, as in the other ⁠assess_*()⁠ functions. Under the recruitment index or with a positive entry_lag in prs_cohort(), this window differs from the time since T1 of cumulative_risk().

Net benefit

Let $w_i$ be the weight of person $i$. Let $e_i$ and $c_i$ be the indicators of an event and a non-event. Let $f_i$ be the indicator of a risk at or above the threshold $p_t$. With $n$ evaluated persons and the prevalence $\rho = n^{-1} \sum_i w_i e_i$:

$\begin{aligned} NB &= \frac{1}{n} \left[\sum_i w_i e_i f_i - \frac{p_t}{1 - p_t} \sum_i w_i c_i f_i\right], \\ NB_\text{all} &= \rho - (1 - \rho) \frac{p_t}{1 - p_t}. \end{aligned}$

Because the weights sum to $n$, the net benefit also equals $\rho\,\text{TPR} - (1 - \rho)\,\text{FPR}\,p_t/(1 - p_t)$. Here, TPR is the sensitivity and FPR is the false_positive_rate of assess_thresholds(). With each weight equal to 1, this net benefit is the net benefit of Vickers and Elkin (2006) for a binary outcome. The weights extend it to censored data. In a simulation, its mean equaled the mean of the Aalen-Johansen form of Vickers et al. (2008) within the Monte Carlo error.

Treating no person has a net benefit of 0 at every threshold. A model is useful at a threshold when its net benefit is larger than that of treating all persons and larger than 0.

The standardized net benefit is $NB / \rho$ (Kerr et al. 2016). Its largest value is 1, for a rule that flags each event and no non-event. Here, "standardized" means the division by $\rho$. It is not the standardized score of prs_cohort(), and not the average over the cohort of standardized_risk().

The change in net benefit of a pair is $NB_\text{full} - NB_\text{reduced}$ at the same threshold. For one threshold, it equals the weighted NRI with two categories (Kerr et al. 2014). The weight is $s_1 = 1$ for events and $s_2 = p_t / (1 - p_t)$ for non-events:

$\Delta NB = \rho\,\text{NRI}_\text{events} + (1 - \rho) \frac{p_t}{1 - p_t}\,\text{NRI}_\text{nonevents},$

with the components nri_events and nri_nonevents of assess_thresholds() at that threshold.

Thresholds

thresholds has no default, because the range of reasonable thresholds depends on the clinical decision (Vickers et al. 2016). A decision curve needs many thresholds, such as seq(0.01, 0.5, 0.01). assess_thresholds() is for a few thresholds, because many thresholds give many categories of the NRI (Kerr et al. 2014).

Standard errors and intervals

Each standard error is the standard deviation of an influence function over the persons, divided by $\sqrt{n}$, as in assess_discrimination(). The influence function includes the change in the Kaplan-Meier censoring weights. The net benefit and the prevalence are means of weighted terms, and their influence functions have the form of that of the Brier score. The standardized net benefit is a ratio of two weighted sums. The net benefit of treating all persons and its standardized value are functions of the prevalence alone. Each 95% interval is the estimate plus or minus 1.96 standard errors, because a net benefit can be negative. When a model flags no person at a threshold, its net benefit and the standard error are 0. No person of the sample can then change the estimate, and the limits are NA. The function does not use the random number generator.

The net benefit on the development persons is optimistic, because fit_ladder() fitted the models on the same persons. For apparent evaluation, the intervals also ignore the uncertainty of the fitted models. Use a plan in validation, or a different cohort in newdata, to get a decision curve without this optimism.

Optimism correction

With a plan from bootstrap_plan() in validation, the function corrects net_benefit, standardized_net_benefit, and delta_net_benefit at each threshold for optimism with the bootstrap (Harrell Jr. et al. 1996; Vickers et al. 2008). bootstrap_plan() describes the method, the random number streams, and the failed replicates. The column evaluation of models and pairs is then "optimism-corrected", and each estimate column holds the corrected estimate, the apparent estimate minus the mean optimism. The element optimism holds the apparent estimate, the optimism, and its Monte Carlo standard error of each estimate. The columns n and n_removed keep their values.

The prevalence and the net benefit of treating all persons do not depend on a model. The element treat_all keeps the values of the development persons, with the label "apparent".

Vickers et al. (2008) recommend repeated 10-fold cross-validation over this correction. In their simulation, the corrected estimates of the two methods differed by less than 0.005 at the thresholds below 60%. The bootstrap had the least bias at 60% and 80%.

The standard error and the interval depend on interval of the plan:

• "analytic": net_benefit_se and standardized_net_benefit_se are the influence-function standard errors of the apparent estimates. The interval is the corrected estimate plus or minus 1.96 of them. When the model flags no person, the apparent standard error is 0, and the limits are NA. The corrected estimate stays, and it can differ from 0, because the refits can flag persons. delta_net_benefit gets an NA standard error and NA limits, as the differences of assess_discrimination() do.

• "two-stage": the standard error is the standard deviation of the outer corrected estimates, and the limits are their 2.5% and 97.5% quantiles.

• "none": the standard errors and limits are NA.

The analytic interval is provisional. A simulation study before the release of the package measures its coverage.

Value

A list with the five elements below. Because each element has the same columns in every mode, the results of two modes stack with base::rbind(). The rows are in the order of horizon, then model as in ladder$models or pair as in pairs, then threshold. When no horizon passes the guard, the data frames have zero rows. Each estimate has three more columns with the suffixes ⁠_se⁠, ⁠_lower⁠, and ⁠_upper⁠, for its standard error and 95% interval (see the section "Standard errors and intervals").

models

One row per horizon, model, and threshold, with these columns:

model

Model name.

horizon

Horizon $h$.

evaluation

"apparent", "optimism-corrected" with validation, or "external" with newdata.

threshold

Risk threshold $p_t$.

n

Number of persons evaluated.

n_removed

Number of persons without a prediction.

net_benefit

Net benefit of the model.

standardized_net_benefit

Net benefit divided by the prevalence.

treat_all

One row per horizon and threshold, with the columns horizon, evaluation, threshold, n, and n_removed of models, and these columns. With validation, evaluation is "apparent".

prevalence

Weighted cumulative incidence of T2 by the horizon, $\rho$. It is the same at every threshold of a horizon.

net_benefit

Net benefit of treating all persons.

standardized_net_benefit

Net benefit of treating all persons divided by the prevalence.

pairs

One row per horizon, pair, and threshold, with these columns:

full, reduced

Model names of the pair.

horizon, evaluation, threshold, n, n_removed

As in models.

delta_net_benefit

Net benefit of full minus net benefit of reduced.

optimism

One row per corrected estimate, with the columns of the element optimism of assess_discrimination() and the column threshold after horizon. The rows are in the order of horizon, then the measures net_benefit, standardized_net_benefit, and delta_net_benefit, then model or pair, then threshold. It has zero rows without validation.

settings

A list with the elements of the element settings of assess_discrimination().

References

Choi SW, Mak TS, O'Reilly PF (2020). “Tutorial: a guide to performing polygenic risk score analyses.” Nature Protocols, 15(9), 2759–2772. doi:10.1038/s41596-020-0353-1.

Harrell Jr. FE, Lee KL, Mark DB (1996). “Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.” Statistics in Medicine, 15(4), 361–387. doi:10.1002/(SICI)1097-0258(19960229)15:4<361::AID-SIM168>3.0.CO;2-4.

Kerr KF, Brown MD, Zhu K, Janes H (2016). “Assessing the clinical impact of risk prediction models with decision curves: guidance for correct interpretation and appropriate use.” Journal of Clinical Oncology, 34(21), 2534–2540. doi:10.1200/JCO.2015.65.5654.

Kerr KF, Wang Z, Janes H, McClelland RL, Psaty BM, Pepe MS (2014). “Net reclassification indices for evaluating risk prediction instruments: a critical review.” Epidemiology, 25(1), 114–121. doi:10.1097/EDE.0000000000000018.

van Geloven N, Giardiello D, Bonneville EF, Teece L, Ramspek CL, van Smeden M, Snell KIE, van Calster B, Pohar-Perme M, Riley RD, Putter H, Steyerberg E (2022). “Validation of prediction models in the presence of competing risks: a guide through modern methods.” BMJ, 377, e069249. doi:10.1136/bmj-2021-069249.

Vickers AJ, Cronin AM, Elkin EB, Gonen M (2008). “Extensions to decision curve analysis, a novel method for evaluating diagnostic tests, prediction models and molecular markers.” BMC Medical Informatics and Decision Making, 8, 53. doi:10.1186/1472-6947-8-53.

Vickers AJ, Elkin EB (2006). “Decision curve analysis: a novel method for evaluating prediction models.” Medical Decision Making, 26(6), 565–574. doi:10.1177/0272989X06295361.

Vickers AJ, Van Calster B, Steyerberg EW (2016). “Net benefit approaches to the evaluation of prediction models, molecular markers, and diagnostic tests.” BMJ, 352, i6. doi:10.1136/bmj.i6.

See Also

validation for the three modes, plot_net_benefit() for the decision curves, and assess_thresholds() for the sensitivity, the predictive values, and the NRI at a few thresholds.

Examples

#----------------------------------------------------------------------------
# assess_net_benefit() example
#----------------------------------------------------------------------------
library(prsladder)

# Fit a ladder on a development cohort.
dat <- simulate_prs_data(n = 2000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(onset = "onset_score", progression = "progression_score")
)
ladder <- fit_ladder(cohort)

# Build an external cohort with the rules of the development cohort.
external <- prs_cohort(
  simulate_prs_data(n = 2000L, seed = 2L, id_prefix = "E"),
  reference = cohort
)

# Net benefit of the 10-year risk at thresholds from 5% to 40%.
net_benefit <- assess_net_benefit(
  ladder,
  horizons = 10,
  thresholds = seq(0.05, 0.4, by = 0.05),
  pairs = list(c("onset+progression", "onset")),
  newdata = external
)
net_benefit$models[
  net_benefit$models$model == "onset+progression",
  c("threshold", "net_benefit", "standardized_net_benefit")
]
net_benefit$treat_all[c("threshold", "prevalence", "net_benefit")]
net_benefit$pairs[c("threshold", "delta_net_benefit")]

---

Assess the Stability of the Predicted Risks

Description

assess_stability() measures how much the predicted risk of each person changes when fit_ladder() fits the ladder on other samples of the same size. It uses the bootstrap samples of a plan from bootstrap_plan(). validation explains the internal validation that gives these samples. For each model and horizon, it reports the mean and the largest mean absolute prediction error over the persons (Riley and Collins 2023).

Usage

assess_stability(
  ladder,
  horizons,
  validation,
  min_events = 10L,
  min_at_risk = 50L,
  min_censoring_survival = 0.02
)

Arguments

ladder	(prs_ladder) A ladder from fit_ladder().
horizons	(numeric(); range ⁠(0, Inf)⁠) Distinct finite times since entry, in the units of the age columns.
validation	(prs_bootstrap_plan) A plan from bootstrap_plan(). The argument has no default.
min_events	(integer(1); default 10L; range ⁠[1, Inf)⁠) Minimum number of T2 events at or before a horizon.
min_at_risk	(integer(1); default 50L; range ⁠[1, Inf)⁠) Minimum number of persons whose time since entry is at least the horizon.
min_censoring_survival	(numeric(1); default 0.02; range ⁠(0, 1]⁠) Minimum censoring survival $G(h)$ at a horizon $h$. A value of 1 requires no censoring at or before the horizon.

Details

Mean absolute prediction error

Let $r_i$ be the risk of person $i$ at horizon $h$ from the model fitted on the development persons $D$. Let $r_{ib}$ be the risk of the same person from the model fitted again on the bootstrap sample $D_b$ of replicate $b$. The mean absolute prediction error of person $i$ is

$\text{MAPE}_i = \frac{1}{B} \sum_{b=1}^{B} |r_{ib} - r_i|,$

over the $B$ replicates that do not fail at $h$ (Riley and Collins 2023). The function reports the mean and the maximum of $\text{MAPE}_i$ over the persons that the other ⁠assess_*()⁠ functions evaluate at $h$. A large value shows that the risk of a person depends on the sample that gave the model. Such instability is larger in a small cohort, and it is a sign of overfitting.

With competing events, the risk is the cumulative incidence of T2 from the models of both events, as in predict_risk().

Replicates

The function uses the replicates and the failure rules of the plan (see bootstrap_plan()). With the same plan and the same guard arguments, the replicates at a horizon are the replicates of the corrected mode of the other ⁠assess_*()⁠ functions. The function ignores interval and n_inner of the plan. When fewer than half of the n_boot replicates remain at a horizon, mape_mean and mape_max are NA.

The horizon guard of horizon_support() applies to the persons with a prediction, as in assess_discrimination(). The function leaves out each horizon that fails the guard and gives one warning for all of them.

Value

A list with the two elements below.

models

A data frame with one row per horizon and model, in the order of horizon, then model as in ladder$models. It has zero rows when no horizon passes the guard.

model

Model name.

horizon

Horizon $h$.

n

Integer number of persons evaluated.

n_boot_used

Integer number of replicates used.

mape_mean

Mean of the mean absolute prediction error over the persons.

mape_max

Maximum of the mean absolute prediction error over the persons.

settings

A list with the elements design (the argument design of fit_ladder()), horizons (the horizons that pass the guard), min_events, min_at_risk, min_censoring_survival, censoring_strata (as in assess_discrimination()), validation (the plan), and replicates. replicates is a data frame with one row per horizon and the integer columns n_boot_used, failed_refit, failed_spline, failed_guard, and failed_risk after the column horizon. These columns count the replicates used and the replicates that fail from each cause (see bootstrap_plan()).

References

Riley RD, Collins GS (2023). “Stability of clinical prediction models developed using statistical or machine learning methods.” Biometrical Journal, 65(8), e2200302. doi:10.1002/bimj.202200302.

See Also

validation for the internal validation, bootstrap_plan() for the replicates, and assess_discrimination() and assess_calibration() for the corrected measures.

Examples

#----------------------------------------------------------------------------
# assess_stability() example
#----------------------------------------------------------------------------
library(prsladder)

dat <- simulate_prs_data(n = 1000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(onset = "onset_score", progression = "progression_score")
)
ladder <- fit_ladder(cohort)

# A small number of replicates keeps the example fast. An analysis needs
# 200 or more.
plan <- bootstrap_plan(n_boot = 20L, seed = 1L)
assess_stability(ladder, horizons = c(5, 10), validation = plan)$models

---

Assess Risk Thresholds and Reclassification

Description

The function reports how each model of a ladder classifies persons at fixed risk thresholds, at each horizon. For each model and threshold, it gives the proportion of persons at or above the threshold, the sensitivity, the false-positive rate, and the positive and negative predictive values. For pairs of models, it gives the categorical net reclassification improvement (NRI) over the risk categories that the thresholds define. The PRS reporting standards ask for these proportions and for the NRI at clinically relevant thresholds (Wand et al. 2021).

The function has three modes. Without newdata and validation, it evaluates the development persons, and the estimates are apparent. With validation, it corrects the apparent estimates for optimism with the bootstrap. With newdata, it evaluates other persons. validation explains when to use each mode.

Usage

assess_thresholds(
  ladder,
  horizons,
  thresholds,
  pairs = NULL,
  newdata = NULL,
  validation = NULL,
  min_events = 10L,
  min_at_risk = 50L,
  min_censoring_survival = 0.02
)

Arguments

ladder	(prs_ladder or prs_export) A ladder from fit_ladder() or its export from export_ladder(). An export holds no persons, and it needs newdata.
horizons	(numeric(); range ⁠(0, Inf)⁠) Distinct finite times since entry, in the units of the age columns.
thresholds	(numeric(); range ⁠(0, 1)⁠) Strictly increasing risk thresholds. The function applies the same thresholds at every horizon. Because predicted risk increases with the horizon, the same threshold flags more persons at a longer horizon. To use different thresholds at different horizons, call the function once per horizon.
pairs	(list of character(2) or NULL; default NULL) Pairs of models to compare, each as c(full, reduced) model names. The names must be in names(ladder$models) and must differ. Pairs need not be nested, and a pair must not occur twice. If NULL, the function compares every nested pair in the order of compare_ladder().
newdata	(prs_cohort or NULL; default NULL) A different cohort, such as another biobank, a later period, or another ancestry group. Build it with prs_cohort(data, reference = ladder$cohort), or with prs_cohort(data, reference = ladder) for an export. Its scores, covariates, and index must be those of the development cohort. If NULL, the function evaluates the development persons. The persons must be separate from the development persons and from the GWAS data that gave the scores (Choi et al. 2020). For a ladder, the function gives a warning with the number of persons of newdata whose ids are also in ladder$cohort. Give newdata or validation, not both.
validation	(prs_bootstrap_plan or NULL; default NULL) A plan from bootstrap_plan() for the optimism correction of the estimates on the development persons (see the section "Optimism correction"). If NULL, the function gives the apparent or the external estimates.
min_events	(integer(1); default 10L; range ⁠[1, Inf)⁠) Minimum number of T2 events at or before a horizon.
min_at_risk	(integer(1); default 50L; range ⁠[1, Inf)⁠) Minimum number of persons whose time since entry is at least the horizon.
min_censoring_survival	(numeric(1); default 0.02; range ⁠(0, 1]⁠) Minimum censoring survival $G(h)$ at a horizon $h$. A value of 1 requires no censoring at or before the horizon.

Details

Persons and weights

At horizon $h$, the function evaluates only the persons with a prediction (see predict_risk()). The column n_removed counts the other persons. The horizon guard of horizon_support() applies to the evaluated persons. The function leaves out each horizon that fails the guard and gives one warning for all of them.

Let $T_i$ be the time since entry of person $i$, and $G$ the Kaplan-Meier estimate of the censoring distribution on the evaluated persons. Each person is in one of three classes:

• An event has T2 at or before $h$, with weight $1 / G(T_i-)$.

• A non-event has $T_i > h$, with weight $1 / G(h)$. With competing events, a person with a competing event at or before $h$ is also a non-event, with weight $1 / G(T_i-)$.

• A person censored at or before $h$ has weight zero.

A competing event is not censoring in the estimate of $G$. Under the recruitment index, the evaluated persons can hold both groups of t1_status. $G$ is then the estimate on the evaluated persons of the group of person $i$. assess_discrimination() describes these groups.

These weights are the inverse probability of censoring weights of Blanche et al. (2013), with the left limit $G(T_i-)$ of Gerds and Schumacher (2006). With one censoring distribution for all persons, they are the weights of riskRegression::Score() with cens.model = "km". That model of riskRegression::Score() uses one censoring distribution also with strat() in the formula. The weights assume that censoring is independent of the scores and the covariates. Under the recruitment index, they also assume that censoring is independent of the entry within each group of t1_status.

Measures at a threshold

A risk equal to a threshold counts as at or above it:

• above is the unweighted proportion of all evaluated persons at or above the threshold. It includes the censored persons.

• sensitivity is the weighted proportion of events at or above the threshold.

• false_positive_rate is the weighted proportion of non-events at or above the threshold.

• ppv is the positive predictive value, the weighted proportion of events among the persons at or above the threshold.

• npv is the negative predictive value, the weighted proportion of non-events among the persons below the threshold.

Let $w_i$ be the weight of person $i$. Let $e_i$ and $c_i$ be the indicators of an event and a non-event. Let $f_i$ be the indicator of a risk at or above the threshold. Then:

$\text{PPV} = \frac{\sum_i w_i e_i f_i}{\sum_i w_i f_i}, \qquad \text{NPV} = \frac{\sum_i w_i c_i (1 - f_i)}{\sum_i w_i (1 - f_i)}.$

Zheng et al. (2008) define the time-dependent PPV and NPV. They estimate them from a Nelson-Aalen, kernel, or Cox estimate of survival, not with censoring weights. This function estimates them with the weights of Blanche et al. (2013), which those authors give for the sensitivity and the specificity. With competing events, a person with a competing event by the horizon is a non-event (Blanche et al. 2013; van Geloven et al. 2022). ppv is then the cumulative incidence of T2 among the persons at or above the threshold (van Geloven et al. 2022). npv is the probability of no T2 by the horizon among the persons below the threshold.

The weights sum to the number of evaluated persons $n$. Let $\rho = n^{-1} \sum_i w_i e_i$ be the prevalence. Then the predictive values follow from the other measures:

$\text{PPV} = \frac{\rho\,\text{TPR}} {\rho\,\text{TPR} + (1 - \rho)\,\text{FPR}},$

where TPR is sensitivity and FPR is false_positive_rate. Because above has no weights, sensitivity * prevalence / above is close to ppv, but not equal.

ppv is NA when the persons at or above the threshold have no weight, for example when the threshold is above every risk. npv is NA when the persons below the threshold have no weight, for example when the threshold is at or below every risk.

The horizon is time since entry, as in the other ⁠assess_*()⁠ functions. Under the recruitment index or with a positive entry_lag in prs_cohort(), a person can enter after T1. A predictive value at a horizon of 10 then covers the 10 time units after entry. The risk at time 10 of cumulative_risk() covers the 10 time units after T1.

Reclassification

For a pair, the thresholds divide the risk into categories. For each person, the function compares the category of the full model with the category of the reduced model. The move is up when the category of the full model is higher, and down when it is lower. The categorical NRI of Pencina et al. (2008) has these terms:

• nri_events is the weighted proportion of events that move up minus the weighted proportion that move down.

• nri_nonevents is the weighted proportion of non-events that move down minus the weighted proportion that move up.

Pencina et al. (2008) define the NRI for a binary outcome without censoring. The censoring weights extend it to censored data and to competing events.

The function does not report the sum of the two components. The sum is not a proportion, and each component is a difference of proportions (Kerr et al. 2014). With one threshold, nri_events is the change in sensitivity. Also, nri_nonevents is the change in false-positive rate with the sign reversed. With three or more categories, the NRI does not account for the clinical importance of each move (Kerr et al. 2014).

The NRI can be positive for an added score with no effect when a model is overfit or poorly calibrated. This bias also occurs on persons outside the development sample. Hilden and Gerds (2014) and Pepe et al. (2015) show this bias for the continuous NRI. Pepe et al. (2015) report the same bias for the categorical NRI. Check calibration with assess_calibration(). On persons outside the development sample, the AUC and the Brier score of assess_discrimination() do not improve because of a poorly fitting model (Hilden and Gerds 2014; Pepe et al. 2015). The function does not report the continuous NRI or the integrated discrimination improvement, which have the same bias (Hilden and Gerds 2014).

Standard errors and intervals

Each standard error is the standard deviation of an influence function over the persons, divided by $\sqrt{n}$, as in assess_discrimination(). The influence function includes the change in the Kaplan-Meier censoring weights. Each measure except above is a ratio of two weighted sums, and its influence function has the form of that of the AUC. above has no weights, and its standard error is that of a proportion.

The proportions are above, sensitivity, false_positive_rate, ppv, and npv. The 95% interval of a proportion is the Wald interval on the logit scale:

$\mathrm{expit}\left\{\mathrm{logit}(\hat p) \pm 1.96 \frac{\widehat{\mathrm{SE}}}{\hat p (1 - \hat p)}\right\}.$

Its limits are NA when $\hat p$ is 0 or 1. The NRI terms can be negative, and their intervals are the Wald intervals on the raw scale, the estimate plus or minus 1.96 standard errors. The standard error and the limits are NA when the estimate is NA. The limits are also NA when the standard error is 0, such as for an NRI term when no person moves. No person of the sample can then change the estimate. The function does not use the random number generator.

The standard errors treat the predicted risks as constant. This assumption holds for persons outside the development sample. On the development persons, the estimates are apparent and tend to be optimistic, and the intervals ignore the uncertainty of the fitted models.

Optimism correction

With a plan from bootstrap_plan() in validation, the function corrects each measure at each threshold for optimism with the bootstrap (Harrell Jr. et al. 1996). bootstrap_plan() describes the method, the random number streams, and the failed replicates. The column evaluation is then "optimism-corrected", and each estimate column holds the corrected estimate, the apparent estimate minus the mean optimism. The element optimism holds the apparent estimate, the optimism, and its Monte Carlo standard error of each estimate. The columns n and n_removed keep their values.

The function corrects a proportion on the raw scale. A logit scale is not possible, because a refit that flags no person gives a sensitivity of 0, with an infinite logit. The raw difference can be below 0 or above 1 when the apparent value is 0, 1, or near these values. This occurs at a threshold near the edge of the risks of a model. The function then truncates the corrected proportion to 0 or 1, and its analytic limits are NA. Because the true value is in $[0, 1]$, the truncation never moves an estimate away from it. The element optimism keeps the raw difference.

A measure that is NA in a replicate, such as ppv when the refit flags no person, does not use that replicate. The column n_boot_used of optimism counts the replicates used for each row. When fewer than half of the replicates remain, the corrected estimate is NA.

The standard error and the interval depend on interval of the plan:

• "analytic": the standard error of a proportion is the influence-function standard error of the apparent estimate. The interval is the logit-scale interval of the apparent estimate. The function moves this interval to the corrected estimate:

  $\mathrm{expit}\left\{\mathrm{logit}(\hat p_\text{corr}) \pm 1.96 \frac{\widehat{\mathrm{SE}}} {\hat p_\text{app} (1 - \hat p_\text{app})}\right\}.$

  Its limits are NA when $\hat p_\text{app}$ or $\hat p_\text{corr}$ is not in $(0, 1)$. The NRI components get NA standard errors and limits. Such an interval of a difference between models ignores the variance of the optimism and of the fitted models (see assess_discrimination()).

• "two-stage": the standard error is the standard deviation of the outer corrected estimates, and the limits are their 2.5% and 97.5% quantiles. The function truncates the limits of a proportion to $[0, 1]$.

• "none": the standard errors and limits are NA.

The analytic interval is provisional. A simulation study before the release of the package measures its coverage.

Value

A list with the four elements below. Because each element has the same columns in every mode, the results of two modes stack with base::rbind(). The rows are in the order of horizon, then model as in ladder$models and threshold, or pair as in pairs. When no horizon passes the guard, the data frames have zero rows. Each measure has three more columns with the suffixes ⁠_se⁠, ⁠_lower⁠, and ⁠_upper⁠, for its standard error and 95% interval (see the section "Standard errors and intervals").

models

One row per horizon, model, and threshold, with these columns:

model

Model name.

horizon

Horizon $h$.

evaluation

"apparent", "optimism-corrected" with validation, or "external" with newdata.

threshold

Risk threshold.

n

Number of persons evaluated.

n_removed

Number of persons without a prediction.

above

Proportion of persons at or above the threshold.

sensitivity

Weighted proportion of events at or above the threshold.

false_positive_rate

Weighted proportion of non-events at or above the threshold.

ppv

Positive predictive value: weighted proportion of events among the persons at or above the threshold.

npv

Negative predictive value: weighted proportion of non-events among the persons below the threshold.

pairs

One row per horizon and pair, with these columns:

full, reduced

Model names of the pair.

horizon, evaluation, n, n_removed

As in models.

nri_events

NRI component for events.

nri_nonevents

NRI component for non-events.

optimism

One row per corrected estimate, with the columns of the element optimism of assess_discrimination() and the column threshold after horizon. The rows are in the order of horizon, then the measures above, sensitivity, false_positive_rate, ppv, npv, nri_events, and nri_nonevents, then model and threshold, or pair. threshold is NA for an NRI component. It has zero rows without validation.

settings

A list with the elements of the element settings of assess_discrimination().

References

Blanche P, Dartigues J, Jacqmin-Gadda H (2013). “Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks.” Statistics in Medicine, 32(30), 5381–5397. doi:10.1002/sim.5958.

Choi SW, Mak TS, O'Reilly PF (2020). “Tutorial: a guide to performing polygenic risk score analyses.” Nature Protocols, 15(9), 2759–2772. doi:10.1038/s41596-020-0353-1.

Gerds TA, Schumacher M (2006). “Consistent estimation of the expected Brier score in general survival models with right-censored event times.” Biometrical Journal, 48(6), 1029–1040. doi:10.1002/bimj.200610301.

Harrell Jr. FE, Lee KL, Mark DB (1996). “Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.” Statistics in Medicine, 15(4), 361–387. doi:10.1002/(SICI)1097-0258(19960229)15:4<361::AID-SIM168>3.0.CO;2-4.

Hilden J, Gerds TA (2014). “A note on the evaluation of novel biomarkers: do not rely on integrated discrimination improvement and net reclassification index.” Statistics in Medicine, 33(19), 3405–3414. doi:10.1002/sim.5804.

Kerr KF, Wang Z, Janes H, McClelland RL, Psaty BM, Pepe MS (2014). “Net reclassification indices for evaluating risk prediction instruments: a critical review.” Epidemiology, 25(1), 114–121. doi:10.1097/EDE.0000000000000018.

Pencina MJ, D'Agostino Sr. RB, D'Agostino Jr. RB, Vasan RS (2008). “Evaluating the added predictive ability of a new marker: from area under the ROC curve to reclassification and beyond.” Statistics in Medicine, 27(2), 157–172. doi:10.1002/sim.2929.

Pepe MS, Fan J, Feng Z, Gerds T, Hilden J (2015). “The net reclassification index (NRI): a misleading measure of prediction improvement even with independent test data sets.” Statistics in Biosciences, 7(2), 282–295. doi:10.1007/s12561-014-9118-0.

van Geloven N, Giardiello D, Bonneville EF, Teece L, Ramspek CL, van Smeden M, Snell KIE, van Calster B, Pohar-Perme M, Riley RD, Putter H, Steyerberg E (2022). “Validation of prediction models in the presence of competing risks: a guide through modern methods.” BMJ, 377, e069249. doi:10.1136/bmj-2021-069249.

Wand H, Lambert SA, Tamburro C, Iacocca MA, O'Sullivan JW, Sillari C, Kullo IJ, Rowley R, Dron JS, Brockman D, Venner E, McCarthy MI, Antoniou AC, Easton DF, Hegele RA, Khera AV, Chatterjee N, Kooperberg C, Edwards K, Vlessis K, Kinnear K, Danesh JN, Parkinson H, Ramos EM, Roberts MC, Ormond KE, Khoury MJ, Janssens ACJW, Goddard KAB, Kraft P, MacArthur JAL, Inouye M, Wojcik GL (2021). “Improving reporting standards for polygenic scores in risk prediction studies.” Nature, 591(7849), 211–219. doi:10.1038/s41586-021-03243-6.

Zheng Y, Cai T, Pepe MS, Levy WC (2008). “Time-dependent predictive values of prognostic biomarkers with failure time outcome.” Journal of the American Statistical Association, 103(481), 362–368. doi:10.1198/016214507000001481.

See Also

validation for the three modes, assess_calibration() to check calibration, assess_discrimination() for the AUC and the Brier score, and compare_ladder() for the likelihood-ratio test of added scores.

Examples

#----------------------------------------------------------------------------
# assess_thresholds() example
#----------------------------------------------------------------------------
# Fit a ladder on a development cohort.
dat <- simulate_prs_data(n = 2000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(onset = "onset_score", progression = "progression_score")
)
ladder <- fit_ladder(cohort)

# Build an external cohort with the rules of the development cohort.
external <- prs_cohort(
  simulate_prs_data(n = 2000L, seed = 2L, id_prefix = "E"),
  reference = cohort
)

# Classify at 10-year risks of 20% and 40%, and compare the model with
# both scores to the model with the onset score only.
thresholds <- assess_thresholds(
  ladder,
  horizons = 10,
  thresholds = c(0.2, 0.4),
  pairs = list(c("onset+progression", "onset")),
  newdata = external
)
thresholds$models
thresholds$pairs

#----------------------------------------------------------------------------
# Optimism correction on the development persons
#----------------------------------------------------------------------------
# A small number of replicates keeps the example fast. An analysis needs
# 200 or more.
corrected <- assess_thresholds(
  ladder,
  horizons = 10,
  thresholds = 0.2,
  pairs = list(c("onset+progression", "onset")),
  validation = bootstrap_plan(n_boot = 20L, seed = 1L)
)
corrected$models[c("model", "sensitivity", "ppv", "ppv_lower", "ppv_upper")]

---

Plan the Bootstrap of an Internal Validation

Description

bootstrap_plan() holds the settings of the bootstrap optimism correction of the ⁠assess_*()⁠ functions. Give the plan to the argument validation of assess_discrimination(), assess_calibration(), assess_thresholds(), assess_net_benefit(), or assess_stability(). validation explains when to use the internal validation. The plan holds no data. The same plan gives the same bootstrap samples in every function.

Usage

bootstrap_plan(
  n_boot = 200L,
  seed = NULL,
  interval = c("analytic", "two-stage", "none"),
  n_inner = 50L
)

Arguments

n_boot	(integer(1); default 200L; range ⁠[2, Inf)⁠) Number of bootstrap replicates $B$. With interval = "two-stage", the number of outer replicates. A double with no fractional part is also accepted.
seed	(integer(1) or NULL; default NULL) Seed of the random number streams. If NULL, the function draws one seed from the current generator (see the section "Random numbers").
interval	(character(1); default "analytic") The interval of the corrected estimates: "analytic", "two-stage", or "none" (see the section "Intervals").
n_inner	(integer(1); default 50L; range ⁠[2, Inf)⁠) Number of inner replicates of each outer replicate with interval = "two-stage". The function checks the value for every interval, and uses it only with "two-stage".

Details

Optimism correction

On the development persons, the measures of a ladder are apparent. They are too optimistic, because fit_ladder() fitted the models on the same persons. The bootstrap optimism correction of Efron (1983); Harrell Jr. et al. (1996) estimates this optimism and subtracts it. The correction uses all persons. For a highly overfitted model, the correction can be too small (Efron 1983).

Let $D$ be the development persons and $f$ the ladder from $D$. Let $\theta(f, E)$ be a measure of $f$ on the persons $E$. In each replicate $b = 1, \dots, B$, the function draws a sample $D_b$ of the persons of $D$ with replacement. It fits every model of the ladder again on $D_b$, which gives $f_b$. The optimism of replicate $b$ and the corrected estimate are

$O_b = \theta(f_b, D_b) - \theta(f_b, D), \qquad \hat\theta_\text{corr} = \theta(f, D) - \frac{1}{B} \sum_{b} O_b.$

The refit repeats every step of fit_ladder(), also the knots of each spline, which come from the values of $D_b$.

Let $n$ be the number of development persons. The corrected value estimates the expected performance on new persons of a ladder that fit_ladder() fits on $n$ persons of the same population. It is a poor estimate of the performance of the one ladder $f$. Bates et al. (2024) prove this for the bootstrap and for cross-validation in the linear model. A report of the corrected estimate for $f$ assumes that $f$ is a typical ladder of that size.

Choices made in advance

The correction repeats only the steps of fit_ladder(). Choose the rung, spline_df, the thresholds, and the pairs before the analysis. A choice made from the results of the same persons is not part of the replicates, and the correction is then too small. The construction of each score is also outside the replicates. A score tuned on the same persons needs other persons for its validation (Choi et al. 2020).

The column heuristic_shrinkage of fit_ladder()$table shows the size of the optimism before any bootstrap. A value near 1 shows a small optimism.

Random numbers

Replicate $b$ draws its persons from stream $b$ of the L'Ecuyer-CMRG generator (L'Ecuyer et al. 2002). Stream 1 comes from base::set.seed() with seed, and stream $b + 1$ from parallel::nextRNGStream() applied to stream $b$. Each function that gets the plan starts again from seed. Because of this, all functions use the same bootstrap samples. Each function sets the generator kinds explicitly, and it restores the generator kinds and state of the session on exit.

When seed is NULL, bootstrap_plan() draws one seed from the current generator and stores it. This draw advances the generator state of the session. Every function that gets the plan then uses the stored seed.

Intervals

interval selects the 95% confidence interval of the corrected estimates:

• "analytic" moves the influence-function interval of the apparent estimate to the corrected estimate. Only the measures of one model with an influence function, such as the AUC and the Brier score, get this interval. Other measures, such as a difference between two models, the observed-to-expected ratio, and the calibration slope, get NA limits. The interval ignores the variance of the optimism and of the fitted models.

• "two-stage" is the interval of Algorithm 2 of Noma et al. (2021). Each outer replicate $b$ repeats the full correction on $D_b$, with n_inner inner replicates from the substreams of stream $b$ (parallel::nextRNGSubStream()). The limits are the 2.5% and 97.5% quantiles of the n_boot outer corrected estimates, and the standard error is their standard deviation. Every corrected measure gets this interval.

• "none" gives no interval.

The limits of "two-stage" are quantiles of the outer estimates. Their Monte Carlo error depends on n_boot. With 200 outer replicates, the Monte Carlo standard deviation of a limit is about 0.18 standard deviations of the outer estimates. With 1000 outer replicates, it is about 0.084. Because of this, use n_boot of 1000 or more with "two-stage" (Noma et al. 2021). Noma et al. (2021) found that the limits changed little from 100 to 2000 inner replicates (their e-Appendix B). The default n_inner = 50L is below that range. In one simulation setting, the analytic interval of the corrected AUC covered 0.93 to 0.94 for the rungs with a score of moderate effect. It covered 0.89 to 0.90 for the rungs with an AUC near 0.5, because the corrected AUC of such a rung is too high. The analytic interval of the Brier score covered 0.96 to 0.97. validation gives the setting. The coverage of the two-stage interval and of the other settings is not measured.

Failed replicates

A replicate fails at every horizon when a spline term or a refit fails on $D_b$. It fails at one horizon when the persons of $D_b$ with a prediction fail the horizon guard of horizon_support(). It also fails at one horizon when $f_b$ gives no risk for a person that $f$ evaluates. The follow-up limits of $f_b$ can be narrower than those of $f$ (see the section "Missing risks" of predict_risk()). The functions drop each failed replicate and give one warning. A measure can also be NA in one replicate. The optimism of a measure then uses the replicates in which both of its terms are not NA. When fewer than half of the n_boot replicates remain, the corrected estimate and its interval are NA.

Time

Each function refits the ladder in each replicate. A test on one computer used 20,000 persons, four models, and two horizons. A replicate took about 0.4 seconds in assess_discrimination() and about 0.7 seconds in assess_calibration(). 200 replicates then take a few minutes. With competing events, the cumulative incidence takes most of the time. A replicate of assess_discrimination() then took about 13 seconds, and 200 replicates take about 45 minutes. assess_stability() needs the same time. A replicate of assess_thresholds() or assess_net_benefit() took about 11 seconds. The corrected measures of assess_calibration() use the model of T2 only. Its replicate then took about 0.4 seconds. With 200 replicates in each of the five functions, the internal validation then takes about 3 hours. Each function repeats the refits and the risks, because each function runs its own replicates. The interval "two-stage" takes n_boot * (n_inner + 1) replicates. It is fast in a small cohort, where it matters most.

Value

An object of class prs_bootstrap_plan, a list with the elements n_boot (integer), seed (integer), interval (character), and n_inner (integer).

References

Bates S, Hastie T, Tibshirani R (2024). “Cross-validation: what does it estimate and how well does it do it?” Journal of the American Statistical Association, 119(546), 1434–1445. doi:10.1080/01621459.2023.2197686.

Choi SW, Mak TS, O'Reilly PF (2020). “Tutorial: a guide to performing polygenic risk score analyses.” Nature Protocols, 15(9), 2759–2772. doi:10.1038/s41596-020-0353-1.

Efron B (1983). “Estimating the error rate of a prediction rule: improvement on cross-validation.” Journal of the American Statistical Association, 78(382), 316–331. doi:10.1080/01621459.1983.10477973.

Harrell Jr. FE, Lee KL, Mark DB (1996). “Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.” Statistics in Medicine, 15(4), 361–387. doi:10.1002/(SICI)1097-0258(19960229)15:4<361::AID-SIM168>3.0.CO;2-4.

L'Ecuyer P, Simard R, Chen EJ, Kelton WD (2002). “An object-oriented random-number package with many long streams and substreams.” Operations Research, 50(6), 1073–1075. doi:10.1287/opre.50.6.1073.358.

Noma H, Shinozaki T, Iba K, Teramukai S, Furukawa TA (2021). “Confidence intervals of prediction accuracy measures for multivariable prediction models based on the bootstrap-based optimism correction methods.” Statistics in Medicine, 40(26), 5691–5701. doi:10.1002/sim.9148.

See Also

validation for the three modes of the ⁠assess_*()⁠ functions, and assess_discrimination(), assess_calibration(), assess_thresholds(), assess_net_benefit(), and assess_stability(), which accept the plan.

Examples

#----------------------------------------------------------------------------
# bootstrap_plan() example
#----------------------------------------------------------------------------
library(prsladder)

dat <- simulate_prs_data(n = 1000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(onset = "onset_score", progression = "progression_score")
)
ladder <- fit_ladder(cohort)

# A small number of replicates keeps the example fast. An analysis needs
# 200 or more.
plan <- bootstrap_plan(n_boot = 20L, seed = 1L)
discrimination <- assess_discrimination(ladder, 10, validation = plan)
discrimination$models
discrimination$optimism

---

Test the Entry Terms and Pick a Design for Absolute Risks

Description

check_design() tests how the hazard of T2 depends on the entry into a cohort under the recruitment index. From these tests, it picks a value of the argument design of fit_ladder() for absolute risks. The tests are for T2 and use the model with all scores.

The pick is for absolute risks only. Fit the hazard ratios and the primary tests with the default design = "t1". Choose this workflow before the analysis, and report the three rows of tests.

Usage

check_design(cohort, spline_df = c(age_t1 = 3L), entry_cuts = c(1, 2, 4, 7))

Arguments

cohort	(prs_cohort) A cohort from prs_cohort() under the recruitment index, with incident and prevalent persons and at least one event.
spline_df	(named integer() or NULL; default c(age_t1 = 3L); range ⁠[1, 8]⁠) Degrees of freedom of the variables with a spline term, as in fit_ladder(). Each name must be age_t1 or a numeric covariate of the cohort. NULL or a numeric vector of length zero gives a linear term to every variable. The function also accepts doubles with no fractional part, such as 3.
entry_cuts	(numeric(); default c(1, 2, 4, 7); range ⁠(0, Inf)⁠) Upper limits of the bands of the time since entry, in the time unit of the cohort, without the last band to infinity. The values must be finite and strictly increasing.

Details

Tests

Let $p$ be 1 for a prevalent person and 0 for an incident person, and let $d$ be t1_duration. Let $t$ be the time since T1, and let $u = t - d$ be the time since entry. The element tests has one row for each test below.

• entry_terms is the likelihood-ratio test of the model of "t1" inside the model of "t1_entry". It tests the terms $\gamma p + \delta d$ of "t1_entry". It has 2 degrees of freedom.

• entry_bands is a score test at the estimates of "t1_entry". The alternative replaces $\gamma p$ by $p \sum_j \gamma_j 1(u \in B_j)$. The bands are $B_j = (c_{j-1}, c_j]$, with $c_0 = 0$ and the values of entry_cuts, and a last band to infinity. The test asks whether the shift of prevalent persons changes with the time since entry. Because the test needs no fit of the alternative, a band with no T2 event needs no special case. The degrees of freedom are the rank of the information matrix of the band terms. This rank is usually one less than the number of bands that hold a prevalent person at the time of a T2 event. A band with no such person adds no information.

• index_duration is the test of survival::cox.zph() with transform = "km" for the duration term of "index" (Grambsch and Therneau 1994). It tests that the effect of $d$ in this model is constant in $u$. Its degrees of freedom are the value of age_t1 in spline_df.

A competing event is a censoring in each test.

To try the tests on data with a known effect of the entry, use the arguments prevalent_hazard_ratio, duration_hazard_ratio, and prevalent_effect_years of simulate_prs_data().

Rule

The function applies this rule at the 5% level:

1. If entry_bands does not reject, the pick is "t1_entry".

2. If entry_bands rejects, the pick is "index".

The tests entry_terms and index_duration are diagnostics, and the rule does not use them. The rule cannot pick "t1". A rejection of entry_terms shows that "t1" is not correct. It does not show that "t1_entry" is adequate. With the birth year in the models and a short recruitment period, entry_terms can have little power against an effect that ends after entry. The birth year and the age at T1 then explain most of t1_status and t1_duration. When the pick is "index", a rejection of index_duration shows that its assumption of a constant effect of $d$ can fail.

Pre-specify entry_cuts. The p-value of entry_bands assumes that the cut points do not depend on the data. Use the same spline_df in check_design() and fit_ladder().

Limits

The rule is for one cohort. Because the pick varies by chance, a study that pools results over cohorts pre-specifies one design for all cohorts. The optimism correction takes the design as fixed. No test of this function targets a score effect that differs by T1 status. Below about 500 events, a result without a rejection says little. One event in a band with very little follow-up can decide entry_bands. The rule was not studied with competing events or with an entry lag. The pick comes from T2, and the competing model uses the same design.

Errors

The function stops with an error in these cases:

• cohort is not under the recruitment index. Every person then has the same entry.

• All persons of the cohort have the same value of t1_status.

• The cohort has no events.

• A value of t1_status has 10 persons or fewer, as in fit_ladder().

• The data do not support a spline term. A smaller value in spline_df can fix the problem (see the section "Errors and warnings" of fit_ladder()).

• survival::coxph() would reduce the time since entry of a person to zero.

• One of the three fits gives a warning or an aliased coefficient, as in fit_ladder(). The error message names the design.

• entry_bands has no degrees of freedom. The band terms then have no information at the T2 event times. An example is a cohort with no prevalent person in follow-up after the first cut at the time of a T2 event.

• survival::cox.zph() fails, as in fit_ladder(). The error message names the design.

Value

A list with the two elements below.

tests

A data frame with one row per test, in the order entry_terms, entry_bands, and index_duration, and the columns below.

test

The name of the test.

null_model

The value of design in fit_ladder() of the model under the null hypothesis.

statistic

The chi-square statistic.

df

Integer degrees of freedom.

p_value

The p-value from the chi-square distribution.

rejected

Logical, TRUE when p_value is below 0.05.

design

The pick of the rule, "t1_entry" or "index". It is a value for the argument design of fit_ladder().

References

Grambsch PM, Therneau TM (1994). “Proportional hazards tests and diagnostics based on weighted residuals.” Biometrika, 81(3), 515–526. doi:10.1093/biomet/81.3.515.

See Also

fit_ladder() for the designs, check_entry() for the test of quasi-independence, and predict_risk() for the risks.

Examples

#----------------------------------------------------------------------------
# check_design() example
#----------------------------------------------------------------------------
library(prsladder)

dat <- simulate_prs_data(n = 2000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(onset = "onset_score", progression = "progression_score"),
  index = "recruitment",
  age_at_recruitment = "age_at_recruitment"
)

# The tests and the pick of the rule.
checks <- check_design(cohort)
checks$tests

# The ladder of the absolute risks.
risk_ladder <- fit_ladder(cohort, design = checks$design)
risk_ladder

---

Test the Quasi-Independence of Entry and Exit

Description

check_entry() tests whether the entry time of a cohort is quasi-independent of its exit time. Delayed entry corrects for the selection of prevalent persons only under this condition (see the section "Time scale and entry" of prs_cohort()). The test is the conditional Kendall's tau of Tsai (1990).

Usage

check_entry(cohort)

Arguments

cohort

(prs_cohort) A cohort from prs_cohort() with two or more distinct entry times, such as a cohort under the recruitment index.

Details

Statistic

Let $L_i$ be the entry time start of person $i$, and $X_i$ the exit time stop, both in time since T1. An exit is the first of T2 and the competing event. A censoring is not an exit. For each exit $k$, the risk set is $R_k = \{j : L_j < X_k \le X_j\}$, which includes $k$. A person who enters at $X_k$ is not in $R_k$, as in the follow-up intervals ⁠(start, stop]⁠ of prs_cohort(). The statistic is

$K = \sum_k \sum_{j \in R_k} \text{sign}(L_j - L_k),$

where the outer sum is over the exits. Each term compares two persons who were both in follow-up at $X_k$. These pairs are the comparable pairs of Tsai (1990). $M$ is the number of comparable pairs, and the estimate of tau is $K / M$. A positive tau means that persons who enter later also exit later.

Two exits at the same time are not a comparable pair. Each is in the risk set of the other, and their two terms cancel in $K$. $M$ does not count them. A pair with the same entry time adds zero to $K$, but it counts in $M$.

Variance

Let $x$ be a time with $d_x$ exits, and let $R_x$ be its risk set, with $r_x$ persons. Under quasi-independence, the $d_x$ persons who exit at $x$ are equally likely to be any $d_x$ persons of $R_x$. The terms of the times are then uncorrelated. The variance of $K$ is

$V = \sum_x \frac{d_x (r_x - d_x)}{r_x - 1}\, \frac{r_x^3 - \sum_g t_g^3}{3 r_x}.$

Here $t_g$ is the size of a group of equal entry times in $R_x$. The factor $(r_x - d_x) / (r_x - 1)$ is the correction for a sample without replacement. It is one when $d_x$ is one. Without equal entry times and without tied exits, each term is $(r_x^2 - 1) / 3$, which is the variance of Tsai (1990). Equal entry times and tied exits make the variance smaller. Under the recruitment index, every incident person has the same entry time. These groups can then be large. The statistic is $z = K / \sqrt{V}$, and the two-sided p-value comes from the standard normal distribution.

Interpretation

A small p-value is evidence that entry and exit are dependent. The absolute risks of cumulative_risk() and predict_risk() can then be biased. A large p-value does not show that the condition holds. The test has low power against a dependence that is not monotone (Chiou et al. 2018). The test is marginal. The condition can hold within groups of persons and fail in the full cohort, or the reverse. Because of this, the test can reject when a covariate, such as the age at T1, changes both the entry and the hazard. The designs "t1" and "t1_entry" of fit_ladder() need the condition only given their terms. The estimates of cumulative_risk() need the condition within each group. check_design() tests the entry given the terms of the models. The test also assumes that censoring is independent of the exit time, given the entry time. Censoring can depend on the entry time, as it does when a biobank ends follow-up on one date. With times in whole years, the test is approximate when censorings fall in the same years as exits. Use ages with decimals (see prs_cohort()).

The time of the computation grows with the number of persons times the number of distinct exit times.

Errors

The function stops with an error when cohort is not a prs_cohort. It also stops when every person has the same entry time, because the cohort then has no delayed entry to test. This case includes the diagnosis index and prevalent = "drop".

Value

A data frame with one row and the columns below.

n

Integer number of persons.

n_exits

Integer number of exits, T2 or the competing event.

n_pairs

Number of comparable pairs $M$.

tau

The estimate $K / M$, or NA when $M$ is zero.

z

The statistic $K / \sqrt{V}$, or NA when $V$ is zero.

p_value

The two-sided p-value of z, or NA when z is NA.

References

Chiou SH, Qian J, Mormino E, Betensky RA (2018). “Permutation tests for general dependent truncation.” Computational Statistics & Data Analysis, 128, 308–324. doi:10.1016/j.csda.2018.07.012.

Tsai W (1990). “Testing the assumption of independence of truncation time and failure time.” Biometrika, 77(1), 169–177. doi:10.1093/biomet/77.1.169.

See Also

prs_cohort() for the entry times of the two indexes, and check_design() for tests of the entry given the terms of the models.

Examples

#----------------------------------------------------------------------------
# check_entry() example
#----------------------------------------------------------------------------
library(prsladder)

dat <- simulate_prs_data(n = 2000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(progression = "progression_score"),
  index = "recruitment",
  age_at_recruitment = "age_at_recruitment"
)

# With its default arguments, the simulator recruits independently of the
# event times.
check_entry(cohort)

---

Compare Nested Models of a Ladder

Description

compare_ladder() tests whether the scores that a full model adds to a nested, reduced model improve the fit. For example, the pair c("onset+progression", "onset") tests whether the progression score adds information to the onset score. The function uses the likelihood-ratio test on the development persons. For each pair, it also reports the change in AIC and in Harrell's concordance.

Usage

compare_ladder(ladder, pairs = NULL)

Arguments

ladder

(prs_ladder) A ladder from fit_ladder().

pairs

(list of character(2), or NULL; default NULL) Pairs of models to test. Each element is c(full, reduced), with the full model first. Each pair must obey these rules: Both names are in names(ladder$models), and neither is NA. The reduced model is nested in the full model. The pair occurs only once in the list. The function ignores the names of the list and of its elements. If NULL, the function tests every nested pair. The pairs are then in order of the full model, then of the reduced model, both in the order of ladder$models.

Details

Nested pairs

Two models form a nested pair when the scores of the reduced model are a proper subset of the scores of the full model. For example, onset is nested in onset+progression, and base is nested in every other model. The models onset and progression do not form a nested pair.

A ladder with $k$ scores has $3^k - 2^k$ nested pairs. For example, two scores give 5 pairs and three scores give 19 pairs. A pair with base as the reduced model gives the gain of the scores over the covariates and age term alone. The PRS reporting standards ask for this gain over a conventional risk model (Wand et al. 2021).

Likelihood-ratio test

For the partial log-likelihoods $\ell_\text{full}$ and $\ell_\text{reduced}$, the statistic is

$\chi^2 = 2(\ell_\text{full} - \ell_\text{reduced}).$

The degrees of freedom equal the number of added scores, because each score has one coefficient. Under the null hypothesis that the added scores have no effect, the statistic has an asymptotic chi-squared distribution. The test is valid because every model of a ladder has the same persons, covariates, age term, and design. For the pre-specified tests, use a ladder with the default design = "t1" of fit_ladder() (see its section "Designs").

The statistic does not change when the analyst rescales a score or reverses its sign. Because of this, the standardization of the scores has no effect on the test. Numerical error in the fits can give a statistic slightly below zero. The function reports that statistic as computed, with a p-value of 1.

With competing events (see prs_cohort()), chisq, delta_aic, and delta_concordance compare the models of the cause-specific hazard of T2. They do not use the competing models of fit_ladder(). The test is then a test of the association of the added scores with the hazard of T2. A score can also change the cumulative incidence of T2 through the hazard of the competing event.

Selection on the first diagnosis can bias the hazard ratio of a score, most of all an onset score. The section on index event bias in fit_ladder() explains this bias. Under this bias, the test applies to the hazard ratio in the selected persons, which can differ from the hazard ratio in the population.

Descriptive measures

The function gives two descriptive differences, full minus reduced:

• delta_aic is the difference in AIC. It equals $-\chi^2 + 2\,\mathrm{df}$. A negative value favors the full model.

• delta_concordance is the difference in Harrell's concordance on the development persons. A positive value favors the full model.

Unlike the AIC, delta_concordance has no penalty for the added scores. Because it uses the development persons, it tends to overstate the gain in new persons. An accurate estimate of the gain needs out-of-sample persons (Choi et al. 2020). Use assess_discrimination() with a plan in validation or a different cohort in newdata for this estimate (see validation).

The function does not test these differences. On the development persons, a test of a difference in AUC between nested models is conservative (Demler et al. 2012; Pepe et al. 2013). Harrell's concordance is a form of the AUC for survival data. Pepe et al. (2013) study a binary outcome. They show that a null hypothesis of no gain in prediction equals a null hypothesis of no association of the added predictor. The association is conditional on the other predictors. The likelihood-ratio test is a test of this association.

Multiple comparisons

The function does not adjust p-values for multiple comparisons. Prespecify the primary pairs and give them in pairs. The default, all nested pairs, is a descriptive overview.

Value

A data frame with one row per pair, in the order of pairs, and the columns below.

full, reduced

Model names of the pair.

added

Scores that the full model adds, joined by +.

df

Integer number of added scores.

chisq

Likelihood-ratio statistic.

p_value

P-value of the test, not adjusted for multiple comparisons.

delta_aic

AIC of the full model minus AIC of the reduced model.

delta_concordance

Harrell's concordance of the full model minus that of the reduced model.

References

Choi SW, Mak TS, O'Reilly PF (2020). “Tutorial: a guide to performing polygenic risk score analyses.” Nature Protocols, 15(9), 2759–2772. doi:10.1038/s41596-020-0353-1.

Demler OV, Pencina MJ, D'Agostino Sr. RB (2012). “Misuse of DeLong test to compare AUCs for nested models.” Statistics in Medicine, 31(23), 2577–2587. doi:10.1002/sim.5328.

Pepe MS, Kerr KF, Longton G, Wang Z (2013). “Testing for improvement in prediction model performance.” Statistics in Medicine, 32(9), 1467–1482. doi:10.1002/sim.5727.

Wand H, Lambert SA, Tamburro C, Iacocca MA, O'Sullivan JW, Sillari C, Kullo IJ, Rowley R, Dron JS, Brockman D, Venner E, McCarthy MI, Antoniou AC, Easton DF, Hegele RA, Khera AV, Chatterjee N, Kooperberg C, Edwards K, Vlessis K, Kinnear K, Danesh JN, Parkinson H, Ramos EM, Roberts MC, Ormond KE, Khoury MJ, Janssens ACJW, Goddard KAB, Kraft P, MacArthur JAL, Inouye M, Wojcik GL (2021). “Improving reporting standards for polygenic scores in risk prediction studies.” Nature, 591(7849), 211–219. doi:10.1038/s41586-021-03243-6.

See Also

fit_ladder() for the models, and assess_discrimination() for the change in AUC and Brier score at a horizon, corrected for optimism or on a different cohort.

Examples

#----------------------------------------------------------------------------
# compare_ladder() example
#----------------------------------------------------------------------------
# Fit a ladder with three scores.
dat <- simulate_prs_data(n = 2000L, seed = 1L, score_correlation = 0.3)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(
    onset = "onset_score",
    progression = "progression_score",
    outcome = "outcome_score"
  )
)
ladder <- fit_ladder(cohort)

# Test the prespecified primary pairs.
compare_ladder(
  ladder,
  pairs = list(
    c("onset+progression", "onset"),
    c("progression+outcome", "outcome")
  )
)

# Describe all 19 nested pairs.
compare_ladder(ladder)

---

Cumulative Risk by Score Group

Description

cumulative_risk() compares the observed risk of the second event (T2) between persons with low and high values of one polygenic risk score. The function divides the cohort into quantile groups of the score. For each group, it estimates the cumulative risk of T2 at given times since the first diagnosis (T1). It also gives the risk difference and the risk ratio to a reference group, and the share of all T2 cases in each group. The other results are the restricted mean time lost to T2 and a test of equal risk across the groups. With score = NULL, it estimates the risk of the whole cohort.

The result is a prs_cumulative_risk object. plot_risk(), plot_stacked(), and plot_gradient() draw it, and describe_risk() makes a table of it. Because these functions use one object, the figures and tables show the same estimates.

The estimates are Kaplan-Meier estimates, or Aalen-Johansen estimates with competing events, with no covariate adjustment. For the risk from a fitted model, use predict_risk().

Usage

cumulative_risk(
  cohort,
  score,
  times,
  groups = 5L,
  cut_probs = NULL,
  reference = 1L,
  min_events = 1L,
  min_at_risk = 50L
)

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

Arguments

cohort	(prs_cohort) A cohort from prs_cohort().
score	(character(1) or NULL) The name of one score of the cohort, as given in the names of scores in prs_cohort(). If NULL, the function estimates the risk of the whole cohort. The argument has no default. Because of this, a missing score is an error.
times	(numeric(); range ⁠(0, Inf)⁠) Distinct finite times since T1, in the units of the age columns. Each time must be after the origin, the entry_lag of the cohort. The order does not matter. The result gives the times in increasing order.
groups	(integer(1); default 5L; range ⁠[2, Inf)⁠) Number of quantile groups of the score. For example, 5L gives quintiles and 4L gives quartiles. Must be missing when score is NULL or when cut_probs is given.
cut_probs	(numeric() or NULL; default NULL; range ⁠(0, 1)⁠) Strictly increasing probabilities at which to cut the score into length(cut_probs) + 1 groups. Each value is a probability, not a score in standard deviation units. For example, 0.95 separates the 5% of persons with the highest scores. The probabilities refer to the standardized score, which includes the sign from prs_cohort(). With a sign of -1, cut_probs = 0.95 separates the 5% of persons with the lowest raw scores. If NULL, the function uses groups groups of near-equal size. Must be NULL when score is NULL. See the section "Groups at cut points" in Details.
reference	(integer(1); default 1L; range ⁠[1, groups]⁠, or ⁠[1, length(cut_probs) + 1]⁠ with cut_probs) The group whose risk the function subtracts from the risk of each group. The default is the group with the lowest scores. Must be missing when score is NULL.
min_events	(integer(1); default 1L; range ⁠[1, Inf)⁠) Minimum number of T2 events in a group at or before the time. The default hides a risk of 0 with an interval of $[0, 0]$, which occurs only with competing events. Without competing events, the default has no effect, because the third condition in Details requires an event. A larger value selects the groups with early events. The risks that pass are then too high, and their intervals cover less than 95%.
min_at_risk	(integer(1); default 50L; range ⁠[1, Inf)⁠) Minimum number of persons of a group at risk at the time and at each earlier event time of the group.
x	(prs_cumulative_risk) A result of cumulative_risk().
...	Not used.

Details

Score groups

The function cuts the score at its sample quantiles into groups groups of near-equal size. The cut uses the standardized score of the cohort, which includes the sign from prs_cohort(). Group 1 has the lowest standardized scores. With a sign of -1, group 1 has the highest raw scores. The function stops with an error when tied quantiles give fewer distinct groups than requested.

Groups at cut points

cut_probs gives groups of other sizes. The breaks are quantile(x, c(0, cut_probs, 1), type = 7), where x is the standardized score. Group $k$ holds the scores in the interval $(q_{k-1}, q_k]$ of the breaks, and group 1 also holds its lower limit. groups = 5L and cut_probs = 1:4 / 5 give identical breaks. The literal c(0.2, 0.4, 0.6, 0.8) also gives identical breaks. seq(0.2, 0.8, 0.2) can give breaks that differ by about $10^{-16}$, because its third element is $0.6 + 10^{-16}$. A score in that small gap goes to another group.

Two examples:

• cut_probs = 0.95 compares the 5% of persons with the highest scores (group 2) with the other 95% (group 1).

• cut_probs = c(0.05, 0.95) with reference = 2L compares each tail with the middle 90%.

The labels of the figures and tables of the result give the percent range of each group, such as ⁠0-95%⁠ and ⁠95-100%⁠. These ranges refer to the standardized score. With a sign of -1, cut_probs = 0.95 separates the 5% of persons with the lowest raw scores. Each label then also gives the range of the raw score, such as ⁠95-100% (raw 0-5%)⁠. The footnote or caption states that the groups run against the raw score.

A small group has fewer persons at risk. Because of this, it fails the condition on min_at_risk at an earlier time. For example, with 2,000 persons, the top 5% is 100 persons, and the risks at later times can be NA. A group with fewer persons than min_at_risk fails the condition at every time. Each of its risks is then NA, and the function gives a warning. Because the case share needs the risks of all groups, each case share is then also NA.

Time scale and risk

The time scale is time since T1, with delayed entry at start. A person is at risk at time $t$ when start $< t \le$ stop. The prediction metrics, such as assess_calibration(), use time since entry. This function uses time since T1, because a progression curve starts at the diagnosis.

The risk at time $t$ is $1 - \hat{S}(t)$, where $\hat{S}$ is the Kaplan-Meier estimate from survival::survfit() on Surv(start, stop, status). The standard error is the Greenwood standard error of $\hat{S}(t)$, which is also the standard error of the risk. With a positive entry_lag in prs_cohort(), the curve is conditional on being free of T2 at the lag.

When the cohort declares the competing status code (see prs_cohort()), the risk at time $t$ is the cumulative incidence of T2. This risk is the probability of T2 by time $t$ before a competing event. The function estimates it with the Aalen-Johansen estimator (Aalen and Johansen 1978) from survival::survfit() on Surv(start, stop, state), where state is 0, 1, or 2. This estimator accepts delayed entry when the entry and censoring times are independent of the event times (Geskus 2011). The standard error is the infinitesimal jackknife standard error of survival::survfit(). With a positive entry_lag, the curve is conditional on being free of both events at the lag. One minus the Kaplan-Meier estimate is at least the cumulative incidence. The Kaplan-Meier estimate treats a competing event as censoring. It assumes that a person with a competing event keeps the hazard of T2 after that event. In fact, that person can no longer have T2. Because of this, one minus the Kaplan-Meier estimate overestimates the cumulative incidence (Putter et al. 2007; Andersen et al. 2012).

Support of a group at a time

A group supports a time when all three conditions below are true:

1. At least min_events events occur in the group at or before the time.

2. At least min_at_risk persons of the group are at risk at the time.

3. At least min_at_risk persons of the group are at risk at each event time of the group up to the time.

The events of the first condition are T2 events. With competing events, the event times of the third condition include the times of competing events, because these events also change the curve. The third condition protects against an early risk set of a few prevalent persons. Each step of the curve at such an early time stays in every later estimate. When a group does not support a time, risk and risk_se are NA, never zero. The competing risk uses the same conditions, with competing events in the first condition. A time after the last follow-up of a group has zero persons at risk, and its risk is NA.

The default min_events = 1L hides only a risk of 0 with an interval of $[0, 0]$. This interval does not include a positive risk. With competing events, a group with competing events and no T2 event by a time has this risk. A larger min_events selects the samples with more early events than expected, because the condition passes only in these samples. The risks that pass are then too high, and their intervals cover the risk less often than stated. A simulation used delayed entry and groups of 200 and 1000 persons. Under the default, the 95% intervals covered the risk in 93.6% to 96.7% of the samples in which the conditions passed. With min_events = 10L, they covered 53% to 89% at the times where the condition passed in 7% to 27% of samples. The default conditions also select samples at a time with about one expected event, because the third condition needs an event. In the simulation, such a risk of 1.6% was too high by 0.7 points on average, and its interval covered the risk in 95% of samples.

These counts differ from those of horizon_support(), which counts on time since entry and also applies a censoring guard.

A row of curves uses only the third condition. The function marks an event row as supported when at least min_at_risk persons are at risk at each event time of the group up to the row. At an event time, the second condition adds nothing to the third. The first condition does not apply to curves, because with an event minimum, a curve starts only in the samples with early events. The origin row has the value of the first event row. The reason is that the risk before the first event is 0 whatever the size of the risk set. The third condition cannot change from false to true. Because of this, the supported rows of a group are its first rows. curves keeps the risks of the rows that are not supported, and supported marks the rows that a figure draws.

The three conditions are statistical. They do not protect the privacy of persons. Each step of a curve is one event, and n_events gives the event counts. The package does no disclosure control. Apply the rules of the data provider before you share a result.

Curves

The element curves holds the step function of the risk of each group. It has one row at the origin and one row at each event time of either cause. The origin is the entry_lag of the cohort, and the risk at the origin is 0. The number at risk at the origin is the number of persons just after it, with start $\le$ origin $<$ stop. The rule start $< t$ gives an empty risk set at the origin. Under the recruitment index, the number at the origin is the number of incident persons. The number at risk can then increase over time, because prevalent persons enter after T1. At a time in times, the last row of curves at or before the time has the risk of risks, when risks shows a risk.

Intervals

The interval of a risk $F$ with standard error $s$ is a pointwise 95% interval on the log scale of $F$:

$\left[F e^{-z s / F},\ \min\left(1,\ F e^{z s / F}\right)\right], \quad z = \Phi^{-1}(0.975).$

Both limits are 0 when $F$ is 0. For example, before the first T2 event of a group, a row of curves has a risk of 0 and the interval $[0, 0]$. The same rule applies to the risk of T2 and to the competing risk, in risks and in curves. The interval at each time holds for that time only, not for the whole curve.

For the Aalen-Johansen estimate, this interval is the default of survival::survfit(). For the Kaplan-Meier estimate, the default of survival::survfit() is on the log scale of $1 - F$. The ggsurvfit and gtsummary packages use that default. Because of this, the Kaplan-Meier intervals of this function differ from the intervals of these packages.

Whole cohort

With score = NULL, all persons form group 1. groups and reference must then be missing, and cut_probs must be NULL. The columns of the difference, the ratio, the case share, and the difference in time lost are NA, and test has zero rows.

Differences and test

The difference is the risk of a group minus the risk of the reference group at the same time. Because the groups hold different persons, the variance of a difference is the sum of the two variances. The interval is a pointwise 95% Wald interval, $d \pm z_{0.975} \text{SE}(d)$, with no adjustment for multiple comparisons. A difference is NA when either of the two risks is NA.

The test is the score test of coxph(Surv(start, stop, status) ~ group), with the group as a factor. It has groups - 1 degrees of freedom. This test accepts delayed entry. It equals the log-rank test when there is no delayed entry and no tied event times. The test uses all follow-up and does not depend on times, reference, or the guards.

With competing events, the test compares the cause-specific hazards of T2, which agrees with the models of fit_ladder(). A competing event counts as censoring in the test. The curves of the groups can also differ because of the hazard of the competing event, which the test does not examine. The test of Gray (1988) compares the cumulative incidence curves. It applies only to right-censored data without delayed entry.

Ratios and case shares

The ratio is the risk of a group divided by the risk of the reference group at the same time:

$\widehat{RR}_g = \hat F_g / \hat F_r, \qquad \mathrm{SE}(\log \widehat{RR}_g) = \sqrt{(s_g / \hat F_g)^2 + (s_r / \hat F_r)^2},$

where $\hat F_g$ is the risk of group $g$, $s_g$ is its standard error, and $r$ is the reference group. This standard error comes from the delta method on the log scale. The covariance is zero, because the groups hold different persons. The interval is a pointwise 95% Wald interval on the log scale, $\widehat{RR}_g \exp(\pm z_{0.975} \mathrm{SE})$. The column log_ratio_se holds the standard error of $\log \widehat{RR}_g$, not of the ratio. A ratio is NA for the reference group, and when either of the two risks is NA. Because min_events is at least 1, each risk that the function shows is positive, and the logarithm of each ratio is finite.

The case share is the share of all T2 cases by the time that occur in a group. With $n_g$ the number of persons in group $g$, $a_g = n_g \hat F_g$, $v_g = n_g^2 s_g^2$, and $S = \sum_k a_k$:

$\hat\pi_g = a_g / S, \qquad \widehat{\mathrm{Var}}(\hat\pi_g) = \frac{(S - a_g)^2 v_g + a_g^2 \sum_{k \ne g} v_k}{S^4}.$

The variance comes from the delta method. The interval is a pointwise 95% interval on the logit scale, $\mathrm{logit}^{-1}\{\mathrm{logit}(\hat\pi_g) \pm z_{0.975} \mathrm{SE} / [\hat\pi_g (1 - \hat\pi_g)]\}$. The shares of the groups sum to 1 at each time. For the persons above the cut point, the case share is an analogue of the proportion of cases followed of Pfeiffer and Gail (2011). Pfeiffer and Gail rank the persons by the predicted risk of a model, and they compute the proportion from the predicted risks. Here, the score ranks the persons, and the share comes from the observed events. At a time where any group has an NA risk, the case shares of all groups are NA.

The case share does not use the raw counts of events. Because censoring removes more of the later events, the raw counts give too large a share to a group whose events occur earlier. The top group is such a group. In a simulation, the raw counts overstated the share of the top quintile by 0.42 percentage points, and $n_g \hat F_g$ by 0.14 points.

The intervals of the ratio and the case share are conditional on the groups of the sample. They treat the cut points and the group sizes as fixed, but the cut points are sample quantiles that change from sample to sample.

The case share describes a population with the group sizes of the cohort. Under the recruitment index, the cohort holds only the persons who were free of T2 at recruitment. Persons with high scores more often had T2 before recruitment. Because of this, the cohort holds fewer of them than the population at T1. The case share of the top group is then not an unbiased estimate of its share in the population at T1. In a simulation with the default hazard ratios of simulate_prs_data(), the share of the top quintile was 0.2646 in the persons of the cohort and 0.2681 in all persons at T1.

Time lost

The element time_lost gives the restricted mean time lost to T2 by each time $\tau$ in times (Andersen 2013; Conner and Trinquart 2021):

$\mathrm{RMTL}_g(\tau) = \int_{t_0}^{\tau} \hat F_g(u)\,du,$

where $t_0$ is the origin and $\hat F_g$ is the risk of group $g$. It is the mean time from T2 to $\tau$, where a person with no T2 by $\tau$ adds 0. The unit is the unit of the age columns. For example, take years as the unit and $\tau = 10$. A time lost of 0.5 then means that, in the first 10 years after T1, a person of the group spends 0.5 years after T2 on average. Without competing events, the time lost is $\tau - t_0$ minus the restricted mean survival time.

With competing events, competing_rmtl is the same area under the cumulative incidence of the competing event. The time lost to T2, the time lost to the competing event, and the restricted mean time free of both events sum to $\tau - t_0$.

The function gets the time lost and its standard error from summary(fit, rmean = tau) on the fit of the risks. The time lost equals the area under the step function of curves from the origin to $\tau$. With competing events, the standard error is the infinitesimal jackknife standard error of survival::survfit(). survival::survfit() gives this standard error at each time of the fit, and summary() interpolates it linearly between these times. In a check, the interpolated value differed from the exact value by less than 0.1%. Without competing events, the standard error is the Greenwood-type standard error of the restricted mean survival time from survival::survfit().

The interval is a pointwise 95% Wald interval, $\mathrm{RMTL} \pm z_{0.975} \mathrm{SE}$, with the lower limit truncated at 0. The difference is the time lost of a group minus the time lost of the reference group. Its variance is $s_g^2 + s_r^2$, and its interval is a 95% Wald interval. The time lost at $\tau$ is NA where the risk at $\tau$ is NA, because the conditions on min_at_risk cover each event time up to $\tau$. For the competing event, the same rule uses the competing risk.

The intervals are conditional on the groups of the sample, as are those of the ratio and the case share. A simulation with competing events used quintile groups, 300 samples of 3,000 persons, and $\tau = 10$. The reference value was the time lost in one sample of 400,000 persons. Without delayed entry, the 95% intervals of the time lost to T2 covered this value in 92.0% to 96.7% of the samples. With delayed entry, they covered it in 94.0% to 96.0% of the samples. The interval of the difference between the top and bottom quintiles covered its reference value in 94.7% and 94.0% of the samples.

Value

An object of class prs_cumulative_risk, a list with the elements below.

risks

A data frame with one row per group and time. The function orders the rows by group, then by increasing time. The columns are below. The columns n_competing and competing_risk to competing_risk_upper occur only when the cohort declares the competing status code.

group

Integer group number, from 1 to groups.

label

Character label of the group, as in the figures and tables. With equal groups, it is the group number, with the tags (lowest) and (highest) for a sign of 1. With cut_probs, it is the percent range of the standardized score, such as ⁠95-100% (highest)⁠, or ⁠95-100% (raw 0-5%)⁠ for a sign of -1. "All persons" for the whole cohort.

time

Time since T1.

n

Integer number of persons in the group.

n_events, n_competing, n_censored

Integer numbers of T2 events, competing events, and censorings in the group at or before time.

n_at_risk

Integer number of persons in the group at risk at time.

min_at_risk_before

Integer smallest number at risk at an event time of the group up to time, for an event of either cause. NA when the group has no event by time.

risk, risk_se

Cumulative risk of T2 at time and its standard error. NA when the group does not support time.

risk_lower, risk_upper

Limits of the 95% interval of the risk.

competing_risk, competing_risk_se, competing_risk_lower, competing_risk_upper

Cumulative incidence of the competing event at time, its standard error, and the limits of its 95% interval.

difference, difference_se

Risk difference from the reference group and its standard error. NA for the reference group.

difference_lower, difference_upper

Limits of the 95% Wald interval of the difference.

ratio, log_ratio_se

Risk ratio to the reference group and the standard error of its logarithm. NA for the reference group.

ratio_lower, ratio_upper

Limits of the 95% Wald interval of the ratio on the log scale.

case_share, case_share_se

Share of all T2 cases by time that occur in the group, and its standard error. NA at a time where any group has an NA risk, and for the whole cohort.

case_share_lower, case_share_upper

Limits of the 95% interval of the case share on the logit scale.

test

A data frame with the columns statistic, df, and p_value of the score test across groups. It has one row, or zero rows when score is NULL.

curves

A data frame with one row at the origin and one row at each event time of either cause, for each group. The function orders the rows by group, then by increasing time. The columns are group, time, n_at_risk, n_events, n_competing, the risk columns of risks, and supported. n_events and n_competing are totals at or before time. The risks are never NA. supported is TRUE for a row that a figure draws.

time_lost

A data frame with the rows of risks, in the same order, and the columns below. The columns competing_rmtl to competing_rmtl_upper occur only when the cohort declares the competing status code. The values are in the unit of the age columns.

group, label, time

As in risks. To merge time_lost with risks, use these three columns as the key.

rmtl, rmtl_se

Restricted mean time lost to T2 from the origin to time, and its standard error. NA where risk is NA.

rmtl_lower, rmtl_upper

Limits of the 95% Wald interval of the time lost, with the lower limit truncated at 0.

competing_rmtl, competing_rmtl_se, competing_rmtl_lower, competing_rmtl_upper

Restricted mean time lost to the competing event, its standard error, and the limits of its 95% interval. NA where competing_risk is NA.

rmtl_difference, rmtl_difference_se

Difference in time lost to T2 from the reference group and its standard error. NA for the reference group.

rmtl_difference_lower, rmtl_difference_upper

Limits of the 95% Wald interval of the difference.

settings

A list with the elements score, sign, groups, cut_probs, reference, times, min_events, min_at_risk, competing, index, origin, delayed_entry, and event_labels. sign is the sign of the score in prs_cohort(), or NULL for the whole cohort. groups is the number of groups, which is length(cut_probs) + 1 with cut_probs. cut_probs is the input cut_probs, or NULL for groups of near-equal size. competing is TRUE when the cohort declares the competing status code. origin is the entry_lag of the cohort. delayed_entry is TRUE when a person enters after the origin. event_labels is the event_labels element of the cohort, which the figures and tables use.

print() shows the guards, the risks of at most 20 rows of risks, a line that names time_lost, and the test. With cut_probs, it first gives the label of each group, one per line. With a sign of -1, it also states that the raw score runs against the groups. It returns x invisibly.

References

Aalen OO, Johansen S (1978). “An empirical transition matrix for non-homogeneous Markov chains based on censored observations.” Scandinavian Journal of Statistics, 5(3), 141–150.

Andersen PK (2013). “Decomposition of number of life years lost according to causes of death.” Statistics in Medicine, 32(30), 5278–5285. doi:10.1002/sim.5903.

Andersen PK, Geskus RB, de Witte T, Putter H (2012). “Competing risks in epidemiology: possibilities and pitfalls.” International Journal of Epidemiology, 41(3), 861–870. doi:10.1093/ije/dyr213.

Conner SC, Trinquart L (2021). “Estimation and modeling of the restricted mean time lost in the presence of competing risks.” Statistics in Medicine, 40(9), 2177–2196. doi:10.1002/sim.8896.

Geskus RB (2011). “Cause-specific cumulative incidence estimation and the Fine and Gray model under both left truncation and right censoring.” Biometrics, 67(1), 39–49. doi:10.1111/j.1541-0420.2010.01420.x.

Gray RJ (1988). “A class of K-sample tests for comparing the cumulative incidence of a competing risk.” The Annals of Statistics, 16(3), 1141–1154. doi:10.1214/aos/1176350951.

Pfeiffer RM, Gail MH (2011). “Two criteria for evaluating risk prediction models.” Biometrics, 67(3), 1057–1065. doi:10.1111/j.1541-0420.2010.01523.x.

Putter H, Fiocco M, Geskus RB (2007). “Tutorial in biostatistics: competing risks and multi-state models.” Statistics in Medicine, 26(11), 2389–2430. doi:10.1002/sim.2712.

See Also

• plot_risk(), plot_stacked(), plot_gradient(), and describe_risk() for figures and tables of the result.

• predict_risk() for the model-based risk from entry.

• horizon_support() for the support of prediction horizons.

Examples

#----------------------------------------------------------------------------
# cumulative_risk() example
#----------------------------------------------------------------------------
library(prsladder)

# Simulate a cohort and standardize its progression score.
dat <- simulate_prs_data(n = 2000L, seed = 1L)
cohort <- prs_cohort(
  dat,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(progression = "progression_score")
)

# Risk by score quartile at 5 and 10 years since the first diagnosis.
res <- cumulative_risk(cohort, "progression", times = c(5, 10), groups = 4L)
res
res$risks[c("group", "time", "risk", "difference")]

# Years lost to T2 in the first 5 and 10 years since the first diagnosis.
res$time_lost[c("group", "time", "rmtl", "rmtl_difference")]

# The top 5% of the score against the other 95%. The top 5% holds only 100
# persons, so the times are early enough for 50 persons to be at risk.
top <- cumulative_risk(
  cohort,
  "progression",
  times = c(2, 5),
  cut_probs = 0.95
)
top$risks[c("group", "time", "n", "risk", "difference")]

# Risk of the whole cohort.
cumulative_risk(cohort, NULL, times = c(5, 10))

---