Characterization of direct and/or indirect genetic associations for multiple traits in longitudinal studies of disease progression

Joint model for one longitudinal and one time-to-event trait

For each individual i, we define y_i = (y_i,1, …, y_i,j, …, y_i,J), as the vector of QT measures collected over the J visit times t_i = (t_i,1, …, t_i,j, …, t_i,J)^T with 1 ≤ j ≤ J and t_i,1 ≤ ⋯ ≤ t_i,j ≤ ⋯ ≤ t_i,J. We denote (T_i, δ_i) as the vector of right-censored event time T_i and event indicator δ_i for the time-to-event trait, and assume , where is the latent (uncensored) event time and C_i is the censoring time (e.g., administrative censoring). We define , with δ_i = 1 if the event occurs during the observation period, and δ_i = 0 otherwise.

Longitudinal sub-model

This is specified by a mixed-effects model for the longitudinal QT, based on the (Laird and Ware 1982) linear mixed model. The model builds on the assumption that for every individual in the sample there exists an underlying smooth trajectory of the longitudinal QT that describes the subject-specific evolution dependent on time, SNP effect, and individual-level random effects b_i. To simplify the presentation, we assume a linear QT trajectory (Equation 1), but the longitudinal sub-model can be adapted for nonlinear trajectories using, for example, higher order polynomials or splines (Rizopoulos 2012): Where:

• g_i, is the genotype of the individual i for the SNP being tested, coded as the number of copies of the minor allele,

• β = (β₀, β₁, β_g)^T is the vector of fixed intercept and slope time effects and fixed genetic effect on the longitudinal QT,

• b_i = (b_i,0, b_i,1)^T are the subject-specific random intercept and slope time effects assuming b_i∼N₂(0, D) and D is the variance-covariance matrix.

This trajectory cannot be observed directly, rather we observe longitudinal measurements y_i collected at discrete time points t_i ; measurements are subject to independent and identically distributed noise contamination variables ε_i∼N_J(0, Σ), where ε_i = (ε_i,1, …, ε_i,j, …, ε_i,J)^T, Σ = σ²I_J, with σ², the residual variance of the QT: We assume that b_i and ε_i are independent (Laird and Ware 1982).

Equation 2 implies that y_i in ℜ^J follows a multivariate normal distribution with: where X_i = (1_J, t_i, g_i1_J) denotes the (J- by-3) design matrix for the fixed intercept, slope, and SNP effects, and Z_i = (1_J, t_i) is the (J- by-2) design matrix for the random intercept and slope effects, with 1_J = (1, …,1, …,1)^T. To increase robustness to misspecification of the variance-covariance matrix D, we adopt an unstructured form for the random-effects variance, defined as , with the added benefit of not requiring additional constraints on the serial dependence between the repeated measurements for each individual. This choice implies that the covariance function between any pair of QT observations for individual i collected at two distinct visit times t_i,j ≠ t_i,s (1 ≤ j ≤ J and 1 ≤ s ≤ J, with j ≠ s) is given by Cov(y_i,j, y_i,s) = t_i,jt_i,sVar(b_i,1) + (t_i,j + t_i,s)Cov(b_i,0, b_i,1) + Var(b_i,0) + σ², with variance function , which is quadratic over time with positive curvature at Var(b_i,1).

Time-to-event sub-model

This is specified by a proportional hazards model (PH model), in which the hazard function of the time-to-event trait is defined as the instantaneous event rate in a small interval around given that the event has not occurred before time t, genetic effect and a function of the history of the true unobserved longitudinal process up to time t that is associated with risk of the event, , with . We specify that the hazard function (Equation 3) depends on the SNP effect adjusted for association of the longitudinal QT risk factor with the time-to-event trait. Where:

• λ₀(t) is a (parametric or non-parametric) baseline hazard function;

• specifies the function of the longitudinal QT trajectory accounting for trajectory values at time t that is associated with risk of the event. In the case of a contemporaneous parametrization, the hazard of an event at a time t depends on the longitudinal trajectory value at the same time . Other functional forms of the QT trajectory can weight earlier QT values according to prior knowledge of the relationship of the QT with the time-to-event trait (Hickey et al. 2016; Mauff et al. 2017; Papageorgiou et al. 2019);

• α is the effect of the longitudinal QT risk factor on the time-to-event trait;

• γ_g denotes the genetic effect on the time-to-event trait accounting for association of the longitudinal QT risk factor with the time-to-event trait.

Interpretation

As depicted in Fig. 1, the joint model parameters characterize relationships among a SNP, an intermediate QT, and a time-to-event trait and decompose possible effects of a SNP on a time-to-event trait into: an indirect effect induced via the SNP effect on the longitudinal QT; and a direct SNP effect independent of the QT (Ibrahim et al. 2010; Hickey et al. 2018b).

Based on the SNP effects, β_g and γ_g, and assuming α ≠ 0, a SNP association with a time-to-event trait can be one of three types:

