R/drexplorer_IAI.R
In lshen1/drexplorer2: Dose-response Explorer2
Defines functions
plotIAI
fitIAI
CI.delta
plot_median_effect
fit_median_efect
detect_ray_design
Documented in
CI.delta
detect_ray_design
fitIAI
fit_median_efect
plotIAI
plot_median_effect
#' Detect if a experiment is fixed-ray design and prepare input for CI computation
#'
#' @param d1 d1.
#' @param d2 d2.
#' @param e e.
#' @param tol tolerance in declaring fixed ratio.
#' @param d2.d1.force the ratio d2/d1 to be forced in fitting the linear model; this is only
#' effective if the data is not a fixed ratio design. Usually, the all-possible combination
#' design; specifying d2.d1.force will use the specified ratio (from available ratios) to fit
#' the model and generate dose-response curve.
#'
#' @return a list
#' @export
#' @examples
#' detect_ray_design(d1 = nl22B2$schd, d2 = nl22B2$hpr, e = nl22B2$y1)
detect_ray_design <- function(d1, d2, e, tol = 0.001, d2.d1.force = NA){
# Here d2.d1.force use ray design method to fix grid data. What is the rationale and how to justify it?
ind_eAbnorm <- e<0 | e>1
if(sum(ind_eAbnorm)>0) {
warning(sprintf('%d responses not in the range [0, 1]; truncated automatically', sum(ind_eAbnorm)))
e[e<0] <- 0
e[e>1] <- 1
}
ind1 <- d1!=0 & d2==0
ind2 <- d1==0 & d2!=0
ind12 <- d1!=0 & d2!=0 # where the dose differs
ratio_d <- d2[ind12]/d1[ind12] # pointwise ratio
isRayDesign <- all.equal(ratio_d, rep(mean(ratio_d), length(ratio_d)), tolerance = tol) # TRUE or a string
if(isRayDesign != TRUE) isRayDesign <- FALSE # make string as false
# d1, e1, d2, e2, d12, e12 added for CI.delta() estimation; d12 is just the sum; d2.d1 is the mean of ratio when not a ray design
if(isRayDesign != TRUE){
# not a ray design
if(is.na(d2.d1.force)){
d2.d1.force <- 1
}
if(!d2.d1.force %in% ratio_d){
stop(sprintf("This is not a fixed ratio design; the specified d2.d1 ratio is not supported in the data!"))
}
# force d2.d1 as the specified ratio or the default
d2.d1 <- d2.d1.force
} else {
d2.d1 <- mean(ratio_d) # this is not necessary for fixed ratio design
}
return(list(drMat = data.frame(d1 = d1, d2 = d2, e = e), isRayDesign = isRayDesign,
ind_eAbnorm = ind_eAbnorm, d1 = d1[ind1], e1 = e[ind1],
d2 = d2[ind2], e2 = e[ind2], d12 = d1[ind12]+d2[ind12], e12 = e[ind12],
ind12 = ind12, d2.d1 = d2.d1, d2.d1_avail = ratio_d, d2.d1.force = d2.d1.force))
}
#' Fit median-effect model for 2-drug combination with a fixed ratio desgin
#'
#' Notice that the median-effect model is a special case of the sigma Emax model;
#' 3 linear models are fitted; fixed ratio thus is required. All 3 linear models are fitted
#' using logit(E) ~ log(d1+d2) using a subset of the data depending on if the data is for durg 1,
#' drug 2 and the mixture
#'
#' @param d1 dose for drug 1.
#' @param d2 dose for drug 2.
#' @param e corresponding response in the range [0, 1].
#' @param name1 drug 1 name.
#' @param name2 drug 2 name.
#' @param d2.d1 d2/d1, the fixed ratio.
#' @param base base of logarithm; disabled since the dose-response curve estimation also
#' needs transformation.
#'
#' @return a list.
#' @export
fit_median_efect <- function(d1, d2, e,
name1 = 'Drug A', name2 = 'Drug B', d2.d1,
base = exp(1)){
dose1 <- d1; dose2 <- d2; fa <- e
fu <- 1-fa
## define an indicator for identifying doses satisfying dose2/dose1=d2.d1
ind_mix <- !(dose1+dose2==0 | fa<=0 | fa>=1) # exclude data outside range of the model
totdose <- rep(NA, length(dose1)) # match with input length,i.e. ind.ratio
totdose[ind_mix] <- (dose1+dose2)[ind_mix] # obtain a total dose to fit 3 linear models
####
resp <- rep(NA, length(fa)) # match with input length as in plot(medianEffect$logd,medianEffect$resp
ind_logit <- fa<=0 | fa>=1 # logit transform for good data; others are NA
resp[!ind_logit] <- log(fa[!ind_logit]/fu[!ind_logit])
logd <- rep(NA, length(totdose))
logd <- log(totdose)
#log10d <- log10(totdose) # log10(dose)
## define an indicator for identifying doses satisfying dose2/dose1=d2.d1
ind.ratio <- dose1!=0 & dose2!=0
ind2 <- abs(dose2-d2.d1*dose1) < 0.0001 # to also include 0 dose for dose-response curve without taking log
## Estimate the parameters using median-effect plot for two single drugs and
## their combination at the fixed ratio (dose of drug 2)/(dose of drug 1)=d2.d1.
lm1 <- lm(resp[dose2==0 & dose1!=0]~logd[dose2==0 & dose1!=0])
dm1 <- exp(-summary(lm1)$coef[1,1]/summary(lm1)$coef[2,1])
lm2 <- lm(resp[dose1==0 & dose2!=0]~logd[dose1==0 & dose2!=0])
dm2 <- exp(-summary(lm2)$coef[1,1]/summary(lm2)$coef[2,1])
lmcomb <- lm(resp[ind.ratio]~logd[ind.ratio])
dmcomb <- exp(-summary(lmcomb)$coef[1,1]/summary(lmcomb)$coef[2,1])
# some non-essential data are attached to facilitate visualization
return(list(lm1 = lm1, dm1 = dm1, lm2 = lm2, dm2 = dm2,
lmcomb = lmcomb, dmcomb = dmcomb, ind.ratio = ind.ratio,
ind2 = ind2, base = base, logd = logd, dose1 = dose1, dose2 = dose2,
resp = resp, totdose = totdose, fa = fa, name1 = name1, name2 = name2,
d2.d1 = d2.d1, d1 = d1, d2 = d2, e = e))
}
#' Generate plot for the fit from median-effect model
#'
#' @param medianEffect the fitted result by fit_median_efect().
#' @param type type of plot, options are from c('medianEffect', 'doseResponseCurve', 'contour').
#' @param contour.level levels (between 0 and 1 representing different response values) to draw
#' contour plot.
#' @param logd whether to plot logd in Dose-response curve; only effective if
#' type=='doseResponseCurve'.
#' @param align whether to align dose range. Imagine that 2 drugs are mixed in fixed ratio,
#' the doses in one drug will have low doses that the other drug and mixture do not have
#' which some biologists deams as a waste of graphic region. An easy fix is to truncate xlim.
#' This is only only effective if type=='doseResponseCurve'.
#' @param ylim limits on y axis. Default is set between 0 and 1.
#' @param cex.legend cex for the legend text'.
#' @param legend.position legend position passed to legend() function.
#'
#' @importFrom graphics contour title
plot_median_effect <- function(medianEffect,
type = c('medianEffect', 'doseResponseCurve', 'contour'),
contour.level = (1:9)/10, logd = FALSE, cex.legend = 1,
align = TRUE, ylim = c(0, 1), legend.position = 'topright'){
if(type == 'medianEffect'){
plot(medianEffect$logd,medianEffect$resp,type='n',
xlab = sprintf('Log(Dose)'), ylab = 'Log(E/(1-E))',
main = sprintf('Median Effect Plot \n-- d2.d1=%.1g', medianEffect$d2.d1))
abline(medianEffect$lm1, lty = 3)
with(medianEffect, points(logd[dose2==0 & dose1!=0], resp[dose2==0 & dose1!=0], pch = 1))
abline(medianEffect$lm2, lty = 2)
with(medianEffect, points(logd[dose1==0 & dose2!=0], resp[dose1==0 & dose2!=0], pch = 2))
abline(medianEffect$lmcomb, lty = 1)
with(medianEffect, points(logd[ind.ratio],resp[ind.ratio], pch = 16))
with(medianEffect, legend('topright', c(name1, name2, 'Mixture'),
pch = c(1, 2, 16), lty = c(3, 2, 1), bty = 'n'))
} else if(type == 'doseResponseCurve'){
if(logd){
dose <- noNA(format_grid(dose1 = with(medianEffect, c(totdose, dose1, dose2)),
n = 300)$xGrid)
if(align){
mindose <- min(log(dose[dose != 0]), na.rm = T)
} else {
dd <- c(dose, with(medianEffect, c(dose1, dose2)))
mindose <- min(log(dd[dd != 0]), na.rm = T)
}
# notice: when plot log dose, at low doses, the predicted response might not be present in medianEffect due to grid selection.
with(medianEffect, plot(log(totdose[ind2]), fa[ind2],
type = 'n', xlab = 'Log(Dose)', ylab = 'Relative viability',
ylim = ylim, xlim = c(mindose, max(log(dose)))))
with(medianEffect, lines(log(dose), with(medianEffect,
1.0/(1.0+(dm1/dose)^(summary(lm1)$coef[2,1]))),
type = "l", lty = 3))
with(medianEffect, points(log(dose1[dose2==0]),fa[dose2==0], pch = 1))
with(medianEffect, lines(log(dose), 1.0/(1.0+(dm2/dose)^(summary(lm2)$coef[2,1])),
type = "l", lty = 2))
with(medianEffect, points(log(dose2[dose1==0]),fa[dose1==0], pch=2))
with(medianEffect, lines(log(dose), 1.0/(1.0+(dmcomb/dose)^(summary(lmcomb)$coef[2,1])),
type = "l", lty = 1))
with(medianEffect, points(log(totdose[ind.ratio]),fa[ind.ratio], pch = 16))
title( "Dose-Response Curves" )
with(medianEffect, legend(legend.position, c(name1, name2, 'Mixture'),
pch = c(1, 2, 16), lty = c(3, 2, 1), bty = 'n', cex = cex.legend))
} else {
dose <- with(medianEffect, seq(0,max(totdose[ind2], na.rm = TRUE),
max(totdose[ind2], na.rm = TRUE)/100))
with(medianEffect, plot(totdose[ind2],fa[ind2], type = 'n',
xlab = 'Dose',ylab = 'Relative viability',
ylim = ylim, xlim = c(0, max(dose))))
with(medianEffect, lines(dose, with(medianEffect,
1.0/(1.0+(dm1/dose)^(summary(lm1)$coef[2,1]))),
type = "l", lty = 3))
with(medianEffect, points(dose1[dose2==0], fa[dose2==0], pch = 1))
with(medianEffect, lines(dose, 1.0/(1.0+(dm2/dose)^(summary(lm2)$coef[2,1])),
type = "l", lty = 2))
with(medianEffect, points(dose2[dose1==0], fa[dose1==0], pch = 2))
with(medianEffect, lines(dose, 1.0/(1.0+(dmcomb/dose)^(summary(lmcomb)$coef[2,1])),
type = "l", lty = 1))
with(medianEffect, points(totdose[ind.ratio],fa[ind.ratio], pch = 16))
title( "Dose-Response Curves" )
with(medianEffect, legend(legend.position, c(name1, name2, 'Mixture'),
pch = c(1, 2, 16), lty = c(3, 2, 1), bty = 'n', cex = cex.legend))
}
} else if(type == 'contour'){
## The contour plot of raw data.
temp1 <- sort(unique(medianEffect$dose1))
temp2 <- sort(unique(medianEffect$dose2))
isAllPossible <- length(temp1)*length(temp2)==length(medianEffect$fa) # all possible combination design
if(isAllPossible){
obsfa <- matrix(medianEffect$fa, length(temp1), length(temp2), byrow = T)
contour(temp1, temp2, obsfa, levels = contour.level,
xlab = sprintf('%s dose', medianEffect$name1),
ylab = sprintf('%s dose', medianEffect$name2))
title("Contour Plot")
} else {
ww <- sprintf('%s has %d doses; %s has %d doses; only %d total responses; This is not a all-possible design and thus no contour plot can be constructed!',
medianEffect$name1, length(temp1),
medianEffect$name2, length(temp2), length(medianEffect$fa))
warning(ww)
}
}
}
#' Estimate the interaction index as well as its confidence interval using delta method for
#' fixed ratio design
#'
#' This code is extracted from the source code distributed at
#' https://biostatistics.mdanderson.org/SoftwareDownload/
#'
#' Two versions are included in the source code: CI_IIV2 2008.SSC and CI_IIV2.SSC. At first
#' we implement CI_IIV2.SSC; this gives wider CI band; so we decide to use CI_IIV2 2008.SSC.
#' Further, CI_IIV2 2008.SSC has a more recent date and thus most updated.
#'
#' @param d1 dose for drug 1 where drug 2 has dose equal to 0.
#' @param e1 corresponding scaled response for drug 1, must between 0 and 1.
#' @param d2 dose for drug 2 where drug 1 has dose equal to 0.
#' @param e2 corresponding scaled response for drug 2, must between 0 and 1.
#' @param d12 the sum of doses from drug 1 and drug 2 when both drugs are administered.
#' @param e12 corresponding scaled response when both drugs are administered.
#' @param d2.d1 the ratio of the ray design. This is the fixed ratio of dose 2 divided by dose 1.
#' @param E a vector of responses (between 0 and 1) where IAI and confidence interval are to be
#' computed from.
#' @param alpha significance level of confidence interval.
#' @param min response values close to 0 and 1 will lead to extreme values when logit
#' transformation is applied. min is used to truncate the response value so that values less than
#' min are truncated at min.
#' @param max response values close to 0 and 1 will lead to extreme values when logit
#' transformation is applied. max is used to truncate the response value so that values larger than
#' max are truncated at max.
#'
#' @return a data frame with columns IAI, IAI.low, IAI.up, E, dx1 (corresponding dose of drug 1),
#' dx2 (corresponding dose of drug 1), dx12 (corresponding dose of combined drug, same as
#' definition of d12).
#' @references Lee, J. J., & Kong, M. (2009). Confidence intervals of interaction index for
#' assessing multiple drug interaction. Statistics in biopharmaceutical research, 1(1), 4-17.
#' @importFrom stats qt
CI.delta <- function(d1, e1, d2, e2, d12, e12, d2.d1, E,
alpha = 0.05, min = 0.02, max = 0.98){
# min and max are used to truncate the effect so that logit transform would not be obsurd; the original code is not very robust!
e1 <- truncate_effect(e1, min = min, max = max)
e2 <- truncate_effect(e2, min = min, max = max)
e12 <- truncate_effect(e12, min = min, max = max)
# if e1, e2 or e12 has value of 0 or 1, the logit would be Inf and leads to error for lm.
# thus use safe log
lm1 <- lm(log(e1/(1-e1))~log(d1))
dm1 <- exp(-summary(lm1)$coef[1,1]/summary(lm1)$coef[2,1])
lm2 <- lm(log(e2/(1-e2))~log(d2))
dm2 <- exp(-summary(lm2)$coef[1,1]/summary(lm2)$coef[2,1])
lmcomb <- lm(log(e12/(1-e12))~log (d12))
dm12 <- exp(-summary(lmcomb)$coef[1,1]/summary(lmcomb)$coef[2,1])
Dx1 <- dm1*(E/(1-E))^(1/summary(lm1)$coef[2,1])
Dx2 <- dm2*(E/(1-E))^(1/summary(lm2)$coef[2,1])
dx12 <- dm12*(E/(1-E))^(1/summary(lmcomb)$coef[2,1])
iix <- (dx12/(1+d2.d1))/Dx1+(dx12*d2.d1/(1+d2.d1))/Dx2
lm1.s <-summary(lm1)
lm2.s <-summary(lm2)
lm12.s <-summary(lmcomb)
c1 <- 1.0/lm1.s$coef[2,1]^2*lm1.s$coef[1,2]^2
temp <- - mean(log(d1))*lm1.s$coef[2,2]^2
c1 <- c1+2.0*(log(E/(1-E))-lm1.s$coef[1,1])/lm1.s$coef[2,1]^3*temp
c1 <- c1+(log(E/(1-E))-lm1.s$coef[1,1])^2/lm1.s$coef[2,1]^4*lm1.s$coef[2,2]^2
c2 <- 1.0/lm2.s$coef[2,1]^2*lm2.s$coef[1,2]^2
temp <- - mean(log(d2))*lm2.s$coef[2,2]^2
c2 <- c2+2.0*(log(E/(1-E))-lm2.s$coef[1,1])/lm2.s$coef[2,1]^3*temp
c2 <- c2+(log(E/(1-E))-lm2.s$coef[1,1])^2/lm2.s$coef[2,1]^4*lm2.s$coef[2,2]^2
c12 <- 1.0/lm12.s$coef[2,1]^2*lm12.s$coef[1,2]^2
temp <- - mean(log(d12))*lm12.s$coef[2,2]^2
c12 <- c12+2.0*(log(E/(1-E))-lm12.s$coef[1,1])/lm12.s$coef[2,1]^3*temp
c12 <- c12+(log(E/(1-E))-lm12.s$coef[1,1])^2/lm12.s$coef[2,1]^4*lm12.s$coef[2,2]^2
var.ii <-((dx12/Dx1)^2*c1+(dx12*d2.d1/Dx2)^2*c2+(1.0/Dx1+d2.d1/Dx2)^2*dx12^2*c12)/(1+d2.d1)^2
t975 <- qt(1-alpha/2,length(d1)+length(d2)+length(d12)-6)
iix.low1 <- iix*exp(-t975*var.ii^0.5/iix)
iix.up1 <- iix*exp(t975*var.ii^0.5/iix)
# add the corresponding dose for visualization
IAI <- data.frame(E = E, IAI = iix, IAI.low = iix.low1, IAI.up = iix.up1,
dx1 = Dx1, dx2 = Dx2, dx12 = dx12)
return(IAI)
}
#' Estimate the interaction index as well as its confidence interval using delta method
#' for fixed ratio design
#'
#' This code is extracted from the source code distributed at
#' https://biostatistics.mdanderson.org/SoftwareDownload/
#'
#' Two papers have been published by Lee et al, one in 2007 (Lee2007) and on in 2009 (Lee2009).
#' The Lee2007 paper described five methods to assess interaction: (1) Lowewe additivity model
#' using interaction index (IAI) (2) Model of Greco et al 1990. This approach uses
#' \deqn{\alpha}{alpha} as the metric and it can be related to IAI (3) Model of Machado and
#' Robinson which uses a metric denoted as \deqn{\eta}{eta} (4) Model of Plummer and Short which
#' can also be linked to IAI through the parameter \deqn{\beta_4}{beta_4} (5) Model of Carter
#' et al that can be linked to IAI through the parameter \deqn{\beta_{12}}{beta_12}. For more
#' details of these models, please refer to Lee2007.
#'
#' The Lee2009 paper provided generalization of IAI to multiple drugs using Lowewe additivity
#' model and assumption of Chou and Talalay's median effect equation. The Chou and Talalay's
#' median effect equation can be expressed as:
#' \deqn{log(E/(1-E))=m(log d - log Dm)}{log(E/(1-E))=m(log d - log Dm)}
#' where E is the effect at dose d for a compound whose median effective dose.
#'
#' Some notes about experiment design. Usually the data is either fixed ratio design (ray design)
#' or grid design which means all-possible combination of drug concentrations between two drugs
#' are available. The Lee2007 paper provided an example of grid design. However, specific fixed
#' ratio is used to fit the median effect model which is the basis to estimate IAI. The Lee2009
#' paper considered with fixed ratio design.
#'
#' @param d1 dose for drug 1.
#' @param d2 dose for drug 2.
#' @param e corresponding response in the range [0, 1].
#' @param E a vector of responses (between 0 and 1) where IAI and confidence interval are to be
#' computed from.
#' @param name1 name of drug 1.
#' @param name2 name of drug 2.
#' @param alpha significance level of confidence interval.
#' @param d2.d1.force a ratio passed to detect_ray_design() function so as to specify a fixed
#' ratio for grid design.
#' @param tol tolerance in declaring fixed ratio.
#'
#' @return a data frame with columns IAI, IAI.low, IAI.up, E, dx1 (corresponding dose of drug 1),
#' dx2 (corresponding dose of drug 1).
#'
#' @references Lee, J. J., & Kong, M. (2009). Confidence intervals of interaction index for
#' assessing multiple drug interaction. Statistics in biopharmaceutical research, 1(1), 4-17.
#' @references Lee, J. Jack, et al (2007). Interaction index and different methods for determining
#' drug interaction in combination therapy. Journal of biopharmaceutical statistics 17.3 461-480.
#'
#' @export
#' @examples
#' fit_allPoss_1 <- fitIAI(d1 = nl22B2$schd, d2 = nl22B2$hpr, e = nl22B2$y1,
#' name1 = 'SCH66336', name2 = '4HPR', d2.d1.force = 1)
fitIAI <- function(d1, d2, e, E = seq(0.05, 0.95, 0.005),
name1 = 'Drug A', name2 = 'Drug B', alpha = 0.05,
d2.d1.force = NA, tol = 0.01){
res_design <- detect_ray_design(d1 = d1, d2 = d2, e = e, d2.d1.force = d2.d1.force, tol = tol)
# fit median-effect model
medianEffect <- fit_median_efect(d1 = d1, d2 = d2, e = e, d2.d1 = res_design$d2.d1,
name1 = name1, name2 = name2)
if(res_design$isRayDesign != TRUE) { # not fixed ratio design, fitCI is NULL
warning("This is not a ray design; computation is done with method derived from fixed ratio!")
} else {
warning(sprintf("Fixed ratio design detected. Fixed ratio is %.3f\n",
res_design$d2.d1))
}
## fit CI
fitCI <- with(res_design,
CI.delta(d1 = d1, e1 = e1, d2 = d2, e2 = e2, d12 = d12,
e12 = e12, d2.d1, E = E, alpha = alpha))
meta <- list(d1 = d1, d2 = d2, e = e, E = E, alpha = alpha,
d2.d1 = res_design$d2.d1, isRayDesign = res_design$isRayDesign)
return(list(CI = fitCI, meta = meta, medianEffect = medianEffect))
}
#' Visualize interaction index
#'
#' @param fit the fitted result from fitIAI()
#' @param type type of plot from c('IAI', 'medianEffect', 'doseResponseCurve', 'contour')
#' @param contour.level contour level. only effective if type='contour'
#' @param ylim y axis limit
#' @param mode specify if to plot against response, dose or both. only effective if type=='IAI'. can be either 'response', 'dose', or 'both'
#' @param logd logd passed to plot_median_effect
#' @param align align passed to plot_median_effect
#' @param cex.legend cex.legend passed to plot_median_effect
#' @param legend.position legend.position passed to plot_median_effect
#'
#' @importFrom graphics mtext
#'
#' @export
plotIAI <- function(fit, type = c('IAI', 'medianEffect', 'doseResponseCurve', 'contour'),
contour.level = (1:9)/10, ylim = NULL, mode = 'both',
logd = FALSE, align = TRUE, cex.legend = 1, legend.position = 'topright'){
if(type == 'IAI' & !is.null(fit$CI)){
resCI <- fit$CI
if(is.null(ylim)) ylim <- c(0.01, 10)
op <- par(mar=c(8, 5, 4, 4) + 0.1)
lty_E <- 4
# CI
ind.CId.l <- resCI$IAI.low >= min(resCI$IAI)
ind.CId.u <- resCI$IAI.up <= max(resCI$IAI)+1
xlim_log10d <- range(log10(resCI$dx12))
col_d <- 'blue2'
lty_d <- 5
if(mode == 'both'){
# IAI ~ E
fa <- fit$meta$e
E <- resCI$E
with(resCI, plot(E, IAI, log = "y", type = "l", xlab = "", xaxt = 'n',
ylab = "Interaction Index",
xlim = c(min(E, fa),max(E, fa)), ylim = ylim,
cex = 0.6, main = 'Interaction plot'))
axis(1, at = pretty(range(E)), col = "black", lwd = 1)
mtext("Relative viability", side = 1, col = "black", line = 2)
abline(h = 1, lty = 4)
lines(E[ind.CId.l], resCI$IAI.low[ind.CId.l], lty = lty_E)
lines(E[ind.CId.u], resCI$IAI.up[ind.CId.u], lty = lty_E)
### add dose info
par(new = TRUE)
## IAI ~ dose
plot(resCI$dx12, resCI$IAI, log = "xy", type = "l", ylim = ylim,
axes = F, xlab = '', ylab = '', lty = lty_d, col = col_d)
axis(1, xlim = xlim_log10d, lwd = 1, line = 3.2,
col = col_d, col.axis = col_d, lty = lty_d)
mtext("Dose", side = 1, col = col_d, line = 5.2)
lines(resCI$dx12[ind.CId.l], resCI$IAI.low[ind.CId.l], lty = lty_d, col = col_d)
lines(resCI$dx12[ind.CId.u], resCI$IAI.up[ind.CId.u], lty = lty_d, col = col_d)
} else if(mode == 'response'){
# IAI ~ E
fa <- fit$meta$e
E <- resCI$E
plot(resCI$E, resCI$IAI, log = "y", type = "l", xlab = "", xaxt = 'n',
ylab = "Interaction Index", xlim = c(min(E, fa),max(E, fa)),
ylim = ylim,cex = 0.6, main = 'Interaction plot')
abline(h = 1, lty = 4)
axis(1, at = pretty(range(resCI$E)), col = "black", lwd = 1)
mtext("Relative viability", side = 1, col = "black", line = 2)
abline(h = 1, lty = 4)
lines(E[ind.CId.l], resCI$IAI.low[ind.CId.l], lty = lty_E)
lines(E[ind.CId.u], resCI$IAI.up[ind.CId.u], lty = lty_E)
} else if(mode == 'dose'){
## IAI ~ dose
plot(resCI$dx12, resCI$IAI, log = "xy", type = "l", xlab = "Dose",
ylab = "Interaction Index", ylim = ylim, cex = 0.6, main = 'Interaction plot')
abline(h = 1, lty = 4)
lines(resCI$dx12[ind.CId.l], resCI$IAI.low[ind.CId.l], lty = lty_d)
lines(resCI$dx12[ind.CId.u], resCI$IAI.up[ind.CId.u], lty = lty_d)
}
par(op)
} else {
medianEffect <- fit$medianEffect
if(type == 'contour'){
plot_median_effect(medianEffect, type = type, contour.level = contour.level)
} else {
plot_median_effect(medianEffect, type = type, logd = logd, align = align,
ylim = ylim, legend.position = legend.position, cex.legend = cex.legend)
}
}
}
lshen1/drexplorer2 documentation built on June 2, 2020, 9:27 p.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 plotIAI fitIAI CI.delta plot_median_effect fit_median_efect detect_ray_design
Documented in CI.delta detect_ray_design fitIAI fit_median_efect plotIAI plot_median_effect
#' Detect if a experiment is fixed-ray design and prepare input for CI computation
#'
#' @param d1 d1.
#' @param d2 d2.
#' @param e e.
#' @param tol tolerance in declaring fixed ratio.
#' @param d2.d1.force the ratio d2/d1 to be forced in fitting the linear model; this is only
#' effective if the data is not a fixed ratio design. Usually, the all-possible combination
#' design; specifying d2.d1.force will use the specified ratio (from available ratios) to fit
#' the model and generate dose-response curve.
#'
#' @return a list
#' @export
#' @examples
#' detect_ray_design(d1 = nl22B2$schd, d2 = nl22B2$hpr, e = nl22B2$y1)
detect_ray_design <- function(d1, d2, e, tol = 0.001, d2.d1.force = NA){
# Here d2.d1.force use ray design method to fix grid data. What is the rationale and how to justify it?
ind_eAbnorm <- e<0 | e>1
if(sum(ind_eAbnorm)>0) {
warning(sprintf('%d responses not in the range [0, 1]; truncated automatically', sum(ind_eAbnorm)))
e[e<0] <- 0
e[e>1] <- 1
}
ind1 <- d1!=0 & d2==0
ind2 <- d1==0 & d2!=0
ind12 <- d1!=0 & d2!=0 # where the dose differs
ratio_d <- d2[ind12]/d1[ind12] # pointwise ratio
isRayDesign <- all.equal(ratio_d, rep(mean(ratio_d), length(ratio_d)), tolerance = tol) # TRUE or a string
if(isRayDesign != TRUE) isRayDesign <- FALSE # make string as false
# d1, e1, d2, e2, d12, e12 added for CI.delta() estimation; d12 is just the sum; d2.d1 is the mean of ratio when not a ray design
if(isRayDesign != TRUE){
# not a ray design
if(is.na(d2.d1.force)){
d2.d1.force <- 1
}
if(!d2.d1.force %in% ratio_d){
stop(sprintf("This is not a fixed ratio design; the specified d2.d1 ratio is not supported in the data!"))
}
# force d2.d1 as the specified ratio or the default
d2.d1 <- d2.d1.force
} else {
d2.d1 <- mean(ratio_d) # this is not necessary for fixed ratio design
}
return(list(drMat = data.frame(d1 = d1, d2 = d2, e = e), isRayDesign = isRayDesign,
ind_eAbnorm = ind_eAbnorm, d1 = d1[ind1], e1 = e[ind1],
d2 = d2[ind2], e2 = e[ind2], d12 = d1[ind12]+d2[ind12], e12 = e[ind12],
ind12 = ind12, d2.d1 = d2.d1, d2.d1_avail = ratio_d, d2.d1.force = d2.d1.force))
}
#' Fit median-effect model for 2-drug combination with a fixed ratio desgin
#'
#' Notice that the median-effect model is a special case of the sigma Emax model;
#' 3 linear models are fitted; fixed ratio thus is required. All 3 linear models are fitted
#' using logit(E) ~ log(d1+d2) using a subset of the data depending on if the data is for durg 1,
#' drug 2 and the mixture
#'
#' @param d1 dose for drug 1.
#' @param d2 dose for drug 2.
#' @param e corresponding response in the range [0, 1].
#' @param name1 drug 1 name.
#' @param name2 drug 2 name.
#' @param d2.d1 d2/d1, the fixed ratio.
#' @param base base of logarithm; disabled since the dose-response curve estimation also
#' needs transformation.
#'
#' @return a list.
#' @export
fit_median_efect <- function(d1, d2, e,
name1 = 'Drug A', name2 = 'Drug B', d2.d1,
base = exp(1)){
dose1 <- d1; dose2 <- d2; fa <- e
fu <- 1-fa
## define an indicator for identifying doses satisfying dose2/dose1=d2.d1
ind_mix <- !(dose1+dose2==0 | fa<=0 | fa>=1) # exclude data outside range of the model
totdose <- rep(NA, length(dose1)) # match with input length,i.e. ind.ratio
totdose[ind_mix] <- (dose1+dose2)[ind_mix] # obtain a total dose to fit 3 linear models
####
resp <- rep(NA, length(fa)) # match with input length as in plot(medianEffect$logd,medianEffect$resp
ind_logit <- fa<=0 | fa>=1 # logit transform for good data; others are NA
resp[!ind_logit] <- log(fa[!ind_logit]/fu[!ind_logit])
logd <- rep(NA, length(totdose))
logd <- log(totdose)
#log10d <- log10(totdose) # log10(dose)
## define an indicator for identifying doses satisfying dose2/dose1=d2.d1
ind.ratio <- dose1!=0 & dose2!=0
ind2 <- abs(dose2-d2.d1*dose1) < 0.0001 # to also include 0 dose for dose-response curve without taking log
## Estimate the parameters using median-effect plot for two single drugs and
## their combination at the fixed ratio (dose of drug 2)/(dose of drug 1)=d2.d1.
lm1 <- lm(resp[dose2==0 & dose1!=0]~logd[dose2==0 & dose1!=0])
dm1 <- exp(-summary(lm1)$coef[1,1]/summary(lm1)$coef[2,1])
lm2 <- lm(resp[dose1==0 & dose2!=0]~logd[dose1==0 & dose2!=0])
dm2 <- exp(-summary(lm2)$coef[1,1]/summary(lm2)$coef[2,1])
lmcomb <- lm(resp[ind.ratio]~logd[ind.ratio])
dmcomb <- exp(-summary(lmcomb)$coef[1,1]/summary(lmcomb)$coef[2,1])
# some non-essential data are attached to facilitate visualization
return(list(lm1 = lm1, dm1 = dm1, lm2 = lm2, dm2 = dm2,
lmcomb = lmcomb, dmcomb = dmcomb, ind.ratio = ind.ratio,
ind2 = ind2, base = base, logd = logd, dose1 = dose1, dose2 = dose2,
resp = resp, totdose = totdose, fa = fa, name1 = name1, name2 = name2,
d2.d1 = d2.d1, d1 = d1, d2 = d2, e = e))
}
#' Generate plot for the fit from median-effect model
#'
#' @param medianEffect the fitted result by fit_median_efect().
#' @param type type of plot, options are from c('medianEffect', 'doseResponseCurve', 'contour').
#' @param contour.level levels (between 0 and 1 representing different response values) to draw
#' contour plot.
#' @param logd whether to plot logd in Dose-response curve; only effective if
#' type=='doseResponseCurve'.
#' @param align whether to align dose range. Imagine that 2 drugs are mixed in fixed ratio,
#' the doses in one drug will have low doses that the other drug and mixture do not have
#' which some biologists deams as a waste of graphic region. An easy fix is to truncate xlim.
#' This is only only effective if type=='doseResponseCurve'.
#' @param ylim limits on y axis. Default is set between 0 and 1.
#' @param cex.legend cex for the legend text'.
#' @param legend.position legend position passed to legend() function.
#'
#' @importFrom graphics contour title
plot_median_effect <- function(medianEffect,
type = c('medianEffect', 'doseResponseCurve', 'contour'),
contour.level = (1:9)/10, logd = FALSE, cex.legend = 1,
align = TRUE, ylim = c(0, 1), legend.position = 'topright'){
if(type == 'medianEffect'){
plot(medianEffect$logd,medianEffect$resp,type='n',
xlab = sprintf('Log(Dose)'), ylab = 'Log(E/(1-E))',
main = sprintf('Median Effect Plot \n-- d2.d1=%.1g', medianEffect$d2.d1))
abline(medianEffect$lm1, lty = 3)
with(medianEffect, points(logd[dose2==0 & dose1!=0], resp[dose2==0 & dose1!=0], pch = 1))
abline(medianEffect$lm2, lty = 2)
with(medianEffect, points(logd[dose1==0 & dose2!=0], resp[dose1==0 & dose2!=0], pch = 2))
abline(medianEffect$lmcomb, lty = 1)
with(medianEffect, points(logd[ind.ratio],resp[ind.ratio], pch = 16))
with(medianEffect, legend('topright', c(name1, name2, 'Mixture'),
pch = c(1, 2, 16), lty = c(3, 2, 1), bty = 'n'))
} else if(type == 'doseResponseCurve'){
if(logd){
dose <- noNA(format_grid(dose1 = with(medianEffect, c(totdose, dose1, dose2)),
n = 300)$xGrid)
if(align){
mindose <- min(log(dose[dose != 0]), na.rm = T)
} else {
dd <- c(dose, with(medianEffect, c(dose1, dose2)))
mindose <- min(log(dd[dd != 0]), na.rm = T)
}
# notice: when plot log dose, at low doses, the predicted response might not be present in medianEffect due to grid selection.
with(medianEffect, plot(log(totdose[ind2]), fa[ind2],
type = 'n', xlab = 'Log(Dose)', ylab = 'Relative viability',
ylim = ylim, xlim = c(mindose, max(log(dose)))))
with(medianEffect, lines(log(dose), with(medianEffect,
1.0/(1.0+(dm1/dose)^(summary(lm1)$coef[2,1]))),
type = "l", lty = 3))
with(medianEffect, points(log(dose1[dose2==0]),fa[dose2==0], pch = 1))
with(medianEffect, lines(log(dose), 1.0/(1.0+(dm2/dose)^(summary(lm2)$coef[2,1])),
type = "l", lty = 2))
with(medianEffect, points(log(dose2[dose1==0]),fa[dose1==0], pch=2))
with(medianEffect, lines(log(dose), 1.0/(1.0+(dmcomb/dose)^(summary(lmcomb)$coef[2,1])),
type = "l", lty = 1))
with(medianEffect, points(log(totdose[ind.ratio]),fa[ind.ratio], pch = 16))
title( "Dose-Response Curves" )
with(medianEffect, legend(legend.position, c(name1, name2, 'Mixture'),
pch = c(1, 2, 16), lty = c(3, 2, 1), bty = 'n', cex = cex.legend))
} else {
dose <- with(medianEffect, seq(0,max(totdose[ind2], na.rm = TRUE),
max(totdose[ind2], na.rm = TRUE)/100))
with(medianEffect, plot(totdose[ind2],fa[ind2], type = 'n',
xlab = 'Dose',ylab = 'Relative viability',
ylim = ylim, xlim = c(0, max(dose))))
with(medianEffect, lines(dose, with(medianEffect,
1.0/(1.0+(dm1/dose)^(summary(lm1)$coef[2,1]))),
type = "l", lty = 3))
with(medianEffect, points(dose1[dose2==0], fa[dose2==0], pch = 1))
with(medianEffect, lines(dose, 1.0/(1.0+(dm2/dose)^(summary(lm2)$coef[2,1])),
type = "l", lty = 2))
with(medianEffect, points(dose2[dose1==0], fa[dose1==0], pch = 2))
with(medianEffect, lines(dose, 1.0/(1.0+(dmcomb/dose)^(summary(lmcomb)$coef[2,1])),
type = "l", lty = 1))
with(medianEffect, points(totdose[ind.ratio],fa[ind.ratio], pch = 16))
title( "Dose-Response Curves" )
with(medianEffect, legend(legend.position, c(name1, name2, 'Mixture'),
pch = c(1, 2, 16), lty = c(3, 2, 1), bty = 'n', cex = cex.legend))
}
} else if(type == 'contour'){
## The contour plot of raw data.
temp1 <- sort(unique(medianEffect$dose1))
temp2 <- sort(unique(medianEffect$dose2))
isAllPossible <- length(temp1)*length(temp2)==length(medianEffect$fa) # all possible combination design
if(isAllPossible){
obsfa <- matrix(medianEffect$fa, length(temp1), length(temp2), byrow = T)
contour(temp1, temp2, obsfa, levels = contour.level,
xlab = sprintf('%s dose', medianEffect$name1),
ylab = sprintf('%s dose', medianEffect$name2))
title("Contour Plot")
} else {
ww <- sprintf('%s has %d doses; %s has %d doses; only %d total responses; This is not a all-possible design and thus no contour plot can be constructed!',
medianEffect$name1, length(temp1),
medianEffect$name2, length(temp2), length(medianEffect$fa))
warning(ww)
}
}
}
#' Estimate the interaction index as well as its confidence interval using delta method for
#' fixed ratio design
#'
#' This code is extracted from the source code distributed at
#' https://biostatistics.mdanderson.org/SoftwareDownload/
#'
#' Two versions are included in the source code: CI_IIV2 2008.SSC and CI_IIV2.SSC. At first
#' we implement CI_IIV2.SSC; this gives wider CI band; so we decide to use CI_IIV2 2008.SSC.
#' Further, CI_IIV2 2008.SSC has a more recent date and thus most updated.
#'
#' @param d1 dose for drug 1 where drug 2 has dose equal to 0.
#' @param e1 corresponding scaled response for drug 1, must between 0 and 1.
#' @param d2 dose for drug 2 where drug 1 has dose equal to 0.
#' @param e2 corresponding scaled response for drug 2, must between 0 and 1.
#' @param d12 the sum of doses from drug 1 and drug 2 when both drugs are administered.
#' @param e12 corresponding scaled response when both drugs are administered.
#' @param d2.d1 the ratio of the ray design. This is the fixed ratio of dose 2 divided by dose 1.
#' @param E a vector of responses (between 0 and 1) where IAI and confidence interval are to be
#' computed from.
#' @param alpha significance level of confidence interval.
#' @param min response values close to 0 and 1 will lead to extreme values when logit
#' transformation is applied. min is used to truncate the response value so that values less than
#' min are truncated at min.
#' @param max response values close to 0 and 1 will lead to extreme values when logit
#' transformation is applied. max is used to truncate the response value so that values larger than
#' max are truncated at max.
#'
#' @return a data frame with columns IAI, IAI.low, IAI.up, E, dx1 (corresponding dose of drug 1),
#' dx2 (corresponding dose of drug 1), dx12 (corresponding dose of combined drug, same as
#' definition of d12).
#' @references Lee, J. J., & Kong, M. (2009). Confidence intervals of interaction index for
#' assessing multiple drug interaction. Statistics in biopharmaceutical research, 1(1), 4-17.
#' @importFrom stats qt
CI.delta <- function(d1, e1, d2, e2, d12, e12, d2.d1, E,
alpha = 0.05, min = 0.02, max = 0.98){
# min and max are used to truncate the effect so that logit transform would not be obsurd; the original code is not very robust!
e1 <- truncate_effect(e1, min = min, max = max)
e2 <- truncate_effect(e2, min = min, max = max)
e12 <- truncate_effect(e12, min = min, max = max)
# if e1, e2 or e12 has value of 0 or 1, the logit would be Inf and leads to error for lm.
# thus use safe log
lm1 <- lm(log(e1/(1-e1))~log(d1))
dm1 <- exp(-summary(lm1)$coef[1,1]/summary(lm1)$coef[2,1])
lm2 <- lm(log(e2/(1-e2))~log(d2))
dm2 <- exp(-summary(lm2)$coef[1,1]/summary(lm2)$coef[2,1])
lmcomb <- lm(log(e12/(1-e12))~log (d12))
dm12 <- exp(-summary(lmcomb)$coef[1,1]/summary(lmcomb)$coef[2,1])
Dx1 <- dm1*(E/(1-E))^(1/summary(lm1)$coef[2,1])
Dx2 <- dm2*(E/(1-E))^(1/summary(lm2)$coef[2,1])
dx12 <- dm12*(E/(1-E))^(1/summary(lmcomb)$coef[2,1])
iix <- (dx12/(1+d2.d1))/Dx1+(dx12*d2.d1/(1+d2.d1))/Dx2
lm1.s <-summary(lm1)
lm2.s <-summary(lm2)
lm12.s <-summary(lmcomb)
c1 <- 1.0/lm1.s$coef[2,1]^2*lm1.s$coef[1,2]^2
temp <- - mean(log(d1))*lm1.s$coef[2,2]^2
c1 <- c1+2.0*(log(E/(1-E))-lm1.s$coef[1,1])/lm1.s$coef[2,1]^3*temp
c1 <- c1+(log(E/(1-E))-lm1.s$coef[1,1])^2/lm1.s$coef[2,1]^4*lm1.s$coef[2,2]^2
c2 <- 1.0/lm2.s$coef[2,1]^2*lm2.s$coef[1,2]^2
temp <- - mean(log(d2))*lm2.s$coef[2,2]^2
c2 <- c2+2.0*(log(E/(1-E))-lm2.s$coef[1,1])/lm2.s$coef[2,1]^3*temp
c2 <- c2+(log(E/(1-E))-lm2.s$coef[1,1])^2/lm2.s$coef[2,1]^4*lm2.s$coef[2,2]^2
c12 <- 1.0/lm12.s$coef[2,1]^2*lm12.s$coef[1,2]^2
temp <- - mean(log(d12))*lm12.s$coef[2,2]^2
c12 <- c12+2.0*(log(E/(1-E))-lm12.s$coef[1,1])/lm12.s$coef[2,1]^3*temp
c12 <- c12+(log(E/(1-E))-lm12.s$coef[1,1])^2/lm12.s$coef[2,1]^4*lm12.s$coef[2,2]^2
var.ii <-((dx12/Dx1)^2*c1+(dx12*d2.d1/Dx2)^2*c2+(1.0/Dx1+d2.d1/Dx2)^2*dx12^2*c12)/(1+d2.d1)^2
t975 <- qt(1-alpha/2,length(d1)+length(d2)+length(d12)-6)
iix.low1 <- iix*exp(-t975*var.ii^0.5/iix)
iix.up1 <- iix*exp(t975*var.ii^0.5/iix)
# add the corresponding dose for visualization
IAI <- data.frame(E = E, IAI = iix, IAI.low = iix.low1, IAI.up = iix.up1,
dx1 = Dx1, dx2 = Dx2, dx12 = dx12)
return(IAI)
}
#' Estimate the interaction index as well as its confidence interval using delta method
#' for fixed ratio design
#'
#' This code is extracted from the source code distributed at
#' https://biostatistics.mdanderson.org/SoftwareDownload/
#'
#' Two papers have been published by Lee et al, one in 2007 (Lee2007) and on in 2009 (Lee2009).
#' The Lee2007 paper described five methods to assess interaction: (1) Lowewe additivity model
#' using interaction index (IAI) (2) Model of Greco et al 1990. This approach uses
#' \deqn{\alpha}{alpha} as the metric and it can be related to IAI (3) Model of Machado and
#' Robinson which uses a metric denoted as \deqn{\eta}{eta} (4) Model of Plummer and Short which
#' can also be linked to IAI through the parameter \deqn{\beta_4}{beta_4} (5) Model of Carter
#' et al that can be linked to IAI through the parameter \deqn{\beta_{12}}{beta_12}. For more
#' details of these models, please refer to Lee2007.
#'
#' The Lee2009 paper provided generalization of IAI to multiple drugs using Lowewe additivity
#' model and assumption of Chou and Talalay's median effect equation. The Chou and Talalay's
#' median effect equation can be expressed as:
#' \deqn{log(E/(1-E))=m(log d - log Dm)}{log(E/(1-E))=m(log d - log Dm)}
#' where E is the effect at dose d for a compound whose median effective dose.
#'
#' Some notes about experiment design. Usually the data is either fixed ratio design (ray design)
#' or grid design which means all-possible combination of drug concentrations between two drugs
#' are available. The Lee2007 paper provided an example of grid design. However, specific fixed
#' ratio is used to fit the median effect model which is the basis to estimate IAI. The Lee2009
#' paper considered with fixed ratio design.
#'
#' @param d1 dose for drug 1.
#' @param d2 dose for drug 2.
#' @param e corresponding response in the range [0, 1].
#' @param E a vector of responses (between 0 and 1) where IAI and confidence interval are to be
#' computed from.
#' @param name1 name of drug 1.
#' @param name2 name of drug 2.
#' @param alpha significance level of confidence interval.
#' @param d2.d1.force a ratio passed to detect_ray_design() function so as to specify a fixed
#' ratio for grid design.
#' @param tol tolerance in declaring fixed ratio.
#'
#' @return a data frame with columns IAI, IAI.low, IAI.up, E, dx1 (corresponding dose of drug 1),
#' dx2 (corresponding dose of drug 1).
#'
#' @references Lee, J. J., & Kong, M. (2009). Confidence intervals of interaction index for
#' assessing multiple drug interaction. Statistics in biopharmaceutical research, 1(1), 4-17.
#' @references Lee, J. Jack, et al (2007). Interaction index and different methods for determining
#' drug interaction in combination therapy. Journal of biopharmaceutical statistics 17.3 461-480.
#'
#' @export
#' @examples
#' fit_allPoss_1 <- fitIAI(d1 = nl22B2$schd, d2 = nl22B2$hpr, e = nl22B2$y1,
#' name1 = 'SCH66336', name2 = '4HPR', d2.d1.force = 1)
fitIAI <- function(d1, d2, e, E = seq(0.05, 0.95, 0.005),
name1 = 'Drug A', name2 = 'Drug B', alpha = 0.05,
d2.d1.force = NA, tol = 0.01){
res_design <- detect_ray_design(d1 = d1, d2 = d2, e = e, d2.d1.force = d2.d1.force, tol = tol)
# fit median-effect model
medianEffect <- fit_median_efect(d1 = d1, d2 = d2, e = e, d2.d1 = res_design$d2.d1,
name1 = name1, name2 = name2)
if(res_design$isRayDesign != TRUE) { # not fixed ratio design, fitCI is NULL
warning("This is not a ray design; computation is done with method derived from fixed ratio!")
} else {
warning(sprintf("Fixed ratio design detected. Fixed ratio is %.3f\n",
res_design$d2.d1))
}
## fit CI
fitCI <- with(res_design,
CI.delta(d1 = d1, e1 = e1, d2 = d2, e2 = e2, d12 = d12,
e12 = e12, d2.d1, E = E, alpha = alpha))
meta <- list(d1 = d1, d2 = d2, e = e, E = E, alpha = alpha,
d2.d1 = res_design$d2.d1, isRayDesign = res_design$isRayDesign)
return(list(CI = fitCI, meta = meta, medianEffect = medianEffect))
}
#' Visualize interaction index
#'
#' @param fit the fitted result from fitIAI()
#' @param type type of plot from c('IAI', 'medianEffect', 'doseResponseCurve', 'contour')
#' @param contour.level contour level. only effective if type='contour'
#' @param ylim y axis limit
#' @param mode specify if to plot against response, dose or both. only effective if type=='IAI'. can be either 'response', 'dose', or 'both'
#' @param logd logd passed to plot_median_effect
#' @param align align passed to plot_median_effect
#' @param cex.legend cex.legend passed to plot_median_effect
#' @param legend.position legend.position passed to plot_median_effect
#'
#' @importFrom graphics mtext
#'
#' @export
plotIAI <- function(fit, type = c('IAI', 'medianEffect', 'doseResponseCurve', 'contour'),
contour.level = (1:9)/10, ylim = NULL, mode = 'both',
logd = FALSE, align = TRUE, cex.legend = 1, legend.position = 'topright'){
if(type == 'IAI' & !is.null(fit$CI)){
resCI <- fit$CI
if(is.null(ylim)) ylim <- c(0.01, 10)
op <- par(mar=c(8, 5, 4, 4) + 0.1)
lty_E <- 4
# CI
ind.CId.l <- resCI$IAI.low >= min(resCI$IAI)
ind.CId.u <- resCI$IAI.up <= max(resCI$IAI)+1
xlim_log10d <- range(log10(resCI$dx12))
col_d <- 'blue2'
lty_d <- 5
if(mode == 'both'){
# IAI ~ E
fa <- fit$meta$e
E <- resCI$E
with(resCI, plot(E, IAI, log = "y", type = "l", xlab = "", xaxt = 'n',
ylab = "Interaction Index",
xlim = c(min(E, fa),max(E, fa)), ylim = ylim,
cex = 0.6, main = 'Interaction plot'))
axis(1, at = pretty(range(E)), col = "black", lwd = 1)
mtext("Relative viability", side = 1, col = "black", line = 2)
abline(h = 1, lty = 4)
lines(E[ind.CId.l], resCI$IAI.low[ind.CId.l], lty = lty_E)
lines(E[ind.CId.u], resCI$IAI.up[ind.CId.u], lty = lty_E)
### add dose info
par(new = TRUE)
## IAI ~ dose
plot(resCI$dx12, resCI$IAI, log = "xy", type = "l", ylim = ylim,
axes = F, xlab = '', ylab = '', lty = lty_d, col = col_d)
axis(1, xlim = xlim_log10d, lwd = 1, line = 3.2,
col = col_d, col.axis = col_d, lty = lty_d)
mtext("Dose", side = 1, col = col_d, line = 5.2)
lines(resCI$dx12[ind.CId.l], resCI$IAI.low[ind.CId.l], lty = lty_d, col = col_d)
lines(resCI$dx12[ind.CId.u], resCI$IAI.up[ind.CId.u], lty = lty_d, col = col_d)
} else if(mode == 'response'){
# IAI ~ E
fa <- fit$meta$e
E <- resCI$E
plot(resCI$E, resCI$IAI, log = "y", type = "l", xlab = "", xaxt = 'n',
ylab = "Interaction Index", xlim = c(min(E, fa),max(E, fa)),
ylim = ylim,cex = 0.6, main = 'Interaction plot')
abline(h = 1, lty = 4)
axis(1, at = pretty(range(resCI$E)), col = "black", lwd = 1)
mtext("Relative viability", side = 1, col = "black", line = 2)
abline(h = 1, lty = 4)
lines(E[ind.CId.l], resCI$IAI.low[ind.CId.l], lty = lty_E)
lines(E[ind.CId.u], resCI$IAI.up[ind.CId.u], lty = lty_E)
} else if(mode == 'dose'){
## IAI ~ dose
plot(resCI$dx12, resCI$IAI, log = "xy", type = "l", xlab = "Dose",
ylab = "Interaction Index", ylim = ylim, cex = 0.6, main = 'Interaction plot')
abline(h = 1, lty = 4)
lines(resCI$dx12[ind.CId.l], resCI$IAI.low[ind.CId.l], lty = lty_d)
lines(resCI$dx12[ind.CId.u], resCI$IAI.up[ind.CId.u], lty = lty_d)
}
par(op)
} else {
medianEffect <- fit$medianEffect
if(type == 'contour'){
plot_median_effect(medianEffect, type = type, contour.level = contour.level)
} else {
plot_median_effect(medianEffect, type = type, logd = logd, align = align,
ylim = ylim, legend.position = legend.position, cex.legend = cex.legend)
}
}
}
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.