Participants

For the present analysis, we included 4377 TrialNet Pathway to Prevention participants [15, 16] who had been genotyped using the Illumina T1DExomeChip array, had at least one islet autoantibody, and had a baseline OGTT performed less than a year after entering TrialNet. If the participant had a single positive antibody, this was confirmed on two successive tests. We excluded 355 participants who had been diagnosed with diabetes or who had hyperglycaemia, defined as fasting plasma glucose ≥7.0 mmol/l (≥126 mg/dl) 2 h plasma glucose ≥11.1 mmol/l (≥200 mg/dl), or HbA1c ≥48 mmol/mol (≥6.5%) at their baseline OGTT. We further removed 55 participants with incomplete data. A flowchart of the cohort selection is provided in electronic supplementary material (ESM) Fig. 1. The dataset included 49% female participants; 89% of participants were of white European ancestry, and socioeconomic information for the participants was unknown. The participants provided informed consent, and the study was approved by the ethics committee at each site.

Procedures Study protocol

Participants were initially screened for GADA, IAA and IA-2A. Only a subset of participants were tested for islet cell antibodies and ZnT8 A, and thus these two autoantibodies were not included in this analysis. Participants were monitored using HbA1c measurements and OGTT at 6- or 12-month intervals depending on estimated risk.

Diagnosis of type 1 diabetes

During follow-up, type 1 diabetes was diagnosed in participants with hyperglycaemia, i.e. fasting plasma glucose ≥7.0 mmol/l, 2 h plasma glucose ≥11.1 mmol/l after 75 g oral glucose, random plasma glucose ≥11.1 mmol/l, or an HbA1c ≥48 mmol/mol (6.5%). The participants were either symptomatic, or, if asymptomatic, met these thresholds on two separate occasions.

Autoantibody assays

GADA, IA-2A and IAA were measured by radioimmunoassay in the TrialNet Core Laboratory at the Barbara Davis Center (Aurora, CO), as previously described [17].

Metabolic measures

Participants underwent an OGTT (oral glucose dose 1.75 g/kg, maximum 75 g) after an overnight fast. C-peptide and glucose measurements were performed in the fasting state and 30, 60, 90 and 120 min later.

SNP genotyping and computation of the type 1 diabetes genetic risk score

