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

Introduction

Despite their known ability to improve inference in clinical and epidemiological studies, particularly in the presence of informative censoring/dropout or when longitudinal traits are measured with error^1–3, joint models for longitudinal and time-to-event outcomes have received limited attention in genetic association studies. Genome-wide association studies (GWAS) of intermediate quantitative traits (QTs) typically require follow-up studies to identify whether SNP associations detected with each of the QT(s), analyzed separately, also affect related outcomes through direct and/or indirect effects induced by those traits. Such intermediate traits may include established risk factors for related clinical outcomes and be measured with errors (e.g. biological variation). By accounting for QT measurement error and dependencies among traits, joint models can improve accuracy and efficiency of effect estimation, as well as detection of SNP associations.

The so-called shared-random-effects joint model^4,5 consists of a sub-model for a single longitudinal trait and a sub-model for a single right-censored time-to-event trait. The longitudinal sub-model describes the QT as an underlying smoothed trajectory that depends on fixed effects of time and baseline covariates, as well as subject-specific random effects. The joint model association structure is induced via the functional dependence between the hazard of an event at any time t and the longitudinal trait trajectory^6,7. This relationship can be based on prior biological knowledge of the link between that longitudinal and time-to-event traits. As elucidated by Ibrahim et al¹, this class of joint models lends itself to interpretations of direct/indirect effects because the relationship between a baseline covariate, such as a SNP genotype, and each of the longitudinal and time-to-event traits, as well as the relationship between the longitudinal and time-to-event traits can be specified via model parameters corresponding to direct, indirect and overall effects. Extensions of joint models for multiple longitudinal and/or multiple time-to-event traits can further improve inference by borrowing information among related traits. Extensions of joint models have been reviewed for multiple longitudinal traits^6,7 or multiple time-to-event traits⁸. A few extensions for both multiple longitudinal and multiple time-to-event traits have been developed (^9–11, among others), but these models are often formulated for a specific study question, and thus can lack generalizability. Such extensions raise computational challenges for maximisation of the likelihood under the multivariate random effect’s distribution. Two-stage approaches^12–15 for joint model fitting appear more computationally efficient and lends itself to more flexible model formulation but inference can be mis-calibrated because parameter estimates and predictions from Stage 1 are obtained from the longitudinal model without consideration of the time-to-event outcome, and the uncertainty in Stage 1 estimates is ignored during Stage 2 estimation, a problem known as propagation of errors¹⁶.

Motivated by the complex genetic architecture of long-term type 1 diabetes complications (T1DC), we propose joint-model extension to evaluate genetic associations with multiple time-to-event traits and multiple longitudinal risk factors. Risk of development of T1DC, including diabetic retinopathy (DR) and diabetic nephropathy (DN), is hypothesized to result from multiple genetic factors with potential direct and/or indirect effects via multiple - shared and/or specific - risk factors¹⁷. In addition to genetic factors, hyperglycemia (measured by Hemoglobin A1c, referred as HbA1c) represents the major QT risk factor for T1DC. Other measured longitudinal traits, that can also be influenced by genetic factors, have postulated associations with T1DC, for example, association of systolic blood pressure (SBP) with DN. Because conventional analysis methods in GWAS ignores dependency between related traits, SNPs can appear associated with multiple traits (including SNPs with direct and/or indirect effects induced by intermediate QT(s) measured or unmeasured). Accurately distinguishing between direct and/or indirect associations is crucial to elucidate the biological pathways through which genetic factors operate into the etiology of those complex traits, like in our motivating example of the T1DC genetic architecture. A naïve approach to disentangle direct from indirect associations can be achieved by including the intermediate QT as a time-dependent covariate in a Cox Proportional Hazard (PH) model. However, measurement errors in intermediate QT(s), presence of informative censoring in longitudinal trait, unmeasured shared risk factors between longitudinal/time-to-event traits can challenge accurate classification of direct and/or indirect SNP associations.

The major contribution of our work is a general formulation for the joint analysis of multiple time-to-event traits and multiple longitudinal QT risk factors in genetic association studies. We develop inference methods for statistical genetic analysis based on joint model parameters estimation, including a hypothesis testing procedure to classify direct and/or indirect SNP associations with each time-to-event trait. A second contribution of the paper is the implementation of a data-informed simulation algorithm to generate multiple causal SNPs with various direct effects on simulated time-to-event traits and/or indirect effects via observed (measured) longitudinal QTs in the Diabetes Control and Complications Trial (DCCT) study^18,19 and unobserved (simulated) longitudinal QTs. This algorithm also provides a general approach to estimate achievable power given sample size. Our simulation results show bias reduction and improved classification of direct and/or indirect SNP associations in comparison with alternative approaches ignoring measurement errors in longitudinal QT(s). Thirdly, we apply the joint model in DCCT data and clarifies two SNPs as having indirect associations via the HbA1c longitudinal risk factor, with consistent conclusions using alternative time-dependent association structures that account for established cumulative and time-weighted effects of HbA1c on T1DC traits^20,21. Example R codes for data simulation and for application of the proposed joint model are available on GitHub.

