A Bayesian Framework for Functional Mapping through Joint Modeling of Longitudinal and Time-to-Event Data

Abstract

The most powerful and comprehensive approach of study in modern biology is to understand the whole process of development and all events of importance to development which occur in the process. As a consequence, joint modeling of developmental processes and events has become one of the most demanding tasks in statistical research. Here, we propose a joint modeling framework for functional mapping of specific quantitative trait loci (QTLs) which controls developmental processes and the timing of development and their causal correlation over time. The joint model contains two submodels, one for a developmental process, known as a longitudinal trait, and the other for a developmental event, known as the time to event, which are connected through a QTL mapping framework. A nonparametric approach is used to model the mean and covariance function of the longitudinal trait while the traditional Cox proportional hazard (PH) model is used to model the event time. The joint model is applied to map QTLs that control whole-plant vegetative biomass growth and time to first flower in soybeans. Results show that this model should be broadly useful for detecting genes controlling physiological and pathological processes and other events of interest in biomedicine.

1. Introduction

To study biology, a classic approach is dimension reduction in which a biological phenomenon or process is dissected into several discrete features over time and space. Most efforts in the past decades have been made to understand biological details of individual features and then use knowledge from each feature to draw an inference about biology as a whole. There has been increasing recognition of the limitation of this approach because it fails to detect a rule that governs the transition from one feature to next, thus leading to a significant loss of information behind the development of a biological trait. More recently, tremendous developments in statistics and computer science have enabled scientists to model and compute the dynamic behavior of a biological phenomenon and construct a comprehensive view of how a cell, tissue, or organ grows and develops across the time-space scale.

A statistical dynamic model, called functional mapping, is one of the products of such developments [1, 2]. The merit of functional mapping lies in its biological relevance to study the tempo-spatial pattern of change for the trait and further predict the physiological or pathological status of trait phenotype. Functional mapping has proven to be powerful for elucidating the dynamic genetic architecture of complex phenotypic traits by identifying when specific genes (known as quantitative trait loci or QTLs) involved turn on and turn off and how long they are expressed in a time course. With the advent of new automatic techniques that collect dynamic data in a cost-effective way, functional mapping can be anticipated to play an increasingly important role in shedding light on the genetic control mechanisms of complex traits or diseases.

The statistical foundation of functional mapping is longitudinal data analysis or functional data analysis. There has been a considerable body of literature on statistical modeling of time-varying mean and covariance structure using various parametric, nonparametric, and semiparametric methods [3–7]. A joint mean-covariance model was proposed by Pourahmadi [8, 9], which shows some advantages over modeling the mean and covariance separately. Since the publication of the pioneering work by Laird and Ware [10], random effects model have been extensively used for longitudinal data analysis [11]. All these statistical approaches have been incorporated into functional mapping [12, 13], aiming to provide the most parsimonious estimates of QTL effects for a given data set. A Bayesian algorithm for functional mapping has been proposed recently by Liu and Wu [14].

The complexity of biology lies in the fact that no biological trait is isolated, rather every trait is affected by other traits through genes and environmental factors. For example, when a plant grows into a particular stage, reproductive behavior, such as flowering, starts to emerge as one of the important events in plant development. The time to first flower is highly associated with the amount of vegetative growth, depending on the environment where the plant is grown. Likewise, the time to recurrence of prostate cancer in humans is related with dynamic changes of prostate specific antigen level. How to jointly model longitudinal and time-to-event data within functional mapping has become an important issue for studying the common genetic basis of these processes and predicting events based on longitudinal traits.

Simultaneous modeling of longitudinal traits and time to events has been an active area in biostatistics during the past twenty years. A linear random effects model and EM estimation approach are proposed by Henderson et al. [15] for joint modeling. Guo and Carlin [16] made a comparative study between separate models and a joint model, showing that a joint model is more powerful when there is a strong correlation between the trait and the event. Wang and Taylor [17] developed a Bayesian method and MCMC algorithm for joint modeling of longitudinal and event time data and applied their algorithm on AIDS data. A review article by Tsiatis and Davidian [18] nicely summarizes the recent developments for such joint modeling.

By simply estimating the correlation between longitudinal traits and event time, Lin and Wu [19] developed a first model that connects these two aspects within functional mapping. However, they developed a likelihood-based framework where the covariance structure for the longitudinal trait was modelled by the known AR(1), and model parameters were estimated using maximum likelihood estimation. Taking advantage of event models, such as semiparametric Cox proportional hazard model, Weibull model, accelerated failure time (AFT) model, we here propose a sophisticated model for joint modeling of longitudinal trait and time to event to locate the QTLs which control the event via a dynamic trait. The detection of those QTLs that are common to these types of traits may help to prevent or accelerate the outcome by genetic approaches. Our model is constructed with a Bayesian paradigm and model parameters are estimated by the MCMC algorithm. Local polynomials are used to model the mean trajectory and generalized-linear-model- (GLM-) based approach is used to model the covariance matrix. The model is validated using a real example in which whole-plant biomass as a longitudinal trait measured at a series of discrete time points and the time to first flower as a time-to-event are jointly modeled through functional mapping. The statistical properties of the model applied to estimate QTL temporal effects in this example and its practical usefulness are investigated by simulation studies.

2. Joint Modeling Framework

2.1. Model for the Longitudinal Trait

