A Joint Model for a Longitudinal Pulse Rate and Respiratory Rate of Congestive Heart Failure Patients: at Ayder Referral Hospital of Mekelle University, Tigray, Ethiopia

Pulse rate(PR) and respiratory rates(RR) are main symptoms of congestive heart failure(CHF) and the abnormal PR and RR are broad indicators of major physiological instabilities. The lower PR and RR are associated with a strong and healthier heart. CHF is a complex clinical syndrome that can result from any structural or functional cardiac disorder that impairs the ability of the ventricle to fill with or eject blood. The main objective of this study is, to investigate the joint evolution of PR and RR of CHF patients and identify the potential risk factors affecting the two end points. The latest data from Sep.2012 up to Aug.2013 have been taken from medical charts of 264 adult CHF patients to model separate and joint linear mixed effect for PR and RR. The baseline mean and standard deviation of PR and RR are 126.11 and 18.98 and 31.64 and 10.99 respectively. The association of the evolution for PR and RR was estimated to be (ρ=0.7054) which is statistically significant with 95% CI of (0.642, 0.769). PR and RR showed a decreasing pattern over time in both joint and separate models. Furthermore, a positive and significant association was observed between the two end points and all covariates except LVEF and time. Finally, to identify associated effect fitting joint model for paired endpoints is recommended. Journal of Biometrics & Biostatistics J o u rn al of Bio metrics & Bistatis t i c s


Background of the study
Abnormal respiratory rates and changes in respiratory rate are broad indicators of major physiological instability, and in many cases RR is one of the earliest indicators of this instability. Therefore, it is critical to monitor RR as an indicator of patient status. RR performs at least as accurately in identifying patients at risk of these adverse events as PR and the SBP. A RR of greater than 24 breaths per minute is able to identify approximately 50% of patients at risk of serious adverse events with 95% specificity. Although the main function of the respiratory system is gas exchange, a broad range of factors can affect ventilation. In patients with CHF, an increase in RR can warn of impending pulmonary edema, or fluid in the lungs, which is a common debilitating symptom of CHF as it was stated in American Heart Association [1].
Heart rate (HR) is among the many vital signs (respiration rate, blood oxygen saturation, arterial blood pressure, etc.), one of the most commonly measured and monitored. Whatever will be the sensing principle or the monitoring method used, data referred to the HR can be considered the primary vital sign information which is needed on a patient approach in both emergency and clinical situations. Gorgas [2] stated that, HR data are used to measure anomalous rate or irregular PR (arrhythmias) or heart block. The HR or PR represents the number of times the heart beats in a certain period of time. It is usually measured in minutes, and normal resting HR is approximately 60 to 80 beats per minute. It can go as high as 100 in a healthy adult and as low as 40 in athletes as it is described in American Heart Association and Gorgas [1,2]. The HR can be measured in various areas of the body, but the two most common sites are the wrist and neck. A lower HR is associated with a stronger and healthier heart. A lower HR means the heart is not pumping or working hard to deliver blood and oxygen to the body. The pulse can be lowered through regular exercise, and there are also breathing exercises to lower the heart. Take slow deep breaths to lower the pulse.
Dennison [3] stated that, heart failure (HF) is a condition in which one or both ventricles cannot pump sufficient blood to meet the metabolic needs of the body. HF, also known as CHF, is a chronic condition that develops over time. In some cases, the heart can't fill with enough blood; in other cases, the heart can't pump blood to the For instance, Thiébauta et al. [9] investigated the bivariate random effects model between the evolution of CD4 and HIV RNA and he reported the bivariate random effects model was significantly better than two separate univariate random effects models with (p-value<0.0001). In addition, the joint mixed effect model on evolution of occurrence and prevalence of antimicrobial resistant zoonotic agents were executed by Ferrari and Cribari-Neto [10]. They used beta-regression to illustrate the joint evolution on both outcome variables and they reported that, the correlation was estimated to be 0.95, with 95% confidence interval [0.414, 0.997] showing that the correlation was positively significant. Thus, there was a strong correlation between percentage resistant and prevalence and that both were increased with time. That correlation however ignores the effect of time.
Furthermore, Lambert and Vandenhende [11] studied that, the hemodynamic effect on diastolic blood pressure (DBP), systolic blood pressure (SBP) and HR. These three responses were measured repeatedly over time on 10 healthy volunteers during the dose escalation. The available covariates included in the study were sex and the concentration of drug in the plasma at time of measurement. The analysis was focused on the safety data, more safety data and more precisely on assessment of drug HR, in beats/min, (DBP) and SBP (in mmHg) for the ten subjects in the treatment arm. These measurements were taken before the first dose on day 1 and 4 hours after the morning dose on days 6-8, 12-14, 18-22. Thus, twelve repeated measurements were recorded per subject for each of the three outcome variables. In addition, the drug concentration (in ng/ml) was measure in plasma at the same times and sex was additional covariate. First, the evolution of DBP, SBP and HR were separately analyzed. And time did not appear explicitly associated with regression parameters. Indeed, time was only used to describe serial dependence between the repeated measurements and serial dependence was only found necessary to model HR profiles. In this dose escalation study, drug concentration tends to increase with time. For this reason, the effect of time appeared indirectly in the model as it was associated with the variation of the drug concentration in plasma. Gamma distribution was selected to fit the evolution of HR and the covariates considered were location drug concentration and sex.
Lambert and Vandenhende [11] reported that, the marginal mean HR was significantly smaller for men than for women but not significantly related to the drug concentration. They suggested that, the choice to the normal copula as the dependence structure could easily be specified through the variance-covariance matrix. The dependence between any two of the three outcomes measured by a parameter ρ with ρ ≤ |1| was again related to Kendall's tau of Lambert and Vandenhende [11]. Two Joint models of HR with SBP and with DBP were modeled. Thus, they reported that there was no significant association between HR and SBP but there was significant positive association between HR and DBP with a fitted Kendall's tau equal to 0.53 before treatment and 0.07 when there was drug in the plasma. There was no significant effect of sex on HR (PR) and DBP. In addition, joint model for SBP and DBP was fitted and there was a significant positive association between two variables with a fitted Kendall's tau equal to 0 and 0.42 for females before and after drug administration respectively and 0.22 for males no significant treatment effect on the association parameter was detected. Then in line to this the joint mixed effect model of two PR and RR is modeled.
Njagi et al. [12] analyzed joint modeling on the risk of rehospitalization and the mean number of times a patient's HR measurements which was classified as "abnormal", with LVEF as a baseline covariate for chronic HF data. He analyzed jointly model rest of the body with enough force. Patients with CHF have a poorer quality of life and shorter life expectancy compared with those of the same age in the general population. CHF is a chronic, debilitating illness, with ever-increasing prevalence in the aged peoples.