Methods

Joint Modelling Approach

Model Formulation

We assume that a set of M SNPs are available with K (1 ≤ k ≤ K) unordered and non-competing time-to-event traits, such as multiple disease complications, and L (1 ≤ l ≤ L) longitudinal QT traits (i.e. intermediate risk factors) possibly measured with error (i.e. biological variation) in N unrelated individuals indexed by i (1 ≤ i ≤ N). For each SNPm (1 ≤ m ≤ M), with minor allele frequency (MAF) p_m, the genotypes are coded as the number of copies of the minor allele, under an additive genetic model. For each individual denotes the vector of quantitative measures collected over scheduled visits at t_ij for j=1,…, n_i for the l^th longitudinal trait (1≤l≤ L) with n_i represents the maximum number of visits recorded. Let (T_i(k), δ_i(k)) be the vector of observed right-censored event time T_i(k) and event indicator δ_i(k) for the k^th time-to-event trait.

We assume , where is the latent (uncensored) time-to-event k and C_i is the censoring time (e.g. administrative censoring). We define , with δ_i(k) = 1 if the event occurs during the observation period and δ_i(k) = 0 otherwise.

To characterize the genetic architecture of a system of multiple time-to-events and multiple longitudinal risk factors, we formulate a shared-random-effects joint model, as specified in Figure 1, combining: (i) longitudinal and (ii) time-to-event sub-models connected by (iii) specified time-dependent association structures.

• Download figure
• Open in new tab

Figure 1. Proposed joint modelling approach for characterization of complex genetic architecture of multiple disease progression

(i) Longitudinal sub-model

The longitudinal sub-model is specified by a multivariate mixed model for the L longitudinal QTs (or risk factors), based on the multivariate extension of the Laird and Ware linear mixed model²². It describes each vector of observed longitudinal QT measures Y_i(l)(t_ij) as noisy observations of a true and unobserved smoothed subject-specific longitudinal trajectory, , with noise terms . Each smooth trajectory describes subject-specific evolution and depends on time, SNP effect, individual-level random effects b_i(l) and can also include effects of other baseline covariates, for example confounding factors or ancestry-related principal components. To simplify the presentation, we assume a linear trajectory, but the longitudinal sub-models can be adapted for nonlinear trajectories using, for example, higher order polynomials or splines²³.

Multivariate mixed model

Where:

• β_0(l), β_1(l), β_g(l) and β_h(l) denote, respectively, the fixed intercept, and baseline effects of time, SNP and other covariate(s) on each longitudinal trait l.

• b_i = ((b_i0(l), b_i1(l))^′, …, (b_i0(L), b_i1(L))^′) ^′ is the 2xL vector of random effects with b_i∼N_2L(0, D).

• denotes the overall variance-covariance matrix for random effects, accounting for serial dependencies within longitudinal traits (D_l,l, unstructured) and cross-dependencies between longitudinal traits (D_l,m for two traits l≠m).

• is the vector of error terms where is assumed independent and identically distributed.

• is the residual variance of the QT l. We further assume that ε_ij(l) and b_i(l) are independent²².

(ii)Time-to-event sub-model

The time-to-event sub-model is specified by a Proportional Hazards (PH) mixed effects model (also known as frailty PH model) where the hazard function of each time-to-event trait k depends on the SNP effect adjusted for a function of the subject’s longitudinal trajectories W_i(k)(t). We introduce a subject-specific random effect (frailty term, u_i) to capture unexplained dependencies (e.g. due to unmeasured shared factors) among the K time-to-event traits. Given the common frailty term u_i, the K time-to-event traits are assumed independent²⁴.

Proportional Hazards mixed effects sub-model

Where:

• λ_0(k)(t) is a (parametric or non-parametric) baseline hazard function for the time-to-event trait k.

• γ_g(k) denotes the SNP effect on the time-to-event k adjusted for the indirect SNP effects via longitudinal QT(s), taken into account by W_i(k)(t).

• γ_v(k) is the effect of the baseline covariate vector(s). Those covariates can be shared or trait-specific covariates. They can also be time-dependent covariates for all or a subset of the K time-to-event traits (in this case, V_i(k) = V_i(k)(t)).

• u_i and b_i, random effects from the time-to-event and longitudinal sub-models respectively, are assumed independent. We assume u_i∼Γ(a, b), with a, b>0.

(iii) Time-dependent association structures

The time-dependent association structures, W_i(k)(t), that connect the longitudinal and time-to-event sub-models (Figure 1), account for the specific temporal relationships between each set of L_k (1 ≤ L_k ≤ L) associated QT(s) and each time-to-event trait k.

