Installation

The MDDC package is available on CRAN and can be installed using the following code. Additionally, the development version can be found on GitHub.

install.packages("MDDC")

We load the MDDC package using the following line:

library(MDDC)

Using Boxplot Method

Our goal is to identify (AE, drug) pairs with abnormally high report counts, specifically those cells with counts significantly exceeding their expected values.

First we perform the analysis using mddc_boxplot(). This function has five argument:

• contin_table: A data matrix of an I × J contingency table with rows representing adverse events and columns representing drugs. We recommend users first check the input contingency table using the function check_and_fix_contin_table().

• col_specific_cutoff: Logical. In step 2 of the algorithm, whether to apply the boxplot method to the standardized Pearson residuals within each drug column (default is TRUE) or to the entire table (FALSE).

• separate: Logical. In step 2 of the algorithm, whether to separate the standardized Pearson residuals for the zero cells and non zero cells and apply boxplot method separately or together. Default is TRUE.

• if_col_cor: Logical. In step 3 of the algorithm, whether to use column (drug) correlation or row (adverse event) correlation. Default is FALSE, indicating the use of adverse event correlation. TRUE indicates the use of drug correlation.

• cor_lim: A numeric value between 0 and 1. Specifies the correlation threshold to select “connected” adverse events in step 3. Default is 0.8.

• coef: A numeric value or a list of numeric values. If a single numeric value is provided, it will be applied uniformly across all columns of the contingency table. If a list is provided, its length must match the number of columns in the contingency table, and each value will be used as the coefficient for the corresponding column.

• num_cores: Number of cores used to parallelize the MDDC Boxplot algorithm. Default is 2.

We now perform the MDDC (boxplot) analysis with the statin49 dataset:

set.seed(42)
test1 <- mddc_boxplot(
  contin_table = statin49,
  col_specific_cutoff = T,
  separate = T,
  if_col_cor = F,
  cor_lim = 0.8,
  coef = 1.5
)

The above function outputs a list with three components:

• boxplot_signal: An I × J data matrix with entries 1 or 0, indicating the signals identified in step 2. A value of 1 indicates signals, 0 indicates no signal.

• corr_signal_pval: An I × J data matrix of p-values for each cell in the contingency table from step 5, when the r_ij values are mapped back to the standard normal distribution.

• corr_signal_adj_pval: An I × J data matrix of the Benjamini-Hochberg adjusted p-values for each cell in step 5. Users can choose whether to use corr_signal_pval or corr_signal_adj_pval, and can set their own p-value threshold (for example, 0.05).

Below, we display the first few rows and columns for each component of test1. We first check the component boxplot_signal:

head(test1$boxplot_signal)[, 1:5]
#>                    Atorvastatin Fluvastatin Lovastatin Pravastatin Rosuvastatin
#> Rhabdomyolysis                1           0          0           1            1
#> Muscle Disorder               0           0          0           0            1
#> Muscle Fatigue                0           0          0           0            0
#> Muscle Haemorrhage            0           0          0           0            0
#> Muscle Necrosis               0           0          0           0            0
#> Muscle Rupture                0           0          0           1            0

This indicates the pairs (Rhabdomyolysis, Atorvastatin), (Rhabdomyolysis, Pravastatin), (Muscle Rupture, Pravastatin), (Rhabdomyolysis, Rosuvastatin), and (Muscle Disorder, Rosuvastatin) are identified as signals in step 2 of MDDC (boxplot). Now we look at the second component corr_signal_pval which shows the p-values of all the cells from step 5:

round(head(test1$corr_signal_pval)[, 1:5], digits = 3)
#>                    Atorvastatin Fluvastatin Lovastatin Pravastatin Rosuvastatin
#> Rhabdomyolysis               NA          NA         NA          NA           NA
#> Muscle Disorder           0.527       0.974      0.178       0.545        0.000
#> Muscle Fatigue            0.556       0.834      0.252       0.061        0.656
#> Muscle Haemorrhage        0.561       0.635      0.521       0.307        0.570
#> Muscle Necrosis           0.534       0.186      0.869       0.835        0.677
#> Muscle Rupture            0.569       0.019      0.999       0.000        0.678

