This retrospective study was based on anonymized data of 189,520 patients with T2DM from SingHealth Diabetes Registry between January 2013 and December 2018. The study was approved by the SingHealth Centralized Institutional Review Board with a waiver of informed consent granted. The board deemed that further ethical deliberation was not required as the study involves analysis of an anonymized dataset. All methods performed in this study were in accordance with the relevant guidelines and regulations. The dataset was split into training data (80% with 152,527 patients) for treatment recommendation models, including three types of treatments (antiglycemic, antihypertensive, and lipid-lowing treatment), and test data (20% of data with 36,993 patients) for evaluating the effectiveness of the three treatment recommendation models.

The EMR data for each patient included demographic information, medical history, physical measurements, laboratory data, and physicians’ prescriptions. Demographic information included age, gender, ethnicity, smoking, and others. Medical history included comorbidities, vascular complications, hospital admissions, emergency department visits, and outpatient visits. Physical measurements included systolic blood pressure (SBP), diastolic blood pressure (DBP), heart rate, weight, height, and BMI. Laboratory data included glycated hemoglobin A_1c (HbA_1c), low-density lipoprotein cholesterol (LDL-c), and fasting plasma glucose. For physicians’ prescriptions, we considered only antiglycemic, antihypertensive, and lipid-lowering drugs and their dosages.

As is the case for all EMR data, our dataset contains errors and missing data. The rate of missingness was generally low, with higher rates for variables under the categories of physical measurements and laboratory data. We handled the errors and missing data using the following strategy. During preprocessing, errors were treated as missing values. For missing physical measurements and laboratory data, we substituted the missing data with the value from the closest preceding data point of the same patient within a 1-year time frame. If data were still missing, we proceeded to impute the missing data using the median of the observed values for that variable for all patients without missing data.

Two types of clinical outcomes were analyzed: short-term and long-term outcomes. Short-term outcomes were evaluated at the patient-visit level, including blood glucose control, blood pressure control, blood lipids control, and hypoglycemia-related admissions. For the long-term outcomes, we evaluated the occurrences of 5 diabetes complications and death in up to 6 years at the patient level, including myocardial infarction, heart failure, stroke (including ischemic and hemorrhagic strokes), diabetic nephropathy, other microvascular complications (diabetic neuropathy, diabetic eye complications, diabetic foot/peripheral angiopathy), and death.

The treatment recommendation models were based on the patient’s clinical information from visits as input to recommend three types of treatments as output: antiglycemic, antihypertensive, and lipid-lowering medications. The input clinical information of a patient contains demographic information, lab data, physical measurements, medical history, and prescriptions currently in use. We utilized three models to recommend the three types of medications, and then combined the output of the three models into a comprehensive treatment recommendation.

Figure 1 illustrates the treatment recommendation approach by combining a knowledge-driven model and a data-driven model. The knowledge-driven model was developed based on the clinical guidelines and expert experiences on managing T2DM [29-34]. For the data-driven model, RL was used to learn the policy of treatment recommendation from real-world data that optimizes a predefined long-term goal via trial and error [35-37]. When integrating these two types of models, the knowledge-driven model was first applied to select the candidate medications, and the data-driven model was used to rank the candidates by the expected clinical outcomes.

We trained the RL model on a set of time-varying data consisting of state s_t (clinical data of the current visit), action a_t (treatment), and reward r_t score (based on clinical outcome). The ultimate goal of RL is to learn a policy π, which for any given state, one can select the action that maximizes cumulative future rewards.

The DQN [38] is a type of RL method that has been recently utilized to solve clinical decision problems with continuous state variables [13,39,40]. Referring to the previous DQN work of sepsis treatment in the intensive care unit [13], we applied deep neural networks to calculate the action-value function Q that estimates the cumulative rewards for each treatment action at the current visit state. To train the DQN model, two neural networks with the same architecture were used: an evaluation network Q(·) and a target network . The evaluation network was used to obtain optimal action max_a Q(s_t, a, θ) and was trained by the loss function L(Q′, Q). The target network was used to estimate the expected action-value Q′ to calculate the loss function L and updated its parameters by slowly tracking the parameters of the evaluation network θ every training iteration: ←τ·θ+θ (1 – τ)· with update parameter τ<1. The loss function L(Q′, Q) is defined in Equation (1):

