<jupyter_start><jupyter_text>Proportionality Assessment using R and the openEBGM packageTo analyse the reported proportionality ratio, we will structure the data as needed by the openEBGM  [[2]](2) and determine the PRR. Additionally, we will use the Gamma-Poission Shrinkage model and use Empricial Bayes Geometric Mean (EBGM) to report the ratio. We will follow the example as presented here to achieve this :(https://cran.r-project.org/web/packages/openEBGM/vignettes/x1_introAndDataPrepVignette.html)<jupyter_code>install.packages('openEBGM')
options(warn=-1)
library(openEBGM)
library(repr)

# Change plot size to 8 x 3
options(repr.plot.width=10, repr.plot.height=10,repr.plot.res = 100)
df<-read.csv('./data/product_event_strat.csv')
head(df)
df[df=='0']<- NA #sex == 0 is unknown, replace with na
df$strat_gender<- df$sex
table(df$strat_gender, useNA="always")
df$strat_age <- ifelse(is.na(df$age), "unknown",
                 ifelse(df$age < 18, "under_18",
                        "18_plus"))
table(df$strat_age, useNA = "always")
vars<-c("id", "var1","var2", "strat_gender", "strat_age") #select the columns
#remove any na
df=na.omit(df)
df=df[vars]
head(df)<jupyter_output><empty_output><jupyter_text>We have processed the data into the right format for using the openEBGM package. We can calculate the RR and PRR for the product-event pairs.- $N$ indicates the actual count for the cell. - $E$ is the expected count under the assumption of independence between rows and columns.  When RR is 1,we observe the exact count you would expect to observe if no association exists between the two variables. When RR>1, we observe a larger count than expected.<jupyter_code>processRaw(data=df, stratify = FALSE,zeroes = FALSE)<jupyter_output><empty_output><jupyter_text>Take for example the first row in the above table. Small counts cause the RR to be unstable. That is when the a single observation for a low expected counted (E=0.15), the RR becomes large. For this reason, we will explore Empirical Bayes approach.   Empirical Bayes ApproachIn order to handle the issues of small counts, we will explore using the Gamma-Poisson Shrinker(GPS) model[[1]](1). The emprical Bayes model shrinks the large RR  from small counts to a value much closer to "null hypothesis" value of 1. This shrinkage is larger for smaller counts and negligible for larger counts. The GPS model uses a Poisson ($\mu_{ij}$) likelihood for actual cell count $N_{ij}$ and the expected counts $E_{ij}$ are treated as constants. We are interested in the ratio $\lambda_{ij}=\frac{\mu_{ij}}{E_{ij}}$. The prior on $\lambda$ is a mixture of two gamma distributions and thus the posterior is also a mixture of two gamma distributions. The Emprical Bayes Geometric Mean(EBGM) of the posterior distribution is normally used in place of RR. Additionally, we can use the quantiles as well. We will explore them below. In the Empirical Bayes approach, the prior is driven by the data and this implies that the known posterior distribution relies on a set of hyper-parameters derived from the prior distribution.  Estimating HyperparametersThe actual counts (N) and expected counts (E) are used to estimate the hyperparameters of the prior distribution. A large contingency table will have many cells, resulting in computational difficulties for the optimization routines needed for estimation. Data squashing [[1]](1). transforms a set of 2-dimensional points to a smaller number of 3-dimensional points.For a given N, "squashing" combines points with similar Es into bins using a specified bin size and uses the average E within each bin as the E for that bin’s “superpoint”. The new superpoints are weighted by bin size. The hyperparameters are estimated by minimizing the negative log-likelihood function.  The hyperparameters are denoted by the vector θ=(α1,β1,α2,β2,P). P is the mixture fraction. Hyperparameters can be calculated by exploring the parameter space of the likelihood function using either the full data set of Ns and Es or the squashed set.<jupyter_code>processed=processRaw(data=df)
squashed <- squashData(processed) #Using defaults
head(squashed)<jupyter_output><empty_output><jupyter_text>**Likelihood Functions**<jupyter_code>theta_init <- c(alpha1 = 0.2, beta1 = 0.1, alpha2 = 2, beta2 = 4, p = 1/3)# These values are taken from the example in the R package documentation
stats::nlm(negLLsquash, p = theta_init,
           ni = squashed$N, ei = squashed$E, wi = squashed$weight, N_star = 1)
theta_init <- data.frame(alpha1 = c(0.2, 0.1, 0.3),
                         beta1  = c(0.1, 0.2, 0.5),
                         alpha2 = c(2,   10,  6),
                         beta2  = c(4,   10,  6),
                         p      = c(1/3, 0.2, 0.5)
)
squashed2 <- squashData(squashed, count = 2, bin_size = 10, keep_pts = 50)
system.time(
hyper_estimates_squashed <- autoHyper(data = squashed2, theta_init = theta_init)
)
hyper_estimates_squashed
exploreHypers(data = squashed2, theta_init = theta_init, std_errors = TRUE)
(theta_hat <- hyper_estimates_squashed$estimates)

#Qn() correspond to the probability that λ came from the first component of the posterior mixture
qn <- Qn(theta_hat, N = processed$N, E = processed$E)#calculates the mixture fractions for the posterior distributions.
head(qn)
identical(length(qn), nrow(processed))# Qn output should match the number of rows in the original data
#EBGM  is a measure of central tendency of the posterior distributions, λij|N=n
processed$ebgm <- ebgm(theta_hat, N = processed$N, E = processed$E, qn  = qn)
head(processed)
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)
suspicious <- processed[processed$QUANT_05 >= 2, ]
nrow(suspicious); nrow(processed); nrow(suspicious)/nrow(processed)<jupyter_output><empty_output><jupyter_text>From the above, less than 1% of product-event pairs are suspect.<jupyter_code>suspicious <- suspicious[order(suspicious$QUANT_05, decreasing = TRUE),
                         c("var1", "var2", "N", "E", "QUANT_05", "ebgm", 
                           "QUANT_95")]
head(suspicious, 5)<jupyter_output><empty_output><jupyter_text>From the above table, most suspicious (top 5) data all belongs "carmustine". Events such as venoocclusive liver disease, basal cell carcinoma etc are recorded much more than expected. These may need further investigation with subject matter experts and will depend on domain-knowledge.<jupyter_code>tabbed <- table(suspicious$var1)
head(tabbed[order(tabbed, decreasing = TRUE)])
proc=processRaw(data=df)
ebout=ebScores(proc, hyper_estimate = hyper_estimates_squashed, quantiles = c(5,95))
ebout
summary(ebout)
plot(ebout, event="Myelodysplastic syndrome",cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09301659988551803, 'fraction_numerical': 0.01760160274756726, 'mean_word_length': 4.46015625, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 27, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '519812', 'license': 'MIT', 'n_tokens_mistral': 2237, 'n_tokens_neox': 2039, 'n_words': 796, 'path': 'Proportional Analysis using R.ipynb', 'repo_name': 'titan-king/openFDA_ADE_Analysis', 'size': '6988'}