Step 0: Loading Libraries and Setting Parameters

library(data.table)

d1_study <- ("RITONAVIR") #Drug 1 of DDI
d2_study <- ("ATORVASTATIN") #Drug 2 of DDI
ae_study <- ("MYOPATHY/RHABDOMYOLYSIS") #AE of DDI

file_path <- ("C://Users//William//Desktop//DDI Codes//data//processed//1.csv")

Step 1: Function to Generate Contingency Table

get_omega_table <- function(d1, d2, ae, merged){
  
  merged$d1i <- ifelse(merged$active_ingredient == d1, 1, 0)
  merged$d2i <- ifelse(merged$active_ingredient == d2, 1, 0)
  merged$aei <- ifelse(merged$pt == ae, 1, 0)
  
  merged <- merged[, .(id, d1i, d2i, aei)]
  
  flat <-  merged[, .(d1i = max(d1i), d2i = max(d2i), aei = max(aei)), by = id]
  
  #individual cell counts
  n111 <- nrow(flat[flat$d1i == 1 & flat$d2i == 1 & flat$aei == 1])
  n110 <- nrow(flat[flat$d1i == 1 & flat$d2i == 1 & flat$aei == 0])
  
  n101 <- nrow(flat[flat$d1i == 1 & flat$d2i == 0 & flat$aei == 1])
  n100 <- nrow(flat[flat$d1i == 1 & flat$d2i == 0 & flat$aei == 0])
  
  n011 <- nrow(flat[flat$d1i == 0 & flat$d2i == 1 & flat$aei == 1])
  n010 <- nrow(flat[flat$d1i == 0 & flat$d2i == 1 & flat$aei == 0])
  
  n001 <- nrow(flat[flat$d1i == 0 & flat$d2i == 0 & flat$aei == 1])
  n000 <- nrow(flat[flat$d1i == 0 & flat$d2i == 0 & flat$aei == 0])
  
  #marginal cell counts
  n..1 <- nrow(flat[flat$aei == 1])
  n..0 <- nrow(flat[flat$aei == 0])
  
  n11. <- nrow(flat[flat$d1i == 1 & flat$d2i == 1])
  n10. <- nrow(flat[flat$d1i == 1 & flat$d2i == 0])
  n01. <- nrow(flat[flat$d1i == 0 & flat$d2i == 1])
  n00. <- nrow(flat[flat$d1i == 0 & flat$d2i == 0])
  
  #total cell count
  n... <- nrow(flat)
  
  #construct the table
  contingency <- rbind(
    c(n111, n110, n11.),  # Drug D1 and Drug D2
    c(n101, n100, n10.),  # Only Drug D1
    c(n011, n010, n01.),  # Only Drug D2
    c(n001, n000, n00.),  # Neither Drug D1 nor Drug D2
    c(n..1, n..0, n...)   # Total
  )
  rownames(contingency) <- c("Drug D1 and Drug D2", "Only Drug D1", "Only Drug D2", "Neither Drug D1 nor Drug D2", "Total")
  colnames(contingency) <- c("Target AE", "All Other AEs", "Total")
  
  return(list(
    contingency_table = contingency,
    cell_counts = list(
      n111 = n111, n110 = n110, n101 = n101, n100 = n100, 
      n011 = n011, n010 = n010, n001 = n001, n000 = n000,
      n..1 = n..1, n..0 = n..0, n11. = n11., n10. = n10., 
      n01. = n01., n00. = n00., n... = n...
    )
  ))
}

Step 2: Function to Calculate Omega Shrinkage Estimate

get_omega_result <- function(n001,n00.,n101,n10.,n011,n01.,n111,n11.){
  
#Compute observed relative reporting rate

f00 <- n001/n00. #Rate of adverse event in reports without usage of both drug
f10 <- n101/n10. #Rate of adverse event in reports with usage of drug 1 but not 2
f01 <- n011/n01. #Rate of adverse event in reports with usage of drug 2 but not 1
f11 <- n111/n11. #Rate of adverse event in reports with usage of both drug 1 and 2

#Compute estimator of expected relative reporting rate, g11, for rate of adverse event in reports with usage of both drug 1 and 2

max1 <- max(f00 / (1 - f00), f10 / (1 - f10))
max2 <- max(f00 / (1 - f00), f01 / (1 - f01))

g11 <- 1 - 1 / (max1 + max2 - f00 / (1 - f00) + 1) #Estimatated value of f11 E[f11]

#When f10<f00 (indicating no risk of AE attributable to D1), this yields the most sensible estimator g11=max(f00, f01) and vice versa when f01< f00.

E111 <- g11*n11. #Expected value of drug1/drug2/ae, Equivalent to E[11]*n[11.]

omega <- log2((n111+0.5)/(E111+0.5)) #alpha = 0.5 is used as default, it provides enough shirinkage. 

#Now, there are both frequentists and bayesian approaches to generate the uncertainty:

#Frequentist 95CI
omega_se <- 1.96/(log(2)*(sqrt(n111))) #Equation 18 on Noguchi et al. 

omega_lb_fr <- omega - omega_se 
omega_ub_fr <- omega + omega_se


#Bayesian 95CRI

alpha <- n111 + 0.5 
beta <- E111 + 0.5

#A Gamma G(a,a) prior for the Poisson rate parameter, a, specifically  a= 0.5 is used in the measure omega 

omega_lb_by <- log2(qgamma(0.025, shape = alpha, rate = beta))
omega_ub_by <- log2(qgamma(0.975, shape = alpha, rate = beta))

result <- data.table(
  omega = omega,
  omega_lb_fr = omega_lb_fr,
  omega_ub_fr = omega_ub_fr,
  omega_lb_by = omega_lb_by,
  omega_ub_by = omega_ub_by
)

return(result)
}

Step 3: Test

Reading in the data

data <- fread(file_path)
setnames(data, old = names(data)[1], new = "id")
data$pt <- toupper(data$pt)

Calling the contingency table function

omega_table <- get_omega_table(
  d1 = d1_study,
  d2 = d2_study, 
  ae = ae_study, 
  merged = data)
list2env(omega_table$cell_counts, envir = .GlobalEnv)
## <environment: R_GlobalEnv>

Calling the function to compute estimator

counts <- data.table(   d1 = d1_study,
                        d2 = d2_study,
                        ae = ae_study,
                        n111 = n111,
                        n110 = n110,
                        n11. = n11.,
                        n011 = n011,
                        n010 = n010,
                        n01. = n10.,
                        n001 = n001,
                        n000 = n000,
                        n00. = n00.,
                        n... = n...)

result <- get_omega_result(n001, n00., n101, n10., n011, n01., n111, n11.)

# Add new columns for labeling the signal
result[, `:=` (d1 = d1_study, d2 = d2_study, ae = ae_study)]

# Reorder columns: `d1`, `d2`, `ae` first
setcolorder(result, c("d1", "d2", "ae"))
DT::datatable(counts, options = list(pageLength = 10, scrollX = TRUE))
DT::datatable(result, options = list(pageLength = 10, scrollX = TRUE))

Signal Determination

DDI signals are detected if lower bound > 0