Genetic mapping should be based on a segregating population, such as the backcross, F₂, or recombinant inbred lines (RILS), initiated with two inbred lines each carrying an alternative allele. An RIL population is generated by self-crossing the hybrids of the two inbred lines continuously for 7-8 successive generations, which leads to two homozygous genotypes for alternative alleles at each locus. Methods for other designs can be derived similarly. Suppose a backcross has 𝑛 progeny which is genotyped to construct a linkage map, aiming at locating putative QTLs that trigger significant effects on a longitudinal trait and its associated event. For each progeny, the trait is measured repeatedly at 𝑇 different time points and a time-to-event is also recorded. At a specific time point 𝑡, the phenotypic value of the trait for progeny 𝑖 affected by a putative QTL can be expressed by a linear model as follows:𝑦_𝑖(𝑡)=²_𝑗=1𝑥_𝑖𝑗𝑢_𝑗(𝑡)+𝑟_𝑖(𝑡)+𝑒_𝑖(𝑡), (1) where 𝑥_𝑖𝑗 is an indicator variable for a possible QTL genotype of progeny 𝑖 and defined as 1 if a particular QTL genotype 𝑗 is indicated and 0 otherwise (𝑗=1 for QTL genotype 𝑄𝑄 and 2 for genotype 𝑞𝑞), 𝑢_𝑗(𝑡) is the mean phenotypic value of QTL genotype 𝑗 for progeny 𝑖 at time 𝑡, 𝑟_𝑖(𝑡) is the subject specific random effect, and 𝑒_𝑖(𝑡) is the residual error assumed to follow a normal distribution with mean zero and covariance matrix 𝚺.

The central theme of functional mapping is to model the mean and covariance structures for the longitudinal trait efficiently. Here, we model the mean vector by polynomial function and the covariance matrix by an approach that guarantees the positive definiteness of the estimated covariance matrix. Without loss of generality, assume the response vector for progeny 𝑖, 𝐲_𝑖=(𝑦_𝑖(1),…,𝑦_𝑖(𝑇)), has mean 0 and covariance matrix 𝚺. The response at time 𝑡, 𝑦_𝑖(𝑡), can be predicted by its predecessors as the follows:𝑦_𝑖(𝑡)=_𝑡^𝑡−1_′=1𝜙_𝑡,𝑡′𝑦_𝑖𝑡^′+𝜖_𝑖(𝑡), (2) where 𝜙_𝑡,𝑡′ is the corresponding regression coefficient, 𝜖_𝑖(𝑡) is the prediction error for progeny 𝑖 with mean=0, and 𝜎²(𝑡) is its variance. Assuming that 𝜖_𝑖(𝑡)’s are uncorrelated (Pourahmadi [8]), we get cov(𝝐_𝑖)=𝐷, a diagonal matrix with 𝜎²(𝑡) being the 𝑡th diagonal element, where 𝝐_𝑖=(𝜖_𝑖(1),…,𝜖_𝑖(𝑇))^′ is the vector of prediction errors. Hence, the matrix representation of the above autoregression becomes𝝐_𝑖=𝐌𝐲_𝑖, (3) where 𝐌 is a lower triangular matrix with 1’s in diagonal elements and −𝜙_𝑡,𝑡′ in the (𝑡,𝑡′)th position. The above equation simply givescov𝝐_𝑖=𝐌cov𝐲_𝑖𝐌^𝑇=𝐌𝚺𝐌^𝑇=𝐃, (4) which is related to the modified Cholesky decomposition of Σ [20].

Equation (4) will be considered as the basis for modeling the covariance structure, since this guarantees the estimated covariance matrix to be positive definite. Following Pourahmadi [8, 9], we model the mean vector, unconstrained variance parameters log𝜎²(𝑡), and dependence parameter 𝜙_𝑡,𝑡′, using a polynomial function of a particular order, expressed as𝑢_𝑗(𝑡)=𝛽_𝑗0+𝛽_𝑗1𝑡+𝛽_𝑗2𝑡²+⋯+𝛽_𝑗𝑟𝑡^𝑟,𝑟_𝑖(𝑡)=𝜃_𝑖0+𝜃_𝑖1𝑡+𝜃_𝑖2𝑡²+⋯+𝜃_𝑖𝑚𝑡^𝑚,log𝜎²_𝑡=𝜂₀+𝜂₁𝑡+𝜂₂𝑡²+⋯+𝜂_𝑔𝑡^𝑔,𝜙_𝑡,𝑡′=𝛿₀+𝛿₁𝑡−𝑡^′+𝛿₂𝑡−𝑡′²+⋯+𝛿_ℎ𝑡−𝑡^′ℎ,𝑡^′=1,2,…,𝑡−1. (5)(6)(7)(8)

The optimal (𝑟,𝑚,𝑔,ℎ) is determined from the information criteria (AIC/BIC). We note that different genotypes are assumed to have the same covariance structure but different means. Note that the above method of modeling the covariance structure for a longitudinal response is more robust than the traditional first-order autoregressive (AR(1)) or compound symmetry (CS) structure since real data might not show a parametric dependence structure. We refer to the proposed approach as GLM-based approach to estimate the covariance matrix.

Denote 𝜷_𝑗=(𝛽_𝑗0,𝛽_𝑗1,…,𝛽_𝑗𝑟) for QTL genotype 𝑗 and 𝜽_𝑖=(𝜃_𝑖0,𝜃_𝑖1,…,𝜃_𝑖𝑚) for subject 𝑖. Then, the conditional mean function for progeny 𝑖 carrying QTL genotype 𝑗 (𝑗=1 or 2) for the given subject specific random effect (𝜽_𝑖) can be expressed as𝝁_𝑖𝑗=𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗+𝐗_𝑖^(𝑚)𝜽^𝑇_𝑖, (9) where𝐗_𝑖^(𝑟)=1𝑡_𝑖1𝑡²_𝑖1⋯𝑡^𝑟_𝑖11𝑡_𝑖2𝑡²_𝑖2⋯𝑡^𝑟_𝑖2⋮⋮⋮⋱⋮1𝑡_𝑖𝜏𝑡²_𝑖𝜏⋯𝑡^𝑟_𝑖𝜏⋮⋮⋮⋱⋮1𝑡_𝑖𝑇𝑡²_𝑖𝑇⋯𝑡^𝑟_𝑖𝑇. (10)

