Modelling Under-Five Mortality among Hospitalized Pneumonia Patients in Hawassa City, Ethiopia: A Cross-Classiﬁed Multilevel Analysis

Community acquired pneumonia refers to pneumonia acquired outside of hospitals or extended health facilities and it is a leading infectious disease. This study aims to model mortality of hospitalized under-5year child pneumonia patients and investigate potential risk factors associated with child mortality due to pneumonia. The study was a retrospective study on 305 sampled under-ﬁve hospitalized patients of community acquired pneumonia. A cross-classiﬁed multilevel logistic regression was employed with resident and hospital classiﬁed at the second level. Bayesian estimation method was applied in which the posterior distribution was simulated via Markov Chain Monte Carlo. The variability attributable to hospital was found to be larger than variability attributable to residence. The odds of dying from the community acquired pneumonia was higher among patients who were; diagnosed in spring season, complicated with malaria, AGE and AFI, in a neonatal age group, diagnosed late (more than a week). The risk of mortality was also found high for lower nurse: patient and physician: patients’ ratios.


Introduction
Pneumonia is a form of acute respiratory tract infection (ARTI) that affects the lungs [1]. It is among the main killers of children under the age five worldwide [2]. Community acquired pneumonia (CAP), which refers to pneumonia acquired outside of B Tariku Tessema tarikutessema1@gmail.com 1 Department of Statistics, College of Computing and Informatics, Haramaya University, Haramaya, Ethiopia hospitals or extended health facilities, is a leading infectious disease requiring hospital admission in both developed and developing countries including Ethiopia. According to World Health Organization (WHO), pneumonia is the second leading cause of death among the under-fives next to preterm birth complication [3]. It killed nearly to one million children under the age of five in 2015, accounting for 16% of all deaths of children under 5 years old [4]. Fifteen countries that are home for 55% of world's children shares 72% of the global burdens of pneumonia and diahrrea child deaths [4]. WHO [2] estimated that, about half of the deaths in children under five occurred in 2015 were caused by infectious diseases and conditions such as pneumonia, diarrhoea, malaria, meningitis, tetanus, HIV and measles.
Ethiopia was among 15 countries contributing to burden of child mortality due to pneumonia and diarrhea between 2015 and 2016 [4]. In Ethiopia, the child mortality rate is 67 deaths per 1000 live births [5]. Pneumonia has been the main killer of children under age of 5. In 2012 UNICEF [6] reported a 21% of the total deaths among children under five was due to pneumonia. A similar Pattern was observed in the 2000 Ethiopian Demographic and Health Survey (EDHS). The infant mortality rate in the 5 years preceding survey is 77 and under-five mortality is 123 deaths per 1000 live births for the same period. This means that one in every thirteen Ethiopian children dies before reaching age one, while one in every eight does not survive to the fifth birthday. However, these figures declined to 68 deaths per 1000 live births where pneumonia accounts for 18% of the total deaths [7].
Various studies have been expended to identify important risk factors of mortality and morbidity from the Pneumonia patients. Some of the literatures have used spatial and non-spatial analytic techniques to examine patterns of the illnesses, and the factors that determine it [8]. There is also literature that have used multivariate logistic regression to determine the factors that predict in-hospital mortality among patients who require hospitalization for the treatment of CAP [9]. These statistical methodologies are not capable to consider variability that could be observed in groups and does not take into account the dependency of observation in the same groups/cluster. Therefore, it is necessary to use different approach that may explore the important determinants of child mortality due to pneumonia at patient, hospital and patient's residence level.
This study is, therefore, intended to model the variability of under-five child mortality (U5M) of hospitalized pneumonia patients observed at resident or hospital, and identify risk factors associated with U5M due to pneumonia using cross-classified multilevel model. A cross-classified model is a multilevel model where individual or lower level units are nested to two different groups [10]. In this model individual is nested into a classification of two groups, that is, patients admitted to one hospital come from different sub-city, likewise patient from one sub-city admitted on different hospitals. Hence, source of variation at individual/patients and group/hospital, subcity levels can be considered. This study also contributes to the existing literature by providing methodological demonstration of cross-classified multilevel models in hospitalized patients of community acquired pneumonia.