• indirect SNP association: the SNP has a non-null effect on the longitudinal QT (β_g ≠ 0), but no effect on the time-to-event trait (γ_g = 0); the overall SNP effect θ depends on the indirect effect (θ = µ_g, with µ_g = αβ_g).

• direct SNP association: the SNP has a non-null effect on the time-to-event trait (γ_g ≠ 0), but no effect on the longitudinal QT (β_g = 0); the overall SNP effect depends only on the direct effect (θ = γ_g).

• both direct and indirect SNP associations: the SNP has non-null effects on the longitudinal risk factor (β_g ≠ 0) and on the time-to-event trait (γ_g ≠ 0). In this case, the overall SNP effect θ aggregates the indirect and direct SNP effects (θ = µ_g + γ_g, with µ_g = αβ_g).

In later subsections, we detail statistical implementation for estimation of the α, β_g and γ_g parameters and associated hypothesis testing in the joint model, which underlie the procedure we then propose to classify the SNP association as indirect, direct, or both direct and indirect.

When an associated longitudinal risk factor is omitted from the time-to-event model, the estimated SNP effect on the time-to-event trait captures the overall SNP effect (θ = µ_g + γ_g). This can occur in GWAS when the time-to-event analysis ignores an intermediate risk factor or when the time-to-event trait is associated with more than one intermediate risk factor. This observation also illustrates one of the limitations of the joint model for one longitudinal with one time-to-event trait, with the consequence that an indirect SNP association can be mistaken as a direct association when other longitudinal risk factors are omitted in the joint model.

Generalization of the joint model to multiple longitudinal and multiple time-to-event traits

To characterize the genetic architecture of a system of multiple longitudinal risk factors and multiple time-to-events, we propose an extension of the joint model to L longitudinal and K time-to-event traits (L>1, K>1), as shown in Fig. 2 and detailed as follows. We define y_i,l = (y_i,l,1, …, y_i,l,j, …, y_i,l,J) as the observed longitudinal measures for each l^th QT, 1 ≤ l ≤ L, collected over the J visit times t_i. We define (T_i,k, δ_i,k) as the vector of observed right-censored event time T_i,k and event indicator δ_i,k for each k^th time-to-event trait for individual i, with . We assume the same censoring time C_i across all K outcomes, but the model can be extended to situations where C_i varies for each time-to-event trait. Again, for ease of presentation, we simplify the model notation with no adjusting covariates and assume linear trajectories for all L longitudinal traits and contemporaneous effects of L longitudinal traits on the K time-to-event traits. The model can be extended to account for non-linear trajectories, cumulative longitudinal effects, and trait-specific and/or shared covariates such as confounding factors or ancestry-related principal components.

• Download figure
• Open in new tab

Fig. 2. Proposed joint modelling approach for characterization of complex genetic architecture of multiple disease progression

Multivariate longitudinal sub-model

In the extension of the longitudinal sub-model to L longitudinal traits, we index subscripts in Equations 1 and 2 for each l^th longitudinal trait (Equations 4 and 5 in Fig 2). The vector of observed trait values becomes: Where:

• X_i,l = (1_J, t_i, g_i1_J) and Z_i,l = (1_J, t_i) are the design matrices for fixed and random effects,

• β_l = (β_0,l, β_1,l, β_g,l)^T and b_i,l = (b_i,0,l, b_i,1,l)^T denote the QT-specific fixed and random effects,

• ε_i,l = (ε_i,1,l, …, ε_i,J,l, …, ε_i,J,l)^T is the vector of residual error terms, with ε_i,l∼N_J(0, Σ_l), where with , the residual variance for the l^th QT; we assume the ε_i,l are independent for all L traits.

To account for dependence among the L longitudinal QTs, we assume the overall random effects vector for all L QTs, b_i = (b_i,1, …, b_i,l, …, b_i,L)^T∼N_2L(0, D), where is the variance-covariance matrix for all L QTs, accounting for serial dependencies within each l^th QT, i.e. , and accounting for cross-dependencies between each pair l, m of QTs with l≠m, that is .

This formulation implies that the vector of stacked repeated measures for all L longitudinal QTs for individual i, y_i = (y_i,1, …, y_i,l, …, y_i,L)^T in ℜ^J×L follows a multivariate normal distribution with mean E[y_i] = X_iβ and variance , where:

• X_i = diag(X_i,1,…, X_i,l, …, X_i,L) and Z_i = diag(Z_i,1,…, Z_i,l, …, Z_i,L) are the overall (JL-by-3L) and (JL-by-2L) design block diagonal matrices for the fixed and random effects respectively;

• β = (β₁, …, β_l, …, β_L)^T is the 3L-length stacked vector of fixed effects,

• D is the (2L- by-2L) covariance matrix for random effects b = (b_i,1, …, b_i,l, …, b_i,L)^T,

• Σ = diag(Σ₁, …, Σ_l, …, Σ_L)^T is the (JL-by-JL) block diagonal matrix of residual variances.