Assume the vectors of subject specific random effects 𝜽_𝑖 follow 𝑚-variate normal distribution with mean 0 and covariance matrix 𝜎²𝐼_𝑚 and they are independent of the residual errors. Note that under this assumption, 𝐲_𝑖|𝜽_𝑖 will follow MVN(𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗+𝐗_𝑖^(𝑚)𝜽^𝑇_𝑖,Σ), and the marginal distribution of 𝐲_𝑖 will be MVN(𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗,Σ+𝜎²𝐗_𝑖^(𝑚)𝐗_𝑖^(𝑚)𝑇).

2.2. Model for the Event Time

We use 𝑠_𝑖 to denote the event time of progeny 𝑖. Since in the current situation, the event time is recorded for all progeny; no progeny is censored. Assuming a Cox proportional hazard model for this event, we get for progeny 𝑖,𝜆_𝑖(𝑡)=𝜆₀(𝑡)exp𝛾𝜇_𝑖𝑗(𝑡), (11) where 𝜆₀(𝑡) denotes the baseline hazard at time 𝑡, 𝜇_𝑖𝑗(𝑡) is the mean longitudinal trait at time 𝑡 for given 𝜽_𝑖 when progeny 𝑖 is of QTL genotype 𝑗 and the regression coefficient 𝛾 represents the effect of the trait on the event time. The survival function for progeny 𝑖 can be expressed in terms of the hazard function as 𝑠_𝑖(𝑡)=exp[−∫^𝑡₀𝜆_𝑖(𝑢)𝑑𝑢].

The longitudinal model described above is linked to the hazard model by 𝛾. If 𝛾 = 0, then the event is independent of the trait, and hence we should better fit separate models for the trait and the event. However, when 𝛾 is different from zero, a joint model performs better than the separate models [17]. For simplicity, the baseline hazard is assumed to be a step function, 𝜆₀(𝑡)=𝜆_0𝑘 over a partition of the observed time scale [0, max(𝑠_𝑖)] into 𝐾 (possibly evenly spaced) intervals, (𝑡^𝜆₀=0,𝑡^𝜆₁,𝑡^𝜆₂,…,𝑡^𝜆_𝐾). The value of 𝐾 is usually not too large, possibly smaller than 10.

2.3. Likelihood for the Joint Model

Since the QTL genotype of a progeny is unknown, we use a mixture model to describe the likelihood of the progeny in terms of its possible underlying QTL genotypes [21]. The joint likelihood of unknown parameters Θ given the longitudinal trait 𝐲=(𝐲_𝑖)^𝑛_𝑖=1 and event time 𝐬=(𝐬_𝑖)^𝑛_𝑖=1 for all 𝑛 progeny can be expressed as𝐿(𝚯∣𝐲,𝐬)=^𝑛_𝑖=1²_𝑗=1𝜔_𝑗|𝑖𝜋𝐲_𝑖,𝑠_𝑖∣𝑄_𝑖=𝑗,𝜽_𝑖=^𝑛_𝑖=1²_𝑗=1𝜔_𝑗|𝑖𝑓𝐲_𝑖∣𝑄_𝑖=𝑗,𝜽_𝑖×𝜆_𝑖𝑠_𝑖exp−^𝑠𝑖₀𝜆_𝑖(𝑢)𝑑𝑢, (12) where 𝜋(⋅) denotes the joint density of the longitudinal trait and event time; 𝑓(.) denotes a multivariate normal with QTL genotype-specific mean 𝝁_𝑖𝑗 modeled as (9) and covariance matrix 𝚺; hazard function 𝜆_𝑖(𝑠_𝑖) is modeled as (11); and 𝜔_𝑗|𝑖 denotes the conditional probability of QTL genotype 𝑗 given that the marker information of projeny 𝑖 and 𝑄_𝑖 is the QTL genotype for the 𝑖-th subject.

The QTL genotype is inferred from marker genotypes of the linkage map. Let 𝑀_𝑖=(𝑀_𝑖1,…,𝑀_𝑖𝑚) be the 𝑚-marker genotypes for progeny 𝑖, 𝐷^∗ the position of the putative QTL measured by its distance from the very first marker of an ordered linkage group, and 𝐷_𝑘 the distances between marker 1 and 𝑘. Assume that the QTL is located between marker 𝑘 and 𝑘+1. Then, the conditional probability of QTL genotype 𝑗 given the genotype of these two markers that flank the QTL is expressed as𝜔_𝑗∣𝑖=Prob𝑄_𝑖=𝑗∣𝐷^∗,𝑀_𝑖𝑘,𝑀_𝑖(𝑘+1),𝐷_𝑘,𝐷_𝑘+1. (13) Note that, given the QTL locations 𝐷^∗, 𝐷_𝑘, and 𝐷_𝑘+1, one can compute 𝑑₁, the distance of the QTL from marker 𝑘 and 𝑑₂, the distance of the QTL from marker 𝑘+1 [22]. Using the Haldane map function, one can compute recombination fractions between marker 𝑘 and QTL (𝑟₁), between QTL and marker 𝑘+1 (𝑟₂), and between markers 𝑘 and 𝑘+1 (𝑟) as follows:𝑟₁=121−𝑒^−2𝑑1,𝑟₂=121−𝑒^−2𝑑2,1𝑟=21−𝑒^−2𝑑, (14) where 𝑑=𝐷_𝑘+1−𝐷_𝑘 is the distance between marker 𝑘 and 𝑘+1. Wu et al. [22] provide a procedure for deriving the conditional probabilities of QTL genotypes given marker interval genotypes for the backcross, F₂, and RIL populations, respectively.