A joint mixed effect model
Repeated observation of multiple outcomes is common in biomedical and public health research. Such experiments result in multivariate longitudinal data, which are unique in the sense that they allow the researcher to study the joint evolution of these outcomes over time. Laird and Ware [4] stated that, in many circumstances, more than one response variable is followed longitudinally, and analyzing all jointly may be beneficial. Molenberghs and Verbeke [5] also stated that, until recently, methods for multiple longitudinal outcomes have largely been based on simple approaches where each outcome is analyzed separately, or by reducing the dimension of the multiple outcomes through a factor analysis or principal components type of approach.
Reducing the dimension of the multiple outcomes is also easy to implement, and can quite often capture much of the correlation between outcomes. Another frequently used method is stated by Gueorguieva and Agresti [6] to introduce random effects, but instead of sharing the random effect across the longitudinal responses, use separate, but correlated random effects in the longitudinal responses. Intrinsically multivariate questions concerning relationships between outcomes and the joint influence of covariates on them may be easily answered by fully exploiting the multivariate nature of the data through joint models.

General objective
The main objective of this study is to investigate the joint evolution of PR and RR of CHF patients and identifying the potential risk factors affecting the two end points.

Specific objectives
The specific objectives of the study are: To explore the evolution of PR and RR of CHF patients over time separately To fit a mixed effect model for the PR and RR and identify the associated factors To fit a joint model for PR and RR of CHF patients and to compare with separated models.

Joint model for longitudinal data
A joint multivariate normal distribution was considered for the corresponding latent variables and each outcome was analyzed with a marginal dose-response model. The covariance matrix takes into account the correlation between outcomes and the correlation due to clustering. That was an important improvement of Catalano and Ryan; Fitzmaurice and Laird [7,8] as model estimates of the correlation between outcomes and evolution of these correlations with dose were available. Hence, in relation to those literatures, the joint model for two symptoms of CHF i.e. PR and RR was considered to assess and identify both estimate of the correlation between two outcomes and the evolution of these correlations with a certain treatment throughout the time.
the recurrent time-to-re-hospitalization and a count version of the dichotomized longitudinal HR by understanding re-hospitalization is important in HF management. HR was first dichotomized into "normal" (50-90; coded 0) and "abnormal" (values higher than 90; coded 1). Values less than 50 were not considered for that analysis. During each period in which a patient was not under hospitalization, the number of times that the patient's HR measurements were classified as "abnormal" was enumerated, generating a count response. Notice that patients who were re-hospitalized and discharged at least once in the course of the study had at least 2 periods in which they were not under hospitalization, separated by a period of hospitalization. As a covariate, they considered the baseline LVEF. LVEF indicates the fraction of blood being pumped out of the ventricle with each contraction.
Njagi et al. [12] first looked at the results from the extended model. The test for a joint effect of ejection status on both processes was not statistically significant (p=0.1650), and therefore they concluded that there is no significant evidence of a joint effect of ejection status on both the mean number of abnormal HR measurements and the risk of re-hospitalization. Based on exponentiation of the relevant parameter estimate, the mean number of abnormal HR measurements in patients with reduced ejection was found to be 3.3531 times that of patients with preserved ejection. That effect was at borderline statistical significant (p=0.0594). The risk of re-hospitalization for patients with reduced ejection was obtained, by also exponentiation the corresponding parameter estimate, as 5.5168 times that of patients with preserved ejection; however, that effect was clearly not significant (p=0.6498).
Njagi et al. [12] then compared the results from the extended and the conventional model. Based on an AIC-based comparison, they observed that their extended model provided improvement to model fit, without compromising parsimony. There was impact on both the point estimates and standard errors. As they noted, the effect of ejection status on the mean number of abnormal HR measurements was borderline significant under the extended model; however, the case was quite different under the conventional model (p=0.0901). There was also a remarkable difference in the scale factor; it was highly significant under the conventional model (p=0.0022), while that was clearly not the case under the extended model, as they mentioned. However, in terms of the hypothesis of a joint effect of ejection status on both processes, the two models had provided close results; (p=0.1650; 0.1648) for the extended and the conventional model respectively.