Description of Study Area and Data
This study had been carried out at Hawassa City, capital of Southern Nations Nationalities and People's Region (SNNPR), Ethiopia. The city, which is the economic and cultural hub of the region, has a total area of about 157.2 km 2 divided in to 8 subcities (Kifle-ketema). These sub-cities (see Table 1) are divided mainly based on the geographical structure. Under Hawassa City Administration, there are a total of five hospitals and six health centres owned by three different sectors (Government, Non-Governmental Organization, and Private Sectors). All hospitals and one health centre under these sectors are currently admitting patients. These hospitals and health centres are listed in Table 2 with their respective sub-cities.
The study used a retrospective study on the under-five patients admitted to hospitals during the period of January 1, 2012 to January 1, 2014. The data were extracted from patients' cards that were confirmed by physician as a community acquired pneumonia patient. All hospitals and health centres, under Hawassa City Administration that admit patients of CAP during the period of January 1, 2012 to January 1, 2014 were taken into account. These Hospitals/Health Centres are; Hawassa Referral Hospital, Bushulo Health Centres, Adare Hospital, Kibru and Asher General Hospitals.
Hawassa city administration has a population of 349,635 people, out of which 169,685 are females and 179,950 are males. Among the population of the city administration, 223,821 people live in urban area, while the remaining 125,814 live in rural area of the administration. Population distributions across the sub-cities are given in Table 1.

Variables Considered in the Study
Variables considered at individual/patient level in the study were selected based on some previous study [9,[12][13][14][15][16] and those that are expected to be the factors/determinants of mortality of under-five from hospitalized patients of CAP.  Patients' discharge status was dependent variable in this study. It is a binary, discharge alive or dead, response which was recorded during discharge of patients from hospital. The independent variables at patient's and hospital levels which could possibly affect this variable are condensed in the following Table 3 with their description.

Multilevel Regression Models
The multilevel regression model often called hierarchical multilevel linear model (HLM) assumes that there is a hierarchical multilevel data set, with one single dependent variable that is measured at lowest level and explanatory variables at all existing levels [17]. Many kinds of data, including observational data collected in the human and biological sciences, have a hierarchical, nested, or clustered structure (Fig. 1a). Traditional linear regression techniques are not capable of taking into account data with a hierarchical multilevel structure [18]. Data with hierarchical structure violates an important assumption of linear regression model. As a result, regression coefficients are not efficient and its standard errors are understated.

Two Level Multilevel Regression Models
Multilevel regression model can be described as the extension of regression model, which adds random effects at multiple levels. For instance, in this study the interest is to know the risk of dying from CAP if one child is in a neonatal age. In this case, both variables are patient level variables; as a result, one can use simple logistic regression. Since, patients are nested with in hospital, using multilevel models allows us to divide the variance in the outcome variable at the group and individual levels. A two-level multilevel model in health area assigns patients to level 1 and hospitals to level 2.
Consider a collection of r independent variables which will be characterized by X = (X 1 , X2, . . .X r ). Then the conditional probability of the outcome of interest in a study is present (in our case death of child due to CAP) can be denoted by P ij = P(Y ij = 1|X = x). The ratio of success probability to failure probability is known as odds of success.
i.e, P ij 1 − P ij , i = 1, 2, . . . , n; j = 1, 2, · · · h Suppose that we have a sample of n independent observations of (X ij , Y ij ) for r+1 variables, where Y ij denotes the value of dichotomous variable, and Y i j is the value Fig. 1 Patients nested to hospitals and hospitals to sub-city (a) and patients nested to cross classification of sub-city with hospital at level 2 (b). a Purely nested data structure. b Cross-classified data structure. Note This figure was adapted for the data analysed here from Fielding et al. [10] of the independent variables for the ith individual (patients). Then, the logistic model can be written using odds of success as follows.
Thus, the data layout can be represented as: The relationship between success probability P i and his/her characteristics are nonlinear. Hence, by applying logit transformation, in order to have a meaningful interpretation, one would have the following logistic regression model equation.
A two level logistic regression models allow to assess variation in a dependent variable at hospital and patient level simultaneously in which patients are nested within hospitals. The response variable is denoted by Y ij : with probabilities P ij = P(Y ij = 1/X lij ), probability of dead for patient 'i' in hospital 'j', and 1-P ij , probability of discharged alive for patient i and hospital j. The standard assumption is that Y ij has a Bernoulli distribution. Let the P ij be modelled using logit link function.
The two-level model is given by: . . , β r j = β r + u r j this becomes: . . , X rij ) represents the first and the second covariates, (β 0 , β 1 , · · · , β r ) are regression coefficients, u 0j , u 1j , u 2j , · · · u kj are the random effect of model parameter at level two. It is assumed that these level two random effects distributed normally with mean 0 and variance σ 2 u .

