The study revealed that, the median recovery time of pneumonia patients admitted at FHRH was 72 hours with minimum and maximum recovery time of 18 hours and 96 hours respectively. Out of the total sampled pneumonia patients, 90 (89.1%) were recovered from pneumonia. When we fit the cox proportional hazards model using the candidate variables: residence, birth order, age of mothers, education of mothers, danger signs, cooking place, comorbidity and severity were significant factors affecting time to recovery of pneumonia patients at 5% level of significance (Table 1).

In this study, three longitudinally measured response variables of pneumonia patients were considered. The linear mixed model was used for all the variables; pulse rate, respiratory rate and oxygen saturation of patients. The study was started by exploring the mean and variance structure of those longitudinally measured response variables. The three longitudinal measures of pneumonia were approximately measured every 6 h a day from admission up to hospital discharge of under-five pneumonia patients. All of 101 sampled under five admitted pneumonia patients were at risk of pneumonia up to the third visit time (t = 12 hour), this tells that, for this study, the minimum follow-up time at which the patient get the event of recovery was the third visit (t = 18 hours) and the number of patients getting the event increases, whereas the number of patients with at risk of pneumonia decreases through visit time.

The study also revealed that, the average values of RR and PR decrease, whereas the average values of oxygen saturation increase through the visit time. At the end of the follow up, the overall average values of RR, PR and Oxygen saturation were 50.55 bpm, 131.20 bpm and 90.18 mmHg with standard deviation of 12.55 bpm, 27.37 bpm and 6.11 mmHg respectively (Table 2).

Checking assumptions of the data is the first step in analyzing longitudinal data. Normal QQ plots in Fig. 1 shows that, the data for the three longitudinal outcomes were approximately normally distributed and then it is better to proceed to the next steps of the analysis.

Normal QQ plots of RR, PR and oxygen saturation

MLMM was fitted using three longitudinal measures of pneumonia (RR, PR and oxygen saturation) for under-five admitted pneumonia patients (Table 3). At 5% level of significance; marital status of mothers, smoking exposure of patients, breast feeding, severity, cooking place, comorbidity and visit time were significant factors related with longitudinal measures of RR. Age, residence, birth order, comorbidity, danger signs, vaccination, severity and visit time were significantly related with longitudinal measures of pulse rate. The variables that significantly related with longitudinal measures of oxygen saturation were; age at the base line, residence, comorbidity, danger signs, age of mothers, severity and visit time. Variables; severity, visit time and comorbidity were simultaneously associated with longitudinal measures of RR, PR and Oxygen saturation of patients.

The random part of MLMM shows the variance and covariance between rate of change and baseline values for the three longitudinal measures of pneumonia (RR, PR and oxygen saturation) were significantly different from zero which tells the existence of a relationship between a patients baseline standing between outcomes, rate of change between outcomes as well as, between baseline standing of one outcome and rate of change of the other outcome through follow-up time.

In the previous sections; determinants of the multivariate longitudinal measures of pneumonia as well as determinants of time to recovery of under-five admitted pneumonia patients were identified. The results of joint model analysis for multivariate longitudinal and survival data found in the Table 4, contains multivariate longitudinal and survival sub models. In the random part of MLMM, estimates of variance and covariance were different from zero, shows the existence of correlation between intercepts of outcomes, between rate of changes of outcomes and correlation between rate of change and baseline values of the three longitudinal measures of pneumonia (RR, PR and oxygen saturation).

Based results of Table 4, the average RR, PR and oxygen saturation of under-five pneumonia patients admitted at FHRH were 47.2660 bpm, 146.7431 bpm and 87.29 mmHg respectively when all categories are at their reference group. As age of patients increased by 1 month, the average RR and PR were significantly decreased by 0.38 bpm and 1.01 bpm respectively. Whereas, age was not a predictor of oxygen saturation. Coming from urban residence increases the average RR and PR by 1.70 bpm and 1.26 bpm respectively, whereas it lowers the average oxygen saturation by1.01 mmHg as compared with rural residency, keeping other variables constant. Being first child significantly rises the average RR and PR by 1.59 bpm and 3.09 bpm respectively; whereas it lowers the average oxygen saturation by 1.72 mmHg as compared with being second or above child; other variables held constant. Being non-exposed by smoking lowers the average RR by 2.27 bpm as compared with patients exposed by smoking; keeping other variables constant, but had no information about PR and oxygen saturation.

