Bayesian Joint Modelling of Survival Time and Longitudinal CD4 Cell Counts Using Accelerated Failure Time and Generalized Error Distributions

Show more

1. Introduction

Survival of HIV/AIDS patients is crucially dependent on comprehensive and targeted medical interventions. Health professionals monitor patients’ health status using such disease markers as CD4 T-cells counts. The disease progression as indicated by the longitudinal CD4 measures may affect the time of an event of interest―death of a patient in this case. The main interest of inference is on the association between the longitudinal and survival processes. Joint models for longitudinal and time-to-event are based on the joint distribution of the two processes [1] . The joint analysis may be appropriate when the longitudinal variable is correlated with patient’s health status and incorporate all information simultaneously so as to provide valid and efficient inferences [2] .

The traditional approach in the analysis of survival data assumes a homogeneous population, where all individuals have same health risks. In practice, individual patients possibly differ in health risks such as their vulnerability to causes of death, responses to treatments, and influences of various risk factors. Joint modelling of the two data often assumes normal distributions in the linear mixed models [2] [3] [4] . It is interesting to look for alternative distributions that can accommodate data that may not be normally distributed.

The current study considers the longitudinal measure in terms of its rate of growth. The rate of growth is an important concept in studying changes. If the level of growth is viewed as the current status of a process at a specific time, the rate of growth measures how fast the process is changing at that time [4] . Studies by [4] extend the usual growth models using the generalized error distribution (GED) and estimate its parameters under Bayesian framework. The author studied such a general form of linear growth model ${Y}_{it}={X}_{t}{\beta}_{i}+{\in}_{it}$, where ${Y}_{it}$ is growth observation for individual i at time t assuming the error term to have the generalized error distribution ${\in}_{it}~\text{GED}\left(0,{\sigma}_{\u03f5}^{2},{\gamma}_{t}\right)$ and independent of the random effects ${\beta}_{i}$ . The GED has a shape parameter ${\gamma}_{t}$ that may possibly vary across time, but here assumed constant. The model can handle data with both leptokurtic and platykurtic errors.

Markos et al. [5] studied joint modelling of survival time and longitudinal CD4 cell counts of HIV/AIDS patients using Bayesian methods. The authors compared various Bayesian joint models involving the Weibull, lognormal and loglogistic AFT distributions and normality assumption for the longitudinal CD4 measure using two data sets. They recommended the Bayesian joint loglogistic model for one data set collected from Shashemene referral hospital and the Bayesian joint lognormal model for the second data set collected from Bale Robe general hospital. These models have same hazard rate functions as that of data sets. In the current study, we further analyze the same data sets with newly defined Bayesian joint models considering the generalized error distribution for the longitudinal measure instead of normal distribution.

2. Methodology

2.1. Description of Data

The study considers two data sets that are collected from two hospitals under similar settings and as considered by [5] [6] [7] . The data are extracted from patients’ charts which contain epidemiological, laboratory and clinical information of the patients. Patients with ages less than 16 years old are also those who started ART before the defined study period were not included in this study.

Data 1: The first data set is obtained for 354 random samples of HIV/AIDS patients who had been under ART follow-up from January 2006 to December 2012 at the Shashemene referral hospital. There are two outcome variables in each data set. The longitudinal measure is the square root of the number of CD4 cell counts per mm^{3} of blood repeatedly measured at approximately every 6 months interval. The survival outcome is the time in months of a patient to the associated death event. Explanatory variables for the longitudinal response are: visit time, square of visit time, sex, functional status, alcohol use, tobacco use, number of opportunistic infections. Predictors for the survival time are: TB infection status at baseline, awareness about ART, condom use, number of opportunistic infections, number of living rooms at home. Data 2: The second data set is obtained for 400 random samples of HIV/AIDS patients who had been under ART follow-up from January 2008 to March 2015 at the Bale Robe general hospital. Outcome variables are: longitudinal measure which is the square root of the number of CD4 cell counts per mm^{3} of blood that repeatedly measured at approximately every 6 months interval, and survival time in months of a patient to death event. Explanatory variables for the longitudinal measure are: visit time, square of visit time, sex, age, weight, number of opportunistic infections, and those for the survival time are: age, weight, functional status, tobacco, condom. Description and codes of the explanatory variables are described in Table 1.

2.2. Linear Mixed Models

The longitudinal data, CD4 T-cell counts, are measurements on the response variable taken from same individuals over several observation times. Thus the set of observations on a subject tends to be inter-correlated [8] [9] . The two sources of variations expected are the within-patient and the between-patients variations.

