---
title: "From type 2 diabetes to coronary artery disease: a prsladder tutorial"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{From type 2 diabetes to coronary artery disease: a prsladder tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)
has_gtsummary <- requireNamespace("gtsummary", quietly = TRUE)
```

```{css, echo = FALSE, eval = !identical(Sys.getenv("IN_PKGDOWN"), "true")}
/*
 * Floating table of contents for html_vignette.
 *
 * On a screen at least 1040px wide, the table of contents is fixed in a
 * sidebar to the left of the text, and it scrolls by itself when it is
 * taller than the window. The body (700px plus 2em of padding on each side)
 * moves right so that the sidebar never covers it, and stays centered on a
 * wide screen. On a narrow screen, the default inline table of contents of
 * vignette.css stays.
 *
 * The chunk does not run in a pkgdown build. The site has its own table of
 * contents, and the margin of `body` would move the page off its grid.
 */
@media (min-width: 1040px) {
  body {
    margin-left: max(270px, calc(50vw - 378px));
    margin-right: auto;
  }

  #TOC {
    position: fixed;
    top: 1em;
    left: max(20px, calc(50vw - 628px));
    width: 230px;
    max-height: calc(100vh - 2em);
    overflow-y: auto;
    margin: 0;
    box-sizing: border-box;
  }

  #TOC ul {
    padding-left: 30px;
  }
}
```

## The question

People with diabetes, including type 2 diabetes (T2D), develop coronary artery disease (CAD) at about twice the rate of people without diabetes (Emerging Risk Factors Collaboration 2010).
Some of that risk is inherited.
A polygenic risk score (PRS) sums many small genetic effects into one number per person.
A research team has three scores for every participant of a biobank:

| Name in the text | Built from a GWAS of | Column of the data | Name in the output |
|---|---|---|---|
| T2D score | T2D in the general population | `t2d_score` | `onset` |
| CAD score | CAD in the general population | `cad_score` | `outcome` |
| Progression score | Time from T2D to CAD, among people with T2D | `t2d_to_cad_score` | `progression` |

The team asks one question:

> Among people who already have T2D, does the progression score tell us who will develop CAD, beyond what the T2D score and the CAD score already tell us?

Wand et al. (2021) call this kind of study disease prognosis: the risk of a further outcome after a diagnosis.
prsladder is built for it.

prsladder starts after the scores exist.
It does not do the genotype quality control or build the scores; see Choi et al. (2020) for those steps.

This tutorial writes the analysis plan first, and then runs the full analysis in twelve steps:

1. Build the cohort, with follow-up from biobank recruitment and death as a competing event.
2. Describe it.
3. Fit a model for every combination of scores.
4. Test whether the progression score adds information.
5. Look at CAD risk by score group.
6. Choose a design for the absolute risks.
7. Predict each person's absolute risk.
8. Check the predictions on the development persons, corrected for optimism (internal validation).
9. Check the predictions in a second biobank (external validation).
10. Compare with an analysis that starts follow-up at the T2D diagnosis.
11. Compare with an analysis that treats deaths as censoring.
12. Answer the question.

It ends with a checklist for the report.

The code below loads the package.
It also sets the ggplot2 theme of every figure in this tutorial.

```{r}
library(prsladder)
ggplot2::theme_set(ggplot2::theme_bw())
```

## Terms and clocks

prsladder works for any first diagnosis and any later event, and it has general names for them.
Throughout the package, the first diagnosis (here T2D) is called T1 and the later event (here CAD) is called T2.

| Term of the package | Meaning in this tutorial | Where it appears |
|---|---|---|
| T1 | The T2D diagnosis | `age_t1`, `t1_status`, `t1_duration`, the designs `"t1"` and `"t1_entry"`, and the axis "Time since T1" |
| T2 | CAD | The help pages and the messages of the package |
| Competing event | Death before CAD | `competing_risk`, `n_competing`, and `event = "competing"` |
| Incident person | A person who joined the biobank before the T2D diagnosis | `t1_status` |
| Prevalent person | A person who joined the biobank after the T2D diagnosis | `t1_status` |
| Entry | The start of the follow-up of a person. In Step 1, it is the diagnosis for an incident person and recruitment for a prevalent person | `start`, and "since entry" in the output and on the axes |
| Horizon | A number of years after entry, for a predicted risk | `horizon` and `horizons` |

The package counts time on two clocks.

| Clock | Starts at | Used by |
|---|---|---|
| Years since the diagnosis | The T2D diagnosis | `times` of `cumulative_risk()` and of `standardized_risk()` |
| Years since entry | The entry of each person | `horizon` of `predict_risk()`, and `horizons` of `horizon_support()` and of every `assess_*()` function |

The two clocks agree for an incident person.
For a prevalent person who joined 3 years after the diagnosis, the risk within 10 years of entry is the risk from year 3 to year 13 after the diagnosis, for a person who is alive and free of CAD at year 3.
For a prevalent person, a risk from the diagnosis and a risk from entry are then two different quantities.
The figure of a function has the clock of that function.

## The data

Real biobank data cannot ship with a package, so the tutorial uses a simulated biobank.
`simulate_prs_data()` creates one row per person with T2D.
Because the true effects of the scores are set by us, we can check that the package recovers them.
We give the three scores hazard ratios for CAD of 1.10, 1.25, and 1.30 per standard deviation (SD).
The three scores are correlated, as real scores for related diseases are.

Each check against the true values starts with the words "Simulation check (not possible with real data)".

Because many people with T2D die before they develop CAD, the simulated biobank has deaths.
The T2D score also raises the rate of death, with a hazard ratio of 1.2 per SD.
The other two scores have no effect on death.

Age and sex change both rates, as they do in real data.
Each year of age at the T2D diagnosis raises the rate of CAD by 4% and the rate of death by 8%.
Women have 0.7 times the rate of men, for CAD and for death.
The principal components and the birth year have no effect.

The simulator sets the level of each rate for a reference person.
This person is a man with a T2D diagnosis at age 55 and average scores.
Without deaths, this man would have a 15% risk of CAD within ten years of the diagnosis.
Without CAD, he would have a 20% risk of death within ten years.
These values were chosen for illustration and are not estimates from any study.

```{r}
simulated <- simulate_prs_data(
  n = 3000L,
  seed = 2026L,
  hazard_ratios = c(t2d = 1.10, t2d_to_cad = 1.25, cad = 1.30),
  score_correlation = 0.3,
  baseline_risk = 0.15,
  incident_probability = 0.3,
  competing_risk = 0.2,
  competing_hazard_ratios = c(t2d = 1.2),
  covariate_hazard_ratios = c(age_at_diagnosis = 1.04, sex = 0.7),
  competing_covariate_hazard_ratios = c(age_at_diagnosis = 1.08, sex = 0.7)
)
```

As in most biobanks, many people had T2D before they joined.
`incident_probability = 0.3` makes 30% of the participants join before their T2D diagnosis (incident), and 70% join after it (prevalent).
With the defaults of the simulator, the biobank recruits from 2006 to 2010, and follow-up ends for every person at the end of the study, in 2022.

The simulator follows every person from the T2D diagnosis.
A real biobank cannot: it recruits only persons who are alive, and their earlier records come from registers.
Some simulated prevalent persons died before they would have joined, and the column `recruited` is `FALSE` for them.
We keep the recruited persons, so that the data are those of a real biobank.
Persons whose CAD came before recruitment are recruited, because the simulator ends follow-up at CAD and so has no deaths after CAD.

```{r}
# Keep the persons who were alive at recruitment.
biobank <- simulated[simulated$recruited, ]
sum(!simulated$recruited)
print(
  head(biobank[c("IID", "t2d_score", "t2d_to_cad_score", "cad_score", "sex",
                 "age_at_diagnosis", "age_at_exit", "exit_status",
                 "age_at_recruitment")], 3),
  digits = 3
)
table(biobank$exit_status)
```

The columns that matter are:

| Column | Meaning |
|---|---|
| `IID` | Participant identifier |
| `t2d_score`, `t2d_to_cad_score`, `cad_score` | The three raw scores |
| `sex`, `PC1` to `PC5`, `birth_year` | Covariates; `sex` is `Male` or `Female`, and the PCs are genetic principal components of ancestry |
| `age_at_diagnosis` | Age at the T2D diagnosis (T1) |
| `age_at_exit` | Age at CAD, at death, or at the end of follow-up |
| `exit_status` | 1 if CAD was observed, 2 if the person died before CAD, 0 if follow-up ended first |
| `age_at_recruitment` | Age at which the person joined the biobank |

The simulated columns have the names of the arguments of `prs_cohort()`, but real data can use other names: each argument takes the name of a column.

The simulator also returns `true_risk`, each person's true risk of CAD within ten years of the T2D diagnosis, which a real data set does not have.
Because a person can die first, this risk is a cumulative incidence (Step 2).

The analysis uses all persons of the biobank.
A model always looks better on the persons it was fitted to than on new persons.
Step 8 corrects for this optimism with the bootstrap, and Step 9 checks the models in a second biobank.
A random split of the biobank into a development part and a test part is not recommended.
The split gives a weaker model and an imprecise check (`?validation`).
The persons must be unrelated, as in any PRS analysis.

## The analysis plan

A choice made after looking at the results biases them, and the bootstrap of Step 8 does not correct for it (`?validation`).
Write the choices down before Step 1.
The choices of this tutorial are:

| Choice | Value | Reason | Steps |
|---|---|---|---|
| Start of follow-up | Recruitment index, on the clock of years since the T2D diagnosis | A biobank holds only persons alive at recruitment | 1, 10 |
| Competing event | Death before CAD | Death ends the chance of CAD | 1, 11 |
| Primary comparisons | `onset+progression` vs `onset`, and `progression+outcome` vs `outcome` | The question | 4, 8, 9 |
| Secondary comparison | `onset+progression+outcome` vs `onset+outcome` | Progression beyond both other scores, in the model of the report | 4, 8 |
| Decision rule | Both primary tests must reject at the 5% level. The secondary test follows only when both reject | The question asks for information beyond each score. This rule needs no adjustment for multiple tests | 4 |
| Prediction model of the report | `onset+progression+outcome` | It uses all three scores | 7 to 9 |
| Horizons | 5 and 10 years after entry | Clinical use, if the data support them (Step 7) | 7 to 9 |
| Risk thresholds | 10% and 20% at 10 years | Risk categories of a guideline (Step 8) | 8 |
| Range of the decision curve | 5% to 30% | Around the thresholds | 8 |
| Age term | Natural spline with 3 degrees of freedom | Default of `fit_ladder()` | 3 |
| Design of the hazard ratios | `"t1"` | The data must not choose the design of a primary test | 3, 4 |
| Design of the absolute risks | The pick of `check_design()`, with its default bands | The hazard after entry can differ between incident and prevalent persons | 6, 7 to 9 |
| Bootstrap | 200 or more replicates, fixed seed (30 here, to run fast) | Monte Carlo error (Step 8) | 8 |

## Step 1: build the cohort

`prs_cohort()` is the only function that reads your data frame, apart from `prediction_data()` for persons without outcome data (Step 7).
You tell it which column holds what.
The names you give the scores (`onset`, `progression`, `outcome`) are used in every later output.

```{r}
cohort <- prs_cohort(
  biobank,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(
    onset = "t2d_score",
    progression = "t2d_to_cad_score",
    outcome = "cad_score"
  ),
  covariates = c("sex", "PC1", "PC2", "PC3", "PC4", "PC5", "birth_year"),
  index = "recruitment",
  age_at_recruitment = "age_at_recruitment",
  status_codes = list(censored = 0, event = 1, competing = 2),
  event_labels = c(event = "CAD", competing = "Death")
)
cohort
```

`prs_cohort()` does three things:

- It checks the columns and stops with a clear message if, for example, an identifier repeats or an age is not numeric.
- It removes persons who cannot be analyzed, such as persons with a missing value, and counts them by reason (`cohort$n_excluded`).
- It standardizes each score to mean 0 and SD 1, so every hazard ratio below is per one SD of the score.

`index = "recruitment"` sets when each person enters follow-up.
The clock of the cohort is years since the T2D diagnosis.
Incident persons enter at the diagnosis.
Prevalent persons enter at recruitment, after the years they have already lived with T2D.
Statisticians call this delayed entry, or left truncation (Putter et al. 2007).
A prevalent person is in the biobank only because the person was alive and free of CAD at recruitment, and delayed entry accounts for this selection.
Prevalent persons whose CAD came before recruitment (historical cases) are removed.
Under either index, every person must also be free of CAD at the T2D diagnosis, because the question is the risk of CAD after T2D.
`prs_cohort()` removes a person whose CAD came at or before the diagnosis, and counts the person under `prior_event`.
Step 10 shows what goes wrong when every person enters at the diagnosis instead.

`status_codes` tells the package what each code of `exit_status` means.
A death ends follow-up, but unlike the end of a study, it also ends the chance of CAD.
Death is therefore a competing event, not censoring.
Step 11 shows what goes wrong when deaths are treated as censoring.

`event_labels` names the two events in the figures and tables.
Without it, they call the events T2 and the competing event.
The labels change no estimate.

The mean and SD used for standardizing are kept, because the persons of another biobank must be standardized with the same values (Step 9):

```{r}
cohort$scores
```

The `sign` column matters when a score was built so that a higher value means lower risk.
Pass `signs = c(outcome = -1)` to flip it.
Choose signs before looking at the outcome, for example from how the score's alleles were aligned.
A sign chosen from the results is a form of overfitting.

### Check the condition of delayed entry

Delayed entry needs a condition.
The times of entry and of censoring must be independent of the times of CAD and death (Geskus 2011).
Part of the condition can be tested.
`check_entry()` tests whether the entry time is quasi-independent of the exit time, the first of CAD and death.
The test is the conditional Kendall's tau of Tsai (1990).

```{r}
entry_test <- check_entry(cohort)
entry_test
```

```{r guard-entry-test, include = FALSE}
# The text below states these results.
stopifnot(abs(entry_test$tau) < 0.05, entry_test$p_value >= 0.05)
```

Tau is near zero, and the p-value is `r sprintf("%.2f", entry_test$p_value)`.
The test gives no evidence that entry and exit are dependent.
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 also marginal.
It can reject when a covariate, such as the age at the T2D diagnosis, changes both the entry and the rate of CAD.
The Cox models of Step 3 need the condition only given their terms, and Step 6 tests the entry in a model with these terms.
The estimates of `cumulative_risk()` need the condition within each score group.

A selection that depends on the future breaks the condition, and delayed entry cannot correct it.
The article [Delayed entry in depth](https://brettklamer.com/work/prsladder/articles/entry-effect.html#a-selection-that-delayed-entry-cannot-correct) shows such a selection.

The test cannot show whether the persons of the biobank are like the population.
For example, UK Biobank participants have lower mortality than the general population of the same age (Fry et al. 2017).
This selection changes the absolute risks, and only data from the population can measure it.

## Step 2: describe the cohort

`summary()` gives the numbers that a report needs: persons, CAD events, deaths, median follow-up, and the distribution of each raw score, overall and by status (0 for no event, 1 for CAD, and 2 for death).

```{r}
summary(cohort)
```

The counts are also split by incident and prevalent persons.
`historical` counts the prevalent persons removed because their CAD came before recruitment.
`prs_cohort()` removed nobody else, because the simulated data have no missing values and the simulator gives CAD only after the T2D diagnosis.
Persons who developed CAD have higher mean scores than persons who did not.
Persons who died have a higher mean T2D score, the score that raises the rate of death.
The next steps measure how strong that difference is and whether it holds after accounting for the other scores.

The PRS reporting standard (PRS-RS) of Wand et al. (2021) asks for the characteristics of each sample, at least age and sex.
`describe_cohort()` gives them as a table of the gtsummary package, with one column per status, and needs that package.
`gtsummary::add_overall()` adds a column for all persons.
The scores are on their raw scale, as they came from the scoring software.

```{r, eval = has_gtsummary}
labels <- list(sex ~ "Sex", birth_year ~ "Year of birth")
year_digits <- birth_year ~
  function(x) gtsummary::style_number(x, big.mark = "")