where the expected Q′ is:

and r_reg is the maximum reward of all r_t that is used in the regularization term to penalize an inexpertly large Q value. Here, we used a double-DQN [36] architecture that calculates in Equation (2) through action from the evaluation network a′=argmax_a [Q(s_t+1, a, θ)] instead of max. Double-DQN leads to a more stable learning target and low-variance action-value estimates. Moreover, we used dueling-DQN [37] that adds a dueling architecture in the network to separate the output of the last hidden layer into two streams to learn state values and state-independent action advantages, respectively. We also used a prioritized experience replay [41] method to speed up the training approach. Each training batch was sampled from the training data according to the importance, which was measured by the samples’ temporal-difference error. The complete training procedure of our DQN model can be found in Algorithm S1 in Multimedia Appendix 1.

To present the specific design of the RL-based data-driven model, we take antiglycemic treatment as an example to explain how we trained and applied the DQN in our models.

The given dataset makes up 80% of the entire dataset, which contained 152,527 patients and 5,166,669 outpatient visits. To define outpatient visit samples from the EMRs, we used lab test results (or physical measurements) within a certain time interval before the visit as state variables of its corresponding sample. The complications that occurred before the current outpatient visit were treated as medical history. The criteria used in the generation of a given dataset are shown in Figure 2. Subsequently, the patients in the given dataset were randomly divided into a training set and validation set at a ratio of 8:2.

The features that were selected to define the states for a patient visit included demographics, medical history, disease risks, previous drugs, lab data, and physical measurements. The detailed state information is presented in Table S1 in Multimedia Appendix 1. Among these features, continuous variables were normalized into a common scale, whereas binary variables were represented using 0 or 1. Other categorical variables were converted into multiple binary variables using one-hot encoding. Finally, we obtained a 49-dimension state vector.

For patients with T2DM, antiglycemic medications are usually prescribed based on the currently used drugs [29,32,33]. Thus, to simplify the action space of the DQN, we defined an action of a visit based on medication changes from the previous prescription. The prescription changes included drug changes at the drug class level and dosage changes of some common drugs. In terms of dosage changes, the increase or decrease in dosage of the three most frequently used drugs, namely metformin, basal insulin, and premixed insulin, was considered. Table 1 lists the actions for antiglycemic treatment used in the DQN model. The most common medication adjustment options in the actions were: changing the dosage of a drug, adding an oral antidiabetic drug (OAD), and changing to insulin treatment. Among the action options, “No prescription change” indicates use of the same drugs and dosages as the previous prescription, and “Using xxx insulin” means changing to the specific insulin or insulin combinations.

The reward function is usually defined to quantify the effectiveness of the action at each time step. In the antiglycemic DQN, for a patient visit at time t, we defined the reward function r_t as shown in Equation (3):

reward_gly=a · sgn(7 – s_t+1^HbA_1c) + (–b) · s_t+1^Hypo + (–c) · s_t+1^Final · s_t+1^CX(3)

where

(4)

In Equation 3, s_t+1^HbA_1c (%) is the HbA_1c of time t+1, s_t+1^Hypo, is a binary (1,0) variable representing whether hypoglycemia occurs before time t+1, s_t+1^Final is a binary (1,0) variable representing whether time t+1 is the final visit of the patient, and s_t+1^CX ∈ {–1,1} indicates whether complications or death occur at time t+1, where –1 represents “No” and 1 represents “Yes.”

The concept underlying the reward function is to give a positive reward when (1) the HbA_1c after 3-6 months reaches the control target (less than 7%), and (2) no complications or death occurred until the last visit of a patient in the next 6 years. A negative reward (ie, penalty) is given when (1) the HbA_1c after 3-6 months is not well controlled, (2) a hypoglycemia event occurs in the next 6 months, and (3) a complication or death occurs after the current visit. Based on the importance of these outcomes, we set the coefficients in Equation (1) as follows: a=1, b=2, c=4. For an intermediate visit of a patient, the DQN model is trained to optimize the cumulative reward, which is equal to the current reward plus the next visit’s expected cumulative reward multiplied by a discount factor, γ=0.9. Therefore, the DQN model is able to estimate the impact of a current action on both short-term and long-term outcomes.