Where:

• 1 ≤ L_k ≤ L is the set of longitudinal risk factors of the time-to-event trait k.

• denotes the functional form of the exposure effect on time-to-event trait k corresponding to the l-th longitudinal risk. In the case of a contemporaneous parametrization, the hazard of an event k at a time t depends on the longitudinal trait value at the same time t (i.e. . Other parametrizations have been described in the literature^6,7,25.

• α_l(k) denotes the association parameter of the function of the smooth longitudinal trajectory with the time-to-event k.

• , with 0 ≤ s ≤ t, 1 ≤ l ≤ L} denotes the history of the underlying longitudinal trajectory for trait l up to time t.

Interpretation

When the joint model is correctly specified, the parameters accurately represent multiple relationships between each SNP and the correlated traits. Note that β_g(l) is the SNP effect on the longitudinal risk factor l, γ_g(k) is the direct SNP effect on the time-to-event trait k via other mechanisms than the indirect SNP association induced by the observed longitudinal risk factors. In this context, θ_k(l) = µ_g(k,l) + γ_g(k) is interpreted as the overall SNP effect on the time-to-event trait k, where µ_g(k,l) = α_l(k)β_g(l) represents the indirect¹ SNP effect on the time-to-event trait k via the longitudinal QT l.

Under the proposed joint model formulation, when a longitudinal risk factor l is associated with a time-to-event trait k (α_l(k) ≠ 0), a SNP that has an effect on the longitudinal QT risk factor l (β_g(l) ≠ 0), but no effect on the time-to-event trait k (γ_g(k) = 0), is interpreted as having an indirect SNP association and its overall effect is equal to its indirect effect (θ_k(l) = µ_g(k,l), with µ_g(k,l) = α_l(k)β_g(l)). A SNP with an effect on the time-to-event trait k (γ_g(k) ≠ 0), but no effect on the longitudinal QT l (β_g(l) = 0), is interpreted as having a direct SNP association and its overall effect is equal to its direct effect (θ_k(l) = γ_g(k)). A SNP with an effect on the longitudinal risk factor l (β_g(l) ≠ 0) and an effect on the time-to-event trait k (γ_g(k) ≠ 0) is interpreted as having both direct and indirect SNP associations via distinct genetic pathways. In this case, its overall effect is an aggregation of both direct and indirect effects (θ_k(l) = µ_g(k,l) + γ_g(k), with µ_g(k,l) = α_l(k)β_g(l)). It is obvious that when an associated longitudinal risk factor is omitted from a time-to-event model, as occurs in separate analyses of longitudinal and time-to-event traits, the marginal SNP effect on the time-to-event trait will reflect partially the overall SNP effect θ_k(l). This is particularly the case in GWAS when intermediate risk factors are ignored in the analysis models.

Effect estimation and test statistic construction

To mitigate the computational complexity involved in the maximization of the joint likelihood and allow more flexibility in the model formulated, we estimate the parameters using a two-stage approach (see Appendix). We work within the framework defined by Tsiatis and Davidian^26,27 which guarantees that our estimators are consistent and asymptotically normal. Specifically, in Stage 1 we fit the multivariate mixed model (Equation 1) using the mvlme() function from the R package JoineRML²⁸ to estimate the parameters of the longitudinal trajectories of the risk factors (Equation 2), and obtain fitted values of the smoothed trajectories. In Stage 2, we fit a Cox PH frailty time-to-event model (Equation 3) adjusting for functions of the smoothed trajectories as time-dependent covariates using the coxph() function from the R survival^29,30 package. We assume a Gamma distribution for the frailty term (u_i) and a separate baseline hazard function for each time-to-event trait using the strata argument in coxph(). To account for propagation of errors, due to the uncertainty in Stage 1 estimates ignored during Stage 2 parameter estimation¹⁶ and empirically estimate the joint covariance matrix of SNP and trajectories effects we apply a nonparametric bootstrap³¹. For each bootstrap sample b (1 ≤ b ≤ B, B is the total number of bootstraps), 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 of the estimated parameters using the B bootstrap parameter estimates. Wald statistics for each β_g(l) are computed as using the empirical bootstrap standard errors , and similarly for each γ_g(k) as .

Procedure to classify direct/indirect SNP associations

In Table 1, we present a practical hypothesis testing procedure to classify SNPs as having direct and/or indirect association with a time-to-event trait k and each associated longitudinal risk factor l. This procedure requires two significance thresholds, and for β_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 one . The latter is applicable for instance, when we want to systematically classify direct/indirect association among a set of M SNPs, and the former when assessing which SNPs, among the ones reported to be associated with the longitudinal risk factor, have a direct effect on a time-to-event trait.

View this table:

• View inline
• View popup
• Download powerpoint

Table 1. Hypothesis testing procedure to classify a SNP as having a direct and/or an indirect association with a time-to-event trait k and an associated longitudinal risk factor l.