Table 1. Explanatory variables with codes.

Analysis of within-patient variation allows studying of changes of the CD4 counts over time, while analysis of between-patients variation allows understanding differences between patients.

Here we assume that the longitudinal CD4 measure has the generalized error distribution for instead of normal. For any variable Y that follows the generalized error distribution, its density function $\text{GED}\left(y;\gamma ,\mu ,{\sigma}^{2}\right)$ with three parameters as adapted by [4] from [10] is:

$f\left(y\right)=a\left(\gamma \right)\frac{1}{\sigma}\mathrm{exp}\left[-c\left(\gamma \right){\left|\frac{y-\mu}{\sigma}\right|}^{\frac{2}{1+\gamma}}\right]$ (1)

with: $a\left(\gamma \right)=\frac{\Gamma {\lceil 3\left(1+\gamma \right)/2\rceil}^{\frac{1}{2}}}{\left(1+\gamma \right)\Gamma {\lceil \left(1+\gamma \right)/2\rceil}^{\frac{3}{2}}}$, $c\left(\gamma \right)={\{\frac{\text{\Gamma}\lceil 3\left(1+\gamma \right)/2\rceil}{\text{\Gamma}\lceil \left(1+\gamma \right)/2\rceil}\}}^{\frac{1}{1+\gamma}}$

Here µ is location and ${\sigma}^{2}$ is scale parameters. And $\gamma $ is the shape parameter of GED that is related to kurtosis of the distribution and characterizes non-normality of Y. The GED can model the error distribution more flexibly than the normal one [4] [10] [11] .

The generalized error distribution generalizes the normal distribution. Normal distribution is a special case of GED when $\gamma =0$ in which case $w\left(0\right)=1/\sqrt{2\text{\pi}}$ and $c\left(0\right)=1/2$ . Other special cases include a Laplace (double exponential) distribution when $\gamma =1$ and when becomes a uniform distribution $\gamma $ approaches-1. It becomes leptokurtic distribution for $\gamma \in \left(0,1\right)$, gets fatter tails when $\gamma >0$, and gets thinner tails than the normal distribution when $\gamma <0$ . Choy and Smith [12] derived the GED as a scale mixture of normal distributions for $\gamma \in \left(0.5,1\right]$ .

The GED in Equation (1) is expressed in a simpler form by [13] as follows:

$f\left(y\right)=\frac{\lambda s}{2\Gamma \left(\frac{1}{s}\right)}\mathrm{exp}\left(-{\lambda}^{s}{\left|y-\mu \right|}^{s}\right)$ (2)

The normal distribution is a special case of this form of GED when s = 2 and so $\lambda =1/\sigma \sqrt{2}$ . But in many situations, data are assumed normal though normality may not be an appropriate assumption. The statistical package such as fitdistrplus developed by [14] can be used to see whether or not measurement errors of data at hand are normally distributed.

The generalized error distribution was first introduced by Subbotin [15] as class of symmetric distributions with variation in kurtosis. The distributions have many structural properties close to a normal distribution. Many researchers have studied GED including its applications but not in the joint models studied here. Nelson [16] developed linear regression and time series models with heavy tails assuming the underlying distribution to be the GED. It can be used in statistical modelling if the observation errors are not necessarily normally distributed. Zhang [4] proposed and studied linear mixed growth models for longitudinal data with the GED so as to handle leptokurtic and platykurtic errors. The author reported that such models fit better to data than the respective models with normality assumptions.

In our case, we first analyzed the CD4 counts data in fitdistrplus package [14] and found that measurement errors of the longitudinal CD4 data seem non-normally distributed. Then we define the linear mixed model using generalized error distribution. Let
${Y}_{it}$ be the longitudinal CD4 measurement of the i^{th} patient
$i=1,2,3,\cdots ,n$ at times
$t={t}_{1},\cdots ,{t}_{T}$ . The linear mixed model for the longitudinal process with assumption of generalized error distribution for the error term is defined as:

${Y}_{it}={\mu}_{i}\left({t}_{it}\right)+{W}_{1i}\left({t}_{it}\right)+{\in}_{it}$ (3)