The Empty Model
The empty two-level model for dichotomous outcome variable refers to population of groups (level-two units) and specifies the probability distribution for the groupdependent probabilities P i j , where P i j characterised by Y i j = P i j +e i j , without taking further explanatory variables into account. The logit transformed models, logit(P i j ), has a normal distribution. Then the empty model is expressed by formula: where γ 0 is the population average of the transformed probabilities and u 0 j the random deviation from this average for hospital j. The residual term associated with deviation, u 0 j is a unique effect of hospital j on the outcome variable; it is assumed to be normally and independently distributed with mean of zero and constant variance, σ 2 that is The level 2 residual variance captures the variation across hospital means. With this model, the amount of variance in the outcome variable that is attributable to within group characteristics (here, patients) and to between group differences (hospitals) can be identified. This model (2.5) does not include a separate parameter for the levelone variance [19]. This is because the level one residual variance of the dichotomous outcome variable follows directly from the success probability is indicated below.
where e i j is individual-dependent residual. Define π 0 -the probability corresponding to the average value by the γ 0 , as follows: For the logit function, the logistic transformation of γ 0 is given by: Now the variance for probability of success (death of child from CAP in-hospital j) can be given by:

Random Intercept Model
The extension of equation 2.5 above in which it incorporates explanatory variables is termed as random intercept model. The success probability is not necessarily the same for all individual in a given group. Therefore, the success probability now depends on the individual as well as on the group, and is denoted by P i j , [19]. Accordingly, the outcome is split into an expected value and a residual i.e.: The logistic random intercept model expresses the log-odds, i.e., the logit of P i j is: As usual, the deviation u 0 j are assumed to be independently distributed with zero mean and variance σ 2 . Equation (2.10) does not include a level one residual because it is an equation for the probability rather than the outcome Y i j . Thus, the level residual has already included in the equation.

Random Slope Model
Another model which is used to assess whether the slope of any of the explanatory variable has a significant variance component between the groups is called the random slope model. Assume that effect of variable X j = (X 1j , X 2j , · · · , X kj ) across the groups, and accordingly has a random slope. Equation (2.10) for the logit of the success probability is extended with the random effect u 1 j x 1i j, (19), which leads to: The characteristic that is unique to this model, compared to others discussed here, is that the random coefficients β are themselves dependent variables in a second regression equation. Of course, in RSM we have: There are now r+1 random group effects, the random intercept u oj and the random slope u l j . It is assumed that both have a zero mean. Their variances are denoted by σ 2 0 , . . . , σ 2 r respectively and their covariance is denoted by σ 0r . This can be described by the covariance matrix .