Assuming, for each l^th QT, an unstructured variance-covariance matrix D_l,l, the variance at each visit time t_i,j is Var(y_i,j,l) and the covariance function Cov(y_i,j,l, y_i,s,l) between two visit times t_i,j ≠ t_i,s are analogous to those defined above for the joint model with a single longitudinal and single time-to-event trait. The multivariate mixed model accounts for dependencies between each QT pair l ≠ m via random-effect covariance functions in the D_l,m matrices where the covariance between observations of two QTs (l, m; l ≠ m) measured at the same visit time t_i,j is which is quadratic over time, and the covariance function between two longitudinal QTs (l ≠ m) measured at different visit times . Thus, joint analysis of correlated longitudinal QTs is expected to improve power over separate analysis of each QT by borrowing information through implied dependency structures among the random effects.

Multivariate time-to-event sub-model

Finally, we extend Equation 3 to a multivariate PH frailty time-to-event sub-model, with a subject-specific random effect (frailty term, u_i) introduced to capture potential unexplained dependencies (e.g. due to unmeasured baseline shared factors) among the time-to-event traits. In Equation 6 (Fig 2), λ_0,k(t) and γ_g,k, correspond to the baseline hazard function, and SNP effect on the k^th time-to-event trait (1 ≤ k ≤ K), accounting for association of each l^th QT with the time-to-event trait k (α_l,k, 1 ≤ l ≤ L). We assume that the subject-specific frailty term follows a gamma distribution, that is u_i∼Γ(a, b) with a, b>0, and represents dependencies among the K traits. Equation 6 (Fig 2) can be expressed as: where α_k = (α_1,k, …, α_l,k, …., α_L,k)^T is the vector of all L QT effects on the k^th time-to-event trait, and w_i,k(t) = (w_i,1,k(t), …., w_i,l,k(t), …, w_i,L,k(t)) specifies the corresponding association profile of each QT with the k^th time-to-event trait. We note , where , denotes the l^th QT trajectory (1 ≤ l ≤ L) at time t, which depends on the fixed and random effects β_l and b_i,l.

Comparisons with joint model of one longitudinal and one time-to-event trait

In the proposed joint model extension for multiple longitudinal and multiple time-to-event traits, the direct, indirect and overall SNP effects defined above for the joint model with one longitudinal and one time-to-event trait are interpreted similarly. However, there are important practical differences between the latter model for a pair of traits and the proposed multi-trait. First, because the joint model extension can account for multiple intermediate longitudinal QT risk factors associated with one (or multiple) time-to-event trait(s), it improves inference for SNP association and accuracy of SNP classification, particularly when a time-to-event trait depends on more than one longitudinal risk factor as illustrated in our numerical experiments. Second, in the multivariate longitudinal sub-model, the variance-covariance matrix D for the random effects specifies non-null covariance terms in D_l,m for each pair of longitudinal QT l, m (1 ≤ l ≤ L and 1 ≤ m ≤ L, l ≠ m). In contrast, under the assumption of null covariance terms in D_l,m for all QT pairs inherent in separate analyses of each QT, the multivariate sub-model reduces to independent sub-models for each longitudinal QT. When longitudinal QTs are correlated, assuming null covariances can fail to make use of information borrowed through the random effects and reduce efficiency of the parameter estimates in the longitudinal trajectories (Shah et al. 1997; Jensen and Ritz 2018). This, in turn, can affect estimation in the time-to-event model. Third, without a frailty term u_i, the time-to-event sub-model (Equations 6 and 7) reduces to separate time-to-event sub-models for each time-to-event trait. Thus, through use of a shared frailty term, the extended joint model accounts for residual dependency between the K time-to-event traits, not explained by the covariates shared by the time-to-event sub-models. Overall, the proposed joint model for multiple longitudinal and multiple time-to-event traits can improve inference by accounting for intermediate longitudinal QT(s) and their dependencies, as well as dependencies among the time-to-event traits, and thereby improve classification accuracy of direct and/or indirect SNP associations.

Effect estimation and test statistic construction

To address computational obstacles involved in the maximization of the joint likelihood and allow more flexible inference, we estimate the parameters using a two-stage approach (see details in the Appendix). We work within the framework originally defined by (Tsiatis et al. 1995; Wulfsohn and Tsiatis 1997; Dafni and Tsiatis 1998; Tsiatis and Davidian 2001) and in the spirit of subsequent authors (Ye et al. 2008; Yuen et al. 2018; Arisido et al. 2019); (Tsiatis and Davidian 2001; Tsiatis and Davidian 2004) specify conditions that guarantee the two-stage estimators are consistent and asymptotically normal. Specifically, in Stage 1 we fit a multivariate mixed model (Equation 5) using the mvlme() function from the R package JoineRML [(Hickey et al. 2018c), version 0.4.2] to estimate the parameters of the longitudinal trajectories of the risk factors, and obtain fitted values of the smoothed trajectories. In Stage 2, we fit a Cox PH frailty time-to-event model (Equations 6 or 7) adjusting for functions of the smoothed trajectories as time-dependent covariates using the coxph() function from the R package survival [(Therneau and Grambsch 2000; Therneau 2020), versions 3.2.7 and 3.2.13]. We assume a Gamma distribution for the frailty term (u_i) and a separate non-parametric baseline hazard function for each time-to-event trait using the strata argument in coxph().