Data source
The data were obtained from Ayder Referral Hospital of Mekelle University of CHF patient Clinic, north of Ethiopia in Tigray region. Mekelle University is one of the higher learning institutions found in Tigray National Regional State. The longitudinal data are extracted from patients' chart which contains epidemiological, laboratory and clinical information of all CHF patients under different drug levels follow-up including a detailed heart failure history.

Study variables
Response variables: Pulse rate and reparatory rate, Covariates (independent variables), Eight covariates are used for both separate and joint analyses. Four of these covariates are continuous while four of them are categorical covariates.

Statistical analysis technique
Different statistical analysis including both descriptive and inferential statistics, such as: summary statistics, data exploring and model comparison have been used in this study. Joint random effects with LMM have been modeled to infer the joint effect of bivariate longitudinal outcomes of CHF patients. Finally, Data were analyzed using SAS and R.

Descriptive statistics and data exploring
Data exploration is a very helpful tool in the selection of appropriate models. Thus, individual profiles plot, the mean profile plot and exploring the random effects and other data exploratory analysis for the data sets have been considered.

Separated linear mixed effect model
A mixed linear model is a generalization of the standard linear model used in the GLM procedure, the generalization being that the data are permitted to exhibit correlation and non-constant variability. The mixed linear model, therefore, provides the flexibility of modeling not only the means of your data but also their variances and covariance. The Linear Mixed Model is also an extension of the Linear Model that allows for incorporation of random effects and is represented in its most general fashions by Molenberghs and Verbeke [13].
where, x i and z i are the fixed and random design of covariates, respectively, β is a vector of unknown fixed effects, i γ is a vector of unknown random effects and i ε is the unknown random error. β Represents parameters that are the same for all subjects; i γ represents parameters that are allowed to vary over subjects.

Joint model for bivariate continuous longitudinal data
In many circumstances, more than one response variable is followed longitudinally, and analyzing both jointly may be beneficial. Until recently, methods for multiple longitudinal outcomes have largely been based on simple approaches where each outcome is analyzed separately, or by reducing the dimension of the multiple outcomes through a factor analysis or principal components type of approach. Bivariate linear mixed models are useful when analyzing longitudinal data of two associated markers. In this paper, a bivariate linear mixed model including random effects and independent measurement error for both PR and RR was presented. Longitudinal data are often collected in epidemiological studies, especially to study the evolution of biomedical markers. Thus, linear mixed models were stated by Laird and Ware [4] which recently available in standard statistical packages (Littell et al. [14] are used to take into account all available information and deal with the intra-subject correlation).
Extension to bivariate Case: Now under bivariate set-up two symptoms of CHF (PR and RR) as outcome variables are observed in each occasion. The two end points were longitudinally measured as a vector of responses, Y i (t), at each occasion with this model: ) 0 Where, 2 2 x ∑ is the variance covariance matrix of 2 endpoints , the response vector for the subject i, with Y ki the n ki vector of the end points k (k=1, 2) with n 1i =n 2i =n i n 1i= n 1i so model for bivariate longitudinal Gaussian data is: refer to the population means at time t. We assume that random effects are jointly distributed as follows: Where, D, the covariance matrix of the random effects, has the following structure: The error components are uncorrelated and not associated with the random effects Clearly, the correlation between the evolution of Y 1 and Y 2 is given by: And the marginal correlation between Y 1 and Y 2 at time t is given by: It is not difficult to comprehend that as the number of response variables (or the dimension of multivariate response) increases, the number of covariance parameters increase exponentially and the problem of estimation of covariance parameters becomes more and more difficult.

Estimation methods
Estimation for separate mixed effect model: Estimation of the parameters in LMM is usually based on maximum likelihood (ML) or restricted maximum likelihood (REML) estimation for the marginal distribution of Y i which can easily be seen to be ( ) Note that model LMM implies a model with very specific mean and covariance structures, which may or may not be valid, and hence needs to be checked for every specific data set at hand. Note also that, when , with I ni equal to the identity matrix of dimension ni, the observations of subject i are independent conditionally on the random effect bi. The model is therefore called the conditional independence model. Even in this simple case, the assumed random-effects structure still imposes a marginal correlation structure for the outcomes Y ij Even if all i Σ equal to 2 ni I σ , the covariance matrix in ( ) is not a diagonal matrix, illustrating that, marginally, the repeated measurements Y i j of subject i are not assumed to be uncorrelated. The marginal mean and marginal variance-covariance matrix of the vector Maximum likelihood estimation: -Suppose a random sample of N observations is obtained from a linear mixed effect model as defined above, then Brady et al. [15] defined that the likelihood of the model parameters, given the vector of N observations, is defined as: Then, the MLE of  β on combining all the information from all the N subjects equals.
Where det refers to the determinant and the elements of the V i matrix are functions of the covariance parameters in ϴ.

