Abstract
Background
The characterization of the relationship between a longitudinal response process and a timetoevent has been a pressing challenge in biostatistical research. This has emerged as an important issue in genetic studies when one attempts to detect the common genes or quantitative trait loci (QTL) that govern both a longitudinal trajectory and developmental event.
Results
We present a joint statistical model for functional mapping of dynamic traits in which the event times and longitudinal traits are taken to depend on a common set of genetic mechanisms. By fitting the Legendre polynomial of orthogonal properties for the timedependent mean vector, our model does not rely on any curve, which is different from earlier parametric models of functional mapping. This newly developed nonparametric model is demonstrated and validated by an example for a forest tree in which stemwood growth and the time to first flower are jointly modelled.
Conclusion
Our model allows for the detection of specific QTL that govern both longitudinal traits and developmental processes through either pleiotropic effects or close linkage, or both. This model will have great implications for integrating longitudinal and event data to gain better insights into comprehensive biology and biomedicine.
Background
Although there has been a upsurge of interest in jointly modelling longitudinal and event data during the last decade [19], no statistical models have been developed to characterize the shared genetic basis for these two types of traits. In biomedicine, the identification of specific genetic variants responsible for an HIV patient's timedependent CD4 count and for the time to onset of AIDS symptoms can help to design individualized drugs to control this patient's progression to AIDS. Similarly, in studies of prostate cancer, a shared genetic basis between prostate specific antigen, repeatedly measured for patients following treatment for prostate cancer, and the time to disease recurrence can be used to make optimal treatment schedules for patients. In plants, knowledge about whether the genetic loci for reproductive behaviors, such as the time to first flower and the time to form seeds, also govern growth rates and sizes of plants helps to understand the etiology of plant's adaptation to the environment in which they are grown.
The genetic mapping of quantitative trait loci (QTL) that are responsible for longitudinal traits has long been a difficult issue because of the dynamic features of these traits. More recently, part of this difficulty has been solved by integrating the statistical analysis of longitudinal data into a QTL mapping framework, leading to a socalled functional mapping strategy [1016]. Statistical models for functional mapping were established on the belief that biological processes can be described by mathematical functions. One of the most significant examples for this is the use of Sshaped logistic curves to model growth trajectories. West et al. [17] indicated from fundamental principles of biophysical processes that logistic forms of growth are biologically crucial for the maintenance of optimal metabolic level and, thereby, the best use of available resources for an organism from birth to adulthood. Because of the embedment of fundamental biological principles within the modelling model, functional mapping provides a quantitative framework for testing biologically relevant hypotheses at the interplay between gene actions and development.
The concept of functional mapping can be further extended to jointly mapping a longitudinal variable and a timetoevent by incorporating statistical theories developed to characterize the relationships between longitudinal response and event processes [19]. However, original functional mapping models for a dynamic trait reply upon explicit mathematical functions that describe the development of the trait. In practice, there are also many situations in which no appropriate curves can be used to describe a biological process. To model an arbitrary shape of curves, a different statistical model based on nonparametric theory should be formulated. Polynomial analyses that can be specified by varying orders have power to fit curves with arbitrary shapes. As shown by Kirkpatrick and Heckman [18], Legendre polynomials have several favorable properties for curve fitting which include: (1) the functions are orthogonal, (2) it is flexible to fit sparse data, (3) higher orders are estimable for high levels of curve complexity and (4) computation is fast because of good convergence.
The purpose of this article is to develop a joint statistical model for nonparametric functional mapping of longitudinal trajectories based on the Legendre polynomials, integrated with timetoevents. This joint model is constructed within the maximum likelihood context, including simultaneously modelling of the mean vector (based on nonparametric approaches) and covariance matrix (based on parametric approaches). By analyzing stem volume growth data in an example of a forest tree, we will demonstrates the implications of our joint model. Lastly, the advantages of our model in general biomedical and biological research and the areas in which the model can be further refined are discussed.
The Model
The likelihood function
Consider a mapping population of size n for which a number of molecular markers are genotyped, aimed to identify QTL for a longitudinal trait and timetoevent. Every individual of the mapping population is measured for the longitudinal trait at multiple (say T) time points (y) and a timetoevent (z). Variable z can be the time to first flower, the timing of cancer malignance, the time of mortality, or the events that happen at a time. The inference of unknown QTL genotypes for the phenotypic traits based on observed marker information (M) can be made due to cosegregation between the QTL and markers.
Suppose there are two segregating QTL for longitudinal and event traits in the mapping population, each with genotypes 2, 1 and 0. These two QTL are assumed to be linked or associated with and, therefore, can be inferred from, markers. The joint likelihood function of the two types of phenotypic data and marker information at the two underlying QTL is written as
where Ω is the unknown vector that defines the QTL positions, timedependent QTL effects ( and ) and covariance matrix (Σ), is the mixture proportion expressed as the conditional probability of a joint genotype j_{1}j_{2 }for the longitudinal and event QTL given marker genotypes for individual i and is the (T + 1)dimensional multivariate normal distribution function with mean vector (, ) and covariance matrix Σ.
Conditional probabilities
There are different descriptions of the conditional probability, depending on the type of the mapping population. If the mapping population is an experimental cross initiated with two contrasting parents, such as the F_{2 }or backcross, the conditional probability is described in terms of the recombination fractions between the markers and QTL [12,19]. If the two QTL are bracketed by different pairs of markers, the conditional probability of joint QTL genotypes given the marker intervals can be expressed as the product of the corresponding conditional probabilities for QTL genotypes given a single marker interval. If the two QTL are located at the same marker interval, the conditional probabilities should be derived using the principle of 4point analysis. For a natural population, the association between the QTL and markers can be described by the coefficients of linkage disequilibria [17]
Modelling the mean vector
The choice of a mean function for a longitudinal trait is based on theory or past experience that suggests a certain mathematical form for the timedependent mean. However, it would be essential to derive a general approach that can fit any kind of curves. By choosing different orders of orthogonal polynomials, the Legendre function has potential to approximate the functional relationships between trait values and times to any specified degree of precision. The Legendre polynomials are solutions to a very important differential equation, the Legendre equation,
The polynomials may be denoted by P_{r}(x), called the Legendre polynomial of order r. The polynomials are either even or odd functions of x for even or odd orders r.
The general form of a Legendre polynomial of order k is given by the sum,
where K = r/2 or (r  1)/2 whichever is an integer. This polynomial is defined over the interval [1, 1]. From Eq. 5, we show the first few polynomials as
In this modelling, independent variable x is expressed as time t, which is adjusted, to rescale the measurement times to the range of the orthogonal function [1, 1], by
where t_{min }and t_{max }are respectively the first and last time points.
Our aim is to model the timedependent genotypic values for different QTL genotypes j_{1}j_{2}, using the orthogonal Lengedre polynomial with a particular order r. A family of such polynomials is denoted by
and a vector of genotypic values, which is timeindependent, denoted by
The timedependent genotypic values (t) can be described as a linear combination of weighted by the family of the polynomials, i.e.,
Substituting the mean vector of the likelihood (1) by the above expression (3), we will need to estimate timeinvariant genotypic values for the longitudinal trait, , and the genotypic mean for the event trait, .
Modelling the covariance matrix
A general form for the covariance matrix among longitudinal trajectories and development event in the likelihood (1) is expressed as
where Σ_{y }and are the covariance matrix and variance for the longitudinal and event traits, respectively, and Σ_{yz }= is the covariance matrix between these two types of traits. The structures of Σ_{y }and Σ_{yz }can be empirically modelled on the basis of prior knowledge or results. Several approaches for parametric modelling of the covariance matrix, reviewed by Zimmerman and NunezAnton [20], can be utilized.
The most common approach for modelling the covariance structure is based on a variancecorrelation specification, in which functions for the responses' variances and correlations are specified. In previous QTL mapping [1014], the covariance structure for longitudinal traits is modelled by the simplest, most parsimonious and most flexible firstorder autoregressive (AR(1)) model in which there are two parameters, stationary variance () and correlation (ρ_{y}). Relaxing the stationary variance assumption for growth data, Wu et al. [15] adopted a transformbothsides (TBS) model to obtain an empirically homogeneous variance. The results from simulation studies suggest that the TBSbased mapping model provides more precise estimates for curve parameters and residual variancecorrelation than the untransformed model.
The TBSbased model displays the potential to relax the assumption of variance stationarity, but the covariance stationarity issue remains unsolved. Zimmerman and NúñezAntón [20] proposed a socalled structured antedependence (SAD) model to model the agespecific change of correlation in the analysis of longitudinal traits. The SAD model has been employed in several studies and displays many favorable properties [21,22].
The emergence of a developmental event (z) at time t* can be correlated with the longitudinal trait. For example, larger tumor sizes may be likely to lead to earlier malignance of cancer than smaller tumor sizes. An AIDS patient would die when his/her HIV load accumulates to a particularly high level. In plants, first flowering only appears after some investment of vegetative growth. All such common knowledge suggests that the correlation between the event trait at time t* and longitudinal trait measured at time t (before t*) decays with time difference (t*  t). In fact, a similar pattern of correlation should also hold for t > t* because of the autocorrelation nature. With all this consideration, the correlation between the event and longitudinal traits can be modelled by the power equation, expressed as
where 0 ≤ η ≤ 1. Equation (5) suggests that the event is correlated with the longitudinal trait, to the same extent, before and after its emergence. The event trait should be individualspecific when it is the timing of development, such as the time to first flower. In this case, Equation (5) and, therefore, the covariance matrix (4), expressed as Σ_{i}, should be individualspecific. If Σ_{y }is modelled by the AR(1) model, one can derive the explicit expressions of the determinant and inverse of Σ_{i }specified by (, ρ_{y}, , η, λ).
Computational algorithms
The unknown parameters (Ω) contained with the mixture model (1) include three types, QTLmarker recombination fractions for a pedigree or QTLmarker linkage disequilibria for a natural population reflected in the conditional probabilities , the curve parameters (, ) that model the mean vector, and the parameters (, ρ_{y}, , η, λ) that model the structure of the covariance matrix. We derived the EM algorithm to estimate these parameters. Using the (prior) conditional probability and the likelihood, we define the posterior probability for individual i to bear on a QTL genotype j_{1}j_{2 }as
The posterior probabilities are then used to derive a closedform maximum likelihood estimates of the QTL locations, expressed as the ratio of recombination fractions, for linkage analysis or QTLmarker haplotype frequencies for linkage disequilibrium analysis [17]. For functional mapping, in which the mean vectors and covariance matrix are modelled by mathematical parameters based on nonlinear equations, it is impossible to derive the closed forms for these parameters, the simplex algorithm, widely used in operations research, is found to provide a fast and precise estimation of the curve parameters and the parameters that model the residual covariances [23]. Thus, we implement the simplex algorithm in the maximization process of the EM algorithm.
For linkage analysis based on an experimental cross, 's are expressed in the recombination fraction between the QTL and two flanking markers. In practical computations, the QTL position parameter can be viewed as a fixed parameter because a putative QTL can be searched at every 1 or 2 cM on a map interval bracketed by two markers throughout the entire genome. The amount of support for a QTL at a particular map position is often displayed graphically through the use of likelihood maps or profiles, which plot the likelihood ratio test statistic as a function of map position of the putative QTL. The peak of the profile corresponds to the position of the QTL over the genome.
For linkage disequilibrium analysis of a natural population, we have derived a closed form for the EM algorithm to estimate QTLmarker haplotype frequencies. From the estimated haplotype frequencies, the allele frequencies of QTL and QTLmarker linkage disequilibria can be estimated. How the markers are associated with the underlying QTL in the population can be tested for the significance of QTLmarker linkage disequilibria.
After the point estimates of parameters are obtained by the EM algorithm, the approximate variancecovariance matrix and the sampling errors of the estimates () can be estimated. The techniques for so doing involve calculation of the incompletedata information matrix which is the negative secondorder derivative of the incompletedata loglikelihood. The incompletedata information can be calculated by extracting the information for the missing data from the information for the complete data [24].
Order selection
For a QTL to be detected, we need to determine the optimal order for the Legendre polynomial that fits the data. We propose using the AIC information criterion to select the best model. The AIC value at a particular order, r, is calculated by
AIC = 2 ln L(r) + 2 dimension(Ωr), (7)
where (r) is the the MLE of parameters for the Legendre polynomial of order r and dimension (Ωr) represents the number of independent parameters under order r.
Also, Bayesian Information Criterion (BIC) [25] is used to determine the optimal order of the Legendre function, which is calculated by
BIC = 2 ln L(r) + 2 dimension(Ωr) ln (nT). (8)
As compared to AIC, BIC adjusts the effects of sample size and the number of time points measured.
Hypothesis tests
Our model allows for a number of hypothesis tests to examine the genetic control of growth processes [14]. All these tests are helpful to address biological questions related to the genetic control mechanisms of growth. Testing whether specific QTL exist to affect the longitudinal and event processes is a first step toward the understanding of the detailed genetic architecture of complex phenotypes. This can be tested by formulating the following hypotheses,
H_{0 }: = u and = v vs. H_{1}: Not all equalities in H_{0 }hold. (9)
The H_{0 }states that there is no QTL affecting longitudinal and event processes (the reduced model), whereas the H_{1 }proposes that such a QTL does exist (the full model). The test statistic for testing the hypotheses is calculated as the loglikelihood ratio of the reduced to the full model:
LR = 2[ln L_{0}(y, z)  ln L_{1 }(y, z, M)], (10)
where and denote the MLEs of the unknown parameters under H_{0 }and H_{1}, respectively. Because the LR calculated by equation (9) may not be asymptotically χ^{2}distributed with eleven degrees of freedom due to violation of regularity conditions, an empirical approach for determining the critical threshold based on permutation tests is used. By repeatedly shuffling the relationships between marker genotypes and phenotypes, a series of the maximum loglikelihood ratios are calculated, from the distribution of which the critical threshold is determined.
After the QTL are detected to be significant for both longitudinal and event traits, we need to test whether the detected QTL are significant separately for each trait. We assume two different genetic settings:
(1) Longitudinal and event traits are under control of the same QTL;
(2) Each process is controlled by different QTL that are linked on the same chromosomal region.
For the first setting, only after it is significant in two separate tests for longitudinal and event traits can the tested QTL be thought to be pleiotropic in affecting both types of traits. For the second setting, we assume that two tested QTL are located in the same interval bracketed by two markers. The comparison of the first and second setting can examine how the detected QTL jointly affect the differentiation in longitudinal and event traits. First, we can test how two genetic mechanisms, pleiotropy or close linkage, contribute to the correlation between these two types of traits. If the two QTL are detected to be significant for both, we then test whether such a correlation is due to pleiotropy or close linkage. Second, when two QTL exist, we can test how they epistatically interact to affect longitudinal trajectories and developmental events. Wu et al. [14] formulated a procedure for testing the epistatic effects on developmental trajectories.
Results
The proposed joint model is used to analyze growth trajectories and flowering behavior in a forest tree. The study material used was derived from the interspecific hybridization of Populus (poplar), P. deltoides and P. eummericana. This hybrid population was planted at a spacing of 4 × 5 m in the complete randomized design in a field trial near Xuzhou City, Jiangsu Province, China. The total stem heights and diameters measured at the end of each of the first 11 growing seasons are used to calculate stem volume indices (y) for QTL analysis. Because the vegetative and reproductive growth processes are generally correlated in plants [26], the ages to first flower (z) was predicted by a regression equation for each of these hybrids. Two genetic linkage maps each based on a different parent were constructed for a subset of hybrids (90) with different types of molecular markers that are segregating in a pattern of pseudotest backcross [27]. Our analysis here will be based on P. deltoides (D)specific linkage map.
Although stem height and diameter for each tree follows a logistic curve [10], the stem volume index derived from these two traits cannot be fit by the growth equation mainly because stem volume has not yet reached its asymptotic growth during this measurement period (Fig. 1A). As shown by Figure 1, the variance of the stem volume index increases markedly with age, but the logtransformation of these indices leads to much parallel curves (Fig. 2B), suggesting that the variance stationarity assumption may be met after the transformation.
Figure 1. Plots of stem volume index growth vs. ages for each of the 90 genotypes used to construct linkage maps in poplar hybrids (Yin et al. 2002). The relationships between growth and age are displayed for untransformed (A) and logtransformed data (B).
Figure 2. The profile of the loglikelihood ratios between the full (there is a QTL) and reduced (there is no QTL) model that combines stem volume index growth trajectories and flower timing across linkage groups in the Populus deltoides parent map. The genomic positions corresponding to the peaks of the curve are the MLEs of the QTL localization. The threshold values for claiming the existence of QTL are given as the horizonal solid lines for the genomewide level and broken lines for the chromosomewide level. Blue color corresponds to the unifying model for jointly mapping growth trajectories and flower trait, whereas red color corresponds to a model for mapping growth trajectories only. The positions of markers on the linkage groups (Yin et al. 2002) are indicated at ticks.
We implemented the Legendre function to model the QTL genotypic mean vector of growth trajectories and the TBSbased AR(1) model to approximate the structure of the covariance matrix. The joint model also allows the estimation of the genotypic means and residual variance for the age to first flower, as well as the correlation of this trait with stem wood growth trajectories. Equation (5) provides a general equation for modelling the correlation between the event and longitudinal traits measured at different ages. In this example, it was observed that there were significant correlations between the age to first flower and volume growth at all different ages (0.29 – 0.82). Thus, for simplicity of computation, we assume that such correlations are consistent across the ages of volume index, denoted as η. Using the adjusted ages, we calculated the coefficients of the Legendre polynomials for the first seven orders (Table 1). These coefficients are used to estimate timedependent genotypic values. The AIC and BIC values calculated consistently suggested an optimal order of 5 to fit stem volume growth (Table 2).
Table 1. Coefficients of the first five Legendre polynomials for adjusted time points (t*) used in the poplar growth study.
Table 2. The AIC and BIC values used to determine the optimal order for the Lengendre polynomials.
While our joint model was derived to detect two QTL at a time, it was reduced to a oneQTL model because of a limited sample size for sufficient estimates of twoQTL model parameters. For the pseudotest backcross there are two genotypes, Qq (j = 1) and qq (j = 0), at each QTL. Figure 2 illustrates the profile of the loglikelihood ratio (LR) values for testing the existence of QTL that control either overall growth curves of stem volume indices from age 0 to 10 years or the ages to first flower, or both, across all of the 19 Dspecific linkage groups. We performed 100 permutation tests to determine critical threshold values for declaring the existence of QTL. By comparing the peaks of the LR profile with the thresholds, three significant QTL were detected, one on linkage group 2 at the 5% genomewide testing level and two on linkage groups 5 and 12 at the 5% chromosomewide testing level (Fig. 2; Table 3). We indicated the positions of these QTL on linkage groups, which correspond to the peaks of the LR profile. If only the stem volume growth is analyzed using traditional functional mapping [10], only the QTL on linkage group 2 is detected, suggesting that the joint model displays better power than a singletrait analysis.
Table 3. The MLEs and their sampling errors (SE, in the parentheses) of the QTL position, timeinvariant QTL effects on growth curves (expressed in the Legendre polynomials), QTL effect on the time to first flower, residual variance and residual correlation under the logtransformed model for the interspecific poplar hybrid mapping population.
Each of the three QTL was tested for their pleiotropic effect on both vegetative growth and reproduction by formulating two independent null hypotheses, one being that the QTL does not affect stem growth and the second being that the QTL does not affect flowering age. The rejection of both the null hypotheses implies that a QTL has a pleiotropic effect on growth and reproduction. As indicated by Table 3, all the detected QTL on linkage groups 2, 5 and 12 only trigger a significant effect on stem volume growth, but neither has an effect on both growth and reproduction.
The MLEs of growth parameters for stem volume indices, covariancestructuring parameters and the parameters dealing with reproductive behaviors, as well as their standard errors estimated from the Fisher information matrix, were tabulated in Table 3. It can be seen that the estimates of all the parameters from our joint model provide reasonable precision, using the estimates of growth curves, we draw two different curves each corresponding to a genotype at each of the detected QTL (Fig. 3). Note that growth curves were first drawn from the estimates of the Legendre parameters (the left panel of Fig. 3) and then transformed back to the normal scale (the right panel of Fig. 3). In general. these QTL are switched on to affect the overall stem growth process after age 4–5 years at which strong intertree competition sets in the stand due to canopy closure. Figure 3 also displays genotypic differences in the age to first flower at each of the growth QTL. But as tested, only QTL on linkage group 12 has a significant impact on the age to first flower (Table 3). At this QTL, the slowergrowing genotype flowers about 0.7 year earlier than the fastergrowing genotype. Through this QTL, the fastgrowing attribute and the capacity to efficiently occupy growth resources can be transmitted to the next generation.
Figure 3. Volume growth curves for two different QTL genotypes for the QTL detected on linkage group 2 by the Legendre polynomialbased model. Left panel: logtransformed curves; Right panel: antetransformed curves. Growth trajectories for all the individuals studied are indicated in yellow background. The effect of the detected QTL on the time to first flower is indicated.
Discussion
A theoretical framework has been constructed for functional mapping of quantitative trait loci (QTL) underlying longitudinal growth [1015]. Functional mapping was grounded on biological reality that every organism follows universal growth laws that can be derived from fundamental principles for the allocation of metabolic energy between maintenance of existing tissue and the production of new biomass [17]. In a couple with linkage disequilibrium mapping, functional mapping has been extended to map host QTL for HIV dynamics for a natural human population [16].
Although functional mapping has proven to be both biologically and statistically advantageous in terms of the estimates of the QTL positions and effects, its practical applications may be limited for two reasons. First, a longitudinal variable, such as HIV dynamics, tumor growth or plant vegetative growth, may be related to timetoevents, like time to onset of AIDS symptoms, time to first malignancy or time to first flower, through a common set of QTL [28,29]. Second, not all longitudinal data measured at a series of discrete time points can be fit by a mathematical function with biological means.
In this article, we have proposed a joint model for functional mapping of longitudinal trajectories and timetoevents with the nonparametric context. Several statistical models have been proposed to jointly analyze longitudinal and event processes [19]. Different from those traditional models, our joint model has been constructed within the mixture model framework, with each mixture component assigned by biological rationale. We incorporated Legendre polynomials to characterize an arbitrary form of growth curves. In a real example for a forest tree, the model has detected a few QTL that affect growth processes and the age to first flower. The detection of the common genetic basis for vegetative and reproductive growth supports the views that any developmental event is not isolated from the growth process [28,29]. Our model provides a complete genetic analysis of growth courses for various organisms at different organization levels. From a statistical perspective, it increases the power of QTL detection and the precision of parameter estimation because the information about growth and development is jointly utilized. Meanwhile, our model allows for the test of several important hypotheses regarding the genetic control of developmental events occurring from fertilized ovum to reproductive maturity.
Our model is based on nonparametric Legendre orthogonal polynomial approaches for growth and development processes. Orthogonal polynomials (including Legendre) have been extensively used in random regression analyses for longitudinal traits with repeated records [18,3032]. There are several favorable properties for Legendre polynomials to be utilized in curve fitting, i.e., (1) the functions are orthogonal, (2) it is flexible to fit sparse data, (3) higher orders are estimable for high levels of curve complexity and (4) computation is fast because of good convergence. Nonparametric regression methods for modelling the mean structure of longitudinal data have been based on more commonly used Bspline basis functions [33]. Brown et al. [9] extended the Bspline basis to model multiple longitudinal variables. As compared to the Bspline approach that constructs curves from pieces of lower degree polynomials smoothed at selected pointed (knots), Legendre polynomials are simpler in which only fewer regression coefficients are needed to model the curve. However, polynomials often overemphasize the observations at the extremes and may be problematic for high orders of fit due to oscillations at the extremes of the curve [34]. It is therefore worthwhile implementing more flexible Bspline basis functions into the nonparametric functional mapping model.
In our joint model, we assumed that the timetoevent is multivariate normally distributed together with longitudinal data (see also [35]). An alternative to model the distribution of longitudinal and event data is to take the product of the normal distribution function of longitudinal trajectories and the distribution of the event trait and sensoring indicator given the trajectory function [9,36]. In addition, our model should be extended to consider multiple longitudinal variables based on a framework by Lin et al. [5], multiple timetoevents [8] and structured covariance matrices among unbalanced repeatedmeasures [37]. In order to unravel the genetic architecture of complex phenotypes that are characterized by a network of biological processes, such extensions will be essential. With appropriate improvements, our joint model will have great power to unlock the genetic secrets hidden in various complicated and biologically realistic life processes.
Conclusion
We have developed a joint statistical model that can detect specific QTL governing both longitudinal traits and developmental processes through either pleiotropic effects or close linkage, or both. This model was integrated by nonparametric approaches that do not rely on mathematical equations to model growth curves. The model will have great implications for integrating longitudinal and event data to gain better insights into comprehensive biology and biomedicine.
Authors' contributions
ML derived the models, programmed the method and performed data analyses. RW conceived the idea and drafted the manuscript.
Figure 4. Volume growth curves for two different QTL genotypes for the QTL detected on linkage group 5 by the Legendre polynomialbased model. See Figure 3 for all the explanations.
Figure 5. Volume growth curves for two different QTL genotypes for the QTL detected on linkage group 12 by the Legendre polynomialbased model. See Figure 3 for all the explanations.
Acknowledgements
We thank the two anonymous referees for their constructive comments on this manuscript. The preparation of this manuscript has been partially supported by NSF grant (0540745) and NIH grant (R01 NS041670) to R. W. The publication of this manuscript is approved as journal series No.10579 by the Florida Agricultural Experimental Station.
References

