{'fraction_non_alphanumeric': 0.031077667001711443, 'fraction_numerical': 0.01347423324549727, 'mean_word_length': 4.879571953362083, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Profiling Multiple Imputation for Missing Electronic Health Records

Electronic health records (EHRs) are increasingly used for clinical and comparative effectiveness research but suffer from usability deficiencies. Motivated by health services research on diabetes care, we seek to increase the quality of EHRs by focusing on missing longitudinal glycosylated hemoglobin (A1c) values. Under the framework of multiple imputation (MI) we propose an individualized Bayesian latent profiling approach to capturing A1c measurement trajectories related to missingness. We combine MI inferences to evaluate the effect of A1c control on adverse health event incidence. We examine different missingness mechanisms and perform model diagnostics and sensitivity analysis. The proposed method is applied to EHRs of adult patients with diabetes who were medically homed in a large academic Midwestern health system between 2003 and 2013. Our approach fits flexible models with computational efficiency and provides useful insights into the clinical setting.


Glycemic testing and control
This work is motivated by critical substantive issues of glycemic testing and control for the approximately 25.2% of American seniors with diabetes (The Centers for Disease Control and Prevention 2017). Diabetes is a management intensive condition and a major cause of morbidity and mortality. The American Diabetes Association (ADA) recommends glycosylated hemoglobin (A1c) testing at least two times a year in patients who are meeting treatment goals and have stable glycemic control, and quarterly in patients whose therapy has changed or who are not meeting glycemic goals (ADA 2018a). Evidence-based guidelines are most applicable to relatively healthy patients with diabetes.
For complex patients with diabetes, those aged 65 years or over, or those with comorbid conditions, adhering tightly to established guidelines might diminish the quality of care and lead to adverse outcomes. Complex patients may be at high risk of complications-such as hypoglycemic coma, seizures, falls, fractures, motor vehicle accidents, cardiovascular events, stroke, or acute renal failure-and also less likely to experience benefits of tight glycemic control, e.g., shown by the findings of the ACCORD and ADVANCE trials (Dluhy and McMahon 2008;Martin et al. 2006). Less stringent goals may be appropriate for complex patients with comorbid conditions, or long-standing diabetes in whom the tight control goal is difficult to achieve (ADA 2018b). Therefore, complex patients would be best suited by individualized treatment plans, but little clinical evidence exists to develop these plans.
To enhance findings from clinical trials, electronic health records (EHRs) are increasingly proposed as a data source for clinical and comparative effectiveness research on improving healthcare (Agency for Healthcare Research and Quality 2013;Cebul et al. 2011). The EHR is a longitudinal record of patient medical information that is maintained by encounters in any care delivery setting, and may include all of the key administrative clinical data relevant to patient care, such as laboratory data (e.g., A1c values), demographics, progress notes, problems, medications, vital signs, past medical history, immunizations, and radiology reports. EHR automates access to information and has the potential to streamline the clinician's workflow and support other care-related activities directly or indirectly through various interfaces, including evidence-based decision support, quality management, and outcomes reporting (U.S. Center for Medicare & Medicaid Services 2012).
We would like to use EHRs to examine the role of patient complexity in the recommendation for glycemic control. However, EHRs suffer from usability deficiencies (Healthcare Information and Management Systems Society 2009). Collected in an unscheduled fashion, typically only when the patient seeks care or the physician orders care, EHRs suffer from a large amount of intermittent missing data.
In this article, we focus on intermittently missing A1c values in EHRs. Multiple factors could influence the degree of missing A1c values for individual patients. Physician nonadherence to guidelines for lab testing will cause missing values. The patient's health conditions potentially affect the propensity to test, leading to informative presence patterns that may cause systematic estimation bias. Incorporation of all available information including patient health status is critical to appropriately handle missing A1c values from EHRs.
We propose a Bayesian profiling approach to incorporate patient characteristics and infer latent groups of A1c measurement trajectories. The measurement trajectories examine both the A1c collection or presence patterns and collected values across time. The Bayesian profiling approach combines with multiple imputation (MI, Rubin (1987)) to produce complete EHR datasets for general analysis purpose. As an illustration of the MI inference, we evaluate the association between A1c levels and the incidence of any acute health events, such as hospitalization, emergency room (ER) visit or death.

Statistical literature on intermittently missing data
Missing A1c values in EHRs present statistical methodology challenges. Intermittent missing data in large-scale, unbalanced longitudinal studies, where subjects reappear after one or more missed visits, call for new imputation approaches especially when the missing values are potentially nonignorable (Little and Rubin 2002) and the observed cases are sparse. Simple methods (e.g., complete case analysis or last observation carried forward) will distort the relationship and are not recommended by National Research Council (2010). Standard missing data methods in longitudinal studies focus on a common set of pre-specified and monotone missing times, where a measurement being missing implies that all followup measurements are also missing or dropped from the analysis (e.g., reviews in Ibrahim and Molenberghs (2009)). Likelihood-based, weighting or extrapolation factorization approaches for missing data mainly apply to monotone missingness. Utilization of the information collected after subjects reappear will potentially correct bias and increase estimation efficiency, and the challenge lies in how to appropriately model sparse data structures with non-monotone missing patterns.
We seek a flexible imputation engine under MI to predict intermittent missing values. MI separates imputation and analysis into two steps and propagates the uncertainty due to missing data. Various MI software packages have been developed assuming data are missing at random (MAR) based on either joint multivariate normal distributions, e.g., PROC MI (SAS Institute Inc. 2017), Amelia (King et al. 2001) and norm (Schafer 1997), or a sequence of fully conditional distributions, e.g., IVEware (Raghunathan et al. 2001), mice (Van Buuren and Oudshoorn 1999), and mi (Gelman et al. 2015). Multilevel models are embedded with MI to handle correlated data, such as the packages REALCOM-IMPUTE (Carpenter and Kenward 2013) and pan (Schafer 2016). However, existing MI software cannot handle high-dimensional data that are subject to high proportions of intermittent missingness in large-scale longitudinal studies, nor nonignorable missing values due to the default MAR assumption.
With not missing at random (NMAR), the general framework needs to jointly model the incomplete variables (i.e., the A1c values) and the missingness indicators (i.e., whether A1c values are present on any given occasion), where the imputation model for the outcome has to be combined with the response propensity model for the presence patterns via a selection model, pattern mixture model, or shared parameter model (Little 1995;Wu and Carroll 1988). Previous work introduces strong and untestable assumptions for parameter identification, and the computation is non-trivial (Ibrahim et al. 2005). The longitudinal A1c collection history in EHRs results in sparse observations across a large number of different presence patterns and calls for flexible modeling strategies that account for the time dependency, borrow information across patterns and adjust for patient heterogeneity.
We jointly model the A1c presence patterns and the lab values taking into account patient characteristics. To achieve parameter identification and dimension reduction, we introduce latent profiles and develop a Bayesian profiling multiple imputation (BPMI) approach for missing A1c data in longitudinal EHRs. The novel latent profiling can handle ignorable and nonignorable missingness.
We assume that latent profiles are primarily determined by the trajectories of the longitudinal A1c measurements, including A1c values and potentially the presence patterns. The latent profiling includes patient characteristics-such as socio-demographics, healthcare utilization measures, medications, complexity, comorbidity and complication-as covariates that affect the profile assignment. This is an improvement over previous work that required joint modeling of the outcome variables, missingness indicators and covariates, leading to the extra need for Monte Carlo integration or approximate inferences and poor performances.
The proposal shares a similar decomposition with Lin et al. (2004) but differs by the MI framework under the Bayesian paradigm. Besides avoiding the extra computation step of Monte Carlo integration, the posterior updating of BPMI has efficiency gains via Gibbs sampling and generates completed datasets for general analysis purposes.
The paper is organized as follows. Section 2 describes the EHR structure and content. The developed BPMI is presented in detail in Section 3 and applied to the EHRs in Section 4 which presents our statistical and substantive findings. Section 5 concludes the contribution and discusses future extensions.

EHR description
We collect the EHR data of adult patients with diabetes who were medically homed at a large academic Midwestern health system between 2003 and 2013. Patients are eligible for inclusion if they have diabetes defined by a qualifying International Classification of Diseases-9 code for diabetes. Because A1c reflects mean blood glucose over the preceding three months and quarterly testing is recommended for complex patients, we construct a longitudinal dataset with one measurement per patient per three-month period (patientquarter). When there are multiple A1c values available during the three months, we use the average of the measurements restricted to be before the first acute health event date if any. We treat the first four quarters as the baseline period and use the average of collected A1c values during the first four quarters as the baseline A1c value. Patients are included if they have baseline A1c values and at least one available followup A1c measurement. Patients' various enrollment dates and lengths of followup result in an unbalanced data structure.
We identified 7372 patients with 113761 quarters after baseline and only 57285 available observations. The intermittent missingness proportions across patients are as high as 0.97, with a median of 0.50. The baseline A1c values of 7372 patients center around mean 6.93% and range from 3.55% to 17.1%, showing heterogeneity in glycemic control.
We randomly select 30 patients and present the spaghetti plot of their collected A1c values in Figure 2.1. The time trends vary across patients with different ranges, and the measurements of the same patients tend to be correlated. We consider mixed effects models accounting for within-patient correlation with random intercepts and slopes with respect to time or functions of time. While the A1c measurement trajectories present patient heterogeneity, subgroups of patients could share similar profiles across time. We assume each subgroup to have different location and scale for the time trends, which results in a mixture distribution overall. EHRs collect rich patient characteristics that could be predictive of the A1c trajectories and likelihood of missingness or presence patterns. The patient-quarter structure creates both time-varying and time-invariant variables. Time-varying variables include age and the count of physician visits that occurred within each quarter. On average patients have 2.39 physician visits every three months with standard deviation (sd) 2.03 across patient-quarters.
The baseline characteristics include socio-demographics, healthcare utilization measures, medications, complexity, comorbid conditions (Elixhauser et al. 1998), and diabetes complications, all of which are carried forward (e.g., chronic conditions will be assumed as present from diagnosis onward). Patient complexity is defined as the presence of congestive heart failure (CHF) and/or Ischemic heart disease (IHD). The comorbid patterns are characterized by Hierarchical Condition Categories (HCC) risk scores. The number of adverse events, such as hospitalization and ER visits, and the number of ER visits that did not lead to hospitalization during the first four quarters, serve as proxy measures of healthcare utilization.
The latent profiling allocation takes into account variables that are predictive of the A1c values or presence patterns into the model. We use the observed A1c values ignoring the patient-quarter structure and fit an ordinary linear regression to select the variables that are highly correlated with A1c. Then we fit a logistic regression model to predict the missingness at each quarter and select the variables that are predictive of the visit pattern. The union of the two sets of selected variables results in the covariate list in Table 2.1. The patient characteristics are potential confounders and moderators of the tight control effect on the health outcomes.

Bayesian profiling multiple imputation
Denote the individual record for patient i by {X i0 , X ij , Y ij , R ij , j ∈ [1, T i ]} for the total number T i of tracked quarters during the followup, where X i0 are the time-invariant covariates, X ij 's are time-varying covariates, and Y ij is the variable to be imputed, A1c values over time. For brevity, we will refer to Y ij as the longitudinal outcomes. Assume that only Y ij 's are subject to missing values, and let R ij be a time-varying binary indicator for its presence, . . , n total number of patients. Assume individuals fall into latent classes C i = 1, . . . , L, where L is the total number of classes. We assume L is a finite and positive known integer. The selection of L and the case of L being an unknown random variable that induces additional uncertainty will be discussed later. We assume that the latent class structure is primarily determined by the trajectories of the longitudinal outcome and potentially its presence patterns, and that the allocation probabilities are affected by the covariates.
The conditional distribution of the observed measurements can be expressed as where f (·) denotes the distribution.
If R i only depends on fully observed (X i0 , X i ), i.e., MAR, the observed data can provide We call this marginal profiling because the latent class structure influences only the marginal distribution of Y .
of its marginal dependence on the observed longitudinal outcomes. The MAR assumption can be relaxed to that is, conditional MAR, where the missingness is independent of the outcome given the latent classes and covariates. However, the unconditional missingness is nonignorable as the latent structure affects both the longitudinal outcome and presence patterns. Joint modeling of (Y i , R i ) is then necessary and will be identified due to the conditional independence assumption given the latent classes. We call this approach joint profiling. Joint profiling is identical to marginal profiling when the estimated parameters in Model (3.1) do not change across profiles, so R i is conditionally independent of C i given (X i0 , X i ). Hence, the joint profiling approach is general with marginal profiling as a special case and able to model both ignorable and nonignorable missing data. We will consider both the marginal and joint profiling for imputation and perform sensitivity analysis of the inference on the effect of A1c levels on acute health outcomes.
The latent classes are at the individual level and thus time invariant. Further, we assume the latent class allocation depends only on time-invariant variables f (C i = l|X i0 , X i ) = f (C i = l|X i0 ). For computational and interpretational convenience, assume the outcome trajectory is independent of time-invariant covariates given the latent class . This assumes the patient characteristics at baseline are represented by the latent profiling.

Marginal profiling
Under MAR we construct latent profile models based on the collected outcomes, i.e., marginal profiling, which is estimated from the pattern of trajectories of the longitudinal outcomes.
Here Y obs i denotes the observed A1c measurement for individual i. The factorization needs two modeling components: the outcome model and the latent profile model. This model specification in (3.2) assumes: 1) the presence pattern is independent of the A1c values and the latent profiles given the covariates R ij ⊥ ⊥ (Y ij , C i )|(X i0 , X ij ); 2) the A1c values are independent of the presence patterns and the baseline characteristics given the latent profiles and time-varying covariates Y ij ⊥ ⊥ (R ij , X i0 )|(C i , X ij ); and 3) the latent profiles are independent of the time-varying covariates given the baseline characteristics C i ⊥ ⊥ X ij |X i0 . Model (3.2) can be generalized by assuming the outcome depends on the baseline characteristics that are assumed to have been captured by the latent profiling in the current specification.