To account for propagation of errors, due to uncertainty in Stage 1 estimates that is not accounted for during Stage 2 parameter estimation (Wulfsohn and Tsiatis 1997), and to empirically estimate the joint covariance matrix of SNP-QT trajectory (β_g,l) and SNP-TTE effects (γ_g,k), we apply a nonparametric bootstrap. The bootstrap also provides reliable standard error estimates in the joint time-to-event sub-model needed when using an unspecified baseline hazard (Hsieh et al. 2006; Lawrence Gould et al. 2015; Furgal et al. 2019). For each bootstrap sample b (1 ≤ b ≤ B, B is the total number of bootstrap repetitions), we generate a new dataset by randomly sampling N individuals with replacement and refitting the joint model on each new dataset b. We compute the empirical joint covariance matrix for all β_g,l, γ_g,k and α_l,k parameter estimates using the B bootstrap parameter vector estimates. Wald statistics for each β_g,l are computed as using the empirical bootstrap standard errors se_βg,l (to test H0: β_g,l = 0 vs H1: β_g,l ≠ 0), and similarly for each γ_g,k as (to test H0: γ_g,k = 0 vs H1: γ_g,k ≠ 0).

In contrast to the two-stage joint model, a conventional one-stage analysis to assess whether a SNP has an association with the time-to-event trait, independent from the QT-TTE association, relies on regression adjustment using observed longitudinal QT values as time-dependent covariates in a Cox-PH model (Paterson et al. 2010; Deng and Pan 2017). This approach, based only on the time-to-event model, does not provide information about SNP-QT association (β_g,l) and interprets the SNP as having a direct association with the time-to-event if the test of SNP-TTE effect (γ_g,k) is declared significant, given the observed QT. Inference for α_l,k under this approach can be biased or inefficient when the QT is measured with random error or high within-subject variability (Faucett and Thomas 1996; Wulfsohn and Tsiatis 1997; Xu and Zeger 2001; Song et al. 2002; Brown and Ibrahim 2003), and inference for γ_g,k may also be affected. Although estimates of the SNP-QT effect (β_g,l) obtained from mixed model QT analysis, fitted separately, may be used to distinguish between direct alone versus both direct and indirect SNP association, unlike the joint model, this two-step conditional approach ignores measurement error in the observed QT values.

Procedure to classify direct and/or indirect SNP associations

In Table 1, we present a practical procedure to classify a SNP as having direct and/or indirect association with a time-to-event trait k, that accounts for the SNP association with a longitudinal risk factor l. This procedure requires two significance thresholds, and for hypothesis tests of each of β_g,l and γ_g,k respectively, to be specified prior to the analysis and adjusted for the number of SNPs tested. Depending on the research question, we can choose different values for and , or the same value . The latter is applicable, for instance, to systematically classify direct and/or indirect association for a set of M SNPs, and the former to assess which SNPs, among those reported to be associated with the longitudinal risk factor, have a direct effect on a time-to-event trait. To our knowledge, no comparable procedure to classify direct and/or indirect SNP association, based on SNP effect estimates from joint models, has been proposed for studies with longitudinal risk factors and time-to-event traits. A key feature of the proposed joint model extension to multiple longitudinal and multiple time-to-event traits is inference for SNP effects on each of the traits in a single integrated statistical model, while accounting for within-subject QT variability and dependencies among the traits. The focus of the simulation study and DCCT data application, which follow, is to evaluate the SNP classification procedure applied in extended joint model analysis.

Algorithm for realistic data generation under a complex genetic architecture

To generate a data structure that combines observed and simulated traits, we formulate a genotype-phenotype multi-trait model including: (i) L=3 linear mixed models linking each SNP with an indirect effect to a longitudinal risk factor, and (ii) K=2 non-independent parametric time-to-event models depending on fitted longitudinal QT trajectories and SNPs with direct effects. For each DCCT individual i with observed longitudinal measures for HbA1c and SBP, and observed baseline covariates (T1D duration, sex), we simulate: genotypes at M causal SNPs with MAF vector p, longitudinal trait values U_i, and time-to-event traits ((T_i,k, δ_i,k), k=1,2 for DR and DN), using the algorithm illustrated in Fig. 4 and detailed in File S2 (sections 1-5). All SNP genotypes are generated under Hardy-Weinberg and linkage equilibrium assumptions. SNPs with indirect effects through longitudinal QTs associated with DR and/or DN are generated from the observed (SNP1, SNP5) or simulated (SNP3) QTs, while SNPs with direct effects (SNP2, SNP4) are generated independently of the longitudinal QTs and are included in the hazard function used to generate each time-to-event trait (Fig. 4).

• Download figure
• Open in new tab

Fig. 4. Illustration of the procedure developed for DCCT-based simulation study under the scenario from Fig. 3