describe_cohort(cohort, label = labels, digits = year_digits) |>
  gtsummary::add_overall()
```

Without `digits`, the year would print with a thousands separator.
In a terminal, print the table with `gtsummary::as_kable()`.
The table does not hide small cells, which can identify persons in biobank data, and the table object holds the data of every person: export the rendered table, never the R object.

How fast does CAD follow T2D in the whole cohort?
`cumulative_risk()` with `score = NULL` estimates the cumulative incidence of CAD.
This incidence is the proportion of persons with CAD by each time since the T2D diagnosis, when death can come first.
It uses the Aalen-Johansen estimator (Aalen and Johansen 1978), with delayed entry: a prevalent person enters the estimate at recruitment.
`competing_risk` is the cumulative incidence of death before CAD.

```{r}
overall_risk <- cumulative_risk(cohort, score = NULL, times = c(5, 10, 15))
overall_risk
```

`risk_lower` and `risk_upper` (printed as `lower` and `upper`) are the limits of a 95% confidence interval.
`plot_stacked()` shows how the risk divides between the two events.
The lower area is the cumulative incidence of CAD, and the area above it is the cumulative incidence of death.

```{r stacked-cohort, fig.height = 3.5, fig.alt = sprintf("Cumulative incidence of CAD and of death in the whole cohort, stacked. By 15 years, about %.0f%% of persons have CAD and about %.0f%% die first.", 100 * overall_risk$risks$risk[3], 100 * overall_risk$risks$competing_risk[3])}
plot_stacked(overall_risk)
```

## Step 3: fit the ladder of models

With three scores there are eight possible models: no score, each score alone, each pair, and all three.
prsladder calls this set the ladder.
`fit_ladder()` fits two Cox proportional-hazards models for every rung of the ladder.
The model for CAD treats a death as censoring, and the model for death treats CAD as censoring.
These models are the cause-specific hazard models (Putter et al. 2007).
Together they give the cumulative incidence of CAD (Step 7).

Every model also contains the covariates and a smooth (spline) term for age at T2D diagnosis.
Adjusting for principal components guards against population structure, which can create spurious associations (Choi et al. 2020).
Under the recruitment index, keep the birth year as well.
With the age at diagnosis, the birth year sets the calendar year of each person on the clock.
The models can then account for a calendar trend (`?prs_cohort`).
Every model uses the same persons, so any two models can be compared fairly.
`spline_df` names the variables with a spline term and gives the degrees of freedom of each.
The default is 3 for age, and `spline_df = NULL` gives a straight line.
A numeric covariate can also have a spline term, for example with `spline_df = c(age_t1 = 3, birth_year = 3)`.

```{r}
ladder <- fit_ladder(cohort)
ladder
```

```{r values-ladder-table, include = FALSE}
range_ph <- range(ladder$table$ph_p_value)
min_epp <- min(ladder$table$events_per_parameter)
```

Each row is a rung, and the columns describe its model for CAD.
The most useful columns are:

- `concordance`: Harrell's C-index, the probability that of two persons, the one who developed CAD first had the higher predicted risk. 0.5 is a coin toss; 1 is perfect. Because it is measured on the persons of the fit, it is optimistic (Step 8).
- `events_per_parameter`: CAD events per model coefficient. `fit_ladder()` warns below 10, where estimates become unreliable.
- `heuristic_shrinkage`: an estimate of how much the model is overfitted, from its likelihood-ratio statistic (van Houwelingen and le Cessie 1990). A value near 1 means little overfitting, and a value near 0 or below 0 means that the fit is mostly noise.
- `ph_p_value`: the global test of the proportional-hazards assumption, over all terms of the model. A small value means that the hazard ratio of one or more terms can change over follow-up.

Here every model has at least `r floor(min_epp)` CAD events per coefficient.
Every `ph_p_value` is below 0.05, with values from `r sprintf("%.3f", range_ph[1])` to `r sprintf("%.3f", range_ph[2])`.
A global test does not show which term causes the result.
`survival::cox.zph()` gives the test of each term, here for the model with all three scores:

```{r}
ph <- survival::cox.zph(
  ladder$models[["onset+progression+outcome"]],
  transform = "km"
)
round(ph$table, 3)
```

```{r guard-ph-test, include = FALSE}
ph_terms <- ph$table[rownames(ph$table) != "GLOBAL", "p"]
ph_scores <- ph_terms[c("onset", "progression", "outcome")]
# The text below states these results.
stopifnot(
  range_ph[2] < 0.05,
  min(ph_scores) >= 0.05,
  setequal(names(ph_terms)[ph_terms < 0.05], c("PC4", "PC5"))
)
```

The three scores have p-values of `r sprintf("%.2f", min(ph_scores))` or more.
The test gives no evidence that their hazard ratios change over follow-up.
The small global value comes from `PC4` and `PC5`.
All eight models have the same covariates, and their global tests then reject together.
With real data, a small p-value for a score means that its reported hazard ratio is an average over follow-up.
`plot(ph)` then shows how each coefficient changes with time.

Simulation check (not possible with real data): the simulator gives the principal components no effect, so the rejection for `PC4` and `PC5` is a chance result.

`summary()` gives the hazard ratio (HR) for CAD of each score in each model, with its 95% confidence interval:

```{r}
hr <- summary(ladder)
print(
  hr[c("model", "term", "hr", "hr_lower", "hr_upper", "p_value")],
  digits = 3
)
```

```{r guard-hr-coverage, include = FALSE}
hr_progression <- hr$hr[
  hr$model == "onset+progression+outcome" & hr$term == "progression"
]
# The text below states that each interval of this model contains its true
# value.
local({
  truth <- c(onset = 1.10, progression = 1.25, outcome = 1.30)
  full_hr <- hr[hr$model == "onset+progression+outcome", ]
  stopifnot(
    full_hr$hr_lower <= truth[full_hr$term],
    full_hr$hr_upper >= truth[full_hr$term]
  )
})
```

### How to read a hazard ratio per SD

In the model with all three scores, the HR of `progression` is `r sprintf("%.2f", hr_progression)`.
It means that at any time after the T2D diagnosis, the rate of CAD is `r round(100 * (hr_progression - 1))`% higher for a person whose progression score is one SD above another's, when the other terms of the model are the same.
This rate is the rate of CAD among persons who are alive and free of CAD.
A hazard ratio for this rate does not give the same effect on the risk of CAD, because the risk also depends on the rate of death (Putter et al. 2007).
Steps 5 and 7 give the risks.

### Why the hazard ratios change between rows

Alone, each score also carries some of the effect of the other two, because the scores are correlated.
In the model with all three scores, each HR is the effect of that score when the other two are held fixed.

Simulation check (not possible with real data): in the model with all three scores, the hazard ratios are close to the true values of 1.10, 1.25, and 1.30, and each 95% interval contains its true value.

### Hazard ratios for death

`competing_coefficients` holds the hazard ratios for death, in the format of `coefficients`:

```{r}
death_hr <- ladder$competing_coefficients
print(
  death_hr[death_hr$model == "onset+progression+outcome" &
             death_hr$term %in% c("onset", "progression", "outcome"),
           c("term", "hr", "hr_lower", "hr_upper")],
  digits = 3
)
```

```{r guard-death-hr, include = FALSE}
death_full <- death_hr[death_hr$model == "onset+progression+outcome", ]
death_onset <- death_full[death_full$term == "onset", ]
death_null <- death_full[death_full$term %in% c("progression", "outcome"), ]
alone <- hr[hr$model == hr$term, ]
# The text below and the alt text of the figure state these results.
stopifnot(
  death_onset$hr_lower <= 1.2,
  death_onset$hr_upper >= 1.2,
  death_null$hr_lower <= 1,
  death_null$hr_upper >= 1,
  death_onset$hr > death_null$hr,
  alone$hr == tapply(hr$hr, hr$term, max)[alone$term]
)
```

Simulation check (not possible with real data): the simulator gives the T2D score a hazard ratio for death of 1.2, and the estimate is `r sprintf("%.2f", death_onset$hr)` (95% CI `r sprintf("%.2f", death_onset$hr_lower)` to `r sprintf("%.2f", death_onset$hr_upper)`).
The simulator gives the other two scores no effect on death, and the interval of each contains 1.

A score that raises the rate of death lowers the cumulative incidence of CAD, because fewer persons live long enough to develop it.
Latouche et al. (2013) recommend that a report give the cause-specific hazard ratios of both events.
`plot_ladder()` with `event = "both"` draws them, with one panel per score and one row per model:

```{r ladder-both, fig.width = 7, fig.height = 6, fig.alt = sprintf("Cause-specific hazard ratio per SD of each score for CAD and for death, in each model of the ladder. For CAD, the hazard ratio of each score is largest when the score is alone. For death, the T2D score has the largest hazard ratio, about %.1f.", death_onset$hr)}
plot_ladder(ladder, event = "both")
```

## Step 4: does the progression score add information?

`compare_ladder()` compares two nested models, a full model and a reduced model that lacks one or more of its scores.
It uses the likelihood-ratio test, whose null hypothesis is that the added scores carry no information.
The test compares the models for CAD.
It asks whether the scores change the rate of CAD among persons alive and free of CAD.
The test uses all persons, which is one reason not to set persons aside for a test sample.

The plan gives two co-primary tests:

- progression beyond the T2D score (`onset+progression` versus `onset`), and
- progression beyond the CAD score (`progression+outcome` versus `outcome`).

The question asks for information beyond each of the other scores.
The answer is then yes only when both tests reject at the 5% level.
Because both tests must reject, the chance of a false yes is at most 5%, and the level needs no adjustment for two tests.
A plan in which one rejection is enough must test each comparison at the 2.5% level.

```{r}
primary <- list(
  c("onset+progression", "onset"),
  c("progression+outcome", "outcome")
)
primary_tests <- compare_ladder(ladder, pairs = primary)
print(primary_tests, digits = 3)
```

```{r guard-primary-tests, include = FALSE}
# The text below and the answer of the last step state these results.
stopifnot(
  primary_tests$p_value < 0.001,
  primary_tests$delta_aic < 0,
  primary_tests$delta_concordance > 0
)
```

Both p-values are below 0.05, and the rule of the plan gives a yes.
The progression score adds information beyond each of the other scores.
`delta_aic` below zero and `delta_concordance` above zero point the same way; these two columns are descriptive and are not tested.

The secondary comparison of the plan adds progression to a model that already has both other scores (`onset+progression+outcome` versus `onset+outcome`).
Its full model is the prediction model of the report.
The plan allows this test, because both primary tests rejected.

```{r}
secondary <- list(c("onset+progression+outcome", "onset+outcome"))
secondary_test <- compare_ladder(ladder, pairs = secondary)
print(secondary_test, digits = 3)
```

```{r guard-secondary-test, include = FALSE}
# The text below and the answer of the last step state this result.
stopifnot(secondary_test$p_value < 0.001)
```

The p-value of this test is below 0.001.

Without `pairs`, `compare_ladder()` runs all 19 nested comparisons of three scores.
Use them to describe the ladder, not to search for the smallest p-value.
The package does not adjust for multiple testing.

## Step 5: CAD risk by score group

Hazard ratios are hard to explain to a clinician.
The proportion of people who develop CAD within a given number of years is easier.
`cumulative_risk()` divides the persons into groups of equal size by one score (here fifths, or quintiles) and estimates the cumulative incidence of CAD by each time since the T2D diagnosis.

```{r}
by_quintile <- cumulative_risk(
  cohort,
  score = "progression",
  times = c(5, 10),
  groups = 5L
)
print(
  by_quintile$risks[c("group", "time", "n", "n_events", "risk", "difference",
                      "difference_lower", "difference_upper")],
  digits = 2
)
by_quintile$test
```

```{r guard-quintile-risks, include = FALSE}
quintile_5 <- by_quintile$risks$risk[by_quintile$risks$time == 5]
quintile_10 <- by_quintile$risks$risk[by_quintile$risks$time == 10]
# The text below states these results.
stopifnot(
  diff(quintile_10) > 0,
  any(diff(quintile_5) < 0),
  which.min(quintile_5) == 1L,
  which.max(quintile_5) == 5L
)
```

`risk` is the Aalen-Johansen estimate of the proportion with CAD by `time` years after the T2D diagnosis.
At 10 years, the risk rises with the score in every step.
At 5 years, the lowest fifth has the lowest risk and the highest fifth has the highest risk, but the risk does not rise in every step between them.
Each fifth holds only `r by_quintile$risks$n[1]` persons, and few of them have CAD by 5 years.
`difference` is the gap to group 1, the lowest fifth of the progression score; change the comparison group with `reference`.
`test` asks whether the groups differ at all.

`by_quintile` holds the curves as well as the table.
Latouche et al. (2013) recommend that a report give the cumulative incidence of both events.
`plot_risk()` with `event = "both"` draws the two cumulative incidences on one scale, one panel per event, with a 95% band around each curve.
The risk table below counts the persons at risk and the events of each kind.
Each curve has a step at every event of its group.
The `times` of the call set where the curves end and the columns of the risk table, so three to five times are enough.
A prevalent person enters the curve at recruitment, so the number at risk in a group can rise after the diagnosis.
An early risk set can hold few persons, and each event in it adds a large step to every later risk.
Because of this, a curve stops at the first event time with fewer than `min_at_risk` persons at risk, and a group whose first risk sets are too small has no curve.
A positive `entry_lag` in `prs_cohort()` then starts the curves after the small risk sets.

```{r incidence-both, fig.width = 7.5, fig.height = 7, fig.alt = "Cumulative incidence of CAD and of death by quintile of the progression score, with a risk table. The incidence of CAD is higher in the upper quintiles, and the incidence of death shows no trend across the quintiles."}
plot_risk(by_quintile, event = "both")
```

`plot_stacked()` shows the same curves stacked, one panel per group:

```{r stacked-quintiles, fig.width = 7, fig.height = 5, fig.alt = "Cumulative incidence of CAD and of death by quintile of the progression score, stacked. The share of CAD is larger in the upper quintiles."}
plot_stacked(by_quintile)
```

`describe_risk()` gives the risks as a table of the gtsummary package, with one row per quintile and one column per time.
Each cell is the risk with its 95% confidence interval.
It gives one table per event, and `gtsummary::tbl_stack()` joins them:

```{r, eval = has_gtsummary}
gtsummary::tbl_stack(
  list(
    describe_risk(by_quintile),
    describe_risk(by_quintile, event = "competing")
  ),
  group_header = c("CAD", "Death")
)
```

The table holds no data of single persons, but it does not hide small counts either.

With ten groups, curves are hard to tell apart.
`plot_gradient()` shows the risk of each group at one time instead:

```{r}
by_decile <- cumulative_risk(
  cohort,
  score = "progression",
  times = 10,
  groups = 10L
)
```

```{r risk-gradient, fig.alt = sprintf("Ten-year cumulative incidence of CAD by decile of the progression score. The risk is about %.0f%% in the lowest decile and about %.0f%% in the highest decile, with an uneven rise between them.", 100 * by_decile$risks$risk[1], 100 * by_decile$risks$risk[10])}
plot_gradient(by_decile, time = 10)
```

The rise between neighboring deciles is uneven, because each decile holds only `r by_decile$risks$n[1]` persons.

```{r values-case-share, include = FALSE}
top_share <- by_quintile$risks$case_share[
  by_quintile$risks$time == 10 & by_quintile$risks$group == 5
]
```

Keep two points in mind (Choi et al. 2020).
A high score raises risk; it does not make CAD certain.
And most people who develop CAD are not in the top fifth, because the other four fifths together hold four times as many people.
Here the top fifth has `r round(100 * top_share)`% of the CAD cases by ten years (`case_share` in `by_quintile$risks`).

The article [Risk by score group in depth](https://brettklamer.com/work/prsladder/articles/score-groups.html) gives three more results of `cumulative_risk()`:

- the share of the CAD cases in each group, with its confidence interval,
- the top 10% of the score against the other persons, with `cut_probs`, and
- the years that each group loses to CAD, the restricted mean time lost.

## Step 6: choose a design for the absolute risks

The models of Step 3 have one clock for all persons, the years since the T2D diagnosis.
At year 5, a prevalent person and an incident person then share one baseline hazard.
This assumption can fail.
A prevalent person was well enough to join the biobank years after the diagnosis, and lived with T2D in an earlier calendar period.
An incident person is followed through the first months after the diagnosis, when the work-up can find CAD.
Such differences matter most for the absolute risks, because the baseline hazard sets those risks.

`fit_ladder()` has three designs:

| `design` | Clock | Hazard after entry | Use |
|---|---|---|---|
| `"t1"` (default) | Years since the diagnosis | One baseline for all persons | Hazard ratios and primary tests |
| `"t1_entry"` | Years since the diagnosis | One baseline, a constant factor for prevalent persons, and a log-linear term for the years from the diagnosis to entry | Absolute risks, when the rule picks it |
| `"index"` | Years since entry | One baseline for incident persons and one for prevalent persons, and a term for the years from the diagnosis to entry | Absolute risks, when the rule picks it |

`?fit_ladder` gives the models and the assumption of each design.
`check_design()` tests the entry in the model with all scores, and it applies a rule that the analysis plan fixed:

```{r}
checks <- check_design(cohort)
checks
```

```{r helper-test-p-value, include = FALSE}
# The p-value of one test of check_design(), selected by its name. A name that
# matches no row stops the build. Without the check, the value has length 0,
# which passes stopifnot() and prints nothing in the text.
test_p_value <- function(checks, test) {
  p_value <- checks$tests$p_value[checks$tests$test == test]
  stopifnot(length(p_value) == 1L)
  p_value
}
```

```{r guard-design-pick, include = FALSE}
# The text below states that the rule picks "t1_entry". Step 7 names this
# design.
stopifnot(
  test_p_value(checks, "entry_bands") >= 0.05,
  checks$design == "t1_entry"
)
```

The three tests ask three questions:

- `entry_terms`: do the two terms of `"t1_entry"` improve on `"t1"`? A rejection shows that `"t1"` is not correct, but not that `"t1_entry"` is adequate.
- `entry_bands`: does the factor of the prevalent persons change with the years since entry? The default bands end at 1, 2, 4, and 7 years. Fix the bands before the analysis, because the p-value assumes that the data did not choose them.
- `index_duration`: does the effect of the years from the diagnosis to entry change with the years since entry? This effect must stay constant under `"index"`.

The rule uses only `entry_bands`.
If the test does not reject at the 5% level, the pick is `"t1_entry"`.
If it rejects, the pick is `"index"`.
The rule cannot pick `"t1"`.
Here `entry_bands` has a p-value of `r sprintf("%.2f", test_p_value(checks, "entry_bands"))`, and the pick is ``r paste0('"', checks$design, '"')``.
With `r sum(cohort$data$status)` CAD events, a result without a rejection is weak evidence.

Simulation check (not possible with real data): the pick is the expected one, because `simulated` has no effect of the entry.
`simulate_prs_data()` adds such an effect only when its arguments `prevalent_hazard_ratio` or `duration_hazard_ratio` differ from their defaults.

The last step of the rule is a new ladder with the pick:

```{r}
risk_ladder <- fit_ladder(cohort, design = checks$design)
hr_risk <- summary(risk_ladder)
in_full <- function(x) x[x$model == "onset+progression+outcome", ]
data.frame(
  term = in_full(hr)$term,
  t1 = round(in_full(hr)$hr, 3),
  pick = round(in_full(hr_risk)$hr, 3)
)
```

```{r values-hr-change, include = FALSE}
hr_change <- max(abs(in_full(hr_risk)$hr - in_full(hr)$hr))
```

The hazard ratios of the two ladders differ by at most `r sprintf("%.3f", hr_change)` here.
The report keeps `ladder` for the hazard ratios and the tests of Steps 3 and 4.
The data then do not choose the design of a primary test.
The absolute risks of Steps 7 to 9 use `risk_ladder`.
Each result holds its design in `settings$design`.
`loglik`, `aic`, and `concordance` do not compare designs (`?fit_ladder`).

`?check_design` gives the limits of the rule.
The pick comes from the CAD model, and the death model uses the same design.
The rule has not been studied with competing events.
A study that pools several biobanks fixes one design for all of them, because the pick varies by chance.

The article [Delayed entry in depth](https://brettklamer.com/work/prsladder/articles/entry-effect.html#a-biobank-with-an-effect-of-the-entry) shows the rule on data in which recruitment changes the hazard of CAD.
In those data, `entry_bands` rejects, and the pick is `"index"`.

## Step 7: predict absolute risk

### Which horizons can be evaluated?

A prediction horizon counts years after entry: after the diagnosis for an incident person, and after recruitment for a prevalent person.
A horizon can only be evaluated if enough persons are followed that long after entry.
`horizon_support()` checks each horizon for at least 10 events, 50 persons still followed, and enough uncensored follow-up.

```{r}
horizon_support(cohort, horizons = c(5, 10, 15, 20))
```

```{r guard-horizon-support, include = FALSE}
# The text below states these results.
local({
  support <- horizon_support(cohort, horizons = c(5, 10, 15, 20))
  stopifnot(
    identical(support$supported, c(TRUE, TRUE, FALSE, FALSE)),
    support$censoring_survival[support$horizon == 15] < 0.02,
    support$n_at_risk[support$horizon == 20] < 50
  )
})
```

The guard uses the smaller censoring survival of the incident and the prevalent persons (see `?horizon_support`).
At 15 years, the censoring survival of the incident persons is below 0.02, so that horizon fails.
Too few persons are followed for 20 years after entry, so that horizon also fails.
Every `assess_*()` function runs this check and skips an unsupported horizon with a warning.
The planned horizons, 5 and 10 years, are supported.

```{r}
horizons <- c(5, 10)
```

### Absolute risk for a person

`predict_risk()` turns a rung into each person's probability of CAD within a horizon after entry, here ten years.
Because each rung holds its model for death, the probability is the cumulative incidence of CAD.
For a prevalent person, it is the risk in the ten years after recruitment, given the years already lived with T2D.

```{r}
report_fit <- risk_ladder$models[["onset+progression+outcome"]]
round(predict_risk(report_fit, newdata = cohort$data[1:5, ], horizon = 10), 3)
```

The persons do not need to be in the data.
Suppose a 55-year-old woman received a T2D diagnosis in 2010 and we know her raw scores.
Her diagnosis year and birth year are within the range of the biobank, because the model uses `birth_year` and cannot be trusted outside its data.
We compare two versions of her, with the progression score one SD above and one SD below the average.
The principal components are set to 0, the average of the simulated ones.
We use the column names of the development data, and `prediction_data()` applies the rules of the cohort: it standardizes each raw score with the development mean, SD, and sign, as `prs_cohort()` did, and sets the entry.

```{r}
persons <- data.frame(
  age_at_diagnosis = 55,
  sex = "Female",
  PC1 = 0, PC2 = 0, PC3 = 0, PC4 = 0, PC5 = 0,
  birth_year = 1955,
  t2d_score = 0.2,
  t2d_to_cad_score = c(-1, 1),
  cad_score = 0.5
)
newdata <- prediction_data(persons, cohort)
newdata[c("start", "age_t1", "t1_status", "t1_duration", "progression")]
woman <- rbind(
  cumulative_incidence = predict_risk(report_fit, newdata, horizon = 10),
  net_risk = predict_risk(report_fit, newdata, horizon = 10, competing = FALSE),
  death = predict_risk(report_fit, newdata, horizon = 10, event = "competing")
)
round(woman, 3)
```

`start = 0` means that she enters follow-up at her diagnosis, and `t1_status` and `t1_duration` show that she counts as an incident person.
Her ten-year risk of CAD is about `r round(100 * woman["cumulative_incidence", 1])`% with a low progression score and `r round(100 * woman["cumulative_incidence", 2])`% with a high one.
This cumulative incidence is the risk to tell a patient.
`competing = FALSE` gives the net risk, the risk in a world where nobody dies before CAD.
The net risk is higher, because it ignores that she can die first.
The last row, from `event = "competing"`, is her risk of death before CAD.
It uses the same two models, with their roles swapped.

Three years later she is 58 and still alive and free of CAD.
Her risk over the next ten years is now a conditional risk: the risk from year 3 to year 13 after her diagnosis, given no CAD and no death by year 3 (Benichou and Gail 1990).
`age_at_entry` names the column with the age from which to predict.

```{r}
persons$age_now <- 58
later <- prediction_data(persons, cohort, age_at_entry = "age_now")
later[c("start", "t1_status", "t1_duration")]
```

`prediction_data()` now calls her prevalent, like a volunteer who joined the biobank three years after the diagnosis.
Under `risk_ladder`, she would then get the factor of those volunteers.
She is a patient of a clinic, and the model cannot show whether she resembles them.
This choice belongs to the analyst (`?prediction_data`).
`t1_status = "incident"` treats her as an incident person who was followed from her diagnosis to year 3:

```{r}
later <- prediction_data(
  persons,
  cohort,
  age_at_entry = "age_now",
  t1_status = "incident"
)
later[c("start", "t1_status", "t1_duration")]
woman_later <- predict_risk(report_fit, later, horizon = 10)
round(woman_later, 3)
```

```{r guard-conditional-risk, include = FALSE}
# The text above and below states these results.
stopifnot(
  woman["net_risk", ] > woman["cumulative_incidence", ],
  woman_later > woman["cumulative_incidence", ]
)
```

`start` stays at 3, the time of the prediction.
`t1_status` and `t1_duration` describe her entry.

Her risk over the next ten years is now about `r round(100 * woman_later[1])`% and `r round(100 * woman_later[2])`%.
It is a little higher than at diagnosis, because in this simulation the rate of CAD increases with time since the T2D diagnosis.


### Risk across the score distribution

A hazard ratio per SD does not give the size of the risk.
`standardized_risk()` gives the absolute risk at percentiles of one score.
For each percentile, it gives every person of the cohort the score at that percentile and keeps the other predictors of each person.
It then averages the predicted risks of the persons (Makuch 1982).

```{r}
by_percentile <- standardized_risk(
  risk_ladder,
  "progression",
  times = 10,
  probs = c(0.1, 0.5, 0.9)
)
print(
  by_percentile$risks[c("prob", "raw_score", "time", "risk", "competing_risk",
                        "difference", "ratio")],
  digits = 2
)
```

```{r guard-percentile-risks, include = FALSE}
median_t1 <- standardized_risk(
  ladder,
  "progression",
  times = 10,
  probs = 0.5
)$risks$risk
percentile_risk <- by_percentile$risks$risk
percentile_death <- by_percentile$risks$competing_risk
hr_alone <- hr$hr[hr$model == "progression" & hr$term == "progression"]
# The text below states that this risk falls across the percentiles.
stopifnot(diff(percentile_death) < 0)
```

At the 10th percentile of the progression score, the ten-year cumulative incidence of CAD is `r round(100 * percentile_risk[1])`%.
At the 90th percentile, it is `r round(100 * percentile_risk[3])`%.
`difference` and `ratio` compare each percentile with the first, here the 10th percentile.
`reference` changes the comparison.
In Step 5, the lowest and the highest fifth have ten-year risks of `r round(100 * quintile_10[1])`% and `r round(100 * quintile_10[5])`%.
These two groups hold the persons with scores near these percentiles, but the results differ, because the scores are correlated.
A group also holds the other scores of its persons, and a standardized risk keeps the other scores of each person of the cohort.
As in `cumulative_risk()`, `times` counts years since the T2D diagnosis.
The average is over the persons of the recruited cohort.
Under `risk_ladder`, each person counts as an incident person who starts at the diagnosis.
The result is then the risk from the diagnosis with the hazard of an incident person (`?standardized_risk`).

For most prevalent persons of this cohort, this risk is an extrapolation.
An incident person has the diagnosis after recruitment starts, and most prevalent persons here have the diagnosis before the first diagnosis of an incident person.
The standardized risk then extrapolates to incident persons diagnosed before that time.
At the median, `risk_ladder` gives `r round(100 * percentile_risk[2])`% and `ladder` gives `r round(100 * median_t1)`%.
Under the design `"t1"` of `ladder`, a risk from the diagnosis does not extrapolate.
The article [Delayed entry in depth](https://brettklamer.com/work/prsladder/articles/entry-effect.html#the-cost-of-the-pick-for-a-risk-from-the-diagnosis) compares the predictions of both ladders with the true risks of the simulator, and it gives the cost of the extrapolation.
This cost is not a reason to use `"t1"` for absolute risks, because an effect of the entry makes `"t1"` biased (`?fit_ladder`).
Give a risk from the diagnosis with an interval, such as the bootstrap interval below.
The risks from entry of Steps 8 and 9 do not extrapolate in this way.

`competing_risk` is the standardized cumulative incidence of death before CAD.
It falls across the percentiles, from `r sprintf("%.1f", 100 * percentile_death[1])`% to `r sprintf("%.1f", 100 * percentile_death[3])`%.
A higher rate of CAD lowers this incidence, because a person with CAD is no longer at risk of death before CAD.
The score is in the models of both events.
Because of this, the contrast between two percentiles also includes any effect of the score on death.

The default model, `model = "progression"`, has the progression score as its only score.
It adjusts for the covariates and the age at diagnosis, but not for the other two scores.
Its HR per SD for CAD is `r sprintf("%.2f", hr_alone)` (Step 3).
A larger model also adjusts for the other scores, which keep their observed values.
`n_boot` adds bootstrap percentile intervals.
As in Step 8, this tutorial uses 30 replicates to run fast, and an analysis needs 200 or more.

```{r boot-standardized-risk}
adjusted <- standardized_risk(
  risk_ladder,
  "progression",
  times = 10,
  model = "onset+progression+outcome",
  n_boot = 30L,
  seed = 2026L
)
print(
  adjusted$risks[c("prob", "risk", "risk_lower", "risk_upper", "difference",
                   "difference_lower", "difference_upper")],
  digits = 2
)
```

```{r guard-interval-width, include = FALSE}
adjusted_t1 <- standardized_risk(
  ladder,
  "progression",
  times = 10,
  model = "onset+progression+outcome",
  n_boot = 30L,
  seed = 2026L
)
top_pick <- adjusted$risks[adjusted$risks$prob == 0.9, ]
top_t1 <- adjusted_t1$risks[adjusted_t1$risks$prob == 0.9, ]
# The text below states that the interval of `ladder` is narrower.
stopifnot(
  nrow(top_pick) == 1L,
  nrow(top_t1) == 1L,
  top_t1$risk_upper - top_t1$risk_lower <
    top_pick$risk_upper - top_pick$risk_lower
)
```

```{r guard-adjusted-contrast, include = FALSE}
# The text below states that the contrast is smaller.
stopifnot(
  adjusted$risks$difference[adjusted$risks$prob == 0.9] <
    by_percentile$risks$difference[by_percentile$risks$prob == 0.9]
)
```

The contrast is smaller, because in this model the HR of the progression score is `r sprintf("%.2f", hr_progression)`.
Each replicate fits the model again, and the intervals include the variance of the extrapolation (`?standardized_risk`).
At the 90th percentile, the interval of the risk is `r round(100 * top_pick$risk_lower)`% to `r round(100 * top_pick$risk_upper)`%.
The same call with `ladder` gives `r round(100 * top_t1$risk_lower)`% to `r round(100 * top_t1$risk_upper)`%.
This interval is narrower, because a risk from the diagnosis does not extrapolate under the design `"t1"`.

## Step 8: internal validation

### Apparent performance

`assess_discrimination()` reports two measures at each horizon:

- the time-dependent AUC, the probability that a person with CAD by the horizon had a higher predicted risk than a person without CAD at the horizon (Blanche et al. 2013), and
- the Brier score, the mean squared difference between predicted risk and outcome (Gerds and Schumacher 2006). Lower is better.

A person who died before CAD and before the horizon counts as a person without CAD (Blanche et al. 2013).
Both measures weight each person by the inverse probability of remaining uncensored.
Without other arguments, the function evaluates the development persons:

```{r}
apparent <- assess_discrimination(risk_ladder, horizons, pairs = primary)
print(
  apparent$models[apparent$models$horizon == 10,
                  c("model", "evaluation", "auc", "auc_lower", "auc_upper",
                    "brier")],
  digits = 3
)
```

The column `evaluation` says `"apparent"`.
These values are optimistic, because the models were fitted to the same persons.
Their intervals also treat the fitted models as fixed.
Do not report them as the performance of the models on new persons.

### One plan for all functions

The bootstrap corrects the apparent values for optimism (Harrell et al. 1996).
Each replicate draws a sample of the persons with replacement and fits the ladder again on it.
The replicate then measures how much better the new fit does on its own sample than on all persons.
The mean of this difference is the optimism, and the corrected value is the apparent value minus the optimism.
The correction uses all persons, and Steyerberg et al. (2001) found it more accurate than a random split.

`bootstrap_plan()` holds the settings.
We give the same plan to every `assess_*()` function, and each function then uses the same bootstrap samples.
This tutorial uses 30 replicates to run fast.
An analysis needs 200 or more.
On one computer, each of the three bootstrap calls of this step took 13 to 19 seconds with 30 replicates.
With 200 replicates, the three calls then take about 5 minutes for this cohort.
`?bootstrap_plan` gives the times for a larger cohort.

```{r}
plan <- bootstrap_plan(n_boot = 30L, seed = 2026L)
```

The bootstrap plan holds only its settings, and no data.
`interval = "analytic"`, the default, gives an interval to the measures of one model that have a standard error from an influence function, such as the AUC.
`interval = "two-stage"` gives an interval to every measure, at the cost of `n_boot * (n_inner + 1)` replicates.

The correction repeats only the steps of `fit_ladder()`.
The analysis plan fixed the pairs, thresholds, rung, and `spline_df` before the analysis.
A choice made from the results makes the correction too small.
The correction takes the design as fixed, and it does not repeat the selection of `check_design()`.

### Discrimination

```{r boot-discrimination}
discrimination <- assess_discrimination(
  risk_ladder,
  horizons,
  pairs = primary,
  validation = plan
)
both <- rbind(apparent$models, discrimination$models)
print(
  both[both$horizon == 10 & both$model == "onset+progression+outcome",
       c("model", "evaluation", "auc", "auc_lower", "auc_upper", "brier")],
  digits = 3,
  row.names = FALSE
)
```

The column `evaluation` of the corrected result says `"optimism-corrected"`.
`optimism` holds the apparent value, the optimism, and its Monte Carlo standard error, the part of the uncertainty that comes from the number of replicates:

```{r}
optimism_auc <- discrimination$optimism
optimism_auc <- optimism_auc[
  optimism_auc$measure == "auc" & optimism_auc$horizon == 10,
]
print(
  optimism_auc[c("model", "apparent", "optimism", "optimism_mcse",
                 "n_boot_used")],
  digits = 3
)
```

```{r guard-optimism, include = FALSE}
optimism_full <- optimism_auc$optimism[
  optimism_auc$model == "onset+progression+outcome"
]
auc_full_se <- apparent$models$auc_se[
  apparent$models$model == "onset+progression+outcome" &
    apparent$models$horizon == 10
]
auc_full <- discrimination$models$auc[
  discrimination$models$model == "onset+progression+outcome" &
    discrimination$models$horizon == 10
]
mcse_ratio <- max(optimism_auc$optimism_mcse / optimism_auc$optimism)
# The text below states that the eight models have about the same optimism.
stopifnot(diff(range(optimism_auc$optimism)) < 0.005)
```

With `r sum(cohort$data$status)` CAD events, the optimism of the model with all three scores is `r sprintf("%.3f", optimism_full)`, against a standard error of its AUC of `r sprintf("%.3f", auc_full_se)`.
The eight models have about the same optimism, from `r sprintf("%.3f", min(optimism_auc$optimism))` to `r sprintf("%.3f", max(optimism_auc$optimism))`.
Each model has the covariates and the age spline, and a score adds only one coefficient.
In a smaller cohort, or with more coefficients, the optimism grows.
With 30 replicates, the Monte Carlo standard error is at most `r round(100 * mcse_ratio)`% of the optimism.
With fewer replicates, it can be as large as the optimism, which is one reason why an analysis needs more replicates.

`plot_discrimination()` draws the corrected values, with one panel per horizon:

```{r discrimination, fig.width = 7, fig.alt = sprintf("Optimism-corrected time-dependent AUC of each model at 5 and 10 years. At 10 years the AUC is about %.2f for the base model and about %.2f for the model with all three scores.", optimism_auc$apparent[1] - optimism_auc$optimism[1], auc_full)}
plot_discrimination(discrimination)
```

```{r guard-base-auc, include = FALSE}
auc_base <- discrimination$models[discrimination$models$model == "base" &
                                    discrimination$models$horizon == 10, ]