Unknown parameters Θ in likelihood (12) contain QTL genotype-specific parameters 𝜷_𝑗, 𝜎², the parameters that model the variance structure and dependence structure 𝜼=(𝜂₀,𝜂₁,…,𝜂_𝑔) and 𝜹=(𝛿₀,𝛿₁,…,𝛿_ℎ), as shown in models (7) and (8), respectively, the effect of the longitudinal trait on the event 𝛾, QTL position 𝐷^∗ and the baseline hazards 𝜆_0𝑘.

2.4. Posterior Distribution and Sampling Procedure

We derive a Bayesian approach for estimating the unknown parameters. This will first need to specify the prior distributions for the parameters and, given the data and the priors, derive the posterior distribution over all the unknown parameters. For 𝛽_𝑗, we place a multivariate normal prior with zero mean and covariance matrix Σ_𝛽. An inverse gamma prior with parameters (𝛼₁,𝛼₂) is considered for 𝜎². We consider a uniform prior on 𝛾 and uniform (0,𝐷_𝑚) prior for the parameter 𝐷^∗. Independent Gamma (𝑎,𝑏) priors are taken for 𝜆_0𝑘, for 𝑘=1,…,𝐾. Priors for 𝜼 and 𝜹 are taken as MVN(0,Σ_𝜂) and MVN(0,Σ_𝛿), respectively.

With the above priors and likelihood function, we have the joint posterior distribution for the parameters. In this case, it is quite straightforward to get the full conditional posterior distributions. Assume that the priors are independent for different parameters. Thus, we get the posterior density of 𝛽, 𝜎², 𝜸, 𝐷^∗, 𝜼, 𝜹, 𝝀_𝟎 as𝜋𝛽,𝜎²,𝛾,𝐷^∗,𝜼,𝜹,𝝀_𝟎∣𝐲,𝐬∝𝜋𝐲,𝐬∣𝜷,𝜎²,𝛾,𝐷^∗,𝜼,𝜹,𝝀_𝟎×𝜋𝛽𝜋𝜎²𝜋𝛾𝜋𝐷^∗𝜋(𝜼)𝜋(𝜹)𝜋𝝀_𝟎. (15)

Assuming that priors for different genotypes are independent, we can express the above posterior distribution as𝜋𝛽,𝜎²,𝛾,𝐷,𝜼,𝜹,𝝀_𝟎∣𝐲,𝐬∝𝜋𝐲,𝐬∣𝜷,𝜎²,𝛾,𝐷^∗,𝜼,𝜹,𝝀_𝟎ײ_𝑗=1𝜋𝜷_𝑗𝜋𝜎²𝜋𝛾𝜋𝐷^∗𝜋(𝜼)𝜋(𝜹)𝜋𝝀_𝟎. (16)

The full conditional distributions for the model parameters, as derived in the Appendix, are used to estimate the parameters using the MCMC algorithm. Note that the full conditional distribution for 𝛽_𝑗 is expressed as a product of a normal distribution term which comes from the longitudinal trait and two other terms from the hazard model. To update 𝛽_𝑗, therefore, we use a Metropolis-Hastings (MH) algorithm with a normal proposal density since it is a part of its posterior distribution. We also note that in the full conditionals of 𝜂 and 𝛿, normal distribution coming from the longitudinal part of the data is the main determinant. Hence for 𝜂 and 𝛿, we consider normal proposals with the current value of the parameter as the mean and covariance matrix as 𝚺_𝜂 and 𝚺_𝛿, respectively. Selection of a good proposal density for 𝛾 is a bit tricky and we follow the recommendation given by Wang and Taylor [17]. By evaluating several choices for a good proposal, we consider a normal distribution with mean as the current state of the parameter and a suitable standard deviation in such a way that the proposed density gets well mixed with the target distribution (acceptance rate between 0.25 to 0.40). Because of conjugacy, we can directly simulate from the full conditional of 𝜆_0𝑘’s.

The parameter 𝐷^∗ which specifies the location of the QTL is updated following the idea of Satagopan et al. [23] by using the MH algorithm. A new value of 𝐷^∗, which we denote by 𝐷^∗new, is generated from Uniform (max(0,𝐷^∗−𝜓),min(𝐷^∗+𝜓,𝐷_𝑚)), where 𝜓 is the tuning parameter. Denote this proposed distribution by 𝑞(𝐷^∗,𝐷^∗new). The proposed value will be considered as the new value of the chain with probability𝛼𝐷^∗,𝐷^∗new𝜋𝐷=min1,^∗new∣𝑦,𝑠,𝛽,𝛾,𝜂,𝛿,𝜆₀𝑞𝐷^∗new,𝐷^∗𝜋𝐷^∗∣𝑦,𝑠,𝛽,𝛾,𝜂,𝛿,𝜆₀𝑞(𝐷^∗,𝐷^∗new). (17) We note that 𝜋(𝐷^∗new∣…)∝𝜋(𝑄∣𝐷^∗new,𝑀)𝜋(𝐷^∗new)=∏^𝑛_𝑖=1𝜋(𝑄_𝑖∣𝐷^∗new,𝑀_𝑖)𝜋(𝐷^∗new) and, similarly, 𝜋(𝐷^∗∣…)=∏^𝑛_𝑖=1𝜋(𝑄_𝑖∣𝐷^∗,𝑀_𝑖)𝜋(𝐷^∗).

Because of the independence among 𝑛 progeny, 𝑄 is updated by updating for each 𝑄_𝑖 separately. For each progeny 𝑖, the full conditional density is in the form of a multinomial with the following cell probabilities:Π𝑄_𝑖=𝑗∣𝑦_𝑖,𝑠_𝑖=𝜋𝑄_𝑖=𝑗𝜋𝑦_𝑖,𝑠_𝑖∣𝑄_𝑖=𝑗∑²_𝑗′=1𝜋𝑄_𝑖=𝑗^′𝜋𝑦_𝑖,𝑠_𝑖∣𝑄_𝑖=𝑗^′. (18) We can sample the QTL genotype directly from this full conditional density at each cycle. Details of the estimation procedure can be found in Satagopan et al. [23].