Gaussian quadrature
The Gaussian Quadrature approximates the integral of a function, with respect to a given kernel, by a weighted sum over predefined abscissas for the random effects. Unlike other numerical integration techniques, the abscissas are spaced unevenly throughout the interval of integration. With a modest number of Quadrature points, along with appropriate centering and scaling of the abscissas, the Gaussian Quadrature approximation can be highly effective see Abramowitz and Stegun for details [16]. Pinheiro and Bates [17] also suggested that, in the particular context of random-effects models, so-called adaptive Quadrature rules, where 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, they 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: bi is q×1 dimensional vector of unknown random effects, b i ~ N (0, D) D is variance-covariance matrix of f(z) and for φ (z) denotes the density of the multivariate standard normal distribution

Model comparison technique
In order to select the best and final model which is appropriately fits with the given longitudinal data, it is necessary to compare the different models by using different techniques and methods. Hence, models are compared with Akaki Information Criteria (AIC), the Bayesian Information Criteria (BIC), and the Likelihood ratio test methods for nested were used at 5% level of significance.

( )
Where, -2 logL is twice the negative log-likelihood value for the model P: number of estimated parameters and npar: the total number of parameters in the model N: Total number of observations. Smaller values of AIC and BIC reflect an overall better fit.

Data descriptions and descriptive statistics
In this cohort study, socio-demographic and clinical data of 264 patients whose age are 15 years and above receiving preferable drugs to improve the symptoms of CHF from Sep. 2012 to Aug.2013 in Ayder referral Hospital of Mekelle University at baseline were considered. The two symptoms of CHF, PR and RR have been used. These longitudinal response variables were measured for at least 4 visits. There were a total of 6494 visits from 264 subjects in the CHF treatment, The number of visits per subject varied from 4 to 46 weeks with a mean follow-up time of 16.103 (SD=10.85), but the time interval for all patients is not equally observed i.e. visits were unequally spaced. The sample sizes at the six consecutive time points were 264, 232, 188, 206, 245 and 184. There is a sharply increasing degree of missing data over time due to deaths, dropouts, missed clinic visits and transferring to other hospital and also there is admitting and readmitting of the patients. Moreover, eight covariates which four continuous and four categorical were encompassed (Table 1).
The baseline characteristics of patients are displayed in Table  2. Among these patients, more than half 151 (57.2%) of them are females and 113 (42.8%) are males; refer Table 2 for the rest details. As indicated in Table 3, the longitudinally measured symptoms of CHF, PR (in beats/minutes) and RR (in breaths/minutes) were considered as bivariate responses. They were measured approximately every week at the study entry, and again a common measuring time is used for all patients. Hence, the baseline mean and standard deviation of PR were 126. 11

Exploratory data analysis
Figures 1 and 2, indicated that there is decreasing trend on the PR and RR of CHF patients throughout the follow up. Both PR and RR that were the heaviest at the beginning tends to be turn down throughout the follow up. That means, the variability of the measurements, at the beginning (baseline) of the follow up were highly decreased relative to the end of the follow up. Likewise, there is variation with in groups throughout the time by decreasing the PR and RR from week to week. According to the average trend line that putted on the individual profile plot with blue color, the PR and RR of the CHF patients were declined throughout the time. Hence, there is the negative evolution on both PR and RR over the time. Furthermore, the average trend line is almost near to straight downward line indicating linear relationship with absence of quadratic effect on the negative evolution of PR and RR of CHF patients.
Besides plotting the PR over follow up time in weeks, it is also useful to include different subgroups on the same graph to illustrate the relationship between the response variable PR and an explanatory variable sex over follow up in weeks as it was shown on Figure 3 of (P.a). Thus the mean profile plot of PR by sex presented on Figure 3 of (P.a) indicated that PR decreased for both men and women over the follow up. However, the slope for the women seems more visibly higher than the slope for the men from baseline up to end of follow up which did not indicated the interaction effect as did not crossed each other.

