Are you susceptible to adverse drug reactions?

Prediction and insights using LightGBM and SHAP — Part two

Chiao-Feng Lin
DNAnexus Science Frontiers
6 min readFeb 16, 2021

--

In part one of this two-part blog, we described how we used UK Biobank (UKB) clinical and genomic data, and the LightGBM framework to build a model for predicting an individual’s susceptibility to Adverse Drug Reactions (ADRs). In this second part, we applied SHAP TreeExplainer to our LightGBM results and gained insights about ADR. We found that both the psychology and physiology of an individual play a role.

SHAP Explanations

We used the SHAP package to understand the trained model and predictions made from it. Installing SHAP is straightforward and it works nicely with LightGBM output.

Summary plots

In SHAP, a model is explained by how much the predicting features impact the prediction from the trained model. The impacts are quantified by SHAP values, which can be positive or negative. Figure 1 shows the top 20 features and the combined magnitude of their impact on predicting all classes. The feature that had the highest impact (the highest mean SHAP value) was “Seen a psychiatrist for nerves, anxiety, tension or depression.” Its margin over the next one was much higher than all the others.

This was surprising as we did not expect that mental state would have anything to do with ADRs. To confirm, we first checked the raw counts in each category (Figure 2). Chi-square tests gave p-values of 0 (class 1 vs. control (class 0)), 0 (class 2 vs. control), and 0 (class 3 vs. control), indicating a strong association between this psychological feature and ADRs. Even though this association is not necessarily causation, it appeared to be a strong predictor for class 1 and class 3. We followed up with a literature search and found that there were related reports. For example, we found a study investigating the connection between mental state and drug allergy that had inconclusive results due to its extremely small sample size (61 cases and 55 controls).

Figure 1: Top 20 impactful features on the ADR prediction model from SHAP analysis
Figure 1: Top 20 impactful features on the ADR prediction model from SHAP analysis. Class 1: individuals having ADR-related ICD10 codes in the main diagnosis (N=1,478). Class 2: individuals having ADR-related ICD10 codes in the secondary diagnosis (N=26,024). Class 3: individuals having ADR-related ICD10 codes in both fields (N=1,828). Class 0 (control): individuals not having ADR-related ICD10 codes in both fields.
Figure 2: Histograms of feature “Seen a psychiatrist for nerves, anxiety, tension or depression” for the four classes.
Figure 2: Histograms of feature “Seen a psychiatrist for nerves, anxiety, tension or depression” for the four classes. Proportions of individuals answering Yes (coding 1) and No (coding 0) are very different between class 0 (control) and the other three classes. (Coding -1: Do not know. Coding -3: Prefer not to answer)

Currently, the SHAP package only supports bar plots for multiclass classification as shown in the above global summary plot. Dot plots can reveal the direction of the impact, but they are only available for binary classification. Therefore, the following analyses were all based on binarized classification. We also looked at how individual features impact the prediction of each class. Consistent with the multiclass global plot, class 1 and class 3 show similar plots. Class 2 and class 3 have overlapping features (Figure 3a-c). Some trends are intuitive and easy to explain, e.g. “Taking other prescription medication” and “Number of treatments/medications taken.” Others are less so. For instance, “Length of time at current address” is interesting. As we know, moving homes is extremely stressful. Presumably, the shorter the stay at one address, the stronger the effect of the lingering stress from the move.

Figure 3a: Top 20 impactful features on the model for predicting class 3
Figure 3a: Top 20 impactful features on the model for predicting class 3
Figure 3b: Top 20 impactful features on the model for predicting class 1
Figure 3b: Top 20 impactful features on the model for predicting class 1
Figure 3c: Top 20 impactful features on the model for predicting class 2
Figure 3c: Top 20 impactful features on the model for predicting class 2

Dependence plots