For each individual i, with {t_i, y_i,1, y_i,2} observed in DCCT, the algorithm simulate: latent longitudinal QT values for U (y_i,u), genetic data at M causal SNPs (with p, the specified vector of MAFs) and time-to-T1DC data (T_i,k, δ_i,k) for K=2 time-to-T1DC traits. The genetic data are simulated under Hardy-Weinberg and linkage equilibrium assumptions. SNPs with indirect effects are simulated from a multinomial distribution with calculated conditional genotype probabilities for individual i (π_gi=0, π_gi=1, π_gi=2) based on y_i,l. Each y_i,l is assumed to follow a multivariate normal distribution with X_i,l and Z_i,l the specified fixed and random effect design matrices in longitudinal trait models and the vector of specified parameter values for each l^th QT. SNPs with direct effects are simulated from the population probabilities, that depend only on the MAF. The specified hazard function for each k^th time-to-event trait depends on the effects of the longitudinal QT trajectories and on the SNPs with direct effects in η_i,k, with for DR and for DN, as well as the effect of the shared latent QT trajectory used to induce some dependencies between the time-to-event traits. We define Γ_k as the vectors of specified parameter values for each k^th time-to-event trait. The uncensored event time , is simulated by calculating the inverse of the cumulative specified hazard function using the Brent univariate root-finding method (Brent 2013; Crowther and Lambert 2013). To simplify the exposition of the simulation procedure, we ignore the effects of the sex on SBP and of the T1D duration on both T1DC traits, but they were included in the data generating model, see File S2 for details. Parameters for the causal genetic scenario are shown in Fig. 3 and File S2 (section 4)

Scenario for DCCT-based complex genetic architecture

Overall, the complex genetic architecture represents multiple types of SNP-trait associations (Fig. 3): direct association with each T1DC trait (SNP2, SNP4), indirect association with both T1DC traits via measured (SNP1) and unmeasured (SNP3) longitudinal QTs; and direct and indirect association via a measured longitudinal QT (SNP5); all longitudinal risk factors exhibit within-subject random variability. Except for SNP3, all other SNP scenarios represent SNP association with a longitudinal QT risk factor (SNP1, SNP5) or a time-to-event trait (SNP2, SNP4, SNP5) testable in a single-trait GWAS. SNP1, SNP3 and SNP5 have indirect effects on T1DC traits, such that their associations with the T1DC traits are detectable in discovery analysis for each TTE trait, using Cox PH time-to-event models fitted separately for each TTE trait and ignoring the longitudinal QT risk factors (File S2, section 8). SNP1 corresponds roughly to rs10810632 and rs1358030 associations reported in the motivating DCCT GWAS of HbA1c (Paterson et al. 2010), while SNP5 represents a strong signal that would be detected in separate GWAS analysis of each longitudinal and time-to-event trait.

Analysis of the simulated data

To evaluate the statistical performance under the complex genetic architecture outlined above, focusing on direct and/or indirect classification, we compare extended joint model analysis of multiple longitudinal QT risk factors and multiple T1DC with alternative analysis methods that do not fully exploit the data structure. These alternative approaches include joint model analyses limited to two longitudinal QT and one time-to-event trait, joint model analyses of one longitudinal QT and one time-to-event trait, and a conditional time-to-event analysis of two time-to-event traits adjusted for observed values of two longitudinal QTs. We assess the impact of model specification on the classification accuracy of direct and/or indirect SNP associations by fitting mis-specified joint models that leave out important covariates or traits. Altogether, the comparisons are designed to assess the merits of extended joint model analysis over simpler available methods. In each replicated dataset, each of the five SNPs is separately analysed for association with the longitudinal QTs and TTE traits, using four alternative analytic approaches:

• JM-cmp: a completely specified joint model analysis that includes observed (HbA1c, SBP) and unobserved (U) longitudinal QT as well as baseline covariates (sex, T1D duration) used in the data simulation. Due to the latent nature of unobserved U, JM-cmp cannot be fitted in practice, but we include it as a benchmark for comparison against the data analysis models JM-mis that are fitted without U.

• JM-mis: includes the same variables as in JM-cmp, but excludes U, the unobserved longitudinal QT.

• JM-sep: joint models of two longitudinal traits and one time-to-event trait that do not account for dependency between the time-to-event traits (where JM-sep (l = 1,2; k=1) denotes the joint model for DR; and JM-sep(l = 1,2; k=2) the joint model for DN), and joint models of one longitudinal and one time-to-event trait that do not account for dependence between the longitudinal traits, nor between the time-to-event traits (referred to as JM-sep(l = 1; k=1), JM-sep(l = 1; k=2), and JM-sep(l = 2; k=2)). Altogether, the JM-sep models assess the merits of the extended joint model methods in comparison to JM-mis and JM-cmp.

