Joint Modeling of Longitudinal CD4 Count and Weight Measurements of HIV/Tuberculosis Co-infected Patients at

As HIV/TB co-infected patients are started to be visited, it is common to measure weight and CD4 repeatedly overtime to determine the health status of patients. Most of the time linear mixed modeling of weight and CD4 count cannot handle the association between the outcomes whereas the joint modeling of multivariate linear mixed model does. Thus, this study was an attempt to model jointly the longitudinal CD4 and weight measurements of HIV/TB co-infected patients. This retrospective study consists of 254 HIV/TB co-infected patients who were 18years old and above, and on ART followup from 1st February 2009 to 1st July 2014 at Jimma University Specialized Hospital. Firstly, weight and square root of CD4 count were analyzed separately. Based on the separate model, the joint models were built to know the correlation between mean change of weight and CD4 count overtime. Finally, appropriate model was selected based on AIC and BIC values. The ﬁt statistics showed that the joint model ﬁtted the data better than the separate model. From the joint model sex, educational level and functional status were the factors contributing to the prediction of HIV/TB co-infected patients weight at baseline. Beside the linear time effect has a positive effect on the mean change of weight whereas the quadratic time change has negative effect. The baseline CD4 count was differ by patient status and functional status. Further, the linear time effect has a positive sign and found to be statistically signiﬁcant at 5% level of signiﬁcance on the mean change of the square root of CD4 count. Nevertheless, the quadratic time effect has a signiﬁcant negative effect. The ﬁnding of the current study revealed that there was a moderate positive association between the mean change of weight and square root of CD4 count overtime.

finding of the current study revealed that there was a moderate positive association between the mean change of weight and square root of CD4 count overtime.
Keywords Co-infection · Longitudinal analysis · Linear mixed model · Multivariate linear mixed model

Background of the Study
Tuberculosis and HIV have been closely linked since the emergence of AIDS and TB is the most common infectious disease affecting HIV-sero positive individuals and causing to their death [1,2]. Globally, the number of TB patients who had been diagnosed with HIV status reached 2.1 million in 2010 which is equivalent to 34 % of notified cases of TB. Of the 8.8 million incident cases worldwide an estimated 1.1 million (13 %) were found to be co-infected with HIV [3]. Today, HIV and TB treatments are common in many societies and the use of drugs has altered the joint dynamics of both diseases. About one third of 39.5 million HIV infected people worldwide were co-infected with TB [4] and up to 50 % of individuals living with HIV are expected to develop TB. Many TB carriers who were infected with HIV are 30-50 times more likely to develop active TB than those without HIV. For individuals infected with HIV the presence of other infections including TB tends to increase the rate of HIV replication. This acceleration may result in higher levels of infection and rapid HIV progression to the AIDS stage [5]. Over the past 20 years, HIV has fuelled TB notification rates which have increased 3-5 fold in many African countries. By 2007, the continent accounted for 79 % of the global burden of HIV associated TB [6]. Worst affected were those countries in the east and south of the continent where HIV prevalence rates are highest. In South Africa and Swaziland approximately 1 % of the population develops TB annually. Notification rates in some poor communities in South Africa have even increased to over 2 % per year rates that are almost unprecedented in the era of short course multi-drug chemotherapy [7,8]. WHO ranked Ethiopia as seventh among the 22 high burden countries with TB where the estimated annual incidence and prevalence, respectively, were 379 and 643 cases per 100,000 populations [9]. The prevalence of HIV among TB patients was estimated as high as 41 % [10,11]. In Ethiopia routine data obtained from 44 locations in the year 2005/2006 showed that about 41 % of TB patients were HIV+. In addition, another routine data set collected in 2006/2007 estimated that some 31 % of TB patients were HIV+ [12]. TB was the cause of 76,000 deaths in Ethiopia out of which 30 % were among HIV+ patients [6]. Beside the rate of mortality caused by the co-infection, the extent of the negative impact on the quality of livelihood resulting from mental disorders was studied [13]. In many medical cases more than one clinical outcome are measured longitudinally at the same time on the same subject where these measured clinical outcomes are correlated. For example, SBP and DBP for hypertensive patients and CD4 count, beta-2-macroglobulin and weight for HIV infected patients were measured longitudinally at same time. Since they are highly related changes in either often affect changes in the other. In such cases the univariate longitudinal analysis does not take into account correlation between observations on different response variables at each time points. Beside this knowing how the evolution of one is related to the evolution of the other, as well as how the association changes or evolves overtime is difficult from univariate longitudinal analysis. Joint modeling of longitudinal data in other way round accounts two types of correlations which are known to be serial correlation and cross correlation. The serial correlation accounts correlation between observations at different time points within a subject and the cross correlation accounts correlation between different responses at each time point. If different types of outcomes are measured at each time point, the correlation structure is more complicated and hence, more difficult for drawing inference [14,15]. Therefore, the primary objective of present study was to model jointly weight which act as a marker for TB patients and CD4 count which act as a biomarker for HIV infected individuals. Furthermore, the study deals with linear mixed modeling of weight and CD4 count measurements to know the factors that affect the mean change of each outcome variable overtime.

