3

I'm trying to figure out how to calculate absolute EC50 values using the LL.3 and LL.4 (3 and 4 parameter) dose response models in the package drc, but I keep getting these errors of "Warning message:In log(exp(-tempVal/parmVec[5]) - 1) : NaNs produced" and the EC50 value is "NA".

Here is an example of the code I'm trying to run

###use rygrass dataset in drc

gr.LL.3 <- drm(ryegrass$rootl ~ ryegrass$conc, fct = LL.3()) # 3 parameter log-logistic model
gr.LL.4 <- drm(ryegrass$rootl ~ ryegrass$conc, fct = LL.4()) # 4 parameter log-logistic model

plot(gr.LL.3) #graph looks fine
plot(gr.LL.4) #graph looks fine

ED (gr.LL.3, respLev = c(50), type = "relative") # this works fine
ED (gr.LL.4, respLev = c(50), type = "relative") # this works fine

ED (gr.LL.3, respLev = c(50), type = "absolute") # this gives me "NA" for EC50 along with warning message  
ED (gr.LL.4, respLev = c(50), type = "absolute") # this gives me "NA" for EC50 along with warning message

It's not due to 0 values for concentrations

### It's not due to 0 values for concentrations
# ryegrass dataset with 0 value concentrations and corresponding rootl removed

rootlength <- c(8.3555556, 6.9142857, 7.75, 6.8714286, 6.45, 5.9222222, 1.925, 2.8857143, 4.2333333, 1.1875, 0.8571429, 1.0571429, 0.6875, 0.525, 0.825, 0.25, 0.22, 0.44)
conc.wo.0 <- c(0.94, 0.94, 0.94, 1.88, 1.88, 1.88, 3.75, 3.75, 3.75, 7.5, 7.5, 7.5, 15, 15, 15, 30, 30, 30)

gro.LL.3 <- drm(rootlength ~ conc.wo.0, fct = LL.3())

plot(gro.LL.3) #graph looks fine

ED (gro.LL.3, respLev = c(50), type = "relative") # this works fine
ED (gro.LL.3, respLev = c(50), type = "absolute") # once again, this gives me "NA" for EC50 along with warning message

It's also not due to the response being in absolute vs relative terms

### It's also not due to the response being in absolute vs relative terms
# ryegrass dataset with response relative to average response with 0 concentration (sorry, I did the absolute to relative conversion in excel, I'm still learning r)

rel.rootl <- c(0.98, 1.03, 1.07, 0.94, 0.95, 1.03, 1.08, 0.89, 1.00, 0.89, 0.83, 0.76, 0.25, 0.37, 0.55, 0.15, 0.11, 0.14, 0.09, 0.07, 0.11, 0.03, 0.03, 0.06)
concentration <- c(0, 0, 0, 0, 0, 0, 0.94, 0.94, 0.94, 1.88, 1.88, 1.88, 3.75, 3.75, 3.75, 7.5, 7.5, 7.5, 15, 15, 15, 30, 30, 30)

rel.gro.LL.3 <- drm(rel.rootl ~ concentration, fct = LL.3())

plot(rel.gro.LL.3) #graph looks fine

ED (rel.gro.LL.3, respLev = c(50), type = "relative") # this works fine
ED (rel.gro.LL.3, respLev = c(50), type = "absolute") # once again, this gives me "NA" for EC50 along with warning message

I'm new to this, so any help is appreciated.

Phillip Martin
  • 33
  • 6

1 Answers1

3
rel.rootl <- c(0.98, 1.03, 1.07, 0.94, 0.95, 1.03, 1.08, 0.89, 1.00, 0.89, 0.83, 0.76, 0.25, 0.37, 0.55, 0.15, 0.11, 0.14, 0.09, 0.07, 0.11, 0.03, 0.03, 0.06)
concentration <- c(0, 0, 0, 0, 0, 0, 0.94, 0.94, 0.94, 1.88, 1.88, 1.88, 3.75, 3.75, 3.75, 7.5, 7.5, 7.5, 15, 15, 15, 30, 30, 30)

rel.gro.LL.3 <- drm(rel.rootl ~ concentration, fct = LL.3())

plot(rel.gro.LL.3) #graph looks fine

ED (rel.gro.LL.3, respLev = c(50), type = "relative") # this works fine
ED (rel.gro.LL.3, respLev = c(50), type = "absolute") # once again, this gives me "NA" for EC50 along with warning message

The problem is because when you are trying to estimate the absolute EC50 the ED function solves for the point on the curve where you want (i.e. the respLev argument) so if your relative response level does not have 50% on the y-axis it will run into an error because your y-axis is proportions.

To fix this issue either multiply your normalized response by 100 to turn it into a percent relative response

rel.gro.LL.3.percent <- drm(rel.rootl*100 ~ concentration, fct = LL.3())

ED (rel.gro.LL.3.percent, respLev = c(50), type = "relative") # same result as above

Estimated effective doses

       Estimate Std. Error
e:1:50  3.26520    0.19915

ED (rel.gro.LL.3.percent, respLev = c(50), type = "absolute") # very similar to relative EC50

Estimated effective doses

   Estimate Std. Error
e:1:50  3.30154    0.20104

Alternatively, you could change the respLev to 0.5 in your original model. 

ED (rel.gro.LL.3, respLev = c(50), type = "relative") # this still works fine

Estimated effective doses

   Estimate Std. Error
e:1:50  3.26520    0.19915

ED (rel.gro.LL.3, respLev = c(0.5), type = "absolute") # Now this works and is the same as we got before with response multiplied by 100 

Estimated effective doses

    Estimate Std. Error
e:1:0.5  3.30154    0.20104
Zachary Noel
  • 46
  • 1
  • 1
    Thanks so much! That solves my problem. I really do appreciate it. For the benefit of other readers, check out Zachary's paper at "Noel, Z. A., Wang, J., & Chilvers, M. I. (2018). Significant Influence of EC 50 Estimation by Model Choice and EC 50 Type. Plant Disease, 102, 708–714." It's a good resource for EC50 calculations. – Phillip Martin Feb 01 '20 at 23:51