3. Application

3.1. Material and Analysis

The new model was applied to analyze a real data set for QTL mapping in soybeans. The mapping population contains 184 RILs derived from two cultivars, Nannong 1138-2 and Kefeng no. 1. A genetic linkage map of this population was first established by Zhang et al. [24] with 452 makers including RFLP, SSR, EST distributed among the 21 linkage group. This map was recently updated by adding some new SSR makers and dumping some unreliable markers. The new map contains 834 molecular makers covering a length of 2,308 cM in 24 linkage groups, with an average genetic distance of 2.85 cM between adjacent markers. Those markers with missing information were excluded from the analysis, leading to a total of 780 markers involved in our analysis.

The plants and their parents were grown in a sample lattice design with two replicates at Jiangpu Soybean Experiment Station, Nanjing Agricultural University, China. After 20 days of seedling emergence, plant biomass (in gms.) were measured once every 5–10 days until most plants stopped height growth. A total of 8 measurements were taken for the biomass and the time to get the first flower in that growing season was also recorded for each plant. Figure 1 shows the raw data, both the trait (biomass) and the event time for 184 plants.

Figure 1 

Whole-plant biomass growth trajectories for 184 soybean RILs. The time to first flower is indicated by a vertical line on each biomass growth curve. The black curve is the mean growth trajectory.

Prior distributions for the model parameters were taken as follows. For genotype specific fixed effect 𝛽_𝑗, a multivariate normal prior was used with zero mean and a diagonal covariance matrix with all diagonal elements 100. For 𝜎², we took 𝐼𝐺(3,1) prior which has mean=0.5 with small variance (0.25). A uniform prior 𝑈(−3.0,−0.1) was taken for 𝛾 following Wang and Taylor [17]. Observed time scale for the event time data was partitioned into 5 parts; that is, we took 𝐾=5 and considered independent gamma (0.04, 1.0) priors for 𝜆₀₁,…,𝜆₀₅, following Wang and Taylor [17]. Uniform prior was taken for 𝐷^∗. For 𝜂 and 𝛿, we considered multivariate normal priors with zero means and diagonal covariance matrices with diagonal elements 30 and 20, respectively. To investigate the effect of prior distributions on the estimation method, we did a sensitivity analysis. Considering different sets of priors, we fitted our model many times and it turned out the estimation is almost insensitive to the choice of priors. Hence the choice of our priors (even with huge variances for 𝛽, 𝜂, and 𝛿) does not affect much the estimation of the model parameters.

We fitted our joint model as described in Section 2, by MCMC sampling. We ran chains 120,000 times. To remove the effect of the starting values, we excluded first 20,000 burn-in iterations. With the remaining 100,000 iterations, we estimated the posterior distributions and the parameters were estimated by the posterior mean and also calculated the sample standard deviations for the posterior densities. For the MH algorithm, the acceptance rates were 0.31, 0.26, 0.25, and 0.30 for 𝜎², 𝛾, 𝜂, and 𝛿, respectively.

Since our model is complex, we perform several standard diagnostic tests to assess the convergence of the Markov chains. First, we use the method proposed by Brooks and Gelman [25]. Considering five different chains with different starting points and discarding the burn-in iterations, we computed multivariate potential scale reduction factors (MPSRF) to assess the convergence of the chains. Starting points for the model parameters were drawn from the respective priors. The computed values of this statistic get stabilized near 1 after 60,000 iterations for our model parameters, which indicates convergence of the chains.

Second, we perform Geweke test which compares the earlier part of the markov chain to the later part for assessing convergence. After deleting the burn-in iterations, from the remaining 100,000 iterations, we take out two subsequences; the first 50,000 and the last 50,000 iterations. Also consistent spectral density estimates at zero frequency are calculated to compute the 𝑧-scores. The calculated 𝑃values for our model parameters are above 0.18, which indicates a small absolute 𝑧-scores assessing the convergence of our chains. A detailed discussion of this method can be found in Geweke [26].

Finally, we perform the Heidelberger and Welch test as proposed by Heidelberger and Welch [27]. This test has two parts: a stationary test and a half-width test. Our chains after deletion of the burn-in iterations, pass the stationary test. To assess whether the number of iterations is adequate to estimate the parameters accurately, we calculate relative half-width (RHW). We consider the default 𝛼 value (0.05) and predetermined tolerance value is taken as 0.1. For all our model parameters, the calculated RHW is in between 0.045 and 0.081. This indicates that we have enough iterations to estimate the model parameters with 95% confidence under tolerance of 0.1.

3.2. Results

Figure 1 is the growth trajectories of whole-plant biomass over time for all RILs, in which the times to first flower for each RIL are also indicated. There are great variability found for these two traits among RILs. In the exploratory data analysis, we computed BIC values for different orders for the triplet (𝑟,𝑚,𝑔,ℎ), showing that the BIC value is smallest for the order (2,2,2,2), that is, second order polynomials can best fit the the mean, variance, and dependence structures (Table 1). By scanning for the existence of QTLs over the genetic linkage map, we obtained the posterior distribution of the model parameters and estimate marginal posterior distributions of the QTL locations (𝐷^∗) for all 24 linkage groups (Figure 2).

Table 1 

BIC values for selecting the optimum (𝑟,𝑚,𝑔,ℎ) to model the mean-covariance structures for whole-plant biomass growth trajectories and the time to first flower in soybeans.

Figure 2 

Marginal posterior plot for QTL locations over 24 linkage groups. Marker locations are indicated by ticks on the 𝑥-axis.

