Solving intercept / equivalence point / crossing point of where two Weibull regressions meet

Question

I'm looking to solve the point where 2 opposite Weibull functions meet. I'm using the drc package, with a type 2 Weibull having 2 parameters in R. I've fit both lines to my data, and now I have the e and b intercepts for the functions defining each line. I would like to now solve for where the x values are equivalent.

The Weibull (type2) 2 parameter function is: $$ f(x)= 1-\exp\left[-\exp\left(b\times (\log(x) - \log(e)\right))\right] $$ exponential decay line:

           Estimate Std. Error   t-value p-value
 b:(Intercept) -2.226194   0.225339 -9.879323       0
 e:(Intercept)  1.209326   0.072042 16.786444       0

exponential increase line:

               Estimate Std. Error   t-value p-value
b:(Intercept)  1.616248   0.145047 11.142956       0
e:(Intercept)  1.837511   0.072107 25.482970       0

Ideally I would also like to get a confidence interval for this point. I have this run in R so if I could input a suggested code that would be great.

This is for a paper I'm writing, and I would add the solver of this question to this publication. In the paper, I'm using the Weibull function to model a PCR reaction.

Here is the code for the exponential increasing instance:

  # first import the raw data:
  # I pull it from Excel

> CCrelative <- read.xlsx('CC relative fold increase.xlsx', 1)
> colnames(CCrelative) <- c("MolOffTarget", "x1", "x2")
> print(CCrelative)
   #here it is after the import

      MolOffTarget    x1    x2
    1      50001.0 8.474 8.372
    2       5001.0 7.795 7.617
    3        501.0 7.090 7.291
    4         56.0 6.258 4.803
    5          6.0 2.093 1.890
    6          1.5 0.911 0.679
    7          1.1 0.480 0.508


  #now set the minimum of each column to 0 by subtracting the lowest value from each
> CClog0 <- transform(CCrelative, x1 = (x1-(CCrelative[7,2])), x2=(x2-(CCrelative[7,3])))
      #now normalize to the maximum value 
> CClogT <- transform(CClog0, x1 = (x1/(CClog0[1,2])), x2=(x2/(CClog0[1,3])))
      # now make the first column log10 scale
> CClogT <- transform(CClogT, MolOffTarget = log10(MolOffTarget))
      #now merge columns x1 and x2 into a new data frame
> CClogT2 <- data.frame(rep(CClogT$MolOffTarget, 2), c(CClogT$x1,CClogT$x2))
> colnames(CClogT2) <- c("MolOffTarget", "FAM")
> CClogT2.W2.2 <- drm(FAM ~ MolOffTarget, data = CClogT2, fct = W2.2())

Here is the code for the exponential decay:

> CCvicRelative <- read.xlsx('CCvicRelative fold increase.xlsx', 1)
> colnames(CCvicRelative) <- c("MolOffTarget", "x1", "x2")
> print(CCvicRelative)

  MolOffTarget     x1     x2
1      50001.0 -0.717 -0.706
2       5001.0 -0.565 -0.567
3        501.0 -0.360 -0.349
4         56.0 -0.001  0.584
5          6.0  1.568  1.582
6          1.5  1.767  1.822
7          1.1  1.802  1.844
      #now set the minimum of each column to 0 by subtracting the lowest value from each
> CCvic0 <- transform(CCvicRelative, x1 = (x1-(CCvicRelative[1,2])), x2=(x2-(CCvicRelative[1,3])))
      #now normalize to the maximum value 
> CCvicT <- transform(CCvic0, x1 = (x1/(CCvic0[7,2])), x2=(x2/(CCvic0[7,3])))
      # now make the first column log10 scale
> CCvicTlog <- transform(CCvicT, MolOffTarget = log10(MolOffTarget))
      #now merge columns x1 and x2 into a new data frame
> CCvicT2 <- data.frame(rep(CCvicTlog$MolOffTarget, 2), c(CCvicTlog$x1,CCvicTlog$x2))
> colnames(CCvicT2) <- c("MolOffTarget", "VIC")
> CCvicT2.W2.2 <- drm(VIC ~ MolOffTarget, data = CCvicT2, fct = W2.2())

Apologies it's coming up a bit funny on this HTML.

It took me a bit to get used to this package. I was using Bioconductor packages before.

I made a mistake with the Weibull function earlier. The correct one is now present.

Answer 1

The intersection point

The intersection point $x^{*}$ of the two functions $f_{1}(x)$ and $f_{2}(x)$ with parameters $b_{1}, e_{1}$ and $b_{2}, e_{2}$ is given by: $$ x^{*}=\exp\left[(b_{1}\log(e_{1})-b_{2}\log(e_{2}))/(b_{1}-b_{2})\right] $$

Here is a small R script depicting your example:

#-----------------------------------------------------------------------------
# The parameters
#-----------------------------------------------------------------------------

b1 <- -2.226194
e1 <- 1.209326

b2 <- 1.616248
e2 <- 1.837511

#-----------------------------------------------------------------------------
# Setting up a plot
#-----------------------------------------------------------------------------

x <- seq(0, 5, by=.01)

f <- function(x,e, b){
  1-exp(-exp(b*(log(x) - log(e))))
}

par(bg="white", cex=1.2)
plot(f(x, e1, b1)~x, type="l", ylim=c(0, 1),
las=1, lwd=2, ylab="Function value")
lines(f(x, e2, b2)~x, col="steelblue", lwd=2)