NO. Variable Description
Value/codes   Similarly, it is also useful to illustrate the relationship between the response variable PR and an explanatory variable NYHA over follow up time in weeks as it was shown on Figure 3 of (P.b). Thus the mean profile plot of PR by NYHA presented on Figure 3 of (P.b) indicated that all categories have decreasing trend on PR over follow up time in weeks. However, the slope for the NYHA class IV seems to be higher than the slope of the others from baseline up to end follow up which does not indicate the interaction effect as they did not crossed each other's.
Generally, lower NYHA class has lower PR whereas higher NYHA class has higher PR throughout the time. As it was shown on Figure 3 at (P.c), the mean interaction plot of PR by residence for CHF patients indicated that, even if there is a down ward trend for both categories, almost both categories have the same effect on PR, since the plot for both rural and urban was almost overlapped and it seems like the absence of significant difference between rural and urban on the evolution of PR. According to Figure 3 at (P.d), the mean interaction plot of PR by diagnosis for CHF patients indicated that, almost all categories have the same effect on PR, since the plot for all levels was overlapped. However, the PR of patients with previous diagnosis history of severe anemia and CHD were relatively higher than others and ACF. Generally, PR in all categories has declining trend as it was shown by Figure 3 at (P.d). Furthermore, similar explanations can be given for the patterns of RR as shown on Figure 4. Thus the mean profile plot of RR sex presented on Figure 4 at (P.1) indicated that RR decreased for both men and women over the follow up. However, the slope for the women is above the slope for the men from baseline up to end of the follow up which did not reflect the interaction effect as plots did not crossed each other ( Figure 4).
As it is indicated on Figures 5 and 6, even though, most of points overlapped for both slope and quadratic slope plot, there is certain variability in the slope which mean that, considering slope random effect for this model is important for both PR and RR models. But as it is clearly depicted on the subject specific intercept plot, there is high variable in the intercept and it became crucial to add both intercept and slope in the random term.

Model selection
A primary goal of model selection is to choose the simplest model that provides the best fit to the observed data. There may be several choices concerning which fixed and random effects should be included in a linear mixed model (LMM). There are also many possible choices of covariance structures for the D and R i matrices.

Model fitting for fixed and random effects
The top-down strategy was used to select statistically significant covariates for the two independent mixed effect models with outcome variables PR and RR. The models based on only fixed effects were selected with constant random effects at first and then after, the significance of random effects was also checked. Hence, passing through procedures, the model which included the covariates time, age, sex, NYHA, LVEF, weight and interaction of time with LVEF and weight for fixed effect with subject specific random intercept and random slope for both models with outcome variables PR and RR were preferred regardless of relatively small values of AIC=42645.03, BIC=42746.71 and log-Likelihood ratio test with P-value of <0.0001 for model with outcome variable PR and AIC=24344.10, BIC=24445.78 and log-Likelihood ratio test with P-value of <0.0001 for model with outcome variable RR. The saturated model including the covariates diagnosis history and place of residence in addition to the covariates that were included in the reduced model were fitted and compared with reduced model. But saturated model did not improve the model. Thus, statistically insignificant covariates such as diagnosis history and place of residence were excluded from the final model. Finally, quadratic fixed effect and random effect did not improve the models and that's why it is also excluded in the model. In addition, the ML method was preferred.

Model selection with covariance structure for the best model
Comparing different covariance structure for the best model for Pulse Rate: According to Final Model in Joint Case: Note: The Notations used in the model: T=time, A=age, SF= sex (female), W= weight, L=LVEF, NCII=NYHA class II, NCIII= NYHA class III and NCIV= NYHA class IV The PR and RR outcomes were modeled with the set of covariates, and the results were described in Table 6. The final model was somewhat complex and included 11 fixed effect parameters for both outcome variables PR and RR including intercept and Time to Weight & LVEF interactions for both separate and joint mixed effect models.

Results of joint mixed effect model
A joint mixed effect model for the two symptoms of CHF syndrome PR and RR was fitted with an unstructured variance-covariance structure. This model is the same as the separate model except the sets of random intercepts and slopes for each response are now correlated rather than independent. This model was fitted allowing for a linear time effect for each covariate that was selected as a fixed effect in the separate linear mixed model. The subject specific random intercepts and random slopes were fitted to account for within-subject correlations.
According to Table 6, the fixed-effect intercept coefficient 10 β = 108.58 (S.E.=4.384) represents an estimate of the average PR at time=0 and excluding all covariates in the model. Likewise, the fixed-effect intercept coefficient 20 β = 36.57 (S.E.=1.78) represents an estimate of the average RR at time=0 and excluding all covariates in the model. All parameters are statistically significant except there is no evidence of a significant relationship between NYHA class II and PR (P=0.0961). Among all covariates, Time, Age, and LVEF were negatively associated with both outcomes that mean the repeatedly follow up made a particular decrease on both outcomes with (P<0.0001).
In addition, sex was significantly associated with both PR and RR outcomes; thus, female patients had 4.528 points higher over evolution of PR (P=0.0008) and 1.734 points higher over evolution of RR (P=0.0024) compared to males. Moreover, there was evidence of a statistically positive relationship between weight and both PR ( In similar way, time-weight interaction was also significantly and positively associated with PR (β = 0.016; P<0.0001) and with RR (β = 0.003; P<0.0001). Generally, as it is indicate in the results in Table  6, both PR and RR have decreasing pattern throughout the follow up with respective clinical treatments. This concept indirectly indicated the improvement on risk of congestive heart failure because the lower value of both symptoms PR and RR is directly related with a stronger and healthier heart.

Variability of error and random effect in joint model
Alike parameter estimation and testing, variability analysis of both fixed and random effects are also another important aspects. High variability is the indicator of less accuracy or high error on prediction of the association of outcome evolutions with respective risk factors. Then as it is shown in ; 95%CI= [1.598, 1.719]) for RR. Thus, the variability due to subject specific random intercepts is higher than that of random slopes for both models. The random effect variability is greater on PR than RR.

