All data were collected from the EHRs at the University of Tokyo Hospital, which included 7034 patients who visited the hospital, had diagnostic codes, and were registered to the Japan diabetes comprehensive database project based on an advanced electronic medical record system (J-DREAMS) cohort [42]. The data were recorded in the EHRs between January 1, 2011 and December 31, 2022. The data, including treatment decisions and outcomes, were reflective of care by endocrinologists. Variables extracted from the EHRs included sex, age, 12 kinds of diabetes-related laboratory tests, and drugs. The laboratory tests included numerical values of HbA_1c, glucose, triacylglycerol, high density lipoprotein cholesterol, total cholesterol [43], urinary albumin creatinine ratio, creatinine [44], alanine transaminase, aspartate transaminase, and γ-glutamyltransferase [45], and categorical values of proteins [45] and glycogen. Drugs were identified using the list of drug price standards [46] provided by the Ministry of Health, Labour and Welfare of Japan, and 44 types prescribed in 2021 were selected.

The records used for training were not used for testing to ensure that the same patients were not included in both groups. A total of 80% (5627/7034) of patients were included in the training group, and the remaining 20% (1407/7034) of patients were included in the testing group.

In order to examine the size of the training data needed to achieve the target prediction accuracy, we extracted three different subsets of training data from the sequences of the 5627 patients in the training group: 2 years from 2020 to 2021 (3013 patients and their 25484 subsequences), 5 years from 2017 to 2021 (4009 patients and their 78,020 subsequences), and 10 years from 2012 to 2021 (4524 patients and their 168,595 subsequences). For testing, we further identified a subset of the testing data that was solely data from 2022 (637 patients and their 2988 subsequences). Thus, there was no overlap in patients or time periods between training and testing. Table 1 shows the characteristics of the patients. The drugs in the table are the top 5 and bottom 5 in terms of number of prescriptions in the 2-year training data. All characteristics are provided in Multimedia Appendix 1.

Patients’ medical conditions and prescription histories are in general irregularly spaced, reflecting variability in patient care appointment dates, with updates to outpatient EHR occurring before and after clinical visits. We organized the data into Monday-to-Sunday weeks and quantized the data to a single value per week, using the average in the case of multiple measurements and treating weeks with no values as having missing values [47]. This approach allowed the ML model to treat irregularly spaced data spanning Y (a natural number) years as regularly spaced data consisting of (Y×365)⁄7 (rounded up to the nearest integer) values, that is we treated all data as weekly data. We did not perform preprocessing, including interpolation, on missing values in the regularly spaced data. No normalization, outlier removal, or dimensionality reduction was performed.

We designed a transformer-based encoder-decoder model (Figure 2) that takes as input a time series of drugs prescribed and laboratory tests over the past 1 year, sex, and age. The model approaches drug selection as a multichoice task. With N types of drugs, there are 2^N potential prescription combinations. By setting the number of units in the output layer of the transformer decoder to N, we implemented N binary classifications. The model outputs a set of scores representing the probability of each drug being prescribed on that day and, for each drug and day, a binary prescription decision based on whether the prescription probability is greater than or equal to 0.5.

We treated age and HbA_1c values as numerical data, sex as a categorical label, and prescription history as a set of categorical labels. In addition, a time series representing the presence or absence of missing data, called an attention mask, was simultaneously generated. Each data type was transformed into a uniform-dimensional vector using an embedding layer specific to that data type, which was then handled by the transformer model. Time information was converted to a uniform-dimensional vector using a positional embedding layer with a periodic function and added to the output value of the embedding layers.

The model incorporates two types of attention layers: self-attention in the encoder, designed to extract relationships in time and meaning from the time series of drugs prescribed and laboratory tests, and cross-attention in the decoder, used to predict the next prescribed drugs based on these relationships. Times in the time series that contain missing data are ignored in the self-attention and mutual attention calculation process using the attention mask. This makes it easy to handle without additional processing for missing data completion. The self-attention weights were optimized through self-supervised learning. This involves the task of predicting both the next laboratory test values and drugs prescribed using a time series of past laboratory tests and drugs prescribed. The cross-attention weights were optimized through supervised learning involving the task of predicting a set of scores representing the probability of each drug being prescribed on that day.

The model consisted of four transformer-encoder layers including four multihead self-attention blocks and four transformer-decoder layers including four multihead cross-attention blocks, along with a hidden layer of dimension 256. Compared to the text data handled by the language model, the prescription drug selection data has a smaller vocabulary (number of drug types) and a shorter time series, so these hyperparameters were set to about one-third the size of the BERT model [34]. The parameters were optimized by Adam with a learning rate of 1e-4, a batch size of 256, and 100 epochs. The loss functions for numerical and categorical data were mean squared error and focal loss [48], respectively. All implementations were written in Python 3.11 (Python Software Foundation) and PyTorch 2.2 (Meta AI).