Data Source
The HIV/TB co-infection data used for this study were obtained from Jimma University Specialized Hospital HIV and TB Outpatient Clinics, South West of Ethiopia. The study population consists of all HIV/TB co-infected patients who were 18 years old and above, and who were on ART treatment any time between 1st February 2009 and 1st July 2014. Of the total 850 co-infected patients from both clinics at the hospital, 254 patients who have at least one CD4 count and weight at the same time were considered for the study. All the patients' epidemiological, laboratory and clinical information were obtained from the patients chart of ART followup retrospectively from patients card.

Ethical Consideration
Ethical clearance was obtained from Department of Statistics of Jimma University. Personal information was kept confidentially without disclosing to others during data collection from patient cards.

Variables of the Study
The two outcome variables considered for this study were the longitudinal CD4 count measurement and weight of co-infected patients which is measured in kg. The CD4 cell counts per mm 3 of blood, which is considered as a biomarker was measured approximately within 6 months interval whereas the weight of the patients were mea-  Table 1 with respective of their categories.

Model Specification
Longitudinal responses may arise in two common situations. One is when the measurements taken on the same subject at different times and the other is when the measurements taken on related subjects. In both cases, the responses are likely to be correlated [14]. For longitudinal data, two sources of variations were considered which were the within-subject which arises during the measurements within each subject, and between subject variation which arises during the measurement between different subjects. Modeling within subject variations help us to study changes overtime while modeling between subject variation help us to understand differences between subjects.

Univariate Linear Mixed Effects Model
Before going to the joint modeling of longitudinally measured weight and CD4 count of HIV/TB co-infected patients, the linear mixed model were employed for each outcome to identify an appropriate covariates that predicts the mean change of weight and CD4 count overtime. Therefore, the linear mixed-effects model handling the two of variation is where y j is the n j × 1 vector of observed response values, β j is the p × 1 vector of fixed-effects parameters, X j is the n j × p observed design matrix corresponding to the fixed-effects, b j is the q × 1 vector of random-effects parameters, Z j is the n j × q observed design matrix corresponding to the random-effects, and ε j is the n j × 1 vector of residuals for the jth response. The corresponding assumptions for model (1) are b j ∼ N q j (0, j ) and ε j ∼ N n j (0, j ); where j and j are the variance-covariance matrix for b j and ε j for each outcome variable respectively.

Methods of Parameter Estimation
Suppose that a random sample of n j for jth response observations is obtained from a univariate normal linear mixed effects model as defined in equation (1). Then, the likelihood of the model parameters given the vector of n j observations is defined as where β j is a vector of fixed-effect parameters and γ j is a vector containing the variance parameters for the jth response. Hence, the log-likelihood of the model parameters is defined as The required parameters of the model were obtained by maximizing the loglikelihood function with respect to β j and γ j . However, the maximum likelihood approach may produce variance parameters that are biased downwards since they are based on the assumption that the fixed-effects parameters are known [16]. To handle this problem the present study used restricted maximum likelihood (REML) method, which treat fixed-effects as parameters rather than constant.

