Overview

This markdown documents results of two different approaches to computing the extended-EBGM signal scores of DDI.

  1. The first approach mirrors the analysis performed in Almenoff 2003. DOI: 10.1002/pds.885, where only the DDIs of interest are re-coded as combinations (ie. “Drug A + Drug B”)

  2. The second is our approach, where we add in records of all drug combinations against each adverse event.

In this markdown, the openEBGM package is used. Both CAERS data and a single quarter of FAERS data is used for validation. CAERS markdown can be found

Step 0: Loading Libraries and Defining Parameters

library(openEBGM)
library(data.table)
library(dplyr)
library(dtplyr) #dplyr for data.table
library(tidyr)
library(purrr)
library(knitr)
library(DT)


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: Data Preperation

Almenoff method: re-code only DDI of interest

recode_1 <- function(data, d1, d2, ae){
  
  data[,role_cod:=NULL] #do not need this column
  
  setnames(data, old = names(data)[1], new = "id") #need to rename these for openEBGM
  setnames(data, old = names(data)[2], new = "var1")
  setnames(data, old = names(data)[3], new = "var2")
    
  common_ids <- intersect(data[var1 %in% d1_study]$id, data[var1 %in% d2_study]$id)

  data[id%in%common_ids & var1 %in% d1_study]$var1 <- paste0(d1_study, " + ", d2_study)

  data <- data[!(id %in% common_ids & var1 %in% d2_study)]

  return(data)
}

Our method: re-code all DDIs

recode_2 <- function(data) {
  # Remove unnecessary column
  data[, role_cod := NULL]

  # Rename columns for openEBGM
  setnames(data, old = names(data)[1:3], new = c("id", "var1", "var2"))

  # Function to generate Cartesian pairs while keeping single records
  expand_pairs <- function(sub_data) {
    id_var2 <- unique(sub_data$var2)  # Get unique var2 per id
    var1_unique <- unique(sub_data$var1)  # Get unique var1 values

    # Generate all possible unique pairs of var1
    var1_combinations <- CJ(var1_a = var1_unique, var1_b = var1_unique, unique = TRUE)[var1_a < var1_b]
    var1_combinations[, var1_combined := paste(var1_a, var1_b, sep = " + ")]
    var1_combinations[, c("var1_a", "var1_b") := NULL]  # Remove unnecessary columns

    # Keep individual records as well
    single_records <- data.table(var1_combined = var1_unique)

    # Combine individual and pair records
    all_records <- unique(rbind(single_records, var1_combinations))

    # Expand with var2 values
    expanded_data <- CJ(var1 = all_records$var1_combined, var2 = id_var2, unique = TRUE)

    return(expanded_data)
  }

  # Apply function by id
  expanded_data <- data[, expand_pairs(.SD), by = id]

  return(expanded_data)
}

Step 2: Writing a function for EBGM procedure

Writing the function

get_ebgm <- function(data){
  
  cat("Calculating Counts, Expected Value, and RR","\n")
  
  processed <- as.data.table(processRaw(data = data, stratify = FALSE, zeroes = FALSE)) 
  
  cat("Squashing","\n")

  squashed <- autoSquash(processed)

  theta_init <- c(alpha1 = 0.2, beta1 = 0.1, alpha2 = 2, beta2 = 4, p = 1/3)
  
  cat("Maximizing","\n")

  nlminb <- stats::nlminb(start = theta_init, objective = negLLsquash,
                            ni = squashed$N, ei = squashed$E, wi = squashed$weight, 
                            N_star = 1)

  theta_hat <- nlminb$par
  
  print(theta_hat)
  
  cat("Calculating EBGM","\n")

  qn <- Qn(theta_hat, N=processed$N, E=processed$E)

  processed$ebgm <- ebgm(theta_hat, N=processed$N, E=processed$E, qn = qn)
  
  processed$EB05 <- quantBisect(5, theta_hat = theta_hat, N=processed$N, E=processed$E, 
                                qn   = qn)
  processed$EB95 <- quantBisect(95, theta_hat = theta_hat, N=processed$N, E=processed$E, 
                                qn = qn)

return(processed)
}

Step 3: Call the Function with Prepared Data

The first approach gives the following:

data <- fread(file_path)

cat("Recoding","\n")
## Recoding
data_1 <- recode_1(data = data, d1 = d1_study, d2 = d2_study, ae = ae_study)

cat("Calculating EBGM","\n")
## Calculating EBGM
ic_result_1 <- get_ebgm(data_1)
## Calculating Counts, Expected Value, and RR 
## Squashing 
## Maximizing 
##       alpha1        beta1       alpha2        beta2            p 
## 1.307008e-08 1.650725e-02 3.974360e+00 2.843270e+00 1.229147e-01 
## Calculating EBGM
DT::datatable(
  ic_result_1[var1== paste0(d1_study," + ",d2_study)][order(var2)],  # Sort by var2
  options = list(pageLength = 10, scrollX = TRUE)
)
## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html
subset <- as.data.table(ic_result_1)
subset <- subset[N == 20][1:10]

DT::datatable(
  subset,
  options = list(pageLength = 10, scrollX = TRUE)
)
## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

The second approach (not running this block cause it takes too long)

#cat("Recoding","\n")

#data_2 <- recode_2(data = data)

#cat("Calculating EBGM","\n")

#ic_result_2 <- get_ebgm(data_2)

Signal Determination:

In Almenoff 2003, Interaction Signal Score is considered positive if:

  1. EB05 (AE, DDI) > 1, and

  2. EB05 (AE, DDI) > EB95 (AE, D1) , and

  3. EB05 (AE, DDI) ) > EB95 (AE, D2)