Bayesian QuantileBased Joint Modelling of Repeated Measurement and TimetoEvent data, with an Application to Lung Function Decline and Time to Infection in Patients with Cystic Fibrosis
{onecolabstract}Background The most widely used approach to joint modelling of repeated measurement and time to event data is to combine a linear Gaussian random effects model for the repeated measurements with a logGaussian frailty model for the timetoevent outcome, linking the two through some form of correlation structure between the random effects and the logfrailty. In this approach, covariates are assumed to affect the mean response profile of the repeated measurement data.
Objectives Some applications raise substantive questions that cannot be captured by this structure. For example, an important question in cystic fibrosis (CF) research is to understand the impact of a patient’s lung function trajectory on their risk of acquiring a variety of infections, and how this varies at different quantiles of the lung function distribution.
Methods Motivated by this question, we develop a joint quantile modelling framework in this paper with an associated Markov Chain Monte Carlo algorithm for Bayesian inference.
Results The translation from the common joint model towards quantile regression succeeds and is applied to CF data from the United Kingdom. The method helps detecting an overall difference in the relation between lung function decline and onset of infection in the different quantiles.
Conclusions Joint modelling without taking into account the special heteroscedastic structure is not sufficient in certain research question and the extensions towards models beyond the mean is necessary.
1 Introduction
Amongst the many statistical methods that have been developed for anlysing data from longitudinal studies, the term joint modelling refers to the statistical modelling of data in which each subject provides data on two qualitatively different kinds of outcome variable: a timesequence of repeated measurements; and a (possibly rightcensored) timetoevent variable. The extensive literature on this topic is reviewed in Tsiatis and Davidian (2004), and in a recent text by Rizopoulos (2012) and Özgür et al. (2015). Most of the joint modelling literature adopts a hierarchical modelling approach in which the repeated measurement and timetoevent outcomes are modelled as conditionally independent linear Gaussian and loglinear proportional hazards models given a latent bivariate stochastic Gaussian process, say whose components affect the conditional mean of the measurement process and conditional loghazard of the timetoevent process, respectively.
In practice, this approach is usually associated with an inferential focus on how a subject’s repeated measurement process affects their prognosis for survival, after adjustment for covariate effects. Quantile regression methods (Koenker and Bassett, 1978), in contrast, are designed to answer questions concerning the relationship between a subject’s covariates and the corresponding quantile of their measured outcomes. Papers that describe quantile regression methods for repeated measurement outcomes include Koenker (2004) and Geraci (2014) in a likelihood based scenario, Fenske et al. (2011) in the boosting context, and Yue and Rue (2011) using a Bayesian approach. To our knowledge, the one paper that extends this to joint modelling as defined above is Farcomeni and Viviani (2015).
In this paper, we develop a novel approach to quantilebased joint modelling, motivated by a question in the epidemiology of cystic fibrosis, a genetic condition that leads to a progressive deterioration of a patient’s lung function throughout their life.
In Section 2 of the paper we give a description of cystic fibrosis, the specific research question that motivated this work and the data that we will use to answer the question. In Section 3 we set out our proposed methodology. We formulate a location scale mixture representation of the asymmetric Laplace distribution whose likelihood can be treated like that of a latent Gaussian distribution, thus allowing the use of tools previously developed for joint modelling based on mean regression. We also develop an MCMC algorithm for Bayesian inference. In Section 4 we describe our analysis of the cystic fibrosis data. Section 5 discusses some limitations of the current methodology and outlines further work to extend its scope. Code is available from the first author.
2 Outline of the problem
Cystic fibrosis (CF) is one of the most common serious genetic diseases in the Western world. It has an impact on a number of organs, primarily the lung, pancreas and liver. The disease is characterized by recurrent lung infection and inflammation with associated longterm lung function decline. Most people with cystic fibrosis die prematurely as a result of respiratory failure.
Most epidemiological studies of lung function in CF have concentrated on investigating how a range of risk factors affect the decline in mean lung function in children and adult populations (Konstan et al., 2007; Salvatore et al., 2011; Konstan et al., 2012; Salvatore et al., 2012; TaylorRobinson et al., 2012, 2013). A number of studies have found that infection with Pseudomonas aeruginosa (PA) has an accelerating effect on the decline in mean lung function (Rosenfeld et al., 2012; TaylorRobinson et al., 2012). However, this acceleration is potentially greater and thus more concerning for patients whose lung function is already relatively poor, i.e. for individuals whose lung function is at lower quantiles of the distribution. A first step to further explore the impact of PA on lung function is to set the problem in a quantile regression context (Koenker and Bassett, 1978). Instead of assessing the impact of a set of covariates on the mean, we analyze their impact on selected quantiles. This provides additional insights without the need to assume a closed form distribution; we give a more detailed description in Section 3.
Quantile regression has been used rarely in research on pulmonary diseases. Examples include van Sickle et al. (2011), who showed differences in lung function decline for groups with different socioeconomic background, Denaro and Bailey (2012), who modelled lung function decline in CF, and Kulich et al. (2005), who estimated reference equations to be able to classify the lung function level of CF patients in comparison to other CF patients, rather than to a healthy population.
In this paper we use data from the UK CF registry, which collects longitudinal data on all patients with CF living in the UK. The registry currently collects data on around 10,000 patients who attend annual examinations examinations where they are assessed for clinical status, pulmonary function and microbiology of lower respiratory tract secretions. Lung function is measured using forced expiratory volume in 1 second as a percentage of predicted (%FEV) which is an agestandardised measure of lung function. We analyse an extract from the registry, containing a total of 16,872 %FEV measures from the 3,200 patients who were seen at least twice between 1995 and 2009. Of the 3,200 patients, 1,237 (39%) were infected with PA at some point during the followup period. We aim to develop a better understanding of the connection between the different levels of decline in the lung function and the risk of PA acquisition.
A quantile regression analysis of the repeated measurements of lung function reveals clear differences in the impact of PA at different quantiles of %FEV. We fitted the following model,
(1) 
where, for each patient, is the vector of the age at which %FEV was measured, contains zeros until the onset of infection and years of infection thereafter, and are random intercept and slope terms, and is the quantile of %FEV. For inference on this model, we used the Bayesian software tool BayesX (Belitz et al. (2013)). Figure 1 shows that the impact of PA on lung function decline, as measured by the parameter , is less pronounced at central and higher quantiles of the distribution. This suggests that infection with PA has a bigger impact on lung function for patients with worse preexisting lung function and constitutes a prima facie case for investigating structural differences in the relationship between lung function and PA infection at different quantiles of %FEV.
Joint modelling of the kind described in Section 1 has previously been used to investigate cystic fibrosis data for the relationship between lung function decline and survival of CF patients (Barrett et al., 2015; Schluchter et al., 2002). To best of our knowledge, the present paper is the first attempt to use quantile regression methods to understand the relationship between lung function decline and risk of PA infection.
3 Quantilebased joint modelling
3.1 Mean Regression Bayesian Joint Modelling
We first set out an example of a mean regression joint model proposed in Faucett and Thomas (1996), which provides the starting point for our qunatilebased methodology and our associated Monte Carlo Markov Chain (MCMC) algorithm. The model specification is
(2)  
where is the th repeated measurement on the th of individuals, is the corresponding measurement time and the are specified functions of , covariates and individuallevel random effects that have an impact on the repeated measurements (subscript ), the hazard function (subscript ) or both (subscript ).
To estimate the distribution of the survival times we use the partial likelihood for the hazard rate. Since the composition of the predictor is rather messy we present the parts in a structured form.