• CM-obs: a Cox PH frailty survival analysis of both time-to-event traits (DR, DN) that includes the same variables as in JM-mis but adjusts for the observed longitudinal QT values as time-dependent covariates; this model corresponds to the conditional analysis approach mentioned above. Here, to classify SNPs as indirect association, we fit a linear mixed model for both QTs to test the SNP effects on the QT(s) adjusted for the same covariates as used for the joint models. Comparisons of estimation, hypothesis testing, and classification results based on CM-obs to those based on JM-mis and JM-cmp allow us to assess the impact of within-subject QT variation/measurement errors on hypothesis testing and classification results for each SNP.

For each of these analyses, we compute empirical covariance matrices for the effect estimates using 500 bootstrap iterations and construct large sample test statistics for each of the SNP effect parameters. Under two-stage JM inference, in stage 1 we fit a bivariate (l=1,2) or univariate (l=1 or l=2) longitudinal QT model and test the β_g,l parameters for indirect SNP association. Then in stage 2, we fit the time-to-event models with the QT trajectories from stage 1 according to the TTE model specification (k=1,2; k=1, or k=2), and test the γ_g,k parameters for direct SNP association. Because stage 1 QT analysis is shared among methods, differences in classification among analytic methods largely arise through differences in the test results in stage 2 TTE analysis.

Under two-stage and conditional independence assumptions described in an Appendix, the SNP association estimates of β_g,l and γ_g,k and corresponding test statistics are expected to be uncorrelated under the global null scenario, but this may not necessarily hold under the genetic alternative scenario when the analysis model is mis-specified or when the time-to-event estimation uses observed longitudinal trait values. We therefore compute empirical correlations in each replicate under both scenarios.

Given hypothesis test results for a pair of QT/TTE traits for each SNP, we apply the procedure defined in Table 1 to classify the SNP-TTE association.

Evaluation criteria

We compare type I error and power of hypothesis tests of SNP-QT (β_g,l) and SNP-TTE (γ_g,k) association among the alternative analytic approaches (JM-cmp, JM-mis, JM-sep, CM-obs) under the global null and causal genetic scenarios for each of the 5 SNPs analysed separately, at P* = 5% and 1% critical values. We assess estimation accuracy using mean bias for β_g,l and γ_g,k estimates, and confidence interval coverage across replicates, and similarly examine the distribution and mean of their bootstrap standard errors and correlation for all the compared models.

For each of the 5 SNPs, we evaluate accuracy of the procedure presented in Table 1 to classify SNP association with each of the TTE traits as direct and/or indirect. Specifically, under the global null and causal genetic scenarios, we compare the empirical classification frequencies to the expected classification frequencies under the assumption of indirect and/or direct association built into the generating model. The empirical frequencies are tabulated from the distribution of simulation replicates in the four classification categories (direct, indirect, direct and indirect, not direct and not indirect) as defined in Table 1, using specified classification thresholds such that , with P;^*=0.05. We calculate expected frequencies under the assumption that the Z statistics, constructed from the estimates of β_g,l and γ_g,k, and their bootstrap variances, follow a bivariate normal distribution with correlation specified by the bootstrap estimates (see File S2, section 9 for details), to allow for potential dependence between the SNP parameter tests. We judge the classification procedure for a SNP association with each QT-TTE trait pair to have high accuracy when the empirical frequencies are consistent with those expected, and we compare accuracy among the different models. We also assess empirical classification frequencies for SNPs with non-null effects across variation in stringency of significance thresholds up to P;^*=10⁻⁵.

SNP association test validity and power

Under the global null scenario of no genetic association with any of the traits, the type I error of each SNP test (β_g,l and γ_g,k) is reasonably well controlled (Table 2, left-hand side), and the P-values P;_βg and P;_γg from the joint models show no departure from the expected large sample distributions (χ² with 1 degree of freedom, see section 7 in File S2).

Under the alternative multi-SNP causal genetic scenario (Table 2, right-hand side), type I errors for tests of each null β_g,l tend to be close to the nominal level of 5% for most analysis models (with exceptions for SNP2 and SNP5), while tests of the two SNPs with effects on intermediate measured longitudinal QT risk factors (SNP1 (β_g,1) and SNP5 (β_g,2)) reach 100% power for all the analysis models (Table 2, right-hand side). Tests of null γ_g,k by JM-cmp, JM-mis and JM-sep(l=1,2; k, with k=1 or 2) show overall good type I error control for all SNPs with direct and/or indirect effects via measured longitudinal QTs (i.e. all SNPs but SNP3). However, for SNP1 with indirect effects on both T1DC traits via HbA1c, tests of null γ_g,2 exhibit liberal type I errors in JM-sep(l=2; k=2), which may be explained by bias in SNP1 effect on DN (, see File S2, section 6) towards the indirect SNP1 effect via HbA1c QT risk factor, which is ignored from JM-sep(l=2; k=2). For SNP2 which has direct effect on DR (γ_g,1), tests of non-null γ_g,1 exhibit equivalent or higher power under JM-cmp and JM-mis (closely followed by JM-sep(l=1,2; k)). However, the power to detect the direct SNP4 effect on DN (γ_g,2) appears reduced in JM-sep(l=1; k=2); which can be explained by the bias in towards the null, due to ignoring SBP QT risk factor in JM-sep(l=1; k=2). For SNP5 with direct and indirect effects on DN via SBP, γ_g,2 appears detected by a larger number of replicates under JM-sep(l=1; k=2) and CM-obs (Table 2, right-hand side); which may be explained by bias in SNP5 effect (, see File S2, section 6) towards the overall SNP5 effect (combining direct and indirect SNP5 effects via SBP), given that SBP is ignored from JM-sep(l=1; k=2) and random SBP variation in CM-obs. Finally, for SNP3 that has indirect effects on both T1DC traits via the unmeasured longitudinal risk factor U, only JM-cmp accounts for indirect SNP3 pathways via U. This translates into elevated type I errors for tests of each null γ_g,k likely explained by the biased SNP3 effects on both T1DC traits (, see File S2, section 6) in JM-mis, all JM-sep, and CM-obs, towards the indirect SNP3 effect via U. Overall, the relative ranking of empirical powers among the analysis models persists for tests of each γ_g,k across P* varying up to 10^-5 (File S2, section 8), with steeper power reduction for detection of SNP4 effect on DN (γ_g,2) in JM-sep(l=1; k=2), and markedly misleading detection of SNP5 effect on DN (γ_g,2) in JM-sep(l=1; k=2) and CM-obs.