Simulation study

Design of the DCCT-data-based simulation study

To assess the performances of the methods, we simulate the data under a genetic causal scenario that imitates the complex genetic architecture of T1DC (Figure 2). The genetic association model involves N=667 DCCT subjects with K=2 non-independent simulated time-to-T1DC traits (DR, DN) that depend on M=5 causal SNPs with direct effects on each time-to-T1DC and/or indirect effects via L=3 longitudinal QTs: two as measured in DCCT (HbA1c, SBP) and one other simulated (U) to induce some unexplained shared dependency between the T1DC traits. We also introduce an effect of sex on SBP, and effects of T1D duration (at baseline) on both T1DC traits, as observed in the original DCCT data.

• Download figure
• Open in new tab

Figure 2. Realistic DCCT-data-based causal genetic scenario

The effect of gender on SBP (noted as H_i(2)), and effects of T1D duration at baseline (noted as V_i(k)) on both time-to-T1DC (k=1, 2) traits are not represented in this figure, but are included in the data generating model, see Supplementary Information SI-1 for details. We simulated R=1000 replicates of N=667 DCCT individuals with M=5 causal variants and K=2 time-to-event traits simulated under this causal genetic scenario and R=1000 replicates of M=5 SNPs (with same MAF as the causal ones) under a global null genetic scenario where none of the SNPs is associated with any traits.

Algorithm for realistic data generation under a complex genetic architecture

To generate such a data structure that combines observed and simulated traits, we formulate a data generating model including: (i) L=3 linear mixed models linking each of the SNPs with an indirect effect to each longitudinal risk factor, and (ii) K=2 non-independent parametric time-to-event models depending on SNPs with direct effects and on functions of the longitudinal QT trajectories. For each DCCT individual i with observed HbA1c (Y_i(1)), SBP (Y_i(2)), visit times (t_i) vectors at n_i time points, and covariates (H_i(l), V_i(k)), we simulate longitudinal trait vector U_i, time-to-event traits and genetic data at M causal SNPs as illustrated in Figure 3 and detailed in Supplementary Information SI-1. All SNPs are simulated under Hardy-Weinberg and linkage equilibrium assumptions. Particularly, SNPs with indirect effects are simulated from the observed (SNP1, SNP5) or simulated (SNP3) longitudinal QTs, while SNPs with direct effects (SNP2, SNP4) on time-to-T1DC traits are simulated independently of the longitudinal QTs and are included in the specified hazard function used to simulate each time-to-event trait k (Figure 3).

• Download figure
• Open in new tab

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

Ω₍₁₎, Ω₍₂₎ and Ω_(U) are specified parameter values for the longitudinal trait models, and Γ₍₁₎, Γ₍₂₎ for the time-to-event trait models, p is the specified vector of minor allele frequencies vector of the M SNPs. Parameters for the causal genetic scenario in Figure 2 are summarized in Table 2 and in Table S1. The DCCT-data-based simulation algorithm uses longitudinal trait/baseline covariates from DCCT and simulated genetic data and time-to-event traits. SNPs with indirect effects are simulated from a multinomial distribution with conditional genotype probabilities (π_0i, π_1i, π_2i) using longitudinal trait values for individual i. Here, X_i and Z_i denote the design matrices of the fixed and random effects of the mixed model assumed for each longitudinal trait and D is the specified covariance matrix for the random effects. SNPs with direct effects are simulated from the population probabilities, that depend on the MAF. Each time-to-event trait k is simulated by calculating the inverse of the cumulative specified hazard function using the Brent univariate root-finding method^34,35. For DR, and for . We further added the effect of the longitudinal trajectory to each η_i(k) in the hazard function of each trait k to induce some unexplained dependencies between the simulated time-to-event traits. Details of the DCCT-based simulation procedure are presented in Supplementary Information SI-1.

For each SNPm, we specify MAF (p_m) and the SNP effects (β_g(l), γ_g(k)), as well as the other parameter values according to the DCCT Genetics Study and the T1DC literature (Table 2, Tables S1). We assume contemporaneous association structures for HbA1c and SBP effects on simulated time-to-T1DC traits. We simulate each time-to-event trait k under a Weibull model, with shape and scale parameters specified to generate ∼54% DR events and ∼25% DN events in each of the R=1000 replicated datasets generated under the causal scenario from Figure 2. Under the global null genetic scenario, where none of the SNPs is associated with any traits, we simulate M SNPs independently of the traits with the same MAF as for the causal SNPs.

View this table:

• View inline
• View popup
• Download powerpoint

Table 2. Simulated scenarios of direct and/or indirect SNP associations, assessed using R=1000 replicates of N=667 DCCT subjects simulated under the causal genetic scenario.

This table presents parameters specified for SNPs and trajectory effects on time-to-event traits for the DCCT-based simulation study and illustrates how each type of SNP would be detected in GWAS discovery analysis based on separate analysis of each trait.

