In Bayesian statistics, the gamma distribution is the conjugate prior distribution for a Poisson likelihood. ‘Conjugate’ means that the posterior distribution will follow the same general form as the prior distribution. DuMouchel (1999) used a model with a Poisson(\(\mu_{ij}\)) likelihood for the counts (for row i and column j of the contingency table). We are interested in the ratio \(\lambda_{ij}=\frac{\mu_{ij}}{E_{ij}}\), where \(E_{ij}\) are the expected counts. The \(\lambda_{ij}\)s are considered random draws from a mixture of two gamma distributions (our prior) with hyperparameter \(\theta=(\alpha_1,\beta_1,\alpha_2,\beta_2,P)\), where \(P\) is the prior probability that \(\lambda\) came from the first component of the prior mixture (i.e., the mixture fraction). The prior is a single distribution that models all the cells in the table; however, there is a separate posterior distribution for each cell in the table. The posterior distribution of \(\lambda\), given count \(N=n\), is a mixture of two gamma distributions with parameters \(\theta=(\alpha_1+n,\beta_1+E,\alpha_2+n,\beta_2+E,Q_n)\) (subscripts suppressed for clarity), where \(Q_n\) is the probability that \(\lambda\) came from the first component of the posterior, given \(N=n\) (i.e., the mixture fraction).
The posterior distribution, in a sense, is a Bayesian representation of the relative reporting ratio, \(RR\) (note the similarity in the equations \(RR_{ij}=\frac{N_{ij}}{E_{ij}}\) and \(\lambda_{ij}=\frac{\mu_{ij}}{E_{ij}}\)). The Empirical Bayes (EB) metrics are taken from the posterior distribution. The Empirical Bayes Geometric Mean \((EBGM)\) is the antilog of the mean of the log2-transformed posterior distribution. The \(EBGM\) is therefore a measure of central tendency of the posterior distribution. The 5% and 95% quantiles of the posterior distributions can be used to create two-sided 90% credibility intervals for \(\lambda_{ij}\), given \(N_{ij}\) (i.e, our “sort of” RR). Alternatively, since we are most interested in the lower bound, we could ignore the upper bound and create a one-sided 95% credibility interval.
Due to Bayesian shrinkage (please see the Background section of the Introduction to openEBGM vignette), the EB scores are much more stable than \(RR\) for small counts.
Once the product/event combinations have been counted and the hyperparameters have been estimated, we can calculate the EB scores:
library(openEBGM)
data(caers) #subset of publicly available HFCS (formerly CAERS) data
processed <- processRaw(caers, stratify = FALSE, zeroes = FALSE)
squashed <- squashData(processed)
squashed <- squashData(squashed, count = 2, bin_size = 10, keep_pts = 50)
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_estimates <- autoHyper(squashed, theta_init = theta_init)
#> 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...
(theta_hat <- hyper_estimates$estimates)
#> alpha1 beta1 alpha2 beta2 P
#> 3.77699506 0.51388202 3.73820989 3.64892019 0.04817177Qn()The Qn() function calculates the mixture fractions for
the posterior distributions. The values returned by Qn()
correspond to the probability that \(\lambda\) came from the first component of
the posterior mixture distribution, given \(N=n\) (recall there is a \(\lambda|N=n\) for each cell in the table,
but that each \(\lambda\) comes from a
common distribution). Thus, the output from Qn() returns a
numeric vector of length equal to the total number of product-symptom
combinations, which is also the number of rows in the data frame
returned by processRaw(). When calculating the \(Q_n\)s, be sure to use the full data set
from processRaw() – not the squashed data set or the raw
data.
ebgm()The ebgm() function calculates the Empirical Bayes
Geometric Mean \((EBGM)\) scores. \(EBGM\) is a measure of central tendency of
the posterior distributions, \(\lambda_{ij}|N=n\). Scores much larger than
one indicate product/adverse event pairs that are reported at an
unusually high rate.
processed$ebgm <- ebgm(theta_hat, N = processed$N, E = processed$E, qn = qn)
head(processed)
#> var1 var2 N
#> 1 7 ZEN 7-HYDROXYMITRAGYNINE CHEWABLE TABLETS ACUTE PSYCHOSIS 1
#> 2 7 ZEN 7-HYDROXYMITRAGYNINE CHEWABLE TABLETS HOSPITALISATION 1
#> 3 7-HYDROXYMITRAGYNINE 7-OH KRATOM DEATH 1
#> 4 ARTRI KING BLOOD PRESSURE ABNORMAL 1
#> 5 ARTRI KING CEREBROVASCULAR ACCIDENT 1
#> 6 ARTRI KING SWELLING 1
#> E RR PRR ebgm
#> 1 0.0003602305 2776.00 Inf 1.27
#> 2 0.1253602305 7.98 8.00 1.19
#> 3 0.0407060519 24.57 24.78 1.24
#> 4 0.0021613833 462.67 555.00 1.27
#> 5 0.0075648415 132.19 138.75 1.26
#> 6 0.0190922190 52.38 53.37 1.25quantBisect()The quantBisect() function calculates quantiles of the
posterior distribution using the bisection method.
quantBisect() can calculate any quantile of the posterior
distribution between 1 and 99%, and these quantiles can be used as
limits for credibility intervals. Below, QUANT_05 is the
5th percentile; QUANT_95 is the 95th
percentile. These form the lower and upper bounds of 90% credibility
intervals for the Empirical Bayes (EB) scores.
processed$QUANT_05 <- quantBisect(5, theta_hat = theta_hat,
N = processed$N, E = processed$E, qn = qn
)
processed$QUANT_95 <- quantBisect(95, theta_hat = theta_hat,
N = processed$N, E = processed$E, qn = qn
)
head(processed)
#> var1 var2 N
#> 1 7 ZEN 7-HYDROXYMITRAGYNINE CHEWABLE TABLETS ACUTE PSYCHOSIS 1
#> 2 7 ZEN 7-HYDROXYMITRAGYNINE CHEWABLE TABLETS HOSPITALISATION 1
#> 3 7-HYDROXYMITRAGYNINE 7-OH KRATOM DEATH 1
#> 4 ARTRI KING BLOOD PRESSURE ABNORMAL 1
#> 5 ARTRI KING CEREBROVASCULAR ACCIDENT 1
#> 6 ARTRI KING SWELLING 1
#> E RR PRR ebgm QUANT_05 QUANT_95
#> 1 0.0003602305 2776.00 Inf 1.27 0.50 3.31
#> 2 0.1253602305 7.98 8.00 1.19 0.48 2.75
#> 3 0.0407060519 24.57 24.78 1.24 0.50 3.05
#> 4 0.0021613833 462.67 555.00 1.27 0.50 3.29
#> 5 0.0075648415 132.19 138.75 1.26 0.50 3.26
#> 6 0.0190922190 52.38 53.37 1.25 0.50 3.18The EB-scores (\(EBGM\) and quantile scores) can be used to look for “signals” in the data. As stated in the Background section of the Introduction to openEBGM vignette, Bayesian shrinkage causes the EB-scores to be far more stable than their \(RR\) counterparts, which allows for better separation between signal and noise. One could, for example, look at all product-symptom combinations where QUANT_05 (the lower part of the 90% two-sided credibility interval) is 2 or greater. This is often used as a conservative alternative to \(EBGM\) since QUANT_05 scores are naturally smaller than \(EBGM\) scores. We can say with high confidence that the “true relative reporting ratios” of product/adverse event combinations above this threshold are much greater than 1, so those combinations are truly reported more than expected. The value of 2 is arbitrarily chosen, and depends on the context. Below is an example of how one may identify product-symptom combinations that require further investigation based on the EB-scores.
suspicious <- processed[processed$QUANT_05 >= 2, ]
nrow(suspicious); nrow(processed); nrow(suspicious)/nrow(processed)
#> [1] 24
#> [1] 13441
#> [1] 0.001785581From above we see that less than 1% of product-symptom pairs are suspect based on the QUANT_05 score. One may look more closely at these product-symptom combinations to ascertain which products may need further investigation. Subject matter knowledge is required to determine which signals might identify a possible causal relationship. The EB-scores find statistical associations – not necessarily causal relationships.
suspicious <- suspicious[order(suspicious$QUANT_05, decreasing = TRUE),
c("var1", "var2", "N", "E", "QUANT_05", "ebgm",
"QUANT_95")]
head(suspicious, 5)
#> var1
#> 3628 EU NATURAL CONCEPTION FOR HER FERTILITY AID AND MULTI
#> 9503 PRESERVISION AREDS 2 FORMULA SOFT GELS
#> 5253 KRATOM
#> 10387 RITUAL SYNBIOTIC PLUS PREBIOTIC PROBIOTIC POSTBIOTIC
#> 4030 FLORASTOR
#> var2 N E QUANT_05 ebgm
#> 3628 ABORTION SPONTANEOUS 9 0.16858790 7.96 13.95
#> 9503 NEOVASCULAR AGE-RELATED MACULAR DEGENERATION 13 0.59510086 7.65 12.20
#> 5253 DEPENDENCE 19 1.42615274 6.84 10.06
#> 10387 GASTROINTESTINAL DISORDER 8 0.17291066 6.80 12.38
#> 4030 FUNGAEMIA 6 0.02269452 5.81 11.96
#> QUANT_95
#> 3628 23.06
#> 9503 18.67
#> 5253 14.39
#> 10387 21.06
#> 4030 22.10
tabbed <- table(suspicious$var1)
head(tabbed[order(tabbed, decreasing = TRUE)])
#>
#> PRESERVISION AREDS 2 FORMULA SOFT GELS AG1 WHOLE FOOD POUCH 360G NON SPECIFIC
#> 7 2
#> CITRACAL MAXIMUM CENTRUM SILVER MEN 50 PLUS
#> 2 1
#> CENTRUM SILVER WOMEN 50 PLUS CITRACAL MAXIMUM PLUS
#> 1 1The output above suggests some products which may warrant further investigation.
Next, the openEBGM Objects and Class Functions vignette will demonstrate the object-oriented features of the openEBGM package.