Pairwise joint modeling of clustered and high-dimensional outcomes with covariate missingness in pediatric pneumonia care

Multiple outcomes reflecting different aspects of routine care are a common phenomenon in health care research. A common approach of handling such outcomes is multiple univariate analyses, an approach which does not allow for answering research questions pertaining to joint inference. In this study, we sought to study associations among nine pediatric pneumonia care outcomes spanning assessment, diagnosis and treatment domains of care, while circumventing the computational challenge posed by their clustered and high-dimensional nature and incompletely recorded covariates. We analyzed data from a cluster randomized trial conducted in 12 Kenyan hospitals. There were varying degrees of missingness in the covariates of interest, and these were multiply imputed using latent normal joint modeling. We used the pairwise joint modeling strategy to fit a correlated random effects joint model for the nine outcomes. This entailed fitting 36 bivariate generalized linear mixed models and deriving inference for the joint model using pseudo-likelihood theory. We also analyzed the nine outcomes separately before and after multiple imputation. We observed joint effects of patient-, clinician- and hospital-level factors on pneumonia care indicators before and after multiple imputation of missing covariates. In both pairwise joint modeling and separate univariate analysis methods, enhanced audit and feedback improved documentation and adherence to recommended clinical guidelines over time in six and five pneumonia care indicators, respectively. Additionally, multiple imputation improved precision of parameter estimates compared to complete case analysis. The strength and direction of association among pneumonia outcomes varied within and across the three domains of pneumonia care