Compared to summary plots, SHAP’s dependence plots give better resolution on how an individual feature impacts the prediction (represented as SHAP values). More importantly, they help examine whether or how one feature interacts with another. For example, Figure 4a clearly shows that answering yes (1 at x-axis) to “Seen doctor (GP) for nerves, anxiety, tension or depression” is more likely (higher SHAP value) to be classified as class 3. Looking at the summary plots, we might suspect that “Seen a psychiatrist for nerves, anxiety, tension or depression” and “Seen doctor (GP) for nerves, anxiety, tension or depression” are related. The dependence plot shows that among the people who answered yes to “Seen doctor (GP) for nerves, anxiety, tension or depression”, those answering yes (red dots) to “Seen a psychiatrist” are less likely to be classified as class 3.

Figure 4a: Dependence plots for the impact of psychology-related features on predicting class 3
Figure 4a: Dependence plots for the impact of psychology-related features on predicting class 3. Answering yes (1 at x-axis) to “Seen doctor (GP) for nerves, anxiety, tension or depression” is more likely (higher SHAP value) to be classified as class 3 than answering no (0 at x-axis). (Coding -1: Do not know. Coding -3: Prefer not to answer)
Figure 4b: Dependence plots for the impact of psychology-related features on predicting class 3
Figure 4b: Dependence plots for the impact of psychology-related features on predicting class 3. Answering yes (1 at x-axis) to “Loneliness, isolation” is more likely (higher SHAP value) to be classified as class 3 than answering no (0 at x-axis). (Coding -1: Do not know. Coding -3: Prefer not to answer)

Conclusions and Discussions

Apart from the most impactful feature — “Seen a psychiatrist for nerves, anxiety, tension or depression” — there are several other features related to mental state that impact the prediction of ADRs, e.g. “Loneliness, isolation” (Figure 4b). This suggests a possible link between ADRs and an individual’s mental state. This connection was proposed and studied before, with no clear conclusion. Our study, powered by the large sample size of UKB and the decision-tree boosting algorithm and SHAP analysis, provides some evidence for a seemingly unlikely factor that should be considered in ADR prediction. Other noteworthy impactful features are biomarkers obtained from blood and urine samples. This is not surprising as protein biomarkers have been used for diagnosing ADRs. Our selected cohorts and the LightGBM modeling suggested that both physiology and psychology are linked to ADRs. SHAP analysis provides possible explanations for the impacts and the dynamic among factors.

We did not find strong individual SNPs and Indels for predicting ADRs. There are several potential reasons for the subpar performance of genomic-based modeling:

  1. Many of the genetic variants that have been demonstrated to be associated with ADRs are located in regions of the genome for which SNParrays are ill suited, such as the incredibly complex HLA regions and the CYP gene family.
  2. SNParrays were designed to investigate common variants, and therefore their data are less informative for rare and low-frequency variants. We used relatively stringent criteria when selecting genomic features to filter them out. But it is possible that some rare and low-frequency variants are informative for ADR prediction.

The fact that phenotypic features are stronger predictors than genomic ones is not surprising. Published GWAS results have found that the genetic component very often only accounts for a small fraction of variance observed in the trait of interest (so called “missing heritability”) meaning that a lot of non-genetic factors also play a role in the said trait.

Future directions/closing thoughts

This study is the first time we applied a decision-tree boosting approach to PGx. There is a lot of room for improvement. In addition to relaxing the selection criteria for genomic variants, we can also explore the use of Identity By Descent (IBD) segments as genomic features. IBD segments consist of multiple variants, and they capture genomic architecture in people, without being subject to sampling error (selecting genomic variants). Previously, detecting IBD segments was computationally intensive, however, recent advancements in methodology have made it efficient for large scale studies.

Moreover, studies on genomic variants for drug-specific ADRs have yielded encouraging results. We may try a narrower-scope approach, such as separating individuals into groups that take the same type of medications (sharing the same target or biological system).

Acknowledgements

This research has been conducted using the UK Biobank Resource under application number 46926.

--

--