(3.3)
For posterior computation we develop a Gibbs sampler by introducing Pólya-Gamma (PG) distributed variables (Polson et al. 2013 Conditional on the PG variables, the posterior distribution of η l will be conjugate with normal prior distributions. The resulting Gibbs sampler improves the posterior fitting without the stickiness that is common with rejection sampling methods. Let c il = I(C i = l) be the latent class indicator and introduce the normal prior distribution N(b l , B l ) on the coefficients η l , the conditional posterior distribution of η l given w is multivariate normal, π(η l | −) ∼ N(m * l , S * l ). Here S * l = (V T Ω l V + B −1 l ) −1 , V is the design matrix with each row X T i0 , Ω l = diag(w 1l , . . . , w nl ), and m * l = S * l (B −1 l b l + S T m l ), where m l is a vector in R n , and the ith component is m For the collected longitudinal and continuous A1c measurements, we assume a linear mixed effects model with profile-specific coefficients and variances (β * C , σ 2 C ). Let D i denote a T i × p covariate matrix, with associated p-vector of coefficients β. The jth row of D i , denoted by D ij , is then a p-vector of covariate values measured at time j. Covariates for profile-specific effects β * C i and for individual-specific random effects b i are denoted by the T i × q matrix D * i and the T i × r matrix D * * i , respectively, both with structures similar to D i . There may be overlap of the covariates in D i , D * i and D * * i , including main effects and high-order interactions of X i . Deterministic functions of time, for example, basis spline functions of time in our EHR application, can be included in the covariates.
which is a linear mixed effects model with a mixture of location and scale parameters varying across profiles. The profile-specific coefficients β * C i can capture the differential trajectories of A1c measurements with time-varying covariates in D * i . Model (3.4) is a location and scale mixture model that is flexible to capture non-Gaussian distributions as shown in Figure 2.1. Here Σ r×r is the covariance-variance matrix of the individual-specific random effects b i , which could be simplified as a diagonal matrix if the components of b i are treated as independent or a scalar if b i only represents random intercepts.
We assign weakly informative and conjugate prior distributions to the parameters (Gelman et al. 2008). The inference is not sensitive to the hyper-parameter values, as expected with large sample sizes. The prior specification and full conditional posterior distributions as efficient Gibbs sampler are presented in the Appendix A.1. After convergence, based on the posterior samples of the parameters, we can impute the missing A1c values and disseminate completed datasets.

Joint profiling
With NMAR, we jointly model the A1c values and presence patterns given the covariates to obtain the joint profiling structure where Y i includes both observed and missing measurements. The joint model has three components: the latent profile model, the longitudinal outcome model and the longitudinal response/presence propensity model. This NMAR assumption move beyond MAR by joint profiling to capture the dependence between the outcome and presence patterns that jointly determine the latent profiles. The presence indicator is conditionally independent of the outcome given the latent profiles and covariates, i.e., conditional MAR: R ij ⊥ ⊥ (Y ij , X i0 )|(C i , X ij ). Conditional MAR makes models identifiable, yet flexible for dimension reduction and interpretation.
We consider a generalized linear mixed effects with a logit link for the longitudinal presence indicator R ij , conditional on latent classes and random effects e i ∼ N (0, E) with E as a covariance-variance matrix.
Where the covariates (D ij , D * ij , D * * ij ) could overlap or be different from those selected in Model (3.4) as subsets of the main effects and high-order interactions in X ij .
We develop the Gibbs sampler with two sets of introduced Pólya-Gamma distributed variables for posterior computation with models (3.3), (3.4) and (3.6) under weakly informative and conjugate prior distributions. Different from marginal profiling, imputation of missing A1c values is nested in the iterative process, where parameter estimation and data augmentation are simultaneously implemented. The posterior computation with prior specification and imputation are presented with details in Appendix A.2.
To facilitate interpretation of the potential profiling structure, we assume the number of latent patterns is fixed with a finite value. This can be extended to allow a random variable under the nonparametric Bayesian modeling framework such as the dependent probit stick-breaking process (Rodríguez and Dunson 2011). However, discrete covariates with low prevalence tend to cause separation problems with a large number of latent profiles, as illustrated in our application study. In this paper, we fix L to be a modest or small integer, and we use diagnostic tools for model selection in Section 4.1.

Imputation and Inference with EHRs
We apply the BPMI to impute missing A1c values in the EHRs. The collected A1c values range from 3.3% to 17.3% with mean 7.1% and are right skewed with potentially multiple modes. We use the location and scale mixture distributions as a flexible strategy to capture the irregular distribution. In the profile allocation model (3.3), the time-invariant covariates X i0 are the baseline characteristics listed in Table 2.1. In the outcome and presence propensity models (3.4) and (3.6), the covariates in D ij with fixed effects include the constant ones and the time-varying variables, age and the count of physician visits. The covariates in D * ij represents the functions of time in years with slopes changing across classes. We use the basis spline functions of time with (1%, 15%, 20%, 50%, 75%, 90%) quantiles as six knots and cubic terms as a piecewise polynomial regression, and the nine spline functions used as covariates in D * ij with profile-specific coefficients as shown in l 's whose starting points are the same across classes. The initial values of η's are 0, representing initial equal probability of class allocation. The scale parameters start at 0.1, and the random effects are initially drawn from normal distributions with mean 0 and the initial scales.
We implement posterior computation for both marginal and joint profiling. The MCMC algorithms efficiently achieve convergence with Gibbs samplers. We run the MCMC chains long enough to obtain 100 multiply imputed datasets and 500 replicated datasets to check imputation plausibility, i.e., the model capability to predict observed outcomes. For classification, we follow Goodman (2007) with a hard partitioning and keep one random assignment based on a random draw from the component-specific probabilities across iterations. Our analysis shows that the results are not sensitive to the classification rules.

Model diagnostics and comparison
We use the Bayesian information criterion (BIC) and the log-pseudo marginal likelihood (LPML) to select the number of classes, L. For BIC, we use the posterior mean estimates of related parameters and obtain the likelihood conditional on the class allocation. LPML is the  sum of logarithms of conditional predictive ordinates and stimated using posterior samples of parameters, θ (t) , t = 1, . . . , T , and LPML=[1/T T t=1 1/f (Y i |θ (t) )] −1 (Gelfand and Dey 1994). This is an approximation of the leave-one-out cross-validation using importance sampling. We use f (Y i |θ (t) ) for the marginal profiling and f (Y i , R i |θ (t) ) for the joint profiling.
We restrict L to be modest or small to avoid computational problems due to separation of the low-frequency predictors, such as dementia and renal failure. Table 4.1 presents the BIC and LPML values for models with different L values. The model with the smallest BIC and largest LPML will be favored. As a tradeoff, we select the model with L = 3 classes both for marginal and joint profiling.
To assess imputation plausibility, we perform a posterior predictive check and generate replicated datasets that are predicted values of observations based on the posterior estimates of model parameters after convergence (He et al. 2010;Meng 1994). Let {R (1) , . . . , R (T 0 ) } be the collection of the T 0 replicated data sets. We then compare statistics of interest in each replicated dataset to those in the observed dataset D. Suppose that S is some statistic of interest, let S R (t) and S D be the values of S computed from R (t) and D, respectively. The quantities S include the mean, 2.5% percentile and 97.5% percentile of A1c levels for every patient and for every quarter. For the patient-level summaries, we look at the two-sided posterior predictive probabilities as diagnostic tools defined as When the value ppp is small, for example, less than 5%, this suggests that the replicated data sets systematically differ from the observed data set, with respect to that statistic.
With not small ppp values, the imputation model preserves the observed characteristics of that statistic. For the 7372 patient-level A1c summaries, the marginal profiling model yields 52 ppp values that are below 1% for the average, and the number of below 1% ppp values for the 2.5% percentile and for the 97.5% percentile are 238 and 375. The number of ppp values that are below 1% for these three statistics model under joint profiling is 66, 464 and 678, respectively. The posterior predictive check does not raise the flag of lack of model fit. However, the performance of joint profiling is inferior in preserving the observed mean values across patients in comparison with marginal profiling. For the quarter-level A1c summaries, we compare the predictions with the empirical values under two profiling approaches. The 95% predictive credible intervals greatly overlap with the empirical 95% CIs. Figure 4.2 shows that the posterior predictive mean estimates are generally close to the empirical mean values, except for two observations, and both approaches perform competitively and similarly. We check the parameter estimates in Model 3.1, and find that the coefficients do not change across profiles. This explains the similar performance of joint profiling and marginal profiling.
To evaluate the robustness against model misspecification for the A1c trajectory, we use the logarithm transformed A1c value as the outcome variable and fit a linear regression model with time in years in D * ij , the slopes of which change across classes. Based on BIC and LPML, we select L = 4 classes for marginal and joint profiling under the linearity assumption. The posterior predictive check shows that imputation performances under both marginal and joint profiling are worse than those with spline functions, with larger bias and variability. Hence, we choose the models with spline functions.

Sensitivity analysis
We aim to address two questions: (1) How using imputed data affects inference on how tight glycemic control is associated with acute diabetes outcomes. (2) How results differ before and after imputation. Since marginal profiling and joint profiling present similar results, we use the results under marginal profiling as an illustration. The overall incidence of any  Figure 4.3: Comparison between BPMI inferences with complete data analysis for the relationship between A1C levels and acute health event incidence.
acute health events is 16% across the patient-quarters. We perform inferences after MI using the combining rules (Rubin 1987) where the variances account for the uncertainty due to imputation. We consider a logistic regression model with a seven-category discrete A1c indicators (< 6%-reference level, 6 ∼ 6.5%, 6.5 ∼ 7%, 7 ∼ 7.5%, 7.5 ∼ 8%, 8 ∼ 9%, ≥ 9%). We fit marginal models using generalized estimating equations (GEE) with exchangeable working correlation to account for the correlation between repeated measures of the same patient. We also fit generalized linear models with maximum likelihood estimation (MLE), and the conclusions are similar. However, MLE suffers from computational problems with large-scale datasets. Figure 4.3 presents the mean and 95% CIs of the regression coefficients. Analysis based on only the observed cases and the BPMI inferences yield different conclusions with different degrees of uncertainty. Generally, the profiling analysis with imputed data shows fewer benefits of A1c values between 6% and 9% than those under the complete case analysis (CCA), in comparison with A1c values below 6%.
Hence, profiling displays a less pronounceable U-shape across A1c values. Specifically, CCA indicates that patients with A1c values between 6% and 8% consistently have significantly lower risks than those with even tighter A1c control (below 6%). However, this is not the case for marginal profiling for A1c values between 7% and 8%. Interestingly, both approaches find that risk levels are not distinguishable between patients with A1c values between 8% and 9% and those with extremely low values (below 6%), and the acute event risk is significantly higher for patients with high A1c values (above 9%).
The BPMI approach supports that the optimal A1c level range for patients with diabetes is 6% to 7%, while CCA shows it to be 6% to 8%. The evidence from BPMI confirms the finding in Kur et al. (2019). Our analysis shows that inferences are sensitive to the choice of method to handle intermittently missing values. The imputed A1c measurements tend to reduce the estimated negative effect of tight A1c control. Figure 4.3 shows that there is a nonlinear or U-shaped relationship between A1c values and the risk of adverse outcomes. Future work would be to study the heterogeneous effects of patent complexity in the U-shaped relationship as stringent A1c goals tend to harm patients with complications. We will further discuss the implications of including adverse event outcomes into the A1c imputations in Section 5.
In this paper we focus on the missing A1c values, where the MAR and NMAR assumptions refer to the model with A1c as the outcome. The adverse event indicators are fully observed. The CCA will exclude the occasions with missing A1c values and then the adverse event outcomes together with other corresponding time-varying information. This could cause different results between BPMI and CCA if the missingness of A1c induces the missingness of adverse events in the observed model that may be MAR or NMAR. However, this assumption cannot be verified and requires to account for adverse events in the A1c imputation model that is not the goal of this paper. Hence we assume the missingness of A1c does not depend on the adverse event indicators, which is reasonable since the baseline adverse event indicator and related patient characteristics are included as covariates. This could explain why the MI estimates differ from CCA may lie in the covariates that BPMI includes. The inferences look at the overall relationship while the BPMI uses rich patient information to impute missing A1c values.

Discussion
We propose a Bayesian latent profiling approach to handling longitudinal intermittent missing A1c values in EHR for patients with diabetes. The latent profiles are primarily determined by the A1c measurement trajectories and potentially the presence patterns, and the patient characteristics affect the assignment tendencies. Different assumptions on the missing data mechanisms result in marginal and joint profiling approaches with location and scale parameters varying across latent profiles. Model performance is evaluated via diagnostic tools, posterior predictive check and sensitivity analyses. The Bayesian profiling model is able to preserve the observed data structure. The inferences after MI are different from CCA. The BPMI model with spline functions presents the best fitting and the most plausible imputation. The performance of marginal profiling and joint profiling is similar in the EHR application, where the MAR assumption tends to hold.
The BPMI has the appealing property for dimension reduction with a sparse visit structure and a large number of covariates and performs robust against model misspecification. It should be noted that currently available software for fully conditional imputation, such as mice and mi, failed to converge in our application, and produced a distorted observed data structure. The shared parameter models are subject to model misspecification and can cause estimation bias. The BPMI uses mixture model with different location and scale parameters that are flexible to capture different distribution structures. For the scenarios where the alternative methods work well, the BPMI yields competitive performances based on our simulation results that are not presented here and will be shared upon request.
Several extensions of BPMI may be of interest. First, incorporation of time-varying health outcomes into modeling can provide additional information to impute the A1c levels. The all-cause health outcome includes a combination of any hospitalization, any ER visit, or death during each quarter. Figure 4.3 demonstrates a nonlinear relationship between A1c values and the outcome. Hospitalizations or ER visits due to possible drug-related morbidity may also be a significant risk for patients with diabetes under tight control. However, the causal pathway makes sequential models more appealing than the joint modeling approach. We would like to generate completed EHRs for general analysis purpose, so incorporating the adverse event outcomes is not practical for comparative effectiveness research. The relationship between A1c and adverse events (i.e. the U-shape at low A1c) could be different for people with complex disease, i.e., the U-shape could be stronger for these people. Future work will be on exploring the heterogeneous effect of patient complexity or across latent profiles in the relationship and then A1c goals.
Other time-varying variables include blood pressure (BP), elevated low-density lipoprotein (LDL), number of primary caregiver office visits and Body Mass Index, shown as associated with tight A1c control (Jackson et al. 2006;McFarlane et al. 2002;Niefeld et al. 2003). However, these variables are also subject to missingness, and can have higher missingness percentages than A1c, such as BP and LDL. With more information available from EHR, variable selection needs further investigation, in particular being integrated with the Bayesian profiling approach as one systematic process.
10. Update e i from [e i |−] ∼ N (m * e i , S * e i ), where m * e i = S * e i ((B * * i ) k e i ), B * * i = (B * i1 ) , ..., (B * iT i ) and k e i is a vector of R T i as the following: where Ω e i = diag(w * i1 , ..., w * iT i ).
11. Update E from IW (n + ν e , n i=1 e i e i + Σ e ).