openEBGM Objects and Class Functions

Creating the Object

As mentioned in the other vignettes, the openEBGM package is capable of calculating \(EBGM\) and quantile scores from the posterior distribution. openEBGM makes it easy to calculate such quantities using a class and object system. While creation of objects of class openEBGM is not necessary (see previous vignette), it provides access to methods for some common generic functions and reduces the number of function calls needed.

To create the object, we first need to calculate the hyperparameter estimates.

library(openEBGM)
data(caers)
proc <- processRaw(caers, stratify = FALSE, zeroes = FALSE)
squashed <- squashData(proc)
theta_init <- data.frame(
  alpha1 = c(1, 2, 3),
  beta1 = c(1, 2, 3),
  alpha2 = c(2, 4, 5),
  beta2 = c(2, 4, 5),
  p = c(.1, 0.2, 0.3)
)
theta_init
#>   alpha1 beta1 alpha2 beta2   p
#> 1      1     1      2     2 0.1
#> 2      2     2      4     4 0.2
#> 3      3     3      5     5 0.3
hyper_estimate <- autoHyper(squashed, theta_init = theta_init, 
  zeroes = FALSE, squashed = TRUE, N_star = 1
)
#> Trying method 'nlminb'...
#>   Working on initial guess number 1 of 3...
#>   Working on initial guess number 2 of 3...
#>   Working on initial guess number 3 of 3...

Once we have the hyperparameter estimates and the processed data, we can calculate the \(EBGM\) scores and any desired quantile(s) from the posterior distribution.

ebout <- ebScores(proc, hyper_estimate = hyper_estimate,
  quantiles = c(5, 95)  #For the 5th and 95th percentiles
)
ebout_noquant <- ebScores(proc, hyper_estimate = hyper_estimate,
  quantiles = NULL  #For no quantiles
)

As seen above, we can calculate the \(EBGM\) scores with or without adding quantiles. If using quantiles, we can specify any number of quantiles.

my_theta <- list()  #autoHyper() returns a list
my_theta$estimates <- c(
  alpha1 = 3.77, beta1 = 0.51, alpha2 = 3.74, beta2 = 3.65, P = 0.048
)
ebout_custom_theta <- ebScores(proc, hyper_estimate = my_theta)

Using the Generic Functions

Once the object has been created, we can use class-specific methods for some of R’s generic functions (namely, print(), summary(), and plot()).

#We can print an openEBGM object to get a quick look at the contents
print(ebout)
#> 
#> There were 23 var1-var2 pairs with a QUANT_05 greater than 2
#> 
#> Top 5 Highest QUANT_05 Scores
#>                                                     var1
#> 3628  EU NATURAL CONCEPTION FOR HER FERTILITY AID AND MU
#> 9503              PRESERVISION AREDS 2 FORMULA SOFT GELS
#> 5253                                              KRATOM
#> 10387 RITUAL SYNBIOTIC PLUS PREBIOTIC PROBIOTIC POSTBIOT
#> 4030                                           FLORASTOR
#>                                               var2  N          E QUANT_05
#> 3628                          ABORTION SPONTANEOUS  9 0.16858790     7.96
#> 9503  NEOVASCULAR AGE-RELATED MACULAR DEGENERATION 13 0.59510086     7.65
#> 5253                                    DEPENDENCE 19 1.42615274     6.84
#> 10387                    GASTROINTESTINAL DISORDER  8 0.17291066     6.80
#> 4030                                     FUNGAEMIA  6 0.02269452     5.81
print(ebout_noquant, threshold = 3)
#> 
#> There were 51 var1-var2 pairs with an EBGM score greater than 3
#> 
#> Top 5 Highest EBGM Scores
#>                                                     var1
#> 3628  EU NATURAL CONCEPTION FOR HER FERTILITY AID AND MU
#> 10387 RITUAL SYNBIOTIC PLUS PREBIOTIC PROBIOTIC POSTBIOT
#> 9503              PRESERVISION AREDS 2 FORMULA SOFT GELS
#> 4030                                           FLORASTOR
#> 5253                                              KRATOM
#>                                               var2  N          E  EBGM
#> 3628                          ABORTION SPONTANEOUS  9 0.16858790 13.96
#> 10387                    GASTROINTESTINAL DISORDER  8 0.17291066 12.38
#> 9503  NEOVASCULAR AGE-RELATED MACULAR DEGENERATION 13 0.59510086 12.20
#> 4030                                     FUNGAEMIA  6 0.02269452 11.96
#> 5253                                    DEPENDENCE 19 1.42615274 10.07

