|Home | About | Journals | Submit | Contact Us | Français|
Biometric Research Branch, Division of Cancer Treatment and Diagnosis, National Cancer Institute, Bethesda, Maryland 20892, USA
In many medical studies, patients are followed longitudinally and interest is on assessing the relationship between longitudinal measurements and time to an event. Recently, various authors have proposed joint modeling approaches for longitudinal and time-to-event data for a single longitudinal variable. These joint modeling approaches become intractable with even a few longitudinal variables. In this paper we propose a regression calibration approach for jointly modeling multiple longitudinal measurements and discrete time-to-event data. Ideally, a two-stage modeling approach could be applied in which the multiple longitudinal measurements are modeled in the first stage and the longitudinal model is related to the time-to-event data in the second stage. Biased parameter estimation due to informative dropout makes this direct two-stage modeling approach problematic. We propose a regression calibration approach which appropriately accounts for informative dropout. We approximate the conditional distribution of the multiple longitudinal measurements given the event time by modeling all pairwise combinations of the longitudinal measurements using a bivariate linear mixed model which conditions on the event time. Complete data are then simulated based on estimates from these pairwise conditional models, and regression calibration is used to estimate the relationship between longitudinal data and time-to-event data using the complete data. We show that this approach performs well in estimating the relationship between multivariate longitudinal measurements and the time-to-event data and in estimating the parameters of the multiple longitudinal process subject to informative dropout. We illustrate this methodology with simulations and with an analysis of primary biliary cirrhosis (PBC) data.
Recently, many studies collect longitudinal data on a panel of biomarkers, and interest is on assessing the relationship between these biomarkers and time to an event. For example, Allen et al. (2007) examined the relationship between five longitudinally collected cytokines measured from serum plasma and survival. Interest focused on whether the values of these multiple cytokines are associated with survival. In another example, patients with primary biliary cirrhosis are followed longitudinally and interest is on examining whether multiple longitudinally biomarkers are prognostic for a poor clinical outcome. Important features in studies of this type are that there may be a relatively large number of biomarkers and that these biomarkers are subject to sizable measurement error due to laboratory error and biological variation.
Various authors have proposed joint modeling approaches for a single longitudinal measurement and time-to-event data [Tsiatis, DeGruttola and Wulfsohn (1995); Wulfsohn and Tsiatis (1997); Tsiatis and Davidian (2004); Henderson, Diggle and Dobson (2000), among others]. There is also limited work on joint models for a few longitudinal measurements and time-to-event data [Xu and Zeger (2001a, 2001b); Huang et al. (2001); Song, Davidian and Tsiatis (2002); Ibrahim, Chen and Sinha (2004); Brown, Ibrahim and DeGruttola (2005); Chi and Ibrahim (2006)]. However, these methods are difficult to implement when the number of longitudinal biomarkers is large since most of these approaches require integrating over the vector of all random effects to evaluate the joint likelihood of the multivariate longitudinal and time-to-event data. This paper proposes an approach for jointly modeling multivariate longitudinal and discrete time-to-event data which easily accommodates many longitudinal biomarkers.
Fieuws and Verbeke (2005) and Fieuws, Verbeke and Molenberghs (2007) have proposed an approach for modeling multivariate longitudinal data whereby all possible pairs of longitudinal data are separately modeled and are then combined in a final step. We use a similar approach along with a recent regression calibration approach for jointly modeling a single series of longitudinal measurements and time-to-event data [Albert and Shih (2009)] to implement the joint modeling approach proposed in this paper. Recently, Fieuws et al. (2008) have proposed a discriminant analysis based approach for using multivariate longitudinal profiles to predict renal graft failure. At the end of their discussion, they mention that a more elegant approach, which has not yet been developed, would involve a joint model for the many longitudinal profiles and time-to-event data. This paper presents such an approach.
We describe the approach in Section 2. We show the advantages of this approach using simulation in Section 3. We illustrate the methodology with an analysis of primary biliary cirrhosis data (PBC) in which we simultaneously examine the relationship between multiple longitudinal biomarker in Section 4. A discussion follows in Section 5.
Define Ti to be a discrete event-time which can take on discrete values tj, j = 1, 2, …, J, and Yij to be a binary indictor of whether the ith patient is dead at time tj. Then , where indicates the number of follow-up measurements before the event or the end of follow-up at time tJ. Longitudinal measurements are measured at times t1, t2, …, tJi. Denote X1i = (X1i1, X1i2, …, X1iJi)′, X2i = (X2i1, X2i2, …, X2iJi)′, …, XPi = (XPi1, XPi2, …, XPiJi)′ as the P biomarkers measured repeatedly at j = 1, 2, …, Ji time points. Further, define as the longitudinal measurements without measurement error for the pth biomarker and . We consider a joint model for multivariate longitudinal and discrete time-to-event data in which the discrete event time distribution is modeled as a linear function of previous true values of the biomarkers without measurement error on the probit scale. Specifically,
where i = 1, 2, …, I, j = 2, 3, …, Ji, Yi1 is taken as 0, α0j governs the baseline discrete event time distribution and αp measures the effect of the pth biomarker (p = 1, 2, …, P) at time tj–1 on survival at time tj. Specifically, (1) allows for examining the effect of multiple “true” biomarker values at time j – 1 on the probability of an event between the (j – 1) th and jth time point.
The longitudinal data is modeled assuming that the fixed and random effect trajectories are linear. Specifically, the multivariate longitudinal biomarkers can be modeled as
where βp0 and βp1 are the fixed effect intercept and slope for the pth biomarker, and γpi0 and γpi1 are the random effect intercept and slope for the pth biomarker on the ith individual. Denote β = (β10, β11, β20, β21, …, βP0, βP1)′ and γi = (γ1i0, γ1i1, γ2i0, γ2i1, …, γPi0, γPi1)′. We assume that γi is normally distributed with mean 0 and variance Σγ, where Σγ is a 2P by 2P dimensional variance matrix, and εpij are independent error terms which are assumed to be normally distributed with mean 0 and variance .
Alternative to (1), where the probability of an event over an interval depends on the true biomarker values at the beginning of the interval, the event-time process could be adapted to depend on the random effects of the multivariate longitudinal process [e.g., γpi1 can replace in (1)].
where rij = P (Yij = 1|Yi(j − 1) = 0), h (Xpi|γpi0, γpi1) is the product of Ji univariate normal density functions each with mean and variance , and f (γ) is a multivariate normal density with mean zero and variance Σγ. When P = 1, (4) can be maximized by numerical integration techniques such as a simple trapezoidal rule or Gaussian quadrature [Abramowitz and Stegun (1974)]. However, these methods are not feasible for even a few longitudinal biomarkers. Alternative Monte Carlo methods such as Monte Carlo EM [Wei and Tanner (1990)] are possible, but these methods do not perform well for even moderately high dimensional random effects (say, P > 2). In the next subsection we develop an alternative approach which is easy to implement with a large number of longitudinal biomarkers.
We propose a two stage regression calibration approach for estimation, which can be described as follows. In the first stage, multivariate linear mixed models can be used to model the longitudinal data. In the second stage, the time-to-event model is estimated by replacing the random effects with corresponding empirical Bayes estimates. There are three problems with directly applying this approach. First, estimation in the first stage is complicated by the fact that simply fitting multivariate linear mixed models results in bias due to informative dropout; this is demonstrated by Albert and Shih (2009) for the the case of P = 1. Second, as described in Section 2.1, parameter estimation for multivariate linear mixed models can be computationally difficult when the number of longitudinal measurements (P) is even moderately large. Third, calibration error in the empirical Bayes estimation needs to be accounted for in the time-to-event model. The proposed approach will deal with all three of these problems.
The bias from informative dropout is a result of differential follow-up whereby the longitudinal process is related to the length of follow-up. That is, in (2)–(3), patients with large values of are more likely to have an early event when αp > 0 for p = 1, 2, …, P. There would be no bias if all J follow-up measurements were observed on all patients. As proposed by Albert and Shih (2009) for univariate longitudinal data, we can avoid this bias by generating complete data from the conditional distribution of Xi = (X1i, X2i, …, XPi) given Ti, denoted as Xi|Ti. Since Xi|Ti under model (2)–(3) does not have a tractable form, we propose a simple approximation for this conditional distribution. Under model (2)–(3), the distribution of Xi|Ti can be expressed as
Since Ti and the values of Xi are conditional independent given γi, h (Xi|γi, Ti) = h (Xi|γi), where . The distribution of Xi|Ti can be expressed as a multivariate linear mixed model if we approximate g(γi|Ti) by a normal distribution. Under the assumption that g(γi|Ti) is normally distributed with mean μTi = μ01Ti, μ11Ti, μ02Ti, μ12Ti, …, μ0PTi, μ1PTi)′ and variance , and by rearranging mean structure parameters in the integrand of (5) so that the random effects have mean zero, Xi|Ti corresponds to the following multivariate linear mixed model:
where i = 1, 2, …, I, j = 1, 2, …, Ji, and p = 1, 2, …, P. The parameters and are intercept and slope parameters for the pth longitudinal measurement and for patients who have an event time at time Ti or who are censored at time tJ. In addition, the associated random effects are multivariate normal with mean 0 and variance , and the residuals are assumed to have an independent normal distribution with mean zero and variance . Thus, this conditional model involves estimating separate fixed effect intercept and slope parameters for each potential event-time and for subjects who are censored at time tJ. Likewise, separate random effects distributions are estimated for each of these discrete time points. For example, the intercept and slope fixed-effect parameters for the pth biomarker for those patients who have an event at time Ti = t3 is and , respectively. Further, the intercept and slope random effects for all P biomarkers on those patients who have an event at time Ti = t3, , is multivariate normal with mean 0 and variance . A similar approximation has been proposed by Albert and Shih (2009) for univariate longitudinal data (P = 1).
Recall that by generating complete data from (6) we are able to avoid the bias due to informative dropout. However, when P is large, direct estimation of model (6) is difficult since the number of elements in grows quadratically with P. For example, the dimension of the variance matrix is 2P by 2P for P longitudinal biomarkers. Fieuws and Verbeke (2005) proposed estimating the parameters of multivariate linear mixed models by formulating bivariate linear mixed models on all possible pairwise combinations of longitudinal measurements. In the simplest approach, they proposed fitting bivariate linear mixed models on all combinations of longitudinal biomarkers and averaging “overlapping” or duplicate parameter estimates. Thus, we estimate the parameters in the fully specified model (6) by fitting bivariate longitudinal models that only include pairs of longitudinal markers. Fieuws and Verbeke (2005) demonstrated with simulations that there is little efficiency loss using their approach relative to a full maximum-likelihood based approach. Fitting these bivariate models is computationally feasible since only four correlated random effects are contained in each model. (That is, is a four by four-dimensional matrix for each discrete event-time Ti.) Duplicate estimates of fixed effects and random-effect variances from all pairwise bivariate models are averaged to obtain final parameter estimates of the fully specified model (6). For example, when P = 4 there are (P − 1) = 3 estimates of for the pth longitudinal biomarker that need to be averaged.
Model (6) is then used to construct complete longitudinal pseudo data sets which in turn are used to estimate the mean of the posterior distribution of an individual’s random effects given the data. Specifically, multiple complete longitudinal data sets can be constructed by simulating Xpij values from the approximation to the distribution of Xi |Ti given by (6) where the parameters are replaced by their estimated values. Since the simulated data sets have complete follow-up on each individual, the bias in estimating the posterior mean of γi caused by informative dropout will be much reduced.
The posterior mean of distribution γi given the data can be estimated by fitting (2)–(3) to the generated complete longitudinal pseudo data. However, similar to fitting the conditional model (6), fitting model (2)–(3) is difficult due to the high dimension of Σγ. Thus, we again use the pairwise estimation approach of Fieuws and Verbeke (2005), whereby we estimate the parameters of (2)–(3) by fitting all pairwise bivariate models and averaging duplicate parameter estimates to obtain final parameter estimates. For each generated complete longitudinal pseudo data set, the estimate of the posterior mean, denoted as i = (1i0, 1i1, …, Pi0, Pi1)′, can be calculated as
where Zi is a PJ × 2P design matrix corresponding to the fixed and random effects in (2)–(3), where , and Vi is the variance of Xi. Estimates of , denoted as , are obtained by substituting (p0, p1, pi0, pi1) for (βp0, βp1, γpi0, γpi1) in (3).
To account for the measurement error in using i as compared with using γi in (1), we note that
In the second stage, α0j (j = 1, 2, …, J) and αp (p = 1, 2, …, P) can be estimated by maximizing the likelihood
where is given by (8). Thus, we propose the following algorithm for estimating α0j and αp (p = 1, 2, …, P):
We choose M = 10 in the simulations and data analysis since this was shown to be sufficiently large for univariate longitudinal modeling discussed in Albert and Shih (2009). Asymptotic standard errors of 0j and p cannot be used for inference since they fail to account for the missing data uncertainty in our procedure. Standard errors and 95% confidence intervals of parameter estimates using the bootstrap [Efron and Tibshirani (1993)] are as follows:
Covariates can be incorporated in (3) by adding them directly into the multivariate linear mixed model (6). Specifically, if , where Zi is a vector of covariates with ηp being parameters for the pth biomarker, then P (Xi|Ti, Zi) = ∫ h(Xi|γi, Zi)g(γi|Ti) dγi and Xi|Ti can be approximated by a multivariate linear mixed model with Zi ηp being added to the right side of (6). Estimation then proceeds as described in Section 2.2 Although more difficult, covariates can also be incorporated into (1). If
then P (Xi|Ti, Zi ) = ∫ h(Xi|γi, Zi )g(γi|Ti, Zi ) dγi, and under the assumption that g(γi|Ti, Zi) is normally distributed with variance not depending on Zi (which we found to be the case in simulations not shown), then Xi|Ti, Zi can be approximated by a multivariate linear mixed model with added to the right-hand side of (6). Extensive simulations showed that the conditional distribution of γi|Ti, Zi is nearly normally distributed over a wide range of parameter values, with slight departures from normality found when the relationship between the longitudinal process and event-time processes is very strong (i.e., αp’s are very large in magnitude). The multivariate linear model approximation is flexible in that it allows the regression parameters to vary with Ti. A more parsimonious model would be to constrain the parameters such that does not vary with Ti.
We conduct a simulation study to examine the statistical properties of the proposed approach. The approach is examined for the situation where P = 3, I = 300, and for p = 1, 2 and 3. The remaining parameters are presented in Table 1. We compare the proposed approach with M = 10 with a model where the ’s are assumed to be known (true model) and with a model where the observed Xpij’s are used in (1) (observed model). Although the true values are never actually observed in practice, we examine the true model as a benchmark in comparing the other models. Table 1 shows that the proposed approach results in nearly unbiased estimates of α1, α2 and α3, whereas the model which uses the observed observations (which are subject to measurement error) has severe bias for estimating α1 and α3 (the two parameters which are nonzero). Although the variability of parameter estimates is larger for the are presented in proposed approach as compared with the observed approach, the root mean squared errors are substantially smaller for α1 and α3 for the proposed approach. For example, the root mean squared errors for α1 is 0.089 for the proposed approach and 0.184 for the observed model. Results when the “true” values for the markers meters Table 1. (markers without measurement errors) are assumed known are also presented in Table 1 (column labeled Truth). As expected, estimates are unbiased for this gold standard case. A comparison of the gold standard case with the proposed approach shows efficiency loss. For example, the relative efficiency for estimating α1, α2 and α3 with the proposed approach versus the gold standard is 0.45 [(0.06/0.089)2], 0.57 and 0.24, respectively. We conducted additional simulations with different parameter values. In all cases tried, the mean squared errors were substantially smaller using the proposed approach as compared with using the observed values (data not shown).
Table 1 also presents the average estimated intercept and slope for each of the three longitudinal biomarkers. The results show that estimates of these fixed effects are nearly unbiased for the proposed approach.
We examine the effect of multiple longitudinal biomarkers on the short-term prognosis for patients with primary biliary cirrhosis (PBC) using the PBC study conducted at the Mayo Clinic from 1974 to 1984 [Murtaugh et al. (1994)]. PBC is a chronic disease characterized by inflammatory destruction of the small bile ducts within the liver, which eventually leads to cirrhosis of the liver, followed by death. Various biomarkers such as biliribin, prothrombin time and albumin were collected longitudinally, and interest is on examining whether these biomarkers relate to the natural history of disease. Of major interest was whether these biomarkers are prognostic for transplantation-free survival (time to either transplantation or death). A total of 312 patients had a baseline measurement and were followed longitudinally at 6 months and at yearly intervals thereafter.
For our application, we focused on the first 4 years of follow-up for a number of reasons. First, individual changes in the biomarkers appeared to be close to linear over this time period. Figure 1 shows 4 examples of complete follow-up in which the changes in log-transformed biliribin appear to be linear over the first 4 years of follow-up and not very linear over the whole range of follow-up. For each panel, the solid line is a least-squares regression line for the first four years of follow-up, while the dashed line is a corresponding line using all the follow-up data. The patterns over the whole follow-up period are nonlinear curves and are not systematic over subjects, and therefore not easily characterized by a simple nonlinear mixed model. The reason why the linear assumption is reasonable over the shorter time interval is that, even though the curves are nonlinear, they can adequately be approximated as linear functions over a short time interval (i.e., a nonlinear function can be locally approximated by a linear function). Second, the methodology makes the assumption that the effect of the biomarkers on prognosis is constant over the follow-up period [i.e., αp parameters in (1) do not vary over time]. This assumption is more reasonable over the shorter 4 year interval rather than the entire follow-up period.
Table 2 shows parameter estimates from fitting model (1) with the observed data as covariates instead of the true values. Although both standard errors and 95% confidence intervals were estimated using the bootstrap, only the 95% confidence intervals are presented since the bootstrap estimates were not normally distributed for many of the parameter estimates. The results demonstrate a statistically significant positive effect of biliribin and prothrombin time and a negative effect of albumin on transplantation-free survival. However, it should be recognized that these parameter estimates may be distorted due to the measurement error in these longitudinal biomarker measurements.
Using the proposed approach, we initially fit model (1)–(3) which incorporated a random intercept and slope term for each of the three biomarkers. However, the random effect for slope for prothrombin time and albumin were estimated as nearly zero. Thus, we re-fit the model without a random slope effect for these two biomarkers. Table 3 shows parameter estimates from the proposed approach with 95% confidence intervals estimated using the bootstrap (as in the analysis with the observed biomarkers presented in Table 2, we do not present the parameter estimates of the standard errors). Except for the effect of biliribin, estimates of the other two biomarkers are substantially larger in magnitude under the proposed approach than when ignoring measurement error and using the observed data (Table 2). This is consistent with the common phenomenon that ignoring measurement error attenuates parameter estimates. In terms of inference, the effect of prothrombin time on short-term prognosis is no longer statistically significant with the proposed approach, while the effects of biliribin and albumin on prognosis are statistically significant with both approaches. In the PBC analysis, estimates of were 0.31, 0.12 and 0.11 for log-transformed values of biliribin, albumin and prothrombin time, respectively. The smaller absolute values for parameter estimates of albumin and prothrombin time using the observed markers (Table 1) relative to estimates for these markers using the proposed approach (Table 2) can be attributed to attenuation due to measurement error, since, in these cases, the residual variances are substantially larger than the between-subject variations.
We also conducted analyses where we adjusted for treatment effect and age in the discrete-time survival model [Zi is treatment group or age in model (10)]. As discussed in Section 2.3, we constrained the parameters so that they did not vary with Ti (results were similar when we did not constrain the parameters). When we adjusted for treatment group (with treatment group coded as 1 for D-penicillamine and 0 for placebo) in (10), we estimated the α coefficient corresponding to treatment as 0.051 (95% CI: −0.84 to 1.14). The estimates of other parameters were almost identical to those presented in Table 3. When we adjusted for age in (10), we estimated the α coefficient corresponding to age as 0.018 (95% CI: 0.01 to 0.129), where age was scaled in units of a year. Although age was statistically significant, the effects of log bilirubin, albumin and prothtime on transplantation-free survival were similar to those for the unadjusted model. Specifically, the regression coefficients (α coefficients) corresponding to three markers are 0.394 (95% CI: 0.07 to 2.21), −10.65 (−83.51 to −5.52) and 5.76 (−4.28 to 30.97).
Table 3 also shows the estimated fixed effect intercept and slope for the three longitudinal biomarkers for the proposed approach. The estimates suggest that biliribin is increasing, while albumin and prothrombin time are nearly constant over time.
The joint modeling approach is important in this application for a number of reasons. First, survival models which use observed biomarkers can result in attenuated estimates of risk. The proposed approach allowed us to account for the measurement error in investigating the effect of multiple biomarker measurements on the short-term prognosis of PBC patients in terms of transplantation-free survival. With the proposed approach, we found that the “true” biliribin and albumin values at the beginning of an interval had a sizable and statistically significant effect on the probability of either a transplantation or death in the subsequent interval. Second, the proposed approach allows us to appropriately make inference about changes in the three “true” biomarkers over time. As stated before, the largest change over time was in biliribin which sizably increased over time. When making these longitudinal inferences, not appropriately modeling the relationship between the multiple biomarkers and survival may lead to bias due to informative dropout [Wu and Carroll (1988)].
We proposed an approach for jointly modeling multivariate longitudinal and discrete time-to-event data. Unlike likelihood-based approaches which require high-dimensional integration to evaluate the joint likelihood, this approach only requires fitting bivariate random effects models. This methodology uses recent methodology for fitting multivariate longitudinal data with bivariate linear mixed models proposed by Fieuws et al. (2005, 2007). They discussed the simple averaging of duplicate parameters estimates as we did in implementing the proposed approach. They also proposed a pseudo-likelihood approach which involves maximizing the sum of likelihoods from bivariate models across all combinations of pairwise longitudinal markers. Although this later approach may provide some minor efficiency gain over simple averaging, it would be substantially more complicated to implement in our setting. Further, one of the advantages of the pseudo-likelihood approach is that it provides an analytic expression for the asymptotic variance of the parameter estimates. Unfortunately, this asymptotic variance is not generalizable to the joint model with time-to-event data, making the pseudo-likelihood approach less attractive in our setting.
There are similarities between our approach and the recent approach by Fieuws et al. (2008) for predicting renal graft failure based on multivariate longitudinal profiles. Both approaches model the conditional distribution of the multivariate longitudinal profiles given the failure time. However, there are major differences between the two approaches. Fieuws et al. model the conditional distribution of the longitudinal measurements given the failure time and then use Bayes rule to estimate the probability of failure given the longitudinal profiles. In our approach, we use an approximation of the conditional distribution of Xi|Ti under the joint model of the multivariate longitudinal and time-to-event data in order estimate the parameters of this joint model.
We demonstrated the feasibility of the proposed approach with three biomarkers (P = 3). However, the approach can easily accommodate a large number of longitudinal profiles since it simply involves fitting bivariate models. The relationship between the multivariate longitudinal and event-time data is governed by expression (1). However, other functional relationships are possible with this approach. For example, we could relate the two processes by averages of “true” longitudinal biomarkers either across time or across different biomarkers. Alternatively, the approach could be formulated so that the event-time process depends on the individual’s intercept and slope for each of the P longitudinal biomarkers.
The multivariate longitudinal profiles are modeled as multivariate linear mixed models in (2) and (3), which was appropriate for the analysis of the PBC data. However, the methodology could be extended to allow for more flexible nonlinear modeling of marker profiles. This would involve approximating the conditional probability Xi|Ti in (5) where h(Xi|γi) follows a multivariate nonlinear mixed model, rather than the linear mixed model discussed in our paper. In the nonlinear case, we could approximate (5) by (6), where (6) would be a nonlinear mixed model with parameters indexed by Ti rather than the linear mixed model presented. However, unless there is biological rational for a particular nonlinear mixed model, it may be difficult to choose a reasonable model in most practical situations.
In our formulation, we assumed that event times are only administratively censored after a fixed follow-up at the end of the study. For the situation in which patients are censored prematurely, dropout times can be imputed based on a model fit using patients who had the potential to be followed over the entire study duration.
In this article methodology was developed using a discrete-time survival model with calibration error being incorporated by using a probit link function. This approach led to an analytically tractable form for incorporating calibration error (8). For the PBC study, little is lost by using a discrete-time model since the probability of an event in each of the five intervals is low (the estimated probability of an event during each of the five time intervals is 0.05, 0.03, 0.08, 0.07 and 0.07). For continuous-time survival models such as the Cox model, incorporating calibration is more difficult since there is no simple analytic solution. That said, we could use a Cox model if we do not account for the calibration error in replacing the random effect by their empirical Bayes estimators. Although not the case in the PBC study, in situations where the within-subject variation is small relative to the between-subject sources of variation, the calibration error will be small and there will be only a small amount of bias induced by not accounting for the calibration error.
We thank the Center for Information Technology, NIH, for providing access to the high-performance computational capabilities of the Biowulf cluster computer system. We thank the Editor, Associate Editor and reviewer for their constructive comments which led to an improved manuscript.
1Supported by the Intramural Research Program of the National Institutes of Health, Eunice Kennedy Shriver National Institute of Child Health and Human Development.