Dose Response - Global curve fitting using R

Question

I have the following dose response data and wish to plot dose response model and global fit curve. [xdata = drug concentration; ydata(0-5) = response values at different concentrations of the drug]. I plotted the Std Curve without an issue.

Std Curve Data fit:

df <- data.frame(xdata = c(1000.00,300.00,100.00,30.00,10.00,3.00,1.00,0.30,
                           0.10,0.03,0.01,0.00),
                 ydata = c(91.8,95.3,100,123,203,620,1210,1520,1510,1520,1590,
                           1620))

nls.fit <- nls(ydata ~ (ymax*xdata / (ec50 + xdata)) + Ns*xdata + ymin, data=df,
               start=list(ymax=1624.75, ymin = 91.85, ec50 = 3, Ns = 0.2045514))

Dose Response Curve Data fit:

df <- data.frame(
        xdata = c(10000,5000,2500,1250,625,312.5,156.25,78.125,39.063,19.531,9.766,4.883,
                 2.441,1.221,0.610,0.305,0.153,0.076,0.038,0.019,0.010,0.005),
        ydata1 = c(97.147, 98.438, 96.471, 73.669, 60.942, 45.106, 1.260, 18.336, 9.951, 2.060, 
                   0.192, 0.492, -0.310, 0.591, 0.789, 0.075, 0.474, 0.278, 0.399, 0.217, 1.021, -1.263),
        ydata2 = c(116.127, 124.104, 110.091, 111.819, 118.274, 78.069, 52.807, 40.182, 26.862, 
                   15.464, 6.865, 3.385, 10.621, 0.299, 0.883, 0.717, 1.283, 0.555, 0.454, 1.192, 0.155, 1.245),
        ydata3 = c(108.410, 127.637, 96.471, 124.903, 136.536, 104.696, 74.890, 50.699, 47.494, 23.866, 
                   20.057, 10.434, 2.831, 2.261, 1.085, 0.399, 1.284, 0.045, 0.376, -0.157, 1.158, 0.281),
        ydata4 = c(107.281, 118.274, 99.051, 99.493, 104.019, 99.582, 87.462, 75.322, 47.393, 42.459, 
                   8.311, 23.155, 3.268, 5.494, 2.097, 2.757, 1.438, 0.655, 0.782, 1.128, 1.323, 0.645),
        ydata0 = c(109.455, 104.989, 101.665, 101.205, 108.410, 101.573, 119.375, 101.757, 65.660, 35.672, 
                   31.613, 12.323, 25.515, 17.283, 7.170, 2.771, 2.655, 0.491, 0.290, 0.535, 0.298, 0.106))

When I tried to get the fit parameters using the R script provided below, I get the following error:

Error in nls(ydata1 ~ BOTTOM + (TOP - BOTTOM)/(1 + 10^((logEC50 - xdata) * :
singular gradient

nls.fit1 <- nls(ydata1 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df,
                start=list(TOP = max(df$ydata1), BOTTOM = min(df$ydata1),hillSlope = 1.0, logEC50 = 4.310345e-08))

nls.fit2 <- nls(ydata2 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df,
                start=list(TOP = max(df$ydata2), BOTTOM = min(df$ydata2),hillSlope = 1.0, logEC50 = 4.310345e-08))

nls.fit3 <- nls(ydata3 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df,
                start=list(TOP = max(df$ydata3), BOTTOM = min(df$ydata3),hillSlope = 1.0, logEC50 = 4.310345e-08))

nls.fit4 <- nls(ydata4 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df,
               start=list(TOP = max(df$ydata4), BOTTOM = min(df$ydata4),hillSlope = 1.0, logEC50 = 4.310345e-08))

nls.fit5 <- nls(ydata0 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df,
                start=list(TOP = max(df$ydata0), BOTTOM = min(df$ydata0),hillSlope = 1.0, logEC50 = 4.310345e-08))

Please advice me on how to fix this issue

G. Grothendieck G. Grothendieck · Accepted Answer · 2017-02-26T17:56:38

First note that the ratio of the largest value of xdata to the smallest is 2 million so we likely want to use log(xdata) in place of xdata.

Now, making this change we get the 4 parameter log-logistic LL2.4 model of the drc package but with a slightly different parameterization than in the question. Assuming that you are ok with these changes we can fit the first model as follows. See ?LL2.4 for the details of the parameterization and see the relevant examples at the bottom of ?ryegrass . Here df is the df shown in the question -- the LL2.4 model itself makes the log(xdata) transformation.

library(drc)

fm1 <- drm(ydata1 ~ xdata, data = df, fct = LL2.4())
fm1
plot(fm1)

Here we fit all 5 models and visually we see from the plot at the end that the fits are pretty good.

library(drc)

fun <- function(yname) {
  fo <- as.formula(paste(yname, "~ xdata"))
  fit <- do.call("drm", list(fo, data = quote(df), fct = quote(LL2.4())))
  plot(fit)
  fit
}

par(mfrow = c(3, 2))
L <- Map(fun, names(df)[-1])
par(mfrow = c(1, 1))

sapply(L, coef)

giving:

                     ydata1   ydata2   ydata3   ydata4   ydata0
    b:(Intercept)  -1.37395  -1.1411  -1.1337  -1.0633  -1.6525
    c:(Intercept)   0.70388   1.9364   1.5800   1.3751   5.7010
    d:(Intercept) 101.02741 122.0825 120.8042 108.2420 107.9106
    e:(Intercept)   6.17225   5.0686   4.3215   3.7139   3.2813

and the following graphical fits (click on the image to expand it):

Dose Response - Global curve fitting using R

2 Answers