Tsiatis AA, DeGruttola V, Wulfsohn MS: Modeling the relationship of survival to longitudinal data measured with error: Applications to survival and CD4 counts in patients with AIDS.
Journal of the American Statistical Association 1995, 90:2737. Publisher Full Text

Wulfsohn MS, Tsiatis AA: A joint model for survival and longitudinal data measured with error.
Biometrics 1997, 53:330339. PubMed Abstract

Henderson R, Diggle P, Dobson A: Joint modeling of longitudinal measurements and event time data.
Biostatistics 2000, 4:465480. Publisher Full Text

Song X, Davidian M, Tsiatis AA: A semiparametric likelihood approach to joint modeling of longitudinal and timetoevent data.
Biometrics 2002, 58:742753. PubMed Abstract  Publisher Full Text

Lin HQ, McCulloch CE, Mayne ST: Maximum likelihood estimation in the joint analysis of timetoevent and multiple longitudinal variables.
Statistics in Medicine 2002, 21:23692382. PubMed Abstract  Publisher Full Text

Lin HQ, Turnbull BW, McCulloch CE, Slate EH: Latent class models for joint analysis of longitudinal biomarker and event process data.
Journal of the American Statistics Association 2002, 457:5365. Publisher Full Text

Tsiatis AA, Davidian M: Joint modeling of longitudinal and timetoevent data: An overview.

