I am new to R and I want to find the ED50 with its confidence interval in the drc package in R. I found the best model is lorentz.4 according to AIC. But I get a warning message for the confidence interval, plus during the summary.
Warning message: In sqrt(diag(varMat)) : NaNs produced
Error in ED.drc(mod.drc, c(10, 50, 90)) : ED values cannot be calculated
Any help is so appreciated. Thank you. Here is the code:
library(devtools)
install_github("onofriandreapg/aomisc")
library(aomisc)
library(drc)
X <- c(25, 50, 100, 200, 300)
Y1 <- c(0.12, 0.14,0.44, 0.14, 0.07)
mod.drc <- drm(Y1 ~ X, fct =DRC.lorentz.4() )
plot(mod.drc, ylim = c(0, 0.8), log = "")
AIC(mod.drc)
summary(mod.drc)
confint(mod.drc)
ED(mod.drc, c(10, 50, 90))
Seems that DRC.lorentz.4 does not work with ED. Calculate it yourself. We assume the same library statements.
We use the fact that for the DRC.lorentz.4 model coef(mod.drc)[[4]]
equals the X value at the maximum point of the fitted curve. Had we not had a parameter for that we could have used x.at.ymax <- optimize(Predict, range(X), maximum = TRUE)$maximum
for that X value.
Note that there is some question of whether we can really rely on the maximum of the fitted curve given that there are no data points near it.
mod.drc <- drm(Y1 ~ X, fct = DRC.lorentz.4() )
AIC(mod.drc)
plot(mod.drc)
Predict <- function(x) predict(mod.drc, data.frame(X = x))
# Take average of 1st fitted value and maximum fitted value
x.at.ymax <- coef(mod.drc)[[4]]
y50 <- (Predict(X[1]) + Predict(x.at.ymax)) / 2
# find the x value that corresponds to y50
opt.ed50 <- uniroot(function(x) Predict(x) - y50, c(X[1], x.at.ymax))
ed50 <- opt.ed50$root
abline(v = ed50, h = y50, lty = 2)