Genotyping was performed using the T1DExomeChip array (Illumina). This is a custom array with more than 90,000 custom SNPs selected from regions of the genome that are robustly associated with autoimmune diseases, the Infinium CoreExome-24 version 1.1 BeadChip. The type 1 diabetes genetic risk score (GRS2) described by Sharp et al [18] includes 67 SNPs, 30 of which were directly genotyped. We imputed 32 SNPs with a median R2 of 0.997 (min 0.858; max 0.999) using whole-genome sequence data from the NHLBI TOPMed Imputation Server (https://imputation.biodatacatalyst.nhlbi.nih.gov/) and the multi-ethnic reference panel that includes 97,256 reference samples and >308 million genetic variants [19]. We also imputed five SNPs in the HLA region (rs72848653, R2=0.999; rs9266268, R2=0.999; rs16899379, R2=0.998; rs2524277, R2=0.995; rs9268500, R2=0.925) using HLA-TAPAS hosted by the Michigan Imputation Server (https://github.com/immunogenomics/HLA-TAPAS/) and the high-resolution HLA reference panel spanning five global populations (n=21,546) based on whole-genome sequencing data [20]. Code to generate the HLA interaction part of the GRS2 is freely available online (https://github.com/t2diabetesgenes/t1dgrs2).

Statistical analysis

We used baseline variables from four categories: clinical/demographic, genetic, metabolic and immunological. In the clinical category, we considered the age (absolute value or natural logarithm) at screening, self-reported gender and BMI (absolute value or age- and gender-adjusted z score). To harmonise adult and paediatric BMI values, BMI z scores were calculated using previously published tables for children [21] and an age of 20 years for adults, as previously reported [22]. In the immunological category, we either included IA-2A (positive, negative) or autoantibody combination (i.e. GADA, IAA, IA-2A, GADA–IAA, GADA–IA-2A, IAA–IA-2A, GADA–IAA–IA-2A). Metabolic variables included fasting C-peptide, AUC C-peptide, AUC glucose and 30 min C-peptide index (C-peptide30). The GRS2 was used to provide an estimated genetic risk for type 1 diabetes. ESM Table 1 summarises the main characteristics of the study participants, both overall and by type 1 diabetes stage at screening.

Primary outcome

The primary outcome was time to type 1 diabetes.

Modelling

We used two modelling approaches. The first is the Cox proportional hazard (CPH) model [23], which is a classical statistical method to measure the effect of a unit increase of a variable with respect to the hazard rate. The second is the survival random forest approach [24], an ensemble learning method that can capture non-linear effects and interactions but may increase the risk of overfitting.

We trained the model using half of the dataset (participants entering the study during or before January 2013) and validated the model on the second half of the dataset (observations entering the study from February 2013). Only validation results are shown.

To compare the modelling approaches, we used measures of predictive power adapted to the survival setting and able to correctly adjust from censoring, i.e. time-dependent receiver operating characteristic (ROC) AUC [25] and the time-dependent Brier score [26]. The formula for the Brier score is given in ESM Methods. Time-dependent ROC AUC measures the discriminative power of a model to predict an event between the present and a chosen future time horizon; a score of 1 indicates perfect discrimination while a score of 0.5 is equivalent to no discrimination. The Brier score measures the agreement between predicted and observed risk; a value closer to 0 indicates a better performance. The time-dependent ROC AUC and Brier score were calculated for each model at 2-, 3-, 5,- 7- and 10-year horizons. We generated a measure of variable importance to assess the relative predictive power of each variable (see [27] and ESM Methods). The results were computed using R version 4.3.2 with the packages survival [25, 28] and randomForestSRC [29]. p values to compare time-dependent ROC AUC were obtained using the approach developed by Blanche et al [25].

Generation of type 1 diabetes predictive models

For the variables used in this study, there were 114,688 combinations that may be used to fit a predictive model. To limit the number of combinations of variables, we imposed rules to avoid strong collinearity (such as BMI and BMI z score not appearing in the same model). A full list of the rules is presented in ESM Methods. In total, we fitted 1943 combinations of variables using the CPH model and survival random forest model for each stage of progression to type 1 diabetes. We compared the performance of these models to that of DPTRS [10], DPTRS60 [11], M120 [12], Index60 [13], CPH [14] and LR [14]. The formulae for each of these models are given in ESM Methods.

Identifying models with global best performance

A Pareto front [30] (see ESM Methods) was used to identify the model with the best trade-offs between competing performance measures (i.e. cost, participant time and predictive performance at 3-year horizon as measured by ROC AUC and Brier score for each stage).

Estimation of financial cost and participant time for each variable

The cost of acquiring data on specific predictors can vary greatly (see Table 1 and ESM Table 2 for a detailed cost breakdown). Some are easily accessible (e.g. age), while others, such as metabolic variables, may require sequential blood sampling during an OGTT, adding to the cost and time required from the participants and the research or healthcare teams. We followed previous guidelines to estimate the cost of each variable [31,32,33]. All costs were inflated to their net present value in 2024 US dollars. Direct costs were calculated using US Medicare reimbursement rates as these are a measure of central tendency on the private market and include costs of the test and healthcare provider time. Indirect costs were calculated using wage data from the US Bureau of Labor Statistics, and accounted for the value of each participant’s time, which was calculated as the lost wages of the participant or the parent/guardian who must accompany paediatric participants [33]. The cost of obtaining the GRS2 was estimated at US$20 based on previous literature [34] and the current cost in research studies. The time needed to obtain the predictors that require a medical procedure (such as an OGTT) was estimated as the duration of the procedure plus 45 min for preparation and discharge tasks. For example, a 120 min OGTT test requires 165 min to account for the time needed before the procedure, such as that required to obtain demographic and anthropometric characteristics and establish an i.v. line, and the time needed after the procedure, such as that required to retrieve the line and discharge the participant. The cost and time ranges apply to a single screening visit. When variables are used simultaneously and can be measured together, then we included assay costs for both but only the time of the longer test. For example, HbA1c can be collected during OGTT, and therefore additional time for this was not included. ESM Table 3 provides an impact inventory (a comprehensive economic assessment of a healthcare intervention [35]).

Table 1 Cost estimation for the acquisition of variable informationInteractive figures

To help the reader to compare predictive models of interest, we developed interactive figures that can be found online at https://lauricf.github.io/interactive_suplement.github.io/. The reader can zoom, identify models with selected variables, or highlight a model of interest across multiples stages.