where ${\mu}_{i}\left({t}_{it}\right)={X}_{it}{\beta}_{i}$ is time dependent mean response of ${Y}_{it}$ as a function of predictors ${X}_{it}$ and coefficients ${\beta}_{i}$, ${W}_{1i}\left({t}_{it}\right)$ is time dependent subject specific random effects having normal distribution with mean zero and variance ${\sigma}^{2}$, and ${\u03f5}_{it}$ is a random error distributed as ${\u03f5}_{it}~\text{GED}\left(\gamma ,0,{\sigma}_{\in}^{2}\right)$ . The shape parameter $\gamma $ is an important parameter to be studied here.

2.3. Survival Models

The survival time is random variable defined on non-negative real numbers. The observed time is taken as the minimum ${T}_{i}=\text{min}\left({t}_{i},{c}_{i}\right)$ of the time to event ${t}_{i}$ and time to censoring ${c}_{i}$ . The time variable is modeled with two AFT distributions (lognormal and loglogistic) as considered in [5] [17] [18] [19] .

We assume that the survival time follows lognormal distribution. Its probability density function $f\left(t\right)$, survival function $S\left(t\right)$ and hazard function $h\left(t\right)$ with parameters $\mu $ and $\sigma $ can be expressed respectively as:

$\begin{array}{l}f\left(t\right)=\frac{1}{\sigma t\sqrt{2\text{\pi}}}\mathrm{exp}\left(-\frac{{\left(\mathrm{log}\left(t\right)-\mu \right)}^{2}}{2{\sigma}^{2}}\right),\\ S\left(t\right)=1-F\left(\frac{\mathrm{log}\left(t\right)-\mu}{\sigma}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}h\left(t\right)=\frac{f\left(t\right)}{S(t)}\end{array}$

The regression model linked with the covariates for each individual i is given as:

${T}_{i}~\text{lognormal}\left({\mu}_{i}\left(t\right),{\sigma}^{2}\right)$ with $\mathrm{log}\left({\mu}_{i}\left(t\right)\right)={X}_{2i}^{\text{T}}{\alpha}_{2i}+{W}_{2i}\left(t\right)$ (4)

Assuming that the survival time follows loglogistic distribution, its probability density function $f\left(t\right)$, survival function $S\left(t\right)$ and hazard function $h\left(t\right)$ with parameters with parameters $\lambda $ and $\rho $ can be expressed respectively as:

$f\left(t\right)=\frac{\lambda \rho {t}^{\rho -1}}{{\left(1+\lambda {t}^{\rho}\right)}^{2}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}S\left(t\right)=\frac{1}{1+\lambda {t}^{\rho}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}h\left(t\right)=\frac{\lambda \rho {t}^{\rho -1}}{1+\lambda {t}^{\rho}}$

The regression model is linked with the covariates for each individual

${T}_{i}~\text{Loglogistic}\left(\rho ,{\lambda}_{i}\left(t\right)\right)$ with $\mathrm{log}\left({\lambda}_{i}\left(t\right)\right)=-\rho \left({X}_{2i}{\alpha}_{i}+{W}_{2i}\left(t\right)\right)$ (5)

The AFT models allow the direct effects of covariates on survival time instead of hazard rate. Given a vector of predictors ${X}_{2i}$, the log-linear form of the AFT model for survival time ${T}_{i}$ of individual patient $i=1,2,3,\cdots ,n$ can be written as:

$\mathrm{log}\left({T}_{i}|{W}_{2i}\right)={X}_{2i}{\alpha}_{i}+{W}_{2i}\left(t\right)+{\epsilon}_{i}$ (6)

where ${\alpha}_{i}$ is a vector of unknown coefficients of ${X}_{2i}$, ${W}_{2i}\left(t\right)$ refers to subject specific random effects having normal distribution, ${\epsilon}_{i}$ is a sequence of mutually independent measurement errors that follows AFT distributions, in this case, lognormal and loglogistic distributions.

3. Bayesian Joint Models

3.1. Likelihood Model

The association between the longitudinal and survival processes is assumed to come through stochastic dependences denoted by ${W}_{1i}\left(t\right)$ and ${W}_{2i}\left(t\right)$ . There are many ways of making the linkages [2] . Here we consider the links used in [5] . Thus the joint models that link the GED based model of longitudinal process in Equation (3) to the AFT based model of survival process in Equations (4)-(6) is given as follows:

${W}_{1i}\left(t\right)={U}_{1i}+{U}_{2i}\ast t$ (7)

${W}_{2i}\left(t\right)={r}_{1}{U}_{1i}+{r}_{2}{U}_{2i}$ (8)