Scenarios for DCCT-based complex genetic architecture

Overall, the simulated complex genetic architecture covers multiple types of SNP-trait associations (Figure 2 and Table 2): direct association with each T1DC trait (SNP2, SNP4), indirect associations with both T1DC traits via measured (SNP1) and unmeasured (SNP3) longitudinal QTs; and direct and indirect associations via a measured longitudinal QT (SNP5); all longitudinal risk factors are subject to measurement error. Except for SNP3, all other SNP scenarios represent GWAS discovery SNPs detectable by separate GWAS of a longitudinal risk factor (SNP1, SNP5) or a time-to-event trait (SNP2, SNP4, SNP5).

SNP1, SNP3 and SNP5 have indirect effects on T1DC traits, such that their associations with the T1DC traits are detectable with the marginal Cox PH time-to-event models (Tables 2 and S2). SNP1 corresponds roughly to rs10810632 and rs1358030 reported in the motivating DCCT GWAS of HbA1c¹⁹, while SNP5 represents a top hit detected by each GWAS (longitudinal and time-to-event traits), analyzed separately. Because of their association with longitudinal and time-to-event traits, SNP1 and SNP5 associations with time-to-event outcomes are conventionally investigated using the Cox PH model adjusting for the observed longitudinal trait as a time-dependent risk factor to determine if SNP association with the time-to-event trait is fully (or partially) explained via the associated longitudinal QT(s).

Results

Simulation results

Estimation accuracy and single-parameter tests validity

When the longitudinal QT(s) are observed and measured with random errors, the simulation results confirm that the proposed joint model reduces the bias of the parameter estimates for all types of SNP association (Figure S1 and Tables S3-S6), even when the analysis model is not fully specified. This applies when: (i) two correlated time-to-event traits are analyzed jointly rather than separately in the joint model (JM-cmp/JM-mis compared to JM-sep(k)); (ii) a SNP has a direct effect on a time-to-event trait with a low event rate (SNP4/SNP5, JM-cmp/JM-mis compared to CM-obs); (iii) a SNP has an indirect effect via a longitudinal QT measured with errors (for example SNP5, JM-cmp/JM-mis compared to CM-obs). Particularly, for SNP5, CM-obs overestimates the direct SNP5 effect on time-to-DN due to its inability to fully account for the indirect SNP5 effect induced via SBP. Direct SNP2 and SNP4 effects on time-to-T1DC traits show overall the larger relative bias reduction using joint models versus CM-obs (Ranges with JM-mis 0.702-0.983, with JM-cmp 0.934-0.966, see Tables S4 and S6).

When a SNP has an indirect effect on both T1DC traits fully explained via an unmeasured longitudinal QT (SNP3), all methods that ignore the longitudinal QT (JM-mis, JM-sep(k) or CM-obs compared to JM-cmp) produce biased direct SNP3 estimates for both T1DC towards the overall SNP3 effect (θ_k(U) = α_U(k)β_g(U) =0.36, Figure S1-c and S1-d, Table S5). In this scenario, the use of the frailty term in JM-mis does not reduce this bias (JM-mis vs JM-sep).

Classification of direct/indirect SNP associations with time-to-event traits

Under the global null genetic scenario of no genetic association with any of the traits, single-parameter SNP test P-values and from the joint model do not show departure from the expected large sample distributions (χ² with 1 degree of freedom (df), see Figure S2) and the type I error of each test under the null is reasonably well controlled (Table S7). The hypothesis testing procedure applied to single-parameter SNP tests, classifies each SNP as direct or indirect association with each T1DC trait with rates close to the nominal level (Table 3 and Figure S3). However, the classification rates for test of direct and indirect SNP association test appears less than nominal (see Table 3 for SNP5 with SBP/DN traits and Figure S3 for the other SNPs with other QT/time-to-event trait pairs).

View this table:

• View inline
• View popup
• Download powerpoint

Table 3. Classification rates of the procedure for each SNP under the global null genetic scenario, assessed using R=1000 replicates of N=667 DCCT subjects.

Under the causal scenario, the hypothesis testing procedure has high correct classification rates for indirect (SNP1) or direct associations (SNP2, SNP4), as shown in Figures 4-a, 4-b and S4. Conclusions are similar for JM-mis, with attenuated correct classification rates compared to JM-cmp for direct SNP associations but still higher than classification rates based on CM-obs (particularly for SNP4). Compared to results based on CM-obs, the largest relative classification improvement is found for direct SNP4 association with DN (ranges: 1.18-251.9 with JM-cmp and 1.18-110.98 with JM-mis, Figure 4).

• Download figure
• Open in new tab

Figure 4.

Performance of the testing procedure to recover underlying generating causal relationships from Figure 2, assessed using R=1000 replicates of N=667 DCCT subjects.