Chi YY, Ibrahim JG: Joint models for multivariate longitudinal and multivariate survival data.

Brown EB, Ibrahim JG, DeGmttola V: A flexible Bspline model for multiple longitudinal biomarkers and survival.
Biometrics 2005, 61:6473. PubMed Abstract  Publisher Full Text

Ma CX, Casella G, Wu RL: Functional mapping of quantitative trait loci underlying the character process: A theoretical framework.
Genetics 2002, 161:17511762. PubMed Abstract  Publisher Full Text

Wu RL, Ma CX, Chang M, Littell RC, Wu SS, Huang M, Wang M, Casella G: A logistic mixture model for characterizing genetic determinants causing differentiation in growth trajectories.
Genetical Research 2002, 19:235245. Publisher Full Text

Wu RL, Ma CX, Zhao W, Casella G: Functional mapping of quantitative trait loci underlying growth rates: A parametric model.
Physiological Genomics 2003, 14:241249. PubMed Abstract  Publisher Full Text

Wu RL, Ma CX, Yang MCK, Chang M, Santra U, Wu SS, Huang M, Wang M, Casella G: Quantitative trait loci for growth in Populus.
Genetical Research 2003, 81:5164. PubMed Abstract  Publisher Full Text

Wu RL, Ma CX, Lin M, Casella G: A general framework for analyzing the genetic architecture of developmental characteristics.
Genetics 2004, 166:15411551. PubMed Abstract  Publisher Full Text

