ED50 in non linear logistic regression in R

162 views Asked by At

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))
1

There are 1 answers

0
G. Grothendieck On BEST ANSWER

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)

screenshot