# The text below states that the scores raise the AUC.
stopifnot(nrow(auc_base) == 1L, auc_base$auc > 0.5, auc_full > auc_base$auc)
```

The `base` model has no score.
It predicts from the covariates and the age at diagnosis.
At 10 years, its corrected AUC is `r sprintf("%.2f", auc_base$auc)`, against `r sprintf("%.2f", auc_full)` for the model with all three scores.

The `pairs` element gives the change in AUC and Brier score when progression is added:

```{r}
print(
  discrimination$pairs[c("full", "reduced", "horizon", "delta_auc",
                         "delta_auc_lower", "delta_auc_upper", "delta_brier")],
  digits = 3
)
```

Each rung has its own model for death, with the same scores as its model for CAD.
A pair of rungs therefore adds the score to both models.
The corrected changes have no confidence interval.
The interval of the apparent change treats the fitted models as fixed.
With an added score of no effect, it covers the true change too rarely (`?assess_discrimination`).
The likelihood-ratio tests of Step 4 answer whether the progression score adds information.
No improvement in performance is the same null hypothesis as a coefficient of zero for the added score (Pepe et al. 2013).
Report `delta_auc` as a size of effect, and use the likelihood-ratio test as the test of added information.
An AUC of about `r sprintf("%.2f", auc_full)` means that the scores separate groups of people, not individuals, which is typical for a PRS.

### Calibration: are the risks right?

A model can rank persons well and still predict risks that are too high or too low.
`assess_calibration()` compares predicted and observed risk.

```{r boot-calibration}
calibration <- assess_calibration(
  risk_ladder,
  horizons,
  groups = 10L,
  validation = plan
)
overall <- calibration$models
print(
  overall[overall$model == "onset+progression+outcome",
          c("horizon", "oe_ratio", "expected_risk", "observed_risk", "slope")],
  digits = 3
)
```

- `oe_ratio` is observed divided by expected CAD events (Crowson et al. 2016). 1 is ideal; above 1 means that the model predicts too few events. It uses the model for CAD only.
- `expected_risk` is the mean predicted cumulative incidence, and `observed_risk` is the Aalen-Johansen estimate. Their comparison checks the cumulative incidence, which uses the models of both events (Wolbers et al. 2009).
- `slope` is the calibration slope of the model for CAD. 1 is ideal; below 1 means that predictions are too extreme, the usual sign of overfitting.

On the development persons, the apparent O/E ratio is close to 1 by construction, and the apparent slope over all follow-up is exactly 1.
Only the corrected values are informative here.
The table holds the corrected O/E ratio and the corrected slope.
The columns `expected_risk` and `observed_risk` keep their apparent values (`?assess_calibration`).

```{r guard-calibration, include = FALSE}
cal_full <- overall[overall$model == "onset+progression+outcome", ]
# The text below states these results, in the order of the horizons.
stopifnot(
  identical(as.numeric(cal_full$horizon), c(5, 10)),
  abs(log(cal_full$oe_ratio)) < 0.05,
  cal_full$slope[1] > 1,
  cal_full$slope[2] < 1
)
```

The corrected O/E ratio is `r sprintf("%.2f", cal_full$oe_ratio[1])` at 5 years and `r sprintf("%.2f", cal_full$oe_ratio[2])` at 10 years.
For new persons of the same population, the model for CAD then predicts about the observed number of events.
The corrected slope is `r sprintf("%.2f", cal_full$slope[1])` at 5 years and `r sprintf("%.2f", cal_full$slope[2])` at 10 years, and the two horizons do not agree.

Because the slope at a horizon uses only the events before the horizon, it is less precise.
Its apparent value is also not exactly 1, because the windows stop at the horizon.
The corrected slope over all follow-up, in the element `shrinkage`, is the better measure of overfitting:

```{r}
print(calibration$shrinkage[c("model", "slope")], digits = 3)
```

```{r values-shrinkage, include = FALSE}
shrink <- calibration$shrinkage
shrink_full <- shrink$slope[shrink$model == "onset+progression+outcome"]
heuristic_full <- risk_ladder$table$heuristic_shrinkage[
  risk_ladder$table$model == "onset+progression+outcome"
]
```

For the model with all three scores, the corrected slope is `r sprintf("%.3f", shrink_full)`.
New persons would need their linear predictors shrunk by this factor toward the mean (van Houwelingen and le Cessie 1990).
For the same model of `risk_ladder`, `heuristic_shrinkage` gave `r sprintf("%.3f", heuristic_full)` without a bootstrap.
The table of Step 3 shows the value of `ladder`, which has another design.

The `groups` element has the predicted and the observed cumulative incidence in each tenth of predicted risk, with apparent values.
`plot_calibration()` draws them:

```{r calibration, fig.alt = "Observed against predicted ten-year cumulative incidence of CAD in tenths of predicted risk on the development persons. The points lie close to the diagonal."}
plot_calibration(calibration, horizon = 10, model = "onset+progression+outcome")
```

The axes give the risk within 10 years of entry, the time scale of the predictions.
Without competing events, `assess_calibration()` also gives a smooth calibration curve, with its summaries `ici`, `e50`, and `e90` (Austin et al. 2020).
With competing events, the package does not estimate this curve yet, and these columns are `NA`.

### Decision curves

Clinical guidelines act on risk categories.
The 2023 ESC guidelines for people with type 2 diabetes define categories of 10-year cardiovascular disease risk (Marx et al. 2023).
A risk of 10% to less than 20% is high, and a risk of 20% or more is very high.
The guidelines estimate this risk with SCORE2-Diabetes (SCORE2-Diabetes Working Group 2023).
This model predicts a composite of cardiovascular events, and it treats death from other causes as a competing event.
Our risk is the cumulative incidence of CAD only.
We use the categories only to show the functions.

A model flags a person whose predicted risk is at or above a threshold.
A flag helps a person who will develop CAD.
It can harm a person who will not, for example through a treatment that the person does not need.
Net benefit weighs the two (Vickers and Elkin 2006).
A threshold $p_t$ sets the weight.
At $p_t$, the harm of one false positive is $p_t / (1 - p_t)$ times the benefit of one true positive.
Net benefit is the number of true positives per person, minus the number of false positives per person times $p_t / (1 - p_t)$.
With censored data, `assess_net_benefit()` weights each person by the inverse probability of remaining uncensored (Blanche et al. 2013).
A model must do better than two simple strategies.
Treat all flags every person.
Treat none flags no person and has a net benefit of 0.

A decision curve shows net benefit across a range of thresholds.
Choose the range from the clinical decision, as the first step of the analysis (Vickers et al. 2016).
Do not choose a threshold from the decision curve (Kerr et al. 2016).
We use 5% to 30%, around the ESC categories at 10% and 20%.
The horizon is 10 years after entry.

```{r boot-net-benefit}
net_benefit <- assess_net_benefit(
  risk_ladder,
  horizons = 10,
  thresholds = (5:30) / 100,
  pairs = c(primary, secondary),
  validation = plan
)
```

`plot_net_benefit()` draws the corrected curves of the report model and of the same model without the progression score, with 95% confidence bands:

```{r net-benefit, fig.alt = "Optimism-corrected net benefit of the model with all three scores and of the model without the progression score, of treat all, and of treat none, for thresholds of ten-year CAD risk from 5% to 30%. Below about 9%, the two models and treat all are close. From about 9% to about 22%, the model with all three scores is above treat all and treat none. Above about 22%, both models are close to treat none. The bands of the two models overlap."}
report_pair <- c("onset+outcome", "onset+progression+outcome")
plot_net_benefit(net_benefit, horizon = 10, model = report_pair)
```

```{r guard-decision-curve, include = FALSE}
# The text below and the alt text above state these results. The limits lie
# between two thresholds of the grid, so no comparison needs an exact double.
local({
  treat_all <- net_benefit$treat_all[net_benefit$treat_all$horizon == 10, ]
  curves <- net_benefit$models[net_benefit$models$horizon == 10, ]
  reduced <- curves[curves$model == "onset+outcome", ]
  reduced <- reduced$net_benefit[match(treat_all$threshold, reduced$threshold)]
  full <- curves[curves$model == "onset+progression+outcome", ]
  full <- full$net_benefit[match(treat_all$threshold, full$threshold)]
  low <- treat_all$threshold < 0.085
  middle <- treat_all$threshold > 0.085 & treat_all$threshold < 0.225
  high <- treat_all$threshold > 0.225
  stopifnot(
    !anyNA(reduced),
    !anyNA(full),
    abs(reduced[low] - treat_all$net_benefit[low]) < 0.005,
    abs(full[low] - treat_all$net_benefit[low]) < 0.005,
    full[middle] > pmax(0, treat_all$net_benefit[middle]),
    mean(full[middle] >= reduced[middle]) > 0.5,
    abs(reduced[high]) < 0.005,
    abs(full[high]) < 0.005
  )
})
```

Below about 9%, the two models and treat all have almost the same net benefit, because the models flag most persons.
From about 9% to about 22%, the report model has a higher net benefit than treat all and treat none at every threshold of the grid.
In this range, it is also higher than the model without the progression score at most thresholds.
Above about 22%, both models are close to treat none.
Each band is the interval of one model.
Because both models use the same persons, their net benefits are correlated.
Because of this, overlapping bands do not show that the models are alike, and the plot cannot compare them.
The change in net benefit compares them directly.
Treat all uses no model, and it has no optimism.
Its curve has the apparent values of the development persons.

`pairs` gives the change in net benefit from the reduced model to the full model of each pair:

```{r}
nb_pairs <- net_benefit$pairs
print(
  nb_pairs[nb_pairs$threshold %in% c(0.1, 0.2),
           c("full", "reduced", "threshold", "delta_net_benefit")],
  digits = 2
)
```

```{r guard-net-benefit-change, include = FALSE}
nb_gain <- nb_pairs$delta_net_benefit[
  nb_pairs$full == "onset+progression+outcome" & nb_pairs$threshold == 0.2
]
nb_gain_10 <- nb_pairs$delta_net_benefit[
  nb_pairs$full == "onset+progression+outcome" & nb_pairs$threshold == 0.1
]
# The text below states these results.
stopifnot(
  length(nb_gain) == 1L,
  length(nb_gain_10) == 1L,
  nb_gain > 0,
  abs(nb_gain_10) < 0.001
)
```

The last two rows are the secondary comparison.
At 20%, adding the progression score to the other two scores gives `r sprintf("%.4f", nb_gain)` more net benefit.
This gain has the value of `r sprintf("%.1f", 1000 * nb_gain)` more true positives per 1000 persons, with no more false positives.
At 10%, the change is `r sprintf("%.4f", nb_gain_10)`, less than 1 true positive per 1000 persons.

These changes are sizes of effect without a confidence interval, as for `delta_auc` (`?assess_net_benefit`).
The likelihood-ratio tests of Step 4 answer whether the progression score adds information.
The decision curve answers how much that information changes decisions.
An interval for `delta_net_benefit` needs a bootstrap plan with `interval = "two-stage"` and at least 1000 outer replicates.
That plan takes several hours for this cohort, and the tutorial does not run it (`?bootstrap_plan`).

The article [Validation in depth](https://brettklamer.com/work/prsladder/articles/validation-depth.html) continues this step with three more results:

- the persons that each model flags at 10% and 20%, with the sensitivity and the predictive values of `assess_thresholds()`,
- the report model against treat all, treat none, and the model without the progression score, and
- the stability of the predicted risk of one person, from `assess_stability()`.


## Step 9: external validation

Internal validation estimates the performance in new persons from the same population.
A different cohort, such as another biobank, a later period, or another ancestry group, can differ in its baseline risk, its case mix, and its measurements.
External validation measures whether the models transport to such a cohort (Debray et al. 2015).

Our second biobank is also simulated, with the same rate of death and the same effects of age and sex.
It differs from the first in two ways, as a biobank of another population can.
Its rate of CAD is much higher.
If nobody died, a man with a T2D diagnosis at age 55 and average scores would have a 40% risk of CAD within ten years, against 15% in the first biobank.
And its scores have weaker effects, with log hazard ratios about half as large, as a score often has in an ancestry that differs from that of its GWAS (Choi et al. 2020).
`id_prefix` gives its persons other identifiers, and we keep its recruited persons, as for the first biobank.

```{r}
second_biobank <- simulate_prs_data(
  n = 4000L,
  seed = 2027L,
  hazard_ratios = c(t2d = 1.05, t2d_to_cad = 1.12, cad = 1.15),
  score_correlation = 0.3,
  baseline_risk = 0.40,
  incident_probability = 0.3,
  competing_risk = 0.2,
  competing_hazard_ratios = c(t2d = 1.2),
  covariate_hazard_ratios = c(age_at_diagnosis = 1.04, sex = 0.7),
  competing_covariate_hazard_ratios = c(age_at_diagnosis = 1.08, sex = 0.7),
  id_prefix = "E"
)
second_biobank <- second_biobank[second_biobank$recruited, ]
```

### The export

A ladder holds the data of every person, and it cannot leave the first biobank.
`export_ladder()` keeps only what the predictions need:

- the coefficients and the knots of the age spline,
- a baseline hazard with steps of at least 10 events, for CAD and for death,
- the follow-up limits, and
- the rules of the cohort.

```{r}
export <- export_ladder(risk_ladder)
names(export)
```

`models` holds one model for each rung, `rules` the rules of the cohort, `development` the summaries for `describe_case_mix()`, and `settings` the settings of the export.
`print(export)` lists each value that leaves the biobank, for a disclosure review.
The list has `r length(capture.output(print(export)))` lines for this ladder, and the tutorial does not show it.
Some values can come from one or two persons, such as a knot of the age spline.
`?export_ladder` gives the disclosure status of each element, for the review of each biobank.

The first biobank saves the export to a file, and the second biobank loads it:

```{r}
path <- tempfile(fileext = ".rds")
saveRDS(export, path)
imported <- readRDS(path)
```

### The external cohort

`prs_cohort()` with `reference = imported` builds the cohort with the rules of the first biobank: the same columns, status codes, settings, and exclusions, and each score standardized with the development mean, SD, and sign.

```{r}
external <- prs_cohort(second_biobank, reference = imported)
external
```

With a second cohort, `describe_cohort()` gives one column per cohort:

```{r, eval = has_gtsummary}
describe_cohort(cohort, external, label = labels, digits = year_digits)
```

The table needs the data of both cohorts.
When the data cannot leave a biobank, each biobank makes its own column.

`describe_case_mix()` compares the linear predictor of each model in the two cohorts.
The development values come from the export:

```{r}
print(describe_case_mix(imported, external), digits = 3)
```

```{r guard-case-mix, include = FALSE}
# The text below states that the two cohorts have almost the same case mix.
local({
  case_mix <- describe_case_mix(imported, external)
  stopifnot(
    abs(case_mix$lp_mean_difference) < 0.05,
    abs(case_mix$lp_sd_ratio - 1) < 0.05
  )
})
```

The two cohorts have almost the same mean and spread of the linear predictor, because the simulator draws the scores and the covariates of both biobanks from the same distributions.
A larger spread in the external cohort usually gives a higher AUC, and a different mean can shift the O/E ratio (Debray et al. 2015).
Here neither can explain a difference.
A loss of performance must then come from a different baseline risk or from different effects of the predictors.

### Performance in the second biobank

Collins et al. (2016) recommend at least 100 events, and preferably 200, for an external validation.
The second biobank has `r sum(external$data$status)` CAD events.

With `newdata`, the `assess_*()` functions evaluate the external cohort, and the column `evaluation` says `"external"`.

```{r}
external_calibration <- assess_calibration(
  imported,
  horizons,
  groups = 10L,
  newdata = external
)
external_overall <- external_calibration$models
print(
  external_overall[external_overall$model == "onset+progression+outcome",
                   c("horizon", "observed_events", "expected_events",
                     "oe_ratio", "oe_ratio_lower", "oe_ratio_upper",
                     "expected_risk", "observed_risk", "slope", "slope_se")],
  digits = 3
)
```

```{r guard-external-calibration, include = FALSE}
ext_10 <- external_overall[
  external_overall$model == "onset+progression+outcome" &
    external_overall$horizon == 10,
]
# The text below states that the models predict far too few events, and that
# the slope is far below 1.
stopifnot(
  ext_10$oe_ratio_lower > 2,
  ext_10$slope + 2 * ext_10$slope_se < 1
)
```

The models predict far too few events.
At ten years, the O/E ratio of CAD events is `r sprintf("%.2f", ext_10$oe_ratio)`.
The mean predicted cumulative incidence is `r round(100 * ext_10$expected_risk)`%, and the observed cumulative incidence is `r round(100 * ext_10$observed_risk)`%.
The tenths of predicted risk lie far above the diagonal:

```{r external-calibration, fig.alt = "Observed against predicted ten-year cumulative incidence of CAD in tenths of predicted risk in the second biobank. The points lie far above the diagonal."}
plot_calibration(
  external_calibration,
  horizon = 10,
  model = "onset+progression+outcome"
)
```

The slope is `r sprintf("%.2f", ext_10$slope)` (SE `r sprintf("%.2f", ext_10$slope_se)`), far below 1.
The predicted risks are too spread out: the scores separate the persons of the second biobank less than the model expects.
With the same case mix, this points to weaker effects of the scores (Debray et al. 2015).
The O/E ratio and the slope show two different problems: a different baseline risk and different effects.

Simulation check (not possible with real data): both findings are correct, because the simulator gives the second biobank a higher rate of CAD and weaker effects of the scores.

The next step is to report this miscalibration, and then to update the model in the external cohort.
An update can estimate a new baseline hazard, or also a new calibration slope of the linear predictor (van Houwelingen 2000; Debray et al. 2015).
Report the performance before and after the update.
The package does not update models.

A different baseline risk alone does not change how the models rank the persons, but weaker effects do.
`rbind()` puts the corrected AUC of the first biobank next to the AUC of the second:

```{r}
external_discrimination <- assess_discrimination(
  imported,
  horizons,
  pairs = primary,
  newdata = external
)
both <- rbind(discrimination$models, external_discrimination$models)
print(
  both[both$model == "onset+progression+outcome" & both$horizon == 10,
       c("evaluation", "n", "auc", "auc_lower", "auc_upper")],
  digits = 3,
  row.names = FALSE
)
```

```{r guard-external-auc, include = FALSE}
auc_external <- external_discrimination$models$auc[
  external_discrimination$models$model == "onset+progression+outcome" &
    external_discrimination$models$horizon == 10
]
# The text below states this result.
stopifnot(auc_external < auc_full)
```

The AUC of the second biobank, `r sprintf("%.3f", auc_external)`, is lower than the corrected AUC of the first, `r sprintf("%.3f", auc_full)`, because the scores have weaker effects there.

Each result has the same elements and columns in every mode, and the same element of two results stacks with `rbind()`.
In the external mode, the change in AUC between two models has a confidence interval again, because the fitted models are fixed for new persons:

```{r}
print(
  external_discrimination$pairs[c("full", "reduced", "horizon", "delta_auc",
                                  "delta_auc_lower", "delta_auc_upper")],
  digits = 3
)
```

```{r guard-external-delta-auc, include = FALSE}
external_delta <- external_discrimination$pairs[
  external_discrimination$pairs$horizon == 10,
]
# The text below states these results, in the order of the primary pairs.
stopifnot(
  identical(as.character(external_delta$reduced), c("onset", "outcome")),
  external_delta$delta_auc > 0,
  external_delta$delta_auc_lower < 0,
  external_delta$delta_auc_upper > 0
)
```

At ten years, both estimates of the change in AUC are above zero, and both intervals include zero.
The interval is `r sprintf("%.3f", external_delta$delta_auc_lower[1])` to `r sprintf("%.3f", external_delta$delta_auc_upper[1])` for progression added to the T2D score, and `r sprintf("%.3f", external_delta$delta_auc_lower[2])` to `r sprintf("%.3f", external_delta$delta_auc_upper[2])` for progression added to the CAD score.
The second biobank then does not show that the progression score improves the AUC of these models.
These intervals describe the models of the first biobank, with their coefficients.

### The primary tests in the second biobank

The results above evaluate predictions.
A second question is whether the association of Step 4 replicates.
The second biobank has its own person-level data, and it can fit its own ladder and repeat the primary tests:

```{r}
external_ladder <- fit_ladder(external)
print(compare_ladder(external_ladder, pairs = primary), digits = 3)
external_hr <- summary(external_ladder)
print(
  external_hr[external_hr$model == "onset+progression+outcome",
              c("term", "hr", "hr_lower", "hr_upper")],
  digits = 3
)
```

```{r guard-external-tests, include = FALSE}
# The text below states these results. The hazard ratios of the simulator for
# the second biobank are 1.05, 1.12, and 1.15.
local({
  tests <- compare_ladder(external_ladder, pairs = primary)
  truth <- c(onset = 1.05, progression = 1.12, outcome = 1.15)
  full_hr <- external_hr[external_hr$model == "onset+progression+outcome", ]
  stopifnot(
    tests$p_value < 0.001,
    full_hr$hr_lower <= truth[full_hr$term],
    full_hr$hr_upper >= truth[full_hr$term],
    full_hr$hr[full_hr$term == "progression"] < hr_progression
  )
})
```

Both primary tests have a p-value below 0.001, and the association replicates.
The hazard ratios are smaller than in the first biobank.
`external` has the scores on the scale of the first biobank.
Each hazard ratio is then per SD of the first biobank, and the two sets of hazard ratios are comparable.
The fit uses the design `"t1"`, the design of the analysis plan for hazard ratios.

Simulation check (not possible with real data): each interval contains the value of the simulator for the second biobank (1.05, 1.12, and 1.15).


## Step 10: start follow-up at the T2D diagnosis

Many analyses start the follow-up of every person at the T2D diagnosis (`index = "diagnosis"`), also for persons who joined the biobank later.
The early follow-up of a prevalent person, between diagnosis and recruitment, was then not observed as it happened.
We know about it only because the person was alive when they joined.
Anyone who died before recruitment cannot be in the biobank.
Under the diagnosis index, this selection removes early deaths, so the risk of death is too low.
The risk of CAD is too high when most of the removed deaths are from other causes, and too low when fatal CAD is common.
The recruitment index of Step 1 removes the bias in both directions.

Settings are easiest to vary by building each cohort from one list of arguments.
The list holds the arguments of Step 1 apart from the index:

```{r}
setup <- list(
  data = biobank,
  id = "IID",
  status = "exit_status",
  age_at_diagnosis = "age_at_diagnosis",
  age_at_exit = "age_at_exit",
  scores = c(
    onset = "t2d_score",
    progression = "t2d_to_cad_score",
    outcome = "cad_score"
  ),
  covariates = c("sex", "PC1", "PC2", "PC3", "PC4", "PC5", "birth_year"),
  status_codes = list(censored = 0, event = 1, competing = 2),
  event_labels = c(event = "CAD", competing = "Death")
)
recruitment <- list(
  index = "recruitment",
  age_at_recruitment = "age_at_recruitment"
)
cohort_diagnosis <- do.call(prs_cohort, c(setup, list(index = "diagnosis")))
```

Simulation check (not possible with real data): a cohort of every simulated person, also those who died before they could join, estimates the true risks.
The table gives the ten-year cumulative incidence of CAD and of death under each analysis:

```{r}
setup_all <- setup
setup_all$data <- simulated
cohort_all <- do.call(prs_cohort, c(setup_all, list(index = "diagnosis")))
ten_year <- do.call(rbind, lapply(
  list(
    "all simulated persons" = cohort_all,
    "recruitment index" = cohort,
    "diagnosis index" = cohort_diagnosis
  ),
  function(x) {
    cumulative_risk(x, score = NULL, times = 10)$risks[
      c("n", "risk", "risk_lower", "risk_upper", "competing_risk",
        "competing_risk_lower", "competing_risk_upper")
    ]
  }
))
print(ten_year, digits = 2)
```

```{r guard-diagnosis-index, include = FALSE}
death_gap <- ten_year["all simulated persons", "competing_risk"] -
  ten_year["diagnosis index", "competing_risk"]
cad_gap <- ten_year["diagnosis index", "risk"] -
  ten_year["all simulated persons", "risk"]
# The text below states these results.
local({
  truth <- ten_year["all simulated persons", ]
  at_recruitment <- ten_year["recruitment index", ]
  at_first_diagnosis <- ten_year["diagnosis index", ]
  stopifnot(
    abs(at_recruitment$risk - truth$risk) < 0.02,
    abs(at_recruitment$competing_risk - truth$competing_risk) < 0.02,
    cad_gap > 0,
    at_first_diagnosis$risk < truth$risk_upper,
    at_first_diagnosis$competing_risk_upper < truth$competing_risk_lower
  )
})
```

The recruitment index agrees with the truth.
The diagnosis index gives a risk of death that is too low by `r round(100 * death_gap)` percentage points.
Its risk of CAD is `r sprintf("%.1f", 100 * cad_gap)` points higher than the truth, the direction expected when the removed deaths are not from CAD.
In a cohort of this size, that difference is within the confidence intervals, and the risk of death shows the bias clearly.
A model fitted under the diagnosis index carries this bias into its absolute risks.
The bias is in the baseline hazards, so the hazard ratios change much less:

```{r}
analyses <- list(
  "recruitment index" = cohort,
  "diagnosis index" = cohort_diagnosis,
  "incident only" = do.call(
    prs_cohort,
    c(setup, recruitment, list(prevalent = "drop"))
  ),
  "1-year landmark" = do.call(
    prs_cohort,
    c(setup, recruitment, list(entry_lag = 1))
  )
)

progression_hr <- do.call(rbind, lapply(names(analyses), function(analysis) {
  fit <- fit_ladder(analyses[[analysis]])
  row <- summary(fit)
  row <- row[
    row$model == "onset+progression+outcome" & row$term == "progression",
  ]
  data.frame(
    analysis = analysis,
    persons = nrow(analyses[[analysis]]$data),
    events = sum(analyses[[analysis]]$data$status),
    row[c("hr", "hr_lower", "hr_upper")],
    row.names = NULL
  )
}))
print(progression_hr, digits = 3)
```

Each fit uses the default design `"t1"`, and the comparison of the indexes stays within one design.
The other designs need the recruitment index, with incident and prevalent persons.
The last two rows answer common reviewer questions.

- `prevalent = "drop"` keeps only incident persons, who were followed from their diagnosis.
- `entry_lag` starts each person's follow-up no earlier than a fixed time after the diagnosis. A lag of one year ignores CAD found during the work-up of a new T2D diagnosis. It gives a landmark analysis of persons still alive and free of CAD one year after diagnosis.

```{r guard-index-hr, include = FALSE}
incident_hr <- progression_hr[progression_hr$analysis == "incident only", ]
# The text below states these results. The true hazard ratio is 1.25.
stopifnot(
  max(progression_hr$hr_lower) < min(progression_hr$hr_upper),
  progression_hr$hr_lower <= 1.25,
  progression_hr$hr_upper >= 1.25,
  which.max(progression_hr$hr_upper / progression_hr$hr_lower) ==
    which(progression_hr$analysis == "incident only")
)
```

The intervals of the HR of the progression score, adjusted for the other two, overlap in the four analyses.
The incident-only estimate, `r sprintf("%.2f", incident_hr$hr)` (95% CI `r sprintf("%.2f", incident_hr$hr_lower)` to `r sprintf("%.2f", incident_hr$hr_upper)`), is the least precise, because it has the fewest events.
With `r sum(analyses[["incident only"]]$data$status)` CAD events, its largest model has fewer than 10 events per coefficient, and `fit_ladder()` warns.

Simulation check (not possible with real data): each of the four intervals contains the true value of 1.25.

Report the recruitment index as the main analysis.
The diagnosis index can be a sensitivity analysis of the hazard ratios, but not of the absolute risks.

## Step 11: deaths treated as censoring

Many analyses treat a death as the end of follow-up, as if the person had left the study.
This step shows what changes.
A death treated as censoring gives estimates of net risks.
A net risk is the risk in a world where nobody dies before CAD.
Net risks are too high for a real patient (Wolbers et al. 2009).
These estimates are net risks only if death is independent of CAD, and the data cannot check this assumption (Tsiatis 1975).

We build a second cohort from the same persons, with death as censoring:

```{r}
# Assign the whole list of status codes. modifyList() merges the two lists,
# keeps `competing = 2`, and prs_cohort() then stops with an error.
setup_censored <- setup
setup_censored$status_codes <- list(censored = c(0, 2), event = 1)
setup_censored$event_labels <- c(event = "CAD")
cohort_censored <- do.call(prs_cohort, c(setup_censored, recruitment))
cohort_censored
```

### Risk by score group

With death as censoring, `cumulative_risk()` gives one minus the Kaplan-Meier estimate.
That estimate keeps a positive risk of CAD for a person who died (Putter et al. 2007).
The table below gives one minus Kaplan-Meier minus the cumulative incidence, in percentage points, by quintile and time:

```{r}
incidence <- cumulative_risk(
  cohort,
  score = "progression",
  times = c(5, 10, 15),
  groups = 5L
)
one_minus_km <- cumulative_risk(
  cohort_censored,
  score = "progression",
  times = c(5, 10, 15),
  groups = 5L
)
gap <- incidence$risks[c("group", "time")]
stopifnot(identical(one_minus_km$risks[c("group", "time")], gap))
gap$points <- 100 * (one_minus_km$risks$risk - incidence$risks$risk)
round(xtabs(points ~ group + time, gap), 1)
```

```{r guard-km-gap, include = FALSE}
# The text below states that one minus Kaplan-Meier is larger in every cell.
stopifnot(gap$points > 0)
```

One minus Kaplan-Meier is larger in every cell.
The gap grows with the number of deaths before the time: it is at most `r sprintf("%.1f", max(gap$points[gap$time == 5]))` points at 5 years and up to `r sprintf("%.1f", max(gap$points[gap$time == 15]))` points at 15 years.

### Absolute risk

The models for CAD are the same in both cohorts, because each treats a death as censoring.
Only the absolute risks differ.

Simulation check (not possible with real data): the code below compares the mean predicted ten-year risk at the diagnosis with the mean true risk, which is a cumulative incidence.
The simulator gives the true risk from the diagnosis only.
`prediction_data()` without `age_at_entry` predicts from the diagnosis, also for prevalent persons.
For a risk from the diagnosis, the rung of `ladder` is the more precise one in these data (Step 7), and this comparison uses it:

```{r}
persons_raw <- biobank[match(cohort$data$id, biobank$IID), ]
at_diagnosis <- prediction_data(persons_raw, cohort)
true_risk <- persons_raw$true_risk
report_fit_t1 <- ladder$models[["onset+progression+outcome"]]
round(c(
  true = mean(true_risk),
  cumulative_incidence = mean(
    predict_risk(report_fit_t1, at_diagnosis, horizon = 10)
  ),
  net_risk = mean(
    predict_risk(report_fit_t1, at_diagnosis, horizon = 10, competing = FALSE)
  )
), 3)
```

```{r guard-net-risk, include = FALSE}
# The text below states these results.
local({
  truth <- mean(true_risk)
  incidence_mean <- mean(
    predict_risk(report_fit_t1, at_diagnosis, horizon = 10)
  )
  net_mean <- mean(
    predict_risk(report_fit_t1, at_diagnosis, horizon = 10, competing = FALSE)
  )
  stopifnot(
    abs(incidence_mean - truth) < 0.01,
    net_mean > truth,
    net_mean - incidence_mean > 0.005
  )
})
```

The mean cumulative incidence is close to the mean true risk, and the mean net risk is too high.

The difference matters at a guideline threshold.
`assess_thresholds()` gives the share of persons with a predicted risk of 20% or more.
These risks are from entry, and they use the design of the analysis plan for absolute risks, `checks$design`.
The first value uses the net risk of a ladder fitted with death as censoring.
The second value uses the cumulative incidence of the same rung of `risk_ladder`, with deaths as a competing event.
Both values are apparent:

```{r}
ladder_censored <- fit_ladder(cohort_censored, design = checks$design)
flagged <- function(ladder) {
  models <- assess_thresholds(ladder, horizons = 10, thresholds = 0.2)$models
  models$above[models$model == "onset+progression+outcome"]
}
above <- c(
  net_risk = flagged(ladder_censored),
  cumulative_incidence = flagged(risk_ladder)
)
round(above, 3)
```

```{r guard-flagged-net-risk, include = FALSE}
# The text below states that the net risk flags more persons.
stopifnot(above[["net_risk"]] > above[["cumulative_incidence"]])
```

The net risk flags `r round(100 * above[["net_risk"]])`% of the persons, and the cumulative incidence flags `r round(100 * above[["cumulative_incidence"]])`%.
A model that ignores death puts too many persons above a treatment threshold (Wolbers et al. 2009).

### Fine-Gray models

The Fine-Gray model is a different regression model for the cumulative incidence (Fine and Gray 1999).
It is not part of the ladder.
The article [Fine-Gray models with the data of a cohort](https://brettklamer.com/work/prsladder/articles/fine-gray.html) gives the reasons, and it shows how to fit a Fine-Gray model for comparison.

## Step 12: answer the question

The question was: among people who already have T2D, does the progression score tell us who will develop CAD, beyond what the T2D score and the CAD score already tell us?
The table collects the results that answer it.

```{r guard-answer, include = FALSE}
answer_hr <- hr[
  hr$model == "onset+progression+outcome" & hr$term == "progression",
]
answer_risk <- adjusted$risks[adjusted$risks$prob %in% c(0.1, 0.9), ]
answer_auc <- discrimination$pairs[discrimination$pairs$horizon == 10, ]
answer_nb <- nb_pairs[nb_pairs$full == "onset+progression+outcome" &
                        nb_pairs$threshold %in% c(0.1, 0.2), ]