The model was trained using prescription records from endocrinologists at the University of Tokyo Hospital. As a result, the model generates outputs aligned with the treatment approaches of these specialists.

We analyzed the characteristics of patients in the dataset using means, SD, and frequency counts. We calculated 95% CI using the bootstrap method where applicable. We performed all statistical analyses using custom Python code.

We compared our model with an established ML method recognized for high accuracy. There were validations on similar T2D prediction tasks favored LightGBM [49,50], making it our chosen reference for comparisons. We compared the predictive accuracy of the two methods with the three different training dataset for the prescription prediction of 44 drugs using both macro- and microaverages [51] of the receiver operating characteristic area under the curve (ROC-AUC) as metrics. The macroaverage is a measure of the average performance across all classes independently. The method computes the metrics of ROC-AUC for each individual drug, then averages them to obtain an overall score. This analysis method gives equal weight to each drug, regardless of its size or imbalance in the dataset. On the other hand, the microaverage is a measure of the aggregate performance that weights all instances in the dataset equally. This analysis method first aggregates the true positives, false positives, true negatives, and false negatives across all drugs, and then computes each metric using these aggregated values. The microaverage treats every prediction equally, without considering the kinds of drugs.

The predictive performance of each drug was also evaluated using ROC-AUC.

This study was approved by the institutional review board of the University of Tokyo School of Medicine (approval number: 10705-(4)) and was conducted in accordance with the Declaration of Helsinki. This was a retrospective, noninterventional database study without patient involvement. Confidentiality was safeguarded by the University of Tokyo Hospital. According to the Guidelines for Epidemiological Studies of the Ministry of Health, Labour and Welfare of Japan, written informed consent was not required. Information about the current study was available to patients on a website, and patients have the right to cease registration of their data at any time [52].

We assessed different sizes of training data (Table 2). The best prediction performance was obtained from training with the 5 years of data from 2017 to 2021. This version achieved a microaverage (95% CI) ROC-AUC of 0.993 (0.992-0.994) and a macroaverage (95% CI) ROC-AUC of 0.988 (0.980-0.993), and we selected it as our model. Performance was similar when trained using the full range of 10 years of data, from 2012 to 2021, producing somewhat worse results. For all sizes of training data, we met the study’s objectives of achieving an ROC-AUC above 0.95 from our macro- and microaverage evaluations. The prediction accuracy of LightGBM with the 5 years of data had a microaverage ROC AUC of 0.988 (0.985-0.990), and the transformer outperformed LightGBM.

We examined the prediction performance for each of the 44 drugs (Table 3). The drugs in the table are the top 5 and bottom 5 in terms of the number of prescriptions in the 2 years of training data. All results are provided in the supplemental file (Multimedia Appendix 2). We achieved an ROC-AUC above our target of 0.95 for 43 of the 44 drugs when trained with the 5 years of data.

The only drug that did not achieve the target value was “Insulin lispro (genetical recombination) [Insulin lispro Biosimilar 1].” The prescription of this drug began in 2020. Therefore, the same instances of prescription were present in all three training data periods. For this drug, the model trained on just 2 years of data had the highest ROC-AUC.

The proposed model performed as well as other transformer-based models considering the ability for extracting relationships in time and meaning [53]. The embedding vectors obtained through training represent the closeness of relationships between vectors as proximity. The embedding vectors in this experiment were the same 256 dimensions as the transformer hidden size. Projecting onto two dimensions using uniform manifold approximation and projection [54] allows visualization (Figure 3). While we did not observe a strong tendency for clustering, several biosimilar drugs, such as insulin glargine (genetical recombination, ie, insulin glargin biosimilar 1), were positioned close to the biosimilar drug.

The proposed model achieved impressive predictive accuracy, with macroaverages (95% CI) of ROC-AUC of 0.988 (0.980-0.993). In previous studies [11,12], the ROC-AUC for drugs was less than 0.95, and our model represents a significant improvement relative to this previous work.

The prediction accuracy was higher when trained with short-term data from 5 years than when trained using data from the past 10 years. We suspect that changes in the treatment environment, such as the introduction of new prescription drugs, are responsible for this difference in prediction accuracy. When training on data from the past N years, the model uniformly learns the prescription selection trends for the past N years. Therefore, when N is large, the model is heavily influenced by the trends of older prescription selections. This suspicion is supported by our results that the prescription drug that did not reach target accuracy, “Insulin lispro (genetical recombination) [Insulin lispro Biosimilar 1],” was only recently approved and prescribed. Introduced in 2020, this drug has remained relatively rare in the dataset, as demonstrated in Table 1. Even when expanding the training data period from 2 to 10 years, the number of prescriptions for this drug has remained stagnant at 59, and its proportion within the overall dataset has decreased from 0.23% to a mere 0.03%. This low frequency and temporal bias in the data have rendered learning this specific drug’s patterns challenging. Although ML models generally perform better with more data, for this application it may be better to gather data from more hospitals rather than a longer time span, since it is desirable to learn more fresh data, including new prescription drugs. Transformer models, known for power-law characteristics, benefit from scale-ups [38], and expanding the study to multiple hospitals could explore potential performance enhancements and test the applicability of the power-law in the medical field.

