Step 0: Loading Libraries and Defining 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: Writing a Function to Generate Counts

get_css_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)
  
  return(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: Writing a Function to Compute CSS

get_css_result <- function(n111, n110, n101, n100, n011, n010, n001, n000){
  
  # Compute marginal counts
  n1r  <- n111 + n110 + n101 + n100  # Drug 1 present
  n0r  <- n011 + n010 + n001 + n000  # Drug 1 absent
  n.1r <- n111 + n011 + n101 + n001  # AE present
  n.0r <- n110 + n010 + n100 + n000  # AE absent
  
  
  # Recode for CSS
  n11r <- n111  # Drug 1 & Drug 2 & AE
  n00r <- n000 + n010 + n100  # Neither Drug 1 nor Drug 2
  n10r <- n110  # Drug 1 & Drug 2, No AE
  n01r <- n001 + n011 + n101  # AE, but no Drug 1 or Drug 2
  
  # PRR Calculations
  PRR_D1D2 <- (n11r / n1r) / (n01r / n0r)  # PRR for Drug 1 + Drug 2
  PRR_D1   <- ((n11r + n10r) / n1r) / ((n01r + n00r) / n0r)  # PRR for Drug 1
  PRR_D2   <- ((n11r + n01r) / n.1r) / ((n10r + n00r) / n.0r)  # PRR for Drug 2
  
  # Confidence Intervals for PRRs
  log_Prr_D1D2 <- log(PRR_D1D2)
  log_Prr_D1   <- log(PRR_D1)
  log_Prr_D2   <- log(PRR_D2)
  
  # Compute standard error for PRRs
  se_log_Prr_D1D2 <- sqrt((1 / n11r) - (1 / n1r) + (1 / n01r) - (1 / n0r))
  se_log_Prr_D1 <- sqrt((1 / (n11r + n10r)) - (1 / n1r) + (1 / (n01r + n00r)) - (1 / n0r))
  se_log_Prr_D2   <- sqrt((1 / (n11r + n01r)) - (1 / n.1r) + 
                            (1 / (n10r + n00r)) - (1 / n.0r))
  
  
  
  # Compute lower (PRR025) and upper (PRR975) confidence bounds
  PRR025_D1D2 <- exp(log_Prr_D1D2 - 1.96 * se_log_Prr_D1D2)
  PRR975_D1D2 <- exp(log_Prr_D1D2 + 1.96 * se_log_Prr_D1D2)

  PRR025_D1   <- exp(log_Prr_D1 - 1.96 * se_log_Prr_D1)
  PRR975_D1   <- exp(log_Prr_D1 + 1.96 * se_log_Prr_D1)
  
  PRR025_D2   <- exp(log_Prr_D2 - 1.96 * se_log_Prr_D2)
  PRR975_D2   <- exp(log_Prr_D2 + 1.96 * se_log_Prr_D2)
  
  # Compute CSS Score
  CSS <- PRR025_D1D2 / max(PRR975_D1, PRR975_D2)
  
   # Store results in a list
  css_result <- data.table(
    PRR_D1D2 = PRR_D1D2,
    PRR025_D1D2 = PRR025_D1D2,
    PRR975_D1D2 = PRR975_D1D2,
    
    PRR_D1 = PRR_D1,
    PRR025_D1 = PRR025_D1, 
    PRR975_D1 = PRR975_D1,
    
    PRR_D2 = PRR_D2,
    PRR025_D2 = PRR025_D2,
    PRR975_D2 = PRR975_D2,
    
    CSS = CSS
  )
  
return(css_result)
}

Step 3: Calling the Functions

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


css_counts <- as.data.table(get_css_table(d1 = d1_study, d2 = d2_study, ae = ae_study, merged = data))
list2env(css_counts, envir = .GlobalEnv)
## <environment: R_GlobalEnv>
css_result <- get_css_result(n111 = n111, 
                             n110 = n110, 
                             n101 = n101, 
                             n100 = n100, 
                             n011 = n011,
                             n010 = n010,
                             n001 = n001,
                             n000 = n000)

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

# Reorder columns: `d1`, `d2`, `ae` first
setcolorder(css_result, c("d1", "d2", "ae"))

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

# Reorder columns: `d1`, `d2`, `ae` first
setcolorder(css_counts, c("d1", "d2", "ae"))

Signal Detection

  1. PRR025 drug D1 and drug D2 > 1

  2. CSS > 1

DT::datatable(css_counts, options = list(pageLength = 10, scrollX = TRUE))
DT::datatable(css_result, options = list(pageLength = 10, scrollX = TRUE))