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 \times 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\times 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\times 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\times 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).
Function for Finding Optimal Boxplot Coefficient
We provide a function that outputs appropriate coefficients for
boxplot using a grid search method to control the FDR. This function
takes the following arguments:
contin_table: A matrix representing the \(I \times J\) contingency table.
n_sim: An integer specifying the number of simulated
tables under the assumption of independence between rows and columns.
Default is 1000.
target_fdr: A numeric value specifying the desired
level of false discovery rate (FDR). Default is 0.05.
grid: A numeric value representing the size of the
grid added to the default value of coef = 1.5 as suggested
by Tukey. Default is 0.1.
col_specific_cutoff: Logical. If TRUE, then a single
value of the coefficient is returned for the entire dataset, else when
FALSE specific values correspoding to each of the columns are
returned.
exclude_small_count: A logical indicating whether to
exclude cells with counts smaller than or equal to five when computing
boxplot statistics. Default is TRUE.
We apply this function to the statin49 dataset:
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
This function outputs a list with the following components:
coef: A numeric vector containing the optimal
coefficient ‘coef’ for each column of the input contingency
table.
FDR: A numeric vector with the corresponding false
discovery rate (FDR) for each column.
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\times 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\times2\) 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).