The Y axis represents the proportion of replicates that accurately classify each causal SNP as: (a.) indirect; (b.) direct and (c.) direct and indirect causal associations for SNP1, SNP2, SNP4 and SNP5. SNP3 is tested as a direct association because the longitudinal risk factor U is assumed unmeasured. Results obtained for direct association of SNP3 and SNP4 with DN/SBP are similar to those presented for DN/HbA1c (Figure S4) and are not shown here. Plot (d.) shows that SNP5 is misclassified as an indirect association in a larger number of replicates compared to (c.) by the joint models compared to CM-obs. This misclassification rate tends to increase with decreasing P^∗. Classification rates of all these SNPs with other time-to-event/longitudinal trait pairs are presented in Figure S4.

On the other hand, SNP3 with an indirect effect via the unobserved longitudinal trait U, is most frequently detected incorrectly as a direct association with each T1DC trait by both JM-mis and CM-obs (Figure 4-b). SNP5, which has both direct and indirect effects on DN induced via SBP, is most frequently misclassified by JM-cmp and JM-mis as an indirect association at significance levels varying from P^∗ = 0.01 to 5×10⁻⁸ (Figure 4-c and 4-d). SNP5 misclassification rate even increases with P^∗ decreasing to 5×10⁻⁸, due to the direct SNP5 . In comparison, the procedure based on CM-obs results most frequently correctly classify SNP5 association in more than 75% of the replicates at the same significance levels (Figures 4-c and 4-d); this is explained by the overestimated direct SNP5 effect, arising from measurement error in SBP as noted before.

In summary, our simulations show that by accounting for measurement errors in the longitudinal QT risk factors and for dependencies between/within the traits, the proposed joint model improves estimation accuracy and correct classification of associations of SNPs directly associated with each time-to-event trait or indirectly associated via a measured longitudinal QT in comparison to classification using the CM-obs approach. However, when a SNP has both direct and indirect effects on a time-to-event trait, the proposed testing procedure can be conservative since it requires the joint significance of the two SNP effects, β_g(l) and γ_g(k), where the power of each test depends on effect size and the trait distribution. As a result, a SNP with a direct and an indirect association can be misclassified as either a direct or an indirect association. Finally, when a SNP has an indirect effect on both T1DC traits via an unmeasured QT, the testing procedure based on JM-mis, that captures some of the unexplained dependency between time-to-event traits through the frailty term, does not prevent misclassification as a direct association. This observation also shows the importance of the correct specification of the joint model - including all the intermediate QT(s) - to avoid misclassification of direct and/or indirect SNP associations.

Application in the DCCT Genetics Study data

We demonstrate the feasibility of the proposed approach by an application in the DCCT Genetics Study data. We use time to mild DR and time to persistent microalbuminuria, for DR and DN outcomes respectively, as previously defined in the motivating GWAS of HbA1c¹⁹ (see Supplementary information SI-2 for details). Genome-wide genotyping in DCCT subjects was performed using HumanCoreExome Bead Array (Illumina, San Diego, CA, USA) and standard quality controls procedures were applied to individuals and genetic markers^19,35. Ungenotyped autosomal SNPs were imputed using 1000 Genomes data phase 3³⁶ (v5) and minimac3³⁷ (v.1.0.13), as previously described³⁵.

Out of the 667 DCCT individuals, we analyze N=516 subjects with genetic data, without mild to moderate non-proliferative retinopathy or without DN event at DCCT baseline. By the time of the DCCT close-out visit, 297 (57.6%) experienced a DR event, 61 (11.8%) a DN event, including 47 subjects (9.1%) that experienced both events. After SNP filtering and pruning on linkage disequilibrium (see Supplementary Information SI-2 for details), we analyze 307 candidate SNPs reported as associated with HbA1c, SBP, and multiple definitions of DR and/or DN^{19,34,38–42} (see Table S8 for the full list of SNPs).

Compared to the DCCT-based simulated datasets, the proportion of DN events observed in the application is lower. We expect that the test power for direct SNP association with DN under the joint model to be reduced; this implies reduce chances to correctly classify a SNP with either a direct or a direct and indirect association with DN. We also observe larger effects of HbA1c on T1DC traits (α_1(k), k=1, 2) and a smaller effect of SBP on DN (α₂₍₂₎), as shown in Tables 2 vs S9. Larger α_l(k) values increase the contribution of the indirect effect (µ_g(l,k)) to the overall SNP effect (θ_k(l)), which can result in a direct SNP effect (γ_g(k)) more severely biased if the longitudinal QT is not accounted properly in the joint model, while lower α_l(k) values reduce µ_g(l,k) and thus γ_g(k) is less severely biased. Overall, in the DCCT application, we expect a SNP with an indirect effect via HbA1c on any T1DC trait to be more subject to be misclassified as having both a direct and indirect association, compared to the simulation results, while a SNP with an indirect effect via SBP on DN is less subject to this misclassification.

