EXTENDED BCPNN
William Bai
2025-02-28
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: Function to Generate Contingency Table
get_ic_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])
#total cell count
n... <- nrow(flat)
#we need to recode this as the BCPNN paper codes the pairs differently, drug-ae-drug instead of drug-drug-ae, basically switching the 2nd index with third
# Recode the counts based on BCPNN paper convention
n111r <- n111 # This remains the same
n110r <- n101 # Drug-AE-No Drug
n101r <- n110 # Drug-No AE-Drug
n011r <- n011 # No Drug-AE-Drug (same)
n100r <- n100 # Drug-No AE-No Drug (same)
n000r <- n000 # No Drug-No AE-No Drug (same)
n010r <- n001 # No Drug-AE-No Drug
n001r <- n010 # No Drug-No AE-Drug
# Marginal totals using recoded values
n1..r <- n111r + n110r + n101r + n100r # Drug 1 = Yes (regardless of Drug 2 and AE)
n0..r <- n010r + n011r + n000r + n001r # Drug 1 = No
n.1.r <- n111r + n011r + n101r + n001r # Adverse Event = Yes (regardless of drugs)
n.0.r <- n110r + n010r + n100r + n000r # Adverse Event = No
n..1r <- n111r + n110r + n011r + n001r # Drug 2 = Yes (regardless of Drug 1 and AE)
n..0r <- n101r + n100r + n010r + n000r # Drug 2 = No
# Pairwise marginals using recoded values
n11.r <- n111r + n110r # Drug 1 = Yes, Drug 2 = Yes
n10.r <- n101r + n100r # Drug 1 = Yes, Drug 2 = No
n01.r <- n011r + n001r # Drug 1 = No, Drug 2 = Yes
n00.r <- n010r + n000r # Drug 1 = No, Drug 2 = No
n1.1r <- n111r + n101r # Drug 1 = Yes, Adverse Event = Yes
n1.0r <- n110r + n100r # Drug 1 = Yes, Adverse Event = No
n0.1r <- n011r + n001r # Drug 1 = No, Adverse Event = Yes
n0.0r <- n010r + n000r # Drug 1 = No, Adverse Event = No
n.11r <- n111r + n011r # Drug 2 = Yes, Adverse Event = Yes
n.10r <- n110r + n010r # Drug 2 = Yes, Adverse Event = No
n.01r <- n101r + n001r # Drug 2 = No, Adverse Event = Yes
n.00r <- n100r + n000r # Drug 2 = No, Adverse Event = No
# Storing all counts in a list with recoded values
counts <- list(
# Individual cell counts
n111 = n111r, n110 = n110r, n101 = n101r, n100 = n100r,
n011 = n011r, n010 = n010r, n001 = n001r, n000 = n000r,
# Marginal totals
n1.. = n1..r, n0.. = n0..r,
n.1. = n.1.r, n.0. = n.0.r,
n..1 = n..1r, n..0 = n..0r,
# Pairwise marginals
n11. = n11.r, n10. = n10.r, n01. = n01.r, n00. = n00.r,
n1.1 = n1.1r, n1.0 = n1.0r, n0.1 = n0.1r, n0.0 = n0.0r,
n.11 = n.11r, n.10 = n.10r, n.01 = n.01r, n.00 = n.00r,
# Total count
n... = n111r + n110r + n101r + n011r + n100r + n010r + n001r + n000r
)
return(counts)
}The cell counts correspond to the following table
Contingency Table Format for Extended
BCPNN
Step 2: Function to Compute Extended IC
get_ic_result <- function(n111, n110, n101, n011, n100, n010, n001, n000, n...) {
#Prior calculation:
#Smoothed marginal Probabilities
q1.. <- (n111 + n110 + n101 + n100 + 0.5) / (n... + 1) #Expected value for drug 1
q0.. <- (n011 + n010 + n001 + n000 + 0.5) / (n... + 1) #Expected value for no drug 1
q.1. <- (n111 + n011 + n101 + n001 + 0.5) / (n... + 1) #Expected value for ae
q.0. <- (n110 + n010 + n100 + n000 + 0.5) / (n... + 1) #Expected value for no ae
q..1 <- (n111 + n110 + n011 + n001 + 0.5) / (n... + 1) #Expected value for drug 2
q..0 <- (n101 + n100 + n010 + n000 + 0.5) / (n... + 1) #Expected value for no drug 2
#Smoothed pairwise Probabilities
q11. <- (n111 + n110 + 0.25) / (n... + 1) #Expected drug 1 and ae
q10. <- (n101 + n100 + 0.25) / (n... + 1) #Expected drug 1 and no ae
q01. <- (n011 + n001 + 0.25) / (n... + 1) #Expected no drug 1 and ae
q00. <- (n010 + n000 + 0.25) / (n... + 1) #Expected no drug 1 and no ae
q1.1 <- (n111 + n101 + 0.25) / (n... + 1) #Expected drug 1 and drug 2
q1.0 <- (n110 + n100 + 0.25) / (n... + 1) #Expected drug 1 and no drug 2
q0.1 <- (n011 + n001 + 0.25) / (n... + 1) #Expected no drug 1 and drug 2
q0.0 <- (n010 + n000 + 0.25) / (n... + 1) #Expected no drug 1 and no drug 2
q.11 <- (n111 + n011 + 0.25) / (n... + 1) #Expected drug 2 and ae
q.10 <- (n110 + n010 + 0.25) / (n... + 1) #Expected drug 2 and no ae
q.01 <- (n101 + n001 + 0.25) / (n... + 1) #Expected no drug 2 and ae
q.00 <- (n100 + n000 + 0.25) / (n... + 1) #Expected no drug 2 and no ae
# Step 2: Compute hyperparameters for moderating prior for third order IC
a... <- 0.5 * ((q1.. * q.1. * q..1) / (q11. * q1.1 * q.11))
a111 <- ((q11. * q1.1 * q.11) / (q1.. * q.1. * q..1)) * a...
a110 <- ((q11. * q1.0 * q.10) / (q1.. * q.1. * q..0)) * a...
a101 <- ((q10. * q1.1 * q.01) / (q1.. * q.0. * q..1)) * a...
a100 <- ((q10. * q1.0 * q.00) / (q1.. * q.0. * q..0)) * a...
a011 <- ((q01. * q0.1 * q.11) / (q0.. * q.1. * q..0)) * a...
a010 <- ((q01. * q0.0 * q.10) / (q0.. * q.1. * q..0)) * a...
a001 <- ((q00. * q0.1 * q.01) / (q0.. * q.0. * q..1)) * a...
a000 <- ((q00. * q0.0 * q.00) / (q0.. * q.0. * q..0)) * a...
# Step 3: Compute posterior probabilities (updating the prior)
gam111 <- a111 + n111
gam110 <- a110 + n110
gam101 <- a101 + n101
gam011 <- a011 + n011
gam100 <- a100 + n100
gam010 <- a010 + n010
gam001 <- a001 + n001
gam000 <- a000 + n000
# Sum of all gamma values
gamma_sum <- sum(gam111, gam110, gam101, gam011, gam100, gam010, gam001, gam000)
# Step 4: Compute joint probabilities
Ep111 <- gam111 / gamma_sum
Ep110 <- gam110 / gamma_sum
Ep101 <- gam101 / gamma_sum
Ep011 <- gam011 / gamma_sum
Ep100 <- gam100 / gamma_sum
Ep010 <- gam010 / gamma_sum
Ep001 <- gam001 / gamma_sum
Ep000 <- gam000 / gamma_sum
# Step 5: Compute marginal probabilities
Ep11. <- (gam111 + gam110) / gamma_sum
Ep1.1 <- (gam111 + gam101) / gamma_sum
Ep.11 <- (gam111 + gam011) / gamma_sum
Ep1.. <- (gam111 + gam110 + gam101 + gam100) / gamma_sum
Ep.1. <- (gam111 + gam110 + gam011 + gam010) / gamma_sum
Ep..1 <- (gam111 + gam101 + gam011 + gam001) / gamma_sum
# Step 6: Compute third-order IC MAP estimate
ic_map <- log2((Ep111 * Ep1.. * Ep.1. * Ep..1) / (Ep11. * Ep1.1 * Ep.11))
# Step 7: Compute ratios for third-order combinations, uncertainty adjustment
r111 <- round(gam111 / pmin(gam111 + gam110 + gam101, gam111 + gam011 + gam101, gam111 + gam110 + gam011), 1)
# Step 8: Look up Ar and Br based on r111
id <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)
Alist <- c(3.09, 2.93, 2.78, 2.62, 2.45, 2.25, 2.03, 1.79, 1.61, 1.13, 0.073)
Blist <- c(2.22, 2.27, 2.26, 2.25, 2.15, 2.12, 2.05, 1.93, 1.89, 1.15, -0.081)
Ar <- Alist[which(id == r111)]
Br <- Blist[which(id == r111)]
# Step 9: Compute Delta
delta <- Ar / sqrt(gam111) + Br * gam111^(-1.5)
# Step 10: Compute Confidence Intervals
IC_LB <- ic_map - delta
IC_UB <- ic_map + delta
# Step 11: Return results as a list
return(list(
IC = ic_map,
IC_LB = IC_LB,
IC_UB = IC_UB
))
}Step 2: Calling the Function to Generate Table and Result
Generating the Table
data <- fread(file_path)
data$pt <- toupper(data$pt)
setnames(data, old = names(data)[1], new = "id")
ic_counts <- as.data.table(get_ic_table(d1 = d1_study, d2 = d2_study, ae = ae_study, merged = data))
ic_counts[, `:=` (d1 = d1_study, d2 = d2_study, ae = ae_study)]
setcolorder(ic_counts, c("d1", "d2", "ae"))
list2env(ic_counts, envir = .GlobalEnv)## <environment: R_GlobalEnv>Computing the IC
ic_result <- as.data.table(get_ic_result(n111 = n111,
n110 = n110,
n101 = n101,
n011 = n011,
n100 = n100,
n010 = n010,
n001 = n001,
n000 = n000,
n... = n...))
ic_result[, `:=` (d1 = d1_study, d2 = d2_study, ae = ae_study)]
setcolorder(ic_result, c("d1", "d2", "ae"))DT::datatable(ic_counts, options = list(pageLength = 10, scrollX = TRUE))DT::datatable(ic_result, options = list(pageLength = 10, scrollX = TRUE))