Deep doubly robust outcome weighted learning

Abstract

Precision medicine is a framework that adapts treatment strategies to a patient’s individual characteristics and provides helpful clinical decision support. Existing research has been extended to various situations but high-dimensional data have not yet been fully incorporated into the paradigm. We propose a new precision medicine approach called deep doubly robust outcome weighted learning (DDROWL) that can handle big and complex data. This is a machine learning tool that directly estimates the optimal decision rule and achieves the best of three worlds: deep learning, double robustness, and residual weighted learning. Two architectures have been implemented in the proposed method, a fully-connected feedforward neural network and the Deep Kernel Learning model, a Gaussian process with deep learning-filtered inputs. We compare and discuss the performance and limitation of different methods through a range of simulations. Using longitudinal and brain imaging data from patients with Alzheimer’s disease, we demonstrate the application of the proposed method in real-world clinical practice. With the implementation of deep learning, the proposed method can expand the influence of precision medicine to high-dimensional abundant data with greater flexibility and computational power.

Appendix

1.1 Appendix A: More on data extraction and preprocessing of the clinical data

Different centers have different data collection methods and policies during different time periods which might not conform to each other. Given that such conformity issues could create noisy heterogeneity in a problematic way, some stringent inclusion and exclusion criteria were applied to keep the multi-center, multi-stage data relatively clean. We found that there were more missing data in earlier form versions and later visits, and the most information was collected at the initial visits. Since we were interested in baseline information, only observations of the first visit of each subject with the latest form version were included. For categorical variables, categories such as unknown (e.g., 9, 99, 999, 8888, 9999 values), not applicable (form submitted did not collect such data or a skip pattern precludes such responses), or left blank were considered missing. For continuous variables, indicators of unknown, not assessed, and not available (e.g., − 4, 888.8) and extreme values outside of the normal range (e.g., height more than 80 inches) were considered missing.

Because the UDS dataset has many forms/variables that contribute to the assessment of the subject’s cognitive status, covariates that have moderate to high estimated Pearson correlations ($> 0.5$) with the outcome variable were excluded to avoid multicollinearity. Covariates with high estimated Pearson correlation ($> 0.8$) with other non-outcome covariates were excluded as well (e.g., height/weight and body mass index, and various CDR^® scores). The ID data, data containing information that links with UDS and MRI, were processed similarly as the UDS data, such as removing severe missing data and multicollinearity.

There were 48 clinical variables included in the analysis (Table 8). They are: demographics (age, gender, smoking status, research center), vitals (blood pressure, resting heart rate), previous and current medication/therapy (antiadrenergic, anxiolytic, anticoagulant, antidepressant, angiotensin converting enzyme inhibitor, antipsychotic, diuretic, estrogen hormone therapy, lipid lowering medication, nonsteroidal anti-inflammatory, vasodilator), previous medical history (anxiety, diabetes, angioplasty, cardiac bypass procedure, vascular brain injury, heart valve replacement, microhemorrhage, incontinence, insomnia, obsessive-compulsive disorder, Parkinson’s disease), number of visits, and time to complete trail making test.

Variables related to smoking (e.g., consumption of tobacco in the past 30 days, number of years smoking, whether or not the subject quit smoking, etc.) were combined into two variables: an indicator of current smoker (has not quit smoking, or has consumed tobacco in the last 30 days) and an indicator of ever smoked (smoked cigarettes in the last 30 days, more than 100 cigarettes in life, non-zero years of smoking, or at least 1 cigarette smoked per day on average). Diabetes was redefined as a binary variables with 1 representing Type 1, Type 2, or other type of diabetes such as diabetes insipidus, latent autoimmune diabetes or Type 1.5, and gestational diabetes, and 0 representing no diabetes reported.

There were a total of 5,616 MRI sessions available. Each session contains multiple DICOM files, with one file representing one MRI slice. The DICOM format was preferred because it contains image information such as slice position and sequence type in the headers. We extracted the MRI slices from each subject’s MRI session in a compressed folder and removed slices without series description or image position because series description informs the sequence type of the MRI scan and image position helps sort the slices in the right order. Since the consecutive slices differed by a matter of milliseconds, we selected every 5th slice from the 150 middle slices to save space and maintain the same reasonable image dimension for every subject. End slices were discarded as they contain less useful information about the brain. We restricted to one MRI per subject because we were only interested in baseline covariates and wanted to keep the input dimension consistent.

The preprocessed ID and preprocessed UDS were merged by subject ID, their center, and visit year. Only complete cases were used because addressing data incompleteness was not our main research focus here, and imputation on such multi-center observations was often unreliable or needed extremely careful manipulation. In the preprocessed data (before merging with MRI scans but after merging UDS with ID), there were 424 observations collected from 12 centers spanning from 2015 to 2019. Among the 424 subjects, 48% had better or maintained normal cognitive status and 48% currently used antihypertensive or blood pressure medication at the initial visit. The preprocessed MRI data were merged with the medical data by file locator information so only subjects who had qualified UDS, ID, and MRI data were included, resulting in a sample size of 186. We further excluded 24 subjects whose MRI date was more than 200 days away from the initial visit date. The final sample size for the preprocessed data was 162. As mentioned above, we applied a relatively strict inclusion criteria to make sure the input data for the DL models were relatively clean and similar. The MRI data were only lightly processed to preserve the original values but could be piped through more systematic image processing tools.

