In gilles-guillot/MDR: Test Loewe additivity with Model Deviation Ratio
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
background=col2rgb("cornsilk3"),
highlight = TRUE
)
library(MDR)
library(ggplot2)
library(gridExtra)
# library(bookdown)
Overview
Loewe additivity
Let $x_e$ and $y_e$ be doses of compounds 1 and 2 producing in combination an effect $e$. We denote by $X_e$ and $Y_e$ the doses of compounds 1 and 2 required to produce effect $e$ alone (assuming this conditions uniquely define them, i.e. that the individual dose-response functions are bijective).
$X_e/Y_e$ quantifies the potency of compound 1 relatively to that of compound 2.
The quantity $y_e X_e/Y_e$ can be interpreted as the dose of compound 2 converted into the corresponding dose of compound 1 after accounting for difference in potencies.
Dose additivity at $ (x_e,y_e)$ is defined as the situation where
\begin{equation}
x_e + y_e X_e/Y_e = X_e
\end{equation}
more coveniently written as
\begin{equation}
x_e / X_e + y_e/Y_e = 1
\end{equation}
If this conditions holds for any $(x_e,y_e)$ eliciting effect $e$, then the isobole at level $e$ is a segment joining the points
$(X_e,0)$ and $(0,Y_e)$ in the domain $(x,y)$.
If we denote by $f()$, $g()$ and $h(,)$ the dose-response functions of compound 1, compound 2 and of the mixture respectively, then dose (Loewe) additivity holds when
\begin{equation} \label{eq:Loewe}
\frac{x}{f^{-1}(h(x,y))} + \frac{y}{g^{-1}(h(x,y))} = 1
\end{equation}
Model deviation ratio
In the MDR method, the departure from Loewe additivity is based on comparing an estimate of the mixture dose eliciting an effect $e$ (typically 50\% max effect) to an estimate of the mixture dose eliciting the same effect under the Loewe additivity model. It involves the following steps:
-
Estimate the mixture dose eliciting effect $e$
-
Estimate the mixture dose that would elicit effect $e$ if Loewe additivity was holding
-
Compute model deviation ratio of previous quantities
Fitting single-variable dose response curves
The MDR package comes with functions for fitting dose-response curves to single-compound experiment data. The main function is Fit_Single_Hill which fits curves from the 2-, 3- or 4-parameter Hill family: $$f(x) = d + \frac{c-d}{ ( 1+ (a/x)^b)}$$
This functions allows users to perform parameter estimation, boostrap compuation to quantify uncertainty attached to parameter estimates and model selection by cross-validation.
Its use is illustrated below:
## Some simulated data are stored in R object of the MDR package named simdat
## subsetting rows for which y==0
subs = simdat$y==0
## simplest instance of call to Fit_Single_Hill
res_scdrx = Fit_Single_Hill(x=simdat$x[subs],
y=simdat$z[subs], ## z contains the response
n=simdat$n[subs],
SCDR_model='Hill2')
res_scdrx$par_hat$Hill2
## same with bootstrap for
res_scdrx_mock = Fit_Single_Hill(x=simdat$x[subs],
y=simdat$z[subs], ## z contains the response
n=simdat$n[subs],
SCDR_model='Hill2',
bootstrap = TRUE,
B=5)
## next lines using res_scdrx pre-computed with B set to 1000
res_scdrx$par_hat$Hill2
# par(mfrow=c(1,2))
hist(res_scdrx_boot$par_boot$Hill2[,1],breaks=30,
main='Boostrap distribution',xlab='a')
abline(v=res_scdrx$par_hat$Hill2[1],col='firebrick3',lwd=2)
hist(res_scdrx_boot$par_boot$Hill2[,2],breaks=30,
main='',xlab='b')
abline(v=res_scdrx$par_hat$Hill2[2],col='firebrick3',lwd=2)
## fit for 2,3 and 4 param model
res_scdrx_234 = Fit_Single_Hill(x=simdat$x[subs],
y=simdat$z[subs], ## z contains the response
n=simdat$n[subs],
SCDR_model=c('Hill2','Hill3','Hill4'))
Defining and computing the MDR for various experimental designs
Ray design and factorial design
d = c(.8,2,6,10,20)
w = 1
xy = cbind(w*d,(1-w)*d)
w = .2
xy = rbind(xy,cbind(w*d,(1-w)*d))
w = .5
xy = rbind(xy,cbind(w*d,(1-w)*d))
w = .8
xy = rbind(xy,cbind(w*d,(1-w)*d))
w = 0
xy = rbind(xy,cbind(w*d,(1-w)*d))
xy = rbind(rep(1e-5,2),xy)
colnames(xy) = c('x','y')
xy = as.data.frame(xy)
p1 <- ggplot(xy,aes(x,y)) + geom_point(col="steelblue", size=3) +
geom_abline(slope=0,intercept=0,col='firebrick2',linetype='solid') +
geom_abline(slope=(1-.2)/.2,intercept=0,col='firebrick2',linetype='solid') +
geom_abline(slope=(1-.5)/.5,intercept=0,col='firebrick2',linetype='solid') +
geom_abline(slope=(1-.8)/.8,intercept=0,col='firebrick2',linetype='solid') +
geom_vline(xintercept=0,col='firebrick2',linetype='solid') +
xlab('') + ylab('Dose compound 2') +
theme(axis.text.x = element_blank(), axis.text.y = element_blank())
d = c(.8,2,6,10,20)
w = 1
xy = cbind(w*d,(1-w)*d)
w = .8
xy = rbind(xy,cbind(w*d,(1-w)*d))
w = 0
xy = rbind(xy,cbind(w*d,(1-w)*d))
xy = rbind(rep(1e-5,2),xy)
colnames(xy) = c('x','y')
xy = as.data.frame(xy)
p3 <- ggplot(xy,aes(x,y)) + geom_point(col="steelblue", size=3) +
geom_abline(slope=0,intercept=0,col='firebrick2',linetype='solid') +
geom_abline(slope=(1-.8)/.8,intercept=0,col='firebrick2',linetype='solid') +
# geom_abline(slope=(1-.5)/.5,intercept=0,col='firebrick2',linetype='dashed') +
# geom_abline(slope=(1-.8)/.8,intercept=0,col='firebrick2',linetype='solid') +
geom_vline(xintercept=0,col='firebrick2',linetype='solid') +
xlab('Dose compound 1') + ylab('Dose compound 2') +
theme(axis.text.x = element_blank(), axis.text.y = element_blank())
d1 = c(0,.8,2,6,10,20)
d2 = c(0,.5,1.5,5,12,20)
x = rep(d1,6)
y = rep(d2,each=6)
xy = cbind(x,y)
xy = as.data.frame(xy)
colnames(xy) = c('x','y')
p2 <- ggplot(xy,aes(x,y)) + geom_point(col="steelblue", size=3) +
geom_abline(slope=0,intercept=0,col='firebrick2',linetype='solid') +
geom_vline(xintercept=d1,col='firebrick2',linetype='solid') +
geom_hline(yintercept=d2,col='firebrick2',linetype='solid') +
geom_vline(xintercept=0,col='firebrick2',linetype='solid') +
xlab('') + ylab('') +
theme(axis.text.x = element_blank(), axis.text.y = element_blank())
d1 = c(0,.8,2,6,10,20)
d2 = c(0,5)
x = rep(d1,2)
y = rep(d2,each=2)
xy = cbind(x,y)
xy = rbind(xy,cbind(rep(0,6),c(0,.5,1.5,5,12,20)))
xy = as.data.frame(xy)
colnames(xy) = c('x','y')
p4 <- ggplot(xy,aes(x,y)) + geom_point(col="steelblue", size=3) +
geom_abline(slope=0,intercept=0,col='firebrick2',linetype='solid') +
# geom_vline(xintercept=d1,col='firebrick2',linetype='solid') +
geom_hline(yintercept=d2,col='firebrick2',linetype='solid') +
geom_vline(xintercept=0,col='firebrick2',linetype='solid') +
xlab('Dose compound 1') + ylab('') + ylim(0,20) +
theme(axis.text.x = element_blank(), axis.text.y = element_blank())
grid.arrange(p1,p2,p3,p4,nrow=2)
Single parralel design
Under this setting, it is assumed that data consist of single compound experiment data for each compound, and mixture experiment data where the dose of one compound (say compound 2) is held fixed while the dose of the other compound varies.
We denote by
-
$e$ some reference effect (typically 50\% max effect),
-
$X_e$ and $Y_e$ the doses of compounds 1 and 2 required to produce effect $e$ alone,
-
$y^0$ the dose of compound 2 held fixed while the dose of compound 1 varies,
-
$x^*_e$ the dose of compound 1 eliciting effect $e$ in presence of dose $y^0$ of compound 2.
-
$x^{DA}_e$ the dose of compound 1 that would elicit effect $e$ in presence of dose $y^0$ of compound 2 if dose addtivity was holding.
Then $x^{DA}_e = X_e(1-y^0 / Y_e)$ and
$$MDR = \frac{x^{DA}_e}{x^_e} = \frac{X_e}{x^}(1-\frac{y^0}{Y_e})$$
In the above, $y^0$ is set by the experimenter while $X_e$, $Y_e$ and $x^*_e$ are estimated from the data with any technique suitable for single-compound dose-response experiment data.
Single ray design
Under this setting, it is assumed that data consist of single compound experiment data for each compound, and mixture experiment data where the dose of both compounds varies while the ratio of their doses is held constant. On such a ray, even though both doses vary, there is only one dimension or degree of freedom. To account for this fact, it is convenient to re-parameterize $(x_i,y_i)$ into
- $(x_i,y_i)= d_i(w,1-w)$ where $d_i=x_i+y_i$ and $w=x_i/(x_i+y_i)$
or
- $(x_i,y_i) = y_i(a,1)$ where $a=x_i/y_i$.
Using one of those parametrizations makes it possible to use methods designed for single compound experiment data.
Parameterization 1
Using the parameterization $(x_i,y_i)=d_i(w,1-w)$, the dose additivity equation $x_e / X_e + y_e/Y_e = 1$ can be written
$$d_e w / X_e + d_e(1-w)/Y_e = 1$$
Denoting by $d^{DA}_e$ the mixture dose that would elicit effect $e$ if DA was holding and by $d_e^$ the mixture dose estimated to elicit effect $e$, we have:
$$d^{DA}_e = \frac{1}{ w/X_e + (1-w)/Y_e }$$
and $$MDR = \frac{1}{ d_e^(w/X_e + (1-w)/Y_e) }$$
Parameterization 2
Using the parameterization $(x_i,y_i) = y_i(a,1)$, we denote by $y_e^{DA}$ the dose of compound 2 that would elicit effect $e$ in presence of a dose $ay_e^{DA}$ of compound 1 if DA was holding and by $y_e^*$ the dose of compound 2 estimated to elicit
effect $e$ in presence of a dose $ay_e^{DA}$ of compound 1.
The DA equation $x_e / X_e + y_e/Y_e = 1$ is rewritten as
$$a y_e / X_e + y_e/Y_e = 1$$
hence
$$y_e^{DA} = \frac{1}{a/X_e + 1/Y_e}$$
and $$MDR = \frac{1}{y_e^*(a/X_e + 1/Y_e)}$$
Testing Loewe additivity with the MDR statistic
Examples of estimation of the MDR
Analysis of simulated data
library(MDR)
str(simdat)
res_mdr = Compute_MDR(x=simdat$x,
y=simdat$y,
z=simdat$z,
n=simdat$n,
SCDR_model=c('Hill2'),
design='single-ray',
lower = c(0.0001,0.0001,0.0001,0.0001) ,
upper = c(10,10,1,1) ,
do.plot=FALSE )
res_mdr$X50_hat # EC50 for first compound
res_mdr$Y50_hat # EC50 for second compound
res_mdr$d50_DA # EC50 for mixture predicted by Loewe additivity equation
## (cf vignette for details)
res_mdr$d50_obs # EC50 for mixture estimated from data
res_mdr$MDR ## res_mdr$d50_DA/res_mdr$d50_obs
## mock call to Resample_MDR with n_resamp=3
## n_resamp should be set to 200 at least
res_resamp_mock = Resample_MDR(x=simdat$x,y=simdat$y,
z=simdat$z,n=simdat$n,
SCDR_model='Hill2',
par_Hill=NULL,
design='single-ray',
lower = c(0.0001,0.0001,0.0001,0.0001) ,
upper = c(10,10,1,1) ,
n_resamp=3,
distribution='binom',
verbose=FALSE)
## lines below using pre-computed output of Resample_MDR
## coming as object res_resamp in MDR package
## with n_resamp=1000
## alpha/2 and 1-alpha/2 quantiles
quantile(res_resamp$MDR_resamp,p=c(0.025,0.975))
res_mdr$MDR
## p-value for observed MDR
cdf_mdr = ecdf(res_resamp$MDR_resamp)
pval = ifelse(res_mdr$MDR > median(res_resamp$MDR_resamp),
2*(1-cdf_mdr(res_mdr$MDR)),
2*cdf_mdr(res_mdr$MDR))
pval
Re-analysis of data of @dalla2014sensitivity
Single-compound dose-response curves
A subset of the numerical data contained as table 1 and 2 in @dalla2014sensitivity is provided as R objects D.magna.sc and D.magna.mix. These data are used below to show how the MDr can be computed together with an estimate of the uncertainty.
## Attempt to reproduce approximately figure 2 in Dalla Bona et al. 2014
# par(mfrow=c(1,3))
for(i in sort((1:10)[D.magna.sc[,1] %in% c('SGD','EFX','CPX')],decreasing=TRUE))
{
r = (D.magna.sc[i,'max'] / D.magna.sc[i,'min'])^(1/4)
r
doses = D.magna.sc[1,'min']*c(1,r,r^2,r^3,r^4)
dr_func = function(doses,EC50,Hill.Slope){ 1/(1+10^((log10(EC50) - log10(doses))*Hill.Slope)) }
dd = seq(1,500,l=1000)
plot(dd,dr_func(doses=dd,EC50=D.magna.sc[i,'EC50'],Hill.Slope = D.magna.sc[i,'Hill.Slope']),
type='l',log='x',
xlab=expression(paste(g.,L ^{-1})),
ylab = expression(paste(Immobilizartion,"%",sep=" ")),
main=D.magna.sc[i,1],lwd=2)
}
Mixture ED50
## Attempt to reproduce approximately panels e-h in figure 2 of Dalla Bona et al. 2014
#par(mfrow=c(2,2),
par(mai=rep(.4,4))
for(i in 1:4)
{
names = unique(D.magna.mix[,1:2])[i,]
name.a = unlist(names[1]) ; name.b = unlist(names[2])
i.a = D.magna.sc[,1] == name.a
r = (D.magna.sc[i.a,'max'] / D.magna.sc[i.a,'min'])^(1/4)
dsc.a = D.magna.sc[i.a,'min'] * c(1,r,r^2,r^3,r^4)
i.b = D.magna.sc[,1] == name.b
r = (D.magna.sc[i.b,'max'] / D.magna.sc[i.b,'min'])^(1/4)
dsc.b = D.magna.sc[i.b,'min'] * c(1,r,r^2,r^3,r^4)
plot(dsc.a,dsc.b,type="n", xlab=name.a, ylab=name.b,cex=1.5,
xlim=c(0,D.magna.sc[i.a,'max']),
ylim=c(0,D.magna.sc[i.b,'max']))
points(dsc.a,rep(0,5),pch=16,cex=1.5,col='chartreuse4') ; abline(h=0,col='chartreuse4',lty=2)
points(rep(0,5),dsc.b,pch=16,cex=1.5,col='chartreuse4') ; abline(v=0,col='chartreuse4',lty=2)
segments(D.magna.sc[D.magna.sc$Compound==name.a,'EC50'],0,
0,D.magna.sc[D.magna.sc$Compound==name.b,'EC50'],col='red',lty=1)
## loop on various mixtures proportions for this particular combination of compounds
for(j in (1:12)[D.magna.mix[,1] == name.a & D.magna.mix[,2] == name.b] )
{
w.max = D.magna.mix[j,'max.b'] / D.magna.mix[j,'max.a']
w.min = D.magna.mix[j,'min.b'] / D.magna.mix[j,'min.a']
r.a = (D.magna.mix[j,'max.a'] / D.magna.mix[j,'min.a'])^(1/4)
dm.a = D.magna.mix[j,'min.a'] * c(1,r.a,r.a^2,r.a^3,r.a^4)
r.b = (D.magna.mix[j,'max.b'] / D.magna.mix[j,'min.b'])^(1/4)
dm.b = D.magna.mix[j,'min.b'] * c(1,r.b,r.b^2,r.b^3,r.b^4)
points(cbind(dm.a,dm.b),pch=16,cex=1.5,col='chartreuse4')
abline(a=0,b=w.max,col='chartreuse4',lty=2)
points(D.magna.mix[j,'EC50.a'], D.magna.mix[j,'EC50.b'],pch=1,cex=1.5,col='red')
segments(D.magna.mix[j,'EC50.a.low'],D.magna.mix[j,'EC50.b.low'],
D.magna.mix[j,'EC50.a.up'],D.magna.mix[j,'EC50.b.up'],col='red',lwd=3)
}
}
n_resamp = 3 ## set to 3 to save time when building package,
## n_resamp should be set to e.g. 1000
## will require around 30 min of computation on a standard laptop
MDR_sim_mock = matrix(nr=nrow(D.magna.mix),nc=n_resamp)
MDR_obs = rep(NA,nrow(D.magna.mix))
for(i in 1:4) ## loop on the four compound combinations
{
names = unique(D.magna.mix[,1:2])[i,]
name.a = unlist(names[1]) ; name.b = unlist(names[2])
i.a = D.magna.sc[,1] == name.a # where to find info in single compound data D.magna.sc
r = (D.magna.sc[i.a,'max'] / D.magna.sc[i.a,'min'])^(1/4)
dsc.a = D.magna.sc[i.a,'min'] * c(1,r,r^2,r^3,r^4)
i.b = D.magna.sc[,1] == name.b
r = (D.magna.sc[i.b,'max'] / D.magna.sc[i.b,'min'])^(1/4)
dsc.b = D.magna.sc[i.b,'min'] * c(1,r,r^2,r^3,r^4)
for(j in (1:12)[D.magna.mix[,1] == name.a & D.magna.mix[,2] == name.b] ) # loop on compound combination x mixture proportions
{
w.max = D.magna.mix[j,'max.b'] / D.magna.mix[j,'max.a']
r.a = (D.magna.mix[j,'max.a'] / D.magna.mix[j,'min.a'])^(1/4)
dm.a = D.magna.mix[j,'min.a'] * c(1,r.a,r.a^2,r.a^3,r.a^4)
r.b = (D.magna.mix[j,'max.b'] / D.magna.mix[j,'min.b'])^(1/4)
dm.b = D.magna.mix[j,'min.b'] * c(1,r.b,r.b^2,r.b^3,r.b^4)
w = D.magna.mix$min.a/(D.magna.mix$min.a + D.magna.mix$min.b)
EC50.DA = 1/(w[j]/D.magna.sc$EC50[i.a] + (1-w[j])/D.magna.sc$EC50[i.b])
EC50.obs = D.magna.mix$EC50.a[j] + D.magna.mix$EC50.b[j]
MDR_obs[j] = EC50.DA / EC50.obs
## put together info on design in form of three vectors (x,y,n)
x = c(dsc.a, rep(0,length(dsc.b)), dm.a)
y = c(rep(0,length(dsc.a)), dsc.b, dm.b)
n = rep(20,length(x)) # OECD 202 (cf section 2.3 of Dallabona)
res = Resample_MDR(x=x,y=y,n=n,
SCDR_model = "Hill2",
par_Hill=c(a1=D.magna.sc$EC50[i.a],
b1=D.magna.sc$Hill.Slope[i.a],
a2=D.magna.sc$EC50[i.b],
b2=D.magna.sc$Hill.Slope[i.b]),
lower = c(0.0001,0.0001,0.0001,0.0001) ,
upper = c(1000,100,1,1),
design='single-ray',
n_resamp=n_resamp,
distribution = "binom",
verbose=FALSE)
MDR_sim_mock[j,] = res$MDR_resamp
}
}
## Plotting some of the simulated MDR values
## MDR_sim has the same structure as MDR_sim_mock but was obtained with n_resamp=1000
cdf_mdr1 = ecdf(MDR_sim[1,])
pval1 = ifelse(MDR_obs[1] > median(MDR_sim[1,]),
2*(1-cdf_mdr1(MDR_obs[1])),
2*cdf_mdr1(MDR_obs[1]))
#
cdf_mdr2 = ecdf(MDR_sim[2,])
pval2 = ifelse(MDR_obs[2] > median(MDR_sim[2,]),
2*(1-cdf_mdr2(MDR_obs[2])),
2*cdf_mdr2(MDR_obs[2]))
# par(mfrow=c(1,2))
hist(MDR_sim[1,],breaks=30,
main='Histogram of resampled \n MDR values',
sub=rownames(MDR_sim)[1],
xlab=paste('pval=',signif(pval1,dig=2)))
abline(v=MDR_obs[1],lwd=2,col=3)
abline(v=quantile(MDR_sim[1,],p=0.025),lwd=2,lty=2,col='orange')
abline(v=quantile(MDR_sim[1,],p=0.975),lwd=2,lty=2,col='orange')
hist(MDR_sim[2,],breaks=30,main='',
sub=rownames(MDR_sim)[2],
xlab=paste('pval=',signif(pval2,dig=2)))
abline(v=MDR_obs[2],lwd=2,col=3)
abline(v=quantile(MDR_sim[2,],p=0.025),lwd=2,lty=2,col='orange')
abline(v=quantile(MDR_sim[2,],p=0.975),lwd=2,lty=2,col='orange')
References
gilles-guillot/MDR documentation built on Jan. 21, 2020, 8:09 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.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", background=col2rgb("cornsilk3"), highlight = TRUE )
library(MDR) library(ggplot2) library(gridExtra) # library(bookdown)
Overview
Loewe additivity
Let $x_e$ and $y_e$ be doses of compounds 1 and 2 producing in combination an effect $e$. We denote by $X_e$ and $Y_e$ the doses of compounds 1 and 2 required to produce effect $e$ alone (assuming this conditions uniquely define them, i.e. that the individual dose-response functions are bijective). $X_e/Y_e$ quantifies the potency of compound 1 relatively to that of compound 2.
The quantity $y_e X_e/Y_e$ can be interpreted as the dose of compound 2 converted into the corresponding dose of compound 1 after accounting for difference in potencies.
Dose additivity at $ (x_e,y_e)$ is defined as the situation where \begin{equation} x_e + y_e X_e/Y_e = X_e \end{equation} more coveniently written as \begin{equation} x_e / X_e + y_e/Y_e = 1 \end{equation}
If this conditions holds for any $(x_e,y_e)$ eliciting effect $e$, then the isobole at level $e$ is a segment joining the points $(X_e,0)$ and $(0,Y_e)$ in the domain $(x,y)$.
If we denote by $f()$, $g()$ and $h(,)$ the dose-response functions of compound 1, compound 2 and of the mixture respectively, then dose (Loewe) additivity holds when
\begin{equation} \label{eq:Loewe}
\frac{x}{f^{-1}(h(x,y))} + \frac{y}{g^{-1}(h(x,y))} = 1
\end{equation}
Model deviation ratio
In the MDR method, the departure from Loewe additivity is based on comparing an estimate of the mixture dose eliciting an effect $e$ (typically 50\% max effect) to an estimate of the mixture dose eliciting the same effect under the Loewe additivity model. It involves the following steps:
-
Estimate the mixture dose eliciting effect $e$
-
Estimate the mixture dose that would elicit effect $e$ if Loewe additivity was holding
-
Compute model deviation ratio of previous quantities
Fitting single-variable dose response curves
The MDR package comes with functions for fitting dose-response curves to single-compound experiment data. The main function is Fit_Single_Hill which fits curves from the 2-, 3- or 4-parameter Hill family: $$f(x) = d + \frac{c-d}{ ( 1+ (a/x)^b)}$$
This functions allows users to perform parameter estimation, boostrap compuation to quantify uncertainty attached to parameter estimates and model selection by cross-validation.
Its use is illustrated below:
## Some simulated data are stored in R object of the MDR package named simdat ## subsetting rows for which y==0 subs = simdat$y==0 ## simplest instance of call to Fit_Single_Hill res_scdrx = Fit_Single_Hill(x=simdat$x[subs], y=simdat$z[subs], ## z contains the response n=simdat$n[subs], SCDR_model='Hill2') res_scdrx$par_hat$Hill2 ## same with bootstrap for res_scdrx_mock = Fit_Single_Hill(x=simdat$x[subs], y=simdat$z[subs], ## z contains the response n=simdat$n[subs], SCDR_model='Hill2', bootstrap = TRUE, B=5) ## next lines using res_scdrx pre-computed with B set to 1000 res_scdrx$par_hat$Hill2 # par(mfrow=c(1,2)) hist(res_scdrx_boot$par_boot$Hill2[,1],breaks=30, main='Boostrap distribution',xlab='a') abline(v=res_scdrx$par_hat$Hill2[1],col='firebrick3',lwd=2) hist(res_scdrx_boot$par_boot$Hill2[,2],breaks=30, main='',xlab='b') abline(v=res_scdrx$par_hat$Hill2[2],col='firebrick3',lwd=2) ## fit for 2,3 and 4 param model res_scdrx_234 = Fit_Single_Hill(x=simdat$x[subs], y=simdat$z[subs], ## z contains the response n=simdat$n[subs], SCDR_model=c('Hill2','Hill3','Hill4'))
Defining and computing the MDR for various experimental designs
Ray design and factorial design
d = c(.8,2,6,10,20) w = 1 xy = cbind(w*d,(1-w)*d) w = .2 xy = rbind(xy,cbind(w*d,(1-w)*d)) w = .5 xy = rbind(xy,cbind(w*d,(1-w)*d)) w = .8 xy = rbind(xy,cbind(w*d,(1-w)*d)) w = 0 xy = rbind(xy,cbind(w*d,(1-w)*d)) xy = rbind(rep(1e-5,2),xy) colnames(xy) = c('x','y') xy = as.data.frame(xy) p1 <- ggplot(xy,aes(x,y)) + geom_point(col="steelblue", size=3) + geom_abline(slope=0,intercept=0,col='firebrick2',linetype='solid') + geom_abline(slope=(1-.2)/.2,intercept=0,col='firebrick2',linetype='solid') + geom_abline(slope=(1-.5)/.5,intercept=0,col='firebrick2',linetype='solid') + geom_abline(slope=(1-.8)/.8,intercept=0,col='firebrick2',linetype='solid') + geom_vline(xintercept=0,col='firebrick2',linetype='solid') + xlab('') + ylab('Dose compound 2') + theme(axis.text.x = element_blank(), axis.text.y = element_blank()) d = c(.8,2,6,10,20) w = 1 xy = cbind(w*d,(1-w)*d) w = .8 xy = rbind(xy,cbind(w*d,(1-w)*d)) w = 0 xy = rbind(xy,cbind(w*d,(1-w)*d)) xy = rbind(rep(1e-5,2),xy) colnames(xy) = c('x','y') xy = as.data.frame(xy) p3 <- ggplot(xy,aes(x,y)) + geom_point(col="steelblue", size=3) + geom_abline(slope=0,intercept=0,col='firebrick2',linetype='solid') + geom_abline(slope=(1-.8)/.8,intercept=0,col='firebrick2',linetype='solid') + # geom_abline(slope=(1-.5)/.5,intercept=0,col='firebrick2',linetype='dashed') + # geom_abline(slope=(1-.8)/.8,intercept=0,col='firebrick2',linetype='solid') + geom_vline(xintercept=0,col='firebrick2',linetype='solid') + xlab('Dose compound 1') + ylab('Dose compound 2') + theme(axis.text.x = element_blank(), axis.text.y = element_blank()) d1 = c(0,.8,2,6,10,20) d2 = c(0,.5,1.5,5,12,20) x = rep(d1,6) y = rep(d2,each=6) xy = cbind(x,y) xy = as.data.frame(xy) colnames(xy) = c('x','y') p2 <- ggplot(xy,aes(x,y)) + geom_point(col="steelblue", size=3) + geom_abline(slope=0,intercept=0,col='firebrick2',linetype='solid') + geom_vline(xintercept=d1,col='firebrick2',linetype='solid') + geom_hline(yintercept=d2,col='firebrick2',linetype='solid') + geom_vline(xintercept=0,col='firebrick2',linetype='solid') + xlab('') + ylab('') + theme(axis.text.x = element_blank(), axis.text.y = element_blank()) d1 = c(0,.8,2,6,10,20) d2 = c(0,5) x = rep(d1,2) y = rep(d2,each=2) xy = cbind(x,y) xy = rbind(xy,cbind(rep(0,6),c(0,.5,1.5,5,12,20))) xy = as.data.frame(xy) colnames(xy) = c('x','y') p4 <- ggplot(xy,aes(x,y)) + geom_point(col="steelblue", size=3) + geom_abline(slope=0,intercept=0,col='firebrick2',linetype='solid') + # geom_vline(xintercept=d1,col='firebrick2',linetype='solid') + geom_hline(yintercept=d2,col='firebrick2',linetype='solid') + geom_vline(xintercept=0,col='firebrick2',linetype='solid') + xlab('Dose compound 1') + ylab('') + ylim(0,20) + theme(axis.text.x = element_blank(), axis.text.y = element_blank()) grid.arrange(p1,p2,p3,p4,nrow=2)
Single parralel design
Under this setting, it is assumed that data consist of single compound experiment data for each compound, and mixture experiment data where the dose of one compound (say compound 2) is held fixed while the dose of the other compound varies. We denote by
-
$e$ some reference effect (typically 50\% max effect),
-
$X_e$ and $Y_e$ the doses of compounds 1 and 2 required to produce effect $e$ alone,
-
$y^0$ the dose of compound 2 held fixed while the dose of compound 1 varies,
-
$x^*_e$ the dose of compound 1 eliciting effect $e$ in presence of dose $y^0$ of compound 2.
-
$x^{DA}_e$ the dose of compound 1 that would elicit effect $e$ in presence of dose $y^0$ of compound 2 if dose addtivity was holding.
Then $x^{DA}_e = X_e(1-y^0 / Y_e)$ and $$MDR = \frac{x^{DA}_e}{x^_e} = \frac{X_e}{x^}(1-\frac{y^0}{Y_e})$$ In the above, $y^0$ is set by the experimenter while $X_e$, $Y_e$ and $x^*_e$ are estimated from the data with any technique suitable for single-compound dose-response experiment data.
Single ray design
Under this setting, it is assumed that data consist of single compound experiment data for each compound, and mixture experiment data where the dose of both compounds varies while the ratio of their doses is held constant. On such a ray, even though both doses vary, there is only one dimension or degree of freedom. To account for this fact, it is convenient to re-parameterize $(x_i,y_i)$ into
- $(x_i,y_i)= d_i(w,1-w)$ where $d_i=x_i+y_i$ and $w=x_i/(x_i+y_i)$
or
- $(x_i,y_i) = y_i(a,1)$ where $a=x_i/y_i$.
Using one of those parametrizations makes it possible to use methods designed for single compound experiment data.
Parameterization 1
Using the parameterization $(x_i,y_i)=d_i(w,1-w)$, the dose additivity equation $x_e / X_e + y_e/Y_e = 1$ can be written
$$d_e w / X_e + d_e(1-w)/Y_e = 1$$ Denoting by $d^{DA}_e$ the mixture dose that would elicit effect $e$ if DA was holding and by $d_e^$ the mixture dose estimated to elicit effect $e$, we have: $$d^{DA}_e = \frac{1}{ w/X_e + (1-w)/Y_e }$$ and $$MDR = \frac{1}{ d_e^(w/X_e + (1-w)/Y_e) }$$
Parameterization 2
Using the parameterization $(x_i,y_i) = y_i(a,1)$, we denote by $y_e^{DA}$ the dose of compound 2 that would elicit effect $e$ in presence of a dose $ay_e^{DA}$ of compound 1 if DA was holding and by $y_e^*$ the dose of compound 2 estimated to elicit effect $e$ in presence of a dose $ay_e^{DA}$ of compound 1.
The DA equation $x_e / X_e + y_e/Y_e = 1$ is rewritten as $$a y_e / X_e + y_e/Y_e = 1$$ hence $$y_e^{DA} = \frac{1}{a/X_e + 1/Y_e}$$ and $$MDR = \frac{1}{y_e^*(a/X_e + 1/Y_e)}$$
Testing Loewe additivity with the MDR statistic
Examples of estimation of the MDR
Analysis of simulated data
library(MDR) str(simdat) res_mdr = Compute_MDR(x=simdat$x, y=simdat$y, z=simdat$z, n=simdat$n, SCDR_model=c('Hill2'), design='single-ray', lower = c(0.0001,0.0001,0.0001,0.0001) , upper = c(10,10,1,1) , do.plot=FALSE ) res_mdr$X50_hat # EC50 for first compound res_mdr$Y50_hat # EC50 for second compound res_mdr$d50_DA # EC50 for mixture predicted by Loewe additivity equation ## (cf vignette for details) res_mdr$d50_obs # EC50 for mixture estimated from data res_mdr$MDR ## res_mdr$d50_DA/res_mdr$d50_obs ## mock call to Resample_MDR with n_resamp=3 ## n_resamp should be set to 200 at least res_resamp_mock = Resample_MDR(x=simdat$x,y=simdat$y, z=simdat$z,n=simdat$n, SCDR_model='Hill2', par_Hill=NULL, design='single-ray', lower = c(0.0001,0.0001,0.0001,0.0001) , upper = c(10,10,1,1) , n_resamp=3, distribution='binom', verbose=FALSE) ## lines below using pre-computed output of Resample_MDR ## coming as object res_resamp in MDR package ## with n_resamp=1000 ## alpha/2 and 1-alpha/2 quantiles quantile(res_resamp$MDR_resamp,p=c(0.025,0.975)) res_mdr$MDR ## p-value for observed MDR cdf_mdr = ecdf(res_resamp$MDR_resamp) pval = ifelse(res_mdr$MDR > median(res_resamp$MDR_resamp), 2*(1-cdf_mdr(res_mdr$MDR)), 2*cdf_mdr(res_mdr$MDR)) pval
Re-analysis of data of @dalla2014sensitivity
Single-compound dose-response curves
A subset of the numerical data contained as table 1 and 2 in @dalla2014sensitivity is provided as R objects D.magna.sc and D.magna.mix. These data are used below to show how the MDr can be computed together with an estimate of the uncertainty.
## Attempt to reproduce approximately figure 2 in Dalla Bona et al. 2014 # par(mfrow=c(1,3)) for(i in sort((1:10)[D.magna.sc[,1] %in% c('SGD','EFX','CPX')],decreasing=TRUE)) { r = (D.magna.sc[i,'max'] / D.magna.sc[i,'min'])^(1/4) r doses = D.magna.sc[1,'min']*c(1,r,r^2,r^3,r^4) dr_func = function(doses,EC50,Hill.Slope){ 1/(1+10^((log10(EC50) - log10(doses))*Hill.Slope)) } dd = seq(1,500,l=1000) plot(dd,dr_func(doses=dd,EC50=D.magna.sc[i,'EC50'],Hill.Slope = D.magna.sc[i,'Hill.Slope']), type='l',log='x', xlab=expression(paste(g.,L ^{-1})), ylab = expression(paste(Immobilizartion,"%",sep=" ")), main=D.magna.sc[i,1],lwd=2) }
Mixture ED50
## Attempt to reproduce approximately panels e-h in figure 2 of Dalla Bona et al. 2014 #par(mfrow=c(2,2), par(mai=rep(.4,4)) for(i in 1:4) { names = unique(D.magna.mix[,1:2])[i,] name.a = unlist(names[1]) ; name.b = unlist(names[2]) i.a = D.magna.sc[,1] == name.a r = (D.magna.sc[i.a,'max'] / D.magna.sc[i.a,'min'])^(1/4) dsc.a = D.magna.sc[i.a,'min'] * c(1,r,r^2,r^3,r^4) i.b = D.magna.sc[,1] == name.b r = (D.magna.sc[i.b,'max'] / D.magna.sc[i.b,'min'])^(1/4) dsc.b = D.magna.sc[i.b,'min'] * c(1,r,r^2,r^3,r^4) plot(dsc.a,dsc.b,type="n", xlab=name.a, ylab=name.b,cex=1.5, xlim=c(0,D.magna.sc[i.a,'max']), ylim=c(0,D.magna.sc[i.b,'max'])) points(dsc.a,rep(0,5),pch=16,cex=1.5,col='chartreuse4') ; abline(h=0,col='chartreuse4',lty=2) points(rep(0,5),dsc.b,pch=16,cex=1.5,col='chartreuse4') ; abline(v=0,col='chartreuse4',lty=2) segments(D.magna.sc[D.magna.sc$Compound==name.a,'EC50'],0, 0,D.magna.sc[D.magna.sc$Compound==name.b,'EC50'],col='red',lty=1) ## loop on various mixtures proportions for this particular combination of compounds for(j in (1:12)[D.magna.mix[,1] == name.a & D.magna.mix[,2] == name.b] ) { w.max = D.magna.mix[j,'max.b'] / D.magna.mix[j,'max.a'] w.min = D.magna.mix[j,'min.b'] / D.magna.mix[j,'min.a'] r.a = (D.magna.mix[j,'max.a'] / D.magna.mix[j,'min.a'])^(1/4) dm.a = D.magna.mix[j,'min.a'] * c(1,r.a,r.a^2,r.a^3,r.a^4) r.b = (D.magna.mix[j,'max.b'] / D.magna.mix[j,'min.b'])^(1/4) dm.b = D.magna.mix[j,'min.b'] * c(1,r.b,r.b^2,r.b^3,r.b^4) points(cbind(dm.a,dm.b),pch=16,cex=1.5,col='chartreuse4') abline(a=0,b=w.max,col='chartreuse4',lty=2) points(D.magna.mix[j,'EC50.a'], D.magna.mix[j,'EC50.b'],pch=1,cex=1.5,col='red') segments(D.magna.mix[j,'EC50.a.low'],D.magna.mix[j,'EC50.b.low'], D.magna.mix[j,'EC50.a.up'],D.magna.mix[j,'EC50.b.up'],col='red',lwd=3) } }
n_resamp = 3 ## set to 3 to save time when building package, ## n_resamp should be set to e.g. 1000 ## will require around 30 min of computation on a standard laptop MDR_sim_mock = matrix(nr=nrow(D.magna.mix),nc=n_resamp) MDR_obs = rep(NA,nrow(D.magna.mix)) for(i in 1:4) ## loop on the four compound combinations { names = unique(D.magna.mix[,1:2])[i,] name.a = unlist(names[1]) ; name.b = unlist(names[2]) i.a = D.magna.sc[,1] == name.a # where to find info in single compound data D.magna.sc r = (D.magna.sc[i.a,'max'] / D.magna.sc[i.a,'min'])^(1/4) dsc.a = D.magna.sc[i.a,'min'] * c(1,r,r^2,r^3,r^4) i.b = D.magna.sc[,1] == name.b r = (D.magna.sc[i.b,'max'] / D.magna.sc[i.b,'min'])^(1/4) dsc.b = D.magna.sc[i.b,'min'] * c(1,r,r^2,r^3,r^4) for(j in (1:12)[D.magna.mix[,1] == name.a & D.magna.mix[,2] == name.b] ) # loop on compound combination x mixture proportions { w.max = D.magna.mix[j,'max.b'] / D.magna.mix[j,'max.a'] r.a = (D.magna.mix[j,'max.a'] / D.magna.mix[j,'min.a'])^(1/4) dm.a = D.magna.mix[j,'min.a'] * c(1,r.a,r.a^2,r.a^3,r.a^4) r.b = (D.magna.mix[j,'max.b'] / D.magna.mix[j,'min.b'])^(1/4) dm.b = D.magna.mix[j,'min.b'] * c(1,r.b,r.b^2,r.b^3,r.b^4) w = D.magna.mix$min.a/(D.magna.mix$min.a + D.magna.mix$min.b) EC50.DA = 1/(w[j]/D.magna.sc$EC50[i.a] + (1-w[j])/D.magna.sc$EC50[i.b]) EC50.obs = D.magna.mix$EC50.a[j] + D.magna.mix$EC50.b[j] MDR_obs[j] = EC50.DA / EC50.obs ## put together info on design in form of three vectors (x,y,n) x = c(dsc.a, rep(0,length(dsc.b)), dm.a) y = c(rep(0,length(dsc.a)), dsc.b, dm.b) n = rep(20,length(x)) # OECD 202 (cf section 2.3 of Dallabona) res = Resample_MDR(x=x,y=y,n=n, SCDR_model = "Hill2", par_Hill=c(a1=D.magna.sc$EC50[i.a], b1=D.magna.sc$Hill.Slope[i.a], a2=D.magna.sc$EC50[i.b], b2=D.magna.sc$Hill.Slope[i.b]), lower = c(0.0001,0.0001,0.0001,0.0001) , upper = c(1000,100,1,1), design='single-ray', n_resamp=n_resamp, distribution = "binom", verbose=FALSE) MDR_sim_mock[j,] = res$MDR_resamp } }
## Plotting some of the simulated MDR values ## MDR_sim has the same structure as MDR_sim_mock but was obtained with n_resamp=1000 cdf_mdr1 = ecdf(MDR_sim[1,]) pval1 = ifelse(MDR_obs[1] > median(MDR_sim[1,]), 2*(1-cdf_mdr1(MDR_obs[1])), 2*cdf_mdr1(MDR_obs[1])) # cdf_mdr2 = ecdf(MDR_sim[2,]) pval2 = ifelse(MDR_obs[2] > median(MDR_sim[2,]), 2*(1-cdf_mdr2(MDR_obs[2])), 2*cdf_mdr2(MDR_obs[2])) # par(mfrow=c(1,2)) hist(MDR_sim[1,],breaks=30, main='Histogram of resampled \n MDR values', sub=rownames(MDR_sim)[1], xlab=paste('pval=',signif(pval1,dig=2))) abline(v=MDR_obs[1],lwd=2,col=3) abline(v=quantile(MDR_sim[1,],p=0.025),lwd=2,lty=2,col='orange') abline(v=quantile(MDR_sim[1,],p=0.975),lwd=2,lty=2,col='orange') hist(MDR_sim[2,],breaks=30,main='', sub=rownames(MDR_sim)[2], xlab=paste('pval=',signif(pval2,dig=2))) abline(v=MDR_obs[2],lwd=2,col=3) abline(v=quantile(MDR_sim[2,],p=0.025),lwd=2,lty=2,col='orange') abline(v=quantile(MDR_sim[2,],p=0.975),lwd=2,lty=2,col='orange')
References
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.