.2 Two Level Cross-classified Multilevel Models
Unlike hierarchical multilevel data with a purely nested structure, data that are crossclassified not only may be clustered into hierarchical multilevel ordered units but also may belong to more than one unit at a given level of a hierarchy [10]. Goldstein [20] labelled multilevel model for purely nested structure as hierarchical multilevel (HM) models. However, multilevel level models for cross-classified data structure is labelled as cross-classified multilevel (CCM) models.
In a cross-classified design, patients at a given hospital might be from several different resident places/sub-cities and one sub-city might have patients who attend a number of different hospitals. In this type of scenario, hospital and sub-city are considered to be cross classified factors, and cross classified random effects modelling (CCREM) should be used to analyse these appropriately.
The consequences of ignoring an important cross-classification are similar to those of ignoring an important hierarchical multilevel classification [21]. Therefore, it is very important to consider the classification at existing level (second level in this study).
In Fig. 1 the first diagram illustrates a purely nested data structure. In this structure, there are three levels; patients, hospitals and sub-city. A typical hierarchical multilevel model has patients nested in hospitals, and hospitals nested within sub-cities. In a cross-classified multilevel structure, however, patients from a given sub-city might visited different hospital and hospital might serve patients from different sub-city. Subsequently the cross classified data structure occurs at level 2.
In this example, we can take into account influences coming from two different factors: sub-cities and hospitals. This may improve the estimation of explanatory variable effects ensuring the structure is properly specified [10].
Assume one has patients nested within a cross-classification of hospital by sub-city, which is the case illustrated in Fig. 1b and P i(hosp,subc) be the probability of child to die from CAP in hospital 'hosp' and came from sub-city 'subc'. v subc(i) is the sub-city effect which are mutually independent random variables and also independent to u hosp(0, j) , which is hospital random effect. All v subc(i) are normally distributed with mean 0 and variance τ 2 v .
Note the notation used in this paper for cross classified model is adapted from Fielding et al. [10].
Like HM model, in a CCM model, before introducing explanatory variables, the variation in outcomes could be examined by looking at the differences between hospitals, between sub-cities and between individual patients after controlling for sub-cities and hospital effects [22]. This gives a basis for extending the model to identify which sub-cities, hospitals, and patients' characteristics might explain some part of these component variances.
After adding explanatory variables, the residual components of variances can be estimated, this will provide more information about the variation in outcomes [21].

Method of Parameter Estimation
Estimation of CCM models by existing maximum likelihood approaches run into important computational limitations, especially when large numbers of units are involved [21]. As a result, the models used in this paper are fitted using Markov Chain Monte Carlo (MCMC) based algorithms as implemented in the MLwiN package [23]. Starting values for the fixed parameters are estimated from simpler models using a maximum likelihood approach, penalized quasi-likelihood (PQL) in MLwiN [23].
Markov Chain Monte Carlo(MCMC) is a simulation techniques that generate random samples from a complex posterior distribution [22]. The simulated posterior distribution is then used to compute a point estimate and a confidence interval. The posterior distribution of θ is defined from Bayes' theorem as p(θ/data)/p(data/θ )p(θ ) Here p(θ ) is the prior distribution for the parameter vector θ and should represent all knowledge we have about θ prior to obtaining the data.
In multilevel models, a non-informative prior for both fixed effect and random effect parameters are recommended [18]. The default prior distribution applied in MLwiN when MCMC is used are 'flat' or 'diffuse' priors for all parameters (23).
Thus, according to Browne et al. [23] the prior for a fixed parameter, β in MLwiN is an improper uniform distribution, p(β) ∝ 1, which is equivalent to normal prior with large variance. The prior for scalar variance, σ 2 is specific as inverse gamma, p(σ 2 ) ∼ invGamma(α, β). In MLwiN, the prior applied for covariance matrices, is an inverse Wishart distribution, p( −1 ) ∼ W ishart p ( p, p, ), where p is the number of rows in the variance matrix and is the estimate for the true value of .
Consider a two-level random slope binary logistic regression in equation ( with Y i j = P i j + R i j is the outcome variables where ε = σ 2 I nj is level one residual is i.i.d. u is a level two covariance matrix. In this study, the prior distribution used for all fixed parameters is non-informative normal distribution with variance 10 4 and different values of α and β in gamma distribution for scalar variance prior was used. Consider Eq. (2.12) above, random slope model, in which each observation is considered as an outcome of a Bernoulli trial (y i j ∼ Bernoulli( p i j )), and hence for the ith patient in jth hospital: Assuming the n observations are independent, the likelihood function is given by: ) where θ stands for the joint parameters, u j , Σ ε , Ω u .
The posterior distribution is obtained by multiplying the prior distribution over all parameters by the full likelihood functions. The posterior distribution is:

Results
This study was carried out to model variability in mortality of under-five hospitalized pneumonia patients, using a retrospective data on the patients admitted between January 2012 and January 2014. The study used a sample of 305 under-5 year patients out of which 50 (16.4%) died in the hospital and the rest of them were discharged alive from hospitals. A hierarchical two level logistic regression and cross-classified logistic regression with the classification of 'sub-city' and hospitals at second level was modelled. Then these two models are compared based on their DIC. It can be seen from Table 4, more death occurred in Bushulo Health Centre which accounted about 10.2% of the total death occurred during the retrospective interval time. Likewise, in Hawassa University Referral Hospital, about 2.6% of the patients of CAP from Hawassa city were expired. In the same manner, the distribution of patients discharge status is displayed in the following Table 5.

Setting for Bayesian Inference
In this study, non-informative priors have been used which express the complete ignorance of the value of the parameter. Accordingly, as we mentioned on chapter three Sect. 3.5.3, we have used normal distribution with mean 0 and variance 10 4 (huge variance) for fixed parameters. However inverse gamma distribution with parameter α and β for single variance was used. The value of parameters for this distribution was chosen based on Bayesian deviance information criteria (BDIC). Table 6 shows the proposed parameter value for prior distribution of single variance. On Table 6, BDIC is smaller when we use α = 100 and β = 100 in each model except on cross classified empty model. Since, there was no significant change on the values of BDIC, the same prior parameter for these models were used too.
The number of iteration used in the analysis for all of the parameters (fixed and random effect) was set up to be 50,000 with 'thining' number of 1; the first 5000 was discarded(burn-in) and the rest of the chains was used to summarize the posterior distribution.

Cross-Classified Multilevel Logistic Regression Results
Now like the case of HM model, an empty model can be extended to random intercept CCM model using explanatory variables those were significantly associated with discharge status of patients. In CCM random intercept model we vary the intercept across both hospitals and sub-city.
The proportion of variation of discharge status explained by explanatory variables is 78.95% in random intercept of hierarchical multilevel logistic regression. In other words, season of diagnosed, CAP complicated disease, duration till diagnosed, age, patient refer status, patient to nurse ratio and patient to physician ration explains about 78.95% of variation in whether patients were dead or alive due to community acquired pneumonia (Tables 7, 8).

Model Diagnostic Check
In MLwiN statistical Packages, MCMC estimation methods can be easily run by taking either PQL result or MQL results as an initial point for simulation. The packages finally display MCMC summery statistics and MCMC plots which is important to test the diagnostic convergence of MCMC chain. Under MCMC summery statistics the package computes chain means (posterior means), the chain standard deviation, etc.
Raftry and Lewis diagnostic [24] is a diagnostic based on a particular quantile of the distribution. The diagnostic that is used to estimate the length of Markov chain required to estimate a particular quintile to a given accuracy. In MLwiN the diagnostic is calculated for the two quintiles (the defaults are the 2.5 and 97.5% quintiles) that will form a central interval estimate.
The dependency among the Markov chain samples are measured by the autocorrelation which is one method of model diagnostic of convergence. MCMC algorithms do not produce independent samples from the distribution. Rather, it generates samples that are auto correlated. As in time series analysis, we can compute the autocorrelation function at each lag (1+) and include the number of iterations we need to skip in order to have a non-significant autocorrelation. High correlations between long lags indicate poor mixing. Thus, low autocorrelation indicates good convergence.
All the coefficients were tested using kernel density plot, Lewis and Reftry test PACF,and ACF. Thus, the results were obtained after all were satisfied which can be attained with increasing number of iteration.

