Joint modeling of longitudinal and time to event data have been widely used for analyzing medical data, where longitudinal measurements is gathered with a time to event or survival data. In most of these studies, distributional assumption for modeling longitudinal response is normal, which leads to vulnerable inference in the presence of outliers in longitudinal measurements and violation of this assumption. Violation of the normality assumption can also make the statistical inference vague. Powerful distributions for robust analyzing and relaxing normality assumption, are skew-elliptical distributions, which include univariate and multivariate versions of the student’s t, the Laplace and normal distributions. In this paper, a linear mixed effects model with skew-elliptical distribution for both random effects and residuals and a Cox’s model for time to event data are used for joint modeling. This strategy allows for the skewness and the heavy tails of random effect distributions and thus makes inferences robust to the violation. For estimation, a Bayesian parametric approach using Markov chain Monte Carlo is adopted. The method is illustrated in a real Intensive Care Unit (ICU) data set and the best model is selected using some Bayesian criteria for model selection.

Keywords: joint models, bayesian approach, cox’s proportional model, longitudinal data, time to event data, markov chain Monte Carlo, skew-elliptical distributions

Data collected in many clinical and epidemiologic studies, contain both the longitudinal measurements and time to event data. In these studies, usually a biological marker, for example CD4 count measurements in AIDS clinical trials, is considered as a predictor of survival. Also, an interest event such as death or disease progression is considered as an important part of the study. These two types of outcomes often are analyzed using joint modeling of longitudinal and time to event data. Joint modeling has become a progressively popular approach to determine the relationship between of these processes.

Development of joint modeling of longitudinal and time to event data have been discussed broadly in the literature. Hogan and Laird,¹ Tsiatis and Davidian,² Ibrahim et al.³ and McCrink et al.⁴ give so excellent reviews of models and methods for joint analysis of this type of responses. The most common approach is to suppose that a random effect underlines both longitudinal and survival outcomes as shared parameter models (Henderson et al.⁵; Hashemi et al.⁵). Also some other researches are done on Bayesian method (Ibrahim et al.⁷; Chi and Ibrahim⁸), nonparametric random effect distributions (Wang and Taylor⁹; Brown and Ibrahim³), cure fractions (Yu et al.¹⁰ Chi and Ibrahim⁸), multiple longitudinal variables (Lin et al.¹¹), count data (Dunson and Herring¹²), zero-inflated outcomes (Rizopoulos et al.¹³), competing risks (Elashoff et al.¹⁴), accelerated failure time (Tseng et al.¹⁵), parametric assumptions by allowing flexible longitudinal trends (Brown et al.¹⁶). In joint modeling of longitudinal and survival data a mixed effects model is often used for analyzing longitudinal response of joint modeling. This part maybe incomplete due to dropout, where dropout mechanism (Diggle and Kenward¹⁷) can be considered or with some reasons ignored.

The traditional treatment, in many statistical models in particular mixed effect models contexts, is that the random effects follow the multivariate normal distribution (Verbeke and Lesaffre,¹⁸ Rosa et al.¹⁹). In practical applications, this assumption is likely to be failed if, for example, potential outliers exist in the data set.

The normality may be a reasonable model for many problems with continuous measurements, in other situations the normality assumption on a transformed scale is, at best, hoped for. Although such methods may give sensible empirical results, they should not be used if a more suitable theoretical model can be found. Some reasons for this are: (i) transforming usually prepares reduced information on an underlying data generation mechanism; (ii) the transformations are commonly applied to each component one by one, and joint normality is not secured; (iii) the interpretation based on transformed variables are more difficult, especially when transformations are various for the different proportions; (iv) multivariate homogeneity often requires a different transformation from the one to get normality; and (v) a model derived for a given data set may not be relevant to subsequent data sets.

Parametric robust analyses are statistical methods based non normal assumption for adjustment of usual statistical methods in the presence of outliers. The use of heavy-tailed distribution in mixed effects model for robust analysis is valuable for adjusting the role of outlier individuals in the sample and is widely applied in the literature, for example: Rosa et al.,¹⁹ Lange et al.²⁰] and Wu.²¹

A common approach for robust inference is to replace the multivariate t-distribution, which has heavier tails than the normal (Lange and Sinsheimer,²² Pinheiro et al.²³). There are, however, many other situations when the underlying distribution no longer satisfies the symmetric property (e.g., Zhang and Davidian).²⁴ Robust joint modeling can be found in Li et al.²⁵ and Huang et al.²⁶ where a student’s t distribution in different structures of joint modeling of longitudinal and survival data is applied. Thus alternative flexible distributions, such as the multivariate skew-elliptical (SE) distribution, are proposed in recent studies (Azzalini and Capitanio,²⁷ Ma et al.,²⁸).

In this paper, we have developed robust inference of joint modeling of longitudinal and time to event data using skew-elliptical distributions (Lange and Sinsheimer).²² These distributions include the skew student’s t (ST), the skew slash (SS), skew Laplace (SL) and the skew normal (SN) distributions. We consider a linear mixed effect model with skew-elliptical distribution assumption for longitudinal response modeling. Also, a Cox proportional hazard model with step baseline hazard in a frailty model structure is chosen for survival response modeling.

