R/IDAPredict.TestvsControl.R
In Alexander-Ling/IDACombo: Predict Drug Combination Efficacy using Independent Drug Action
Defines functions
IDAPredict.TestvsControl
Documented in
IDAPredict.TestvsControl
#' Predicts and compares IDA efficacies for a pair of control and test treatments
#'
#' This function creates efficacy predictions for a pair of control and test treatments, each treatment consisting of a combination of one or more drugs, using monotherapy efficacy data and the assumptions of independent drug action. Concentrations must be specified for each drug in each treatment.
#'
#' @importFrom stats complete.cases rnorm sd aggregate
#'
#' @param Monotherapy_Data A data frame where each row contains information about the response of a single cell line to a single drug at a single concentration. Must minimally include columns containing the following information: cell line name, drug name, drug concentration, and measured drug efficacy. May optionally include a column recording the standard error (SE) of the measured drug efficacy.
#' @param Cell_Line_Name_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains cell line names.
#' @param Drug_Name_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains drug names.
#' @param Drug_Concentration_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains drug concentrations.
#' @param Efficacy_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains measured drug efficacies (i.e. percent Viability, percent Cell Growth, etc.).
#' @param LowerEfficacyIsBetterDrugEffect A logic vector of length 1 indicating whether or not lower values in Efficacy_Column indicate a more effective drug effect (i.e. for percent viability). Set TRUE if so. Otherwise, set FALSE if higher values in Efficacy_Column indicate a more effective drug response (i.e. for percent cell death).
#' @param Efficacy_Metric_Name A character vector of length 1 indicating the name of the efficacy metric being used (i.e. Percent_Viability, Percent_Growth, etc.). Used to correctly label column names in output. Defaults to "Efficacy".
#' @param Control_Treatment_Drugs A character vector of length > 0 containing the names of the drugs in the control drug treatment for which efficacy predictions are to be made.
#' @param Control_Treatment_Drug_Concentrations A vector of drug concentrations for Control_Treatment_Drugs with the first concentration in Control_Treatment_Drug_Concentrations corresponding to the first drug in Control_Treatment_Drugs etc. Only one concentration may be specified for each drug in the control treatment, but, if a drug is included in both the control and test treatments, there is no need for the same concentration of that drug to be used in both treatments.
#' @param Test_Treatment_Drugs A character vector of length > 0 containing the names of the drugs in the control drug treatment for which efficacy predictions are to be made.
#' @param Test_Treatment_Drug_Concentrations A vector of drug concentrations for Test_Treatment_Drugs with the first concentration in Test_Treatment_Drug_Concentrations corresponding to the first drug in Test_Treatment_Drugs etc. Only one concentration may be specified for each drug in the test treatment, but, if a drug is included in both the control and test treatments, there is no need for the same concentration of that drug to be used in both treatments.
#' @param Calculate_Uncertainty A logic vector of length one indicating whether or not a semi-parametric bootstrap should be performed to estimate uncertainties in the efficacy predictions based on uncertainties in the monotherapy efficacy measurements. Set TRUE if you wish to calculate uncertainties. Defaults to FALSE.
#' @param Efficacy_SE_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains the standard errors of measured drug efficacies. Must be specified if Calculate_Uncertainty is set to TRUE.
#' @param n_Simulations A positive, integer vector of length 1 with a value >= 40 indicating the number of random samples to be drawn when calculating output efficacy prediction uncertainties. Defaults to 1000.
#' @param Calculate_Hazard_Ratio A logic vector of length 1 indicating whether or not a Hazard Ratios (HR) should be calculated between the control and test treatments. Set TRUE if so. Should only be set to TRUE for efficacy metrics that range between 0 and 1 (i.e. percent viability). Defaults to FALSE.
#' @param Average_Duplicate_Records A logic vector of length 1 indicating whether or not duplicated records (where a cell line has multiple records for being tested with a given drug at a given concentration) should be averaged. If TRUE, Efficacy values are averaged, and, if Calculate_Uncertainty is also TRUE, Efficacy_SE values are added in quadrature and divided by the number of duplicate records for that cell line/drug/concentration set.
#' @param Return_Bootstrap_Values A logic vector of length 1 indicating whether or not the function should return the Control Efficacies, Test Efficacies, and, if Calculate_Hazard_Ratio = TRUE, Hazard Ratios (HRs) simulated in the semi-parametric bootstrap used to estimate the uncertainties of those values. Defaults to FALSE. This parameter can only be set to TRUE if Calculate_Uncertainty = TRUE.
#'
#'@details
#'Uncertainty estimates for values calculated by this function are generated using a semi-parametric bootstrap approach. This is performed in several steps.\enumerate{
#'\item Control efficacies for each compound/concentration are simulated by random sampling from normal distributions with means equal to the provided calculated efficacies and standard deviations equal to the provided efficacy standard errors.
#'\item Test efficacies are simulated in the same fashion as the control efficacies, except in cases when a test drug also exists in the control therapy. In such cases, it is assumed that the efficacy values for this drug are derived from the same dose-response curve for the both the control and test therapies, so each simulated efficacy for that test drug is matched to the corresponding simulated efficacy from that drug in the control therapy using a standard normal deviate.
#'\item Efficacy predictions are made for the control and test therapies for each cell line and set of simulated efficacies using the assumptions of independent drug action.
#'\item Cell lines are randomly sampled with replacement for each simulation as many times as there are original cell lines. The simulated control and test therapy efficacies are then sampled according to the sampled cell lines for each simulation.
#'\item Mean efficacies are calculated for the control and test therapies for each simulation. If specified to do so, these values are then used to calculate simulated HRs.
#'\item The simulated distributions of each efficacy metric are used to estimate uncertainties for those metrics.
#'}
#'
#'@return \itemize{
#'\item If Return_Bootstrap_Values = FALSE, this function returns a list with 4 elements: 1) Either a data frame with the calculated efficacy predictions, or, if an error occurred, a character vector of length one with the error message. 2) A data frame listing the control treatment drug names and concentrations. 3) A data frame listing the test treatment drug names and concentrations. 4) A character vector containing the names of the cell lines used to make the efficacy predictions.
#'\item If Return_Bootstrap_Values = TRUE & Calculate_Uncertainty = TRUE & Calculate_Hazard_Ratio = FALSE, this function returns a list with 6 elements: the first 4 elements are the same as when Return_Bootstrap_Values = FALSE and the fifth and sixth elements being numeric vectors of, respectively, the control and test viabilities simulated during the semi-parametric bootstrap used to estimate uncertainties.
#'\item If Return_Bootstrap_Values = TRUE & Calculate_Uncertainty = TRUE & Calculate_Hazard_Ratio = TRUE, this function returns a list with 7 elements: the first 4 elements are the same as when Return_Bootstrap_Values = FALSE and the fifth, sixth, and seventh elements being numeric vectors of, respectively, the control viabilities, test viabilities, and HRs simulated during the semi-parametric bootstrap used to estimate uncertainties.
#'}
#'
#' @examples
#' #Loading Package
#' library(IDACombo)
#'
#' #Making fake monotherapy dataset
#' CellLineNames <- rep(c("CL1", "CL2", "CL3", "CL4", "CL5", "CL6"), 6)
#' DrugNames <- c(rep("D1", 12), rep("D2", 12), rep("D3", 12))
#' Concentrations <- c(rep(1, 6), rep(2, 6), rep(1.5, 6), rep(3, 6), rep("A", 6), rep("B", 6))
#' Viability <- c(sample(seq(0.4,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.2,0.8,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.4,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.2,0.6,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.9,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.2,0.6,length.out = 10), 6, replace = TRUE))
#' Viability_SE <- Viability * sample(seq(0,0.1,length.out = 100), 36, replace = TRUE)
#' Fake_Data <- data.frame(CellLineNames, DrugNames, Concentrations, Viability, Viability_SE)
#'
#' #Creating efficacy predictions for control and test treatments and comparing without
#' #uncertainty calculations
#' #For case where drugs in test treatment are at reduced concentrations from
#' #those used in the control treatment due to the addition of a third drug.
#' #Note that this may mean that the test treatment is less effective than
#' #the control treatment, such that the Hazard Ratio is > 1.
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(1, 1.5, "B"),
#' Calculate_Uncertainty = FALSE,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE)
#'
#' #For case where drugs in test treatment are at same concentrations as
#' #those used in the control treatment.
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = FALSE,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Creating efficacy predictions for control and test treatments and comparing with
#' #uncertainty calculations but without returning simulated values that are generated
#' #using a semi-parametric bootstrap to estimate uncertainties.
#'
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Viability_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Creating efficacy predictions for control and test treatments and comparing with
#' #uncertainty calculations and with returning simulated values that are generated
#' #using a semi-parametric bootstrap to estimate standard errors.
#'
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Viability_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE,
#' Return_Bootstrap_Values = TRUE)
#'
#' #Converting Viabilty to reduction in viability and redoing calculations
#' #without returning simulated values from semi-parametric boostrap. Note the change
#' #in the LowerEfficacyIsBetterDrugEffect flag from TRUE to FALSE
#' Reduction_in_Viability <- 1-Viability
#' Reduction_in_Viability_SE <- Viability_SE
#' Fake_Data <- data.frame(CellLineNames,
#' DrugNames,
#' Concentrations,
#' Reduction_in_Viability,
#' Reduction_in_Viability_SE)
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Reduction_in_Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Reduction_in_Viability_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = FALSE,
#' Efficacy_Metric_Name = "Reduction_in_Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Changing efficacy metric to percent growth (range -1 to 1)
#' #Note that calculating Hazard Ratios is no longer valid, so
#' #Calculate_Hazard_Ratio is set to FALSE.
#' Percent_Growth <- c(sample(seq(0.4,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(-0.4,0.2,length.out = 10), 6, replace = TRUE),
#' sample(seq(-0.2,0.3,length.out = 10), 6, replace = TRUE),
#' sample(seq(-1,0.2,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.8,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(-1,-0.2,length.out = 10), 6, replace = TRUE))
#' Percent_Growth_SE <- abs(Percent_Growth * sample(seq(0,0.1,length.out = 100), 36, replace = TRUE))
#' Fake_Data <- data.frame(CellLineNames,
#' DrugNames,
#' Concentrations,
#' Percent_Growth,
#' Percent_Growth_SE)
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Percent_Growth",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Percent_Growth_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Percent_Growth",
#' Calculate_Hazard_Ratio = FALSE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Adding duplicate records for each cell line, and showing behavior with
#' #Average_Duplicate_Records = FALSE. Should produce warning messages that
#' #duplicates were found and removed.
#' Percent_Growth <- c(sample(seq(0.4,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(-0.4,0.2,length.out = 10), 6, replace = TRUE),
#' sample(seq(-0.2,0.3,length.out = 10), 6, replace = TRUE),
#' sample(seq(-1,0.2,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.8,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(-1,-0.2,length.out = 10), 6, replace = TRUE))
#' Percent_Growth_SE <- abs(Percent_Growth * sample(seq(0,0.1,length.out = 100), 36, replace = TRUE))
#' Fake_Data_to_add <- data.frame(CellLineNames,
#' DrugNames,
#' Concentrations,
#' Percent_Growth,
#' Percent_Growth_SE)
#' Fake_Data <- rbind(Fake_Data, Fake_Data_to_add)
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Percent_Growth",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Percent_Growth_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Percent_Growth",
#' Calculate_Hazard_Ratio = FALSE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Now setting to average duplicate values.
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Percent_Growth",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Percent_Growth_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Percent_Growth",
#' Calculate_Hazard_Ratio = FALSE,
#' Average_Duplicate_Records = TRUE)
#'
#' @export
IDAPredict.TestvsControl <- function(Monotherapy_Data, Cell_Line_Name_Column, Drug_Name_Column, Drug_Concentration_Column, Efficacy_Column, LowerEfficacyIsBetterDrugEffect, Efficacy_Metric_Name = "Efficacy", Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations, Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations, Calculate_Uncertainty = FALSE, Efficacy_SE_Column = NULL, n_Simulations = 1000, Calculate_Hazard_Ratio = FALSE, Average_Duplicate_Records = FALSE, Return_Bootstrap_Values = FALSE){
#Checking that all input variables are in correct format
if(! is.data.frame(Monotherapy_Data)){
stop("Monotherapy_Data is not a data frame.")
}
if(! is.vector(Cell_Line_Name_Column) | ! is.character(Cell_Line_Name_Column) | ! length(Cell_Line_Name_Column) == 1){
stop("Cell_Line_Name_Column is not a character vector of length 1.")
}
if(! Cell_Line_Name_Column %in% colnames(Monotherapy_Data)){
stop("Cell_Line_Name_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(Drug_Name_Column) | ! is.character(Drug_Name_Column) | ! length(Drug_Name_Column) == 1){
stop("Drug_Name_Column is not a character vector of length 1.")
}
if(! Drug_Name_Column %in% colnames(Monotherapy_Data)){
stop("Drug_Name_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(Drug_Concentration_Column) | ! is.character(Drug_Concentration_Column) | ! length(Drug_Concentration_Column) == 1){
stop("Drug_Concentration_Column is not a character vector of length 1.")
}
if(! Drug_Concentration_Column %in% colnames(Monotherapy_Data)){
stop("Drug_Concentration_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(Efficacy_Column) | ! is.character(Efficacy_Column) | ! length(Efficacy_Column) == 1){
stop("Efficacy_Column is not a character vector of length 1.")
}
if(! Efficacy_Column %in% colnames(Monotherapy_Data)){
stop("Efficacy_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(LowerEfficacyIsBetterDrugEffect) | ! is.logical(LowerEfficacyIsBetterDrugEffect) | ! length(LowerEfficacyIsBetterDrugEffect) == 1){
stop("LowerEfficacyIsBetterDrugEffect is not a logical vector of length 1.")
}
if(! is.vector(Efficacy_Metric_Name) | ! is.character(Efficacy_Metric_Name) | ! length(Efficacy_Metric_Name) == 1){
stop("Efficacy_Metric_Name is not a character vector of length 1.")
}
if(! is.vector(Control_Treatment_Drugs) | ! is.character(Control_Treatment_Drugs) | ! length(Control_Treatment_Drugs) > 0){
stop("Control_Treatment_Drugs is not a character vector with length > 0.")
}
if(! all(Control_Treatment_Drugs %in% Monotherapy_Data[,Drug_Name_Column])){
missing.control.drugs <- Control_Treatment_Drugs[! Control_Treatment_Drugs %in% Monotherapy_Data[,Drug_Name_Column]]
stop(paste0("No data for the following Control_Treatment_Drugs found in Monotherapy_Data: ", paste(missing.control.drugs, collapse = ", ")))
}
if(! is.vector(Control_Treatment_Drug_Concentrations) | ! length(Control_Treatment_Drug_Concentrations) == length(Control_Treatment_Drugs)){
stop("Control_Treatment_Drug_Concentrations is not a vector with length = length(Control_Treatment_Drugs).")
}
if(! is.vector(Test_Treatment_Drugs) | ! is.character(Test_Treatment_Drugs) | ! length(Test_Treatment_Drugs) > 0){
stop("Test_Treatment_Drugs is not a character vector with length > 0.")
}
if(! all(Test_Treatment_Drugs %in% Monotherapy_Data[,Drug_Name_Column])){
missing.control.drugs <- Test_Treatment_Drugs[! Test_Treatment_Drugs %in% Monotherapy_Data[,Drug_Name_Column]]
stop(paste0("No data for the following Test_Treatment_Drugs found in Monotherapy_Data: ", paste(missing.control.drugs, collapse = ", ")))
}
if(! is.vector(Test_Treatment_Drug_Concentrations) | ! length(Test_Treatment_Drug_Concentrations) == length(Test_Treatment_Drugs)){
stop("Test_Treatment_Drug_Concentrations is not a vector with length = length(Test_Treatment_Drugs).")
}
if(! is.vector(Calculate_Uncertainty) | ! is.logical(Calculate_Uncertainty) | ! length(Calculate_Uncertainty) == 1){
stop("Calculate_Uncertainty is not a logical vector of length 1.")
}
if(Calculate_Uncertainty == TRUE){
if(is.null(Efficacy_SE_Column)){
stop("Calculate_Uncertainty is TRUE but Efficacy_SE_Column is not specified. Please either set Calculate_Uncertainty to FALSE or specify Efficacy_SE_Column.")
}
if(! is.vector(Efficacy_SE_Column) | ! is.character(Efficacy_SE_Column) | ! length(Efficacy_SE_Column) == 1){
stop("Calculate_Uncertainty is TRUE, but Efficacy_SE_Column is not a character vector of length 1.")
}
if(! Efficacy_SE_Column %in% colnames(Monotherapy_Data)){
stop("Calculate_Uncertainty is TRUE, but Efficacy_SE_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(n_Simulations) | ! is.numeric(n_Simulations) | ! length(n_Simulations) == 1){
stop("Calculate_Uncertainty is TRUE, but n_Simulations is not a numeric vector of length 1.")
}
if(! n_Simulations%%1==0 | ! n_Simulations >= 40){
stop("Calculate_Uncertainty is TRUE, but n_Simulations is not a positive integer >= 40.")
}
}
if(! is.vector(Return_Bootstrap_Values) | ! is.logical(Return_Bootstrap_Values) | ! length(Return_Bootstrap_Values) == 1){
stop("Return_Bootstrap_Values is not a logical vector of length 1.")
}
if(Return_Bootstrap_Values == TRUE & ! Calculate_Uncertainty == TRUE){
stop("Return_Bootstrap_Values is TRUE, but Calculate_Uncertainty is not TRUE.")
}
if(! is.vector(Calculate_Hazard_Ratio) | ! is.logical(Calculate_Hazard_Ratio) | ! length(Calculate_Hazard_Ratio) == 1){
stop("Calculate_Hazard_Ratio is not a logical vector of length 1.")
}
if(! is.vector(Average_Duplicate_Records) | ! is.logical(Average_Duplicate_Records) | ! length(Average_Duplicate_Records) == 1){
stop("Average_Duplicate_Records is not a logical vector of length 1.")
}
#Organizing data into standard format based on column names provided for each desired set of information
#Also subsetting to only include data pertaining to drugs in control or test treatment
if(Calculate_Uncertainty == TRUE){
Data <- Monotherapy_Data[Monotherapy_Data[,Drug_Name_Column] %in% c(Control_Treatment_Drugs, Test_Treatment_Drugs),c(Cell_Line_Name_Column, Drug_Name_Column, Drug_Concentration_Column, Efficacy_Column, Efficacy_SE_Column)]
colnames(Data) <- c("CellLine", "Drug", "Conc", "Efficacy", "Efficacy_SE")
} else {
Data <- Monotherapy_Data[Monotherapy_Data[,Drug_Name_Column] %in% c(Control_Treatment_Drugs, Test_Treatment_Drugs),c(Cell_Line_Name_Column, Drug_Name_Column, Drug_Concentration_Column, Efficacy_Column)]
colnames(Data) <- c("CellLine", "Drug", "Conc", "Efficacy")
}
rm(Monotherapy_Data)
#Making sure all columns are in correct formats
Data$CellLine <- as.character(Data$CellLine)
Data$Drug <- as.character(Data$Drug)
Data$Conc <- as.character(Data$Conc)
Data$Efficacy <- as.numeric(as.character(Data$Efficacy))
if(Calculate_Uncertainty == TRUE){
Data$Efficacy_SE <- as.numeric(as.character(Data$Efficacy_SE))
}
#Removing rows that are missing information
#Note: missing SE information is ignored if Calculate_Uncertainty == FALSE
Data <- Data[complete.cases(Data),]
#Subsetting into control and test group data with specified concentrations
ControlData <- list(NULL)
for(i in 1:length(Control_Treatment_Drugs)){
ControlData[[i]] <- Data[Data$Drug %in% Control_Treatment_Drugs[i],]
#Checking that all provided concentrations are available for their respective drugs
if(! Control_Treatment_Drug_Concentrations[i] %in% ControlData[[i]]$Conc){
stop(paste0(Control_Treatment_Drug_Concentrations[i], " concentration is unavailable for ", Control_Treatment_Drugs[i], " in control treatment."))
} else {
ControlData[[i]] <- ControlData[[i]][ControlData[[i]]$Conc %in% Control_Treatment_Drug_Concentrations[i],]
}
}
names(ControlData) <- Control_Treatment_Drugs
TestData <- list(NULL)
for(i in 1:length(Test_Treatment_Drugs)){
TestData[[i]] <- Data[Data$Drug %in% Test_Treatment_Drugs[i],]
#Checking that all provided concentrations are available for their respective drugs
if(! Test_Treatment_Drug_Concentrations[i] %in% TestData[[i]]$Conc){
stop(paste0(Test_Treatment_Drug_Concentrations[i], " concentration is unavailable for ", Test_Treatment_Drugs[i], " in test treatment."))
} else {
TestData[[i]] <- TestData[[i]][TestData[[i]]$Conc %in% Test_Treatment_Drug_Concentrations[i],]
}
}
names(TestData) <- Test_Treatment_Drugs
rm(Data)
#Finding cell line overlap between all drugs
ControlCellLines <- list(NULL)
for(i in 1:length(ControlData)){
ControlCellLines[[i]] <- sort(unique(ControlData[[i]]$CellLine))
}
TestCellLines <- list(NULL)
for(i in 1:length(TestData)){
TestCellLines[[i]] <- sort(unique(TestData[[i]]$CellLine))
}
All_Drug_CellLines <- c(ControlCellLines, TestCellLines)
Usable_CellLines <- sort(unique(unlist(All_Drug_CellLines)))
for(i in 1:length(All_Drug_CellLines)){
Usable_CellLines <- Usable_CellLines[Usable_CellLines %in% All_Drug_CellLines[[i]]]
}
rm(ControlCellLines, TestCellLines, All_Drug_CellLines)
#Checking again if at least 2 cell lines remain for all drugs. If not, exiting with no
#predictions and a warning.
if(! length(Usable_CellLines) >= 2){
#Returning NA predictions with warning due to too few cell lines.
warning(paste0("<2 overlapping cell lines available for comparison of (", paste(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations, sep = "_", collapse = " + "), ") vs. (", paste(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations, sep = "_", collapse = " + "), ")"))
if(Return_Bootstrap_Values == FALSE){
Return_Object <- list("Less than 2 overlapping cell lines available.", as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used")
return(Return_Object)
} else if(Return_Bootstrap_Values == TRUE & Calculate_Hazard_Ratio == TRUE){
Return_Object <- list("Less than 2 overlapping cell lines available.", as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines, NULL, NULL, NULL)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used", "Bootstrap_Mean_Control_Efficacies", "Bootstrap_Mean_Test_Efficacies", "Bootstrap_HRs")
return(Return_Object)
} else if(Return_Bootstrap_Values == TRUE & Calculate_Hazard_Ratio == FALSE){
Return_Object <- list("Less than 2 overlapping cell lines available.", as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines, NULL, NULL)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used", "Bootstrap_Mean_Control_Efficacies", "Bootstrap_Mean_Test_Efficacies")
return(Return_Object)
}
}
#Subsetting drug data to only include overlapping cell lines
for(i in 1:length(ControlData)){
ControlData[[i]] <- ControlData[[i]][ControlData[[i]]$CellLine %in% Usable_CellLines,]
}
for(i in 1:length(TestData)){
TestData[[i]] <- TestData[[i]][TestData[[i]]$CellLine %in% Usable_CellLines,]
}
#Checking that cell lines aren't duplicated in each drug dataset
#If Average_Duplicate_Records == FALSE, removing cell line duplicates with warning if duplicates are found.
#If Average_Duplicate_Records == TRUE, averaging duplicate records without warning.
for(i in 1:length(Control_Treatment_Drugs)){
CL_Conc <- paste(ControlData[[i]]$CellLine, ControlData[[i]]$Conc, sep = "_")
Dups <- CL_Conc[duplicated(CL_Conc)]
if(length(Dups) > 0 & Average_Duplicate_Records == FALSE){
warning(paste0("Duplicated information found for the following cell lines and ", Control_Treatment_Drugs[i], " concentrations in the control treatment. Average_Duplicate_Records = FALSE so duplicates removed: ", paste(Dups, collapse = ", ")))
ControlData[[i]] <- ControlData[[i]][! duplicated(CL_Conc),]
} else if(length(Dups) > 0 & Average_Duplicate_Records == TRUE){
if(Calculate_Uncertainty == FALSE){
#Simply averaging efficacy values
colnames <- colnames(ControlData[[i]])
ControlData[[i]] <- aggregate(ControlData[[i]]$Efficacy, by = list(ControlData[[i]]$CellLine, ControlData[[i]]$Drug, ControlData[[i]]$Conc), FUN = mean)
colnames(ControlData[[i]]) <- colnames
} else if(Calculate_Uncertainty == TRUE){
#Averaging efficacy
Efficacy_Average <- aggregate(ControlData[[i]]$Efficacy, by = list(ControlData[[i]]$CellLine, ControlData[[i]]$Drug, ControlData[[i]]$Conc), FUN = mean)
colnames(Efficacy_Average) <- c("CellLine", "Drug", "Conc", "Efficacy")
#Calculating uncertainty in averaged efficacy by adding efficacy uncertainties in quadrature and dividing by number of values used in average
Efficacy_Average_Uncertainties <- aggregate(ControlData[[i]]$Efficacy_SE, by = list(ControlData[[i]]$CellLine, ControlData[[i]]$Drug, ControlData[[i]]$Conc), FUN = function(x){sqrt(sum(x^2))/length(x)})
colnames(Efficacy_Average_Uncertainties) <- c("CellLine", "Drug", "Conc", "Efficacy_SE")
#Combining results
ControlData[[i]] <- merge(Efficacy_Average, Efficacy_Average_Uncertainties)
}
}
}
for(i in 1:length(Test_Treatment_Drugs)){
CL_Conc <- paste(TestData[[i]]$CellLine, TestData[[i]]$Conc, sep = "_")
Dups <- CL_Conc[duplicated(CL_Conc)]
if(length(Dups) > 0 & Average_Duplicate_Records == FALSE){
warning(paste0("Duplicated information found for the following cell lines and ", Test_Treatment_Drugs[i], " concentrations in the test treatment. Average_Duplicate_Records = FALSE so duplicates removed: ", paste(Dups, collapse = ", ")))
TestData[[i]] <- TestData[[i]][! duplicated(CL_Conc),]
} else if(length(Dups) > 0 & Average_Duplicate_Records == TRUE){
if(Calculate_Uncertainty == FALSE){
#Simply averaging efficacy values
colnames <- colnames(TestData[[i]])
TestData[[i]] <- aggregate(TestData[[i]]$Efficacy, by = list(TestData[[i]]$CellLine, TestData[[i]]$Drug, TestData[[i]]$Conc), FUN = mean)
colnames(TestData[[i]]) <- colnames
} else if(Calculate_Uncertainty == TRUE){
#Averaging efficacy
Efficacy_Average <- aggregate(TestData[[i]]$Efficacy, by = list(TestData[[i]]$CellLine, TestData[[i]]$Drug, TestData[[i]]$Conc), FUN = mean)
colnames(Efficacy_Average) <- c("CellLine", "Drug", "Conc", "Efficacy")
#Calculating uncertainty in averaged efficacy by adding efficacy uncertainties in quadrature and dividing by number of values used in average
Efficacy_Average_Uncertainties <- aggregate(TestData[[i]]$Efficacy_SE, by = list(TestData[[i]]$CellLine, TestData[[i]]$Drug, TestData[[i]]$Conc), FUN = function(x){sqrt(sum(x^2))/length(x)})
colnames(Efficacy_Average_Uncertainties) <- c("CellLine", "Drug", "Conc", "Efficacy_SE")
#Combining results
TestData[[i]] <- merge(Efficacy_Average, Efficacy_Average_Uncertainties)
}
}
}
#Ordering cell lines the same for each drug
for(i in 1:length(ControlData)){
ControlData[[i]] <- ControlData[[i]][order(ControlData[[i]]$Conc, ControlData[[i]]$CellLine),]
}
for(i in 1:length(TestData)){
TestData[[i]] <- TestData[[i]][order(TestData[[i]]$Conc, TestData[[i]]$CellLine),]
}
#Creating object to store prediction results in
Prediction_Results <- as.data.frame(matrix(rep(NA, 11), nrow = 1))
colnames(Prediction_Results) <- c("Mean_Control_Treatment_Efficacy", "Mean_Control_Treatment_Efficacy_SE", "Mean_Control_Treatment_Efficacy_95%_Confidence_Interval", "Mean_Test_Treatment_Efficacy", "Mean_Test_Treatment_Efficacy_SE", "Mean_Test_Treatment_Efficacy_95%_Confidence_Interval", "HR_Test_vs_Control_Treatment", "HR_Test_vs_Control_Treatment_SE", "HR_Test_vs_Control_Treatment_95%_Confidence_Interval", "p_HR>=1", "p_HR=1")
#Performing combination efficacy predictions for cases where lower efficacy values
#indicate a more effective drug effect (i.e. efficacy = percent viability, percent growth, etc.)
if(LowerEfficacyIsBetterDrugEffect == TRUE){
#Calculating expected combination efficacy for each cell line
#using Independent Drug Action for both Control and Test treatment
Control_Efficacy <- do.call(pmin, lapply(ControlData, function(x){return(x$Efficacy)}))
Test_Efficacy <- do.call(pmin, lapply(TestData, function(x){return(x$Efficacy)}))
#Calculating average efficacy across all cell lines
Prediction_Results$Mean_Control_Treatment_Efficacy <- mean(Control_Efficacy)
Prediction_Results$Mean_Test_Treatment_Efficacy <- mean(Test_Efficacy)
if(Calculate_Hazard_Ratio == TRUE){
Prediction_Results$HR_Test_vs_Control_Treatment <- Prediction_Results$Mean_Test_Treatment_Efficacy / Prediction_Results$Mean_Control_Treatment_Efficacy
}
rm(Control_Efficacy, Test_Efficacy)
#If Calculate_Uncertainty == TRUE, doing Monte Carlo simulation to estimate
#uncertainties in output parameters based on uncertainties in monotherapy efficacies.
if(Calculate_Uncertainty == TRUE){
#Looping through each drug for each treatment and simulating efficacies based on
#measured efficacies and SE's
Control_MC_Efficacies <- as.list(NULL)
for(i in 1:length(ControlData)){
Control_MC_Efficacies[[i]] <- apply(ControlData[[i]][,c("Efficacy", "Efficacy_SE")], 1, function(x){rnorm(n = n_Simulations, mean = x[1], sd = x[2])})
colnames(Control_MC_Efficacies[[i]]) <- ControlData[[i]]$CellLine
}
names(Control_MC_Efficacies) <- names(ControlData)
Test_MC_Efficacies <- as.list(NULL)
for(i in 1:length(TestData)){
#Checking if this test drug is also in the control therapy.
#If so, not randomly sampling for test viabilities. Instead, calculating how many SE's from
#the measured value each simulated control viability fell, and matching that distance in
#the simulated test viabilities. This is done because the drug in the control and test
#therapies are not independent in this case--they will have come from the same dose-response curve.
if(names(TestData)[i] %in% names(ControlData)){
#Calculating the number of SEs away from the measured value each simulated value is in the control therapy for this drug
Measured_Control_Data <- ControlData[[names(TestData)[i]]][,c("CellLine", "Efficacy", "Efficacy_SE")]
Measured_Control_Data <- Measured_Control_Data[match(colnames(Control_MC_Efficacies[[names(TestData)[i]]]), Measured_Control_Data$CellLine),]
Measured_Control_Efficacies <- matrix(Measured_Control_Data$Efficacy, ncol = length(Measured_Control_Data$Efficacy), nrow = n_Simulations, byrow = TRUE)
Measured_Control_Efficacy_SEs <- matrix(Measured_Control_Data$Efficacy_SE, ncol = length(Measured_Control_Data$Efficacy_SE), nrow = n_Simulations, byrow = TRUE)
SEs_deviated <- (Control_MC_Efficacies[[names(TestData)[i]]] - Measured_Control_Efficacies) / Measured_Control_Efficacy_SEs
#Calculating simulated test viabilities based on the SE deviations from the simulated control viabilities for this drug
Measured_Test_Data <- TestData[[names(TestData)[i]]][,c("CellLine", "Efficacy", "Efficacy_SE")]
Measured_Test_Data <- Measured_Test_Data[match(colnames(SEs_deviated), Measured_Test_Data$CellLine),]
Measured_Test_Efficacies <- matrix(Measured_Test_Data$Efficacy, ncol = length(Measured_Test_Data$Efficacy), nrow = n_Simulations, byrow = TRUE, dimnames = list(NULL, Measured_Test_Data$CellLine))
Measured_Test_Efficacy_SEs <- matrix(Measured_Test_Data$Efficacy_SE, ncol = length(Measured_Test_Data$Efficacy_SE), nrow = n_Simulations, byrow = TRUE, dimnames = list(NULL, Measured_Test_Data$CellLine))
Test_MC_Efficacies[[i]] <- Measured_Test_Efficacies + (SEs_deviated * Measured_Test_Efficacy_SEs)
} else {
#If test drug is not in control therapy, randomly sampling
Test_MC_Efficacies[[i]] <- apply(TestData[[i]][,c("Efficacy", "Efficacy_SE")], 1, function(x){rnorm(n = n_Simulations, mean = x[1], sd = x[2])})
}
}
names(Test_MC_Efficacies) <- names(TestData)
if(exists("Measured_Control_Data")){
rm(Measured_Control_Data, Measured_Control_Efficacies, Measured_Control_Efficacy_SEs, Measured_Test_Data, Measured_Test_Efficacies, Measured_Test_Efficacy_SEs)
}
#Calculating expected combination efficacies for each cell line
#using Independent Drug Action for both Control and Test treatment
Control_MC_Combo_Efficacies <- do.call(pmin, Control_MC_Efficacies)
Test_MC_Combo_Efficacies <- do.call(pmin, Test_MC_Efficacies)
#Resampling cell lines for each simulation with replacement to account for the effect of random variation in cell line selection
Selected_CCLs_Per_Sim <- t(apply(Control_MC_Combo_Efficacies, 1, function(x){return(sample(1:length(x), replace = TRUE))}))
Control_MC_Combo_Efficacies <- t(mapply(function(x,y){return(x[y])}, x = split(Control_MC_Combo_Efficacies, row(Control_MC_Combo_Efficacies)), y = split(Selected_CCLs_Per_Sim, row(Selected_CCLs_Per_Sim))))
Test_MC_Combo_Efficacies <- t(mapply(function(x,y){return(x[y])}, x = split(Test_MC_Combo_Efficacies, row(Test_MC_Combo_Efficacies)), y = split(Selected_CCLs_Per_Sim, row(Selected_CCLs_Per_Sim))))
#Calculating average efficacy across all cell lines
Mean_Control_Efficacies_MC <- rowMeans(Control_MC_Combo_Efficacies)
Mean_Test_Efficacies_MC <- rowMeans(Test_MC_Combo_Efficacies)
#Calculating standard errors
Prediction_Results$Mean_Control_Treatment_Efficacy_SE <- sd(Mean_Control_Efficacies_MC)
Prediction_Results$Mean_Test_Treatment_Efficacy_SE <- sd(Mean_Test_Efficacies_MC)
#Calculating 95% confidence intervals for mean treatment efficacies
index_2.5 <- floor(0.025*n_Simulations)
index_97.5 <- ceiling(0.975*n_Simulations)
Prediction_Results$`Mean_Control_Treatment_Efficacy_95%_Confidence_Interval` <- paste(sort(Mean_Control_Efficacies_MC, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
Prediction_Results$`Mean_Test_Treatment_Efficacy_95%_Confidence_Interval` <- paste(sort(Mean_Test_Efficacies_MC, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
#Calculating HR if specified to do so
if(Calculate_Hazard_Ratio == TRUE){
MC_HRs <- Mean_Test_Efficacies_MC / Mean_Control_Efficacies_MC
Prediction_Results$HR_Test_vs_Control_Treatment_SE <- sd(MC_HRs)
#Calculating HR 95% confidence interval
Prediction_Results$`HR_Test_vs_Control_Treatment_95%_Confidence_Interval` <- paste(sort(MC_HRs, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
#Calculating HR p values
n_equal_to_1 <- sum(MC_HRs == 1)
n_less_than_1 <- sum(MC_HRs < 1)
n_greater_than_1 <- sum(MC_HRs > 1)
lower_tail <- (n_less_than_1 + 0.5*n_equal_to_1)/length(MC_HRs)
upper_tail <- (n_greater_than_1 + 0.5*n_equal_to_1)/length(MC_HRs)
#Calculating one-sided p-value with null hypothesis that HR >= 1
Prediction_Results$`p_HR>=1` <- sum(MC_HRs >= 1) / length(MC_HRs)
#If p-value is 0, setting as p < minimum p value that can be estimated using this many simulations
if(Prediction_Results$`p_HR>=1` == 0){Prediction_Results$`p_HR>=1` <- paste0("<", 1/n_Simulations)}
#Calculating two-sided p-value with null hypothesis that HR = 1
Prediction_Results$`p_HR=1` <- min(lower_tail, upper_tail)*2
#If p-value is 0, setting as p < minimum p value that can be estimated using this many simulations
if(Prediction_Results$`p_HR=1` == 0){Prediction_Results$`p_HR=1` <- paste0("<", 1/n_Simulations)}
}
#Cleaning up
rm(Control_MC_Combo_Efficacies, Test_MC_Combo_Efficacies, Selected_CCLs_Per_Sim, index_2.5, index_97.5)
if(Return_Bootstrap_Values == FALSE){
rm(Mean_Test_Efficacies_MC, Mean_Control_Efficacies_MC)
}
if(Calculate_Hazard_Ratio == TRUE){
rm(n_equal_to_1, n_less_than_1, n_greater_than_1, lower_tail, upper_tail)
if(Return_Bootstrap_Values == FALSE){
rm(MC_HRs)
}
}
}
}
#Performing combination efficacy predictions for cases where lower efficacy values
#indicate a less effective drug effect (i.e. efficacy = percent cell death, etc.)
if(LowerEfficacyIsBetterDrugEffect == FALSE){
#Calculating expected combination efficacy for each cell line
#using Independent Drug Action for both Control and Test treatment
Control_Efficacy <- do.call(pmax, lapply(ControlData, function(x){return(x$Efficacy)}))
Test_Efficacy <- do.call(pmax, lapply(TestData, function(x){return(x$Efficacy)}))
#Calculating average efficacy across all cell lines
Prediction_Results$Mean_Control_Treatment_Efficacy <- mean(Control_Efficacy)
Prediction_Results$Mean_Test_Treatment_Efficacy <- mean(Test_Efficacy)
if(Calculate_Hazard_Ratio == TRUE){
Prediction_Results$HR_Test_vs_Control_Treatment <- (1 - Prediction_Results$Mean_Test_Treatment_Efficacy) / (1 - Prediction_Results$Mean_Control_Treatment_Efficacy)
}
rm(Control_Efficacy, Test_Efficacy)
#If Calculate_Uncertainty == TRUE, doing Monte Carlo simulation to estimate
#uncertainties in output parameters based on uncertainties in monotherapy efficacies.
if(Calculate_Uncertainty == TRUE){
#Looping through each drug for each treatment and simulating efficacies based on
#measured efficacies and SE's
Control_MC_Efficacies <- as.list(NULL)
for(i in 1:length(ControlData)){
Control_MC_Efficacies[[i]] <- apply(ControlData[[i]][,c("Efficacy", "Efficacy_SE")], 1, function(x){rnorm(n = n_Simulations, mean = x[1], sd = x[2])})
colnames(Control_MC_Efficacies[[i]]) <- ControlData[[i]]$CellLine
}
names(Control_MC_Efficacies) <- names(ControlData)
Test_MC_Efficacies <- as.list(NULL)
for(i in 1:length(TestData)){
#Checking if this test drug is also in the control therapy.
#If so, not randomly sampling for test viabilities. Instead, calculating how many SE's from
#the measured value each simulated control viability fell, and matching that distance in
#the simulated test viabilities. This is done because the drug in the control and test
#therapies are not independent in this case--they will have come from the same dose-response curve.
if(names(TestData)[i] %in% names(ControlData)){
#Calculating the number of SEs away from the measured value each simulated value is in the control therapy for this drug
Measured_Control_Data <- ControlData[[names(TestData)[i]]][,c("CellLine", "Efficacy", "Efficacy_SE")]
Measured_Control_Data <- Measured_Control_Data[match(colnames(Control_MC_Efficacies[[names(TestData)[i]]]), Measured_Control_Data$CellLine),]
Measured_Control_Efficacies <- matrix(Measured_Control_Data$Efficacy, ncol = length(Measured_Control_Data$Efficacy), nrow = n_Simulations, byrow = TRUE)
Measured_Control_Efficacy_SEs <- matrix(Measured_Control_Data$Efficacy_SE, ncol = length(Measured_Control_Data$Efficacy_SE), nrow = n_Simulations, byrow = TRUE)
SEs_deviated <- (Control_MC_Efficacies[[names(TestData)[i]]] - Measured_Control_Efficacies) / Measured_Control_Efficacy_SEs
#Calculating simulated test viabilities based on the SE deviations from the simulated control viabilities for this drug
Measured_Test_Data <- TestData[[names(TestData)[i]]][,c("CellLine", "Efficacy", "Efficacy_SE")]
Measured_Test_Data <- Measured_Test_Data[match(colnames(SEs_deviated), Measured_Test_Data$CellLine),]
Measured_Test_Efficacies <- matrix(Measured_Test_Data$Efficacy, ncol = length(Measured_Test_Data$Efficacy), nrow = n_Simulations, byrow = TRUE, dimnames = list(NULL, Measured_Test_Data$CellLine))
Measured_Test_Efficacy_SEs <- matrix(Measured_Test_Data$Efficacy_SE, ncol = length(Measured_Test_Data$Efficacy_SE), nrow = n_Simulations, byrow = TRUE, dimnames = list(NULL, Measured_Test_Data$CellLine))
Test_MC_Efficacies[[i]] <- Measured_Test_Efficacies + (SEs_deviated * Measured_Test_Efficacy_SEs)
}else {
#If test drug is not in control therapy, randomly sampling
Test_MC_Efficacies[[i]] <- apply(TestData[[i]][,c("Efficacy", "Efficacy_SE")], 1, function(x){rnorm(n = n_Simulations, mean = x[1], sd = x[2])})
}
}
names(Test_MC_Efficacies) <- names(TestData)
if(exists("Measured_Control_Data")){
rm(Measured_Control_Data, Measured_Control_Efficacies, Measured_Control_Efficacy_SEs, Measured_Test_Data, Measured_Test_Efficacies, Measured_Test_Efficacy_SEs)
}
#Calculating expected combination efficacies for each cell line
#using Independent Drug Action for both Control and Test treatment
Control_MC_Combo_Efficacies <- do.call(pmax, Control_MC_Efficacies)
Test_MC_Combo_Efficacies <- do.call(pmax, Test_MC_Efficacies)
#Resampling cell lines for each simulation with replacement to account for the effect of random variation in cell line selection
Selected_CCLs_Per_Sim <- t(apply(Control_MC_Combo_Efficacies, 1, function(x){return(sample(1:length(x), replace = TRUE))}))
Control_MC_Combo_Efficacies <- t(mapply(function(x,y){return(x[y])}, x = split(Control_MC_Combo_Efficacies, row(Control_MC_Combo_Efficacies)), y = split(Selected_CCLs_Per_Sim, row(Selected_CCLs_Per_Sim))))
Test_MC_Combo_Efficacies <- t(mapply(function(x,y){return(x[y])}, x = split(Test_MC_Combo_Efficacies, row(Test_MC_Combo_Efficacies)), y = split(Selected_CCLs_Per_Sim, row(Selected_CCLs_Per_Sim))))
#Calculating average efficacy across all cell lines
Mean_Control_Efficacies_MC <- rowMeans(Control_MC_Combo_Efficacies)
Mean_Test_Efficacies_MC <- rowMeans(Test_MC_Combo_Efficacies)
#Calculating standard errors
Prediction_Results$Mean_Control_Treatment_Efficacy_SE <- sd(Mean_Control_Efficacies_MC)
Prediction_Results$Mean_Test_Treatment_Efficacy_SE <- sd(Mean_Test_Efficacies_MC)
#Calculating 95% confidence intervals for mean treatment efficacies
index_2.5 <- floor(0.025*n_Simulations)
index_97.5 <- ceiling(0.975*n_Simulations)
Prediction_Results$`Mean_Control_Treatment_Efficacy_95%_Confidence_Interval` <- paste(sort(Mean_Control_Efficacies_MC, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
Prediction_Results$`Mean_Test_Treatment_Efficacy_95%_Confidence_Interval` <- paste(sort(Mean_Test_Efficacies_MC, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
#Calculating HRs if specified to do so
if(Calculate_Hazard_Ratio == TRUE){
MC_HRs <- (1 - Mean_Test_Efficacies_MC) / (1 - Mean_Control_Efficacies_MC)
Prediction_Results$HR_Test_vs_Control_Treatment_SE <- sd(MC_HRs)
#Calculating HR 95% confidence interval
Prediction_Results$`HR_Test_vs_Control_Treatment_95%_Confidence_Interval` <- paste(sort(MC_HRs, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
#Calculating HR p values
n_equal_to_1 <- sum(MC_HRs == 1)
n_less_than_1 <- sum(MC_HRs < 1)
n_greater_than_1 <- sum(MC_HRs > 1)
lower_tail <- (n_less_than_1 + 0.5*n_equal_to_1)/length(MC_HRs)
upper_tail <- (n_greater_than_1 + 0.5*n_equal_to_1)/length(MC_HRs)
#Calculating one-sided p-value with null hypothesis that HR >= 1
Prediction_Results$`p_HR>=1` <- sum(MC_HRs >= 1) / length(MC_HRs)
#If p-value is 0, setting as p < minimum p value that can be estimated using this many simulations
if(Prediction_Results$`p_HR>=1` == 0){Prediction_Results$`p_HR>=1` <- paste0("<", 1/n_Simulations)}
#Calculating two-sided p-value with null hypothesis that HR = 1
Prediction_Results$`p_HR=1` <- min(lower_tail, upper_tail)*2
#If p-value is 0, setting as p < minimum p value that can be estimated using this many simulations
if(Prediction_Results$`p_HR=1` == 0){Prediction_Results$`p_HR=1` <- paste0("<", 1/n_Simulations)}
}
#Cleaning up
rm(Control_MC_Combo_Efficacies, Test_MC_Combo_Efficacies, Selected_CCLs_Per_Sim, index_2.5, index_97.5)
if(Return_Bootstrap_Values == FALSE){
rm(Mean_Test_Efficacies_MC, Mean_Control_Efficacies_MC)
}
if(Calculate_Hazard_Ratio == TRUE){
rm(n_equal_to_1, n_less_than_1, n_greater_than_1, lower_tail, upper_tail)
if(Return_Bootstrap_Values == FALSE){
rm(MC_HRs)
}
}
}
}
#Returning Outputs
#If Calculate_Uncertainty == FALSE, removing SE columns
if(Calculate_Uncertainty == FALSE){
Prediction_Results <- Prediction_Results[,-which(colnames(Prediction_Results) %in% c("Mean_Control_Treatment_Efficacy_SE", "Mean_Control_Treatment_Efficacy_95%_Confidence_Interval", "Mean_Test_Treatment_Efficacy_SE", "Mean_Test_Treatment_Efficacy_95%_Confidence_Interval", "HR_Test_vs_Control_Treatment_SE", "HR_Test_vs_Control_Treatment_95%_Confidence_Interval", "p_HR>=1", "p_HR=1"))]
}
#If Calculate_Hazard_Ratio == FALSE, removing HR columns
if(Calculate_Hazard_Ratio == FALSE){
Prediction_Results <- Prediction_Results[,-which(colnames(Prediction_Results) %in% c("HR_Test_vs_Control_Treatment", "HR_Test_vs_Control_Treatment_SE", "HR_Test_vs_Control_Treatment_95%_Confidence_Interval", "p_HR>=1", "p_HR=1"))]
}
#Replacing "Efficacy" with Efficacy_Metric_Name in column names of Dose_Comparisons
colnames(Prediction_Results) <- gsub("Efficacy", Efficacy_Metric_Name, colnames(Prediction_Results))
#Constructing Return_Object
if(Return_Bootstrap_Values == FALSE){
Return_Object <- list(Prediction_Results, as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used")
} else if(Return_Bootstrap_Values == TRUE & Calculate_Hazard_Ratio == FALSE & Calculate_Uncertainty == TRUE){
Return_Object <- list(Prediction_Results, as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines, Mean_Control_Efficacies_MC, Mean_Test_Efficacies_MC)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used", "Bootstrap_Mean_Control_Efficacies", "Bootstrap_Mean_Test_Efficacies")
} else if(Return_Bootstrap_Values == TRUE & Calculate_Hazard_Ratio == TRUE & Calculate_Uncertainty == TRUE){
Return_Object <- list(Prediction_Results, as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines, Mean_Control_Efficacies_MC, Mean_Test_Efficacies_MC, MC_HRs)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used", "Bootstrap_Mean_Control_Efficacies", "Bootstrap_Mean_Test_Efficacies", "Bootstrap_HRs")
}
#Returning output
return(Return_Object)
}
Alexander-Ling/IDACombo documentation built on July 31, 2020, 7:31 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions IDAPredict.TestvsControl
Documented in IDAPredict.TestvsControl
#' Predicts and compares IDA efficacies for a pair of control and test treatments
#'
#' This function creates efficacy predictions for a pair of control and test treatments, each treatment consisting of a combination of one or more drugs, using monotherapy efficacy data and the assumptions of independent drug action. Concentrations must be specified for each drug in each treatment.
#'
#' @importFrom stats complete.cases rnorm sd aggregate
#'
#' @param Monotherapy_Data A data frame where each row contains information about the response of a single cell line to a single drug at a single concentration. Must minimally include columns containing the following information: cell line name, drug name, drug concentration, and measured drug efficacy. May optionally include a column recording the standard error (SE) of the measured drug efficacy.
#' @param Cell_Line_Name_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains cell line names.
#' @param Drug_Name_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains drug names.
#' @param Drug_Concentration_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains drug concentrations.
#' @param Efficacy_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains measured drug efficacies (i.e. percent Viability, percent Cell Growth, etc.).
#' @param LowerEfficacyIsBetterDrugEffect A logic vector of length 1 indicating whether or not lower values in Efficacy_Column indicate a more effective drug effect (i.e. for percent viability). Set TRUE if so. Otherwise, set FALSE if higher values in Efficacy_Column indicate a more effective drug response (i.e. for percent cell death).
#' @param Efficacy_Metric_Name A character vector of length 1 indicating the name of the efficacy metric being used (i.e. Percent_Viability, Percent_Growth, etc.). Used to correctly label column names in output. Defaults to "Efficacy".
#' @param Control_Treatment_Drugs A character vector of length > 0 containing the names of the drugs in the control drug treatment for which efficacy predictions are to be made.
#' @param Control_Treatment_Drug_Concentrations A vector of drug concentrations for Control_Treatment_Drugs with the first concentration in Control_Treatment_Drug_Concentrations corresponding to the first drug in Control_Treatment_Drugs etc. Only one concentration may be specified for each drug in the control treatment, but, if a drug is included in both the control and test treatments, there is no need for the same concentration of that drug to be used in both treatments.
#' @param Test_Treatment_Drugs A character vector of length > 0 containing the names of the drugs in the control drug treatment for which efficacy predictions are to be made.
#' @param Test_Treatment_Drug_Concentrations A vector of drug concentrations for Test_Treatment_Drugs with the first concentration in Test_Treatment_Drug_Concentrations corresponding to the first drug in Test_Treatment_Drugs etc. Only one concentration may be specified for each drug in the test treatment, but, if a drug is included in both the control and test treatments, there is no need for the same concentration of that drug to be used in both treatments.
#' @param Calculate_Uncertainty A logic vector of length one indicating whether or not a semi-parametric bootstrap should be performed to estimate uncertainties in the efficacy predictions based on uncertainties in the monotherapy efficacy measurements. Set TRUE if you wish to calculate uncertainties. Defaults to FALSE.
#' @param Efficacy_SE_Column A character vector of length 1 containing the name of the column in the Monotherapy_Data data frame which contains the standard errors of measured drug efficacies. Must be specified if Calculate_Uncertainty is set to TRUE.
#' @param n_Simulations A positive, integer vector of length 1 with a value >= 40 indicating the number of random samples to be drawn when calculating output efficacy prediction uncertainties. Defaults to 1000.
#' @param Calculate_Hazard_Ratio A logic vector of length 1 indicating whether or not a Hazard Ratios (HR) should be calculated between the control and test treatments. Set TRUE if so. Should only be set to TRUE for efficacy metrics that range between 0 and 1 (i.e. percent viability). Defaults to FALSE.
#' @param Average_Duplicate_Records A logic vector of length 1 indicating whether or not duplicated records (where a cell line has multiple records for being tested with a given drug at a given concentration) should be averaged. If TRUE, Efficacy values are averaged, and, if Calculate_Uncertainty is also TRUE, Efficacy_SE values are added in quadrature and divided by the number of duplicate records for that cell line/drug/concentration set.
#' @param Return_Bootstrap_Values A logic vector of length 1 indicating whether or not the function should return the Control Efficacies, Test Efficacies, and, if Calculate_Hazard_Ratio = TRUE, Hazard Ratios (HRs) simulated in the semi-parametric bootstrap used to estimate the uncertainties of those values. Defaults to FALSE. This parameter can only be set to TRUE if Calculate_Uncertainty = TRUE.
#'
#'@details
#'Uncertainty estimates for values calculated by this function are generated using a semi-parametric bootstrap approach. This is performed in several steps.\enumerate{
#'\item Control efficacies for each compound/concentration are simulated by random sampling from normal distributions with means equal to the provided calculated efficacies and standard deviations equal to the provided efficacy standard errors.
#'\item Test efficacies are simulated in the same fashion as the control efficacies, except in cases when a test drug also exists in the control therapy. In such cases, it is assumed that the efficacy values for this drug are derived from the same dose-response curve for the both the control and test therapies, so each simulated efficacy for that test drug is matched to the corresponding simulated efficacy from that drug in the control therapy using a standard normal deviate.
#'\item Efficacy predictions are made for the control and test therapies for each cell line and set of simulated efficacies using the assumptions of independent drug action.
#'\item Cell lines are randomly sampled with replacement for each simulation as many times as there are original cell lines. The simulated control and test therapy efficacies are then sampled according to the sampled cell lines for each simulation.
#'\item Mean efficacies are calculated for the control and test therapies for each simulation. If specified to do so, these values are then used to calculate simulated HRs.
#'\item The simulated distributions of each efficacy metric are used to estimate uncertainties for those metrics.
#'}
#'
#'@return \itemize{
#'\item If Return_Bootstrap_Values = FALSE, this function returns a list with 4 elements: 1) Either a data frame with the calculated efficacy predictions, or, if an error occurred, a character vector of length one with the error message. 2) A data frame listing the control treatment drug names and concentrations. 3) A data frame listing the test treatment drug names and concentrations. 4) A character vector containing the names of the cell lines used to make the efficacy predictions.
#'\item If Return_Bootstrap_Values = TRUE & Calculate_Uncertainty = TRUE & Calculate_Hazard_Ratio = FALSE, this function returns a list with 6 elements: the first 4 elements are the same as when Return_Bootstrap_Values = FALSE and the fifth and sixth elements being numeric vectors of, respectively, the control and test viabilities simulated during the semi-parametric bootstrap used to estimate uncertainties.
#'\item If Return_Bootstrap_Values = TRUE & Calculate_Uncertainty = TRUE & Calculate_Hazard_Ratio = TRUE, this function returns a list with 7 elements: the first 4 elements are the same as when Return_Bootstrap_Values = FALSE and the fifth, sixth, and seventh elements being numeric vectors of, respectively, the control viabilities, test viabilities, and HRs simulated during the semi-parametric bootstrap used to estimate uncertainties.
#'}
#'
#' @examples
#' #Loading Package
#' library(IDACombo)
#'
#' #Making fake monotherapy dataset
#' CellLineNames <- rep(c("CL1", "CL2", "CL3", "CL4", "CL5", "CL6"), 6)
#' DrugNames <- c(rep("D1", 12), rep("D2", 12), rep("D3", 12))
#' Concentrations <- c(rep(1, 6), rep(2, 6), rep(1.5, 6), rep(3, 6), rep("A", 6), rep("B", 6))
#' Viability <- c(sample(seq(0.4,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.2,0.8,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.4,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.2,0.6,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.9,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.2,0.6,length.out = 10), 6, replace = TRUE))
#' Viability_SE <- Viability * sample(seq(0,0.1,length.out = 100), 36, replace = TRUE)
#' Fake_Data <- data.frame(CellLineNames, DrugNames, Concentrations, Viability, Viability_SE)
#'
#' #Creating efficacy predictions for control and test treatments and comparing without
#' #uncertainty calculations
#' #For case where drugs in test treatment are at reduced concentrations from
#' #those used in the control treatment due to the addition of a third drug.
#' #Note that this may mean that the test treatment is less effective than
#' #the control treatment, such that the Hazard Ratio is > 1.
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(1, 1.5, "B"),
#' Calculate_Uncertainty = FALSE,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE)
#'
#' #For case where drugs in test treatment are at same concentrations as
#' #those used in the control treatment.
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = FALSE,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Creating efficacy predictions for control and test treatments and comparing with
#' #uncertainty calculations but without returning simulated values that are generated
#' #using a semi-parametric bootstrap to estimate uncertainties.
#'
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Viability_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Creating efficacy predictions for control and test treatments and comparing with
#' #uncertainty calculations and with returning simulated values that are generated
#' #using a semi-parametric bootstrap to estimate standard errors.
#'
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Viability_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE,
#' Return_Bootstrap_Values = TRUE)
#'
#' #Converting Viabilty to reduction in viability and redoing calculations
#' #without returning simulated values from semi-parametric boostrap. Note the change
#' #in the LowerEfficacyIsBetterDrugEffect flag from TRUE to FALSE
#' Reduction_in_Viability <- 1-Viability
#' Reduction_in_Viability_SE <- Viability_SE
#' Fake_Data <- data.frame(CellLineNames,
#' DrugNames,
#' Concentrations,
#' Reduction_in_Viability,
#' Reduction_in_Viability_SE)
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Reduction_in_Viability",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Reduction_in_Viability_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = FALSE,
#' Efficacy_Metric_Name = "Reduction_in_Viability",
#' Calculate_Hazard_Ratio = TRUE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Changing efficacy metric to percent growth (range -1 to 1)
#' #Note that calculating Hazard Ratios is no longer valid, so
#' #Calculate_Hazard_Ratio is set to FALSE.
#' Percent_Growth <- c(sample(seq(0.4,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(-0.4,0.2,length.out = 10), 6, replace = TRUE),
#' sample(seq(-0.2,0.3,length.out = 10), 6, replace = TRUE),
#' sample(seq(-1,0.2,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.8,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(-1,-0.2,length.out = 10), 6, replace = TRUE))
#' Percent_Growth_SE <- abs(Percent_Growth * sample(seq(0,0.1,length.out = 100), 36, replace = TRUE))
#' Fake_Data <- data.frame(CellLineNames,
#' DrugNames,
#' Concentrations,
#' Percent_Growth,
#' Percent_Growth_SE)
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Percent_Growth",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Percent_Growth_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Percent_Growth",
#' Calculate_Hazard_Ratio = FALSE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Adding duplicate records for each cell line, and showing behavior with
#' #Average_Duplicate_Records = FALSE. Should produce warning messages that
#' #duplicates were found and removed.
#' Percent_Growth <- c(sample(seq(0.4,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(-0.4,0.2,length.out = 10), 6, replace = TRUE),
#' sample(seq(-0.2,0.3,length.out = 10), 6, replace = TRUE),
#' sample(seq(-1,0.2,length.out = 10), 6, replace = TRUE),
#' sample(seq(0.8,1,length.out = 10), 6, replace = TRUE),
#' sample(seq(-1,-0.2,length.out = 10), 6, replace = TRUE))
#' Percent_Growth_SE <- abs(Percent_Growth * sample(seq(0,0.1,length.out = 100), 36, replace = TRUE))
#' Fake_Data_to_add <- data.frame(CellLineNames,
#' DrugNames,
#' Concentrations,
#' Percent_Growth,
#' Percent_Growth_SE)
#' Fake_Data <- rbind(Fake_Data, Fake_Data_to_add)
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Percent_Growth",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Percent_Growth_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Percent_Growth",
#' Calculate_Hazard_Ratio = FALSE,
#' Average_Duplicate_Records = FALSE)
#'
#' #Now setting to average duplicate values.
#' IDAPredict.TestvsControl(Monotherapy_Data = Fake_Data,
#' Cell_Line_Name_Column = "CellLineNames",
#' Drug_Name_Column = "DrugNames",
#' Drug_Concentration_Column = "Concentrations",
#' Efficacy_Column = "Percent_Growth",
#' Control_Treatment_Drugs = c("D1", "D2"),
#' Control_Treatment_Drug_Concentrations = c(2, 3),
#' Test_Treatment_Drugs = c("D1", "D2", "D3"),
#' Test_Treatment_Drug_Concentrations = c(2, 3, "B"),
#' Calculate_Uncertainty = TRUE,
#' Efficacy_SE_Column = "Percent_Growth_SE",
#' n_Simulations = 1000,
#' LowerEfficacyIsBetterDrugEffect = TRUE,
#' Efficacy_Metric_Name = "Percent_Growth",
#' Calculate_Hazard_Ratio = FALSE,
#' Average_Duplicate_Records = TRUE)
#'
#' @export
IDAPredict.TestvsControl <- function(Monotherapy_Data, Cell_Line_Name_Column, Drug_Name_Column, Drug_Concentration_Column, Efficacy_Column, LowerEfficacyIsBetterDrugEffect, Efficacy_Metric_Name = "Efficacy", Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations, Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations, Calculate_Uncertainty = FALSE, Efficacy_SE_Column = NULL, n_Simulations = 1000, Calculate_Hazard_Ratio = FALSE, Average_Duplicate_Records = FALSE, Return_Bootstrap_Values = FALSE){
#Checking that all input variables are in correct format
if(! is.data.frame(Monotherapy_Data)){
stop("Monotherapy_Data is not a data frame.")
}
if(! is.vector(Cell_Line_Name_Column) | ! is.character(Cell_Line_Name_Column) | ! length(Cell_Line_Name_Column) == 1){
stop("Cell_Line_Name_Column is not a character vector of length 1.")
}
if(! Cell_Line_Name_Column %in% colnames(Monotherapy_Data)){
stop("Cell_Line_Name_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(Drug_Name_Column) | ! is.character(Drug_Name_Column) | ! length(Drug_Name_Column) == 1){
stop("Drug_Name_Column is not a character vector of length 1.")
}
if(! Drug_Name_Column %in% colnames(Monotherapy_Data)){
stop("Drug_Name_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(Drug_Concentration_Column) | ! is.character(Drug_Concentration_Column) | ! length(Drug_Concentration_Column) == 1){
stop("Drug_Concentration_Column is not a character vector of length 1.")
}
if(! Drug_Concentration_Column %in% colnames(Monotherapy_Data)){
stop("Drug_Concentration_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(Efficacy_Column) | ! is.character(Efficacy_Column) | ! length(Efficacy_Column) == 1){
stop("Efficacy_Column is not a character vector of length 1.")
}
if(! Efficacy_Column %in% colnames(Monotherapy_Data)){
stop("Efficacy_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(LowerEfficacyIsBetterDrugEffect) | ! is.logical(LowerEfficacyIsBetterDrugEffect) | ! length(LowerEfficacyIsBetterDrugEffect) == 1){
stop("LowerEfficacyIsBetterDrugEffect is not a logical vector of length 1.")
}
if(! is.vector(Efficacy_Metric_Name) | ! is.character(Efficacy_Metric_Name) | ! length(Efficacy_Metric_Name) == 1){
stop("Efficacy_Metric_Name is not a character vector of length 1.")
}
if(! is.vector(Control_Treatment_Drugs) | ! is.character(Control_Treatment_Drugs) | ! length(Control_Treatment_Drugs) > 0){
stop("Control_Treatment_Drugs is not a character vector with length > 0.")
}
if(! all(Control_Treatment_Drugs %in% Monotherapy_Data[,Drug_Name_Column])){
missing.control.drugs <- Control_Treatment_Drugs[! Control_Treatment_Drugs %in% Monotherapy_Data[,Drug_Name_Column]]
stop(paste0("No data for the following Control_Treatment_Drugs found in Monotherapy_Data: ", paste(missing.control.drugs, collapse = ", ")))
}
if(! is.vector(Control_Treatment_Drug_Concentrations) | ! length(Control_Treatment_Drug_Concentrations) == length(Control_Treatment_Drugs)){
stop("Control_Treatment_Drug_Concentrations is not a vector with length = length(Control_Treatment_Drugs).")
}
if(! is.vector(Test_Treatment_Drugs) | ! is.character(Test_Treatment_Drugs) | ! length(Test_Treatment_Drugs) > 0){
stop("Test_Treatment_Drugs is not a character vector with length > 0.")
}
if(! all(Test_Treatment_Drugs %in% Monotherapy_Data[,Drug_Name_Column])){
missing.control.drugs <- Test_Treatment_Drugs[! Test_Treatment_Drugs %in% Monotherapy_Data[,Drug_Name_Column]]
stop(paste0("No data for the following Test_Treatment_Drugs found in Monotherapy_Data: ", paste(missing.control.drugs, collapse = ", ")))
}
if(! is.vector(Test_Treatment_Drug_Concentrations) | ! length(Test_Treatment_Drug_Concentrations) == length(Test_Treatment_Drugs)){
stop("Test_Treatment_Drug_Concentrations is not a vector with length = length(Test_Treatment_Drugs).")
}
if(! is.vector(Calculate_Uncertainty) | ! is.logical(Calculate_Uncertainty) | ! length(Calculate_Uncertainty) == 1){
stop("Calculate_Uncertainty is not a logical vector of length 1.")
}
if(Calculate_Uncertainty == TRUE){
if(is.null(Efficacy_SE_Column)){
stop("Calculate_Uncertainty is TRUE but Efficacy_SE_Column is not specified. Please either set Calculate_Uncertainty to FALSE or specify Efficacy_SE_Column.")
}
if(! is.vector(Efficacy_SE_Column) | ! is.character(Efficacy_SE_Column) | ! length(Efficacy_SE_Column) == 1){
stop("Calculate_Uncertainty is TRUE, but Efficacy_SE_Column is not a character vector of length 1.")
}
if(! Efficacy_SE_Column %in% colnames(Monotherapy_Data)){
stop("Calculate_Uncertainty is TRUE, but Efficacy_SE_Column does not match any column names in Monotherapy_Data.")
}
if(! is.vector(n_Simulations) | ! is.numeric(n_Simulations) | ! length(n_Simulations) == 1){
stop("Calculate_Uncertainty is TRUE, but n_Simulations is not a numeric vector of length 1.")
}
if(! n_Simulations%%1==0 | ! n_Simulations >= 40){
stop("Calculate_Uncertainty is TRUE, but n_Simulations is not a positive integer >= 40.")
}
}
if(! is.vector(Return_Bootstrap_Values) | ! is.logical(Return_Bootstrap_Values) | ! length(Return_Bootstrap_Values) == 1){
stop("Return_Bootstrap_Values is not a logical vector of length 1.")
}
if(Return_Bootstrap_Values == TRUE & ! Calculate_Uncertainty == TRUE){
stop("Return_Bootstrap_Values is TRUE, but Calculate_Uncertainty is not TRUE.")
}
if(! is.vector(Calculate_Hazard_Ratio) | ! is.logical(Calculate_Hazard_Ratio) | ! length(Calculate_Hazard_Ratio) == 1){
stop("Calculate_Hazard_Ratio is not a logical vector of length 1.")
}
if(! is.vector(Average_Duplicate_Records) | ! is.logical(Average_Duplicate_Records) | ! length(Average_Duplicate_Records) == 1){
stop("Average_Duplicate_Records is not a logical vector of length 1.")
}
#Organizing data into standard format based on column names provided for each desired set of information
#Also subsetting to only include data pertaining to drugs in control or test treatment
if(Calculate_Uncertainty == TRUE){
Data <- Monotherapy_Data[Monotherapy_Data[,Drug_Name_Column] %in% c(Control_Treatment_Drugs, Test_Treatment_Drugs),c(Cell_Line_Name_Column, Drug_Name_Column, Drug_Concentration_Column, Efficacy_Column, Efficacy_SE_Column)]
colnames(Data) <- c("CellLine", "Drug", "Conc", "Efficacy", "Efficacy_SE")
} else {
Data <- Monotherapy_Data[Monotherapy_Data[,Drug_Name_Column] %in% c(Control_Treatment_Drugs, Test_Treatment_Drugs),c(Cell_Line_Name_Column, Drug_Name_Column, Drug_Concentration_Column, Efficacy_Column)]
colnames(Data) <- c("CellLine", "Drug", "Conc", "Efficacy")
}
rm(Monotherapy_Data)
#Making sure all columns are in correct formats
Data$CellLine <- as.character(Data$CellLine)
Data$Drug <- as.character(Data$Drug)
Data$Conc <- as.character(Data$Conc)
Data$Efficacy <- as.numeric(as.character(Data$Efficacy))
if(Calculate_Uncertainty == TRUE){
Data$Efficacy_SE <- as.numeric(as.character(Data$Efficacy_SE))
}
#Removing rows that are missing information
#Note: missing SE information is ignored if Calculate_Uncertainty == FALSE
Data <- Data[complete.cases(Data),]
#Subsetting into control and test group data with specified concentrations
ControlData <- list(NULL)
for(i in 1:length(Control_Treatment_Drugs)){
ControlData[[i]] <- Data[Data$Drug %in% Control_Treatment_Drugs[i],]
#Checking that all provided concentrations are available for their respective drugs
if(! Control_Treatment_Drug_Concentrations[i] %in% ControlData[[i]]$Conc){
stop(paste0(Control_Treatment_Drug_Concentrations[i], " concentration is unavailable for ", Control_Treatment_Drugs[i], " in control treatment."))
} else {
ControlData[[i]] <- ControlData[[i]][ControlData[[i]]$Conc %in% Control_Treatment_Drug_Concentrations[i],]
}
}
names(ControlData) <- Control_Treatment_Drugs
TestData <- list(NULL)
for(i in 1:length(Test_Treatment_Drugs)){
TestData[[i]] <- Data[Data$Drug %in% Test_Treatment_Drugs[i],]
#Checking that all provided concentrations are available for their respective drugs
if(! Test_Treatment_Drug_Concentrations[i] %in% TestData[[i]]$Conc){
stop(paste0(Test_Treatment_Drug_Concentrations[i], " concentration is unavailable for ", Test_Treatment_Drugs[i], " in test treatment."))
} else {
TestData[[i]] <- TestData[[i]][TestData[[i]]$Conc %in% Test_Treatment_Drug_Concentrations[i],]
}
}
names(TestData) <- Test_Treatment_Drugs
rm(Data)
#Finding cell line overlap between all drugs
ControlCellLines <- list(NULL)
for(i in 1:length(ControlData)){
ControlCellLines[[i]] <- sort(unique(ControlData[[i]]$CellLine))
}
TestCellLines <- list(NULL)
for(i in 1:length(TestData)){
TestCellLines[[i]] <- sort(unique(TestData[[i]]$CellLine))
}
All_Drug_CellLines <- c(ControlCellLines, TestCellLines)
Usable_CellLines <- sort(unique(unlist(All_Drug_CellLines)))
for(i in 1:length(All_Drug_CellLines)){
Usable_CellLines <- Usable_CellLines[Usable_CellLines %in% All_Drug_CellLines[[i]]]
}
rm(ControlCellLines, TestCellLines, All_Drug_CellLines)
#Checking again if at least 2 cell lines remain for all drugs. If not, exiting with no
#predictions and a warning.
if(! length(Usable_CellLines) >= 2){
#Returning NA predictions with warning due to too few cell lines.
warning(paste0("<2 overlapping cell lines available for comparison of (", paste(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations, sep = "_", collapse = " + "), ") vs. (", paste(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations, sep = "_", collapse = " + "), ")"))
if(Return_Bootstrap_Values == FALSE){
Return_Object <- list("Less than 2 overlapping cell lines available.", as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used")
return(Return_Object)
} else if(Return_Bootstrap_Values == TRUE & Calculate_Hazard_Ratio == TRUE){
Return_Object <- list("Less than 2 overlapping cell lines available.", as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines, NULL, NULL, NULL)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used", "Bootstrap_Mean_Control_Efficacies", "Bootstrap_Mean_Test_Efficacies", "Bootstrap_HRs")
return(Return_Object)
} else if(Return_Bootstrap_Values == TRUE & Calculate_Hazard_Ratio == FALSE){
Return_Object <- list("Less than 2 overlapping cell lines available.", as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines, NULL, NULL)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used", "Bootstrap_Mean_Control_Efficacies", "Bootstrap_Mean_Test_Efficacies")
return(Return_Object)
}
}
#Subsetting drug data to only include overlapping cell lines
for(i in 1:length(ControlData)){
ControlData[[i]] <- ControlData[[i]][ControlData[[i]]$CellLine %in% Usable_CellLines,]
}
for(i in 1:length(TestData)){
TestData[[i]] <- TestData[[i]][TestData[[i]]$CellLine %in% Usable_CellLines,]
}
#Checking that cell lines aren't duplicated in each drug dataset
#If Average_Duplicate_Records == FALSE, removing cell line duplicates with warning if duplicates are found.
#If Average_Duplicate_Records == TRUE, averaging duplicate records without warning.
for(i in 1:length(Control_Treatment_Drugs)){
CL_Conc <- paste(ControlData[[i]]$CellLine, ControlData[[i]]$Conc, sep = "_")
Dups <- CL_Conc[duplicated(CL_Conc)]
if(length(Dups) > 0 & Average_Duplicate_Records == FALSE){
warning(paste0("Duplicated information found for the following cell lines and ", Control_Treatment_Drugs[i], " concentrations in the control treatment. Average_Duplicate_Records = FALSE so duplicates removed: ", paste(Dups, collapse = ", ")))
ControlData[[i]] <- ControlData[[i]][! duplicated(CL_Conc),]
} else if(length(Dups) > 0 & Average_Duplicate_Records == TRUE){
if(Calculate_Uncertainty == FALSE){
#Simply averaging efficacy values
colnames <- colnames(ControlData[[i]])
ControlData[[i]] <- aggregate(ControlData[[i]]$Efficacy, by = list(ControlData[[i]]$CellLine, ControlData[[i]]$Drug, ControlData[[i]]$Conc), FUN = mean)
colnames(ControlData[[i]]) <- colnames
} else if(Calculate_Uncertainty == TRUE){
#Averaging efficacy
Efficacy_Average <- aggregate(ControlData[[i]]$Efficacy, by = list(ControlData[[i]]$CellLine, ControlData[[i]]$Drug, ControlData[[i]]$Conc), FUN = mean)
colnames(Efficacy_Average) <- c("CellLine", "Drug", "Conc", "Efficacy")
#Calculating uncertainty in averaged efficacy by adding efficacy uncertainties in quadrature and dividing by number of values used in average
Efficacy_Average_Uncertainties <- aggregate(ControlData[[i]]$Efficacy_SE, by = list(ControlData[[i]]$CellLine, ControlData[[i]]$Drug, ControlData[[i]]$Conc), FUN = function(x){sqrt(sum(x^2))/length(x)})
colnames(Efficacy_Average_Uncertainties) <- c("CellLine", "Drug", "Conc", "Efficacy_SE")
#Combining results
ControlData[[i]] <- merge(Efficacy_Average, Efficacy_Average_Uncertainties)
}
}
}
for(i in 1:length(Test_Treatment_Drugs)){
CL_Conc <- paste(TestData[[i]]$CellLine, TestData[[i]]$Conc, sep = "_")
Dups <- CL_Conc[duplicated(CL_Conc)]
if(length(Dups) > 0 & Average_Duplicate_Records == FALSE){
warning(paste0("Duplicated information found for the following cell lines and ", Test_Treatment_Drugs[i], " concentrations in the test treatment. Average_Duplicate_Records = FALSE so duplicates removed: ", paste(Dups, collapse = ", ")))
TestData[[i]] <- TestData[[i]][! duplicated(CL_Conc),]
} else if(length(Dups) > 0 & Average_Duplicate_Records == TRUE){
if(Calculate_Uncertainty == FALSE){
#Simply averaging efficacy values
colnames <- colnames(TestData[[i]])
TestData[[i]] <- aggregate(TestData[[i]]$Efficacy, by = list(TestData[[i]]$CellLine, TestData[[i]]$Drug, TestData[[i]]$Conc), FUN = mean)
colnames(TestData[[i]]) <- colnames
} else if(Calculate_Uncertainty == TRUE){
#Averaging efficacy
Efficacy_Average <- aggregate(TestData[[i]]$Efficacy, by = list(TestData[[i]]$CellLine, TestData[[i]]$Drug, TestData[[i]]$Conc), FUN = mean)
colnames(Efficacy_Average) <- c("CellLine", "Drug", "Conc", "Efficacy")
#Calculating uncertainty in averaged efficacy by adding efficacy uncertainties in quadrature and dividing by number of values used in average
Efficacy_Average_Uncertainties <- aggregate(TestData[[i]]$Efficacy_SE, by = list(TestData[[i]]$CellLine, TestData[[i]]$Drug, TestData[[i]]$Conc), FUN = function(x){sqrt(sum(x^2))/length(x)})
colnames(Efficacy_Average_Uncertainties) <- c("CellLine", "Drug", "Conc", "Efficacy_SE")
#Combining results
TestData[[i]] <- merge(Efficacy_Average, Efficacy_Average_Uncertainties)
}
}
}
#Ordering cell lines the same for each drug
for(i in 1:length(ControlData)){
ControlData[[i]] <- ControlData[[i]][order(ControlData[[i]]$Conc, ControlData[[i]]$CellLine),]
}
for(i in 1:length(TestData)){
TestData[[i]] <- TestData[[i]][order(TestData[[i]]$Conc, TestData[[i]]$CellLine),]
}
#Creating object to store prediction results in
Prediction_Results <- as.data.frame(matrix(rep(NA, 11), nrow = 1))
colnames(Prediction_Results) <- c("Mean_Control_Treatment_Efficacy", "Mean_Control_Treatment_Efficacy_SE", "Mean_Control_Treatment_Efficacy_95%_Confidence_Interval", "Mean_Test_Treatment_Efficacy", "Mean_Test_Treatment_Efficacy_SE", "Mean_Test_Treatment_Efficacy_95%_Confidence_Interval", "HR_Test_vs_Control_Treatment", "HR_Test_vs_Control_Treatment_SE", "HR_Test_vs_Control_Treatment_95%_Confidence_Interval", "p_HR>=1", "p_HR=1")
#Performing combination efficacy predictions for cases where lower efficacy values
#indicate a more effective drug effect (i.e. efficacy = percent viability, percent growth, etc.)
if(LowerEfficacyIsBetterDrugEffect == TRUE){
#Calculating expected combination efficacy for each cell line
#using Independent Drug Action for both Control and Test treatment
Control_Efficacy <- do.call(pmin, lapply(ControlData, function(x){return(x$Efficacy)}))
Test_Efficacy <- do.call(pmin, lapply(TestData, function(x){return(x$Efficacy)}))
#Calculating average efficacy across all cell lines
Prediction_Results$Mean_Control_Treatment_Efficacy <- mean(Control_Efficacy)
Prediction_Results$Mean_Test_Treatment_Efficacy <- mean(Test_Efficacy)
if(Calculate_Hazard_Ratio == TRUE){
Prediction_Results$HR_Test_vs_Control_Treatment <- Prediction_Results$Mean_Test_Treatment_Efficacy / Prediction_Results$Mean_Control_Treatment_Efficacy
}
rm(Control_Efficacy, Test_Efficacy)
#If Calculate_Uncertainty == TRUE, doing Monte Carlo simulation to estimate
#uncertainties in output parameters based on uncertainties in monotherapy efficacies.
if(Calculate_Uncertainty == TRUE){
#Looping through each drug for each treatment and simulating efficacies based on
#measured efficacies and SE's
Control_MC_Efficacies <- as.list(NULL)
for(i in 1:length(ControlData)){
Control_MC_Efficacies[[i]] <- apply(ControlData[[i]][,c("Efficacy", "Efficacy_SE")], 1, function(x){rnorm(n = n_Simulations, mean = x[1], sd = x[2])})
colnames(Control_MC_Efficacies[[i]]) <- ControlData[[i]]$CellLine
}
names(Control_MC_Efficacies) <- names(ControlData)
Test_MC_Efficacies <- as.list(NULL)
for(i in 1:length(TestData)){
#Checking if this test drug is also in the control therapy.
#If so, not randomly sampling for test viabilities. Instead, calculating how many SE's from
#the measured value each simulated control viability fell, and matching that distance in
#the simulated test viabilities. This is done because the drug in the control and test
#therapies are not independent in this case--they will have come from the same dose-response curve.
if(names(TestData)[i] %in% names(ControlData)){
#Calculating the number of SEs away from the measured value each simulated value is in the control therapy for this drug
Measured_Control_Data <- ControlData[[names(TestData)[i]]][,c("CellLine", "Efficacy", "Efficacy_SE")]
Measured_Control_Data <- Measured_Control_Data[match(colnames(Control_MC_Efficacies[[names(TestData)[i]]]), Measured_Control_Data$CellLine),]
Measured_Control_Efficacies <- matrix(Measured_Control_Data$Efficacy, ncol = length(Measured_Control_Data$Efficacy), nrow = n_Simulations, byrow = TRUE)
Measured_Control_Efficacy_SEs <- matrix(Measured_Control_Data$Efficacy_SE, ncol = length(Measured_Control_Data$Efficacy_SE), nrow = n_Simulations, byrow = TRUE)
SEs_deviated <- (Control_MC_Efficacies[[names(TestData)[i]]] - Measured_Control_Efficacies) / Measured_Control_Efficacy_SEs
#Calculating simulated test viabilities based on the SE deviations from the simulated control viabilities for this drug
Measured_Test_Data <- TestData[[names(TestData)[i]]][,c("CellLine", "Efficacy", "Efficacy_SE")]
Measured_Test_Data <- Measured_Test_Data[match(colnames(SEs_deviated), Measured_Test_Data$CellLine),]
Measured_Test_Efficacies <- matrix(Measured_Test_Data$Efficacy, ncol = length(Measured_Test_Data$Efficacy), nrow = n_Simulations, byrow = TRUE, dimnames = list(NULL, Measured_Test_Data$CellLine))
Measured_Test_Efficacy_SEs <- matrix(Measured_Test_Data$Efficacy_SE, ncol = length(Measured_Test_Data$Efficacy_SE), nrow = n_Simulations, byrow = TRUE, dimnames = list(NULL, Measured_Test_Data$CellLine))
Test_MC_Efficacies[[i]] <- Measured_Test_Efficacies + (SEs_deviated * Measured_Test_Efficacy_SEs)
} else {
#If test drug is not in control therapy, randomly sampling
Test_MC_Efficacies[[i]] <- apply(TestData[[i]][,c("Efficacy", "Efficacy_SE")], 1, function(x){rnorm(n = n_Simulations, mean = x[1], sd = x[2])})
}
}
names(Test_MC_Efficacies) <- names(TestData)
if(exists("Measured_Control_Data")){
rm(Measured_Control_Data, Measured_Control_Efficacies, Measured_Control_Efficacy_SEs, Measured_Test_Data, Measured_Test_Efficacies, Measured_Test_Efficacy_SEs)
}
#Calculating expected combination efficacies for each cell line
#using Independent Drug Action for both Control and Test treatment
Control_MC_Combo_Efficacies <- do.call(pmin, Control_MC_Efficacies)
Test_MC_Combo_Efficacies <- do.call(pmin, Test_MC_Efficacies)
#Resampling cell lines for each simulation with replacement to account for the effect of random variation in cell line selection
Selected_CCLs_Per_Sim <- t(apply(Control_MC_Combo_Efficacies, 1, function(x){return(sample(1:length(x), replace = TRUE))}))
Control_MC_Combo_Efficacies <- t(mapply(function(x,y){return(x[y])}, x = split(Control_MC_Combo_Efficacies, row(Control_MC_Combo_Efficacies)), y = split(Selected_CCLs_Per_Sim, row(Selected_CCLs_Per_Sim))))
Test_MC_Combo_Efficacies <- t(mapply(function(x,y){return(x[y])}, x = split(Test_MC_Combo_Efficacies, row(Test_MC_Combo_Efficacies)), y = split(Selected_CCLs_Per_Sim, row(Selected_CCLs_Per_Sim))))
#Calculating average efficacy across all cell lines
Mean_Control_Efficacies_MC <- rowMeans(Control_MC_Combo_Efficacies)
Mean_Test_Efficacies_MC <- rowMeans(Test_MC_Combo_Efficacies)
#Calculating standard errors
Prediction_Results$Mean_Control_Treatment_Efficacy_SE <- sd(Mean_Control_Efficacies_MC)
Prediction_Results$Mean_Test_Treatment_Efficacy_SE <- sd(Mean_Test_Efficacies_MC)
#Calculating 95% confidence intervals for mean treatment efficacies
index_2.5 <- floor(0.025*n_Simulations)
index_97.5 <- ceiling(0.975*n_Simulations)
Prediction_Results$`Mean_Control_Treatment_Efficacy_95%_Confidence_Interval` <- paste(sort(Mean_Control_Efficacies_MC, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
Prediction_Results$`Mean_Test_Treatment_Efficacy_95%_Confidence_Interval` <- paste(sort(Mean_Test_Efficacies_MC, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
#Calculating HR if specified to do so
if(Calculate_Hazard_Ratio == TRUE){
MC_HRs <- Mean_Test_Efficacies_MC / Mean_Control_Efficacies_MC
Prediction_Results$HR_Test_vs_Control_Treatment_SE <- sd(MC_HRs)
#Calculating HR 95% confidence interval
Prediction_Results$`HR_Test_vs_Control_Treatment_95%_Confidence_Interval` <- paste(sort(MC_HRs, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
#Calculating HR p values
n_equal_to_1 <- sum(MC_HRs == 1)
n_less_than_1 <- sum(MC_HRs < 1)
n_greater_than_1 <- sum(MC_HRs > 1)
lower_tail <- (n_less_than_1 + 0.5*n_equal_to_1)/length(MC_HRs)
upper_tail <- (n_greater_than_1 + 0.5*n_equal_to_1)/length(MC_HRs)
#Calculating one-sided p-value with null hypothesis that HR >= 1
Prediction_Results$`p_HR>=1` <- sum(MC_HRs >= 1) / length(MC_HRs)
#If p-value is 0, setting as p < minimum p value that can be estimated using this many simulations
if(Prediction_Results$`p_HR>=1` == 0){Prediction_Results$`p_HR>=1` <- paste0("<", 1/n_Simulations)}
#Calculating two-sided p-value with null hypothesis that HR = 1
Prediction_Results$`p_HR=1` <- min(lower_tail, upper_tail)*2
#If p-value is 0, setting as p < minimum p value that can be estimated using this many simulations
if(Prediction_Results$`p_HR=1` == 0){Prediction_Results$`p_HR=1` <- paste0("<", 1/n_Simulations)}
}
#Cleaning up
rm(Control_MC_Combo_Efficacies, Test_MC_Combo_Efficacies, Selected_CCLs_Per_Sim, index_2.5, index_97.5)
if(Return_Bootstrap_Values == FALSE){
rm(Mean_Test_Efficacies_MC, Mean_Control_Efficacies_MC)
}
if(Calculate_Hazard_Ratio == TRUE){
rm(n_equal_to_1, n_less_than_1, n_greater_than_1, lower_tail, upper_tail)
if(Return_Bootstrap_Values == FALSE){
rm(MC_HRs)
}
}
}
}
#Performing combination efficacy predictions for cases where lower efficacy values
#indicate a less effective drug effect (i.e. efficacy = percent cell death, etc.)
if(LowerEfficacyIsBetterDrugEffect == FALSE){
#Calculating expected combination efficacy for each cell line
#using Independent Drug Action for both Control and Test treatment
Control_Efficacy <- do.call(pmax, lapply(ControlData, function(x){return(x$Efficacy)}))
Test_Efficacy <- do.call(pmax, lapply(TestData, function(x){return(x$Efficacy)}))
#Calculating average efficacy across all cell lines
Prediction_Results$Mean_Control_Treatment_Efficacy <- mean(Control_Efficacy)
Prediction_Results$Mean_Test_Treatment_Efficacy <- mean(Test_Efficacy)
if(Calculate_Hazard_Ratio == TRUE){
Prediction_Results$HR_Test_vs_Control_Treatment <- (1 - Prediction_Results$Mean_Test_Treatment_Efficacy) / (1 - Prediction_Results$Mean_Control_Treatment_Efficacy)
}
rm(Control_Efficacy, Test_Efficacy)
#If Calculate_Uncertainty == TRUE, doing Monte Carlo simulation to estimate
#uncertainties in output parameters based on uncertainties in monotherapy efficacies.
if(Calculate_Uncertainty == TRUE){
#Looping through each drug for each treatment and simulating efficacies based on
#measured efficacies and SE's
Control_MC_Efficacies <- as.list(NULL)
for(i in 1:length(ControlData)){
Control_MC_Efficacies[[i]] <- apply(ControlData[[i]][,c("Efficacy", "Efficacy_SE")], 1, function(x){rnorm(n = n_Simulations, mean = x[1], sd = x[2])})
colnames(Control_MC_Efficacies[[i]]) <- ControlData[[i]]$CellLine
}
names(Control_MC_Efficacies) <- names(ControlData)
Test_MC_Efficacies <- as.list(NULL)
for(i in 1:length(TestData)){
#Checking if this test drug is also in the control therapy.
#If so, not randomly sampling for test viabilities. Instead, calculating how many SE's from
#the measured value each simulated control viability fell, and matching that distance in
#the simulated test viabilities. This is done because the drug in the control and test
#therapies are not independent in this case--they will have come from the same dose-response curve.
if(names(TestData)[i] %in% names(ControlData)){
#Calculating the number of SEs away from the measured value each simulated value is in the control therapy for this drug
Measured_Control_Data <- ControlData[[names(TestData)[i]]][,c("CellLine", "Efficacy", "Efficacy_SE")]
Measured_Control_Data <- Measured_Control_Data[match(colnames(Control_MC_Efficacies[[names(TestData)[i]]]), Measured_Control_Data$CellLine),]
Measured_Control_Efficacies <- matrix(Measured_Control_Data$Efficacy, ncol = length(Measured_Control_Data$Efficacy), nrow = n_Simulations, byrow = TRUE)
Measured_Control_Efficacy_SEs <- matrix(Measured_Control_Data$Efficacy_SE, ncol = length(Measured_Control_Data$Efficacy_SE), nrow = n_Simulations, byrow = TRUE)
SEs_deviated <- (Control_MC_Efficacies[[names(TestData)[i]]] - Measured_Control_Efficacies) / Measured_Control_Efficacy_SEs
#Calculating simulated test viabilities based on the SE deviations from the simulated control viabilities for this drug
Measured_Test_Data <- TestData[[names(TestData)[i]]][,c("CellLine", "Efficacy", "Efficacy_SE")]
Measured_Test_Data <- Measured_Test_Data[match(colnames(SEs_deviated), Measured_Test_Data$CellLine),]
Measured_Test_Efficacies <- matrix(Measured_Test_Data$Efficacy, ncol = length(Measured_Test_Data$Efficacy), nrow = n_Simulations, byrow = TRUE, dimnames = list(NULL, Measured_Test_Data$CellLine))
Measured_Test_Efficacy_SEs <- matrix(Measured_Test_Data$Efficacy_SE, ncol = length(Measured_Test_Data$Efficacy_SE), nrow = n_Simulations, byrow = TRUE, dimnames = list(NULL, Measured_Test_Data$CellLine))
Test_MC_Efficacies[[i]] <- Measured_Test_Efficacies + (SEs_deviated * Measured_Test_Efficacy_SEs)
}else {
#If test drug is not in control therapy, randomly sampling
Test_MC_Efficacies[[i]] <- apply(TestData[[i]][,c("Efficacy", "Efficacy_SE")], 1, function(x){rnorm(n = n_Simulations, mean = x[1], sd = x[2])})
}
}
names(Test_MC_Efficacies) <- names(TestData)
if(exists("Measured_Control_Data")){
rm(Measured_Control_Data, Measured_Control_Efficacies, Measured_Control_Efficacy_SEs, Measured_Test_Data, Measured_Test_Efficacies, Measured_Test_Efficacy_SEs)
}
#Calculating expected combination efficacies for each cell line
#using Independent Drug Action for both Control and Test treatment
Control_MC_Combo_Efficacies <- do.call(pmax, Control_MC_Efficacies)
Test_MC_Combo_Efficacies <- do.call(pmax, Test_MC_Efficacies)
#Resampling cell lines for each simulation with replacement to account for the effect of random variation in cell line selection
Selected_CCLs_Per_Sim <- t(apply(Control_MC_Combo_Efficacies, 1, function(x){return(sample(1:length(x), replace = TRUE))}))
Control_MC_Combo_Efficacies <- t(mapply(function(x,y){return(x[y])}, x = split(Control_MC_Combo_Efficacies, row(Control_MC_Combo_Efficacies)), y = split(Selected_CCLs_Per_Sim, row(Selected_CCLs_Per_Sim))))
Test_MC_Combo_Efficacies <- t(mapply(function(x,y){return(x[y])}, x = split(Test_MC_Combo_Efficacies, row(Test_MC_Combo_Efficacies)), y = split(Selected_CCLs_Per_Sim, row(Selected_CCLs_Per_Sim))))
#Calculating average efficacy across all cell lines
Mean_Control_Efficacies_MC <- rowMeans(Control_MC_Combo_Efficacies)
Mean_Test_Efficacies_MC <- rowMeans(Test_MC_Combo_Efficacies)
#Calculating standard errors
Prediction_Results$Mean_Control_Treatment_Efficacy_SE <- sd(Mean_Control_Efficacies_MC)
Prediction_Results$Mean_Test_Treatment_Efficacy_SE <- sd(Mean_Test_Efficacies_MC)
#Calculating 95% confidence intervals for mean treatment efficacies
index_2.5 <- floor(0.025*n_Simulations)
index_97.5 <- ceiling(0.975*n_Simulations)
Prediction_Results$`Mean_Control_Treatment_Efficacy_95%_Confidence_Interval` <- paste(sort(Mean_Control_Efficacies_MC, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
Prediction_Results$`Mean_Test_Treatment_Efficacy_95%_Confidence_Interval` <- paste(sort(Mean_Test_Efficacies_MC, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
#Calculating HRs if specified to do so
if(Calculate_Hazard_Ratio == TRUE){
MC_HRs <- (1 - Mean_Test_Efficacies_MC) / (1 - Mean_Control_Efficacies_MC)
Prediction_Results$HR_Test_vs_Control_Treatment_SE <- sd(MC_HRs)
#Calculating HR 95% confidence interval
Prediction_Results$`HR_Test_vs_Control_Treatment_95%_Confidence_Interval` <- paste(sort(MC_HRs, decreasing = FALSE)[c(index_2.5, index_97.5)], collapse = "_")
#Calculating HR p values
n_equal_to_1 <- sum(MC_HRs == 1)
n_less_than_1 <- sum(MC_HRs < 1)
n_greater_than_1 <- sum(MC_HRs > 1)
lower_tail <- (n_less_than_1 + 0.5*n_equal_to_1)/length(MC_HRs)
upper_tail <- (n_greater_than_1 + 0.5*n_equal_to_1)/length(MC_HRs)
#Calculating one-sided p-value with null hypothesis that HR >= 1
Prediction_Results$`p_HR>=1` <- sum(MC_HRs >= 1) / length(MC_HRs)
#If p-value is 0, setting as p < minimum p value that can be estimated using this many simulations
if(Prediction_Results$`p_HR>=1` == 0){Prediction_Results$`p_HR>=1` <- paste0("<", 1/n_Simulations)}
#Calculating two-sided p-value with null hypothesis that HR = 1
Prediction_Results$`p_HR=1` <- min(lower_tail, upper_tail)*2
#If p-value is 0, setting as p < minimum p value that can be estimated using this many simulations
if(Prediction_Results$`p_HR=1` == 0){Prediction_Results$`p_HR=1` <- paste0("<", 1/n_Simulations)}
}
#Cleaning up
rm(Control_MC_Combo_Efficacies, Test_MC_Combo_Efficacies, Selected_CCLs_Per_Sim, index_2.5, index_97.5)
if(Return_Bootstrap_Values == FALSE){
rm(Mean_Test_Efficacies_MC, Mean_Control_Efficacies_MC)
}
if(Calculate_Hazard_Ratio == TRUE){
rm(n_equal_to_1, n_less_than_1, n_greater_than_1, lower_tail, upper_tail)
if(Return_Bootstrap_Values == FALSE){
rm(MC_HRs)
}
}
}
}
#Returning Outputs
#If Calculate_Uncertainty == FALSE, removing SE columns
if(Calculate_Uncertainty == FALSE){
Prediction_Results <- Prediction_Results[,-which(colnames(Prediction_Results) %in% c("Mean_Control_Treatment_Efficacy_SE", "Mean_Control_Treatment_Efficacy_95%_Confidence_Interval", "Mean_Test_Treatment_Efficacy_SE", "Mean_Test_Treatment_Efficacy_95%_Confidence_Interval", "HR_Test_vs_Control_Treatment_SE", "HR_Test_vs_Control_Treatment_95%_Confidence_Interval", "p_HR>=1", "p_HR=1"))]
}
#If Calculate_Hazard_Ratio == FALSE, removing HR columns
if(Calculate_Hazard_Ratio == FALSE){
Prediction_Results <- Prediction_Results[,-which(colnames(Prediction_Results) %in% c("HR_Test_vs_Control_Treatment", "HR_Test_vs_Control_Treatment_SE", "HR_Test_vs_Control_Treatment_95%_Confidence_Interval", "p_HR>=1", "p_HR=1"))]
}
#Replacing "Efficacy" with Efficacy_Metric_Name in column names of Dose_Comparisons
colnames(Prediction_Results) <- gsub("Efficacy", Efficacy_Metric_Name, colnames(Prediction_Results))
#Constructing Return_Object
if(Return_Bootstrap_Values == FALSE){
Return_Object <- list(Prediction_Results, as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used")
} else if(Return_Bootstrap_Values == TRUE & Calculate_Hazard_Ratio == FALSE & Calculate_Uncertainty == TRUE){
Return_Object <- list(Prediction_Results, as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines, Mean_Control_Efficacies_MC, Mean_Test_Efficacies_MC)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used", "Bootstrap_Mean_Control_Efficacies", "Bootstrap_Mean_Test_Efficacies")
} else if(Return_Bootstrap_Values == TRUE & Calculate_Hazard_Ratio == TRUE & Calculate_Uncertainty == TRUE){
Return_Object <- list(Prediction_Results, as.data.frame(cbind(Control_Treatment_Drugs, Control_Treatment_Drug_Concentrations), stringsAsFactors = F), as.data.frame(cbind(Test_Treatment_Drugs, Test_Treatment_Drug_Concentrations), stringsAsFactors = F), Usable_CellLines, Mean_Control_Efficacies_MC, Mean_Test_Efficacies_MC, MC_HRs)
names(Return_Object) <- c("Efficacy_Predictions", "Control_Treatment", "Test_Treatment", "Cell_Lines_Used", "Bootstrap_Mean_Control_Efficacies", "Bootstrap_Mean_Test_Efficacies", "Bootstrap_HRs")
}
#Returning output
return(Return_Object)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.