where the in the sum can be different kinds of functions of the covariates and the time . The parameter vector includes all parameters necessary for the corresponding functions. The index indicates that this part of the predictor is assumed to be influential for both parts of the model, the longitudinal outcomes and the survival times. In the survival part of the model this part of the predictor is multiplied by the association parameter .

. The structure of this function is the same as above, but it contains the covariates being assumed to only be relevant for the longitudinal part of the model.

. The structure of this function is the same as above, but it contains the covariates being assumed to only be relevant for the survival part of the model.
To achieve better legibility we present the functions without the corresponding arguments (i.e. instead of ), unless the arguments differ from the and . This is the case e.g. in the second line of the following formula, where the function is evaluated at the survival times or the third line, where the argument of the formula is the integration variable. The likelihood for the model of structure (2) results in:
where is the number of individuals, the number of observations of the th individual, is the time individual enters the study and is either the time of event or last time observation. The dependent variable is the th observation of individual . To cope with the baseline hazard , the time is split into a grid for which is approximated by a piecewise constant function with values for . The integral in equation (3.1) can hence be approximated by a sum of integrals. The association parameter links the mixed model to the survival analysis and quantifies the strength of the influence. In the following we will give a short overview over the full conditionals and refer to the appendix for more detailed calculation. We present the full conditionals exemplary for linear effects (e.g. ), for possible extensions see section 5.
Predictor part only concerning longitudinal data:
The predictor is independent of the survival part of the model. When choosing a Gaussian prior for the parameters the resulting full conditional is also of Gaussian type. Just as in standard Bayesian longitudinal modelling, inference on random effects can be implemented by inducing the correlation within the individual/group by using an appropriate prior covariance matrix.
Predictor part only concerning survival data:
The approach is reverse to the above presented. The likelihood for the parameters collected in can be presented independently of the longitudinal process. The full conditional can hence be constructed as if we were only doing a Cox regression. Unfortunately there is no closed form distribution for the full conditional and we follow Faucett and Thomas (1996) doing adaptive rejection sampling (ARS).
Predictor part concerning both longitudinal and survival data:
The full conditional for the predictor part is even more complicated than the survival specific part, since it appears in both parts of the likelihood. We hence follow Faucett and Thomas (1996) again and use an ARS step.
Model variance:
The model variance is not linked to the survival part of the model, the full conditional is hence conform with the result in Bayesian mixed models.
Baseline Hazard:
The relevant part of the likelihood (3.1) for the full conditional of the is the second part, which reduces to a gamma distribution.
Association Parameter:
The full conditional for consists of the second term of the likelihood and is also generated by ARS.
3.2 Bayesian Quantile Regression
Our proposed adaptation of the above framework to quantilebased joint modelling is as follows:
where everything is just as in model 2, except for the error not being distributed with a Gaussian distribution and the outcome not being the mean, but a quantile. The function denotes the cumulative distribution function of . We hence have one model for each quantile of interest . Before we explain, how this is included into the joint modelling framework, we sketch out the basic idea of inference in quantile regression. Quantile regression can be estimated by minimising the so called check function:
which puts weights on the distance of the observations to the estimated line with the value . Optimising this criterion however is not straightforward, given the non differentiability due to the absolute value characteristics. The Bayesian alternative most commonly used, is to use the asymmetric Laplace distribution (ALD) as an auxiliary likelihood distribution:
which when maximising leads to the same point estimators as minimising the checkfunction. This detour obviously inherits the problems from the checkfunction as it still contains the non differentiable part of the latter. Kozumi and Kobayashi (2011) suggested to use the location scale mixture presentation of the ALD to circumvent those problems. To this end define weights that follow a exponential distribution with rate : and a standard Gaussian random variable . Define the auxiliary variables and , then Y constructed in the following way:
follows an . This leads to the possibility to construct a regression model of the form:
which then renders possible to conduct MCMC in the same way as done for mean regression with an extra sampling step for the above defined weights. Here a summary of the standard approach to longitudinal modelling:
Fixed and Random effects
The priors for the model parameters are chosen to be Gaussian with the appropriate covariance matrices. The resulting full conditional is also Gaussian, with slight changes in the parameters compared to mean regression. Hyperpriors (e.g. to perform selection algorithm) do usually not have a different impact than in mean regression.
Model Variance
If choosing an inverse gamma distribution as a prior for the model variance, the full conditional is an inverse gamma distribution, too, again, with slight parameter changes, caused by offset and weights.
Weights
The key difference to the mean regression MCMC approach is the extra sampling step, necessary for the weights. Resulting from the location scale mixture representation of the ALD, the weights are a priori exponentially distributed with rate . The full conditional is an inverse Gaussian distribution.
Given the limited number of changes, having to consider in the transition from mean to quantile regression, we considered the transition from Gaussian joint modelling to Quantile regression joint modelling.
3.3 Bayesian Quantile Regression in Joint Modelling
The above described data set shows strong differences in the covariates when measuring the impact on different quantiles. The structure of the presented model will hence be extended to a model describing the association between a mixed model quantile regression and the survival model. The ultimate goal of this experiment is to model the survival function based on a set of quantiles mimicking the whole distribution and thus obtaining a more complete information on the relation between conditional distribution and survival time.
To achieve this goal we use the above described location scale mixture construction of the ALD. This makes the necessary adjustments in the likelihood (3.1) where the first line has to be replaced by
4 Analysing the cystic fibrosis registry data
4.1 Model formulation
As described in Section 2, our repeated measurement outcome is %FEV and our timetoeventoutcome is the time of first infection with PA. The longitudinal predictor was chosen to contain an intercept and a linear time effect as well as random intercept and time slope. The shared predictor was chosen to include the random effects. The list of predictors is:
This is a simple model formulation and by adding the linear fixed time effect to the longitudinal submodel, but not to the survival part, we make sure that the differences in the temporal progression captured but the magnitude of the quantiles do not interfere with .
The model was implemented in R Development Core Team (2011), MCMC chains consisted of a sample of with a burnin period of and a thining of nine. The sampling paths showed no inconsistencies.
4.2 Results
We ran the joint Bayesian quantile model for the nine percentiles . Figure 2 displays the different in the association parameter for the nine different quantiles. The association parameter was considered significant, if 95% of the values in the sample are on one side of and the corresponding boxes colored grey in the figure.
To aid interpretation of these values Figure 3 shows a selection of the longitudinal submodel for four different quantiles. The black lines display the overall quantile specific regression estimation, while the grey lines show the deviation from this overall estimation for individual patients. We can see that the individual lines scatter rather differently below and above the black lines. The impact of this deviation on the hazard rate is quantified by . Negative values indicate that the patients which deviate negatively from the overall trend result in having a multiplicative constant bigger than one and hence have and increased risk of infection. If we look at the example of the negative (see Figure 2) implies that risk of infection is higher for patients, whose lines lie in the group below the dark line in the upper left figure in 3. Positive have the opposite effect: for the more central quantiles, we can see that the values indicate that the individuals above the line have a higher risk of infection whereas the risk is lower for individuals below the line. Please note for the direct interpretation that the individuals are denser around the mean trajectory for the central quantiles and the fact that the values are higher does not mean that the influence is higher. Note that the individuals are not necessarily in the same group for the different quantiles. The overall conclusion is that there are differences in the relation between lung function decline and onset of infection in the different quantiles. Note that the purpose of this paper is to outline the utility of our statistical approach. Future analyses of CF data using the full dataset and taking into account a complete range of covariates will clarify the relationship between lung function and onset of infection with PA.
5 Summary and Outlook
We have shown how joint models for repeated measurement and timetoevent outcomes can be adapted to include ideas from quantile regression. The proposed model could be extended in various ways. One such extension would be to use to use lagged effects between the repeated measurement and timetoevent components of the model. A second extension would be a model in which the hazard for the timetoevent outcome depends on the rate of change in the repeated measurement process. For a random intercept and slope model, this is straightforward, we simply include the random slope as shared random effect. More generally, the model would need to include a differentiable stochastic process, say, amongst the random effects in the repeated measurement component, with the derivative of included in the timetoevent component. Estimation of effects of this kind is difficult without relatively long, frequently observed series of repeated measurements on individual patients. A third extension, also straightforward in principle but difficult in practice without extensive followup data, would be to allow nonlinear timevarying effects.
Including more elaborate effects should be straight forward, since the location scale mixture allows for extending the predictor by numerous effects (Waldmann et al., 2013). Besides nonlinear effect also spatial information could be included. Modelling the random effects with a Dirichlet process mixture and then including the resulting cluster information into the survival model could also lead to interesting insights into the structure of infection.
A further interesting extension would be to replace the quantile model by a Bayesian distributional regression as done by Klein et al. (2015) or with any other kind of non mean regression in order to get a deeper understanding of the connection between a longitudinal process and time to event observations.
Acknowledgments: DCTR was supported by an MRC Population Health Scientist Fellowship (G0802448) for this work. We thank the UK Cystic Fibrosis Trust for access to the UK cystic fibrosis registry and all the centre directors for the input of data to the registry. EW was supported by the Interdisciplinary Center for Clinical Research (IZKF) of the FriedrichAlexanderUniversity ErlangenNürnberg (Project J61). Special thanks go to Peter J Diggle whose comments and insights in the topic that greatly improved the manuscript.
Appendix A Full Conditionals
a.1 Bayesian Quantile Regression
This is the likelihood setup resulting from location scale mixture explained in section 3.2:
We demonstrate the calculation for the example of a linear effect , which forms part of the assumed additive predictor . We assign a Gaussian prior define the sub predictor . Then we get a Gaussian full conditional for parameter vector . This result is generasible towards all kinds of effects with Gaussian prior (e.g. PSplines, Gaussian Markov random fields or random effects).
We assign an inverse gamma distribution as prior to the variance parameter : .
The full conditionals for the weights are calculated separately for each :
With the density transformation theorem follows:
a.2 Bayesian Joint Mean Regression
This is the likelihoodsetup from section 3.1:
To show the full conditional for effects on the longitudinal sub predictor we again demonstrate the calculation for the example of a linear effect , which forms part of the assumed additive predictor . We assign a Gaussian prior define the sub predictor . Then we get a Gaussian full conditional for parameter vector . This result is generasible towards all kinds of effects with Gaussian prior (e.g. PSplines, Gaussian Markov random fields or random effects). Please note that is merely an example and indices are thus suppressed.
We assign a gamma distribution as prior to the pieces of the piecewise constant baseline hazard : .
We assign an inverse gamma distribution as prior to the variance parameter : .
a.3 Bayesian Joint Quantile Regression
Everything is analogous to the above. The results gets either obvious, when extending what is described in A.1 by the extra subpredictor or when replacing the Gaussian distribution in A.2 by the location scale mixture.
References
 Barrett et al. [2015] Jessica Barrett, Peter Diggle, Robin Henderson, and David TaylorRobinson. Joint modelling of repeated measurements and timetoevent outcomes: flexible model specification and exact likelihood inference. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77(1):131–148, 2015. ISSN 14679868. doi: 10.1111/rssb.12060. URL http://dx.doi.org/10.1111/rssb.12060.
 Belitz et al. [2013] Christiane Belitz, Andreas Brezger, Thomas Kneib, Stefan Lang, and Nikolaus Umlauf. BayesX: Software for Bayesian Inference in Structured Additive Regression Models, 2013. URL http://www.BayesX.org/. Version 2.1.
 Denaro and Bailey [2012] Kameryn Denaro and Barbara Bailey. Bayesian quantile regression for cystic fibrosis lung function data. at conference JSM 2012, 2012.
 Farcomeni and Viviani [2015] Alessio Farcomeni and Sara Viviani. Longitudinal quantile regression in the presence of informative dropout through longitudinal survival joint modeling. Statistics in Medicine, 34(7):1199–1213, 2015. ISSN 10970258. doi: 10.1002/sim.6393. URL http://dx.doi.org/10.1002/sim.6393.
 Faucett and Thomas [1996] Cheryl L. Faucett and Duncan C. Thomas. Simultaneously modelling censored survival data and repeatedly measured covariates: A gibbs sampling approach. Statistics in Medicine, 15(15):1663–1685, 1996. ISSN 10970258. doi: 10.1002/(SICI)10970258(19960815)15:15¡1663::AIDSIM294¿3.0.CO;21. URL http://dx.doi.org/10.1002/(SICI)10970258(19960815)15:15<1663::AIDSIM294>3.0.CO;21.
 Fenske et al. [2011] Nora Fenske, Thomas Kneib, and Torsten Hothorn. Identifying risk factors for severe childhood malnutrition by boosting additive quantile regression. Journal of the American Statistical Association, 106:494–510, 2011.
 Geraci [2014] Marco Geraci. Linear quantile mixed models: The lqmm package for laplace quantile regression. Journal of Statistical Software, 57(1):1–29, 2014. ISSN 15487660.
 Klein et al. [2015] Nadja Klein, Thomas Kneib, and Stefan Lang. Bayesian generalized additive models for location, scale and shape for zeroinflated and overdispersed count data. Journal of the American Statistical Association, 110:405–419, 2015.
 Koenker and Bassett [1978] R Koenker and G Bassett. Regression quantiles. Econometrica, 46:33–50, 1978.
 Koenker [2004] Roger Koenker. Quantile regression for longitudinal data. Journal of Multivariate Analysis, 91:74–89, 2004.
 Konstan et al. [2007] M.W. Konstan, W.J. Morgan, S.M. Butler, D.J. Pasta, M.L. Craib, S.J. Silva, D.C. Stokes, M.E. Wohl, J.S. Wagener, W.E. Regelmann, and C.A. Johnson. Risk factors for rate of decline in forced expiratory volume in one second in children and adolescents with cystic fibrosis. J Pediatr, 151:134–9, 139, 2007.
 Konstan et al. [2012] M.W. Konstan, Wagener J.S., D. R. Vandevanter, D.J. Pasta, A. Yegin, L. Rasouliyan, and W.J. Morgan. Risk factors for rate of decline in fev1 in adults with cystic fibrosis. J Cyst Fibros, 11:405–11, 2012.
 Kozumi and Kobayashi [2011] Hideo Kozumi and Genya Kobayashi. Gibbs sampling methods for bayesian quantile regression. Journal of Statistical Computation and Simulation, 81:1565–1578, 2011.
 Kulich et al. [2005] Michal Kulich, Margaret Rosenfeld, Jonathan Campbell, Richard Kronmal, Ron L. Gibson, Christopher H. Goss, and Bonnie Ramsey. Diseasespecific reference equations for lung function in patients with cystic fibrosis. American Journal of Respiratory and Critical Care Medicine, 172:7:885–891, 2005.
 Özgür et al. [2015] Asar Özgür, J. Ritchie, P. Kalra, and Diggle PJ. Joint modelling of repeated measurement and timetoevent data: an introductory tutorial. International Journal of Epidemiology, 44:334–344, 2015.
 R Development Core Team [2011] R Development Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2011. URL http://www.Rproject.org. ISBN 3900051070.
 Rizopoulos [2012] D. Rizopoulos. Joint Models for Longitudinal and TimetoEvent Data: With Applications in R. Chapman & Hall/CRC Biostatistics Series. Taylor & Francis, 2012. ISBN 9781439872864. URL https://books.google.de/books?id=xotIpb2duaMC.
 Rosenfeld et al. [2012] M. Rosenfeld, J. Emerson, S. Mcnamara, V. Thompson, B. W. Ramsey, W. Morgan, and R. L. Gibson. Risk factors for age at initial pseudomonas acquisition in the cystic fibrosis epic observational cohort. J Cyst Fibros, 11:446–53, 2012.
 Rue et al. [2009] Havard Rue, Sara Martino, and Nicolas Chopin. Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations. Journal of the Royal Statistical Society Series B, 71(2):319–392, 2009.
 Salvatore et al. [2011] D. Salvatore, R. Buzzetti, E. Baldo, M. P. Forneris, V. Lucidi, D. Manunza, I. Marinelli, B. Messore, A. S. Neri, V. Raia, M. L. Furnari, and G. Mastella. An overview of international literature from cystic fibrosis registries. part 3. disease incidence, genotype/phenotype correlation, microbiology, pregnancy, clinical complications, lung transplantation, and miscellanea. J Cyst Fibros, 10:71–85, 2011.
 Salvatore et al. [2012] D. Salvatore, R. Buzzetti, E. Baldo, M. L. Furanri, V. Lucidi, D. Manunza, I. Marinelli, B. Messore, A. S. Neri, V. Raia, and G. Mastella. An overview of international literature from cystic fibrosis registries. part 4: update 2011. J Cyst Fibros, 11, 2012.
 Schluchter et al. [2002] M. D. Schluchter, M. W. Konstan, and P. B. Davis. Jointly modelling the relationship between survival and pulmonary function in cystic fibrosis patients. Stat Med, 21:1271–87, 2002.
 TaylorRobinson et al. [2012] D. TaylorRobinson, M. Whitehead, F. Diderichsen, H. V. Olesen, T. Pressler, R. L. Smyth, and P. Diggle. Understanding the natural progression in with cystic fibrosis: a longitudinal study. Thorax, 67:860–6, 2012.
 TaylorRobinson et al. [2013] D. TaylorRobinson, R. L. Smyth, P. Diggle, and M. Whitehead. The effect of social deprivation on clinical outcomes and the use of treatments in the uk cystic fibrosis population: A longitudinal study. Lancet Respir Med, 1:121–128, 2013.
 Tsiatis and Davidian [2004] A.A. Tsiatis and M.A. Davidian. Joint modelling of longitudinal and timetoevent data: an overview. Statistica Sinica, 14:809–834, 2004.
 van Sickle et al. [2011] David van Sickle, Sheryl Magzamen, and John Mullahy. Understanding socioeconomic and racial differences in adult lung function. American Journal of Respiratory and Critical Care Medicine, 2011.
 Waldmann et al. [2013] Elisabeth Waldmann, Thomas Kneib, Stefan Lang, Yu Yue, and Claudia Flexeder. Bayesian semiparametric additive quantile regression. Statistical Modelling, 13:223–252, 2013.
 Yue and Rue [2011] Yu Yue and Harvard Rue. Bayesian inference for additive mixed quantile regression models. Computational Statistics and Data Analysis, 55:84–96, 2011.