We have used Bayesian approach and the available software OpenBUGS (Spiegelhalter et al.)²⁹ for implementation of the models, where Bayesian criteria of DIC (deviance information criterion has been used for model selection. In application section, we have reanalyzed an ICU data set and we present a strategy, in a Bayesian perspective, and show that these skewed distributions are more appropriate than the commonly used distributions.

The rest of the article is organized as follows: In Section 2, some SE distributions are briefly introduced. In Section 3, the specification of joint models as a shared parameter model and Bayesian hierarchical joint models when the random effect terms follow SE distributions. This is done by the Gibbs sampling implementation. In Section 4, we describe the real data set on the hospitalized patients in ICU and subsequently present a Bayesian analysis of the data set. The McMC scheme and statistical interpretation are given in Section 5. Finally, discussions and concluding remarks are illustrated in Section 6.

Several versions of the multivariate SE distributions have been offered in the literature and some of them have recently been used in the univariate mixed effect models (Gosh et al.,³⁰ Jara et al.³¹). Here, we utilize the multivariate SN and ST distributions introduced by Sahu et al.³² and extended by others (Arellano-Valle and Genton).³³ We also make use of the multivariate Laplace distribution and its skewed version proposed by Arsla.³⁴ These distributions are shown to be a scaled mixture of normal that produces flexible distributions with, for example, heavier tails than the multivariate normal distribution.

The multivariate SN distribution

In general, the p-variate Y follows the SN distribution with location vector μ, scale matrix ∑ , and $p\times p$  skewness matrix $D=\text{diag}\left({\text{δ}}_1,{\text{δ}}_2,\dots ,{\text{δ}}_{\text{p}}\right)$ if the probability density function (pdf) of Y is given by

$f\left(y\text{|}\mu ,{{\displaystyle \sum }}^{\text{}},D\right)=2^p{}_p{\left\{y\text{|}\mu ,\right\}}_p\left\{{\left(I-D^{-1}D\right)}^{-1/2}{D}^1\left(y\mu \right)\right\}$  (1)

where $\Omega ={{\displaystyle \sum }}^{\text{}}+D^2$ , ${\varphi }_p\left\{y\text{|}\mu ,\Omega \right\}$  is the pdf of $N_p\left(\mu ,\Omega \right)$ evaluated at y and ${\Phi }_p$  denotes the cumulative density function (cdf) of the p-variate standard normal distribution (Sahu et al. [32]). We denote the density by $S{N}_p\left(\mu ,{{\displaystyle \sum }}^{\text{}},D\right)$ . If D = 0 density (2) reduces to the usual symmetric multivariate normal distribution. The mean and variance of Y are given by $E\left(Y\right)=\mu +\delta \sqrt{2/\pi }$ and $Var\left(Y\right)={{\displaystyle \sum }}^{\text{}}+\left(I-2/\pi \right)D^2$ ,

where ${\delta }^{\prime }=\left({\delta }_1,\dots .,{\delta }_p\right)$ .

The multivariate ST distribution

The p-variate Y follows the ST distribution with location vector μ, scale matrix ∑ , degree of freedom ν, and $p\times p$ skewness matrix $D=\text{diag}\left({\text{δ}}_1,{\text{δ}}_2,\dots ,{\text{δ}}_{\text{p}}\right)$ if the pdf of Y is given by

$f\left(y|\mu ,{{\displaystyle \sum }}^{\text{}},D,v\right)=2^p{t}_{p,v}(y|\mu ,)T_{p,v+p}\left\{{\left(\frac{v+g}{v+p}\right)}^{-1/2}{\left(I-D^{-1}D\right)}^{1/2}{D}^1\left(y\mu \right)\right\},$     (2)      

Where $g={\left(y-\mu \right)}^{\text{'}}{}^{-1}\left(y-\mu \right)$ , $t_{p,v}(y|\mu ,)$ is the pdf of student’s evaluated at y, and $T_{r,v+r}$  denotes the cdf of the p-variate standard t distribution with degrees of freedom $v+r$  (Sahu et al. [32]). We denote the density by $S{T}_p\left(\mu ,{{\displaystyle \sum }}^{\text{}},D,v\right)$ . As SN, If D = 0, then density (3) reduces to the usual symmetric multivariate student’s t distribution. As ν tends to infinity and D = 0, then the distribution reduces to the multivariate normal. The mean and the variance of Y are shown to be

$E(Y) = \mu +{(\frac{v}{\pi })}^{(1/2)\frac{\Gamma ((v-1)/2)}{\Gamma (v/2)}}\delta$                                and

$Var\left(Y\right) = \frac{v}{v-2}\frac{v}{\pi }{\left[\frac{\text{Γ}(\left(v-1\right)/2)}{\text{Γ}(v/2)}\right]}^2{D}^2 for n>2.$

The multivariate SL distribution

The p-variate Y follows the SL distribution with location vector μ, scale matrix ∑ , and $p\times p$ skewness vector $\delta =\left({\delta }_1,{\delta }_2,\dots ,{\delta }_p\right)\text{'}$  if the pdf of Y is given by

$f\left(y\text{|}\mu , {{\displaystyle \sum }}^{\text{}},\delta \right) = \frac{{\left|{{\displaystyle \sum }}^{\text{}}\right|}^{-1/2}}{2^p{\pi }^{\left(p-1\right)/2}\alpha \text{Γ}(\left(p+1\right)/2)}\text{exp} \left[{\left(y-\mu \right)}^{\text{'}}{{{\displaystyle \sum }}^{\text{}}}^{-1}\delta \alpha \sqrt{{\left(y\mu \right)}^{\text{'}}{{{\displaystyle \sum }}^{\text{}}}^{-1}\left(y\mu \right)}\right]$      (3)           

Where $\alpha =\sqrt{1+{\delta }^{\prime }{{{\displaystyle \sum }}^{\text{}}}^{-1}\delta }$  (Arslan, [31]). We denote the density by $S{L}_p\left(\mu ,{{\displaystyle \sum }}^{\text{}},D,v\right)$ . If D = 0, then (4) reduces to the density function of the multivariate symmetric Laplace distribution. The mean and the variance of Y are shown to be $E\left(Y\right)=\mu +\left(p+1\right)\delta$ and $Var\left(Y\right)=(p+1)(\Sigma +22\delta {\delta }^{\prime })).$ .

The longitudinal part of the joint model can be described as follows. Let $y_i\left(t\right)$ denote the value of longitudinal measurement, at time point t for the $i^{th}$ individual $i:1,2,\dots ,m$ . The observed times are $t_{ij}$ , $i:1,2,\dots ,m$ . Thus, the observed longitudinal data for the $i^{th}$ individual consist of measurements $y_i=\left\{y_i\left(t_{ij}\right),j:1,2,\dots ,n_i\right\}$ . For the longitudinal process, we consider the following linear mixed effect model

$y_{ij}=x_1^{\text{'}}\left(t_{ij}\right){\beta }_1+z_1^{\text{'}}\left(t_{ij}\right)b_{1i}+{\epsilon }_{ij}; i:1,2,\dots ,m, j:1,2,\dots ,n_i,$       (4)

where components of ${\epsilon }_i=\left({\epsilon }_{i1},\dots ,{\epsilon }_{i{n}_i}\right)\text{'}$ are measurement errors, ${\beta }_1=\left({\beta }_{11},\dots ,{\beta }_{1p_1}\right)\text{'}$ is a $p_1$ -dimensional vector of longitudinal fixed-effect parameters. $b_{1i}=\left(b_{1i1},\dots ,b_{1i{q}_1}\right)\text{'}$ is a $q_1$ -dimensional vector of random effects and is independent of ${\epsilon }_i$ . $x_1^{\text{'}}\left(t_{ij}\right)$  and $z_1^{\text{'}}\left(t_{ij}\right)$  are, respectively, $p_1$ -dimensional and $q_1$ -dimensional explanatory variables. Let $U_i$ denote the observed survival time for the $i^{th}$ individual, $i:1,2,\dots ,m$ , which is taken as the minimum of the true event time $U_i^{*}$ and the censoring time $C_i$ , i.e., $U_i=\text{min}\left(U_i^{*},C_i\right)$ . We define a censoring indicator, ${\delta }_{it}=I\left(U_i^{*}\le C_i\right)$ , which is 0 for right-censored and 1 for complete observed individuals. Therefore, the observed data for the time to event consist of the pairs $\left\{\left({\text{U}}_i,{\text{ δ}}_{it}\right),\text{ i }=1,2,\mathrm{...},\text{m}\right\}$ .

For the time to event process, we consider a semiparametric proportional hazard model in a frailty structure.

The hazard function in our proposed model is given by:

$h\left(u_i\text{|}{x}_{2i},z_{2i},b_{2i}\right)=h_0\left(u_i\right)\exp \left\{x_{2i}^{\text{'}}{\beta }_2+z_{2i}^{\text{'}}{b}_{2i}\right\},$ (5)

where $h_0\left(u_i\right)$  is the baseline hazard function. Thus, the density function of survival time for $i^{th}$  individual can be written in the form:

$h^{{\delta }_{it}}(u_i|x_{2i},z_{2i},b_{2i}) \times \text{exp} \left\{-H_0\left(u_i\right)\exp \left\{x_{2i}^{\text{'}}{\beta }_2+z_{2i}^{\text{'}}{b}_{2i}\right\}\right\}$  (6)

Where $H_0\left(w\right)=\underset{0}{\overset{w}{{\displaystyle \int }}}{h}_0\left(u\right)du$ , $x_2^{\text{'}}$ and $z_2^{\text{'}}$ are $p_2$ and $q_2-$ dimensional vectors of explanatory variables, respectively. ${\beta }_2=\left({\beta }_{21},\dots ,{\beta }_{2p_2}\right)\text{'}$ is a $p_2$ -dimensional vector of time to event fixed effect parameters, and $b_{2i}=\left(b_{2i1},\dots ,b_{2i{q}_2}\right)\text{'}$  is a $q_2-$ dimensional of random effects of time to event process. It is important to note that some elements of $b_{1i}$  and $b_{2i}$ are shared between two models and joint modeling is formed using this structure.

Suppose that $b_1=\left(b_{11}^{\text{'}},b_{12}^{\text{'}},\dots ,b_{1m}^{\text{'}}\right)\text{'}$ , $b_2=\left(b_{21}^{\text{'}},b_{22}^{\text{'}},\dots ,b_{2m}^{\text{'}}\right)\text{'}$ , $t=\left(t_1,t_2,\dots ,t_m\right)\text{'}$ .

We further assume that the baseline hazard is a step function, $h_o\left(u\right)=h_k$ , for $s_{k-1}<u<s_k,\text{ k }=1,2,\mathrm{...},\text{K}$ , where $0<s_1<{\text{s}}_2<\dots <s_K<\infty$ is a partition of (0,∞) and K indicates the number of steps for the baseline hazard, also $s_0 = o$ . Therefore the cumulative baseline hazard is given by

$H_0\left(u\right)=\left(h_j\left(u-s_{j-1}\right)+{\displaystyle \sum }_{i=1}^{j-1}{h}_i\left(s_i-s_{i-1}\right)\right)I\left(u\in \left(s_{i-1},s_j\right]\right).$ (7)

We implement the Bayesian methodology using MCMC techniques for the joint modeling of longitudinal and survival data with SE distributional assumption. ST distribution was used for random effects terms distribution in tow longitudinal and survival models because of more parameters than other distribution. Inference based on the another distribution will be similar. Multivariate SE distributions will be reconstructed by the mixing strategy approach, therefore they are very convenient to implement in the Bayesian framework. A key characteristic of this model, which allows writing BUGS codes, is that it can be expressed in a flexible hierarchical representation. For example by using multivariate ST distribution for $b_{ik}~ST\left(0,{\Psi }_k,D_{b_k}\right), k=1,2$ , hierarchical representation can be followed:

$b_{ik}|{\Psi }_k,D_{b_k},w_{b_{ik}}{~}^{iid}{N}_{q_k}\left(D_{b_k}{w}_{b_{ik}},s_{b_{ik}}^{-2}{\Psi }_k\right)$

$w_{b_{ik}}|s_{b_{ik}}^{-2}{~}^{iid}{N}_{q_k}\left(0,s_{b_{ik}}^{-2}{I}_{q_k}\right)I\left(w_{b_{ik}}>0\right)$

$s_{b_{ik}}^2|{\nu }_{b_k}~Gamma\left(\frac{{\nu }_{b_k}}{2},\frac{{\nu }_{b_k}}{2}\right)$ ……………. (8)

Bayesian specification of the model needs to consider prior distribution for all the unknown parameters. In our modeling $\theta =\left({\beta }_1^{\text{'}},{\beta }_2^{\text{'}},{\Psi }_1,{\Psi }_2,{\sigma }_e^2,{\nu }_{b_1},{\nu }_{b_2}\right)\text{'}$ is the unknown vector of parameters and because of not having any prior information from historical data or previous experiment, we try to assign non-informative prior distributions for the parameters. Assuming elements of the parameter vector to be independent, the prior distributions are given by: 

${\beta }_1~N_{p_1}\left({\beta }_{01},S_{{\beta }_1}\right), {\beta }_2~N_{p_2}\left({\beta }_{02},S_{{\beta }_2}\right), {\sigma }_e^2~IG\left(\frac{{\tau }_e}{2},\frac{T_e}{2}\right)$ …. (9)

${\Psi }_1~I{W}_{\tau b_1}\left(T_{b_1}\right), {\Psi }_2~I{W}_{\tau b_2}\left(T_{b_2}\right),$

${\nu }_{b_1}~\exp \left({\lambda }_{b_1}\right)I\left\{{\nu }_{b_1}>0\right\}, {\nu }_{b_2}~\exp \left({\lambda }_{b_2}\right)I\left\{{\nu }_{b_2}>0\right\},$

${\delta }_{b_1}~N_{q_1}\left({\mu }_{b_1},{\gamma }_{b_1}\right)I\left\{{\delta }_{b_1}>0\right\}, {\delta }_{b_2}~N_{q_2}\left({\mu }_{b_2},{\gamma }_{b_2}\right)I\left\{{\delta }_{b_2}>0\right\}$

where ${\sigma }_e^2$  follows inverted gamma, ${\beta }_1$ and ${\beta }_2$  are assumed to be normally distributed, ${\Psi }_1$ and ${\Psi }_2$  were denoted by inverse-Wishart distribution, ${\delta }_{b_1}$ and ${\delta }_{b_2}$ distributions are assumed truncated exponential distribution and in final truncated normal distribution for ${\delta }_{b_1}$ and ${\delta }_{b_2}$ are determined. The hyper parameters of these priors are selected so that they lead to the non-informative prior distributions.

According to $y_{\text{i}}|{\beta }_1,b_{i1},{\sigma }_e^2{~}^{\text{i}i\text{d}}N\left(x_{1i}^{\text{'}}{\beta }_1+z_{1i}^{\text{'}}{b}_{1i},{\sigma }_e^2{I}_i\right),$  the joint posterior density of all unobservable is given by :

$\begin{array}{l}\left(\theta ,b_{1i},b_{2i},W_{b_{i1}},W_{b_{i2}}|y,t\right)\\ \propto \left(\varphi \left({\text{y}}_{\text{i}}|x_{1i}^T{\beta }_1+z_{1i}^T{b}_{1i},{\text{σ}}_{\text{e}}^2{\text{I}}_{\text{i}}\right)\times {\text{φ}}_{q_1} \left(b_{1i}|{\mu }_{b_1}+D_{b_1}{\text{w}}_{b_1},{\text{s}}_{b_{1i}}^{-2}{\Psi }_1\right)\times \varphi \left({\text{w}}_{b_1}|0,{\text{s}}_{b_{1i}}^{-2}{\text{I}}_{q_1}\right)\text{I}\left({\text{w}}_{b_1}>0\right)\right)\\ \times {\displaystyle \prod }_{i=1}^m\left\{h^{{\delta }_{it}}\left(t_i|x_{2i},z_{2i},b_{2i}\right) \times \exp \left\{-H_0\left(t_i\right)\exp \left\{x_{2i}^T{\beta }_2+z_{2i}^T{b}_{2i}\right\}\right\}\right\}\\ \times {\text{φ}}_{q_2}\left(b_{2i}|{\mu }_{b_2}+D_{b_2}{\text{w}}_{b_2},{\text{s}}_{b_{2i}}^{-2}{\Psi }_2\right)\times \varphi \left({\text{w}}_{b_2}|0,{\text{s}}_{b_{2i}}^{-2}{\text{I}}_{q_2}\right)\text{I}\left({\text{w}}_{b_2}>0\right) \times \pi \left(\theta \right) \end{array}$ …. (11)       

The above described joint posterior distribution is analytically intractable but MCMC methods such as the Gibbs sampler and Metropolis-Hastings algorithm can be used to draw samples, from which features of marginal posterior distribution of interest can be inferred. The Gibbs sampler works by drawing samples iteratively from conditional posterior distributions deriving from (11). For this purpose, we have needed all full conditional distribution. Let $\vartheta =\left({\theta }^{\prime },b_1^{\text{'}},b_2^{\text{'}},,W_{b_{i1}},W_{b_{i2}}\right)\text{'}$ , and let $\xi$ is one of the component of it, we define ${\vartheta }_{\left(-\xi \right)}$  for the above-mentioned vector when v is omitted from it. The full conditional distributions for ${\beta }_1$ , ${\beta }_2$   ${\delta }_{b_1}$ and ${\delta }_{b_2}$ , respectively, have the form:

$\begin{array}{l}\pi ({\beta }_2|{\theta }_{\left(-{\beta }_2\right)},y,t) \propto {\displaystyle \prod }_{i=1}^m\left\{h^{{\delta }_{it}}\left(t_i|x_{2i},z_{2i},b_{2i}\right) \times \exp \left\{-H_0\left(t_i\right)\exp \left\{x_{2i}^T{\beta }_2+z_{2i}^T{b}_{2i}\right\}\right\}\right\} \times \\ \text{exp}\left\{\frac{-1}{2}\left({\beta }_2-{\beta }_{02}\right)S_{{\beta }_2}^{-1}\left({\beta }_2-{\beta }_{02}\right)\right\}\end{array}$

and

$\begin{array}{l}\pi (b_{2i}|{\theta }_{\left(-{\beta }_2\right)},y,t)\\ \propto {\displaystyle \prod }_{i=1}^m\left\{h^{{\delta }_{it}}\left(t_i|x_{2i},z_{2i},b_{2i}\right) \times \exp \left\{-H_0\left(t_i\right)\exp \left\{x_{2i}^T{\beta }_2+z_{2i}^T{b}_{2i}\right\}\right\}\right\}\\ \times \exp \left\{\frac{-s_{b_{i2}}^2}{2}{\left(b_{i2}-D_{b_2}{w}_{b_{i2}}+{\mu }_{b_2}\right)}^T{\Psi }_2{}^{-1}\left(b_{i2}-D_{b_2}{w}_{b_{i2}}+{\mu }_{b_2}\right)\right\}\end{array}$

and

$\begin{array}{l}\pi ({\delta }_{b_k}|{\theta }_{\left(-S_{b_{ki}}^2\right)},y,t)\propto f\left(b_i|{\Psi }_1,{\delta }_{b_k},w_{b_{ki}}\right)T\left({\delta }_{b_k}|{\mu }_{b_k},{\gamma }_{b_k}\right)\\ \propto \exp \left\{\frac{-s_{b_i}^2}{2}{\left(b_{ki}-D_{b_k}{W}_{b_{ki}}+{\mu }_{b_k}\right)}^T{\Psi }_1{}^{-1}\left(b_{i1}-D_{b_k}{W}_{b_{ki}}+{\mu }_{b_k}\right)\right\}\\ \exp \frac{-1}{2}{\left({\delta }_{b_k}-{\mu }_{b_k}\right)}^{\text{T}}{\gamma }_{b_k}^{-1}\left({\delta }_{b_k}-{\mu }_{b_k}\right)I\left\{{\delta }_{b_k}>0\right\}, k=1,2\end{array}$

For ${\beta }_1$  and $b_{1i}$ , i = 1, 2, ...,n the full conditional distributions have closed form, and are given
by:

$\pi ({\beta }_1\text{|}{\vartheta }_{\left(-{\beta }_1\right)},y,t{\text{)~N}}_{p_1}\left(A_{{\beta }_1}^{-1}{a}_{{\beta }_1},A_{{\beta }_1}^{-1}\right),$

where   $A_{{\beta }_1}=\underset{}{\overset{}{{\displaystyle {\sum }_{i=1}^m}}}\frac{1}{{\text{σ}}_{\text{e}}^2}{x}_{1i}^{\text{'}}{I}_i^{-1}{x}_{1i}+S_{\beta }^{-1})$  and $a_{{\beta }_1}=S_{\beta }^{-1}{\beta }_{01}+\underset{}{\overset{}{{\displaystyle {\sum }_{i=1}^m}}}{\text{σ}}_{\text{e}}^2{x}_{1i}^{\text{'}}\left(y_i-z_{1i}^{\text{'}}{b}_{i1}\right)$ ,

also $\pi (b_{i1}|{\theta }_{\left(-b_{i1}\right)},y,t)~N_{q_1}\left(A_{b_{i1}}^{-1}{a}_{b_{i1}},A_{b_{i1}}^{-1}\right)$

where, $A_{(b_i)}^{(-1)}=s_{(b_{i}1)}^2{\psi }_1{}^{(-1)}+{\sigma }_e^{(-2)}{z}_i{1}^T{I}_i^{(-1)}{z}_{i}1$  and $a_{(b_{i}1)}={\sigma }_e^{(-2)}{z}_i{1}^T{I}_i^{(-1)}(y_i-x_{i}1{\beta }_1)+{\Psi }_1{}^{(-1)}(D_{(b_1)}{w}_{(b_{i}1)}-{\mu }_{(b_1)}).$
The full conditional distributions of ${\nu }_{b_1},{\nu }_{b_2}$ , $b_{2i},{\delta }_{b_1}$  ,and ${\delta }_{b_2}$ are not in closed forms of known distributions. Thus, the Gibbs sampling can readily be implemented. For these parameters a Metropolis-Hastings can be embedded in the Gibbs sampling scheme to obtain draws for them. The full conditional distributions of ${\Psi }_k,S_{{\text{b}}_{k2}}^2$ , $k=1,2$  and ${\sigma }_e^2$  are presented in Appendix1.

We analyzed the data from a retrospective cohort study. All consecutive adult patients (N=380) were considered who, from October 2010 through October 2012, underwent cardiac surgical procedure with cardiopulmonary bypass at Masih Daneshvari Hospital,Tehran, Iran as a referral center. Patients who underwent more than one cardiac surgical procedure during hospitalization (n=40), were excluded from the study. Additional exclusion criteria were surgeries performed off-pump (n=15) and preoperative renal failure requiring dialysis (n=25). Overall, 300 patients could be included in the study. After surgery, patients were admitted to the intensive care unit (ICU) and were followed from the day of ICU admission until ICU discharge or end of study. Urine output (UO) and other physiological variables were repeatedly measured in 2 hours in the first 8 hours of admission in ICU. Patients’ demographics and laboratory data were gathered repeatedly. Main survival endpoint was the time of occurrence of acute kidney injury (AKI) after cardiac surgery that was defined, the amount of urine output less than 0.5 ml/kg per hour in first 6 hours. Patients were followed since ICU admission until AKI occurrence.

We have used the joint modeling which fitted by Guo and Carlin,³⁵ but with some various distributional assumptions. The aim of fitting joint models is to additionally study the relationship between two longitudinal and time to event responses in terms of available covariates. Parameter  in our joint model is determined to assessment relationship strength between UO and AKI occur hazard. The joint modeling with random intercept and slop is given by:

$U{O}_{ij}={\beta }_{11}+{\beta }_{12}{t}_{ij}+{\beta }_{13}ag{e}_i+{\beta }_{14}Se{x}_i+{\beta }_{15}DB{P}_{ij}+{\beta }_{16}infectio{n}_{ij}+b_{1i}+b_{2i}{t}_{ij}+{\epsilon }_{ij}$

For the time to event process, we have applied a Cox proportional hazard model; the hazard function for this model is given by:

$h\left(t_i\right)=h_0\left(t_i\right)\exp \exp \left\{\begin{array}{l}{\beta }_{21}+{\beta }_{22}APACHEI{I}_i+{\beta }_{23}Se{x}_i+{\beta }_{24}DB{P}_i+{\beta }_{25}hospital\_sta{y}_i\\ +{\beta }_{26}emergency\_surger{y}_i+\gamma \left(b_{1i}+b_{2i}\right)\end{array}\right\}$
the random effects, $b_i=\left(b_{1i},b_{2i}\right)\text{'}$ , are shared between two models. Random effects are assumed to have a bivariate skew elliptical distribution, that is, $b_i~S{E}_2\left(0,\Sigma ,D\right)$ . In this model, $U{O}_{ij}$  is the $j^{t}h$ measurement urine output on the $i^{th}$  individual in the study, $j:1,2,\mathrm{...},n_i$ and i : 1, 2, ...,m. $Se{x}_i$ is a gender indicator (0=male, 1=female), also other three explanatory variables are $infectio{n}_{ij}$  (0=no, 1= yes), $DB{P}_{ij}$ (0=no, 1=yes), and $emergency\_surger{y}_i$  (0=no, 1=yes). $\text{APACHE II}$ score and DBP mean Acute Physiology and Chronic Health Evaluation II score and diastolic blood pressure.

In the Bayesian MCMC method, we ran two parallel MCMC chains with different starting values for each 100,000 iteration. Then, we discarded the first 40, 000 iterations as pre-convergence burn-in and retained 50, 000 as the posterior analysis.

 In all models, ${\beta }_k\sim N_{pk}\left(0,10000\right)$  , ${\sigma }_e\sim \Gamma \left(0.1,0.1\right)$ and $\gamma \sim N\left(0,100\right)$ where, k = 1, 2, p1 = 6, p2 = 6, and ${\Psi }_k$  ∼ IW(I2, 2). The hyper parameter of student’s t, ${\text{ν}}_k$ , is assumed to have U(1, 7) distribution. For the piecewise baseline hazard function [ $h_l$ , $l=1,2,3,4$ (the number of piecewise baseline=4)] the gamma (2, 2) prior distribution is considered.

Hyperparameters are chosen such that the priors of the parameters tend to be weakly informative. The implementation of this method is relatively easy in the publicly available software Open BUGS (Spiegelhalter et al. 2003). The parameter estimates are not sensitive to the choice of hyperparameters and initial values. In order to see how stable the final estimates are, multiple parallel chains with different initial values should be carried out (we considered two chains).

Figure 1 shows the observed longitudinal measures of UO plotted against time for the patients included in the analysis based on the AKI occurrence. UO decreases over time and it seems censoring mechanism be randomly, hence tow graphs have similar pattern.

In total, our results were organized on 10 different models, 9 models in all combinations of three distribution SN, ST and SL for for $b_i$  and ${\epsilon }_i$ , In addition to normal assumption for both of them table 1. In this table for example SN,SL model shows skew elliptical normal distribution for shared random effects and skew elliptical Laplace distribution for the residuals. The estimated random effects distribution surface and estimated residuals distribution histogram in N,N model, are depicted in Figure 2 (a) , (b) respectively. The plots clearly show deviation from bivariate normality distribution in shared random effects and normally distributed residuals. According to DIC criteria the model with bivariate skew t and skew Laplace distribution for shared random effects and residuals respectively (ST,SL), is the best-fitted model for the data, while (ST,ST) model is the second-best model.

Figure 2 (a) Estimated random effects distribution under the multivariate normal distribution fitted to the ICU data. (b) Histogram of residual estimates.

Table 2 includes parameter estimates and %95 highest posterior density intervals for parameters of (ST, SL) model for the ICU data set. According to the results of the mixed model of joint analysis, a decrease of UO occurred more in infected patients and also in low DBP, significantly. The total number of UO longitudinal measurements was 1540 and the average of measurements was 3.85 per patients. Low DBP (HR=1.2) and higher hospital length of stay (HR=2.0) had significant effects on the risk of AKI after cardiac surgery, similar results were obtained by other researchers such as: Kim and et al.,³⁶ Khwaja.³⁷ After using joint modeling, our findings showed a negative significant association between the risk of AKI following cardiac surgery and UO in patients ( $\gamma$ = -0.31), that expressed the suitable selection of the joint modeling for this research.

Joint modeling of longitudinal and survival data provides a suitable tools for the analysis of relationship between tow responses. In this paper, flexible joint modeling of these responses is proposed to deal with the situations in which the underlying distributions follow asymmetric structures. We fitted some flexible models for an empirical example taken from a retrospective cohort study to illustrate the suggested strategy. In most studies, the random effects and residuals terms vector is taken to be normally distributed. However, in the analysis of ICU data, we showed that this assumption was violated, leading to inappropriate results. Thus, we applied suitable non-normal distributions to avoid model mis-specifications.

We have used the Cox proportional hazard model for time to event process and a linear mixed effect model for longitudinal measurements. These two models have joined by a vector of random effects which shared in some components together.

In this study, we assessed the association between UO and AKI using the joint modeling of longitudinal and survival data for first time.

After using joint modeling, our findings showed a negative significant association between the risk of AKI following cardiac surgery and UO, which expressed the suitable selection of the joint modeling for this research. Also significant skewness parameter for random intercept implied accurate choice in skew distribution for random effects. According to the results of the mixed model of joint analysis, a decrease of UO occurred more in female and infected patients and also in low DBP, significantly.

In survival part, Female gender, older age, high APACHE II score, low DBP, long hospitalization and emergency surgery had significant effects on the risk of AKI after cardiac diseases. Similar results were obtained by other researchers such as: Kumar,³⁸ Bagshaw³⁹ and Kim.³⁶ Sepsis was also introduced as considerable risk factors of AKI. Majority of infected patients were septic in this study, Hashemian.⁴⁰ Many other factors are associated with an increased risk of AKI but their influence will be highly dependent on the specific nature of the population, Brochard et al.⁴¹

However, the occurrence of AKI after cardiac surgery is often not rapid and other factors concerning intra operative and postoperative management of patients could be relevant. We recommend definite prevention programs in the ICU target patients with traditional risks of AKI such as older age, sex, cardiovascular surgery, chronic pulmonary disease etc.^42–55

The authors are grateful to the reviewers for their constructive comments that improved the manuscript.

None.

1. Hogan JW, Laird NM. Model-based approaches to analyzing incomplete longitudinal and failure time data. Stat Med. 1997;16(1–3):259–272.
2. Tsiatis AA, Davidian M. A semi parametric estimator for the proportional hazards model with longitudinal covariates measured with error. Biometrika. 2001;88(2):447–458.
3. Brown ER, Ibrahim JG. A Bayesian semi parametric joint hierarchical model for longitudinal and survival data. Biometrics. 2003;59(2):221–228.
4. Mccrink K, Spelke ES, Dehaene S, et al. Non-symbolic halving in an Amazonian indigene group. Developmental Science. 2013;16(3):451–462.
5. Henderson R, Diggle P, Dobson A. Joint modeling of longitudinal measurements and event time data. Biostatistics. 2000;1(4):465–480.
6. Hashemi RH, Jacqmin-Gadda, Commenges D. A latent process model for joint modeling of events and markers. Lifetime Data Anal. 2003;8(4):335–348.
7. Ibrahim JG, Chen MH, Sinha D. Bayesian Survival Analysis. Springer- Verlag, US; 2002.
8. Chi YY, Ibrahim JG. Joint models for multivariate longitudinal and multivariate survival data. Biometrics. 2006;62(2):432-445.
9. Wang Y, Taylor JMG. Jointly modeling longitudinal and event time data with application to acquired immunodeficiency syndrome. Journal of the American Statistical Association. 2001;96(455):895–905.
10. Yu M, Law NJ, Taylor JMG. Joint longitudinal survival cures models and their application to prostate cancer. Statistica Sinica. 2004;14:835–862.
11. Lin H, McCulloch CE, Mayne ST. Maximum likelihood estimation in the joint analysis of time-to-event and multiple longitudinal variables. Stat Med. 2002;21(16):2369–2382.
12. Dunson DB, Herring AH. Bayesian latent variable models for mixed discrete outcomes. Biostatistics. 2005;6(1):11–25.
13. Rizopoulos D, Verbeke G, Lesaffre E, et al. A two-part joint model for the analysis of survival and longitudinal binary data with excess zeros. Biometrics. 2008;64(2):611–619.
14. Elashoff RM, Li G, Li N. A joint model for longitudinal measurements and survival data in the Presence of multiple failure types. Biometrics. 2008;4(3):762–771.
15. Tseng YK, Hseih F, Wang JL. Joint modeling of accelerated failure time and longitudinal data. Biometrika. 2005;92(3):587–603.
16. Brown ER, Ibrahim JG, DeGruttola V. A flexible B-spline model for multiple longitudinal biomarkers and survival. Biometrics. 2005;61(1):64–73.
17. Diggle PJ, Kenward MG. Informative dropout in longitudinal data analysis (with discussion). Journal of the Royal Statistical Society Series C (Applied Statistics). 1994;43(1):49–93.
18. Verbeke G, Lesaffre E. The effect of misspecifying the random-effects distribution in linear mixed models for longitudinal data. Computational Statistics & Data Analysis. 1997. p. 541–556.
19. Rosa GJM, Padovani CR, Gianola D. Robust linear mixed models with normal/independent distributions and Bayesian MCMC implementation. Biometrical Journal. 2003;45 (5):573–590.
20. Lange KL, Little JA, Taylor MGJ. Robust statistical modeling using the t distribution. Journal of the American Statistical Association. 1989;84:881–896.
21. Wu L. Mixed Effects Models for Complex Data. Chapman, Hall, editors. Taylor & Francis, CRC Press, USA; 2010.
22. Lange KL, Sinsheimer JS. Normal/independent distributions and their applications in robust regression. Journal of Computational and Graphical Statistics. 1993;2(2):175–198.
23. Pinheiro JC, Liu CH, Wu YN. Efficient algorithms for robust estimation in linear mixed-effects models using the multivariate t distribution. Journal of Computational and Graphical Statistics. 2001;10(2):249–276.
24. Zhang D, Davidian M. Linear mixed models with flexible distributions of random effects for longitudinal data. Biometrics. 2001;57(3):795–802.
25. Li N, Elashoff RM, Li G. Robust Joint Modeling of Longitudinal Measurements and Competing Risks Failure Time Data. Biom J. 2009;51(1):19–30.
26. Huang X, Li G, Elashoff RM. A joint model of longitudinal and competing risks survival data with heterogeneous random effects and outlying longitudinal. Stat Interface. 2010);3(2):185–195.
27. Azzalini A, Capitanio A. Statistical applications of the multivariate skew normal distribution. Series B Statistical Methodology. 1999;61(3):579–602.
28. Ma J, Plesken H, Treisman JE, et al. Lightoid and Claret: a rab GTPase and its putative guanine nucleotide exchange factor in biogenesis of Drosophila eye pigment granules. Proc Natl Acad Sci U S A. 2004;101(32):11652–11657.
29. Spiegelhalter DJ, Thomas A, Best N, et al. Win BUGS Examples, MRC Biostatistics Unit, Institute of Public Health and Department of Epidemiology and Public Health, Imperial College School of Medicine, UK, 2003.
30. Ghosh SK, Hajra S, Jayaram M. Faithful segregation of the multicopy yeast plasmid through cohesin-mediated recognition of sisters. Proc Natl Acad Sci U S A. 2007;104(32):13034–13039.
31. Jara A, Mezzano S. Vascular damage in chronic kidney disease. Rev Med Chil. 2008;136(11):1476–1484.
32. Sahoo R, Sengupta R, Ghosh S. Nitrosative stress on yeast: inhibition of glyoxalase-I and glyceraldehyde-3-phosphate dehydrogenase in the presence of GSNO. Biochem Biophys Res Commun. 2003;302(4):665–670.
33. Arellano-Valle RB, Genton MG. On fundamental skew distributions. Journal of Multivariate Analysis. 2005;96:93–116.
34. Arslan E, Mack W, Eraslan AN. Transient heating of a rotating elastic-plastic solid shaft. 2008;8(1)::0257–10258.
35. Guo X, Carlin BP. Separate and Joint Modeling of Longitudinal and Event Time Data Using Standard Computer Packages. The American Statistician. 2004;58(1):16–24.
36. Kim MY, Park JH, Kang NR, et al. Increased risk of acute kidney injury associated with higher infusion rate of mannitol in patients with intracranial hemorrhage. J Neurosurg. 2014;120(6):1340–1348.
37. Khwaja A. KDIGO clinical practice guidelines for acute kidney injury. Nephron Clin Pract. 2012;120(4):c179–c184.
38. Kumar AB, Suneja M. Cardiopulmonary bypass-associated acute kidney injury. Anesthesiology. 2011;114(4): 964–970.
39. Bagshaw SM, George C, Gibney RN, et al. A multi-center evaluation of early acute kidney injury in critically ill trauma patients. Ren Fail. 2008;30(6):581–589.
40. Hashemian S, Jamaati HR, Malekmohammad M, et al. Assessing the performance of two clinical severity scoring systems in the ICU of a tertiary respiratory disease center. Tanaffos. 2010;9(3):58–64.
41. Brochard L, Abroug F, Brenner M, et al. An official ATS/ERS/ESICM/SCCM/SRLF statement: prevention and management of acute renal failure in the ICU patient: an international consensus conference in intensive care medicine. Am J Respir Crit Care Med. 2010;181(10):1128–1155.
42. Acuna E, Rodriguez CA. Meta-analysis study of outlier detection methods in classification, Technical paper, University of Puerto Rico at Mayaguez, Retrieved from academic.uprm.edu/ eacuna/ paperout.pdf. In proceedings IPSI 2004, Italy. 2004.
43. Brooks SP. Discussion on the paper by Spiegelhalter Best Carlin and vande Linde, Journal of the Royal Statistical Society Series B 2002;64:3616–3618.
44. Carlin P, Louis TA. Bayes and Empirical Bayes Methods for Data Analysis. second edition. Chapman & Hall, CRC Press, New York; 2000.
45. Chen MH, Ibrahim JG, Sinha D. A new joint model for longitudinal and survival data with a cure fraction. Journal of Multivariate Analysis. 2004;91:18–34.
46. Chi YY, Ibrahim JG. Bayesian approaches to joint longitudinal and survival models accommodating both zero and nonzero cure fractions. Statistica Sinica. 2007;17:445–462.
47. Gelman A, Rubin DB. Inference from Iterative Simulation Using Multiple Sequences. Statistical Science. 1992;7(4):457–511.
48. Geweke J. Evaluating the accuracy of sampling-based approaches to calculating posterior moments. In: Bernado JM, editors. Bayesian Statistics, Clarendon Press, UK; 1992 p.
49. Geweke J. Bayesian treatment of the independent Student-t linear model. Journal of Applied Economic. 1993;8:19–40.
50. Goldman AI, Carlin BP, Crane LR, et al. Response of CD4+ and Clinical Consequences to Treatment Using ddI or ddC in Patients with Advanced HIV Infection. J Acquir Immune Defic Syndr Hum Retrovirol. 1996;11(2):161–169.
51. Heidelberger P, Welch PD. Simulation run length control in the presence of an initial transient. Operations Research. 1983;31(6):1109–1144.
52. Ibrahim JG, Chen MH, Sinhav D. Bayesian methods for joint modeling of longitudinal and survival data with application to cancer vaccine trials. Statistica Sinica. 2004;14:863–883.
53. Iglewics B, Martinez J. Outlier Detection using robust measures of scale. Journal of Statistical Computation and Simulation. 1982;15(4):285–293.
54. Little RJA. Robust Estimation of the mean and covariance matrix from data with missing values. Journal of the Royal Statistical Society. Series C (Applied Statistics). 1988;37(1):23–38.
55. Spiegelhalter DJ, Best NG, Carlin BP, et al. Bayesian measures of model Complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2002;64(4):583–639.