In this output, we observe that the first row, corresponding to the adverse event “Rhabdomyolysis”, does not have associated p-values. This is because, in step 2 of the algorithm, “Rhabdomyolysis” was already identified as an AE signal for Atorvastatin, Pravastatin, Rosuvastatin, and Simvastatin. Consequently, the standardized Pearson residual values for these four drugs were replaced with NA. With only two residual values remaining in the first row, it was not possible to find connected AEs for “Rhabdomyolysis”. Therefore, this adverse event was excluded from the subsequent steps of the analysis. Note that for computing the Pearson correlation in step 3, at least three values are required in the matching positions. Applying a p-value threshold of 0.05, we identify the following pairs as signals by considering AE correlations: (Muscle Rupture, Fluvastatin), (Muscle Rupture, Pravastatin), and (Muscle Disorder, Rosuvastatin).

The third component, corr_signal_adj_pval, provides the Benjamini-Hochberg adjusted p-values. Users can choose whether to use corr_signal_pval or corr_signal_adj_pval and can set their own p-value threshold (for example, 0.05).

set.seed(42)
find_optimal_coef(
  contin_table = statin49,
  n_sim = 1000,
  target_fdr = 0.05,
  grid = 0.1,
  col_specific_cutoff = TRUE,
  exclude_small_count = TRUE
)
#> $coef
#> [1] 2.5 3.2 2.8 2.7 2.6 2.4 2.0
#> 
#> $FDR
#> [1] 0.050 0.044 0.050 0.049 0.047 0.038 0.044

Using the Monte Carlo Method

Next, we introduce another primary function of this package, mddc_mc(), which implements the MDDC (MC) algorithm. This function has the following arguments:

• contin_table: A data matrix of an I × J contingency table with rows representing adverse events and columns representing drugs. We recommend users first check the input contingency table using the function check_and_fix_contin_table().

• quantile_mc: In step 2 of the algorithm, this specifies the quantile of the null distribution obtained via the Monte Carlo (MC) method to use as a threshold for identifying cells with high values of standardized Pearson residuals. The default is 0.95.

• mc_num: The number of Monte Carlo replications to perform in step 2. The default is 10,000.

• exclude_same_drug_class: In step 2, when applying Fisher’s exact test to cells with a count less than six, a 2 × 2 contingency table is constructed. This argument specifies whether to exclude other drugs in the same class as the drug of interest. The default is TRUE.

• col_specific_cutoff: Logical. Specifies whether to apply the MC method to the standardized Pearson residuals of the entire table or within each drug column in step 2. The default is TRUE, indicating column-specific cutoff. FALSE applies the MC method to the residuals of the entire table.

• separate: Logical. In step 2 of the algorithm, indicates whether to separate the standardized Pearson residuals for zero cells and non-zero cells, applying the MC method separately. The default is TRUE.

• if_col_cor: Logical. In step 3 of the algorithm, specifies whether to use column (drug) correlation or row (adverse event) correlation. The default is FALSE, indicating the use of adverse event correlation. TRUE indicates the use of drug correlation.

• cor_lim: A numeric value between 0 and 1. Specifies the correlation threshold to use in step 3 for selecting “connected” adverse events. The default is 0.8.

• num_cores: Number of cores used to parallelize the MDDC Boxplot algorithm. Default is 2.

• seed: An optional integer to set the seed for reproducibility. If NULL, no seed is set.

We now apply MDDC (MC) algorithm to statin49 using the following code:

set.seed(42)
test2 <- mddc_mc(
  contin_table = statin49,
  quantile = 0.95,
  rep = 10000,
  exclude_same_drug_class = T,
  col_specific_cutoff = T,
  separate = T,
  if_col_cor = F,
  cor_lim = 0.8
)

This function outputs a list with five components:

• mc_pval: Returns the p-values for each cell in step 2. For cells with counts greater than five, the p-values are obtained via the Monte Carlo (MC) method. For cells with counts less than or equal to five, the p-values are obtained via Fisher’s exact tests.

• mc_signal: Indicates signals for cells with counts greater than five, identified in step 2 by the MC method. A value of 1 indicates a signal, while 0 indicates no signal.

• fisher_signal: Indicates signals for cells with counts less than or equal to five, identified in step 2 by Fisher’s exact tests. A value of 1 indicates a signal, while 0 indicates no signal.

• corr_signal_pval: Returns the p-values for each cell in the contingency table in step 5, where the r_ij values are mapped back to the standard normal distribution.

• corr_signal_adj_pval: Returns the Benjamini-Hochberg adjusted p-values for each cell in step 5. Users can choose whether to use corr_signal_pval or corr_signal_adj_pval, and select an appropriate p-value threshold (for example, 0.05).