How to do 4-parametric regression for ELISA data in R

Question

I am a biology student. We do many Enzyme Linked Immunosorbent Assay (ELISA) experiments and Bradford detection. A 4-parametric logistic regression (reference) is often used for regression these data following this function: $$ F(x) = \left(\frac{A-D}{1+(x/C)^B}\right) + D $$ How can I do this in R? I want to get the $A$, $B$, $C$ and $D$ values and plot the curve.

PS. If I have some data, how can I use the calculated function $F(x)$ to get the value? I mean how do I go from "data -> F(x) -> value"?

Answer 1

The 2nd answer to a Google search for 4 parameter logistic r is this promising paper in which the authors have developed and implemented methods for analysis of assays such as ELISA in the R package drc. Specifically, the authors have developed a function LL.4() which implements the 4 paramater logistic regression function, for use with the general dose response modeling function drm.

Christian Ritz, Jens Streiberg. Bioassay Analysis Using R. Journal of Statistical Software, 2005, Vol. 12, No. 5.

Ritz et al have published a new paper that covers improvements to the 'drc' R package. Dose Response Analysis using R (PLOS ONE, 2015)

Answer 2

After days, I record my found here:

http://www.bioassay.dk/index-filer/start/DraftDrcManual.pdf gives me the current manual of drc package in R. For example:

library(drc) model1 <- drm(SLOPE~DOSE, CURVE, fct=LL.4(names=c("Slope", "Lower", "Upper", "ED50")),data=spinach) summary(model1) plot(model1) If I wanna predict the dose from observation.

model2 <- drm(DOSE~SLOPE, CURVE, fct=LL.4(names=c("Slope", "Lower", "Upper", "ED50")),data=spinach) predict(model2, newdata, type="response") newdata is a dataframe

'predict' is not the best way to estimate the DOSE from SLOPE in this case, because you have to reverse them in your model2, which doesn't work in this example.

If you want to estimate the DOSE from SLOPE, or 'Concentration' from 'OD' in case of an ELISA, just use the ED function of the 'drc' package

EXAMPLE:

library(drc)
model1 <- drm(SLOPE~DOSE, CURVE, 
              fct=LL.4(names=c("Slope", "Lower", "Upper", "ED50")),data=spinach)
plot(model1)

# the ED function is used to give the EDx value. For example, the ED50 is the 
# DOSE value for the 50% response
ED(model1,50)

# check ?ED
?ED

# The result is a matrix, from which the Estimate values can be extracted using
# the index (display=F is a good option also) 
ED(model1,50,display=F)[1:5]

# type="absolute" gives you the ability to use absolute values for the response, to 
# estimate the DOSE
response<-0.5   #lets use 0.5 for the response
DOSEx<-ED(model1,response,type="absolute",display=F)[1:5] # the estimated DOSE
points(y=rep(response,5),x=DOSEx,col="blue",pch=1:5)

Answer 3

You can find the least-square estimate of the parameters using nonlinear regression. Example:

f=function(B,x)
  (B[1]-B[4])/(1+(x/B[3])^B[2])+B[4]

LS=function(B,y,x)
  sum((y-f(B,x))^2)

x=runif(100,0,5)
B=c(1,5,2.5,5)

y=f(B,x)
plot(x,y)

### Estimate should be very close to B
nlm(LS,c(1,1,1,1),x=x,y=y)

Answer 4

There is an excellent R tutorial on fitting the 4 parameter logistic model for calibration purposes (e.g. on ELISA data) here: http://weightinginbayesianmodels.github.io/poctcalibration/over_tutorials.html http://weightinginbayesianmodels.github.io/poctcalibration/calib_tut4_curve_ocon.html http://weightinginbayesianmodels.github.io/poctcalibration/calib_tut5_precision_ocon.html http://weightinginbayesianmodels.github.io/poctcalibration/AMfunctions.html#sdXhat

They just use the nlsLM function in the minpack.lm package. I.e. they fit a a model of the form

x <- c(478, 525, 580,  650,  700,  720,  780,  825,  850,  900,  930,  980, 1020, 1040, 1050, 1075, 1081, 1100, 1160, 1180, 1200)
y <- c(1.70, 1.45, 1.50, 1.42, 1.39, 1.90, 2.49, 2.21, 2.57, 2.90, 3.55, 3.80, 4.27, 4.10, 4.60, 4.42, 4.30, 4.52, 4.40, 4.50, 4.15)