where ${r}_{1},{r}_{2}$ are association parameters. Note that ${U}_{1i},{U}_{2i}$ are latent variables that are independent subject-specific random effects having bivariate Gaussian distribution with mean zeros and constant variances. These effects are assumed to be induced by the longitudinal process to the time-to-event process through the random intercept and random slope terms in the linear mixed model.

We assume Y and T are conditionally independent given the random effects $w=\left({w}_{1},{w}_{2}\right)$ and model parameters $\theta =\left({\theta}_{1},{\theta}_{2}\right)$ . The two sets of parameters are one for the linear mixed model ${\theta}_{1}=\left(\beta ,\gamma ,{\sigma}_{\u03f5}^{2},{\sigma}_{u1}^{-2},{\sigma}_{u2}^{-2}\right)$ and those for the survival model ${\theta}_{2}=\left(\alpha ,{\sigma}_{LN}^{-2},{r}_{1},{r}_{2}\right)$ in the lognormal or ${\theta}_{2}=\left(\alpha ,\rho ,{r}_{1},{r}_{2}\right)$ in the loglogistic case. The joint likelihood function of the data from the two processes can be given as:

$\begin{array}{l}f\left(y,t,\delta |{\theta}_{1},{\theta}_{2},{w}_{1},{w}_{2}\right)\\ ={\displaystyle \int f\left(y|{\theta}_{1},{w}_{1}\right)f\left(t,\delta |{\theta}_{2},{w}_{2}\right)f\left({w}_{2}|{w}_{1}\right)f\left({w}_{1}\right)\text{d}{w}_{2}\text{d}{w}_{1}}\end{array}$ (9)

where each $\delta $ is an indicator for a patient’s survival with $\delta =1$ if death event occurs and $\delta =0$ if censored.

3.2. Prior Distributions

Non-informative joint prior distribution of the parameters $\pi \left(\theta ,{w}_{1},{w}_{2}\right)$ is considered. Individual parameters β’s and α’s are assumed to be independently and identically normally distributed with mean zero and large variance 1000. The association parameters ${r}_{1},{r}_{2}$ are each assumed to have normal distribution with mean zero and variance 1000. The shape parameter $\gamma $ of GED, the shape parameter of loglogistic distribution $\rho $, and precision parameters $\tau ={\sigma}^{-2}$ all are assumed to follow Gamma(2, 0.5).

3.3. Posterior Distribution

The Bayesian model [10] is defined by the posterior distribution $\pi \left(\theta ,w|y,t,\delta \right)$ of the model parameters $\theta $ and random effects $w$ given the data and is expressed by:

$\pi \left(\theta ,w|y,t,\delta \right)=\frac{f\left(y,t,\delta |\theta ,w\right)\pi \left(\theta ,w\right)}{{\displaystyle \int f\left(y,t,\delta |\theta ,w\right)\pi \left(\theta ,w\right)\text{d}\theta \text{d}w}}$ (10)

where $f\left(y,t,\delta |\theta ,w\right)$ is the likelihood function, $\pi \left(\theta ,w\right)$ is the prior probability distribution, and $\int f\left(y,t,\delta |\theta ,w\right)\pi \left(\theta ,w\right)\text{d}\theta \text{d}w$ is the normalizing constant. The main challenge here is computational difficulty. The standard maximum likelihood method involves integrating out latent variables from the log likelihood function which is difficult when the parameters are of high dimensional. Simulation simplifies the computational challenges. Here the Bayesian model in Equation (10) is computed using Markov chain Monte Carlo methods with the Gibbs sampler algorithm that is based on full conditional distributions of the parameters [7] [13] [20] . The Gibbs sampler algorithm is implemented in the WinBUGS software version 1.4 [20] . The final inferences are made based on independent samples taken from the posterior distribution after convergence of the realizations. Time series plots, auto-correlations and Gelman-Rubin statistics are used to assess and confirm convergences.

4. Results and Discussion

The objective of this study is to model the longitudinal CD4 measurement and the associated time to death data using Bayesian joint modelling approach. The generalized error distribution is assumed for the square root of the CD4 T-cell counts, while lognormal and loglogistic distributions are assumed for the survival time. Two data sets collected from two hospitals are analyzed using four Bayesian joint models. The findings from the models are all interpreted as they are important in many ways.

4.1. Descriptive Analysis

For Data 1 taken from Shashemene referral hospital, the average baseline CD4 cell counts is estimated to be 156.9 with standard deviation of 92.5 per mm^{3} of blood sample. By the end of study period, percentage of death event is about 5.9%. The average survival time of the patients is estimated to be 48.7 with standard deviation of 21.3 in months.