Wu RL, Ma CX, Lin M, Wang ZH, Casella G: Functional mapping of growth QTL using a transformbothsides logistic model.
Biometrics 2004, 60:729738. PubMed Abstract  Publisher Full Text

Wang ZH, Wu RL: A statistical model for highresolution mapping of quantitative trait loci determining human HIV1 dynamics.
Statistics in Medicine 2004, 23:30333051. PubMed Abstract  Publisher Full Text

West GB, Brown JH, Enquist BJ: A general model for ontogenetic growth.
Nature 2001, 413:628631. PubMed Abstract  Publisher Full Text

Kirkpatrick M, Heckman N: A quantitative genetic model for growth, shape, reaction norms, and other infinitedimensional characters.
Journal of Mathematical Biology 1989, 27:429450. PubMed Abstract  Publisher Full Text

Lander ES, Botstein D: Mapping Mendelian factors underlying quantitative traits using RFLP linkage maps.
Genetics 1989, 121:185199. PubMed Abstract  Publisher Full Text

Zimmerman DL, NúñezAntón V: Parametric modeling of growth curve data: An overview (with discussions).

Jaffrézic F, Thompson R, Hill WG: Structured antedependence models for genetic analysis of repeated measures on multiple quantitative traits.
Genetical Research 2003, 82:5565. PubMed Abstract  Publisher Full Text