| CORRELATED RANDOM-EFFECTS JOINT MODEL
Let Y rij denote the r th r ¼ 1, 2, …, p ð Þ outcome for the i th i ¼ 1, 2,…,N ð Þ subject in cluster j j¼ 1, 2,…, n i ð Þ . The corresponding univariate random effects model for the r th outcome can be defined as where h À1 Á ð Þ is an appropriate link function depending on the type of outcome (i.e., whether continuous, binary, count, etc.), 12 X rij is a vector of known covariates with fixed effects β r , and Z rij is a vector of covariates with random effects b ri . The univariate random effects model can be extended to jointly model all the outcomes simultaneously, by imposing a joint multivariate distribution on the random effects. 14,15,17 Moreover, the number of random effects can vary among the outcomes of interest. Conditional on the vector of random effects b ri ð ), the outcomes are assumed to be independent 8 and the corresponding log-likelihood contribution for subject i equals l i y 1i , y 2i ,…, y pi jΘ Ã ¼ ln The vector Θ Ã contains all parameters of the full joint model (i.e., fixed parameters denoted by β Ã and covariance parameters denoted by Σ Ã ), while f ri y ri jb ri ,θ r ð Þis the density of y ri conditional on the random effects for the r th outcome on subject i. The vector of random-effects b i is assumed to follow a multivariate normal distribution with mean zero and covariance matrix D , that is, The elements D rs in D correspond to blocks of random effects variance-covariance between the r th and the s th outcomes r, s ¼ 1, 2, …,p ð Þ . For example, assuming that each outcome has a random intercept (b 0 ) and a random slope (b 1 ), then D rs is given by The elements of the variance covariance matrix D can be used to measure the strength of association between any two outcomes of interest. As mentioned earlier, the dimension of the random effects vector b i in the full joint model, increases with an in increase in the number of outcomes. This leads to computational challenges for high dimensional vectors of outcomes. 8,10,14 2.1 | The pairwise modeling approach In light of computational challenges highlighted above, Fieuws and Verbeke 14 proposed a pairwise approach within the pseudo-likelihood framework to handle high-dimensional vectors of outcomes. With a vector of p outcomes, the pairwise approach maximizes the likelihood for all Q ¼ p pÀ 1 ð Þ=2 pairwise models separately, instead of maximizing the full joint multivariate likelihood. 14,24 Precisely, this produces a so-called pseudo-likelihood (pl) of the following form: For a given pair of responses r, s ¼ 1, 2::,p ð Þ , l Y r , Y s jΘ rs ð Þdenotes the likelihood, while Θ rs is the vector of all parameters encountered in a pairwise joint model. 14 The corresponding pseudo-log likelihood function (pll) is given by where Y q and Θ q contain all the observations and parameters, respectively, in the q th response pair q ¼ 1, 2, …,Q ð Þ . All Q pair-specific parameter vectors Θ q q ¼ 1, 2, …, Q ð Þare stacked together into Θ with fixed parameters denoted by β. It is clear that if b Θ q maximizes l Y q jΘ q À Á , then b Θ, the stacked vector combining all b Θ q , maximizes pll Θ ð Þ: 24 The asymptotic distribution of b Θ is multivariate normal given by where H À1 GH À1 is a sandwich estimator and H and G are based on cluster-wise Hessians and gradients of the logpseudo-likelihood function, respectively. 8,10,18,24 The vector of all parameters in the full joint model (Θ Ã Þ and stacked vector from pairwise models (ΘÞ are not equivalent. Specifically, some parameters in Θ Ã have a single counterpart in Θ, while other elements in Θ Ã have multiple counterparts in Θ: 8 A set of fixed effects (β Ã ), for the full joint model, are obtained by averaging duplicate parameter estimates from the pairwise joint models. 8,14 This can be achieved by multiplying the stacked vector of regression parameters (βÞ with an appropriate weight matrix A as below The standard errors follow as the square root of diagonal elements of variance-covariance estimator Further details on estimation of fixed effects and corresponding standard errors are presented in the application section.

| PNEUMONIA TRIAL DATA
In this study, we analyzed routine pediatric data collected in a cluster randomized trial in 12 Kenyan hospitals between March and November 2016. 2,25 The trial's objective was to investigate the level of uptake of pediatric pneumonia treatment guidelines recommended by the World Health Organization (WHO) in 2013. 26 Details on the trial are contained in the trial report. 2 In brief, hospitals were randomly allocated to the intervention arm or control arm. Six hospitals in the intervention received an enhanced monthly audit and feedback (A&F) report on assessment, diagnosis and treatment of pneumonia cases, a bi-monthly standard A&F report assessing performance and adherence to general inpatient pediatric care guidelines at facility level. Besides A&F reports, the trial intervention package contained network intervention strategies such as peer learning among clinicians across study facilities, workshops and follow-up emails and phone calls by the trial pediatrician. On the other hand, six control hospitals received a bi-monthly standard A&F report and network intervention strategies. 2,25 During the trial period, 2299 children aged 2 to 59 months were admitted in general pediatric wards with childhood pneumonia in 12 study hospitals. However, this analysis excluded 172/2299 (7.5%) case records lacking admitting clinician's information. The remaining 2127/2299 (92.5%) patients were admitted by 378 clinicians. Of the 2127 pneumonia cases, 953 (44.8%) were admitted to six intervention hospitals. On average, there were 32 clinicians per hospital, and the number of patients per clinician ranged between 3 and 46. Data were extracted by trained data clerks from pediatric admission record (PAR) (a structured paper based medical record/form used in pediatric wards in CIN hospitals) after discharge from hospital. The data were entered into an open source data capture tool (Research Electronic Data Capture, [REDCap]) 27 using a standard operating procedure manual.

| Pediatric pneumonia care indicators
The outcomes of interest were nine pneumonia care indicators spanning three domains of care (Table 1). These are 1 = cough, 2 = difficult breathing, 3 = respiratory rate, 4 = oxygen saturation, 5 = level of consciousness measured on the 'Alert', 'Verbal response', 'response to Pain', and 'Unresponsive' (AVPU) scale, 6 = lower chest wall indrawing (signs and symptoms in the assessment domain), 7 = pneumonia severity classification (diagnosis domain), 8 = oral amoxicillin prescription to treat pneumonia, and 9 = oral amoxicillin dosage and frequency of administration (treatment domain). While these indicators were measured on different scales reflecting different aspects of care, we created a binary variable for each one of them as appropriate (Table 1). For each case record, we assessed documentation of cough and difficult breathing (primary pneumonia signs and symptoms required for identification of pneumonia cases), respiratory rate, oxygen saturation, AVPU, lower chest wall indrawing (secondary signs and symptoms required for classification of pneumonia severity). 26 For each sign and symptom, we created a binary variable with the value one representing documentation in pediatric admission record (PAR) (e.g., cough assessed at point of admission and marked in a check box as present or absent) and zero representing lack of documentation of a sign and symptom in the medical record by the admitting clinician (Table 1). In the diagnosis domain, we assessed whether clinical pneumonia diagnosis and the severity classification documented in a patient's PAR by the admitting clinician was in line with the diagnosis and the severity implied by presenting signs and symptoms. Here, we created a binary variable with value one representing correct diagnosis and severity classification and zero representing misclassification of pneumonia severity (Table 1).
In the treatment domain, we had two binary indicators, one assessing adherence to prescription guidelines and the other assessing adherence to dosing guidelines. For the prescription indicator, the value one represented prescription of oral amoxicillin to treat pediatric pneumonia as per the guidelines while zero represented deviation from ideal care (i.e., lack of evidence in a patient's medical record that oral amoxicillin was prescribed) ( Table 1).
To determine correctness of dose among oral amoxicillin recipients, we considered actual dose prescribed, patient's weight and frequency of administration as documented in a patient's medical record. We created a binary indicator with value one representing oral amoxicillin correct dosage and correct frequency of administration (i.e., 32-48 international units per Kilogram [IU/Kg] every 12 h). Inappropriate oral amoxicillin dosing was considered as: lack of documentation of actual oral amoxicillin dose prescribed, lack of documentation of patient's weight, undocumented/wrong frequency of oral amoxicillin administration, oral amoxicillin underdosing (<32 IU/Kg every 12 h) or overdosing (>48 IU/Kg every 12 h) ( Table 1).

| Covariates
In this analysis, the covariates of intertest included intervention arm, follow-up time (in months) and their interaction, hospital malaria prevalence status and pediatric admission workload. Five out of 12 hospitals were drawn from high malaria endemic regions while the remaining seven hospitals were drawn from regions with low malaria endemicity in Kenya. 28 Hospitals with less than 1000 pediatric admissions per year were categorized as low admission workload while those with 1000 or more pediatric admissions per year were categorized as high admission workload hospitals. This categorization allowed us to assess the impact of admission workload on quality of inpatient pediatric pneumonia care. This is considering that public hospitals in LMICs are often characterized by a shortage in workforce, potentially impeding delivery of health care services. [29][30][31] At clinician level, gender and cadre were considered (here cadre refers to clinician's level of training, that is, clinical officers with diploma-level training and medical officers with bachelors' degree level training). Among 295 clinicians with observed cadre, majority were clinical officer interns at 62.4% (n = 184) followed by medical officer interns at 33.4% (n = 99). Clinical officer and medical officers accounted for 2.0% (n = 6) each. Among 296 clinicians with observed gender, 43.2% (n = 128) were females.
At patient level, we considered gender, number of comorbid illnesses and age in months at point of admission. While the WHO pediatric pneumonia treatment guidelines apply for children aged 2-59 months 26 we categorized patients into two age groups, (i.e., 2-11 months and 12-59 months) in order to assess whether pneumonia care administered varied between infants and older children. This is considering that older children tend to have better outcomes compared to infants. 32 Approximately, 42.5% (903/2127) of the patients were aged between two and 11 months and 57.5% (1224/2127) were females and among 2114 patients with observed gender, 55.1% (n = 1164) were males. Regarding comorbidities, we determined the total number of clinical diagnoses documented in patient medical records. The diagnoses of interest in the comorbidity variables included malaria, malnutrition, asthma, tuberculosis (TB), rickets, anemia, diarrhea and dehydration. For each patient, we created separate binary variables for each diagnosis above with value 1 denoting the presence of a disease and 0 denoting absence of a disease. We then created a categorical variable which consisted of a count of comorbidities defined as (0 = none, 1 = one, 2 = two, 3 = three or more comorbid illnesses). The above categorization was to assess whether care among pneumonia patients varied with an increase in the number of comorbid illness. Clinically, 46.8% (995/2127) of the patients had no comorbidities, 29.8% (633/2127) had one comorbidity, 17.9% (381/2127) had two comorbidities, and 5.5% (118/2127) had at least three comorbidities. Note: AVPU a :-Alert, Verbal response, Pain response, Unresponsive, *Pneumonia diagnosis for patients with history of cough and/or difficult breathing (primary signs) in combination with signs of lower chest wall indrawing and/or respiratory rate (RR) ≥50 (≥40) for patients aged 2-11 (12-59 months), in the absence of danger any sign (inability to drink/breastfeed, cyanosis, grunting or oxygen saturation < 90% or AVPU = 'V', 'P' or 'U').

| Missingness in the trial data
Missing data occurred in patient-and clinician level covariates. Approximately, 21.9% (83/378) and 21.7% (82/378) clinicians had missing data on the gender and cadre variables respectively, while patient's gender was missing in 0.7% (17/2127) case records. An assessment of the missing data pattern revealed that nearly all clinicians with missing cadre had gender missing as well.

| Multiple imputation
Before fitting the analyses models of interest, we imputed partially observed covariates assuming a missing at random (MAR) mechanism. MI was conducted within joint modeling (JM) framework where imputation values are drawn from a multivariate normal distribution in a single step. 23,33 We used the latent normal approach to impute incomplete categorical variables of interest. 23 Multiple imputation was implemented in the jomo 34 and mitml 35 packages in R (version 3.5.4) which allow imputation of categorical variables with more than two levels in the second and higher levels of the multilevel structure. For the i th i ¼ 1, …,2127 ð Þpatient nested within the j th clinician j ¼ 1, …,378 ð Þin hospital l l¼ 1, …,12 ð Þ , we defined a two-level JM imputation model corresponding to jl are vectors of partially observed level 1 variables (patient's sex) and level two variables (clinician's sex and cadre) respectively. Predictor variables X 1 ð Þ ijl for missing patient's sex were fully observed variables (i.e., follow-up time, interacted with feedback arm, hospital admission workload and hospital malaria prevalence status, patient's age and the number of comorbid illnesses). Predictor variables X 2 ð Þ jl for missing clinician's sex and cadre at the second level of the imputation model included follow-up time interacted with feedback arm, hospital admission workload and hospital malaria prevalence status. We also included all the nine binary response variables in both levels of the imputation model. A random intercept b jl À Á was included to account for clustering at clinician level. Missing values were imputed 20 times with a burn-in of 500, and 500 updates between each imputed data set. Imputed values were assessed as appropriate while trace plots were used to assess convergence of the imputation model. 36

| Separate univariate analyses
First, we analyzed the nine outcomes separately under complete case analysis and after multiple imputation of missing covariates. For each outcome r ¼ 1, 2, …,9 ð Þ , we fitted a generalized linear mixed model defined by.
where β r1 , β r2 …, β r12 are regression parameters associated with known fixed covariates for the r th outcome. Due to relatively low numbers of clinical and medical officers, we grouped clinicians into two cadres from the initial four. That is, clinical officers (CO) combine clinical officers and clinical officer interns and medical officers (MO) combine medical officers and medical officer interns, respectively. The vector of random clinicians' intercepts b ijl is assumed to follow a normal distribution with mean zero and variance σ 2 b :

| Full multivariate joint model
To analyze the nine pneumonia outcomes jointly, a full multivariate joint model was considered: where X i denotes a vector of known covariates and β 1 , β 2 ,…, β 9 are vectors of regression parameters to be estimated for each of the nine outcomes. The random clinicians' intercepts were assumed to follow a joint zero-mean normal distribution denoted by where D, the covariance matrix of the random effects has the following structure:

| Pairwise joint modeling
To circumvent computational burden associated with model (9), we applied the pairwise approach to jointly model the probability of documentation among nine pneumonia outcomes under complete case analysis and after multiple imputation of missing covariates. We fitted 36 pairwise models where each pairwise model was defined by.
where Y rijl and Y sijl denote the r th and the s th outcomes, r ≠ s for (r,s = 1,2,..., 9) for patient i admitted by clinician j in hospital l. Each outcome occurred in eight specific pairs and we included a random clinicians' intercept in each model. For each pairwise joint model, the random effects were assumed to follow a bivariate normal distribution denoted by We fitted all pairwise joint models using the JMbayes package 37 using a server with the following specification: 40 GB memory, Intel Xeon E5-4650 (2.70GHz) processor (12 cores/24 threads), Gnu/Linux Ubuntu 14.04 OS, and R (version 3.4.4) programming language. For verification purposes, complete case analysis was also conducted in SAS version 9.4 using a SAS macro provided by. 10 Under complete case analysis, regression estimates, and standard errors were averaged across 36 pairwise models using the pseudo-likelihood approach presented in Section 2. Likewise, regression parameters were averaged across the various pairwise models for each imputed data set. Variance-covariance estimators 6 were also obtained for each imputed data set. This step resulted in 20 sets of averaged regression parameters and variance-covariance estimators respectively. Thereafter, Rubin's rules 38 were applied to obtain final estimates while accounting for within and between imputation variability. More details on the two-step procedure are as follows.

| Inference for fixed regression parameters
Each bivariate model in the m th imputed dataset had a vector of 26 regression coefficients (i.e., 13 regression coefficients for each outcome) denoted by b β qm , q ¼ 1, 2, …, 36, m ¼ 1, 2, …,20. We stacked the 36 pairwise parameter estimate vectors resulting into a column vector with 936 rows, that is, Any two pairwise joint models with a common outcome (e.g., l Y r , Y s ð Þand l Y r ,Y s 0 ð Þs ≠ s 0 Þ shared the parameters for the r th outcome. [8][9][10]24 To account for duplicate parameter estimates, we pre-multiplied b β m with an appropriate weight matrix A as follows, The weight matrix A had 117 rows (i.e., 13 regression parameters for each of the nine outcomes) and 936 columns and it was constructed such that, it averaged all duplicate parameter estimates of an outcome across the eight pairwise models in which it occurred. The resulting vector, b β Ã m was a stacked column vector of 117 parameter estimates for all nine outcomes. Each outcome had 13 regression parameters denoted by b β Ã mr . This step was repeated for all 20 imputed data sets.

| Inference for standard errors
The corresponding standard errors were obtained using the pseudo-likelihood approach introduced above. For each bivariate pair ,Y mq , q ¼ 1, 2, …, 36, in the m th imputed dataset, we estimated the variance-covariance matrix, H À1 GH À1 . Since H and G depend on the unknown parameters in Θ, 8,24 estimation proceeded as follows. N indicates the total number of subjects.
Step 1: We obtained b J mq and b K mq for each pairwise model using.
where X imq corresponds to the i th subject's contribution in the design matrix for the fixed effects, b where Z imq is the i th subject's contribution in the design matrix for random effects 24 and D mq is the variance-covariance matrix for the random effects for the q th pair in the m th imputed data set. N indicates the number of subjects.
Step2: We combined b J mq and b K mq estimated across all the 36 pairs, Step 3: We estimated H m and G m as follows. b where N is defined above.
Step 4: We obtained a variance-covariance matrix, b Σ Ã m for each imputed dataset using.

| Pooling final estimates
In the final step, we pooled the final estimates using Rubin's rules 38 for each of the 9 outcomes. This was based on the set of pairwise regression parameters and the estimated variance covariance matrices b Σ Ã m estimated in. 10 The pooled MI estimator for β is given by with variance estimator is the between imputation variance. Final MI estimates were compared to those obtained under complete case analysis.

| Wald test for joint covariates effects under complete case analysis and after multiple imputation
To test for the joint effects of covariates on the outcomes, we used a Wald-type test under complete case analysis and after multiple imputation of missing covariates. The general linear hypothesis corresponded to.
The alternative hypothesis stated that at least one of the parameters differs from zero. Under complete case analysis, the test statistic for the joint interaction term was calculated using where L 12 is a matrix of zeros and ones defined to eliminate all parameter estimates except those associated with the interaction term (i:e:, β r,12 r ¼ 1, 2, …,9 ð Þ ), b β Ã is a stacked vector of parameter estimates averaged across 36 pairwise models and b Σ Ã is the variance-covariance matrix estimated using pseudo-likelihood. Wald-test statistics for the other variables were calculated in a similar manner but adjusting the L matrix appropriately.
For imputed datasets, the joint null hypotheses were tested using linear systems of equations like those defined under complete case analysis. Nonetheless, the test statistics were calculated differently. For instance, the test statistic for the joint interaction term effect on the nine outcomes after multiple imputation was calculated using where β Ã m is a stacked vector of pooled parameter estimates for all the nine outcomes and b V m is the variancecovariance matrix based on Rubin's rules. Wald-test statistics for the other variables were calculated in a similar manner but adjusting the L matrix appropriately. In each case, the test statistic was multiplied by nine (removing the fraction in front) resulting in test statistics that were distributed according to chi-squared distribution with nine degrees of freedom. A 5% level of significance was considered in all statistical tests.

| Association among pneumonia outcomes
The strength of association among documentation of pneumonia care indicators was evaluated using the variance-covariance matrix of the random-effects. Since the covariance matrix D defined in (10) was not estimated directly at analysis stage, we constructed it using blocks of random-effects variance-covariance matrices in (12) estimated in the pairwise joint models. Under multiple imputation, we first averaged duplicate variance across 36 pairwise random intercept models for each of the 20 imputed data set. Specifically, we extracted the random intercepts variance-covariance matrix for all 36 pairwise joint models, that is, We then created an overall variance-covariance matrix D m for each imputed data set accounting for overlapping information. For example, in each imputed data set, (m ¼ 1, 2, …,20), documentation of cough occurred in the variancecovariance matrices of the first eight pairs, that is, We extracted the random intercept variances of each outcome from the pairs it occurred in and averaged them into a single random intercept variance estimate of Y r (e.g., σ 2 b 1 denoting the random intercept variance for cough). On the other hand, unique off-diagonal elements corresponding to covariance between any two outcomes were also mapped into D m . Thereafter, we averaged all the 20 D m matrices, m ¼ 1, 2, …, 20 to obtain the overall 9 Â 9 variance covariance matrix D for all the nine outcomes. We used the same procedure to construct the random-intercepts variance-covariance under complete case analysis where we averaged duplicate variances across 36 pairwise random intercept models. The strength of association between any 2 outcomes, say Y r and Y s was calculated using We performed principal component analysis (PCA) on random clinicians' intercepts variance-covariance matrices obtained under complete case analysis and after multiple imputation. This was to help visualize factor loadings among pneumonia outcomes of interest and how they correlated with one another.

| WALD-TYPE TESTS FOR JOINT COVARIATES EFFECTS
After multiple imputation of missing clinician-and patient level covariates, the Wald-type test suggested a significant joint interaction effect between the trial arm and follow-up time on documentation and adherence to recommended clinical guidelines on all the nine pediatric pneumonia outcomes of interest (P-value <0.05). Likewise, pediatric admission workload and malaria prevalence status at hospital level also exhibited significant joint effects on all the nine outcomes (Table 2). At clinician level, gender and cadre had significant joint effect on documentation and adherence to recommended pediatric pneumonia care guidelines (Table 2). At patient level, age and comorbidity had significant joint effect on documentation of all the nine outcomes. On the other hand, patient's gender did not have a significant joint effect on the outcomes of interest ( Table 2). The Wald-type test results under complete case analysis were consistent with those from imputed datasets for all the covariates. That is, all the covariates had significant joint effects on the nine outcomes except patient's gender (Table 2). Figure 1 and Supplementary Table A1-A2 present the odds ratios and their 95% confidence intervals estimated from the pairwise joint model under complete case analysis and after multiple imputation of missing covariates. Separate univariate analyses results are presented in Figure 2 and Supplementary Table A3-A4. Under pairwise joint modeling, the magnitude and direction of covariates effects varied among pneumonia outcomes of interest. Over time, documentation and adherence to recommended clinical guidelines improved in six out of nine pneumonia care indicators among children admitted to six intervention hospitals (i.e., enhanced audit and feedback arm). That is, for a unit increase in follow-up month, the change in the adjusted odds of oxygen saturation, respiratory rate, lower wall indrawing documentation (in the assessment domain), correct pneumonia diagnosis, oral amoxicillin prescription and correct dosage among patients admitted to intervention hospitals (i.e., enhanced A&F arm) were significantly more positive in comparison to the change among patients admitted to control hospitals. These observations were made under complete case analysis and after multiple imputation of missing clinician-and patient level covariates (Figure 1). Nevertheless, the estimated 95% confidence intervals estimated were consistently narrower after multiple imputation. On the other hand, there was no significant difference in the documentation of cough, difficult breathing and AVPU over time between patients admitted to six intervention hospitals (enhanced A&F arm) and patients admitted to six control hospitals (standard A&F arm).
We also observed a few instances of contrasting results. For example, after multiple imputation, the adjusted odds of AVPU documentation were significantly lower among patients admitted to hospitals with low pediatric admission workload (Figure 1 and supplementary Table A1). Under complete case analysis, however, there was no significant difference (Figure 1, supplementary Table A2). Similarly, the adjusted odds of documentation of difficult breathing and correct oral amoxicillin dose among patients admitted in low malaria prevalence hospitals were lower compared to the odds of patients admitted in high malaria hospitals. However, under complete case analysis, the difference was not statistically significant (Figure 1, supplementary Table A2). F I G U R E 1 Odds ratios (dots) and 95% confidence intervals (horizontal bars) under complete case analysis and after multiple imputation of missing covariates: Pairwise joint modeling of nine pneumonia care outcomes With regards to separate univariate analysis, the direction and magnitude of effects of most of the covariates across the nine outcomes were by and large consistent with those observed under a pairwise joint model. Additionally, it was found that documentation and adherence to recommended clinical guidelines improved over time in five out of nine pneumonia care indicators among children admitted to six hospitals in the enhanced A&F arm (intervention arm). To be specific, for a unit increase in follow-up month, the change in the adjusted odds of oxygen saturation, respiratory rate, correct pneumonia diagnosis, oral amoxicillin prescription and correct dosage among patients admitted to intervention hospitals were significantly more positive in comparison to the change among patients admitted to control hospitals. These observations were made under complete case analysis and after multiple imputation. However, multiple imputation improved precision of the estimated odds ratios compared to complete case analysis (Figure 2, Supplementary Tables A3-A4). The estimated variance among admitting clinicians (i.e., variance between random clinicians' intercepts) varied across the nine pneumonia outcomes, both under complete case analysis and after multiple imputation of missing covariates (Tables A3, A4). Tables 3 and 4 present variance-correlation matrices of random clinicians' intercepts among 9 pneumonia outcomes under complete case analysis and after multiple imputation, respectively. Generally, the magnitude of correlation estimated among outcomes was consistently larger under multiple imputation compared to complete case analysis. Moreover, the strength and direction of association among outcomes varied within and across domains of care. For instance, the strength of association between documentation of oxygen saturation and respiratory rate was somewhat high, compared to association with other indicators in the assessment domain. To be specific, correlation between oxygen saturation and respiratory rate documentation increased from 0.69 (Table 3) under complete case analysis to 0.89 (Table 4) after multiple imputation of missing covariates. In the treatment domain, prescription of oral amoxicillin and correct dosage exhibited a strong positive association with a correlation coefficient of 0.73 under complete case analysis (Table 3) and 0.80 after multiple imputation of missing covariates (Table 4).
Across domains of care, correct pneumonia diagnosis was strongly associated with prescription of oral amoxicillin and correct dosage both in the treatment domain. We also observed that documentation of oxygen saturation, respiratory rate, and lower wall chest wall indrawing, in the assessment domain were positively associated with correct pneumonia diagnosis, amoxicillin prescription and correctness of the dose. These observations were made under complete case analysis (Table 3) and after multiple imputation (Table 4). On the other hand, documentation of cough and  difficult breathing (primary pneumonia signs and symptoms) and AVPU in the assessment domain were negatively associated with documentation of other pneumonia care indicators. Under complete case analysis, a principal component analysis (PCA) on the correlation matrix of the random intercepts showed that the first and second principal components explained 57.6% and 24.6% of the variation respectively ( Figure 3, panel a). After multiple imputation, the first and second principal components explained 60.3% and 26.2% of the variation respectively ( Figure 3, panel b). Vectors of two positively correlated outcomes in the loading plots were close, forming a small angle between them (e.g., oxygen saturation and respiratory rate). On the other hand, vectors of negatively correlated outcomes (e.g., cough and treatment) were diverging forming a large angle between them. The direction of vectors for all the outcomes was consistent under complete case analysis and after multiple imputation.

| DISCUSSION
In this study we sought to estimate the joint and separate effects of covariates on nine pediatric pneumonia outcomes from a routine data set collected during a cluster randomized trial conducted in Kenyan hospitals. We also estimated the strength of association among the outcomes using a correlated random-effects joint model. 8,14 Missing data in covariate across two level of hierarchy were handled using multiple imputation.
During the trial period, documentation and adherence to recommended pediatric pneumonia guidelines by clinicians depended on individual quality of care indicators. For instance, documentation of pneumonia care indicators, that did not require a lot of cognitive effort, were highly documented (e.g., cough, difficult breathing) compared to indicators that required more cognitive effort on the part of the clinician (e.g., prescribing the right treatment in the right dosage). These variations in delivery of quality care could also be due to hospital level factors, such as lack of or broken medical devices, impeding delivery of recommended care (e.g., pulse oximeter to measure oxygen saturation).
From Wald type test, we observed significant joint effects of all covariates of interest except patient's gender and these observations were consistent between complete case analysis and after multiple imputation of missing patient and clinician level covariates. After fitting pairwise joint model, results showed that documentation and adherence to recommended clinical guidelines improved over time in six out of nine pneumonia care indicators among children admitted to six hospitals in the intervention arm. In separate analysis, documentation and adherence to recommended clinical guidelines improved over time in five out of nine pneumonia care indicators among children admitted to six hospitals in the intervention arm.
In both analysis approaches (i.e., pairwise joint modeling and separate univariate analysis), multiple imputation led to more precise estimates compared to those from complete case analysis. These observations were attributed to loss of information under complete case analysis resulting in larger standard errors hence wider 95% confidence intervals.
Further results revealed that the strength and direction of association among pneumonia outcomes varied within and across domains of care. Thus, an assumption of common random-effects among all outcomes would be too restrictive and unrealistic for pneumonia trial data analyzed in this study.
In the pairwise modeling approach, estimates obtained by averaging over several auxiliary estimates (from the various pairs) do not maximize the full multivariate likelihood. However, Fieuws and Verbeke 39 demonstrated with simulations that the loss of efficiency is small in the pairwise approach relative to a full maximum-likelihood based approach. Moreover, the averaged estimates are consistent and asymptotically normal, 8 a property which holds for imputed data sets thus, ensuring valid within imputation estimates. Validity of within imputation estimates is a prerequisite for the application of Rubin's rules which then account for between imputation variability. 38 Although we did not evaluate computational complexity explicitly, combining pairwise joint model fitting and multiple imputation comes with its computational expense as demonstrated in this study. At imputation stage, the level of complexity is compounded when missing data occur in more than one level of clustering. In such occurrences, it is paramount to account for the hierarchical structure present in the analysis model of interest in the imputation model as well. This is because incompatibility between imputation and analysis model may lead to biased estimates, underestimated cluster level variances and overestimated individual level variances. 23,33 In the current study, missing covariates were imputed using the latent normal approach within the joint model imputation framework while accounting for clustering at clinician level. Additionally, the outcomes of interest, all fully observed were included in the imputation model as auxiliary variables. Nonetheless, there is need for further research possibly through a simulation study to evaluate compatibility between imputation and substantive model or the lack thereof, in the high dimensional joint modeling context. At analysis stage, complexity stems from calculating parameters of interest (e.g., obtaining variance-covariance matrices for each imputed data set using the pseudolikelihood approach before applying Rubin's rules). Besides, constructing the overall variance covariance matrix for the random effects is not straight forward, hence the need for greater care to avoid incorrect inferences due to miscalculations. Therefore, future studies can consider developing and incorporating generic functions and packages into standard statistical software to handle missing data and other computational aspects (e.g., Wald-type tests to test for joint covariate effects after multiple imputation) more efficiently when the substantive model of interest entails joint modeling of clustered and high-dimensional vectors of outcomes.
The correlated random-effects joint model fitted using the pairwise approach has been previously used to jointly analyze clustered binary data 10 as well as continuous longitudinal outcomes. 14 However, there is essentially no example in the literature on how to account for missing covariates in a high-dimensional joint modeling context. Additionally, we extended and exemplified Wald-type tests for joint covariate effects after multiple imputation in a high-dimensional joint modeling context. To our knowledge, there are no examples in the literature demonstrating application of Waldtype tests for joint covariate effects tests in high-dimensional joint modeling after multiple imputation, hence the novelty of this study.
Besides estimating the joint effects of covariates after multiple imputation, we estimated the strength of association among quality-of-care outcomes, aspects that are largely ignored in routine pediatric care studies. In previous analysis of the trial data, for instance, diagnosis and classification of pneumonia cases was the primary outcome of interest. 2 In yet another study, pneumonia quality of care indicators were combined into a single ordinal composite outcome known as the pediatric quality of care indicator (PAQC) score. 5 Therefore, when there is need for joint inference, we recommend this study as a practical example for handling high-dimensional vector of outcomes using a pairwise fitting approach and at the same time performing multiple imputation to account for missing covariates. However, if the research question does not necessitate joint inference, then univariate mixed models as tools for analysis suffice. 14 Evidently, this study has several limitations. Firstly, we imputed missing covariates assuming a missing at random (MAR) mechanism, an assumption that cannot be verified using the observed data alone. 8,23,40 Therefore, sensitivity analysis is recommended to explore the robustness of the inferences to the MAR assumptions.
As already noted, fitting pairwise joint models on multiply imputed data sets was time intensive. Future studies may consider multiple outputation, an approach suggested by 18 as alternative to the pairwise joint modeling using a sandwich-type robust variance estimator.
In conclusion, there were significant joint effects of covariates on nine pneumonia outcomes before and after multiple imputation of missing covariates. In both pairwise joint modeling and separate univariate analysis approaches, enhanced audit and feedback improved documentation and adherence to recommended clinical guidelines over time in six and five out of nine pneumonia care outcomes of interest. Irrespective of the analysis approach, multiple imputation of missing covariates improved precision of parameter estimates compared to complete case analysis. The strength and direction of association estimated using clinicians' random intercepts estimated from the pairwise joint model varied among pneumonia outcomes within and across the three domains of pneumonia care. Across domains of care, pneumonia diagnosis was strongly correlated with oral amoxicillin prescription and dosage.