The network architecture and training settings are provided below. We adopted a fully connected neural network with 2 hidden layers of 64, with 32 units for the Q networks. Each hidden layer contained batch normalization and Leaky-ReLU activation. The input layer was 49 dimensions and the output layer was 14 dimensions, which were the same as the sizes of the state vector and the action space. The learning rate η was 0.001, the batch size was 256, and the target network update parameter τ was set to 0.01. For regulation, we set the reward threshold r_reg=4 and λ=0.5. We trained the DQN model for a maximum of 100,000 iterations using the Adam optimizer [42].

For antihypertensive and lipid-lowering treatments, actions and reward functions of DQNs are shown in Table 2 and Table 3, respectively.

In the antihypertensive DQN, for a patient visit at time t, we defined the reward function as in Equation (5):

s_t+1^SBP(mmHg) is the SBP of time t+1, s_t+1^DBP (mmHg) is the DBP of time t+1, and other terms are defined as in the reward function of the antiglycemic DQN.

In the lipid-lowering DQN, for a patient visit at time t, we defined the reward function as in Equation (7):

Reward_lip = a · sgn(2.6 – s_t+1^LDL-c) + (–c) · s_t+1^Final · s_t+1^CX (7)

where s_t+1^LDL-c (mmol/L) is the LDL-c of time t+1, and other terms are defined as in the reward function of the antiglycemic DQN.

Of the 36,993 patients in the test data, 18,878 were men (51%). With respect to ethnicity, the majority of the patients were Chinese (69%), followed by Malay (15%) and Indian (11%). By 2019, the median age was 67 years (IQR 59-76) and the median duration of diabetes was 10 years (IQR 6-16).

An overview of the short-term and long-term outcome evaluation cohorts is shown in Figure 3. Further information is provided in the following two subsections.

To evaluate the short-term outcomes, we took model concordance as the exposure variable to evaluate the effect of our treatment recommendation model at the patient-visit level. We generated the test datasets separately for different short-term outcomes, namely the percentages of patient visits with well-controlled parameters (HbA_1c<7% [53 mmol/mol], SBP/DBP<140/90 mmHg, LDL-c<2.6 mmol/L [100 mg/dl]) after 3-6 months of therapy. For a hypoglycemia event, the occurrence rate in the following 6 months was compared between two groups. In short-term evaluation, potential confounding factors were adjusted by stratification, the PS weighting method, and multivariate regression, such as age, gender, and ethnicity.

Specifically, the model concordance was defined at the level of the standard drug class, since the treatment recommendation models output the standard drug class rather than the specific brand name. For example, if the model recommends alpha-glucosidase inhibitors (AGI), the patient visit with a prescription for Acarbose (a type of drug belonging to the drug class AGI) is model-concordant. Only the top-ranking recommended medication for antiglycemic, antihypertensive, and lipid-lowering therapy was considered for evaluation.

After meeting all exclusion criteria, a total of 178,489 visits were included to evaluate the short-term clinical outcomes of HbA_1c control. Of the total samples, 78,670 patient visits (44.08%) were model-concordant and 99,819 (55.92%) were nonconcordant. The characteristics of the model-nonconcordant and model-concordant groups in the test data are shown in Table 4 for short-term blood glucose control evaluation, and patient characteristics for other short-term outcomes (eg, blood pressure control) are shown in Tables S2-S4 in Multimedia Appendix 1.

The evaluation results of short-term outcomes, including HbA_1c, SBP/DBP, LDL-c control, and hypoglycemia event, were based on test samples with corresponding outcome data during the follow-up period. After confounder adjustment for patient characteristics, the model-concordant treatments were associated with good blood glucose control, good blood pressure control, and good blood lipid control compared with model-nonconcordant treatments. There was no significant difference in the occurrences of hypoglycemia events between model-concordant treatments and model-nonconcordant treatments (Table 5).