Multivariate Linear Mixed Model
The mixed model can be easily extended to include multiple response variables by further stacking the data and defining a specific variance-covariance structure for the random effects. In our case, consider modeling weight and CD4 count measurements of co-infected patients jointly overtime by incorporating random intercepts and slopes in order to model the correlations overtime between the two responses. As we defined earlier, the primary aim of this study was to model jointly the longitudinally measured CD4 count and weight of HIV/TB co-infected patients which was modeled by multivariate linear mixed model.
Let y i jk represent the ith observation from the jth subject for the kth response variable, where i = 1, 2, . . . , n jk ; j = 1, 2, . . . , s; k = 1, 2. Also define N k = s i=1 n jk and N = N 1 + N 2 . The vector y jk = [y 1 jk , y 2 jk , . . . , y n jk jk ] represents the n jk observations of the kth response variable from the jth subject and the vector y k = [y 1k , y 2k , . . . , y sk ] represents the N k observations for the kth response variable across all response variables and subjects. Finally, the vector y = [y 1 , y 2 , . . . , y k ] represents the N observations across all response variables and subjects. In the context of modeling our response variables, the linear mixed-effects models for each response variable for subject j taken at time t can be specified as where μ k (t) refers to the average evolution (of the kth response over time) and is a function of the fixed effects. The subject specific random intercepts a jk and slopes b jk (t) describe how the subject specific profiles deviate from the average profile for the kth response. The two response trajectories are joined together by assuming a joint distribution for the vector of random-effects, b j , such as where the variance-covariance matrix for the random effects, , has the following structure: The error components for each response, which are independent of the random effects can be taken to be correlated or uncorrelated (σ 12 = 0), such that the error components are defined as

Special Cases for the Random Effect Variance Matrix
We obtain special cases for the variance-covariance matrix of the random effects by making specific assumptions for the variance-covariance matrix . The first assumption is when the two outcome variable could be taken to be completely independent at any point in time, we impose their variance-covariance matrix which has the following special form given by Within a response variable, the random intercept and slope induce within-subject correlations in the repeated measures overtime, while assuming independence between subjects. Moreover, this model assumes that the two responses are completely independent. The results for the model would be identical to fitting two separate random-effect models.
The second assumption is when the two response variables could be taken to be completely dependent. In this case, the two responses essentially "share" the same set of random effect parameters (intercept and slope). When two parameters are completely dependent, the correlation between them is equal to one. This occurs when the covariance between the parameters is equal to the square root of the product of their respective variances. Most notation, however, define the model with a 2 × 1 vector of random effects, such as Clearly, the aforementioned structure imposes strong assumptions on the relationship between the two response variables. It is very unlikely that the two responses would exhibit complete dependence in the association between the random slopes and between the random intercepts. One advantage of this model, when the assumption is tenable, is that it drastically reduces the number of random effects that must be estimated when the number of response variables is large. For models with a large number of response variables, estimation would likely be impossible if the shared-parameters (or alternative approach) were not used.

Association of the Evolutions (AOE)
One important question that may be addressed with a joint mixed-effects model is how the evolution of one response is associated with the evolution of another response ("association of the evolutions"). By definition, the correlation between the evolutions for the two random slopes is given by It may be noted that the above expression is produced using those values from the matrix defined in (5).

Evolution of the Association (EOA)
A similar idea that may be investigated using a joint mixed effects model is how the association between the responses evolves overtime ("evolution of the association"). Assuming uncorrelated errors, the marginal correlation between the two responses as a function of time is given by Assuming correlated errors, the marginal correlation between the two responses as a function of time is given by The delta method could be used to obtain 95 % confidence bounds for r m (t) at any particular point in time. Two observations can be made from the uncorrelated errors by noticing t = 0 the marginal correlation reduces to which is essentially the correlation between the two random intercepts. In fact, when the error components are small, the closer the marginal correlation at t = 0 approximates the correlation between the random intercepts. Also, as t increases r m (t) converges to r e for the case with uncorrelated errors, and to r m (t) = σ a 1 a 2 + σ 12 for the case of correlated errors, which indicates that the absolute value of the marginal correlation at t = 0 cannot be higher than the correlation between the random intercepts. It may also be noted that as t increases the marginal correlation converges to the correlation between the random slopes, while the variance-covariance parameters of the random effects determine the shape of the marginal correlation function [17].

Joint Model Estimation Techniques
In the particular context of random-effects models, so-called adaptive quadrature rules can be used [18], were the numerical integration is centered on the estimates of the random effects, and the number of quadrature points is then selected in terms of the desired accuracy. To illustrate the main ideas, we consider Gaussian and adaptive Gaussian quadrature, designed for the approximation of integrals of the form for a known function f (z) and for φ(z) the density of the multivariate standard normal distribution. Therefore, first standardize the random effects such that they get the identity covariance matrix. Then, the likelihood contribution for subject i equals where b i is q × 1 dimensional vector of unknown random effects and b i ∼ N (0, ), β is a vector of fixed-effects parameters and γ is a vector containing the variance parameters. f (z) and for φ(z) denotes the density of the multivariate standard normal distribution.

