<figure><img src="image_1.png" /></figure>

HAL
open science

# Detecting adverse high-order drug interactions from individual case safety reports using computational statistics on disproportionality measures

Jules Bangard, Einar Holsbø, Kristian Svendsen, Vittorio Perduca, Etienne Birmelé

► To cite this version:

Jules Bangard, Einar Holsbø, Kristian Svendsen, Vittorio Perduca, Etienne Birmelé. Detecting adverse high-order drug interactions from individual case safety reports using computational statistics on disproportionality measures. 2025. hal-05012416

HAL Id: hal-05012416

https://hal.science/hal-05012416v1

Preprint submitted on 31 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Detecting adverse high-order drug interactions from
individual case safety reports using computational
statistics on disproportionality measures

Jules Bangard¹, Einar Holsbø², Kristian Svendsen³, Vittorio Perduca⁴,
Étienne Birmelé¹

¹Institut de Recherche Mathématique Avancée,
UMR 7501 Université de Strasbourg et CNRS 7 rue René-Descartes,
67000 Strasbourg, France

²Department of Computer Science,
Faculty of Science and Technology, UiT-The Arctic University of Norway,
PO, Box 6050 Stakkevollan, N-9037 Tromsø, Norway

³Department of Pharmacy, Faculty of Health Sciences,
UiT the Arctic University of Norway, Tromsø, Norway

⁴CNRS, MAP5, Université Paris Cité, F-75006 Paris, France.

March 31, 2025

Abstract

Adverse drug interactions are a critical concern in pharmacovigilance, as both
clinical trials and spontaneous reporting systems often lack the breadth to detect
complex drug interactions. This study introduces a computational framework for
adverse drug interaction detection, leveraging disproportionality analysis on individual
case safety reports. By integrating the Anatomical Therapeutic Chemical classification,
the framework extends beyond drug interactions to capture hierarchical pharmacological
relationships. To address biases inherent in existing disproportionality measures, we
employ a hypergeometric risk metric, while a Markov Chain Monte Carlo algorithm
provides robust empirical p-value estimation for the risk associated to cocktails. A
genetic algorithm further facilitates efficient identification of high-risk drug cocktails.
Validation on synthetic and FDA Adverse Event Reporting System data demonstrates
the method's efficacy in detecting established drugs and drug interactions associated
with myopathy-related adverse events. Implemented as an R package, this framework
offers a reproducible, scalable tool for post-market drug safety surveillance.

**Keywords**– Disproportionality Analysis; Drug-Interactions; Genetic Algorithm; Markov Chain Monte Carlo; Pharmacovigilance; Spontaneous Reporting Systems
# 1 Introduction

There are inherent limitations of clinical trials (RCTs) done before a medication is authorised in terms of cost, surveillance time, size and a lack of diversity of patients included (e.g. that pregnant woman, children and elderly often are not included) [1]. Another issue is that RCTs are limited in their ability to detect rare adverse drug reactions (ADRs) and potential drug interactions (DDIs) since they usually only include a few thousand participants at most. Moreover, the use of multiple medications is often an exclusion criterion in RCTs, leaving gaps in the assessment of polypharmacy risks. As a result, tested drug combinations tend to be limited to those where physicians and pharmacologists have existing clinical experience or hypotheses about potential interactions [2]. This means that monitoring of the safety of medications and the potential for DDIs after marketing is essential. This monitoring is often referred as post-authorization pharmacovigilance.

Elderly individuals are commonly prescribed multiple medications due to the increasing prevalence of comorbidities associated with aging. One study suggested that patients aged 75 or older in Austria took on average 7.5 drugs [3]. Among this demographic, 10% of hospitalizations were attributed to ADRs, and 18.7% of these ADRs were caused by a combination of more than one drug, a drug-drug interaction. More broadly, a meta-analysis attributed approximately 22.2% of hospitalizations caused by adverse drug reactions to DDIs [4].

Pharmacovigilance plays a key role in identifying adverse DDIs that went undetected during the pre-clinical and clinical trial phases. One key part of post-marketing pharmacovigilance is spontaneous reporting systems (SRSs), such as the Food & Drug (FDA) Adverse Event Reporting System (FAERS) [5]. Reports contained in SRSs are called Individual Case Safety Reports (ICSRs), and the content of these is harmonized through international guidance delivered by the International Council on Harmonisation (ICH E2B) [6]. SRSs always contain information about the ADR experienced and medications suspected of causing these reactions. However other information such as simultaneous use of other medication by patients is not mandatory to report and is often lacking.

Disproportionality analysis (DA) encompass widely used methods in pharmacovigilance to detect ADRs by assessing the frequency of adverse events associated with specific drug consumption relative to what would be expected. Most established disproportionality methods like the proportional reporting ratio (PRR) [7], reporting odds ratio (ROR) [8], Bayesian Confidence Propagation Neural Network (BCPNN) [9], and the Gamma Poisson Shrinker (GPS) [10] allow for the identification of signals that warrant further investigation in the context of single drug assessment. Fewer DAs methods are available to use for multiple drugs consumption, such as the Ω shrinkage method [11]. This method allows the use of a DA measure in order to detect sets of the type {Drug-Drug-Adverse Event} and stands as the standard measure when it comes to detecting the interactions of two drugs. PRR have also been adapted to DDIs [12]. The chi-square method [13] allows to control the false positive rate when few patients represent the event in the database.