#-----------------------------------------------------------------------------
# Calculate the intersection point
#-----------------------------------------------------------------------------

x.meet <- exp((b1*log(e1)-b2*log(e2))/(b1-b2))
x.meet
[1] 1.442003

#-----------------------------------------------------------------------------
# Check if the intersection point corresponds with the graphic
#-----------------------------------------------------------------------------

abline(v=x.meet)
axis(side=1, at=x.meet, label=round(x.meet,2),
font=3, col.axis="purple")
f(x.meet, e1, b1)
[1] 0.49129
f(x.meet, e2, b2)
[1] 0.49129

---

Confidence intervals

Here is the code I've come up with to perform the non-parametric bootstrap in R(note: I have copied your data into csv-files):

#-----------------------------------------------------------------------------
# Load packages
#-----------------------------------------------------------------------------

library(drc)
library(boot)

#-----------------------------------------------------------------------------
# Load data and prepare them
#-----------------------------------------------------------------------------

path <- "F:/" # change path according to your system

CCrelative <- read.table(paste(path, "CC relative fold increase.csv", sep=""), header=T, sep=";")

CClog0 <- transform(CCrelative, x1 = (x1-(CCrelative[7,2])), x2=(x2-(CCrelative[7,3])))
CClogT <- transform(CClog0, x1 = (x1/(CClog0[1,2])), x2=(x2/(CClog0[1,3])))
CClogT <- transform(CClogT, MolOffTarget = log10(MolOffTarget))

CClogT2 <- data.frame(rep(CClogT$MolOffTarget, 2), c(CClogT$x1,CClogT$x2))

colnames(CClogT2) <- c("MolOffTarget", "FAM")
CClogT2.W2.2 <- drm(FAM ~ MolOffTarget, data = CClogT2, fct = W2.2())

CCvicRelative <-read.table(paste(path, "CCvicRelative fold increase.csv", sep=""), header=T, sep=";")

#now set the minimum of each column to 0 by subtracting the lowest value from each

CCvic0 <- transform(CCvicRelative, x1 = (x1-(CCvicRelative[1,2])), x2=(x2-(CCvicRelative[1,3])))

#now normalize to the maximum value 
CCvicT <- transform(CCvic0, x1 = (x1/(CCvic0[7,2])), x2=(x2/(CCvic0[7,3])))

# now make the first column log10 scale
CCvicTlog <- transform(CCvicT, MolOffTarget = log10(MolOffTarget))

#now merge columns x1 and x2 into a new data frame
CCvicT2 <- data.frame(rep(CCvicTlog$MolOffTarget, 2), c(CCvicTlog$x1,CCvicTlog$x2))
colnames(CCvicT2) <- c("MolOffTarget", "VIC")
CCvicT2.W2.2 <- drm(VIC ~ MolOffTarget, data = CCvicT2, fct = W2.2())

#-----------------------------------------------------------------------------
# Combine the data into one data frame
#-----------------------------------------------------------------------------

combined.data <- data.frame(
  MolOffTarget=c(CClogT2$MolOffTarget, CCvicT2$MolOffTarget),
  values=c(CClogT2$FAM, CCvicT2$VIC), 
  ind=c(rep("FAM", 14), rep("VIC", 14)))

#-----------------------------------------------------------------------------
# The bootstrap function
#-----------------------------------------------------------------------------

intersection_boot_fun <- function(data, indices) {

  data.temp <- data[indices,]

  fitFAM <- drm(data.temp$values[data.temp$ind=="FAM"] ~ data.temp$MolOffTarget[data.temp$ind=="FAM"] ,  fct = W2.2())
  fitVIC <- drm(data.temp$values[data.temp$ind=="VIC"] ~ data.temp$MolOffTarget[data.temp$ind=="VIC"] ,  fct = W2.2())

  coefFAM <- coef(fitFAM)
  coefVIC <- coef(fitVIC)

  b1 <- as.numeric(coefFAM[1])
  e1 <- as.numeric(coefFAM[2])

  b2 <- as.numeric(coefVIC[1])
  e2 <- as.numeric(coefVIC[2])

  x.intersect <- exp((b1*log(e1)-b2*log(e2))/(b1-b2))

  x.intersect

}

#-----------------------------------------------------------------------------
# Bootstrapping and inspecting the results
#-----------------------------------------------------------------------------

results <- boot(data=combined.data, statistic=intersection_boot_fun, R=50000, strata=combined.data$ind)

results

Bootstrap Statistics :
    original     bias    std. error
t1* 1.442003 0.03151577   0.1268457

#jack.after.boot(results) # influential observations

plot(results)

#-----------------------------------------------------------------------------
# BCa confidence intervals
#-----------------------------------------------------------------------------

boot.ci(results, conf=c(0.95, 0.99), type="bca")

Intervals : 
Level       BCa          
95%   ( 1.303,  1.708 )   
99%   ( 1.258,  1.844 )

The output of the bootstrap gives the original estimate, which is about 1.44 in this case. It also gives the bias, which is the difference between the mean of the bootstrap samples (which is about 1.474 in this case) and the original estimate. Technical note: The bootstrapped samples are stored in results$t and the mean could be calculated by mean(results$t). The standard deviation of the bootstrap samples is also provided, it is about 0.13. The bootstrap bias-corrected and accelerated (BCa) gives a 95%-CI from 1.30 to 1.71. We can't use a normal approximation to calculate the confidence interval, as the bootstrap samples are clearly not normal: