From type 2 diabetes to coronary artery disease: a prsladder tutorial

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:

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.

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.

The package counts time on two clocks.

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.

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.

# Keep the persons who were alive at recruitment.
biobank <- simulated[simulated$recruited, ]
sum(!simulated$recruited)
#> [1] 177
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
)
#>     IID t2d_score t2d_to_cad_score cad_score    sex age_at_diagnosis
#> 1 S0001     0.492           -0.720     0.147 Female             61.5
#> 2 S0002    -0.889            0.214    -0.439 Female             68.8
#> 3 S0003     0.222            1.584     1.451   Male             46.7
#>   age_at_exit exit_status age_at_recruitment
#> 1        75.7           2               64.8
#> 2        74.2           0               61.6
#> 3        66.2           0               52.4
table(biobank$exit_status)
#> 
#>    0    1    2 
#> 1667  550  606

The columns that matter are:

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:

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.

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> recruitment index
#> Persons: 2711  Events: 438  Competing events: 606
#> Prevalent persons: keep
#> Scores:
#>  name        column           sign
#>  onset       t2d_score        1   
#>  progression t2d_to_cad_score 1   
#>  outcome     cad_score        1   
#> Covariates: sex, PC1, PC2, PC3, PC4, PC5, birth_year
#> Excluded: historical = 112

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

cohort$scores
#>          name           column sign         mean        sd
#> 1       onset        t2d_score    1 -0.022290906 0.9832644
#> 2 progression t2d_to_cad_score    1 -0.013592420 0.9866685
#> 3     outcome        cad_score    1 -0.003607009 1.0055165

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.

entry_test <- check_entry(cohort)
entry_test
#>      n n_exits n_pairs        tau        z   p_value
#> 1 2711    1044 1575897 0.02326992 1.272315 0.2032612

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

summary(cohort)
#> Cohort, recruitment index
#> 
#> Persons, events, and median follow-up since entry:
#>      group    n n_events n_competing median_followup
#>        all 2711      438         606       13.242408
#>   incident  863      103         152        9.684496
#>  prevalent 1848      335         454       13.963115
#> 
#> Raw scores:
#>         name      mean     sd     q25    median    q75
#>        onset -0.022291 0.9833 -0.7038 -0.011270 0.6385
#>  progression -0.013592 0.9867 -0.6684 -0.006664 0.6504
#>      outcome -0.003607 1.0055 -0.6869 -0.014481 0.6601
#> 
#> Raw scores by status of the second event:
#>         name status    n      mean     sd
#>        onset      0 1667 -0.114312 0.9779
#>        onset      1  438  0.146540 0.9947
#>        onset      2  606  0.108816 0.9599
#>  progression      0 1667 -0.069970 0.9798
#>  progression      1  438  0.250548 0.9605
#>  progression      2  606 -0.049420 0.9949
#>      outcome      0 1667 -0.084866 0.9982
#>      outcome      1  438  0.298670 1.0471
#>      outcome      2  606  0.001444 0.9534
#> 
#> Excluded persons:
#>              missing          prior_event   no_recruitment_age 
#>                    0                    0                    0 
#>           historical    prevalent_dropped nonpositive_followup 
#>                  112                    0                    0

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.

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()
#> Found litedown! Enabling r-universe template

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.

overall_risk <- cumulative_risk(cohort, score = NULL, times = c(5, 10, 15))
overall_risk
#> <prs_cumulative_risk> whole cohort, Aalen-Johansen
#> Guards: min_events = 1, min_at_risk = 50
#>  group time n_at_risk n_events n_competing  risk lower upper competing_risk
#>      1    5      1756       85         117 0.061 0.049 0.075          0.082
#>      1   10      1774      244         335 0.125 0.110 0.141          0.169
#>      1   15      1085      373         518 0.182 0.165 0.200          0.250
#> Restricted mean time lost in `x$time_lost`

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.

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

ladder <- fit_ladder(cohort)
ladder
#> <prs_ladder> 8 models of 3 scores, recruitment index, t1 design, with competing events
#>                      model n_scores    n n_events    loglik      aic
#>                       base        0 2711      438 -3133.552 6287.104
#>                      onset        1 2711      438 -3121.678 6265.356
#>                progression        1 2711      438 -3114.111 6250.222
#>                    outcome        1 2711      438 -3106.294 6234.588
#>          onset+progression        2 2711      438 -3109.304 6242.607
#>              onset+outcome        2 2711      438 -3102.451 6228.902
#>        progression+outcome        2 2711      438 -3096.907 6217.813
#>  onset+progression+outcome        3 2711      438 -3095.375 6216.749
#>  concordance concordance_se events_per_parameter heuristic_shrinkage ph_p_value
#>    0.6162337     0.01403476             43.80000           0.8525664 0.02031953
#>    0.6320331     0.01371023             39.81818           0.8798806 0.01898256
#>    0.6418475     0.01345075             39.81818           0.8969158 0.02219477
#>    0.6507156     0.01340799             39.81818           0.9100894 0.01575074
#>    0.6469502     0.01335119             36.50000           0.8968400 0.01473862
#>    0.6554349     0.01340875             36.50000           0.9077128 0.01172659
#>    0.6628911     0.01320114             36.50000           0.9149648 0.02220847
#>    0.6643486     0.01324871             33.69231           0.9098361 0.01174094

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 33 CAD events per coefficient. Every ph_p_value is below 0.05, with values from 0.012 to 0.022. 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:

ph <- survival::cox.zph(
  ladder$models[["onset+progression+outcome"]],
  transform = "km"
)
round(ph$table, 3)
#>                                 chisq df     p
#> onset                           1.730  1 0.188
#> progression                     0.864  1 0.353
#> outcome                         0.094  1 0.759
#> sex                             0.000  1 0.996
#> PC1                             0.251  1 0.616
#> PC2                             2.393  1 0.122
#> PC3                             3.657  1 0.056
#> PC4                             4.175  1 0.041
#> PC5                             6.054  1 0.014
#> birth_year                      0.893  1 0.345
#> natural_spline(age_t1, df = 3)  3.212  3 0.360
#> GLOBAL                         27.183 13 0.012

The three scores have p-values of 0.19 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:

hr <- summary(ladder)
print(
  hr[c("model", "term", "hr", "hr_lower", "hr_upper", "p_value")],
  digits = 3
)
#>                        model        term   hr hr_lower hr_upper  p_value
#> 1                      onset       onset 1.26    1.151     1.39 1.12e-06
#> 2                progression progression 1.35    1.228     1.48 4.79e-10
#> 3                    outcome     outcome 1.43    1.300     1.57 1.03e-13
#> 4          onset+progression       onset 1.17    1.059     1.29 1.95e-03
#> 5          onset+progression progression 1.29    1.164     1.42 6.99e-07
#> 6              onset+outcome       onset 1.15    1.042     1.27 5.55e-03
#> 7              onset+outcome     outcome 1.37    1.240     1.51 4.48e-10
#> 8        progression+outcome progression 1.24    1.126     1.37 1.50e-05
#> 9        progression+outcome     outcome 1.35    1.219     1.48 3.58e-09
#> 10 onset+progression+outcome       onset 1.10    0.989     1.21 8.01e-02
#> 11 onset+progression+outcome progression 1.21    1.097     1.34 1.73e-04
#> 12 onset+progression+outcome     outcome 1.32    1.189     1.46 1.12e-07

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
)
#>           term    hr hr_lower hr_upper
#> 80       onset 1.248    1.144     1.36
#> 81 progression 0.929    0.853     1.01
#> 82     outcome 1.008    0.923     1.10

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.

primary <- list(
  c("onset+progression", "onset"),
  c("progression+outcome", "outcome")
)
primary_tests <- compare_ladder(ladder, pairs = primary)
print(primary_tests, digits = 3)
#>                  full reduced       added df chisq  p_value delta_aic
#> 1   onset+progression   onset progression  1  24.7 6.53e-07     -22.7
#> 2 progression+outcome outcome progression  1  18.8 1.47e-05     -16.8
#>   delta_concordance
#> 1            0.0149
#> 2            0.0122

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.

secondary <- list(c("onset+progression+outcome", "onset+outcome"))
secondary_test <- compare_ladder(ladder, pairs = secondary)
print(secondary_test, digits = 3)
#>                        full       reduced       added df chisq  p_value
#> 1 onset+progression+outcome onset+outcome progression  1  14.2 0.000169
#>   delta_aic delta_concordance
#> 1     -12.2           0.00891

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.

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
)
#>    group time   n n_events  risk difference difference_lower difference_upper
#> 1      1    5 543       12 0.038         NA               NA               NA
#> 2      1   10 543       32 0.080         NA               NA               NA
#> 3      2    5 542       16 0.059      0.021          -0.0155            0.057
#> 4      2   10 542       35 0.097      0.017          -0.0252            0.060
#> 5      3    5 542       15 0.052      0.014          -0.0206            0.049
#> 6      3   10 542       57 0.136      0.056           0.0119            0.100
#> 7      4    5 542       17 0.063      0.025          -0.0126            0.062
#> 8      4   10 542       57 0.143      0.063           0.0176            0.109
#> 9      5    5 542       25 0.089      0.050           0.0099            0.091
#> 10     5   10 542       63 0.165      0.085           0.0374            0.133
by_quintile$test
#>   statistic df      p_value
#> 1   32.7004  4 1.375634e-06

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

plot_risk(by_quintile, event = "both")

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

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:

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:

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

plot_gradient(by_decile, time = 10)

The rise between neighboring deciles is uneven, because each decile holds only 272 persons.

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 27% of the CAD cases by ten years (case_share in by_quintile$risks).

The article Risk by score group in depth 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:

?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:

checks <- check_design(cohort)
checks
#> $tests
#>             test null_model statistic df   p_value rejected
#> 1    entry_terms         t1  1.448719  2 0.4846350    FALSE
#> 2    entry_bands   t1_entry  3.763538  4 0.4389523    FALSE
#> 3 index_duration      index  0.661233  3 0.8822822    FALSE
#> 
#> $design
#> [1] "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 0.44, and the pick is "t1_entry". With 438 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:

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)
)
#>          term    t1  pick
#> 1       onset 1.096 1.094
#> 2 progression 1.214 1.214
#> 3     outcome 1.316 1.314

The hazard ratios of the two ladders differ by at most 0.002 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 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

horizon_support(cohort, horizons = c(5, 10, 15, 20))
#>   horizon    n n_events n_at_risk censoring_survival supported
#> 1       5 2711      181      2254        0.972117588      TRUE
#> 2      10 2711      345      1531        0.458958669      TRUE
#> 3      15 2711      434       253        0.006080652     FALSE
#> 4      20 2711      438         0        0.000000000     FALSE

horizons <- c(5, 10)

report_fit <- risk_ladder$models[["onset+progression+outcome"]]
round(predict_risk(report_fit, newdata = cohort$data[1:5, ], horizon = 10), 3)
#> [1] 0.151 0.122 0.216 0.118 0.052

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")]
#>   start age_t1 t1_status t1_duration progression
#> 1     0     55  incident           0  -0.9997355
#> 2     0     55  incident           0   1.0272877
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)
#>                       [,1]  [,2]
#> cumulative_incidence 0.096 0.141
#> net_risk             0.104 0.151
#> death                0.146 0.124

persons$age_now <- 58
later <- prediction_data(persons, cohort, age_at_entry = "age_now")
later[c("start", "t1_status", "t1_duration")]
#>   start t1_status t1_duration
#> 1     3 prevalent           3
#> 2     3 prevalent           3

later <- prediction_data(
  persons,
  cohort,
  age_at_entry = "age_now",
  t1_status = "incident"
)
later[c("start", "t1_status", "t1_duration")]
#>   start t1_status t1_duration
#> 1     3  incident           0
#> 2     3  incident           0
woman_later <- predict_risk(report_fit, later, horizon = 10)
round(woman_later, 3)
#> [1] 0.108 0.157

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
)
#>   prob raw_score time  risk competing_risk difference ratio
#> 1  0.1   -1.2840   10 0.097           0.21         NA    NA
#> 2  0.5   -0.0067   10 0.140           0.20      0.042   1.4
#> 3  0.9    1.2459   10 0.196           0.19      0.098   2.0

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
)
#>   prob risk risk_lower risk_upper difference difference_lower difference_upper
#> 1  0.1 0.10      0.070       0.14         NA               NA               NA
#> 2  0.5 0.13      0.095       0.18      0.029            0.016            0.046
#> 3  0.9 0.17      0.122       0.22      0.064            0.033            0.107

Step 8: internal validation

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
)
#>                        model evaluation   auc auc_lower auc_upper brier
#> 9                       base   apparent 0.589     0.555     0.622 0.114
#> 10                     onset   apparent 0.600     0.567     0.633 0.113
#> 11               progression   apparent 0.626     0.593     0.658 0.112
#> 12                   outcome   apparent 0.638     0.605     0.670 0.112
#> 13         onset+progression   apparent 0.628     0.596     0.661 0.112
#> 14             onset+outcome   apparent 0.639     0.607     0.672 0.112
#> 15       progression+outcome   apparent 0.654     0.622     0.686 0.111
#> 16 onset+progression+outcome   apparent 0.654     0.622     0.686 0.111

plan <- bootstrap_plan(n_boot = 30L, seed = 2026L)

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
)
#>                      model         evaluation   auc auc_lower auc_upper brier
#>  onset+progression+outcome           apparent 0.654     0.622     0.686 0.111
#>  onset+progression+outcome optimism-corrected 0.632     0.600     0.664 0.112

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
)
#>                        model apparent optimism optimism_mcse n_boot_used
#> 21                      base    0.589   0.0211       0.00298          30
#> 22                     onset    0.600   0.0204       0.00287          30
#> 23               progression    0.626   0.0204       0.00256          30
#> 24                   outcome    0.638   0.0207       0.00301          30
#> 25         onset+progression    0.628   0.0211       0.00239          30
#> 26             onset+outcome    0.639   0.0213       0.00285          30
#> 27       progression+outcome    0.654   0.0211       0.00282          30
#> 28 onset+progression+outcome    0.654   0.0221       0.00270          30

plot_discrimination(discrimination)

print(
  discrimination$pairs[c("full", "reduced", "horizon", "delta_auc",
                         "delta_auc_lower", "delta_auc_upper", "delta_brier")],
  digits = 3
)
#>                  full reduced horizon delta_auc delta_auc_lower delta_auc_upper
#> 1   onset+progression   onset       5    0.0165              NA              NA
#> 2 progression+outcome outcome       5    0.0124              NA              NA
#> 3   onset+progression   onset      10    0.0273              NA              NA
#> 4 progression+outcome outcome      10    0.0163              NA              NA
#>   delta_brier
#> 1   -0.000197
#> 2   -0.000115
#> 3   -0.000964
#> 4   -0.000544

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
)
#>    horizon oe_ratio expected_risk observed_risk slope
#> 8        5    0.984        0.0669        0.0668 1.022
#> 16      10    0.986        0.1309        0.1324 0.933

print(calibration$shrinkage[c("model", "slope")], digits = 3)
#>                       model slope
#> 1                      base 0.888
#> 2                     onset 0.897
#> 3               progression 0.901
#> 4                   outcome 0.914
#> 5         onset+progression 0.898
#> 6             onset+outcome 0.909
#> 7       progression+outcome 0.910
#> 8 onset+progression+outcome 0.905

plot_calibration(calibration, horizon = 10, model = "onset+progression+outcome")

net_benefit <- assess_net_benefit(
  risk_ladder,
  horizons = 10,
  thresholds = (5:30) / 100,
  pairs = c(primary, secondary),
  validation = plan
)

report_pair <- c("onset+outcome", "onset+progression+outcome")
plot_net_benefit(net_benefit, horizon = 10, model = report_pair)

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
)
#>                         full       reduced threshold delta_net_benefit
#> 6          onset+progression         onset       0.1           9.9e-05
#> 16         onset+progression         onset       0.2           2.7e-03
#> 32       progression+outcome       outcome       0.1           3.5e-03
#> 42       progression+outcome       outcome       0.2           1.7e-03
#> 58 onset+progression+outcome onset+outcome       0.1          -9.7e-05
#> 68 onset+progression+outcome onset+outcome       0.2           1.6e-03

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.

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, ]

export <- export_ladder(risk_ladder)
names(export)
#> [1] "models"      "rules"       "development" "settings"

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

external <- prs_cohort(second_biobank, reference = imported)
external
#> <prs_cohort> recruitment index
#> Persons: 3405  Events: 1326  Competing events: 620
#> Prevalent persons: keep
#> Scores:
#>  name        column           sign
#>  onset       t2d_score        1   
#>  progression t2d_to_cad_score 1   
#>  outcome     cad_score        1   
#> Covariates: sex, PC1, PC2, PC3, PC4, PC5, birth_year
#> Excluded: historical = 412

describe_cohort(cohort, external, label = labels, digits = year_digits)

print(describe_case_mix(imported, external), digits = 3)
#>                       model lp_mean_development lp_sd_development
#> 1                      base              -0.191             0.411
#> 2                     onset              -0.172             0.475
#> 3               progression              -0.188             0.506
#> 4                   outcome              -0.158             0.546
#> 5         onset+progression              -0.172             0.528
#> 6             onset+outcome              -0.144             0.564
#> 7       progression+outcome              -0.168             0.586
#> 8 onset+progression+outcome              -0.157             0.592
#>   lp_mean_external lp_sd_external lp_mean_difference lp_sd_ratio
#> 1           -0.184          0.412           0.006842       1.001
#> 2           -0.164          0.479           0.007968       1.009
#> 3           -0.188          0.508          -0.000352       1.004
#> 4           -0.156          0.545           0.001856       0.998
#> 5           -0.170          0.528           0.001774       1.000
#> 6           -0.141          0.565           0.003140       1.002
#> 7           -0.170          0.585          -0.002590       0.998
#> 8           -0.158          0.590          -0.001194       0.998

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
)
#>    horizon observed_events expected_events oe_ratio oe_ratio_lower
#> 8        5             646             206     3.13            2.9
#> 16      10            1106             359     3.08            2.9
#>    oe_ratio_upper expected_risk observed_risk slope slope_se
#> 8            3.38         0.066         0.190 0.663   0.0677
#> 16           3.26         0.129         0.335 0.641   0.0522

plot_calibration(
  external_calibration,
  horizon = 10,
  model = "onset+progression+outcome"
)

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
)
#>          evaluation    n   auc auc_lower auc_upper
#>  optimism-corrected 2711 0.632      0.60     0.664
#>            external 3405 0.591      0.57     0.613

print(
  external_discrimination$pairs[c("full", "reduced", "horizon", "delta_auc",
                                  "delta_auc_lower", "delta_auc_upper")],
  digits = 3
)
#>                  full reduced horizon delta_auc delta_auc_lower delta_auc_upper
#> 1   onset+progression   onset       5   0.00481        -0.00830          0.0179
#> 2 progression+outcome outcome       5   0.00482        -0.00487          0.0145
#> 3   onset+progression   onset      10   0.00638        -0.00624          0.0190
#> 4 progression+outcome outcome      10   0.00474        -0.00441          0.0139

external_ladder <- fit_ladder(external)
print(compare_ladder(external_ladder, pairs = primary), digits = 3)
#>                  full reduced       added df chisq  p_value delta_aic
#> 1   onset+progression   onset progression  1  20.3 6.72e-06     -18.3
#> 2 progression+outcome outcome progression  1  15.9 6.72e-05     -13.9
#>   delta_concordance
#> 1           0.00269
#> 2           0.00227
external_hr <- summary(external_ladder)
print(
  external_hr[external_hr$model == "onset+progression+outcome",
              c("term", "hr", "hr_lower", "hr_upper")],
  digits = 3
)
#>           term   hr hr_lower hr_upper
#> 10       onset 1.06    0.996     1.12
#> 11 progression 1.11    1.051     1.18
#> 12     outcome 1.10    1.044     1.17

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:

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:

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)
#>                          n risk risk_lower risk_upper competing_risk
#> all simulated persons 3000 0.12       0.11       0.14           0.18
#> recruitment index     2711 0.12       0.11       0.14           0.17
#> diagnosis index       2823 0.13       0.12       0.14           0.12
#>                       competing_risk_lower competing_risk_upper
#> all simulated persons                 0.16                 0.19
#> recruitment index                     0.15                 0.19
#> diagnosis index                       0.11                 0.14

The recruitment index agrees with the truth. The diagnosis index gives a risk of death that is too low by 5 percentage points. Its risk of CAD is 0.8 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:

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
  )
}))
#> Warning in fit_ladder(analyses[[analysis]]): The largest model has 7.92 events
#> per coefficient, fewer than 10.
print(progression_hr, digits = 3)
#>            analysis persons events   hr hr_lower hr_upper
#> 1 recruitment index    2711    438 1.21     1.10     1.34
#> 2   diagnosis index    2823    550 1.25     1.14     1.37
#> 3     incident only     863    103 1.47     1.19     1.83
#> 4   1-year landmark    2686    425 1.20     1.08     1.33

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.

The intervals of the HR of the progression score, adjusted for the other two, overlap in the four analyses. The incident-only estimate, 1.47 (95% CI 1.19 to 1.83), is the least precise, because it has the fewest events. With 103 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:

# 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
#> <prs_cohort> recruitment index
#> Persons: 2711  Events: 438
#> Prevalent persons: keep
#> Scores:
#>  name        column           sign
#>  onset       t2d_score        1   
#>  progression t2d_to_cad_score 1   
#>  outcome     cad_score        1   
#> Covariates: sex, PC1, PC2, PC3, PC4, PC5, birth_year
#> Excluded: historical = 112

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)
#>      time
#> group   5  10  15
#>     1 0.2 0.9 2.1
#>     2 0.3 1.0 2.7
#>     3 0.2 1.5 3.0
#>     4 0.2 1.1 2.8
#>     5 0.4 1.6 4.1

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)
#>                 true cumulative_incidence             net_risk 
#>                0.124                0.123                0.138

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)
#>             net_risk cumulative_incidence 
#>                0.225                0.143

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.

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:

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

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.