Given prior evidence for cumulative effects of HbA1c on T1DC traits^20,21, we compare joint model results obtained with contemporaneous, updated cumulative mean, and time-weighted cumulative HbA1c effects on T1DC (Supplementary Information SI-2 for details). We present results under the latter association structure since it has a stronger prior association with T1DC and in the DCCT individuals analyzed here (Table S9). We obtain similar conclusions with alternative association structures (Figure S5-a). As shown in Figure 5-a, rs10810632 and rs1358030 are classified as indirect associations with both T1DC traits via their association with HbA1c shared risk factor ( and , P^∗=1.7×10⁻⁴ after Bonferroni correction for the 289.02 effective SNPs tested⁴³). These two SNPs have calculated indirect effects (µ_g(1,k)) and 95% bootstrap confidence intervals that do not include 0 (Figure 5-b). As noted before, rs10810632 and rs1358030 were discovered previously in a GWAS of HbA1c in DCCT¹⁹. Classification of these two SNPs as direct and/or indirect associations may be affected by the Winner’s curse bias, due to the modest sample size in this analysis, and the same individuals contributing to discovery and classification^44,45. For the other candidate SNPs which were discovered in larger and/or independent studies^{19,34,38–42}, we find little evidence of association with the discovery traits in marginal analyses (Table S8).

• Download figure
• Open in new tab

Figure 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) (X axis) and γ_g(k) (Y axis) for DR/HbA1c, DN/HbA1c and DN/ SBP trait pairs. Significance levels P^∗ =1.7×10⁻⁴ and P^∗ = 0.05 are indicated by red and grey horizontal and vertical dashed lines. (b.) Association results for the two SNPs that exibit indirect association with DR/DN outcomes at P^∗ =1.7×10⁻⁴. Left panel present results from separate analysis of each trait; and right panel presents results from the joint model with bootstrap 95% confidence intervals for the direct and indirect SNP effects. Results are presented using the time-weighted cumulative HbA1c effects on T1DC traits.

Conclusions regarding classification of direct and/or indirect associations based on the procedure using from CM-obs are similar to those based on the joint model results for all SNPs (Figure S5-b), although we noticed some differences in ; that could be explained by the attenuated association of HbA1c and SBP with T1DC traits due to measurement errors ignored from CM-obs (Tables S10 vs S9).

Discussion

We present new methods for statistical genetic analysis under a joint model specification of multiple time-to-event traits and multiple longitudinal risk factors designed to characterize the complex genetic architecture of related traits in longitudinal studies of disease progression. The proposed model accounts for dependencies within and between traits and can account for trait-specific and shared covariates, measurement errors in the longitudinal traits, as well as effects of unobserved baseline confounding factors between the time-to-event traits through the subject-specific frailty term.

Evaluation by realistic data-informed simulation study of complex T1DC genetic architecture shows that the proposed joint model improves estimation accuracy and efficiency compared to the Cox PH frailty model adjusted for observed longitudinal trait values as time-dependent covariates. This improvement holds particularly when the longitudinal risk factors are measured with errors or the time-to-event outcome has a low event rate. Indeed, measurement error attenuates the association of the longitudinal risk factor with the time-to-event trait, and as a result, if a SNP has an indirect effect induced by the QT, the estimated direct SNP effect on a time-to-event is biased away from zero. The magnitude of this bias tends to be worse with larger measurement error^14,15. Identification of direct and/or indirect SNP associations is also facilitated under joint modelling in the presence of risk factor measurement error. When, however, a SNP has an indirect effect via an unmeasured longitudinal risk factor, both the proposed and Cox PH model approaches produce biased estimates of the direct SNP effects. The extent of this bias depends on the magnitude of the effect sizes β_g(l) and α_l(k). Although the frailty term in the joint model time-to-event sub-model captures some of the unexplained shared dependency between the time-to-event traits, we did not observe bias reduction of the direct SNP estimate in the case of an unmeasured risk factor, which can lead to misclassification of the SNP as a direct association. This may be a limitation of using a time-invariant frailty term that does not fully account for the indirect SNP effect via the unmeasured longitudinal risk factor. Because joint significance of both SNP effects is required to classify a SNP as having both direct and indirect associations, a SNP with both direct and an indirect association can be misclassified as having either a direct or an indirect association if one of the SNP effects is too small to be detected at the specified level. The proposed procedure can be adapted to other research questions, for example to identify SNPs that have direct association among the SNPs found associated with a longitudinal risk factor. In this case, it may be desirable to apply a less stringent significance threshold for the test of the direct SNP effect while controlling the overall type I error. Alternatively, intersection-union tests, as based on the likelihood ratio test⁴⁶, could be used to assess joint significance of β_g(l) and γ_g(k) but would require computation of the full joint likelihood.