When quantiles are present, simply printing the object shows, by default, how many var1-var2 pairs exist that have QUANT\(>x\), where \(x\) is the minimum quantile threshold used for the data (default 2). In the absence of quantiles, it simply outputs the number of var1-var2 pairs that have an \(EBGM\) score greater than the specified threshold. In both cases, it also shows a quick look at the var1-var2 pairs with the highest \(x\) or \(EBGM\), depending on whether quantiles were calculated or not.

One can also use the summary() function on an openEBM object to get further information about the calculations.

summary(ebout)

#> 
#> Summary of the EB-Metrics
#>       EBGM           QUANT_05       QUANT_95     
#>  Min.   : 0.290   Min.   :0.11   Min.   : 0.600  
#>  1st Qu.: 1.220   1st Qu.:0.49   1st Qu.: 2.920  
#>  Median : 1.250   Median :0.50   Median : 3.160  
#>  Mean   : 1.263   Mean   :0.51   Mean   : 3.153  
#>  3rd Qu.: 1.260   3rd Qu.:0.50   3rd Qu.: 3.270  
#>  Max.   :13.960   Max.   :7.96   Max.   :23.070

As seen above, by default the summary() function, when called on an openEBGM object, outputs some descriptive statistics on the \(EBGM\) and quantile scores, and a histogram of the \(EBGM\) scores. There are options to disable plot output, or to calculate the log₂ transform of the scores, which provides a Bayesian information statistic (when applied to the \(EBGM\) score).

summary(ebout, plot.out = FALSE, log.trans = TRUE)
#> 
#> Summary of the EB-Metrics
#>       EBGM            QUANT_05          QUANT_95     
#>  Min.   :-1.7859   Min.   :-3.1844   Min.   :-0.737  
#>  1st Qu.: 0.2869   1st Qu.:-1.0291   1st Qu.: 1.546  
#>  Median : 0.3219   Median :-1.0000   Median : 1.660  
#>  Mean   : 0.3179   Mean   :-0.9953   Mean   : 1.626  
#>  3rd Qu.: 0.3334   3rd Qu.:-1.0000   3rd Qu.: 1.709  
#>  Max.   : 3.8032   Max.   : 2.9928   Max.   : 4.528

Finally, openEBGM provides a method for the plot() function that can produce a variety of different plots. These are shown below.

plot(ebout)

As seen, by default, the plot() function shows the top \(EBGM\) scores by var1-var2 combinations (only var1 is shown for space preservation) and “error bars” using the lowest and highest quantiles calculated. The sample size for each var1-var2 combination is also plotted.

A specific event from var2 may also be selected, and only the var1-var2 combinations that include this particular event will be shown. An example is shown below.

plot(ebout, event = "CHOKING")
#> Warning in plot.openEBGM(ebout, event = "CHOKING"): 2 or more matches found for
#> event specified

In addition to the bar chart, the plot() function can also create a histogram of the \(EBGM\) scores.

plot(ebout, plot.type = "histogram")
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Again, one may choose an event from var2 by which to subset the data when plotting.

plot(ebout, plot.type = "histogram", event = "CHOKING")
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Finally, the last type of plot included with the plot() function shows the shrinkage performed by the algorithm. It is called the “Chirtel Squid Plot”, titled after its creator, Stuart Chirtel.

plot(ebout, plot.type = "shrinkage")

While a specific event may be selected by which to subset the data, it can lead to a less informative plot due to smaller sample size.

plot(ebout, plot.type = "shrinkage", event = "CHOKING")

Conclusion

openEBGM was designed to give the user a high level of control over data analysis choices (stratification, data squashing, etc.) using DuMouchel’s (1999, 2001) Gamma-Poisson Shrinkage (GPS) method. The GPS method applies to any large contingency table, so openEBGM can be used to mine a variety of databases in which the rate of co-occurrence of two variables or items is of interest (sometimes known as the “market basket problem”). U.S. FDA products and adverse events is just one of many possible applications.