For Data 2 taken from Bale Robe general hospital, the average number of baseline CD4 counts is about 177.6 with standard deviation of 104.8 per mm^{3} of blood sample. Among the sample of patients considered, percentage of death event is about 11.5%. The average survival time is about 55.1 with standard deviation of 21.8 in months.

The baseline CD4 counts reveal same variabilities in the two studies which is about 59% as measured by coefficient of variation. However, it seems that there a slight difference between variabilities for time to event data with 44% for Data 1 and 40% for Data 2.

To understand the relationship between the longitudinal measure and follow-up time, mean structures are plotted in Figure 1. The plots show that the average of square root of CD4 counts may have a quadratic relationship with patient’s follow-up time. We thus include both observation time and its square in the linear mixed models as predictor variables.

4.2. Inferential Analysis in the Case of Data 1

In the analysis of Data 1, twenty one parameters are estimated using the two defined Bayesian joint models based on GED-lognormal and GED-loglogistic distributions. The results of analysis are displayed in Table 2 & Table 3.

Table 2. Parameter estimations for the Bayesian Joint GED Lognormal Model in the case of Data 1.

*Significant at 5% significant level; SD: Standard deviation; CI: Credible interval.

Table 3. Parameter estimations for the Bayesian Joint GED Loglogistic Model in the case of Data 1.

*Significant at 5% significant level; SD: Standard deviation; MC: Monte Carlo; CI: Credible interval.

(a)(b)

Figure 1. Plots of mean of square roots of CD4 Counts over observed time for the two studies: (a) Data 1 and (b) Data 2.

4.2.1. Bayesian GED Lognormal Analysis

The results of analysis are displayed in Table 2. They reveal that the shape parameter of the GED of the longitudinal measure is significantly different from zero $\left(\stackrel{^}{\gamma}=1.021,\text{CI}\left(0.8291,1.227\right)\right)$ at 5% level of significance. This is because the 95% credible interval does not contain zero. The standard deviation of the GED is as well significant $\left({\stackrel{^}{\sigma}}_{\in}=3.051,\left(2.911,3.206\right)\right)$ .

The subject-specific random effects ${U}_{1}$ and ${U}_{2}$ are found to be significant since their respective precision parameters are significant $\left({\stackrel{^}{\tau}}_{u1}=0.097,\left(0.079,0.117\right)\right)$ and $\left({\stackrel{^}{\tau}}_{u2}=2.363,\left(1.796,3.133\right)\right)$ . For AFT model, the precision parameter in the lognormal distribution is significant $\left({\stackrel{^}{\tau}}_{LN}=7.869,\left(6.743,9.081\right)\right)$ . However, the association parameters ${r}_{1}$ and ${r}_{2}$ in the joint model are not significant $\left({\stackrel{^}{r}}_{1}=0.001\left(-0.003,0.004\right)\right)$ and $\left({\stackrel{^}{r}}_{2}=0.003\left(-0.018,0.025\right)\right)$ . Thus the Bayesian joint GED lognormal model is not significant.

4.2.2. Bayesian GED Loglogistic Analysis

Analysis results of the Bayesian GED loglogistic model in the case of Data 1 are displayed in Table 3. The parameters of the generalized error distribution are both significant at significant 5% level based on estimations of its shape parameter $\left(\stackrel{^}{\gamma}=1.017,\left(0.823,1.229\right)\right)$ and scale parameter $\left({\stackrel{^}{\sigma}}_{\in}=3.047,\left(2.910,3.198\right)\right)$ . The random effects ${U}_{1}$ and ${U}_{2}$ are significant as well since their respective precision parameters are significant $\left({\stackrel{^}{\tau}}_{u1}=0.097,\left(0.0795,0.117\right)\right)$ and $\left({\stackrel{^}{\tau}}_{u2}=2.326,\left(1.775,3.016\right)\right)$ . In the survival sub-model, the shape parameter of the loglogistic distribution is significant $\left(\stackrel{^}{\rho}=6.450,\left(5.878,7.034\right)\right)$ . Similar to the lognormal case, the Bayesian joint GED loglogistic model is not significant since the association parameters ${r}_{1}$ and ${r}_{2}$ are both insignificant $\left({\stackrel{^}{r}}_{1}=0.002,\left(-0.007,0.012\right)\right)$ and $\left({\stackrel{^}{r}}_{2}=0.001,\left(-0.055,0.058\right)\right)$ .