We observed posterior peaks in linkage groups 1, 4, 15, 19, 20, 21, and 23. To draw inference about the existence of a putative QTL in each of these groups, we computed the Bayes Factor (BF), defined as𝑃()BF=𝐘,𝐒∣𝜅=0𝑃(),𝐘,𝐒∣𝜅=1 (19) where 𝜅 denotes the number of QTLs in that particular group. Following Jeffrey’s scale, the BFS with value smaller than 1 gives strong evidence against the null hypothesis and higher than 10 gives enough evidence for the null hypothesis. Note that in this case, we are testing the existence of no QTL (null) versus the existence of one QTL (alternative) for each linkage group. Here, no QTL in a group means 𝛽₁=𝛽₂, for that group. So, in order to compute BF we run our MCMC twice, first under the null and then under the alternative. The calculated BF for the above 7 linkage groups were 0.2781, 11.274, 13.493, 10.610, 0.5913, 11.475, and 0.7953, respectively, implying the existence of QTL in groups 1, 20, and 23. Table 2 provides genotype-specific mean parameters for the QTLs located on linkage groups 1, 20, and 23, along with their 95 percent credible intervals (C.I.).

Table 2 

Estimates of the parameters that describe genotype-specific biomass growth trajectories and QTL locations on linkage groups 1, 20, and 23, with 95% credible intervals.

Since the nature of our model is complex and our estimation is based on MCMC, we perform posterior predictive check for the aforementioned 7 linkage groups. We simulate observations to get the posterior predictive distribution. Let Θ be the set of all model parameters and 𝐃^rep=(𝐲^rep,𝐬^rep) be the replicated data. Then given the data 𝐃=(𝐲,𝐬), the posterior predictive distribution of 𝐃^rep is given by 𝑝(𝐃^rep∣𝐃)=∫𝑝(𝐃^rep∣Θ)𝑝(Θ∣𝐃)𝑑Θ.

One can simulate from the posterior predictive distribution using the following two steps. First from the posterior distributions of the model parameters, simulate 𝑚 values (vectors) of Θ. Next for each value of Θ, simulate a value (vector) 𝐃^rep from the likelihood. The 𝑚 values (vectors) of 𝐃^rep drawn in this way will essentially come from posterior predictive distribution 𝑝(𝐃^rep∣𝐃).

We simulate 100 draws (𝑚=100) from the posterior predictive distribution and then apply proposed joint analysis and estimate the model parameters using MCMC as described earlier. We compute BF for each of those 7 linkage groups by considering the problem of testing the existence of no QTL (null) versus the existence of one QTL (alternative). Table 3 shows the average BF (with estimated SE) for each linkage group and the existence of QTL in groups 1, 20, and 23 is quite evident.

Table 3 

Posterior predictive check for 7 linkage groups.

Heritability (broad-sense) for the traits is estimated from the data, as the proportion of phenotypic variance attributable to genetic variance. The estimated heritability in our case is 32.6%. Also we compute the percentage of variance explained by three identified QTLs. It turns out the QTLs identified in linkage groups 1, 20, and 23 explain 6.8%, 14.3% and 11.4% of the total variance, respectively.

We show the marginal posterior plots with 95% credible intervals for the parameter 𝛾 in Figure 3. Note that for all three groups, the estimates and the confidence intervals are in the negative part which indicates a negative relationship between the trait and the event time. Biologically this is sensible since the plants with higher body mass will take less time to have the first flower compared to the plants with lower body masses.

(a)

(b)

(c)

Figure 3 

Marginal posterior plot for 𝛾 for linkage groups 1 (a), 20 (b), and 23 (c).

Figure 4 illustrates genotypic differences in whole-plant biomass trajectory and the time to first flower for the three QTLs detected. At the QTL on linkage group 1, the allele (𝑄) inherited from parent Nannong 1138-2 leads to increased biomass growth and earlier flowering than the allele (𝑞) from parent Kefeng no. 1. The inverse pattern is observed for the QTL on linkage group 20. The QTL on linkage group 23 alters its direction of genetic effect. Affected by the first two QTLs, poor vegetative biomass growth in plants stimulates early flowering (Figures 4(a) and 4(b)). Yet, the QTL on linkage group 23 makes the fast-growing genotype to flower earlier than the slow-growing genotype (Figure 4(c)).

(a)

(b)

(c)

Figure 4 

Whole-plant biomass growth trajectories and times to first flower (T1 and T2) for two different genotypes at each of the QTLs detected on linkage groups 1 (a), 20 (b), and 23 (c). Genotypes 𝑄𝑄 inherit two alleles from parent Nannong 1138-2, whereas genotype 𝑞𝑞 inherits two alleles from parent Kefeng no. 1.

4. Simulation Study

We performed simulation studies to study the statistical properties of the joint model. We assumed an RIL design of 200 progeny and simulated 11 evenly spaced markers on a linkage group of length 100 cM. A QTL is located at 43 cM from the very first marker of the linkage group. To reflect a practical problem, we used parameter estimates of the soybean QTL detected in linkage group 20 as true values to simulate the data, allowing the covariance structure. Time-dependent phenotypic values were assumed to follow a multivariate normal distribution and the event times were taken the same as the soybean data. To make a comparison, we analyzed the simulated data using our nonparametric GLM-based covariance structure and the traditional AR(1) and CS covariance structures.

The prior distributions for the model parameters were taken in the same way as discussed in Section 3.1. A uniform prior on (0,100) was considered for 𝐷^∗. For each situation, we ran Markov chains 120,000 iterations and initial 20,000 burn-in iterations were discarded. Model parameters were estimated from the posterior distributions on the basis of remaining 100,000 iterations. The computed BF was 0.649 giving a strong evidence against the null hypothesis. Table 4 shows the means of the Bayesian estimates of model parameters with their respective Monte Carlo standard errors, (MCSE). It can be seen that the estimates are quite close to the actual values with a reasonably small standard errors which justifies the accuracy and precision of our estimation procedure. However, our GLM-based approach provides better estimation of parameters than AR(1) and CS-based approaches.