As noted above, the improvement of the type I errors control and powers in tests of each γ_g,k under the alternative simulation scenario by JM-cmp (closely followed by JM-mis and JM-sep(l=1,2; k=1 or 2)) compared to JM-sep(l=1 or 2; k=1 or 2) for one longitudinal and one time-to-event trait and CM-obs can be explained by more efficient estimation accuracy in , but also in the QT effect estimates on the TTE as illustrated by: bias reduction, coverage probabilities closer to the nominal level (95%), and empirical 95% confidence intervals of the parameter estimates narrower around the true parameter (File S2, section 6). Moreover, we find little evidence for correlation between and . The average bootstrap correlation across replicates is low for each QT/TTE trait pair (l; k) under the alternative genetic (Tables 3-7) and global null scenarios (File S2, section 9) for all joint model analyses . However, we see larger values in CM-obs, particularly for the SBP/DN trait pair; this may be explained by larger random variation in SBP which is ignored by CM-obs.

Joint model fitting

We fit the joint model for each SNP one at a time, including baseline covariates (age at diagnosis, T1D duration, cohort, sex, and year of DCCT study entry), longitudinal HbA1c and SBP, and T1DCs DR and DN (File S3, section 3). Given prior evidence for long term HbA1c effects on T1DC, we present results under the time-weighted cumulative for HbA1c association with time-to-T1DC traits, which exhibits stronger prior association with T1DC in the DCCT individuals analyzed here, but we obtained similar classification results under the two alternative specifications considered for HbA1c association (i.e. contemporaneous value and updated mean, see File S3, section 4). We apply diagnostic tools in the joint model analysis, including residual analysis in the longitudinal and time-to-event sub-models, finding little evidence for model misspecification (File S3, section 5). In the Cox PH frailty time-to-event sub-model, PH assumptions are well-satisfied when the baseline hazard is stratified on the cohort variable; overall conclusions are equivalent when cohort is included as a covariate. As shown in File S3 (section 5), the martingale residuals of the Cox PH frailty model are consistent with the assumption of a linear relationship between the QTs and each time-to-T1DC outcome.

In Fig. 5 (Panel A), we show that rs1358030 and rs10810632, which were discovered in a previous GWAS of HbA1c in DCCT (Paterson et al. 2010), are classified as indirect associations with both T1DC traits via their association with the HbA1c shared risk factor (P;_βg≤ P;^* and P;_γ > P;^*, P;^*=1.7×10^-4 after Bonferroni correction for the 289.02 effective SNPs tested (Li and Ji 2005)). Although the other candidate SNPs were selected from larger meta-analysis in T1D and/or independent studies in general populations (Paterson et al. 2010; Grassi et al. 2011; Sandholm et al. 2012; Hosseini et al. 2015; Wheeler et al. 2017; Evangelou et al. 2018; Pollack et al. 2019), the SNP-trait association tests did not reach Bonferroni-corrected significance thresholds in our analysis (File S3, section 4), potentially due to effect heterogeneity and low power to detect these variants or by some inherent heterogeneity factors in the study designs and/or phenotypic definitions. Thus, classification for these candidate SNPs is uninformative.

• Download figure
• Open in new tab

Fig. 5. Classification of direct and/or indirect SNP associations in the DCCT Genetics Study data.

(A) Scatter plots of the P-values (-log10) for tests of β_g,l (H0: β_g,l = 0 vs H1 β_g,l ≠ 0) on the X axis and γ_g,k (H0: γ_g,k = 0 vs H1: γ_g,k ≠ 0) on the Y axis for HbA1c/DR, HbA1c/DN and SBP/DN trait pairs. Significance levels P;^* =1.7×10^-4 and P;^* = 0.05 are indicated by red and grey horizontal and vertical dashed lines. (B) and (C) represent association results for rs10810632 and rs1358030 detected as indirect associations at P;^* =1.7×10^-4. Left panels present results from separate analysis of each trait (ie. longitudinal model for each QT and Cox PH time-to-event model without adjusting for the longitudinal traits) as used in naïve discovery GWAS; and right panels show results from the joint model with bootstrap 95% confidence intervals for the direct and indirect SNP effects. Results are presented using time-weighted cumulative HbA1c effects on T1DC traits.