Model Selection Criteria
Several model selection procedures exists but none of which were the best. To have an appropriate model for the univariate LMM and multivariate LMM most commonly known model selection criterions; Akaike Information Criterion (AIC) [19] and the Bayesian Information Criterion (BIC) [14] were considered for this study.
where −2 log L is twice the negative log-likelihood value for the model, p is the number of estimated parameters, npar denotes the total number of parameters in the model and N is the total number of observations used to fit the model. Smaller values of AI C and B I C reflect an overall better fit.

Descriptive Statistics
Some demographic information and baseline characteristics of all patients disaggregated by patient status were presented in Table 2. Regarding the sex composition of patients, out of total of 254 co-infected patients 139 (54.72 %) of them were males and 18 (56.25 %) death were also occurred in male group in comparison with female group. More than half 146 (57.48 %) of the co-infected patients belongs to orthodox religious group whereas 17 (6.69 %) belongs to protestant religious group. Of the total deaths occurred in these categories small number 2 (6.25 %) of deaths were occurred in protestant religious group. When we look at the educational level category of the co-infected patients, 108 (42.52 %) attended primary education while only 17 (6.69 %) attended tertiary education.   working categories of patients functional status respectively. To this end, 37 (66.07 %) and 12 (21.43 %) were missed to followup from these two groups respectively in comparison with bedridden group at baseline during the time of co-infection periods.
Mean with the corresponding standard deviation at each time points with respective sample sizes for both outcome variables was shown in Table 3. The average number of square root of CD4 count was 11.965 at baseline. There was a general increment in the mean value up to 30 months starting from baseline and starts declining. However, when we look at the standard deviations up to 36 months there was slight variation among the measurement times where smaller variation in square root of CD4 measurement was observed at 48 month. Also it can be observed from Table 3, smaller mean weight of co-infected patients were observed at the baseline of co-infection period and larger mean weight was observed at 48 month with larger weight variation among the co-infected patients.

Exploring Individual Profile and the Mean Structure
To underpin the model building and visualize the pattern of weight and CD4 measurements of the patients overtime, the overall individual plots were considered.
The graph (Fig. 1) demonstrates the variability (within and between patients) in square root of CD4 count and weight measurements of HIV/TB co-infected patients. The red line loess smoothing technique suggests that the mean structure growth of both variables was in a linear and quadratic pattern overtime. Nevertheless, the rate of increment was high in square root of CD4 count than weight.

The Separate Longitudinal Analysis
The separate longitudinal analysis were started from the fixed-effect modeling to select an appropriate covariates that predicts weight and square root of CD4 count where this linear model only considers the between source of variation among co-infected patients. First, we have used stepwise automatic variable selection method to select appropriate fixed-effects from Table 1. After arriving at an optimal subset, the selected  fixed-effects model was fitted with different random effects for both responses to have an appropriate linear mixed-effects models. The summary of the linear mixed-effects models which were modeled by considering different random effects were shown in Table 4. From Table 4, for both responses only considering the quadratic time effect to the model worsen the LMM. However, considering the linear time effects to the quadratic time effect improves the LMM with only quadratic time effects since its inclusion lowers the AIC and BIC values. Finally, when we look at the improvement of the model with inclusion of random intercept to that of random linear and quadratic time effects of the LMM, there was an improvement of the model and this model have lower AIC and BIC values than the remaining six LMMs. Hence, we consider the LMM with random intercept linear and quadratic time effects as an appropriate model for the separate models.    As we can observe from the fitted LMM for weight of co-infected patients (Table 5), there was significance difference weight at baseline by sex, patients status and educational level since the estimated coefficient for male, missed the follow up and secondary education categories were significant ( p < 0.05). But there was no significance difference weight at baseline by the functional status category of co-infected patients.
When we look at the mean change in weight overtime, the linear time have positive significant effects whereas the quadratic time have negative significant effects. And also, the mean change overtime in weight of co-infected patients were differ by the functional status of the patients since the interaction of linear time with working functional status have significant effects on weight at 5 % level of significance. When we look at the fitted LMM for square root of CD4 count, there was a significance difference in square root of CD4 measurement at baseline by the functional status of the patients since the estimated coefficient of bedridden category was significant at 5 % level of significance. Likewise the case of weight, the linear time have positive effect on mean change of square root of CD4 count overtime whereas the quadratic time change have negative effect at 5 % level of significance.