Other methods than DA measures exist in order to detect DDIs. Among them, well known methods encompass logistic models [14, 15] and association rules methods [16, 17]. One can refer to multiple reviews for a more detailed overview of existing methods [18, 19].

In this study, a novel framework is presented that extends traditional pharmacovigilance
methods. The approach enables detection of DDIs associated to a given adverse event for drug cocktails of arbitrary size, moving beyond the typical pairwise interaction analysis by exploring higher order interactions. Notably, the Anatomical Therapeutic Chemical (ATC) classification tree is leveraged to explore pharmacological relationships between families of drugs at different hierarchical levels, offering a broader perspective beyond simple drug interactions. A key contribution is the ability to compute empirical p-values for risk scores, ensuring more robust thresholds for screened signals compared to classical DAs methods. An R package is provided to implement these methods, facilitating further research, reproducibility, and practical applications in pharmacovigilance.

# 2 Methods

## 2.1 ATC Tree and Cocktail Definition

Multiple classifications of drug active ingredients exist, one being the ATC classification. This system is structured hierarchically and can be represented as a tree with five levels containing, at the time, 6809 nodes. The leaves of this tree are the active ingredients, while the first level consists of nodes representing organs or systems affected by the descendant active ingredients. The remaining nodes correspond to therapeutic or pharmacological families [20].

By applying a Depth-First Search algorithm to enumerate each node of the ATC tree, we can define a drug cocktail as a set of integer corresponding to specific nodes in the tree. More formally, a cocktail $C$ is defined by $C = (x_1, ..., x_k) \in \Delta^k$ where $k \in \mathbb{N}$, and $\Delta$ represents the set of numbered nodes of the ATC tree $T$. Figure 1 shows examples in a simplified tree, the green nodes denoting the considered cocktail.

Considered cocktails can include internal nodes of the tree (thus representing families of active ingredients), allowing for the detection of more general signals. For example, paracetamol might send a weak signal, while if we move up the tree, analgesics as a whole could represent a stronger and more general signal. Therefore, all patients taking at least one drug from this drug family will be considered in the computation of the risk.

## 2.2 Cocktail Risk Characterization

In the field of pharmacology, accurately characterizing the risk associated with drug administration is a hard task.

Highly used methods to evaluate drug interactions, like the PRR [7,12] are biased toward drugs or drugs cocktails taken by few patients, assigning high risk to them (also see Figure 2). In order to mitigate this phenomenon, we propose to use the hypergeometric distribution to evaluate the hazard of medication cocktails. Consider a dataset of $N$ patients, among which $K$ experience the adverse event AE. Let $n_C$ be the number of people taking a cocktail $C$ and $x$ the number of patients taking $C$ and experiencing AE. We define the risk to be

$$H(C) = -\log(\mathbb{P}(X \geq x))$$

where $X \sim H(n_C, K, N)$. This measure is the log p-value under the null hypothesis that the number of people with AE in a uniform sample of $n_C$ out of $N$ people follow the
hypergeometric distribution. It has the advantage to take into account the relative number of patients taking the cocktail *C*. For example, $H(C)$ will be greater for a cocktail taken by 100 patients, among whom 50 experience AE, than for a cocktail taken by 10 patients, among whom 5 experience AE. In contrast, the PRR [7] makes no distinction between the two.

Such measures are commonplace in bioinformatics, especially when working on functional enrichment analysis [21].

## 2.3 Cocktail Identification

### 2.3.1 High-risk Drug Cocktails Identification

Due to combinatorial complexity, obviously it is not possible to compute $H(C)$ for each possible cocktail $C$. Instead, to explore the space of cocktails and locate those at high-risk of AE, we use a genetic algorithm [22]. The main idea is to simulate the evolution of a population of cocktails according to the principle of natural selection. The goal is to search for the most performing individuals with respect to an evaluation criterion based on $H$. The steps of the algorithm are the following, the algorithm repeating them until a user-defined number of iterations is reached.

**Initialization.** The genetic algorithm's population consists of $m$ cocktails. These cocktails are randomly initialized and can vary in size.

**Evaluation and selection.** At each iteration, the population undergoes an evaluation and selection phase. The evaluation computes the score $H(C)$ for each cocktail $C$ in the population. In order to avoid the convergence of the population to a single high-score cocktail, scores of similar cocktails are penalized as further explained in Section 2.3.2. The scores enable tournament selections, where $d$ individuals are drawn from the population, and the one with the highest $H(C)$ among these $d$ is retained for the reproduction phase; where $d$ is a hyperparameter. Such tournaments are conducted until the desired number of individuals for reproduction is reached.

**Modification and population replacement.** The evolution of the population occurs in two stages, each aiming to improve performance according to the evaluation criterion. First, a crossover operation allows two cocktails to exchange information. Here, the crossover involves exchanging sub-trees between two cocktails as follows:

* an internal node $v$ of the tree is randomly selected;

* the nodes of the subtree rooted at $v$ are exchanged between the two sequences.

After performing the crossover, a mutation is applied to the resulting individuals, chosen from two types. The first type is a local mutation which changes a randomly selected node of the cocktail to one of its free neighboring nodes. This mutation is further explained in Section 2.4. The second type is an addition/deletion mutation which operates as follows, with $k$ being the cocktail length and $\alpha$ a chosen hyperparameter:
• with probability min(1, α/k ), a node uniformly drawn from Δ is added to the sequence;