Model Selection
Model was selected based on deviance information criteria (DIC) found for each model, where model with small deviance information criteria was selected for interpretation of the results. We have used deviance information criteria to choose better model for this analysis. The DIC for cross-classified random intercept model was small as compared to the other models.  All of the variables significantly predicted discharge status of patients in both random intercept models except referred status of patients. This variable was significant in CCM random intercept model. However, it was insignificant in HM random intercept model. In addition, the explained proportion of variation (R 2 ) in outcome variable is large for CCM random intercept model.
The model comparison result suggests CCM logistic model as the best fit for discharge status of hospitalized pneumonia patients. Thus, the effects of each covariate can be interpreted using the estimated odds observed under cross-classified multilevel logistic. The odds of dying due to community acquired pneumonia are the probability that the patient died divided by the probability that the patient discharged alive.
Most of the patient level variables included in the model significantly affect the discharge status of hospitalized pneumonia patients. Season diagnosed, CAP complication, duration of onset to diagnosis and age are patient level variable that significantly affect the discharge status of hospitalized pneumonia patients. Whereas, patient to physician ratio and patient to nurse ratio are hospital level variable that significantly affect the outcome of response variable.
The value of Exp(β) for spring season is 7.54 which implies that the odds of being at risk to death during spring season were 7.54 times more likely than patients diagnosed in winter season. In CAP complication variable, there was high probability of death on patients with AGE, AFI, and malaria than patients only with CAP. The odds of these coefficients are 6.64, 7.52 and 7.15, respectively. There was high probability of dying from CAP if the patients are complicated with these diseases. There is less likely to die from CAP, if patients are postnatal and child age than those patients in neonatal age group. Accordingly, the risk of dying from CAP for postnatal age group and child age group patients was less by 91 and 92% respectively. On the result, it was found that patients who late diagnosis has high probability of being at risk than those admitted before 3 days. The odds of being at risk for patients who diagnosed after a week is 8.71 times than those patients diagnosed before three days. The odds of being at risk for patients who diagnosed with 4 to 7 days was also 4.74, which means the risk of dying from pneumonia if patients diagnosed between these days was 4.74 times than those diagnosed before three days keeping another risk factors constant.
Patient to nurse ratio, patient to physician ratio were variables significantly predicts the discharge status of patients in the chosen CCM logistic model among hospital level variables. The odds of being at risk increases by 1.39 for a unit increases of patient to nurse ratio keeping the others risk factors constant. In the same manner, the odds of being at risk increases by 1.01 for a unit increase of patient to physician ratio keeping the risk factors constant. These means number of patients increased by the number physician/nurses keeping number of physician/nurse constant. Generally, patients diagnosed during patient to nurse ratio and/or patient to physician ratio was large had high probability of dying from CAP.

