A Birnbaum-Saunders Model for Joint Survival and Longitudinal Analysis of Congestive Heart Failure Data

We consider a parametric joint modelling of longitudinal measurements and survival times, motivated by a study conducted at the Heart Institute (Incor), São Paulo, Brazil, with the objective of evaluating the impact of B-type Natriuretic Peptide (BNP) collected at different instants on the survival of patients with Congestive Heart Failure (CHF). We employ a linear mixed model for the longitudinal response and a Birnbaum-Saunders model for the survival times, allowing the inclusion of subjects without longitudinal observations. We derive maximum likelihood estimators of the joint model parameters and conduct a simulation study to compare the true survival probabilities with dynamic predictions obtained from the fit of the proposed joint model and to evaluate the performance of the method for estimating the model parameters. The proposed joint model is applied to the cohort of 1609 patients with CHF, of which 1080 have no BNP measurements. The parameter estimates and their standard errors obtained via: i) the traditional approach, where only individuals with at least one measurement of the longitudinal response are included and ii) the proposed approach, which includes survival information from all individuals, are compared with those obtained via marginal (longitudinal and survival) models.


Introduction
In many studies, repeated measurements of one or more variables (longitudinal responses), time until the occurrence of one or more events (survival responses) and additional observations on explanatory variables are collected on a set of subjects in order to characterize their relationship. This is the case of a study conducted at the Heart Institute (Incor), São Paulo, Brazil, where data related to: i) longitudinal measurements of B-type Natriuretic Peptide (BNP) levels (pg/mL), ii) the time between admission to the study and the date of death or censoring (in months), as well as iii) the values of basal covariates, were collected on patients with Congestive Heart Failure (CHF) to identify prognostic factors for the time to death.
In practical situations, data with this nature are often analyzed considering the longitudinal and survival responses separately as noted in Rizopoulos (2010), among others. However, there are two scenarios in which it is more appropriate to perform a joint modelling: when interest is to analyze the behavior of the longitudinal response, considering a possible dependence of time to dropout, potentially informative, treated as the survival response (Hogan & Laird 1997a, Hogan & Laird 1997b, Diggle, Farewell & Henderson 2007 and when interest is to analyze the time-to-event considering the effect of the longitudinal response measurements (Wulfsohn & Tsiatis 1997, Henderson, Diggle & Dobson 2000, Rizopoulos 2010, Crowther, Abrams & Lambert 2012. Different authors suggest that in these cases joint modelling can facilitate the understanding of the mechanisms underlying the phenomenon under investigation and can improve the properties of parameter estimators, being an appealing alternative that has attracted the attention of recent research (Tsiatis & Davidian 2004, Yu, Law, Taylor & Sandler 2004, Diggle, Sousa & Chetwynd 2008, Rizopoulos, Verbeke & Molenberghs 2010, Wu, Liu, Yi & Huang 2012, Rizopoulos 2012b, Gould, Boye, Crowther, Ibrahim, Quartey, Micallef & Bois 2015.
A naive approach when the interest lies exclusively in the survival component, is to consider the longitudinal response as a time-dependent covariate in the Cox model, which requires that the time-dependent covariate values be known exactly at each instant of failure and further, that the time-dependent covariates are external, as decribed by Kalbfleisch & Prentice (2002). This approach may not be appropriate because the longitudinal observations are usually measured intermittently and subject to errors. Additionally, it may be influenced by the occurrence of the event under investigation (Hu, Tsiatis & Davidian 1998, Greene & Cai 2004, Rizopoulos 2010). An alternative is the two-stage approach, where a model for longitudinal response is initially fitted to the data and, using the values of the first stage estimated individual longitudinal, are subsequently incorporated as a time-dependent covariate in the Cox model (Tsiatis, DeGruttola & Wulfsohn 1995, Ye, Lin & Taylor 2008, Albert & Shih 2010. Despite the advantage of its simple computational implementation, this method has the limitation of not considering the effect of the survival response on the modelling of the longitudinal data. Another alternative is to estimate the model parameters by maximizing the likelihood function corresponding to the joint distribution of longitudinal and survival responses (Wulfsohn & Tsiatis 1997, Henderson et al. 2000, Hsieh, Tseng & Wang 2006, Crowther et al. 2012, Rizopoulos 2010, Rizopoulos 2012b. Although the computational implementation of this approach is more complex, it has the advantage of using longitudinal and survival data simultaneously in the process of estimating model parameters. This approach is adopted in this paper. In this context, most authors consider the Cox model to describe survival times (see Wulfsohn & Tsiatis 1997, Henderson et al. 2000, Rizopoulos 2010, Rizopoulos 2012b, among others) although log-normal and Weibull parametric models have also been considered for such purposes in Schluchter (1992), Pawitan & Self (1993), DeGruttola & Tu (1994). Linear mixed models are commonly employed to represent the longitudinal component. The usual methods, however, only use data for subjects that have at least one measurement of the longitudinal response. In studies where such measurements are not recorded for some participants, the corresponding estimates may be biased or less efficient. Among many alternatives to describe survival times, Birnbaum-Saunders distributions seem appropriate in the context of CHF because in chronic cardiac diseases, a cumulative damage caused by several risk factors may lead to a degradation and to a consequent failure, an inherent feature of such models, as described in Galea, Leiva & Paula (2004), Leiva, Barros, Paula & Galea (2007), Barros, Paula & Leiva (2008), Balakrishnan, Leiva, Sanhueza & Vilca (2009) or Leiva, Athayde, Azevedo & Marchant (2011). With this in mind, we propose a joint modelling methodology that considers a linear mixed model to describe the longitudinal response and a Birnbaum-Saunders model to describe the survival response, allowing the inclusion of the survival data of subjects without longitudinal observations.
In Section 2, we present the model and discuss inferential aspects including the dynamic prediction of the survival probability based on the available data up to the instant when we want to make the prediction. In Section 3, we summarize the results of a simulation study designed to compare the true survival probabilities with dynamic predictions obtained from the fitted model and to evaluate the performance of the method for estimating the model parameters. In Section 4 we analyze the data that motivated our research. We conclude with a discussion and suggestions for future research in Section 5.

Proposed Methodology
Consider a set of n subjects followed over the time interval [0, τ ), τ > 0, and suppose that for the i-th subject (i = 1, . . . , n) we observe: i) a sequence of measurements of a longitudinal response y ij = {y i (t ij ), j = 1, . . . , n i } summarized in y i = (y i1 , . . . , y ini ) ⊤ and occurring at times t ij represented in t i = (t i1 , . . . , t ini ) ⊤ , ii) a record of the time between admission to the study and the occurrence of the event of interest (T i ) or censoring (C i the first a h of which are independent of time. The subscript h = 1, 2 indicates whether they correspond to the longitudinal or to the survival components, respectively.
The longitudinal response for the i-th subject at time t ≥ 0 is modelled as where m i (t) = µ 1i (t) + w 1i (t) denotes the true value of the longitudinal response, specified as a function of a mean response µ 1i (t) = x 1i (t) ⊤ β 1 , with β 1 representing the fixed effects corresponding to p 1 explanatory variables in x 1i (t) and the process w 1i (t), characterized in terms of a specific time invariant random intercept for the i-th subject, b 0i ∼ N (0, σ 2 0 ) and e i (t) ∼ N (0, σ 2 e ) denotes the measurement error, considered independent of b 0i , for all t ≥ 0.
We assume that the survival or censoring time observed for the i-th subject follows the log-linear Birnbaum-Saunders regression model where β 2 contains the fixed effects corresponding to p 2 explanatory variables in x 2i and the model errors ε i ∼ SinhN(α, 0, 2), with SinhN(α, ψ, σ) denoting the Normal hyperbolic sine distribution with α, ψ and σ representing the shape, location and scale parameters, respectively. The associated density and survival functions are, respectively and For details on the relation between the Birnbaum-Saunders and the SinhN distributions, see Leiva et al. (2007).
In this set-up we develop the likelihood function corresponding to the joint distribution of longitudinal and survival responses. Assuming that the random effects b 0i (i = 1, . . . , n) take into account both the association between the longitudinal and survival responses, and the correlation between the longitudinal observations (conditional independence), that the censoring mechanism and the observation times of the longitudinal response are not informative and considering independence between survival and censoring times [see Rizopoulos (2012b)], the joint likelihood function is where ω i = 1 if the i-th subject has at least one observation of the longitudinal response and 2 ) ⊤ denoting the vectors containing the parameters for the survival responses for subsets of subjects that have or have not measurements of the longitudinal outcome, respectively, θ y = (β ⊤ 1 , σ 2 e ) ⊤ denotes the vector containing the longitudinal response parameters and θ b = σ 2 0 . Additionally, letting M i (t) = {m i (u), 0 ≤ u ≤ t} denote the history of the true unobserved longitudinal process up to time t, we assume that where f i (·) and S i (·) respectively denote the probability density and survival functions of the SinhN distribution with shape parameter α > 0, scale parameter σ = 2 and location parameter representing Normal probability density functions for the longitudinal response and the random effects, respectively. The parameter γ measures the association between the longitudinal and survival processes. Furthermore, we assume that where f i (·) and S i (·) respectively denote the probability density and survival functions of the SinhN distribution with shape parameter α > 0, scale parameter σ = 2 and location parameter ψ 0i = x ⊤ 2i β 2 .
Explicit expressions for the terms composing the likelihood function (5) are given in the Appendix.
Maximum Likelihood (ML) estimates of the joint model parameters are obtained by direct maximization of (5) using a quasi-Newton algorithm implemented via the PORT routines in the R optimizer nlminb (see Gay 1990), one of the algorithms employed in the JM library (Rizopoulos 2012a). Numerical integration is required because the integrals with respect to the random effects b 0i , i = 1, . . . , n in (5) have no analytical solution. For such purposes, we use Gauss-Hermite quadrature methods as suggested by Wulfsohn & Tsiatis (1997), Henderson et al. (2000) or Rizopoulos (2010), among others, for situations where the random effects vector for each subject has low dimension. In particular, this is the default method employed in the JM library (Rizopoulos 2012a), except in cases where a large number of random effects per subject is available. In such cases, Rizopoulos, Verbeke & Lesaffre (2009) . For such purpose we consider the structure of the estimator of π i (u | t) using the empirical Bayes estimate for b i proposed by Rizopoulos (2011), namely where θ denotes the maximum likelihood estimates of the model parameters θ and with n i (t) denoting the number of longitudinal measurements of the i-th unit up to time t. The performance of (9) for finite samples depends on the quality of the ML estimates of θ and on the prediction of the random effects b i .
The methodology proposed in this paper was fully implemented in R (R Development Core Team 2013).

Simulation
We conducted an extensive simulation study to compare the true survival probabilities with dynamic predictions obtained from the fitted joint model, based on the longitudinal data collected up to different instants and to evaluate the performance of the method for estimating the model parameters.
The longitudinal response for the i-th subject at time t ≥ 0, y i (t), was generated by (1) under Normal distributions for both the random effects and error terms considering observation time and dichotomized CHF etiology (Chagas disease = 0 or Other cardiomyopathies = 1) as covariates. The survival response was generated by a log-linear Birnbaum-Saunders regression model (2) including only CHF etiology as a covariate. To predict the survival probabilities, we considered only individuals with at least one measurement of the longitudinal response, i.e., those for which ω i = 1 in (5).
The location parameter for the Birnbaum-Saunders model can be expressed as and the variance components are σ 2 0 and σ 2 e . The parameter values were taken as the estimates obtained by fitting the joint models to the 529 patients with at least one longitudinal observation in the Incor data.
First we generated the longitudinal observations for each patient and then considered the corresponding mean observation as well as the CHF etiology covariate to generate the survival data, inducing an association between the two types of response. To mimic the set-up in the Incor data, the survival times were right censored either by specifying a Type I censoring scheme with a maximum followup time of τ = 180 or by randomly selecting censoring times from a Uniform distribution in the [0, τ ] interval.
Sixteen different scenarios resulting from the combination of 4 sample sizes (n = 100, 250, 500, 1000) and 4 censoring percentages (p c = 0%, 25%, 50%, 75%) were considered. For each scenario we generated 500 samples and, for each sample, we randomly selected 95% of the subjects to fit the Birnbaum-Saunders timeto-event model; the remaining 5% were considered to estimate the conditional survival probabilities based on the estimator π i (t + ∆t | t) in (9), as a function of the longitudinal response observation times t = 0, 24, 48, 72, 96 and the time increments ∆t = 12, 24, 36. The true probabilities π i (t + ∆t | t) were computed using the true parameter and the true (generated) random effects values. The comparison was carried out in terms of the absolute differences between estimated conditional survival probabilities and their true values. The performance of estimators was evaluated via the bias ( θ − θ) with θ = 500 −1 ∑ 500 l=1 θ l , the relative bias [( θ − θ)/θ] as well as via the square root of the mean squared error (MSE), [500 −1 ∑ 500 l=1 ( θ l − θ) 2 ] 1/2 . For each scenario, we constructed histograms for such differences based on the 500 samples and computed the corresponding 95-th percentile for each combination of t and ∆t.
We show the results obtained for a sample of size n = 1000. Results for the remaining scenarios may be obtained in https://www.ime.usp.br/∼acarlos/Diana/ ApendiceF.pdf In Table 1 we present these values for the case without random censoring and 7% Type I censoring. The number of differences corresponds to the available longitudinal measurements for the 50 subjects selected for prediction in the 500 samples. The results indicate that for a fixed time t, the differences between the true and estimated probabilities increase as ∆t increases. This is justified by the increasing distance between the time up to which longitudinal data are available and the instant for which the prediction is made. Furthermore, for ∆t fixed, there is a decrease in these differences as t increases. A possible explanation is that availability of more longitudinal measurements for each subject improves the predictions of the random effects involved in the proposed estimator. To verify this, we computed 95-th percentile of the empirical distribution (based on the 500 generated samples) of the differences between the true and predicted random effects. The results for the scenario described above are displayed in Table 2 and confirm our conjecture. We also computed the mean and median of the maximum (among patients) differences between the true and estimated survival probabilities for the 500 generated samples. The results, displayed in Table 3 also confirm our previous conclusions. In all cases, the bias, the relative bias and the square root of the MSE of the parameter estimates decrease as the sample size increases or the percentage of censoring decreases.
In Table 4 we present the average mean and median of the maximum (among patients) differences between the true and estimated survival probabilities for the 500 generated samples for all combinations of t and ∆t values, weighted by the number of maxima in each case. In Table 5, we present the bias, relative bias and square root of the MSE for the case where n = 1000, no random censoring and 7% Type I censoring time. Comparing the results from the different scenarios, we observed that the longitudinal component parameter estimates are more stable than those associated to the survival component which have a poor performance for small (n = 100) sample sizes. In summary, we conclude that the performance of the proposed model is better when: i) the sample size is larger and ii) the right censoring is smaller. In particular when the right censoring percentage is large, the parameter estimates have a poor performance and the prediction of the random effects are not as good as those obtained when more longitudinal measurements are available. This leads to a decrease in the quality of the prediction of the dynamic survival probabilities. Therefore the Birnbaum-Saunders joint model is recommended for situations where the underlying rationale for such model is reasonable, the sample size is large and there is little right censoring.

Analysis of the Incor Data
The proposed joint model was applied to the cohort of 1609 patients with CHF, of which 1080 have no BNP measurements.
Initially, the adequacy of the Birnbaum-Saunders model for the survival times description in comparison with the exponential, Weibull and log normal models was verified, considering the analysis of the Cox-Snell residuals (see Klein & Moeschberger 2003).
Then, we conducted an analysis involved a selection of "acceptable" models for the longitudinal and for the survival data fitted separately via standard techniques as suggested by Wu et al. (2012). The 529 patients with at least one longitudinal observation were considered for the former and the complete set of 1609 patients was used for the latter. In this process, each of 24 covariates were fitted individually along with CHF etiology (as suggested by the physicians) and the significant ones were subsequently fitted simultaneously to either the longitudinal or the survival component. The separate longitudinal and survival component models were sequentially refitted with the removal of the non-significant variables or grouping levels with non-significant effects in each step. The (significant) variables in the last step were chosen to compose the joint model. Observation time, CHF etiology and left atrium diameter were used as covariates for the longitudinal component of the model whereas CHF etiology and left ventricular ejection fraction were considered as covariates for the survival response.
The longitudinal component was modelled via (1) with where x 1i1 (t) denotes the time at which the response was observed, x 1i2 , x 1i3 , x 1i4 , x 1i5 , x 1i6 are dummy variables associated to the categories of CHF etiology (dilated, ischaemic, hypertensive, alcoholic or other cardiopathies), x 1i7 , x 1i8 are dummy variables associated to the categories of left atrium diameter (augmented or missing) and w 1i (t) consisting of a random intercept b 0i ∼ N (0, σ 2 0 ). The survival component was modelled via (2) with where x 2i1 , x 2i2 , x 2i3 , x 2i4 , x 2i5 are dummy variables associated to the categories of CHF etiology, x 2i6 , x 2i7 , x 2i8 are dummy variables associated to the categories of left ventricular ejection fraction (very low, low or missing) and ε i ∼ SinhN(α, 0, 2).
Finally, the association between both components was imposed by (7).
Maximum likelihood parameter estimates and their standard errors obtained via: i) the proposed joint model approach, which accommodates survival information of all individuals [likelihood given by (5)] and ii) the traditional joint model approach, where only individuals with at least one measurement of the longitudinal response are included (likelihood given by the first component of (5), i.e., when ω i = 1), were compared with those obtained via marginal (longitudinal and survival) models in each set-up. The results are summarized in Tables 6-9.
In case i), although standard errors are slightly smaller under the joint model, no relevant differences between the joint and marginal models longitudinal parameter estimates were detected with the exception of the time coefficient for which non-significance is more evident under the former model (p = 0.1100 versus p = 0.0586). Estimates of the survival parameters and corresponding standard errors are comparable for both models and the association parameter estimate is positive and marginally significant (p = 0.0609).
In case ii), the standard errors of the longitudinal parameter estimates are consistently smaller under the joint model, enhancing the significance of the time coefficient (p < 0.0001 versus p = 0.0586). Survival parameter estimates and corresponding standard errors obtained under the marginal and joint models are relatively different, leading to changes in the significance in some cases. The association parameter estimate is negative and highly significant (p < 0.0001).
To evaluate the effect of the percentage of patients with no longitudinal observations we refitted the joint model to the data of the 529 patients with at least one longitudinal measurements eliminating at random these measurements for 10%, 30% and 50% of them. The results obtained, displayed in Table 10, suggest that the longitudinal components are less affected than the survival ones, when the percentage of subjects without measurements of the longitudinal response increases. Although there are some changes in the sign of some coefficients, this is probably due to the lack of significance in most cases. In particular, the estimates of the coefficient of association between the longitudinal and survival components of the model are affected as the proportion of patients without at least one measurement of the longitudinal response increases, reaching a positive value for the case where 50% of subjects with no longitudinal measurements are included. Furthermore, we note that the standard errors tend to be greater as the proportion of patients without at least one measurement of the longitudinal response increases.