Being non-comorbid significantly lowers the average RR and PR by 3.98 bpm and 3.64 bpm respectively, while it rises the average oxygen saturation by 2.33 mmHg as compared with being comorbid, keeping other variables constant. Having sever pneumonia at the baseline increases the average values of RR and PR by 5.46 bpm and 1.30 bpm respectively, whereas it lowers the average oxygen saturation by 1.03 mmHg as compared with those having non-sever pneumonia, other factors held constant. Having literate mother increases the average oxygen saturation by 2.70 mmHg as compared with those from illetrate mothers, held other variables as constant. Cooking food inside the living room lowers the average oxygen saturation by 2.11 mmHg as compared with those whose parents cook their food out of living room, keeping remaining factors constant. A unit increase in the number of visits lowers the average RR and PR by 0.19 bpm and 0.16 bpm respectively, whereas it rises the average oxygen saturation by 0.90 mmHg, keeping other predictors constant.

Getting vaccination lowers the average PR by 8.59 bpm as compared with unvaccinated by remaining other variables constant. Feeding exclusive breast within first 6 months decreases the average RR by 1.85 bpm as compared with no breast feeding. The estimated hazard ratio of patients from urban area relative to patients from rural area was 0.61 indicates, patients from urban residence were 0.547 times less likely to recover from pneumonia than patients from rural residence, other variables held constant. Patients without comorbidity were about 2.296 times more likely to experience the event of recovery compared to patients without comorbidity. Patients at the first birth were 0.284 times less likely to get the chance of recovery compared to patients at the second and above births, keeping other variables constant. As age of mothers increase by 1 year, experiencing the event of recovery increases about 2.462 times, other variables held constant. Exclusively breast feed patients within first 6 months of life were about 4.06 times more likely to get recovery as compared with patients having no breast feed.

Patients with severe pneumonia were about 0.206 times less likely to experience the event of recovery compared to patients with non-sever pneumonia, keeping other variables constant. The estimated values of association parameters γ_1 = − 0.297 (p-value = 0.0021), γ_2 = − 0.121 (p-value< 0.001) and γ_3 = 0.545 (p-value = 0.006) indicates; RR and PR were negatively associated with time to recovery, whereas oxygen saturation was positively associated with time to recovery of under-five admitted pneumonia patients.

Model comparison: The multivariate longitudinal sub-model was consistent with the results from the multivariate longitudinal analysis of RR, PR and oxygen saturation. The differences in magnitudes of the parameter estimates were negligible and there were some parameter difference in terms of statistical significance in separate MV longitudinal and separate survival model. But, longitudinal sub-model had narrow confidence interval which indicates that standard error is small for all significant predictors as compared to separate model in MV longitudinal and survival model. When evaluating the overall performance of both the separate and joint models in terms of model parsimonious and goodness of fit, the joint model was preferred as it has smaller standard error than the separate model. This result also supports the study done by [25, 26].

As Table 4 revealed, under MV joint model, estimate of the association parameters in the survival sub model was significantly different from zero (γ_1 = − 0.297, γ_2 = − 0.121 and γ_3 = 0.5452), this indicates that three longitudinal outcomes were correlated with time to recovery of under-five admitted pneumonia patients supported by [27,28,29], stats that the longitudinal and survival data are correlated. The joint model was more parsimonious fit than the separate model. Therefore, the joint model found preferable and parsimonious to fit the data better than the separate one [24] when the association parameter of the joint model is significant. Therefore, the final model for this study was joint model of MLMM and cox PH model.