1.2 Appendix B: More on transfer learning and ResNet34

Transfer learning can be regarded as a feature selection tool because images often contain a large amount of nuisance pixels. The outputs of ResNet34 prior to the dense layers have a lowered dimension of 1000, much smaller and more extracted than the original dimension. If we fed them together with medical data directly into a DL architecture, the MRI data would dominate the dimension. An alternative to transfer learning is applying unsupervised learning such as autoencoder (AE) to the MRI data. AE is a good dimension reduction method but the encoded outputs are sometimes not necessarily good representations of the original input.

In addition to the reasons mentioned in the main text, ResNet34 was chosen as the pretrained model because it has lower model complexity and relatively low top-1 and top-5 errors on the ImageNet data compared wth other famous deep learning architectures such as AlexNet or VGG. Top-1 error means the proportion of test images whose true label does not match with the prediction class with the highest estimated probability. Top-5 error means the proportion of test images whose true label is not among the 5 prediction classes with the top 5 highest estimated probabilities. We applied a pretrained model instead of training our own structure because the lower-level representations extracted from the earlier layers of existing models are generally transferrable across images.

1.3 Appendix C: Definition of Value Functions

We used the same CV definition of value function estimator and its variance estimator as in Jiang et al. (2021). Let $j=1,\ldots , MK$ denote all MK tuning folds regardless of repetition across M repetitions and K CV folds and $i = 1, \ldots , n_j$ be the ith observation in the jth overall fold. We applied cross validation on a dataset of size $n_{tr}$ which will be split into training and validation sets. The CV estimated value function was used to compare tuning performance:

$$\begin{aligned} \widehat{V}^{cv} (\hat{d}_{n_{tr}}^{(-j)}) = \frac{\sum _{j=1}^{MK}\sum _{i=1}^{n_j} U_{ji} }{\sum _{j=1}^{MK}\sum _{i=1}^{n_j} W_{ji}}, \end{aligned}$$

where $W_{ji} = \frac{1\{A_{ji} = \hat{d}_{n_{tr}}^{(-j)} ({\varvec{X}}_{ji}) \}}{\hat{P}^{(-j)}(A_{ji} \vert {\varvec{X}}_{ji})}$, $U_{ji} = Y_{ji}W_{ji}$, $\hat{d}_{n_{tr}}^{(-j)}$ is the decision rule estimated from the dataset of size $n_{tr}$ with the jth fold left out, and $\hat{P}^{(-j)}(A_{ji} \vert {\varvec{X}}_{ji})$ is the estimated propensity score. Its variance is used to measure the estimation uncertainty

$$\begin{aligned} \widehat{\text {Var}}[ \widehat{V}^{cv} (\hat{d}_{n_{tr}}^{(-j)})] = \frac{1}{K(MK-1)} \sum _{j=1}^{MK} \sum _{i=1}^{n_j} R_{ji}^2, \end{aligned}$$

where $R_{ji} = \frac{1}{\bar{W}_j} U_{ji} - \frac{\bar{U}_j}{\bar{W}_j^2} W_{ji}$ is an influence function-inspired form of the value function with $\bar{U}_j = \sum _{i=1}^{n_j}U_{ji}$ and $\bar{W}_j = \sum _{i=1}^{n_j} W_{ji}$. By definition, $\sum _{j=1}^{MK}\sum _{i=1}^{n_j} R_{ji} = 0$.

For testing results, there is no CV so $j = 1, \ldots , M$. The estimated value function becomes

$$\begin{aligned} \widehat{V}(\hat{d}_{n_{te}}) = \frac{\sum _{i=1}^{n_{te}} U_i}{\sum _{i=1}^{n_{te}} W_i}, \end{aligned}$$

where $n_{te}$ is sample size of the testing set, and $U_i, W_i$ are defined similarly as $U_{ji}$ and $W_{ji}$ but with $i=1,\ldots , n_{te}$ and decision rule $\hat{d}_{n_{te}}$. Its variance is given as

$$\begin{aligned} \widehat{\text {Var}} [\widehat{V}(\hat{d}_{n_{te}}) ] = \frac{\sum _{j=1}^{M} (\widehat{V}(\hat{d}_{n_{te}, j}) - \bar{V}(\hat{d}_{n_{te}, M}))^2}{M-1}, \end{aligned}$$

where $\hat{d}_{n_{te}, j}$ is the single estimated decision rule from the jth repetition and $\bar{V}(\hat{d}_{n_{te}, M}) = \sum _{j=1}^{M} \widehat{V}(\hat{d}_{n_{te}, j})$ is the average estimated value functions over M single estimated decision rules. The SD is the square root of the variance.

The performance of the model during tuning and testing was determined by higher estimated value functions and lower SDs.