We further stratified the patient visits by the confounder that was most strongly correlated to the clinical outcome. For the outcome of glucose control (ie, HbA_1c after 3-6 months), the current HbA_1c of the patient was used to stratify the patient visits into three groups of low (<7%), medium (7-9%), and high (>9%) levels. The short-term evaluation was performed on each group separately (Table 6), showing that model-concordant treatments were associated with improved short-term HbA_1c outcomes (ie, higher HbA_1c goal-achieving rate) in each group.

Figure 4 illustrates the relationship between the patient’s model-concordant rate and the occurrence rate of long-term clinical outcomes for all patients with respect to antiglycemic, antihypertensive, and lipid-lowering therapy. Specifically, the patients were divided into different groups according to the patient’s model-concordant rate (eg, every 20% as a group), and the occurrence rate of complications or death in each group was computed. In general, the curves show a downward trend. In other words, there is a negative correlation between the model-concordant rate and the occurrence rate of complications or death; the higher the patient’s model-concordant rate, the lower the occurrence rate of complications or death. Table 7 shows the slope of each curve by fitting to a straight line using all data points, which indicates the extent of the downward trend. In addition, the number of patients in each long-term outcome evaluation cohort is shown in Table S5 in Multimedia Appendix 1.

Furthermore, for combined treatments with antiglycemic, antihypertensive, and lipid-lowering drugs, we evaluated if the patient’s model-concordant rate for the three types of treatments had a positive impact on the reduction of the complication or death risk by multivariate regression. Only patients with all three model-concordant rates of antiglycemic, antihypertensive, and lipid-lowering treatment were included in the multivariate regression model. Table 8 shows the multivariate regression results of each long-term outcome with confounder adjustment for the corresponding risk score generated by the prediction model. All of the prediction models, based on XGBoost, outperformed the clinical baseline models [49-53] and demonstrated good prediction capability, with an area under the receiver operating characteristic curve ranging from .711 to .874. The model-concordant rate for antiglycemic treatment had a negative correlation with the occurrence of major complications and death, with coefficients ranging from –1.12 to –.33. A similar result was found for the model-concordant rate for antihypertensive treatment (coefficient range –1.44 to –.40) and the model-concordant rate for lipid-lowering treatment (coefficient range –1.17 to –.52). All of these coefficients were significant (P<.05), except for the coefficients of antiglycemic treatment and antihypertensive treatment in the evaluation of stroke outcome. This implies that the patients whose treatments were more concordant with the model recommendation were more likely to be associated with a lower risk of diabetes complications (including both macrovascular and microvascular complications) and death. All of the coefficients for the risk score had positive values, which also validated the soundness of the risk prediction models. In addition, the number of patients with or without the corresponding long-term outcomes in multivariate regression are shown in Table S6 in Multimedia Appendix 1.

In all patient visits, the percentages of model-concordant visits were 43.30%, 51.25%, and 58.93% for antiglycemic, antihypertensive, and lipid-lowering therapy, respectively. To analyze the distribution characteristics for model concordance, we compared the antiglycemic medication patterns between physicians’ prescriptions and model recommendations.

Figure 5 depicts the medication patterns used by physicians (top panels) and the model (bottom panels) for the three different groups. All patient-visit samples were categorized into three groups based on the current HbA_1c as follows: low (<7%), medium (7-9%), and high (>9%). For each group, the medication patterns of physicians’ prescriptions and model recommendations are visualized by 2D histograms, in which the x-axis represents the number of OADs and the y-axis represents the intensity of insulin (a value of 0 indicates no insulin used, 1 indicates single use of basal insulin or prandial insulin, and 2 indicates a premix or combination use of basal and prandial insulin, or others). The color indicates the usage number of corresponding medication patterns. First, Figure 5 shows that medication patterns of model recommendations are consistent with clinical knowledge, as most patients in the low group were prescribed with a single OAD, whereas patients in the medium and high groups showed increased use of multi-OAD, insulin, and insulin plus OAD. Second, the medication patterns of model recommendations are visually similar to those of physicians’ prescriptions in the low group, whereas in the medium and high groups, the patterns of model recommendations are more vigorous than those of physicians’ prescriptions. This indicates that the model learns from the data, showing that active adjustment of the medication for the patients in medium and high groups may be associated with a better clinical outcome.