answer_external <- external_discrimination$pairs[
  external_discrimination$pairs$horizon == 10,
]
primary_full <- c("onset+progression", "progression+outcome")
# Each object holds the rows of the table, in the order of the table. The
# other conditions are the results that the text below states.
stopifnot(
  nrow(answer_hr) == 1L,
  nrow(answer_risk) == 2L,
  nrow(answer_nb) == 2L,
  identical(as.character(answer_auc$full), primary_full),
  identical(as.character(answer_external$full), primary_full),
  answer_hr$hr_lower <= 1.25,
  answer_hr$hr_upper >= 1.25,
  answer_auc$delta_auc > 0,
  identical(answer_nb$threshold, c(0.1, 0.2)),
  abs(answer_nb$delta_net_benefit[1]) < 0.001,
  answer_nb$delta_net_benefit[2] > 0,
  answer_external$delta_auc_lower < 0,
  answer_external$delta_auc_upper > 0
)
```

| Result | Estimate | Step |
|---|---|---|
| HR for CAD per SD of the progression score, adjusted for the other two scores | `r sprintf("%.2f", answer_hr$hr)` (95% CI `r sprintf("%.2f", answer_hr$hr_lower)` to `r sprintf("%.2f", answer_hr$hr_upper)`) | 3 |
| Likelihood-ratio tests of the two primary comparisons | Both p-values below 0.001 | 4 |
| Likelihood-ratio test of the secondary comparison | p-value below 0.001 | 4 |
| Ten-year risk of CAD at the 10th and the 90th percentile of the progression score, report model | `r round(100 * answer_risk$risk[1])`% (95% CI `r round(100 * answer_risk$risk_lower[1])`% to `r round(100 * answer_risk$risk_upper[1])`%) and `r round(100 * answer_risk$risk[2])`% (`r round(100 * answer_risk$risk_lower[2])`% to `r round(100 * answer_risk$risk_upper[2])`%) | 7 |
| Change in AUC at ten years, for progression added to the T2D score and to the CAD score | `r sprintf("%.3f", answer_auc$delta_auc[1])` and `r sprintf("%.3f", answer_auc$delta_auc[2])` | 8 |
| Change in net benefit at 10% and 20%, for progression added to the other two scores | `r sprintf("%.1f", 1000 * answer_nb$delta_net_benefit[1])` and `r sprintf("%.1f", 1000 * answer_nb$delta_net_benefit[2])` true positives per 1000 persons | 8 |
| Second biobank: O/E ratio of the report model at ten years | `r sprintf("%.2f", ext_10$oe_ratio)` (95% CI `r sprintf("%.2f", ext_10$oe_ratio_lower)` to `r sprintf("%.2f", ext_10$oe_ratio_upper)`) | 9 |
| Second biobank: change in AUC at ten years, for the same two additions | `r sprintf("%.3f", answer_external$delta_auc[1])` (`r sprintf("%.3f", answer_external$delta_auc_lower[1])` to `r sprintf("%.3f", answer_external$delta_auc_upper[1])`) and `r sprintf("%.3f", answer_external$delta_auc[2])` (`r sprintf("%.3f", answer_external$delta_auc_lower[2])` to `r sprintf("%.3f", answer_external$delta_auc_upper[2])`) | 9 |
| Second biobank: likelihood-ratio tests of the two primary comparisons, in its own ladder | Both p-values below 0.001 | 9 |

The answer is yes, and the size of the gain is modest.
The progression score is associated with the rate of CAD beyond the T2D score, beyond the CAD score, and beyond both scores together.
The added information changes the predictions by a small amount.
The net benefit rises at the threshold of 20%, and it is almost unchanged at 10%.
The changes in AUC and in net benefit of the first biobank are sizes of effect without a confidence interval (Step 8).
The second biobank replicates the two primary tests in its own ladder.
For the predictions of the models of the first biobank, the intervals of both changes in AUC include zero there.
The models predict too few events there, and they need an update before use in that population (Step 9).

Simulation check (not possible with real data): the simulator gave the progression score a hazard ratio of 1.25, and the interval of the first row contains this value.

All values of this tutorial come from simulated data.
They show the form of a report, and they are not results about T2D and CAD.

## Reporting checklist

The PRS-RS reporting standard (Wand et al. 2021) lists what a PRS risk-model paper must report.
prsladder gives the items that come from the analysis:

| PRS-RS item | Where to find it |
|---|---|
| Study design and recruitment | The analysis plan, and `summary(cohort)` for the incident and prevalent persons and the median follow-up |
| Participant demographics and clinical characteristics | `describe_cohort(cohort, external)`, with one column per cohort |
| Non-genetic variables | `cohort$covariates`, the age spline, and the terms of the design |
| Outcome of interest | `summary(cohort)` for the number of events of each kind |
| Missing data | `cohort$n_excluded`, the persons removed for each reason |
| Risk model type | Cox models on the years since the first diagnosis, or on the years since entry under the design `"index"`, and the absolute risk at each horizon from entry. With a competing event, the models are cause-specific and the risk is the cumulative incidence |
| Integrated risk model(s) description and fitting | `fit_ladder()`: each model has its scores, the covariates, and the age spline. `formula()` of a model gives its terms, and `ladder$coefficients` gives every coefficient |
| PRS distribution | `summary(cohort)`, `cohort$scores`, and `describe_cohort()` |
| Risk model predictive ability | `summary(ladder)` for the HR per SD, `compare_ladder()`, `plot_ladder()`, and `standardized_risk()` for the absolute risk at score percentiles |
| Risk model discrimination | `assess_discrimination()`, `plot_discrimination()`, and `assess_thresholds()` |
| Risk model calibration | `assess_calibration()`, `plot_calibration()`, and `shrinkage` for the corrected slope |
| Subgroup analyses | The package has no function for subgroups. Name the subgroups in the analysis plan, build one cohort for each subgroup with `prs_cohort()`, and fit a ladder in each cohort |
| Risk model intended uses | `assess_net_benefit()` and `plot_net_benefit()` for the net benefit at the clinical thresholds |

A report of this kind of study also needs results that PRS-RS does not name:

| Other result | Where to find it |
|---|---|
| Risk by score group, figure | `plot_risk()` with its risk table, `plot_stacked()`, and `plot_gradient()` for many groups |
| Risk by score group, table | `describe_risk()`, and `ratio` and `case_share` in `cumulative_risk()$risks` |
| Time lost by score group | `cumulative_risk()$time_lost` |
| Competing events | The hazard ratios and the cumulative incidence of both events (Latouche et al. 2013): `ladder$competing_coefficients`, `plot_ladder(ladder, event = "both")`, `plot_risk(..., event = "both")`, and `describe_risk()` |
| Design of the absolute risks | The three rows of `check_design()$tests`, and the statement that the analysis plan fixed the rule |
| Mode of each evaluation | The `evaluation` column of every result: apparent, optimism-corrected, or external |
| Internal validation | The number of replicates and the seed of `bootstrap_plan()`, and the Monte Carlo standard errors of the optimism |
| External validation | The settings of `export_ladder()`, such as `min_cell_count`, and `describe_case_mix()` for the two cohorts |

You must supply the rest: how the scores were built and their PGS Catalog IDs, the genotyping and quality control, the definitions of T2D, CAD, and death (for example ICD codes), and the ancestry of every sample.
State whether the GWAS behind a score included participants of your biobank.
Any overlap inflates the results (Choi et al. 2020).

## Limitations to report

- Index event bias: everyone in the cohort has T2D. A score that raises the risk of T2D, such as the T2D score, can then show a distorted HR for CAD (Dudbridge et al. 2019). This bias occurs when T2D and CAD share causes that the model does not include. The package does not correct this.
- Censoring: the Kaplan-Meier, Aalen-Johansen, and censoring-weighted estimates (Gerds and Schumacher 2006) assume that the end of follow-up is unrelated to the scores and covariates, and, under the recruitment index, to the time of entry. The observed-to-expected ratio does not need the last assumption.
- Competing risks: with death declared as a competing event, the absolute risks are cumulative incidences from cause-specific Cox models. A hazard ratio for CAD then describes the rate of CAD among persons alive and free of CAD, and a score can also change the cumulative incidence through the rate of death. With death declared as censoring, the absolute risks are net risks, which assume that death is unrelated to CAD (Step 11).
- Design: the absolute risks depend on an assumption about the hazard after entry (`?fit_ladder`). The rule of `check_design()` chooses the design from the data, and the bootstrap does not repeat this choice. The rule has not been studied with competing events or with an entry lag (`?check_design`).
- Ancestry: a score usually predicts less well in ancestries that differ from its GWAS (Choi et al. 2020). Report results for each ancestry group, and do not assume that they transfer.
- Model choice: do not choose the best score, threshold, or rung of the ladder on the same persons whose results you then report. The bootstrap repeats only the steps of `fit_ladder()`, and it does not correct for such a choice (`?validation`).
- Intervals of the corrected mode: the analytic interval does not include the error of the estimated optimism. One simulation check measured its coverage for the AUC and the Brier score (`?validation`). For a model with an AUC near 0.5, the 95% interval of the corrected AUC covered the true value in 89% to 90% of the samples. The coverage of the other intervals is not measured.

## References

Aalen OO, Johansen S. An empirical transition matrix for non-homogeneous Markov chains based on censored observations. Scandinavian Journal of Statistics 1978;5:141-150.

Austin PC, Harrell FE Jr, van Klaveren D. Graphical calibration curves and the integrated calibration index (ICI) for survival models. Statistics in Medicine 2020;39:2714-2742.

Benichou J, Gail MH. Estimates of absolute cause-specific risk in cohort studies. Biometrics 1990;46:813-826.

Blanche P, Dartigues J-F, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Statistics in Medicine 2013;32:5381-5397.

Chiou SH, Qian J, Mormino E, Betensky RA. Permutation tests for general dependent truncation. Computational Statistics and Data Analysis 2018;128:308-324.

Choi SW, Mak TS-H, O'Reilly PF. Tutorial: a guide to performing polygenic risk score analyses. Nature Protocols 2020;15:2759-2772.

Collins GS, Ogundimu EO, Altman DG. Sample size considerations for the external validation of a multivariable prognostic model: a resampling study. Statistics in Medicine 2016;35:214-226.

Crowson CS, Atkinson EJ, Therneau TM. Assessing calibration of prognostic risk scores. Statistical Methods in Medical Research 2016;25:1692-1706.

Debray TPA, Vergouwe Y, Koffijberg H, Nieboer D, Steyerberg EW, Moons KGM. A new framework to enhance the interpretation of external validation studies of clinical prediction models. Journal of Clinical Epidemiology 2015;68:279-289.

Dudbridge F, Allen RJ, Sheehan NA, et al. Adjustment for index event bias in genome-wide association studies of subsequent events. Nature Communications 2019;10:1561.

Emerging Risk Factors Collaboration; Sarwar N, Gao P, Seshasai SR, et al. Diabetes mellitus, fasting blood glucose concentration, and risk of vascular disease: a collaborative meta-analysis of 102 prospective studies. Lancet 2010;375:2215-2222.

Fine JP, Gray RJ. A proportional hazards model for the subdistribution of a competing risk. Journal of the American Statistical Association 1999;94:496-509.

Fry A, Littlejohns TJ, Sudlow C, Doherty N, Adamska L, Sprosen T, Collins R, Allen NE. Comparison of sociodemographic and health-related characteristics of UK Biobank participants with those of the general population. American Journal of Epidemiology 2017;186:1026-1034.

Gerds TA, Schumacher M. Consistent estimation of the expected Brier score in general survival models with right-censored event times. Biometrical Journal 2006;48:1029-1040.

Geskus RB. Cause-specific cumulative incidence estimation and the Fine and Gray model under both left truncation and right censoring. Biometrics 2011;67:39-49.

Harrell FE Jr, Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine 1996;15:361-387.

Kerr KF, Brown MD, Zhu K, Janes H. Assessing the clinical impact of risk prediction models with decision curves: guidance for correct interpretation and appropriate use. Journal of Clinical Oncology 2016;34:2534-2540.

Latouche A, Allignol A, Beyersmann J, Labopin M, Fine JP. A competing risks analysis should report results on all cause-specific hazards and cumulative incidence functions. Journal of Clinical Epidemiology 2013;66:648-653.

Makuch RW. Adjusted survival curve estimation using covariates. Journal of Chronic Diseases 1982;35:437-443.

Marx N, Federici M, Schuett K, et al. 2023 ESC Guidelines for the management of cardiovascular disease in patients with diabetes. European Heart Journal 2023;44:4043-4140.

Pepe MS, Kerr KF, Longton G, Wang Z. Testing for improvement in prediction model performance. Statistics in Medicine 2013;32:1467-1482.

Putter H, Fiocco M, Geskus RB. Tutorial in biostatistics: competing risks and multi-state models. Statistics in Medicine 2007;26:2389-2430.

SCORE2-Diabetes Working Group and the ESC Cardiovascular Risk Collaboration. SCORE2-Diabetes: 10-year cardiovascular risk estimation in type 2 diabetes in Europe. European Heart Journal 2023;44:2544-2556.

Steyerberg EW, Harrell FE Jr, Borsboom GJ, Eijkemans MJ, Vergouwe Y, Habbema JD. Internal validation of predictive models: efficiency of some procedures for logistic regression analysis. Journal of Clinical Epidemiology 2001;54:774-781.

Tsai W-Y. Testing the assumption of independence of truncation time and failure time. Biometrika 1990;77:169-177.

Tsiatis A. A nonidentifiability aspect of the problem of competing risks. Proceedings of the National Academy of Sciences 1975;72:20-22.

van Houwelingen HC. Validation, calibration, revision and combination of prognostic survival models. Statistics in Medicine 2000;19:3401-3415.

van Houwelingen JC, le Cessie S. Predictive value of statistical models. Statistics in Medicine 1990;9:1303-1325.

Vickers AJ, Elkin EB. Decision curve analysis: a novel method for evaluating prediction models. Medical Decision Making 2006;26:565-574.

Vickers AJ, Van Calster B, Steyerberg EW. Net benefit approaches to the evaluation of prediction models, molecular markers, and diagnostic tests. BMJ 2016;352:i6.

Wand H, Lambert SA, Tamburro C, et al. Improving reporting standards for polygenic scores in risk prediction studies. Nature 2021;591:211-219.

Wolbers M, Koller MT, Witteman JCM, Steyerberg EW. Prognostic models with competing risks: methods and application to coronary risk prediction. Epidemiology 2009;20:555-561.