Associated (common) effect parameters
By referring Table 6, based on 6494 pair symptoms of CHF assessments from 264 subjects, a substantial correlation (ρ=0.7054, S.E.=0.032) with 95% CI: [0.642, 0.769] between the PR and RR within the same subjects is noted. From the random effects, it may be seen that variability is relatively higher for PR than RR. The same may be said of the covariance for subject specific random intercept of PR and RR with (  Table 6. With the joint mixed effect model for the two endpoints of CHF, it is possible to investigate how the evolution of PR is associated with RR. Hence, the association of the evolution (AOE) is to be estimated 0.7054(S.E.=0.032, p-value<0.0001). Not only that but also it is possible to determine how the association between the two symptoms of CHF (PR and RR) evolves over time; thus, the evolution of the association (EOA). For instance, at baseline the evolution of the association was 0.45029 and at first, second and third weeks follow up it increased into 0.4505118, 0.4508186 and 0.4523747 respectively indicating the evolution of association between PR and RR over the time. In addition to that, the evolution of the association (EOA) throughout the time is well visualized as it is shown on the marginal association plot of Figure 7; there is the positive evolution of the association between two outcomes PR and RR. Thus, the association positively evolved over the time. Generally, there is evidence that time has reasonable effect on association of evolution of both outcomes.

Results of separate mixed effect model
Technically, the separate models were fitted for the two outcomes, PR and RR together but assuming that ρ=0, which is entirely equivalent to fitting the models separately or independently as their results were shown in Table 7. Hence, as interpretations for the models those modeled independently for PR and RR is entirely equivalent to that of separate models by assuming ρ=0.

Variability of error and random effect in separate model
By referring Table 6, even if there is slight difference variability's in joint mixed model results, there is almost similar results are computed for separate one. Thus, the subject specific random intercept variance is estimated to be 141.27(S.E.=13.14) with 95% CI of (118.68, 171.03) for PR and 31.63 (S.E.=2.915) with 95% CI of (26.612, 38.222) for RR. The subject specific random slope variance is estimated to be 0.182(S.E.=0.02) with 95% CI of (0.148, 0.229) for PR and 0.066 (S.E.=0.007) with 95% CI of (0.054, 0.083) for RR. The estimated variance of the random error is ( Finally, similar to that of joint mixed model results, the variability due to subject specific random intercepts was higher than that of random slopes for both models. The random effect variability is greater on PR than RR. Comparison of separate and joint or shared mixed effect models LRT=-2LL(PR&RR)-(-2LL(PRRR)=66929.1-66802.9=126.2 (15) Here, both separate and joint mixed effect models have been considered and parameter estimates for the separate and joint models are summarized in Table 6. Technically, the separate models were fitted for two outcomes together, but assuming that ρ=0, which entirely equivalent to fitting the two independent models separately as results were shown in Table 7. It also allows for a single likelihood for the model parameters enabling direct comparison with the correlated bivariate model fitted subsequently. Clearly, PR and RR show a strong positive relationship as evidenced by the correlation of the random effects in joint mixed models. In addition, likelihood comparison shows a convincing improvement in model fit, when random effects are allowed to correlate. Comparing the separate and joint models, although parameter estimates for both outcomes are nearly equivalent, small changes are observed in parameter of some covariate. When comparing the results from the separate settings to the results from joint settings, there are several points of interest. The -2log-likelihood value corresponding to the two separate models (i.e. fitted as a joint model but assuming ρ=0) was equal to 66929.1 and the -2loglikelihood value for the joint model was 66802.9. Hence, the joint random effect model of the two symptoms of CHF, PR and RR was significantly better than two separate random effect models of PR and RR (-2LL=66802.9 vs. 66929.1; LRT=, 126.2 DF=4, P-value<0.0001). With regards to AIC, the joint model (AIC=66870.9) is also indicated as a better fit than the separate model (AIC=66989.1). Notice how the joint model of two symptoms of CHF i.e. PR and RR seem to decrease the variability in the random effects, this may be seen in Table 6. Taking into account the standard errors for the variance and covariance estimates, the joint model in general allowed for more accurate prediction (small errors) of the variability in the random effects, though just slightly.
Comparing the fixed effects for the separate and joint mixed models, some important things may be considered for the two symptoms of CHF patients. First, and foremost, there is the question of whether the different models reached the same bottom line conclusion. Comparing the covariates between two types of models will yield further information of interest. Both separate and joint models found a significant relationship between weight and PR and RR. Weight was positively associated with PR (β =0.538 compared to 0.521), hence, the 95%CI (0.474, 0.603 vs 0.458, 0.584) and with RR (β =0.054 compared to 0.059); 95%CI (0.038, 0.07 vs 0.043, 0.075); therefore, the 95% CI is equally tighter for both models. Sex (females) was positively associated with RR (=1.605 compared to 1.734), hence, the 95%CI (0.47, 2.739    ; therefore, the 95% CI is tighter for joint models relative to males. Both models also concluded a nominal decrease with regards to LVEF for both PR (β =-0.307 compared to -0.295), 95%CI ((-0.339,-0.275 vs -0.327,-0.263) and RR (β =-0.075 compared to -0.075), 95%CI (-0.082,-0.067 vs -0.082,-0.067) has no difference in both models as it was shown in Table 6. All parameters except NYHA class II on PR are statistically significant in both models. Finally, similar explanation could be given through the rest covariates displayed in Table 6.

Model diagnostic checking
Diagnostic checking and Residual plot for fixed effects: Different diagnostic checking plots for the final separate mixed linear models of PR and RR are presented in Figures 8-13 Thus; the result is explained as shown below. According to the Figures 8 and 9, plot of fitted versus standardized residuals, even if there are some outliers, it was indicated that the variability of the errors in both PR and RR were almost nearly constant. That means the errors did not far deviate from each other. Distances of individual residuals were equally far from the horizontal lines. Furthermore, according to the probability plots those were shown on Figures 10 and 11, even if the points were compacted at the two end tails for both outcomes PR and RR, the normality assumption was supported through the upward nearly straight line of normal plots. Similarly, based on the normal probability plots of random effects with subject (MRN) specific random intercepts and random slopes those are shown on Figures 12 and 13, even if it seems a slight deviation of normality at the bottom tail on the random slope (Time) for RR that is not that much worse deviation. Hence, there is no problem with normality assumptions of both random intercepts and random slopes for both PR and RR models and the normality assumption are almost fulfilled.

Discussions
Based on different well organized literatures and analysis that were included in this paper, some discussions and review of works is organized as follows.
This study was conducted on the title of a joint model for a longitudinal PR and RR of CHF patients in Ayder Referral Hospital of Mekelle University. In summary, a joint mixed effect model for paired outcomes with the sets of both continuous and categorical covariates and the interaction of time with weight and LVEF is presented. This model extends previous work by accommodating longitudinally measured two main symptoms of CHF as outcome variables. According to some related works that were reviewed, even if old age, sex, smoking, hypertension, diabetes, obesity, valvular heart disease, and CHD were considered as the important risk factors for CHF, but as a result of the absence of those particular covariates, only some of those covariates were included. For implementation a necessary computational procedure is developed. Using the proposed methods, the influence of different covariates which were listed earlier in this paper is examined. With the proc mixed statistical methods, the influences of the covariates on longitudinally measured bivariate outcomes PR and RR is executed. Since joint model building usually starts from separate models for each component, initially each data are analyzed separately. Such separate analysis is preferred for several reasons. Firstly, it helps to specify the mean response of the model. Secondly, the random effects to be included in the longitudinal model can be easily determined, and thirdly initial values to be provided for the joint models can be obtained.
The finding provides direct evidence that decreasing in LVEF (in %) is the primarily driver of the risk of CHF by causing reasonable increase on both PR and RR longitudinally throughout the follow up. The finding is consistent with the latest literature which suggested by Njagi et al. [12] based on exponentiation of the relevant parameter estimate, the mean number of abnormal HR measurements in patients with reduced ejection was found to be 3.3531 times that of patients with preserved ejection. That effect was at borderline statistical significant (p=0.0594). The test for a joint effect of ejection status on both processes was not statistically significant (p=0.1650) but in contrast to that in this finding there is statistically negative significant association between the LVEF and both PR and RR. As they noted, the effect of ejection status on the mean number of abnormal HR measurements was border-line significant under the extended model; however, the case was quite different under the conventional model (p=0.0901), this statement contrasts the finding of this study.
Lambert and Vandenhende [11] reported that there was no significant association between HR and SBP but there was significant positive association between HR and DBP with a fitted Kendall's tau equal to 0.53 before treatment and 0.07 when there was drug in the plasma. The finding is consistent with it because PR and RR have significant positive association. Furthermore, there was significant association between sex and both PR and RR in contrast to Lambert and Vandenhende [11] that there was no significant effect of sex on HR and DBP.
The finding provides direct evidence of strong correlation between two symptoms of CHF (PR and RR) estimated to be 0.7054(70.54%) with 95% CI of (0.642, 0.769). Thus, the joint mixed effect model was better fit than two separate random effect models. This finding is consistent with the previous literatures that was studied by Thiébauta [9] on bivariate mixed effect model or first-order autoregressive process and independent measurement error for both markers of CD4 and HIVRNA in HIV patients (p 10 . Similarly the finding is also consistent with the previous literatures of Ferrari and Cribari-Neto [10] studied on application of joint models for resistance and prevalence a strong correlation between percentage resistant and prevalence and that both increase with time. The correlation is estimated to be 0.95, with 95% confidence interval [0.414, 0.997] showing that the correlation is significant. That correlation however ignores the effect of time. In contrast to this finding, statistically significant marginal correlations over time have increased throughout the time. Finally, joint mixed model was preferred to find and identify joint evolutions in this finding and this is consistent to Njagi et al. [12] who compared the results from the extended and the conventional model. Based on an AICbased comparison, they observed that their extended model provided improvement to model fit, without compromising parsimony. There was an impact on both the point estimates and standard errors.

Conclusions
The main aim of this thesis was to develop joint mixed effects model for paired symptoms of CHF (i.e. PR measured in beats/minutes and RR measured in breaths/ minute) as outcome variables. Toward this goal, the previously introduced joint model allows the joint modeling of mixed model for PR and RR with specification of subject specific random intercepts and slopes. Then it can be generalized the joint model to the longitudinal data, which necessitates the modeling of association between the continuous outcomes (PR and RR) considered very important. This is accomplished with incorporation of random effects (i.e. subject specific random intercepts and random slopes (time), by excluding quadratic random slopes) in individual linear mixed effect models for PR and RR. The unstructured covariance structure was preferred to fit both separate and joint mixed effect model. Estimation of the fixed and random effects was described, along with formal definitions of the association in the evolution (AOE) of the two responses and the evolution in the associations (EOA). Thus, the question of AOE and the EOA of the PR and RR were clearly addressed.
After passing many procedures, among all covariates diagnosis history and place of residence were excluded in final models because of their insignificant effect on both outcomes but the rest covariates such as time, age, sex, weight LVEF, NYHA and interaction of time with weight and LVEF were included in final models. Out of those covariates, time, age and LVEF were found to be negatively associated with both outcomes in both separate and joint mixed model. Moreover, among all the covariates included in separated and joint mixed models, only NYHA class II were statistically insignificant on the evolution of PR. Non-zero covariance of random intercepts and random slopes explained the statistical significance of association between two outcomes. Likewise, it can be generalized that, the two outcomes have a strong positive correlation and the correlation was statistically significant. Thus, the joint mixed effect model was preferred because the joint mixed effects model is more flexible in allowing separate fixed and random effects for each response i.e. PR and RR through appropriate choice of potential risk factors (covariates) or fixed effect and random effects, while accommodating dependence in the longitudinal trajectories through dependence in the random effects. The baseline mean of the two symptoms PR and RR were out of the normal range for CHF patients but throughout the consecutive follow up of the clinical treatment, decreasing values of PR and RR has been shown. That decreasing trend on PR and RR indirectly indicated the reduction of the risk of congestive heart failure.
Finally, it is concluded that, joint modeling of longitudinal bivariate responses is necessary to explore the association between paired response variables like PR and RR. A usual problem with the joint modeling is failing to convergence because of large number of association parameter to estimate. Gradually, for future work, one might want to look at modeling the joint mixed model with correlated measurement errors which may violate the result of mixed effect when uncorrelated error is considered. Moreover, some one also might want to look at modeling more than two response variables over time. This issue typically can be implemented using modern computing methods for a joint model in which there are more than two response variables. However, with increasing response variables, there is an exponential increase in the amount of computing power necessary to produce estimates and the complexity is high.

Recommendation
As the selection of an appropriate statistical model is directly related to the qualities and nature of the data, in the case of limited quality data, the associations of factors or covariates with outcome variables could not assessed. Therefore, special intention should be given to the quality of the data.
A lot of investigators doing longitudinal research used to model repeated outcomes separately, to assess the evolutions of the outcomes through time by ignoring the associated effects. But it is recommended to check the associated evolutions in some case as the outcomes might have the association of the evolutions. Even if almost equivalent questions answered through joint model and separate model, joint model is able to address the same questions with better accuracy. And also address the additional and important concepts of AOE and the EOA of the outcomes. Thus, fitting joint model is recommended. In many cases including in this study, uncorrelated error is considered in modeling joint mixed models, but in some cases it is crucial to consider correlated error in models because using uncorrelated error model may display less accurate results if there is suspicion of correlated measurement errors in the data.

Acknowledgement
First and foremost I must give thanks to almighty God for all the blessings he bestowed upon me till today.
Secondly, Then I want to express my sincere gratitude and appreciation to my co-advisor Geremew Muleta (MSc.) for his advice and encouragement from inception up to completion.
Thirdly, I would like to take this opportunity to express my heartiest appreciation to Abba Thomas Weldesilassie and Mr. Tesfaye Debella (PHD Scholar) for their insightful guidance, consistent support, and impressive encouragement in grammatical and writing skill.
Finally, and most importantly my deepest gratitude goes to my mother, my brothers and sisters and all my friends who were greatly involved in my work in any means. Last but not the least my special thanks goes to Fitiwi Weldesilassie (Nurse) for he helped me considerably in obtaining necessary knowledge about various medical factors.

Dedication
This paper is dedicated, praying for their eternal life, to my beloved father, the late Hailu Fissuh, and to my cousin and my friend, the late Gaim Massa. My father was God fearing, loving, hardworking and outstanding in his concern for the common good. He is my inspiration, aspiration and the reason I am the man I am today. Daddy, I thank you and I love you. You are my hero, my teacher and my companion in spirit. Gaim was a brave young man, full of life and with great sense of humor. He was friendly, kind and magnanimous. Dear Gaim, I am sorry for your untimely death, but I remember and cherish all the good times with you.