Comparing their total DIC values, the Bayesian GED loglogistic model fits better to Data 1 than the lognormal case. From the analysis of Bayesian GED loglogistic model, the covariates that are found to affect the longitudinal CD4 measure are: observed time, squared observed time and gender. They imply that the disease marker improves over time but later reaches maximum and then declines. Female patients gain more CD4 counts as compared to the males. Effects of the functional status, alcohol use, tobacco use, and opportunistic infection status of patients are not statistically significant. From the survival sub-model, the results show that TB infection status at baseline, awareness about ART and condom use have significant effects on the survival time of a patient.

4.3. Inferential Analysis in the Case of Data 2

Here we analyze Data 2 and estimate twenty parameters of the same two defined models: Bayesian joint GED-lognormal and Bayesian GED-loglogistic models. The results are given in Table 4 & Table 5.

4.3.1. Bayesian GED Lognormal Analysis

The results of analysis in Table 4 show that the shape parameter of the generalized growth distribution is significantly different from zero $\left(\stackrel{^}{\gamma}=0.05244,\left(0.00209,0.1442\right)\right)$ . The standard deviation of the GED is significant $\left({\stackrel{^}{\sigma}}_{\in}=5.006,\left(4.846,5.176\right)\right)$ and relatively larger than those obtained from Data 1 case. The subject-specific random effects ${U}_{1}$ and ${U}_{2}$ are significant due to the fact that their respective precision parameters are significant $\left({\stackrel{^}{\tau}}_{u1}=4.764,\left(0.547,15.47\right)\right)$ and $\left({\stackrel{^}{\tau}}_{u2}=5.231,\left(0.866,12.12\right)\right)$ .

For the survival sub-model, the precision of the lognormal distribution is significant $\left({\stackrel{^}{\tau}}_{LN}=7.43,\left(5.678,10.53\right)\right)$ . Again, the association parameters ${r}_{1}$ and ${r}_{2}$ are not significant ${\stackrel{^}{r}}_{1}=\left(-0.0047,\left(-0.1297,0.101\right)\right)$ and $\left({\stackrel{^}{r}}_{2}=-0.010,\left(-0.144,0.178\right)\right)$ and hence the Bayesian joint GED lognormal model is not significant in the case of Data 2.

The covariates observed time, square of observed time, sex, age, weight and number of opportunistic infection have statistically significant effects on CD4 counts of the HIV/AIDS patients. For survival part, age, functional status, tobacco use, and condom use have effects on survival time of the patients. Patient’s weight is not significant in this case.

4.3.2. Bayesian GED Loglogistic Analysis

The results are displayed in Table 5. The shape parameter of the generalized

Table 4. Parameter estimations for the Bayesian Joint GED Lognormal Model in the case of Data 2.

*Significant at 5% significant level; SD: Standard deviation; MC: Monte Carlo; CI: Credible interval.

Table 5. Parameter estimations for the Bayesian Joint GED Loglogistic Model in the case of Data 2.

*Significant at 5% significant level; SD: Standard deviation; MC: Monte Carlo; CI: Credible interval.

growth distribution is significant $\left(\stackrel{^}{\gamma}=0.093,\left(0.089,0.190\right)\right)$ and its standard deviation is also significant $\left({\stackrel{^}{\sigma}}_{\in}=4.778,\left(4.778,4.95\right)\right)$ . The subject-specific random effects ${U}_{1}$ and ${U}_{2}$ are significant because their respective precision parameters are so $\left({\stackrel{^}{\tau}}_{u1}=4.705,\left(3.939,15.05\right)\right)$ and $\left({\stackrel{^}{\tau}}_{u2}=4.605,\left(4.356,10.15\right)\right)$ .

For the survival sub-model, the precision parameter of the loglogistic distribution is significant $\left(\stackrel{^}{\rho}=5.559,\left(5.526,6.539\right)\right)$ . Only the slope term ${r}_{2}$ of the association parameters is significant $\left({\stackrel{^}{r}}_{2}=0.627,\left(0.643,1.005\right)\right)$ but not the intercept $\left({\stackrel{^}{r}}_{1}=-0.021,\left(-0.005,0.238\right)\right)$ . The Bayesian joint GED loglogistic model is significant implying that the joint model is important for this data set.