M.4pl <- function(x, lower.asymp, upper.asymp, inflec, hill){
      f <- lower.asymp + ((upper.asymp - lower.asymp)/
                            (1 + (x / inflec)^-hill))
      return(f)
    }

require(minpack.lm)
nlslmfit = nlsLM(y ~ M.4pl(x, lower.asymp, upper.asymp, inflec, hill),
                     data = data.frame(x=x, y=y),
                     start = c(lower.asymp=min(y)+1E-10, upper.asymp=max(y)-1E-10, inflec=mean(x), hill=1),
                     control = nls.control(maxiter=1000, warnOnly=TRUE) )
    summary(nlslmfit)
    # Parameters:
    #   Estimate Std. Error t value Pr(>|t|)    
    #   lower.asymp   1.5371     0.1080   14.24 7.06e-11 ***
    #   upper.asymp   4.5508     0.1497   30.40 2.93e-16 ***
    #   inflec      889.1543    14.0924   63.09  < 2e-16 ***
    #   hill         13.1717     2.5475    5.17 7.68e-05 ***
    require(investr)
    xvals=seq(min(x),max(x),length.out=100)
    predintervals = data.frame(x=xvals,predFit(nlslmfit, newdata=data.frame(x=xvals), interval="prediction"))
    confintervals = data.frame(x=xvals,predFit(nlslmfit, newdata=data.frame(x=xvals), interval="confidence"))
    require(ggplot2)
    qplot(data=predintervals, x=x, y=fit, ymin=lwr, ymax=upr, geom="ribbon", fill=I("red"), alpha=I(0.2)) +
      geom_ribbon(data=confintervals, aes(x=x, ymin=lwr, ymax=upr), fill=I("blue"), alpha=I(0.2)) +
      geom_line(data=confintervals, aes(x=x, y=fit), colour=I("blue"), lwd=2) +
      geom_point(data=data.frame(x=x,y=y), aes(x=x, y=y, ymin=NULL, ymax=NULL), size=5, col="blue") +
      ylab("y")

They also show how one could use weights and iteratively refitted least squares to allow for non-homogeneous variance. And it also covers how to do inverse prediction and calculating derived statistics like determining the limit of detection, limit of quantification and working range.

Myself I had more luck using a constrained strictly monotone P spline fit though, fitted using the scam package, to do calibration curves, as that resulted in much narrower 95% confidence intervals and prediction intervals than using the four parameter logistic model... I.e. a model of the form

require(scam)
nknots = 20 # desired nr of knots
fit = scam(y~s(conc,k=nknots,bs="mpi",m=2), family=gaussian, data=data)

If your data are discrete counts as opposed to some continuous measure you could also use family=poisson with a log link, or work with a log(y+1) transformed dependent variable. You could also use log(conc+1) in your formula as well, as concentration can never go negative.

Answer 5

After days, I record my found here:

http://www.bioassay.dk/index-filer/start/DraftDrcManual.pdf gives me the current manual of drc package in R. For example:

library(drc)
model1 <- drm(SLOPE~DOSE, CURVE, fct=LL.4(names=c("Slope", "Lower", "Upper", "ED50")),data=spinach)
summary(model1)
plot(model1)

If I wanna predict the dose from observation.

model2 <- drm(DOSE~SLOPE, CURVE, fct=LL.4(names=c("Slope", "Lower", "Upper", "ED50")),data=spinach)
predict(model2, newdata, type="response")

newdata is a dataframe

If I wanna check the regression, R square is not good for nonlinear regression

RSD <- abs(sqrt(summary(model1)$"resVar") / mean(fitted(model1)))

Thanks for Christian Ritz and my father's help.

Answer 6

For what it's worth, below is an example comparing drc::drm() and gnlm::gnlr():

library(drc)
spinach1 <- subset(spinach, CURVE==1)
model.drm <- drm(SLOPE~DOSE, CURVE,
              fct=LL.4(names=c("B", "D", "A","C")),data=spinach1)
summary(model.drm)


library(gnlm)
attach(spinach1)
model.gnlr <- gnlr(y = SLOPE, 
                   mu =~ (A-D)/(1+(DOSE/C)^B) + D,
                   pmu = list(A=0.1,B=0.1,C=0.1,D=0.1),
                   pshape=log(0.05))
model.gnlr