Discussion
We propose a methodology for joint modelling of longitudinal and survival data, which differs from the methods proposed in the literature by considering a Birnbaum-Saunders model to describe the survival response and incorporating the survival information of subjects without observations of the longitudinal response.
The results of the simulation and practical application to the Incor data when only individuals with at least one measurement of the longitudinal response are included in the joint model, suggest that the inclusion of longitudinal measurements of an appropriate response may be employed to improve the analysis of survival data. In particular, an increase in the number of subjects with measurements of the longitudinal response can improve the evidence of the association between the longitudinal and survival responses and can lead to an increase in the precision of parameter estimates. Further, an increase in the number of observations of the longitudinal response collected in a subject can improve the quality of the prediction of survival. Despite the need of joint modelling in this case, according to the results obtained previously, the results were not so evident in the practical application to the Incor data when survival information of all subjects is considered in the joint model, probably because of the observational nature of the study, carried out during 9 years with no fixed protocol for data collection. The large proportion (67%) of patients with no measurements of the longitudinal response may have masked the association between the two components of the joint model. Despite these results, the proposed methodology can be useful if there are cases of subjects without longitudinal information that should be included in the analysis.
In addition, the results obtained show the usefulness of the Birnbaum-Saunders model to describe the survival times in situations in which a cumulative damage caused by several risk factors may lead to a degradation and to a consequent failure, as in the case of CHF.
Future research is needed before this approach can be routinely used in practical problems. In particular, we mention diagnostic techniques and simulation studies to determine the effect of the proportion of units without longitudinal data on the joint model parameter estimates, considering the context of balanced studies where all subjects have the same amount of measurements of the longitudinal response, because the performance of the estimation method depends not only on the number of units with measurements of the longitudinal response, but also on the number of longitudinal observations collected in each subject, as previously discussed. ( . The Gauss-Hermite quadrature approximation of the logarithm of (5), namely, is given by , with s k denoting the k-th root of the Q-th order Hermite polynomial and w k the corresponding weight, } .