The Bayesian joint GED loglogistic model fits better to the data than the lognormal case based on the total DIC values. From analysis of Data 2 with the selected loglogistic model, the covariates that are found to significantly affect longitudinal CD4 measure of an HIV/AIDS patient are: observed time, square of observed time, sex, age, weight and number of opportunistic infection. With increasing ART follow-up time, CD4 counts of a patient can get parabolic mean growth with reaching maximal point. Female patients achieve higher CD4 counts on average than males. Moreover, the mean CD4 counts of a patient declines as patient’s age increases, weight decreases and number of opportunistic infection increases.

For survival sub-model, the results reveal that improving a patient’s weight improves her/his survival time. Note that weight is an important variable in explaining both the longitudinal and survival processes of a patient. Healthy functional status and condom use during sexual intercourse have also positive effects on survival time. Age and tobacco use are not significant for this data case.

4.4. Assessment of Convergence

For each of the Bayesian models, three parallel sampling chains of 60000 iterations with different starting values are generated. Some plots are given in Figure 2. Inferences are made based on samples of the posterior distributions that are taken with thinning of 10 after burn-in of 25000. Time series plots of the history of the simulations show a reasonable degree of randomness and they may convergence to same values. Auto-correlations and Gelman-Rubin statistics are also used to assess convergences. Finally independent samples are taken from the

(a)(b)(c)

Figure 2. Plots from analysis of Data 1 using the Bayesian GED Loglogistic. (a) Time Series plots of simulations of shape parameter γ of GED and shape parameter ρ of loglogistic distribution; (b) Autocorrelation plots; (c) Plots of Gelman-Rubin statistics.

posterior distribution after convergence of the realization with specified burn-in and thinning values, and then all inferences are made using those samples.

Assessment plots are displayed in Figure 2 from the analysis of Data 1 using the Bayesian GED loglogistic model. Results of two parameters: shape parameter $\gamma $ of GED and shape parameter $\rho $ of loglogistic distribution are illustrated. They show that the simulations converge.

The findings reveal that the generalized error distribution for both data sets has positive estimate of the shape parameter and so it is of fatter tails than normal distribution. There is higher variation on the CD4 T-cell counts for the data from Bale Robe general hospital (about 49%) than that obtained from Shashemene referral hospital (about 10%).

The posterior distributions estimated under the selected models, Bayesian GED loglogistic models, are the solutions required in this analysis. Though the association parameter in the joint model is significant for one data but not for the other data case, both fitted models are still important to consider as they are newly defined models implementing the generalized error distribution. Thus the respective results are used to report findings on how the longitudinal CD4 counts and survival time of a patient are related and parameter estimations.

The findings from this study are fairly consistent with the studies by [5] except for the type of models selected. They suggested different models for the two data sets while only one is recommend here. This may be expected as the GED is involved in this study and as indicated by [4] the GED model can gain more insight on the error distributions than the normal growth curve models.

5. Conclusions

The current study focuses on developing Bayesian joint models with the assumption of generalized error distribution for the longitudinal CD4 observations and of two AFT distributions for the survival time of HIV/AIDS patients. Analyses of two different data sets show that measurement errors of the longitudinal CD4 variable are not normally distributed and so are modelled by the generalized error distribution. The distributions have fatter tails than the normal distribution. The Bayesian joint GED loglogistic models are found to be important models in fitting to the data sets. Fairly consistent estimates of parameters and more insights of the data are obtained from the models. In one of the data sets, it is found that survival time of a patient is affected by the latent variables generated from the longitudinal CD4 T-cell counts.

Covariates with significant effects are identified from analysis of the Bayesian GED loglogistic models. The findings reveal that under ART follow-up, patients’ health can be improved over time. Female patients gain more CD4 counts as compared to the males. Survival time of a patient is negatively affected by TB infection. Awareness about ART is an important factor as well.

Increase in number of opportunistic infection implies decline of CD4 counts. Age of a patient negatively affects the disease marker but no effects on the survival time. Weight loss is related with decline of CD4 counts and shortening of survival time on average. Improving a patient’s weight may improve her/his survival time. Condom use is positively related with survival time of a patient.

Bayesian joint models with GED and AFT distributions are found to be useful in modelling the longitudinal and survival processes. Bayesian computations of such complex models are well handled in the WinBUGS software by providing simulations and estimations of the posterior distributions. We recommend Bayesian joint models with generalized error distributions when measurement errors of the longitudinal data are not necessarily normally distributed. Further studies may investigate the models with various types of shared random effects, more covariates and prediction of future observations.

Acknowledgements

The authors would like to sincerely thank the Shashemene Referral and Bale Robe General Hospitals for providing us the data sets used in this study. We acknowledge the anonymous reviewers for their detailed comments and suggestions.