Joint Analysis of Square Root of CD4 Count and Weight
To have an appropriate joint model that represent the longitudinally measured CD4 count and weight of the co-infected patients, different candidate joint model with different random effect for the joint modeling were considered. The AIC and BIC were used as a guideline in selecting covariates for the model. A smaller AIC and BIC values were generally indicates a better model. The results of the fitted models information criterion were displayed in Table 6. Table 6 elucidates that the joint model with random effects was considered as the worst model compared to the joint model with inclusion of random linear and quadratic slopes since it have larger AIC and BIC values. Moreover, inclusion of random linear slope to the random intercept for the joint modeling improves the joint model with random intercept only. However, when we include quadratic slope time effects for the joint modeling to that of with random intercept and linear slope, the joint model becomes worst since it has larger AIC and BIC. Therefore, we consider the joint model with random intercept and linear slope as an appropriate joint model since it have smaller AIC and BIC values than the remaining models. Table 6 also shows that the value of fit statistics (AIC and BIC) for all joint models were less than that of the separate model. This indicates that joint model was fitting better the data as compared to the univariate model. Therefore, our joint analysis of these two outcomes was justified and the separate analysis results were likely to be biased. The result of selected appropriate joint model for the longitudinally measured CD4 count and weight of HIV/TB co-infected patients was shown in Table 7.
As can be seen from the fitted joint model (Table 7), covariates like sex, educational level and functional status were the factors significantly contributing to the prediction of HIV/TB co-infected patients weight at baseline. Table 7 also suggests that the linear time effect has a positive effect on the mean change of weight whereas the quadratic time change has negative effect overtime at 5 % level of significance. To this end, weight change of the patients were differ by functional status since the estimated coefficient for ambulatory functional status group with linear time interaction have significant effect.
From the fitted joint model, the baseline CD4 cell count was differ by patient status and functional status since the estimated coefficients for missed patient status category in comparison with transferred category, ambulatory and bedridden functional status category in comparison with working functional status category have significant effects at 5 % level of significance. With regard to the mean change in square root of CD4 measurements overtime, the linear time effect has a positive sign and found to be statistically significant. Nevertheless, the quadratic time effect has a significant negative effect on the mean change of CD4 count overtime. Further, the mean change in CD4 cell count of the patients differ by functional status category since the estimated coefficient of ambulatory functional status with time interaction in comparison with working functional status category with time interaction have a significant effects on the mean change of CD4 count overtime at 5 % level of significance.
The estimated variance-covariance matrix for the final joint model waŝ  1ρ a 1 b 1ρ a 1 a 2ρ a 1  From the estimated variance-covariance matrix (ˆ ) of the estimated joint model for the random effects, it can be seen that variability was lower for square root of CD4 count (σ 2 a 1 = 78.95) than weight (σ 2 a 2 = 907.20). Beside, the correlation matrix shows that (R) the baseline subject specific baseline CD4 measurement and weight were negatively correlated (ρ a 1 a 2 = −0.841). On the other hand, there was a moderate positive association between subject specific change in square root of CD4 and weight of the patient overtime (ρ b 1 b 2 = 0.541). There was a negative cross correlation between the baseline weight and time slope of square root of CD4 measurement (ρ b 1 a 2 = −0.303) overtime. Similarly, there was a negative correlation between subject specific baseline weight and the patient specific change in square root of CD4 overtime (ρ b 2 a 1 = −0.214) overtime. It can thus be concluded that there was a moderate positive association between the mean change of weight and square root of CD4 counts overtime.

Conclusion
In this study univariate and bivariate linear mixed methods were considered for fitting two continuous response variables measured longitudinally. Based on separate analysis, the evolution of weight was significantly differ with respect to sex (female), baseline age, time, quadratic time, missed patient status category, secondary education category and the interaction of linear time with working functional status. Like wise, the baseline CD4 measurement was differ by patient status and functional status since the estimated coefficients for missed patient status category in comparison with transferred category, ambulatory and bedridden functional status category in comparison with working functional status category have significant effect. The results of this study also demonstrates that the joint modeling of longitudinally CD4 count and weight measurements fits the data better than those obtained from the separate model. To sum up, the joint model suggests that there was a moderate positive association between the mean change of weight and square root of CD4 count overtime.