Application of the proposed joint model approach in longitudinal studies of disease progression, such as in the DCCT Genetics Study, improves classification of direct and/or indirect SNP association which can help to elucidate the genetic architecture of complex traits. In the context of mediation analysis, Liu et al⁴⁷ discussed various formulations and interpretations of joint models, with shared-random-effects accounting for potential unmeasured baseline confounding factors between one longitudinal and one time-to-event traits. Using applications in datasets from two clinical trials, they illustrate interpretation of sensitivity analysis to unmeasured baseline confounders. Adaptation of the joint model we propose for multiple longitudinal and multiple time-to-event traits for mediation analysis requires extension of the mediation assumptions^48,49 to the case of multiple mediators and multiple time-to-event traits. Specific evaluations of the proposed model under these assumptions are also warranted.

Although our primary aim in this report is to develop statistical methods to accurately distinguish among direct and/or indirect SNP associations with each time-to-event trait, the multi-trait aspect of the joint model lends itself to development of multi-trait SNP association testing for SNP discovery. In Supplementary Information SI-3, we present a joint-parameter test based on a generalized Wald statistic. In application to the simulated DCCT-based complex genetic architecture, we observe good type I error control under the global genetic null scenario, improved power for SNP discovery when a SNP has multiple trait effects, and power maintenance in other SNP association scenarios.

We acknowledge several features of the proposed joint model approach that warrant examination in further work. Firstly, to reduce computational complexity and improve model flexibility, we used two-stage parameter estimation. In some circumstances, this approach can produce biased estimates and/or underestimated standard errors¹⁶. Biased estimates can result from non-random censoring of the longitudinal trait values due to the occurrence of an event or from informative dropout^50,51. Our simulation results show minimal biases in the absence of informative censoring, even when the joint model is mis-specified. In the DCCT application, characterized by administrative censoring and a high completion rate, these biases are of minimal concern because longitudinal trait values continued to be recorded regardless of the occurrence of any T1DC events; we estimated the trajectories using all the available measurements. Furthermore, we obtain robust bootstrap estimates of the covariance matrix, and simulation results under the null did not show deviation from expected distributions. Use of the bootstrap also facilitates joint testing of parameters from both stages. In the presence of informative censoring, we recommend sensitivity analysis using existing implementations joint likelihood estimation. To our knowledge those implementations only exist for simpler joint model formulation with either one longitudinal and one time-to-event trait⁵² or multiple longitudinal traits and one time-to-event outcome^28,53. Secondly, because the joint model integrates longitudinal and time-to-event sub-models, model misspecification can occur in multiple ways and lead to invalid inference⁵⁴. Therefore, we recommend a careful assessment of the model assumptions using joint model diagnostic tools²³. Thirdly, patient visits were scheduled with high frequency in DCCT, so we ignored the modest degree of interval censoring in the current implementation of the joint model; when there are longer gaps between visits, extended methods are needed to account for interval censoring with additional simulation studies to assess impact on joint model estimates. Lastly, joint models are very computational demanding, particularly for genetic association studies that test millions of variants. In DCCT application, it took ∼60 seconds to fit the joint model for each SNP and ∼1090.32 more seconds (∼ 18 minutes) to estimate the covariance matrix with 500 bootstraps run in parallel on 4 nodes (each node with 40 CPU and 202 GB RAM). While analysis at the genome level involving, up to 85 million 1000G-imputed SNPs in DCCT, seems computationally unrealistic, a screening approach without bootstrap to select the SNPs based on their association P-values, followed by the bootstrap refinement would reduce the computational burden. Recent computationally efficient algorithms have been developed to improve the efficiency of the linear mixed model⁵⁵ and Cox PH model^56,57 for genetic association studies, but to date, they remain to be implemented for multivariate outcomes.

With the increasing development of national biobanks^58–60 consisting of a large collection of phenotypes, environmental factors, biomarkers and other risk factors collected over time and linked with genetic data, we anticipate that joint model methods can be applied to characterize the genetic architecture of complex traits using these data resources⁶¹. In addition, the results of such analysis could contribute towards the translation of human genetic findings to personalized medicine by providing more efficient SNP effect estimates for increased precision in polygenic risk score development⁶², and causal inference using mediation and mendelian randomization studies. Lastly, the joint model framework can also enable dynamic prediction beneficial for dynamic risk assessment^7,63 and optimization of intervention strategies^64,65.

Contributions

All named authors affirm that authorship is merited based on the International Committee of Medical Journal Editors (ICMJE) authorship criteria. MB, SBB and RVC designed the study. MB conducted the simulation study and data analysis and drafted the manuscript. SBB, RVC and OEG contributed to the data analysis and interpretation of the results and reviewed/edited the manuscript. ADP contributed to the acquisition of the data and reviewed/edited the manuscript. All authors approve the final version to be published. SBB and MB are the guarantors of this work.