• with probability 1 − min(1, α/k ), a uniformly drawn node from the cocktail is removed.

An example of crossover and mutations can be seen on Figure 1 (a), (c) and (d).

<figure><img src="image_3.png" /><figcaption>Figure 1: Cocktail modifications used for the genetic algorithm (a, c and d) and the MCMC algorithm (b and c). Green nodes are part of considered cocktails. In Crossover, the orange node represents the selected internal node whose subtrees are being swapped. In local mutations, orange nodes represent legal moves.</figcaption></figure>

2.3.2 Distance between drug cocktails definition and cocktail penalisation

When a high-scoring cocktail is identified, the genetic algorithm may converge to a population consisting entirely of that cocktail. To avoid this uniformization phenomenon, similar cocktails are penalized in the evaluation phase as follows,

$$
H_{pen}(C) = \frac{H(C)}{\sum_{C_i \in C} \mathrm{Sim}(C, C_i)}
$$
The computation of the similarity $Sim(C, C')$ is based on a distance inspired by the Levenshtein distance [23]. However, unlike the traditional Levenshtein distance, sequences are treated as unordered sets. For two drug cocktails, $C_1$ and $C_2$, of sizes $n_1$ and $n_2$, the distance $d(C_1, C_2)$ is defined as the minimal cost required to transform $C_1$ into $C_2$ using three elementary operations.

* $Ins_a(C)$ consists of adding $a$ to the cocktail $C$.

* $Del_a(C)$ consists of deleting $a$ from the cocktail $C$.

* $Sub_{a,b}(C)$ consists of substituting $a \in C$ by $b$.

An associated cost is defined for each operation.

**Substitution.** The cost associated with the substitution operation is chosen to be consistent with the conceptual similarity of cocktails. If $a$ is a drug belonging to $C$, the cost should increase as the drug $b$ diverges further from drug $a$. For example, if $a$ is a drug, and $b$ a drug family that contains $a$, the cost should be moderate. Conversely, if $b$ is a drug family not containing $a$ the cost should be higher. This distance is thus defined by the maximal distance between $a$ and $b$ to their Lowest Common Ancestor.

**Insertion, Deletion.** The deletion and insertion cost are chosen as $\frac{depth(T)}{2}$. This choice implies that a substitution always costs less than a deletion followed by an insertion. The latter are used only when the two cocktails do not have the same length.

A transformation $f$ from $C_1$ to $C_2$ is a composition of elementary operations that go from $C_1$ to $C_2$. The associated cost $cost(f)$ is defined as the sum of the cost of the operations used in $f$. Finally, $d(C_1, C_2) = \min_{f} cost(f: f(C_1) = C_2)$.

Finally, the maximum distance between $C_1$ and $C_2$ being $(n_1 + n_2)\frac{depth(T)}{2}$, we define the similarity as

$$Sim(C_1, C_2) = 1 - \frac{2D(C_1, C_2)}{(n_1 + n_2)depth(T)}$$

The computation of the similarity is achievable in $O(n \times m \times depth(T) + |\Delta|)$ operations in the worst case, where $n$ and $m$ denotes the size of the cocktails $C_1$ and $C_2$, respectively. The algorithm to compute the similarity is detailed in Supplementary Materials (SM).

### 2.3.3 Output clustering

Despite the diversity mechanisms integrated into the genetic algorithm, not all drug cocktails generated within the population are unique. It is moreover common to encounter solutions that are merely variations of others, differing only by transformations such as changing a node in the tree to its parent or child. To streamline analysis and enhance the efficiency, a post-treatment clustering of similar solutions is implemented. This method allows to focus on the most risky cocktails within each cluster or to interpret pharmaceutically clusters rather than individual cocktails.
To do so, drug cocktails are embedded into a two-dimensional space using the UMAP algorithm [24], which aims to preserve similarity in the latent space. This representation enables the effective use of conventional machine learning clustering algorithms in $\mathbb{R}^2$. Specifically, the DBSCAN algorithm [25] is applied to identify clusters of similar drug cocktails with an intuitive way of choosing hyperparameters.

## 2.4 P-Value Assignment for Drug Cocktails

Once a list of high-score cocktails has been found, an important step of the analysis is the attribution of a p-value to their scores to decide if they are significant. Such a p-value measures the probability that a randomly chosen cocktail has a score greater than the observed score. However, a naive sampling of cocktails leads to essentially only cocktails taken by no patient. We therefore use a Metropolis-Hastings MCMC algorithm to sample under a null distribution of scores of cocktails of a given size. Indeed, such a method can be used by conditioning on the fact that all visited cocktails are present in the dataset [26].

To employ such an algorithm, it is necessary to define a state space $C = \{C_1, \dots, C_p\}$, a computable target measure $f(C_i)$, and conditional laws $q(\cdot|C_i)$ under which simulation is possible and new states can be proposed.

**State set.** The state set is made of all cocktails of $k$ drugs for a fixed $k$.

**Proposal law.** The proposal law is defined as a mixture of two mutation laws of the current cocktail. They operate as follows:

* Random mutation consists of a completely random movement in the cocktail space.

* Local mutation involves a local movement relative to the structure of the drug tree. Here, a node $x_p$ of the state $C_i$ is changed to one of its free neighboring nodes.

At each iteration, the random and local mutations have probability $p_R$ and $1 - p_R$, where $p_R$ is a hyperparameter.

Figure 1 (b, c) presents examples of a random and a local mutations.

**State evaluation.** The evaluation of a drug cocktail is based on the score $H(C)$. The chosen target measure is then:

$$f_T(C_i) = \frac{1}{Z(T)} \times e^{\frac{H(C_i)}{T}}$$

where $Z(T) = \sum_C e^{\frac{H(C)}{T}}$. $T$ is a parameter known as temperature, which modulates space exploration by more readily accepting cocktails with moderate scores (high $T$) or, conversely, by strongly favoring combinations of drugs with high scores (low $T$).

The acceptance probability of cocktail $C_{i+1}$ from cocktail $C_i$ is given by:

$$\min \left( 1, \frac{f_T(C_{i+1})}{f_T(C_i)} \times \frac{q(C_i|C_{i+1})}{q(C_{i+1}|C_i)} \right)$$
The theory related to the Metropolis-Hastings algorithm [27] ensures that the empirical distribution of $f_T(C_i)$ for the constructed cocktail chain converges to the distribution of $f_T(C)$. A very long realization of such a walk, therefore, allows for the approximation of the distribution that can be used to determine an empirical p-value for the score of a cocktail of interest. This enable the possibility to say whether or not a high-risk is truly significant (e.g. among the top 5% of scores).

## 2.5 Datasets

### 2.5.1 Simulated data

Multiple datasets were simulated to evaluate the method performance against known outcomes. The datasets, designed to challenge the algorithm, simulate various patient scenarios. Each patient record includes prescribed medications and the corresponding occurrence of an adverse event AE.

The first dataset is composed of 200,000 patients, and has the following characteristics:

* 1% of the patients take a size-3 drug cocktail $C_1$ and have a $\frac{1}{100}$ chance of having AE.

* 1% take a size-3 drug cocktail $C_2$ and have a $\frac{1}{200}$ chance of having AE.

* 1% take a size-2 drug cocktail $C_3$ and have a $\frac{1}{100}$ chance of having AE.

* 1% take a size-2 drug cocktail $C_4$ and have a $\frac{1}{200}$ chance of having AE.

A small percentage of the dataset (1.5% per combination) are combinations of 2 out of the 3 drugs from $C_1$ and $C_2$, but with no risk of AE. This helps to mitigate the false identification of sub-cocktails of $C_1$ and $C_2$ as high risk cocktail because those who take two drugs of these cocktails will almost surely take the remaining drug of the cocktail.

The remaining 87% of the datasets consists of patients assigned with random cocktails drawn uniformly. The size $s$ of each cocktail is drawn according to a Poisson distribution with $\lambda = 4$ (mean size of drugs cocktail taken by patients in the dataset). For each cocktail, $s$ nodes of the ATC tree are selected uniformly, with each combination assigned an adverse event with probability $\frac{1}{15000}$.

Three others datasets were similarly constructed, the difference lying in the size of the cocktails inducing AE. One has only size-two cocktails, other only size-three cocktails and the last one, size-two, three and four cocktails. Reader may refer to the SM for more details on these three datasets.

### 2.5.2 FAERS Data

The method was also assessed using the FAERS dataset, which consists of ICSRs submitted by healthcare professionals, consumers, and manufacturers. These reports include details on patient drug intake and the side effects experienced. We deployed our methods on FAERS data from the second quarter of 2013 to the second quarter of 2015.

Significant refinement of the FAERS data was required. Duplicate reports were removed, retaining only the report with the most recent ID. Subsequently, a link was established
between the prescribed drugs and the ATC codes of each active ingredient. This process involved matching drug names to their respective active ingredients and converting these ingredients to their corresponding ATC codes. The DiANA dictionary [28] facilitated the standardization of FAERS drug names to ATC codes. Reports with unmatchable drug names in the DiANA dictionary were excluded, reducing the dataset from 2,043,231 to 1,612,931 patients.

For this study, we focused on myopathy as the selected adverse event outcome. It is a clinically concerning condition with a sufficient number of reported cases in the dataset (536 cases). To validate our results, we compared the identified drug-myopathy associations with known drugs already established to cause myopathy [29–31]. Code for data refinement is available on GitHub (See SM).

# 3 Results

## 3.1 Score Comparison

To support the use of the hypergeometric score $H(C)$ introduced in Section 2.2, we compared it to other risk scores on simulated data. Figure 2 illustrates the performance of various risk scoring methods for detecting high-risk drug combinations, including the Relative Risk (RR), PRR, our Hypergeometric Score $H(C)$, CSS, and Omega Shrinkage. Each subplot displays the score values for cocktails representing true solutions (green) and cocktails not representing true solutions (red). The bottom-right panel presents the Precision-Recall (PR) curve, comparing the detection power of the scores in identifying high-risk cocktails. PRR is a special score because its value is either one or zero, which allows the computation of only one value of precision and recall. PRR is therefore represented by a single point rather than a PR curve.

The scores that assigned higher risk to solutions known to give the AE (true solutions) were the hypergeometric score and the Omega Shrinkage measure. Other scores appeared to assign higher risk to cocktails that are not true solutions. It must be emphasized that the Omega Shrinkage measure is difficult to compare using a threshold, as the original article [11] suggests signaling a cocktail when the score exceeds zero. However, in the simulations, the computed score never exceeded zero, as depicted by the x-axis of its jitter plot. Despite these considerations, these two methods are the least biased toward cocktails taken by only a few patients.

Moreover among the listed scores, the hypergeometric score is the only one that can be easily generalized to larger sized cocktails. Another method has been proposed in order to compute scores of larger size cocktails [32], but the computational cost constrains the authors to limit the study to the 20 most distributed drugs.

## 3.2 Application to the Simulated Dataset

### 3.2.1 Estimation of Risk Distribution

Risk distribution was estimated for size-two drug cocktails on Section 2.5.1 dataset. Estimation of risk distribution for higher cocktails sizes are possible but it is nearly impossible to compare
<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison of scores for cocktail of size-two. Each dot denotes the risk of a cocktail computed on the formerly presented synthetic dataset (Section 2.5.1). Red dots represent cocktails that do not induce adverse event, green ones represent cocktail inducing adverse event. The bottom-right corner shows the Precision-Recall curves for each score. The perfect classification corresponds to the upper right corner. The areas under the curves allow to compare different methods.</figcaption></figure>
it to the true distribution as it is computationally prohibitive to obtain. The distribution
estimated by the MCMC algorithm, is compared to the true risk distributions in Figure 3.

The left panels display the distributions of risk scores for both the estimated (top-left) and true (bottom-left) risk values. Both distributions share a similar shape, with the majority of risk scores concentrated at low values, as 95% of the scores fall below 11. However, some differences can be observed, particularly in the tail of the distribution, where cocktails of higher risks are under-represented in the estimated distribution.

<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison of estimated and true risk distributions for size-two drug cocktails. Left panels show comparison of risk distribution among size-two cocktails, right panels allows to compare probabilities and quantiles of both distributions</figcaption></figure>

The right panel of Figure 3 presents a Probability-Probability (PP) plot (top) and a Quantile-Quantile (QQ) plot (bottom), comparing the quantiles and probabilities of the estimated and true risk distributions. While the right panel of the figure demonstrates good agreement at lower risks, where most of the data lie, deviations at higher risk values suggest that the estimated distribution slightly underrepresents the risk for more extreme values.

Results indicate that the method performs well in estimating risk scores for the majority of cocktails, capturing the overall risk distribution with reasonable accuracy. The slight underestimations which are present for high-risk cocktails are not a problem since the interest of the method is to assign p-values. P-values are still reliable as shown in the PP-plot, Figure 3. Consequently, we can assign robust empirical p-values to any size-two cocktails based on their risks.
<figure><img src="image_6.png" /><figcaption>Figure 4: Clustering of high-risk drug cocktails identified by the genetic algorithm on the simulated dataset. True cocktails solutions are the center of square shape</figcaption></figure>

3.2.2 Genetic algorithm output and clustering

The genetic algorithm was applied to the simulated dataset to identify high-risk drug cocktails.
Multiple runs of the algorithm were conducted using different hyperparameter sets to ensure
robustness, and the results were subsequently combined to create a comprehensive list of
high-risk cocktails.

The genetic algorithm successfully identified nearly all high-risk size-two and size-three
drug cocktails in the simulated dataset (Figure 4). For size-two cocktails, the algorithm
consistently found the exact high-risk combinations. However, for size-three cocktails, the
algorithm sometimes identified cocktails that were very close to the true high-risk combinations,
missing only one drug from the correct cocktail in a few cases (oftentimes, choosing parent
nodes instead of the actual drugs).

To streamline the analysis of the large set of results, significant cocktails were filtered using
the empirical p-value by setting a threshold of 5%. Clustering techniques were then applied
to group similar cocktails together. As discussed in Section 2.3.3, the UMAP algorithm has
been used for dimensionality reduction, followed by the DBSCAN clustering method. This
post-processing step allows to reduce redundancy by grouping cocktails that differed only
slightly, such as by substituting a drug for another within the same pharmacological family.
Such cocktail would have similar medical interpretations.

Figure 4 illustrates the results of the clustering process. Each point represents a drug cocktail, and different clusters are colored distinctly to highlight groups of similar cocktails. The clustering of solutions effectively captured the diversity of the drug cocktails identified
by the genetic algorithm. We can see that the algorithm effectively found expected solutions or very similar solutions. In addition, the number of clusters corresponding to cocktails that are not supposed to induce an adverse event is limited, as only clusters 3, 4 and 5 correspond to cocktails that are not meant to be found. It is possible to argue that the second cluster may be further split into two distinct clusters, one corresponding to the real solution while the other would correspond to no real solutions.

## 3.3 Application to the FAERS Spontaneous Reporting Data

### 3.3.1 Estimation of Risk Distribution

The risk estimation method was applied to the FAERS spontaneous reporting dataset presented in Section 2.5.2. Figure 5 presents a comparison between the estimated risk distribution and the true risk distribution for size-two drug cocktails. The left panel of the figure shows the distributions of risk scores, while the right panel presents the QQ plot and the PP plot, comparing the quantiles and probabilities of both distributions.

<figure><img src="image_7.png" /><figcaption>Figure 5: Comparison of estimated and true risk distributions for size-two drug cocktails on the FAERS dataset presented Section 2.5. Left panels show comparison of risk distribution among size-two cocktails, right panels allows to compare probabilities and quantiles of both distributions</figcaption></figure>

The histogram reveals that the estimated distribution aligns well with the true distribution for the majority of cocktails with lower risk scores. However, deviations begin to emerge in the tail of the distribution. Specifically, 15 of the riskiest cocktails in the true distribution were not captured by the estimated distribution as we see in the QQ plot. This explains the observed shift in the QQ plot at the highest quantiles since higher risk cocktails have not been found by the MCMC algorithm, highlighting the need of a complementary method like the genetic algorithm in order to find riskiest cocktails.
Despite this slight deviation, the empirical p-values remain robust for both lower and higher-risk cocktails as shown by the PP-plot in the Figure 5.

### 3.3.2 Genetic algorithm output and clustering

The genetic algorithm was applied on the FAERS data focusing on the myopathy AE by running 180 parallel executions of the genetic algorithm on varying population sizes (from 100 to 1000 cocktails per generation). The whole procedure run in less than 8 hours on a 24-core server.

Cocktail sizes present in the merged final populations vary from 1 to 6. The MCMC algorithm was run for each cocktail size in this range to assign an empirical p-value to each solution. All solutions having a p-value lower than 0.05 were kept. No multiple testing correction was made in order to avoid false negatives, that is disregarding interesting cocktails, even if this may inflate the number of false positives. 682 drug cocktails composes the final list. The previous clustering procedure was applied, leading to 15 clusters.

To validate the results, some drugs or drug families known to be linked to myopathy adverse events were considered, that is hypolipemic drugs, Colchicine, corticosteroids, Ciclosporine, Beta blocking agents, Fluoroquinolones, and anti-malaria drugs [29–31]. The first 150 solutions were tagged with the families they correspond to, if any.

The final result is a table of 682 rows, one per selected cocktail, indicating its composition, the number of patients taking it, the number of patient taking it and facing myopathy, hypergeometric and RR score. It also contains the cluster it belongs to and the tagged families. The entire table is available as a supplementary file. Table 1 summarizes the cluster assignments for the top 150 cocktails identified in the genetic algorithm. Among these, eleven cocktails could not be assigned to any drug or drug family listed in the table headers. Note that a cocktail can be associated with more than one drug or drug family.

<table><caption>Table 1: Summary of clustering applied to the identified solutions. Each row corresponds to a cluster. The 150 cocktails with the highest risk were analyzed. Box (i, j) in the table represents the number of cocktails in cluster i that include the drug or drug family j.</caption><thead><tr><th>Cluster</th><th>Hypolipemic Drugs</th><th>Colchicine</th><th>Steroids</th><th>Cyclosporine</th><th>Beta blocking agents</th><th>Domperidone</th><th>Fluoroquinolones</th></tr></thead><tbody><tr><td>1</td><td>31</td><td>0</td><td>3</td><td>0</td><td>0</td><td>0</td><td>4</td></tr><tr><td>2</td><td>0</td><td>11</td><td>0</td><td>4</td><td>0</td><td>0</td><td>0</td></tr><tr><td>3</td><td>7</td><td>0</td><td>0</td><td>0</td><td>19</td><td>0</td><td>0</td></tr><tr><td>4</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>50</td><td>0</td></tr><tr><td>5</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>19</td><td>0</td></tr><tr><td>6</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>4</td></tr></tbody></table>

This table highlights the ability of the clustering method to group cocktails with similar pharmacological interpretations and showcase the method's capability to detect known signals from the FAERS dataset. Specifically, cluster 1 corresponds to cocktails containing Hypolipemic Drugs, cluster 2 to Colchicine and Ciclosporine, cluster 3 to beta blocking agents and cluster 6 to Fluoroquinolones. These clusters align with previously reported pharmacological associations [29]. Similarly, clusters 4 and 5 correspond to Domperidone-containing cocktails which have been associated with adverse effects, including potential cardiac complications, as reported by the British Medicines and Healthcare products Regulatory Agency [33].
Interestingly, additional drug families reported in [29], such as anti-malarial drugs, were also identified but at later ranks. For instance, anti-malarial drugs are primarily grouped in cluster 9, as shown in the complete solution table.

# 4 Conclusion

As co-medication becomes increasingly common, there is a growing need for methods capable of detecting signals of harmful drug combinations from the available large databases. The proposed method addresses this need by identifying signals and assigning them a p-value using a hypergeometric disproportionality analysis measure. Additionally, the method enables the identification of broader signals within the ATC hierarchy by proposing "cocktails" of not only active substances but also chemical, therapeutic, and anatomical families, leveraging the hierarchical classification of active substances.

Application on synthetic datasets demonstrated that using the hypergeometric score reduces the false positives from cocktails taken by a small number of patients, enhancing the robustness of the measure. The results on these datasets of our MCMC algorithm to estimate the distributions of cocktail risks were encouraging, as the estimated distributions closely aligned with the true distributions, indicating reliable p-value assignment. Furthermore, the genetic algorithm effectively identified the majority of the harmful cocktails, with a high success rate, highlighting its efficiency in navigating the solution space.

Applying this method to previous FAERS data for the myopathy adverse event yielded promising results. A literature review confirmed the intersection between the identified signals and drugs known to have a higher likelihood of causing myopathy, demonstrating the effectiveness of the proposed methodology. Results also indicate that certain drug combinations are more strongly associated with myopathy than individual medications. Notably, the cyclosporine/colchicine combination exhibits a higher hypergeometric score (73.1 vs. 63.4) and PRR (879.4 vs. 57.1) compared to colchicine alone, reinforcing the importance of analyzing drug interactions in pharmacovigilance. This combination is known to be more likely to induce myopathy [34].

Furthermore, our approach identified a size-four drug cocktail (metformin, prasugrel, bisoprolol, simvastatin) associated with an increased risk signal (hypergeometric score = 72.2, RR = 3060.6). This combination was observed in nine patients, all of whom experienced the adverse event. While this finding demonstrates the feasibility of detecting higher-order drug interactions, we emphasize that no clinical validation is currently available for this specific combination. This underscores both the potential of the method in identifying complex drug interactions and the need for further validation through complementary studies.

This approach can also be extended to other settings where adverse events are explored. For instance, the method could be applied using the ICD diagnosis classification or the MedDRA system, both of which are hierarchical classifications. Such an application would facilitate the identification of symptoms associated with the consumption of drug combinations.

The proposed method is implemented as an R package *emcAdr*, available on GitHub (see SM) and on the CRAN [35].
# Funding

This work was partially supported by the “PHC AURORA” programme (project number: 49704QC), funded by the French Ministry for Europe and Foreign Affairs, the French Ministry for Higher Education and Research and the Norwegian Council for Research.

# References

[1] Robert William Sanson-Fisher, Billie Bonevski, Lawrence W Green, and Cate D'Este. Limitations of the randomized controlled trial in evaluating population-based health interventions. *American journal of preventive medicine*, 33(2):155-161, 2007.

[2] Peter GM van der Heijden, Eugène P van Puijenbroek, Stef van Buuren, and Jacques W Van der Hofstede. On the assessment of adverse drug reactions from spontaneous reporting systems: the influence of under-reporting on odds ratios. *Statistics in medicine*, 21(14):2027-2044, 2002.

[3] Jochen Schuler, Christina Dückelmann, Wolfgang Beindl, Erika Prinz, Thomas Michalski, and Max Pichler. Polypharmacy and inappropriate prescribing in elderly internal-medicine patients in austria. *Wiener klinische Wochenschrift*, 120, 2008.

[4] Supinya Dechanont, Sirada Maphanta, Bodin Butthum, and Chuenjid Kongkaew. Hospital admissions/visits associated with drug-drug interactions: a systematic review and meta-analysis. *Pharmacoepidemiology and drug safety*, 23(5):489-497, 2014.

[5] U.S. Food and Drug Administration. Fda adverse event reporting system (faers) latest quarterly data files. https://www.fda.gov/drugs/questions-and-answers-fdas-adverse-event-reporting-system-faers/fda-adverse-event-reporting-system-faers-latest-quarterly-data-files, 2024. Accessed: 2024-10-07.

[6] International Council for Harmonisation of Technical Requirements for Pharmaceuticals for Human Use (ICH). E2b(r3) individual case safety report (icsr) specification and related files. https://ich.org/page/e2b3-individual-case-safety-report-icsr-specification-and-related-files, 2024. Accessed: 2024-10-07.

[7] Stephen JW Evans, Patrick C Waller, and S Davis. Use of proportional reporting ratios (prrs) for signal generation from spontaneous adverse drug reaction reports. *Pharmacoepidemiology and drug safety*, 10(6):483-486, 2001.

[8] Eugene P. van Puijenbroek, Andrew Bate, Hubert G. M. Leufkens, Marie Lindquist, Roland Orre, and Antoine C. G. Egberts. A comparison of measures of disproportionality for signal detection in spontaneous reporting systems for adverse drug reactions. *Pharmacoepidemiology and Drug Safety*, 11, 2002.
[9] Andrew Bate, Marie Lindquist, I Ralph Edwards, Sten Olsson, Roland Orre, Anders Lansner, and R Melhado De Freitas. A bayesian neural network method for adverse drug reaction signal generation. *European journal of clinical pharmacology*, 54:315–321, 1998.

[10] William DuMouchel. Bayesian data mining in large frequency tables, with an application to the fda spontaneous reporting system. *The American Statistician*, 53(3):177–190, 1999.

[11] G Niklas Norén, Rolf Sundberg, Andrew Bate, and I Ralph Edwards. A statistical methodology for drug–drug interaction surveillance. *Statistics in medicine*, 27(16):3057–3070, 2008.

[12] Xueying Wang, Lang Li, Lei Wang, Weixing Feng, and Pengyue Zhang. Propensity score-adjusted three-component mixture model for drug-drug interaction data mining in fda adverse event reporting system. *Statistics in Medicine*, 39(7):996–1010, 2020.

[13] Masahiko Gosho, Kazushi Maruo, Keisuke Tada, and Akihiro Hirakawa. Utilization of chi-square statistics for screening adverse drug-drug interactions in spontaneous reporting systems. *European journal of clinical pharmacology*, 73:779–786, 2017.

[14] Eugène P van Puijenbroek, Antoine CG Egberts, Eibert R Heerdink, and Hubert GM Leufkens. Detecting drug–drug interactions using a database for spontaneous adverse drug reactions: an example with diuretics and non-steroidal anti-inflammatory drugs. *European journal of clinical pharmacology*, 56:733–738, 2000.

[15] Eugène P Van Puijenbroek, Antoine CG Egberts, Ronald HB Meyboom, and Hubert GM Leufkens. Signalling possible drug–drug interactions in a spontaneous reporting system: delay of withdrawal bleeding during concomitant use of oral contraceptives and itraconazole. *British journal of clinical pharmacology*, 47(6):689–693, 1999.

[16] Yoshihiro Noguchi, Anri Ueno, Manami Otsubo, Hayato Katsuno, Ikuto Sugita, Yuta Kanematsu, Aki Yoshida, Hiroki Esaki, Tomoya Tachi, and Hitomi Teramachi. A new search method using association rule mining for drug-drug interaction based on spontaneous report system. *Frontiers in pharmacology*, 9:197, 2018.

[17] Heba Ibrahim, Amr Saad, Amany Abdo, and A Sharaf Eldin. Mining association patterns of drug-interactions using post marketing fda’s spontaneous reporting data. *Journal of biomedical informatics*, 60:294–308, 2016.

[18] Heba Ibrahim, A Abdo, Ahmed M El Kerdawy, and A Sharaf Eldin. Signal detection in pharmacovigilance: a review of informatics-driven approaches for the discovery of drug-drug interaction signals in different data sources. *Artificial intelligence in the life sciences*, 1:100005, 2021.

[19] Manfred Hauben. Artificial intelligence and data mining for the pharmacovigilance of drug–drug interactions. *Clinical Therapeutics*, 45(2):117–133, 2023.
[20] Yoshihiro Noguchi, Keisuke Aoyama, Satoaki Kubo, Tomoya Tachi, and Hitomi Teramachi. Improved detection criteria for detecting drug-drug interaction signals using the proportional reporting ratio. *Pharmaceuticals*, 14(1):4, 2020.

[21] Steffen Grossmann, Sebastian Bauer, Peter N Robinson, and Martin Vingron. Improved detection of overrepresentation of gene-ontology annotations with parent–child analysis. *Bioinformatics*, 23(22):3024–3031, 2007.

[22] Alain Pétrowski and Sana Ben-Hamida. *Evolutionary algorithms*. John Wiley & Sons, 2017.

[23] Vladimir I Levenshtein et al. Binary codes capable of correcting deletions, insertions, and reversals. In *Soviet physics doklady*, volume 10, pages 707–710. Soviet Union, 1966.

[24] Leland McInnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction. *arXiv preprint arXiv:1802.03426*, 2018.

[25] Martin Ester, Hans-Peter Kriegel, Jörg Sander, Xiaowei Xu, et al. A density-based algorithm for discovering clusters in large spatial databases with noise. In *kdd*, volume 96, pages 226–231, 1996.

[26] Siu-Kui Au and James L Beck. Estimation of small failure probabilities in high dimensions by subset simulation. *Probabilistic engineering mechanics*, 16(4):263–277, 2001.

[27] Christian P. Robert and George Casella. *The Metropolis–Hastings Algorithm*, pages 267–320. Springer New York, New York, NY, 2004.

[28] Michele Fusaroli, Valentina Giunchi, Vera Battini, Stefano Puligheddu, Charles Khouri, Carla Carnovale, Emanuel Raschi, and Elisabetta Poluzzi. Enhancing transparency in defining studied drugs: The open-source living diana dictionary for standardizing drug names in the faers. *Drug Safety*, pages 1–14, 2024.

[29] Sebastian Miernik, Agata Matusiewicz, and Marzena Olesińska. Drug-induced myopathies: A comprehensive review and update. *Biomedicines*, 12(5):987, 2024.

[30] Mederic M Hall, Jonathan T Finnoff, and Jay Smith. Musculoskeletal complications of fluoroquinolones: guidelines and precautions for usage in the athletic population. *PM&R*, 3(2):132–142, 2011.

[31] Ritu Valiyil and Lisa Christopher-Stine. Drug-related myopathies of which the clinician should be aware. *Current rheumatology reports*, 12:213–220, 2010.

[32] Pengyue Zhang, Lei Du, Lei Wang, M Liu, L Cheng, C-W Chiang, H-Y Wu, SK Quinney, L Shen, and L Li. A mixture dose–response model for identifying high-dimensional drug interaction effects on myopathy using electronic medical record databases. *CPT: pharmacometrics & systems pharmacology*, 4(8):474–480, 2015.
[33] Medicines and Healthcare products Regulatory Agency. Domperidone: risks of cardiac side effects. *Drug Safety Update*, 7(10):A1, May 2014. Accessed: 2025-02-24.

[34] D Ducloux, V Schuller, C Bresson-Vautrin, and JM Chalopin. Colchicine myopathy in renal transplant recipients on cyclosporin. *Nephrology, dialysis, transplantation: official publication of the European Dialysis and Transplant Association-European Renal Association*, 12(11):2389-2392, 1997.

[35] Jules Bangard. *The emcAdr package*, 2025. The Comprehensive R Archive Network. https://cran.r-project.org/web/packages/emcAdr/.