Table 4 

Simulation results for genotypic-mean parameters and QTL locations under different covariance structures, AR(1)-, CS- and GLM-based approach.

Figure 5 elucidates the marginal posterior plot for the QTL location under three different covariance structures. It is found that both AR(1)- and CS-based models provide the peaks at wrong locations, whereas GLM-based nonparametric covariance structure locates QTL more accurately in which case the length of the credible interval is narrower than those obtained from the former two structures. This provides numerical evidence that the proposed GLM-based model has better precision of QTL localization.

(a)

(b)

(c)

Figure 5 

The Bayesian estimate of QTL location (indicated by dash vertical lines) from simulation studies under different covariance structures, AR(1) (a), CS (b), and nonparametric (c). The true QTL location is at 43 cM from the very first marker of the linkage group.

We perform further simulation studies to assess the reliability of BF in our data application. For each of the linkage groups 1, 20, and 23, we simulate data under the “null” model. As mentioned earlier, under the “null” model, 𝛽_𝟏=𝛽_𝟐. For group 1, we consider the null model 𝛽_𝟏=𝛽_𝟐=(−0.4762,3,3214,0.1548), for group 20 it is 𝛽_𝟏=𝛽_𝟐=(−1.1524,2.6829,0.2295), and for group 23, our null model is given as 𝛽_𝟏=𝛽_𝟐=(0.4190,0.5971,0.5438). The computed BFS for these three groups are 34.55, 46.79, and 41.86, respectively, suggesting strong evidence for the null.

5. Discussion

Tools to reveal the secret of life should reflect the dynamic nature of life. More recently, a series of statistical models have been developed to map quantitative trait loci (QTLs) that control the dynamic process of a complex trait [1, 2, 22, 28]. These so-called functional mapping models integrate mathematical aspects of biological processes into a statistical framework derived to map complex trait QTLs and have proved to be useful for detecting and identifying genes and genetic interactions involved in quantitative genetic variation for plant height, plant rooting ability, and animal body mass. Functional mapping is also flexible to incorporate complex biological phenomena, such as genotype-environment interactions and allometric scaling providing powerful means for addressing biological questions of fundamental importance.

In this paper, we develop a new version of functional mapping that can map QTLs for developmental events affected by organismic growth trajectories in time. This version is benefited from recent statistical developments for joint modeling of longitudinal traits and event time [16–18, 29]. In the presence of strong correlation between a longitudinal trait and event, a joint model performs better than submodels separately for a single trait. In the joint modeling framework for these two types of traits, we applied a GLM-based approach to model the covariance structure and local polynomials for the mean curves. Bayesian estimation method using the MCMC algorithm was used since it is computationally much simpler than a likelihood-based approach. Simulation results show the effectiveness of GLM-based covariance model compared to traditional parametric compound symmetry or autoregressive structure.

Our joint model, embedded within functional mapping, promotes the study of testing how QTLs pleiotropically affect different biological processes and how one trait is predicted by other traits through genetic information. The application of the new model to soybean mapping data does not only validate its usefulness and utilization, but also gains new insight into the genetic and developmental regulation of trait correlations in plants. There is no doubt that the new model can be modified to study the genetic associations between HIV dynamics and the time to death as well as prostate specific antigen change and the time to recurrence of prostate. However, there is much room for modifying this model. First, to clearly describe our idea, we assume one QTL at a time for trait control. Epistatic interactions between multiple QTLs may play an important role in trait development as well as in correlations between longitudinal traits and events. Second, from a dynamic systems perspective, we need to model dynamic correlations among multiple longitudinal traits and multiple events. Third, with the availability of efficient genotyping techniques, our model should accommodate a high-dimension model selection scheme to identify significant genetic variants from a flood of marker data.

Appendix

Denote 𝛽_−𝑗=(𝛽_𝑗′: 𝑗′=1,2;𝑗′≠𝑗), 𝑆(𝑡)=exp(−∫^𝑡₀𝜆(𝑢)𝑑𝑢), ℎ(𝑡)=(1,𝑡,𝑡²,…,𝑡^𝑟), and 𝑔(𝑡)=(1,𝑡,𝑡²,…,𝑡^𝑚). Let 𝑚_𝑗 be the number of progeny that carries QTL genotype 𝑗. We derive the full conditional distributions for the unknown parameters as𝜋𝜷_𝑗∣𝐲,𝐬,𝜎²,𝛾,𝐷,𝜼,𝜹,𝝀_𝟎,𝜷_−𝑗∝𝜋𝜷_𝑗𝜋𝐲,𝐬∣𝜷,𝜎²,𝛾,𝐷,𝜼,𝜹,𝝀_𝟎∝𝜋𝜷_𝑗𝜋𝐬∣𝐲,𝜷,𝜎²,𝛾,𝐷,𝜼,𝜹,𝝀_𝟎×𝜋𝐲∣𝜷,𝜎²,𝛾,𝐷,𝜼,𝜹,𝝀_𝟎−1∝exp2𝜷_𝑗𝚺_𝛽⁻¹𝜷^𝑇_𝑗−12^{𝑚𝑗𝑖=1}𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗^𝑇×𝚺⁻¹𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗×exp𝛾^𝑚𝑗_𝑖=1ℎ_𝑖𝑠_𝑖𝜷^𝑇_𝑗×^𝑚𝑗_𝑖=1𝑆_𝑖𝑠_𝑖−1=exp2𝜷_𝑗𝚺_𝛽⁻¹𝜷^𝑇_𝑗−12^{𝑚𝑗𝑖=1}𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗^𝑇×𝚺⁻¹𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗×exp𝛾^𝑚𝑗_𝑖=1ℎ_𝑖𝑠_𝑖𝜷^𝑇_𝑗×^𝑚𝑗_𝑖=1𝑆_𝑖𝑠_𝑖∝exp𝜷_𝑗𝚺_𝛽⁻¹𝜷^𝑇_𝑗+^𝑚𝑗_𝑖=1𝜷_𝑗𝐗_𝑖^(𝑟)𝑇𝚺⁻¹𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗−2^𝑚𝑗_𝑖=1𝐲^𝑇_𝑖𝚺⁻¹𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗exp𝛾^𝑚𝑗_𝑖=1ℎ_𝑖𝑠_𝑖𝜷^𝑇_𝑗×^𝑚𝑗_𝑖=1𝑆_𝑖𝑠_𝑖=exp𝜷_𝑗Σ_𝛽⁻¹+^𝑚𝑗_𝑖=1𝐗_𝑖^(𝑟)𝑇𝚺⁻¹𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗−2^𝑚𝑗_𝑖=1𝐲^𝑇_𝑖𝚺⁻¹𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗×exp𝛾^𝑚𝑗_𝑖=1ℎ_𝑖𝑠_𝑖𝜷^𝑇_𝑗^𝑚𝑗_𝑖=1𝑆_𝑖𝑠_𝑖. (A.1)