Physicians select drugs considering insulin secretion and insulin resistance, age, obesity, severity of chronic complications, and liver and kidney function. In this experiment, we used age, sex, 12 kinds of diabetes-related laboratory tests, and past prescription history. By adding comprehensive items likely used in medication selection to the model input, performance may be improved.

Our ultimate goal is to improve the treatment outcomes of diabetes. Merely predicting drug selection alone cannot achieve this goal. Expanding the scope to predict the impact of prescribed drugs could further enhance the model’s utility in diabetes treatment.

Our study has notable limitations. First, the model uses only age, sex, 12 kinds of diabetes-related laboratory tests, and past prescription history as inputs. It is desirable to fully account for the various constraints and contraindications that physicians consider in real-world clinical practice. For example, when considering patient characteristics, factors such as age (older adults or pediatric), pregnancy, lactation, BMI, type of diabetes, c-peptide, renal and hepatic function, comorbidities, allergies, urinary tract infection history, diabetic ketoacidosis history, cardiac history, and hypoglycemia must be considered. Regarding medications, considerations include adverse effects, drug resistance, and cost. However, despite these limitations, this study demonstrates that it is possible to narrow down treatment options to some extent using a limited set of variables. We believe that this result holds promise for providing primary care physicians with some guidance and direction. Regarding warnings of constraints and contraindications in drug selection, CDSS have traditionally excelled in this [20,21] and can effectively complement this model. In addition, recent advancements in natural language processing within ML models [33,34] have enabled the extraction of important information on constraints and contraindications from fresh sources such as research papers published daily. We believe there is substantial potential for future research to integrate these important factors into ML models, and we intend to work on this improvement.

Second, the data were sourced from a single hospital, limiting generalizability. Potential biases in the results include the race and geography of the patients and the prescribing tendencies of a limited number of endocrinologists at the hospital. ML models tend to learn biases present in the training data, resulting in predictions that reflect those biases. Prediction accuracy is significantly degraded when different biases exist between the training data and the test data. Therefore, in order to generalize the prediction results of the model, a large dataset of data from multiple hospitals containing data on patients with diverse backgrounds is required. However, in reality, it is difficult to build an ideal dataset due to constraints such as privacy protection and data collection costs. As one solution, we believe that it is good to collect training data that resembles the patterns of patients treated by the primary care physicians who are the users of the model, thereby eliminating biase differences between the training data and the test data as much as possible. Specifically, we collect EHRs of patients who have been involved in prescriptions by endocrinologists at multiple hospitals in a specific region, and use this to build the model that selects prescription drugs for patients treated by primary care physicians in the same region. We believe that this approach will reduce the bias differences in regional characteristics between the training data and the test data, and improve generalization performance with realistic data collection. We would like to verify the validity of this solution in the future.

Third, ML reflects majority characteristics, potentially limiting applicability to diverse patient populations. In the dataset used in the experiment, over 40% (297/637) of patients were prescribed metformin hydrochloride, and patient characteristics are biased. There is a risk that truly effective treatments may not be prescribed correctly to a minority of patients. Prediction failure analysis needs to be further scrutinized, including versus patient characteristics. We should examine this issue by comparing prediction accuracy for each patient cluster.

Fourth, this was a backward-looking study, using past data, and the essential next phase is to assess the model’s predictive capabilities in clinical practice. There is a need for a careful exploration of the model’s effectiveness in real clinical scenarios.

Fifth, while the proposed model requires hyperparameters, we determined these values based on previous studies. Although fine-tuning hyperparameters is generally desirable, transformer models are typically computationally expensive to train, and many previous studies have likewise not fully tuned their parameters. There is potential for performance improvement through hyperparameter tuning, and we intend to investigate this further in future work.

Sixth, the model has poor interpretability. We investigated the proximity relationship between embedding vectors, but no strong tendency was found. It would be better to use other information obtained from the transformer, such as attention weights, to further improve interpretability.

These limitations raise important ethical considerations in the development and application of medical artificial intelligence. Physicians who use the model must be aware of these limitations and exercise appropriate clinical judgment. Obtaining appropriate informed consent from patients is important.

The proposed model addresses the challenge of predicting the next prescribed drugs. This model, trained using past prescriptions of endocrinologists, has the potential to improve treatment outcomes for nonspecialists by assisting them in making prescription decisions. Future efforts should focus on improving accuracy by incorporating disease state information beyond current inputs and validating the model on large clinical datasets across multiple hospitals.