Sample size and power considerations

The DCCT-based simulation study shows that multiple-trait joint model analysis can achieve high classification rates for direct and/or indirect SNP associations under SNP effect sizes in the specified range (SNP effects: on the QTs 0.7-7, on T1DC 0.4-0.8; and MAFs 10-40%) and classification rates under the global null hypothesis (null SNP effects on all traits) controlled to the nominal levels for direct or indirect SNP association but conservative for direct and indirect association. Compared to the simulation, the number of DN events is lower in the DCCT data application. We thus expect lower power to detect direct SNP association with DN; this implies reduced accuracy to classify a SNP as having either a direct, or a direct and indirect association.

In light of the simulation results, we acknowledge that a larger sample size would be required to maintain performance in DCCT for genetic associations with lower MAFs/effect sizes than the ones specified. In particular, lower effect sizes and MAFs can lead to larger variances of estimated effects and larger 95% confidence intervals, as observed for the direct effects of rs10810632 (MAF=7%) compared to rs1358030 (MAF=36%) for both time-to-T1DC traits (Fig. 5, Panels B and C). SNP effects with larger variances have lower power to detect association. Although we expect the effect sizes estimated by the joint model to be relatively insensitive to the study sample size, reduction in variances by increasing study sample size would be expected to improve classification accuracy for these two SNPs (particularly for rs10810632 that has a lower MAF).

To assess how genetic association and classification results for rs10810632 and rs1358030 may be affected by increasing sample size, we applied parametric resampling (see File S3, section 6 for details) to draw datasets with sample size up to five times the DCCT sample of N=516, and we then extrapolated the classification beyond N=2580. This analysis demonstrates decreasing variances of the SNP effects for both SNPs as sample size increases and narrowing of the 95% and 99% confidence intervals (File S3, section 6); panel A in Fig. 6 illustrates the corresponding shift in classification of rs10810632 association with DR as indirect via HbA1c towards classification as both indirect and direct. On the other hand, given that both SNPs were previously discovered in a GWAS of HbA1c in DCCT, SNP effect estimates for HbA1c association may be overestimated due to Winner’s curse bias (Kraft 2008; Sun et al. 2011), which would impact the classification results. We therefore repeated the sample size analyses specifying an adjusted SNP effect size for the HbA1c association equal to 50% of its effect estimate in our DCCT analysis. In this scenario the test of SNP association with HbA1c falls just below the P* threshold in the resample size N=516, although power improves in larger sample sizes (Fig. 6). Here again classification tends to shift from indirect towards both direct and indirect SNP associations with both T1DC traits (see File S3, section 6 for complete results for both SNPs); however, a much larger sample size is needed to cross the P* threshold for test of direct association.

• Download figure
• Open in new tab

Fig. 6. Change in classification results for rs10810632 with HbA1c/DR to increasing sample size (A) and to the Winner’s curse bias for the SNP effect on HbA1c (B) investigated using parametric resampling.

We used parametric resampling (File S3, section 6) to draw datasets with sample size up to five times the DCCT sample size of N=516, and then extrapolated the classification results beyond N=2580. The X axes of (A) and (B), show the P-values (-log10) for the test of rs10810632 effect on HbA1c (H0: β_g,1 = 0 vs H1 β_g,1 ≠ 0)), while the Y axes show the P-values (-log10) for the test of rs10810632 effect on DR (H0: γ_g,1 = 0 vs H1: γ_g,1 ≠ 0) for the sample sizes investigated (shown by different colors). We fitted a regression line on each plot to project the trend of the classification beyond N=2580 individuals. These plots illustrate the corresponding shift in classification of rs10810632 association with DR as indirect via HbA1c towards classification as both indirect and direct association. Complete results with HbA1c/DN are shown in File S3 (section 6).

In the DCCT application, we take advantage of the largest available clinical study of T1D complications with long-term follow-up and high-density longitudinal QT measurements. This highlights the dearth of longitudinal studies with both a large number of individuals and long-term clinical follow-up, as well as the related challenges in joint model analysis. In prospective study designs with both longitudinal and time-to-event traits, an inherent imbalance can exist among traits in detection of SNP associations, in that power for QT(s) depends on the number of measurements while power for time-to-event traits depends on the number of events and duration of follow-up. Although there is currently a trade-off between epidemiological studies with large sample sizes but low density of longitudinal measurements, and clinical studies with more modest sample sizes but high measurement densities, we anticipate informative application of joint model analysis in large biobanks, for example (Scholtens et al. 2015; Bycroft et al. 2018; Dummer et al. 2018), now assembling longitudinal measures jointly with binary outcomes and genetic data.