Hence, we have𝜋𝜷_𝑗∣.∶MVN𝛽_𝑗,𝑀_𝛽,𝑉_𝛽×exp𝛾^𝑚𝑗_𝑖=1ℎ_𝑖𝑠_𝑖𝜷^𝑇_𝑗^𝑚𝑗_𝑖=1𝑆_𝑖𝑠_𝑖, (A.2) where 𝑀_𝛽=𝚺_𝛽⁻¹+^𝑚𝑗_𝑖=1𝐗_𝑖^(𝑟)𝑇𝚺⁻¹𝐗_𝑖^{(𝑟)−1𝑚𝑗}_𝑖=1𝐗_𝑖^(𝑟)𝑇𝚺⁻¹𝐲_𝑖,𝑉_𝛽=𝚺_𝛽⁻¹+^𝑚𝑗_𝑖=1𝐗_𝑖^(𝑟)𝑇𝚺⁻¹𝐗_𝑖^(𝑟)−1. (A.3)

Similarly, the full conditional distributions for the other parameters can be derived as follows:𝜋𝛾∣⋅∝exp𝛾²_𝑗=1^𝑚𝑗_𝑖=1ℎ_𝑖𝑠_𝑖𝜷^𝑇_𝑗+𝑔_𝑖𝑠_𝑖𝜽^𝑇_𝑖×^𝑛_𝑖=1𝑆_𝑖𝑠_𝑖.𝜋𝜂∣⋅∝^𝑛_𝑖=1det(𝚺)^−(1/2)f−1×exp2𝜂^′𝚺_𝜂⁻¹1𝜂−2^{2𝑗=1𝑚𝑗𝑖=1}𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗×−𝐗_𝑖^(𝑚)𝜽^𝑇_𝑖^𝑇𝚺⁻¹𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗−𝐗_𝑖^(𝑚)𝜽^𝑇_𝑖.𝜋(𝜹∣⋅)∝^𝑛_𝑖=1det(𝚺)^−(1/2)−1×exp2𝛿^′𝚺_𝛿⁻¹1𝛿−2^{2𝑗=1𝑚𝑗𝑖=1}𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗×−𝐗_𝑖^(𝑚)𝜽^𝑇_𝑖^𝑇𝚺⁻¹𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗−𝐗_𝑖^(𝑚)𝜽^𝑇_𝑖.𝜋𝜎²∣⋅∝^𝑛_𝑖=1det𝚺+𝜎²𝐗_𝑖^(𝑚)𝐗_𝑖^{(𝑚)𝑇−(1/2)}×𝜎^{2−𝛼1−1}−𝛼exp²𝜎²−1×exp2^{2𝑗=1𝑚𝑗𝑖=1}𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗^𝑇×𝚺+𝜎²𝐗_𝑖^(𝑚)𝐗_𝑖^{(𝑚)𝑇−1}𝐲_𝑖−𝐗_𝑖^(𝑟)𝜷^𝑇_𝑗.𝜋𝜆_0𝑘∣⋅∝Gamma𝑑_𝑘+1,𝜙_𝜆⁻¹_0𝑘𝜋𝜆_0𝑘=Gamma𝑎+𝑑_𝑘+1,𝜙_𝜆0𝑘+1𝑏⁻¹, (A.4) where 𝑑_𝑘 is the number of events in the interval (𝑡^𝜆_𝑘−1,𝑡^𝜆_𝑘] and 𝜙_𝜆0𝑘=_{𝑖∶𝑠𝑖≥𝑡𝜆𝑘}^𝑡𝜆𝑘_{𝑡𝜆𝑘−1}exp𝛾ℎ(𝑡)𝛽^𝑇_𝑗+𝑔(𝑡)𝜃^𝑇_𝑖+𝑑𝑡_{𝑖∶𝑡𝜆𝑘−1<𝑠𝑖≤𝑡𝜆𝑘𝑠𝑖𝑡𝜆𝑘−1}exp𝛾ℎ(𝑡)𝛽^𝑇_𝑗+𝑔(𝑡)𝜃^𝑇_𝑖𝑑𝑡. (A.5)

Acknowledgments

This work is partially supported by NSF/IOS-0923975, Changjiang Scholars Award, “Thousand-Person Plan" Award, the China National Key Basic Research Program (2006CB1017, 2009CB1184, 2010CB125906), the China National Hightech R&D Program (2006AA100104), the Natural Science Foundation of China (30671266), and the China MOE 111 Project (B08025). R.Li was supported by an NIDA Grant P50-DA10075 and an NNSF of China Grant 11028103. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIDA, or the NIH.