Zhao W, Chen YQ, Casella G, Cheverud JM, Wu RL: A nonstationary model for functional mapping of complex traits.
Bioinformatics 2005, 21:24692477. PubMed Abstract  Publisher Full Text

Zhao W, Wu RL, Ma CX, Casella G: A fast algorithm for functional mapping of complex traits.
Genetics 2004, 167:21332137. PubMed Abstract  Publisher Full Text

Louis TA: Finding the observed information matrix when using the EM algorithm.
Journal of the Royal Statistics Society Series B 1982, 44:226233.

Kozlowski TT, Pallardy SG: Acclimation and adaptive responses of woody plants to environmental stresses.

Yin TM, Zhang XY, Huang MR, Wang MX, Zhuge Q, Tu SM, Zhu LH, Wu RL: The molecular linkage maps of the Populus genome.
Genome 2002, 45:541555. PubMed Abstract  Publisher Full Text

Ambros V: Control of developmental timing in Caenorhabditis elegans.
Current Opinion in Genetics and Development 2000, 10:42833. Publisher Full Text

Rougvie AE: Control of developmental timing in animals.
Nature Reviews Genetics 2001, 2:690701. PubMed Abstract  Publisher Full Text

Schaeffer LR: Application of random regression models in animal breeding.
Livestock Production Science 2004, 86:3545. Publisher Full Text

Meyer K: Estimates of genetic covariance functions for growth of Angus cattle.
Journal of Animal Breeding and Genetics 2005, 122:7385. Publisher Full Text

Meyer K: Random regression analyses using Bsplines to model growth of Australian Angus cattle.
Genetics Selection Evolution 2005, 37:473500. Publisher Full Text

Rice JA, Wu CO: Nonparametric mixed effects models for unequally sampled noisy curves.
Biometrics 2001, 57:253259. PubMed Abstract  Publisher Full Text

de Boor C: A Practical Guide to Splines. 2nd edition. SpringerVerlag; 2001.

Degmttola V, Tu XM: Modeling progression of CD4lymphocyte count and its relationship to survivaltime.
Biometrics 1994, 50:10031014. PubMed Abstract

JacqminGadda H, Thiebaut R, Chene G, Commenges D: Analysis of leftcensored longitudinal data with application to viral load in HIV infection.
Biostatistics 2000, 1:355368. PubMed Abstract  Publisher Full Text

Jennrich RI, Schluchter MD: Unbalanced repeatedmeasures models with structured covariance matrices.
Biometrics 1986, 42:805820. PubMed Abstract