References

[1] Sousa, I. (2011) A Review on Joint Modelling of Longitudinal Measurements and Time-to-Event. Statistical Journal, 9, 57-81.

[2] Henderson R., Diggle, P. and Dobson, A. (2000) Joint Modeling of Longitudinal Measurements and Event Time Data. Biostatistics, 4, 465-480.

https://doi.org/10.1093/biostatistics/1.4.465

[3] Hatfield, A., Hodges, S. and Carlin, P. (2012) Combining Longitudinal and Survival Information in Bayesian Joint Models: When Are Treatment Estimates Improved? Biostatistics, 1-31.

[4] Zhang, Z.Y. (2013) Bayesian Growth Curve Models with the Generalized Error Distribution. Journal of Applied Statistics, 40, 1779-1795.

https://doi.org/10.1080/02664763.2013.796348

[5] Erango, M.A., Goshu, A.T., Buta, G.B. and Dessiso, A.H. (2017) Bayesian Joint Modeling of Survival of HIV/AIDS Patients Using AFT Data and Longitudinal CD4 Cell Counts. British Journal of Medicine & Medical Research, 20, 1-12.

https://doi.org/10.9734/BJMMR/2017/32123

[6] Buta, G.B., Goshu, A.T. and Worku, H.M. (2014) Bayesian Joint Modelling of Disease Progression Marker and Time to Death Event of HIV/AIDS Patients under ART Follow-Up. British Journal of Medicine & Medical Research, 5, 1034-1043.

https://doi.org/10.9734/BJMMR/2015/12907

[7] Dessiso, A.H. and Goshu, A.T. (2017) Bayesian Joint Modelling of Longitudinal and Survival Data of HIV/AIDS Patients: A Case Study at Bale Robe General Hospital, Ethiopia. American Journal of Theoretical and Applied Statistics, 6, 182-190.

https://doi.org/10.11648/j.ajtas.20170604.13

[8] Wu, L. (2009) Mixed Effects Models for Complex Data. CRC Press, Boca Raton.

https://doi.org/10.1201/9781420074086

[9] Viviani, S. (2012) Mixed Effect Joint Models for Longitudinal Responses with Dropout: Estimation and Sensitivity Issues. Thesis, Department of Statistics, Sapienza—University of Rome, Rome, 1-83.

[10] Box, G.E.P. and Tiao, G.C. (1973) Bayesian Inference in Statistical Analysis. Wiley, New York.

[11] Nadarajah, S. (2005) A Generalized Normal Distribution. Journal of Applied Statistics, 32, 685-694.

https://doi.org/10.1080/02664760500079464

[12] Choy, S.T.B. and Smith, A.F.M. (1997) On Robust Analysis of a Normal Location Parameter. Journal of the Royal Statistical Society: Series B, 59, 463-474.

https://doi.org/10.1111/1467-9868.00079

[13] Purczynski, J. (2003) The Use of Computer Simulations in the Estimation of Selected Econometric and Statistical Models. Publishing Series Rozprawy I Studia (DXXV) 451, University of Szczecin, Szczecin.

[14] Delignette-Muller, M.L. and Dutang, C. (2015) Fitdistrplus: An R Package for Fitting Distributions. Journal of Statistical Software, 64, 1-34.

https://doi.org/10.18637/jss.v064.i04

[15] Subbotin, M.T. (1923) On the Law of Frequency of Error. Matematicheski Sbornik, 31, 296-301.

[16] Nelson, D.B. (1991) Conditional Heteroskedasticity in Asset Returns: A New Approach. Econometrical, 59, 347-370.

https://doi.org/10.2307/2938260

[17] Klein, J.P. and Moeschberger, M.L. (2003) Survival Analysis Techniques for Censored and Truncated Data. Second Edition, Statistics for Biology and Health, Springer, Berlin.

[18] Guo, X. and Carlin, B.P. (2004) Separate and Joint Modeling of Longitudinal and Event Time Data Using Standard Computer Packages. The American Statistician, 58, 16-24.

https://doi.org/10.1198/0003130042854

[19] Tseng, Y., Hsieh, F. and Wang, J.L. (2005) Joint Modelling of Accelerated Failure Time and Longitudinal Data. Biometrika, 92, 587-603.

https://doi.org/10.1093/biomet/92.3.587

[20] Spiegelhalter, D., Thomas, A., Best, N. and Lunn, D. (2003) WinBUGS Version 1.4 User Manual.

http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/manual14.pdf