Discussion
CCM logistic random intercept model revealed that more of the variation of discharge status of patients (alive/dead) was due to hospitals as compared to variation attributable to resident/sub-city. During the period of our study, the pandemic strain of community acquired pneumonia resulted in the death of 16.34% patients in the sample. Some of previous study had found approximate result, for example, 18.79% in [12] 14.7% in Cui et al. [14], 23.8% in Kuo-Tung [25].
On the result, it was also observed that patients admitted due to the case (CAP) during spring season was high in all hospitals; this is may be due to the fact that there was high prevalence of CAP in that season.
In this study discharge status of hospitalized patients due to CAP was not associated with sex as reported earlier by Teka [12] and Kuo-Tung [25]. There was also no association between treatment given by the physician and patients discharge status. Since all patients were hospitalized which means they were under the care of nurse and physicians there might be a proper usage of treatments. Thus, the non-significance of association between discharge status and treatment taken might be due to this proper usage.
In the study, also it was found that length of hospital stay (LOS) is not significantly associated with discharge status of CAP under-five patients. The median LOS was 4 days which means 50% of patients admitted due pneumonia stayed in the hospital were less than 4 days. Similar study in Spain during 2003 had found that 9% of LOS due CAP, (16).
In this model, some variation of discharge status of under-five age CAP patient due to hospitals and residence/sub-city were 19.3 and 18.9% respectively. The previous studies [12] that used some other statistical methods like binary logistic regression in which it is assumed that all patients were independent also found a significant association between the residence and discharge status of patients. Thus, this results in Type I error [26].
Finally, it was investigated that hospitals and sub-city also contributed to a variation among discharge status of CAP patients. The variation due to hospitals might be due to service given, quantity and quality of physician/nurse, in-efficient management, etc. In the same manner inaccessibility of health facility, neatness of the environment which results in other complicated disease such as malaria, high population density which results in high prevalence of pneumonia disease, could be the reason for the variation due to residence/sub-city.
In the study, it was found that patients admitted in summer and spring season had high risk of dying from CAP as compared with other seasons. According to Lieberman's [15] finding there is high incidence of CAP during spring and summer seasons. This could be the reason for the high risk of death from CAP in these seasons.
Result showed that patients who were delayed to be diagnosed had high risk of dying from CAP than those patients who diagnosed early. An appropriate usage of treatment as well as delaying to start treatment for one disease leads to high risk of dying. This is may be due to the fact that; the expansion of causing organisms increases through time in patients' body. As a result, duration of onset to diagnosis was a significant predictor of discharge status of patients. This result match with results of Menendez et al. [16] and Teka [12].
Among hospital level variables patients to nurse ratio and patients to physician ratio significantly predict discharge status of CAP patients. Patients who admitted during patients to nurse ratio was high had high risk of dying from pneumonia. Since patients in hospital are nurtured by nurses and this has a positive impact on the recovery from their illness. Fortunately, patients who admitted during ratio of patient to nurse is high has less chance to survive as it is compared to others patients. Thus, the result obtained had a complete agreement with this logic. Another important hospital level variable which again significantly predicted the risk of dying from pneumonia was patient to physician ratio. Patients who admitted during patients to nurse ratio was high had high probability of dying from CAP.
Bed occupancy rate which is one measure of quality of hospitals was insignificantly associated with discharge status of under patients. The reason for this result is not clear but it might be due to the fact that, in this study, there was no patients admitted during the beds were crowded.

Conclusions
sThis study was aimed to model the mortality of hospitalized under-five pneumonia patients and identify risk factors associated to discharge status of patients. Hierarchical multilevel and cross-classified multilevel were used. It was interesting to find out that cross-classified multilevel model was better than simple hierarchical multilevel modelling in fitting to the data and explore variability attributed to hospital and residence. The variability attributable to hospital was found to be larger than variability attributable to residence. From the analysis, it was found that age of patients, season of diagnosis, duration of onset to diagnosis, complicated disease, patient to physician ratio, patient to nurse ratio. Hospitals and place of residence/sub-city were found to contribute significantly to the variability of patient's status (dead/alive).

Recommendations
Scaling up of interventions such as health promotions on appropriate and effective treatment in home or community based care should be emphasized across all sub-cities. And also, the awareness of early diagnosis to the community should be expanded via health extensions, clinicians or any concerned body in collaboration with governmental and non-governmental organizations. Health facilities based on quality parameter of hospitals like bed occupancy rate, patient to nurse ratio, and patient to physician ratio should be improved in hospitals and be handled with effective management. Finally, researchers who are interested on similar area are recommended to choose cross-classified models for the analysis in the case where patients are nested to the classification of hospitals and resident to explore variation of response variable at all existing levels.