2

I am trying to estimate IC50 values with the drc package. Plotting the model works fine. But when I use the ED I get unreasonable results.

Plot of DRC Modell

In my understanding the IC50 should by around the inflection point of the model. In this case I estimate around 0.25. But the ED calculate it as -1.43. Do I use ED wrong or did I missed something?

Thanks

Minial example with data: 

# Create Data
Conc <- c(0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000)
Response <- c(167.11246201, 53.96960486, 128.42857143, 43.67173252, 4.51975684, 0.34042553, 120.10334347, 101.14589666, 155.17629179, 35.31306991, 8.56534954, 1.71124620, 146.34954407, 108.50151976, 163.60182371, 64.70212766, 2.88145897, 0.50759878, 82.92401216, 109.80547112, 116.69300912, 26.85410334, 3.01519757, 0.37386018, 87.06990881, 84.82978723, 118.36474164, 27.52279635, 2.34650456, 0.10638298, 89.47720365, 109.47112462, 85.43161094, 17.69300912, 2.31306991, 0.07294833)

df <- data.frame(Conc = Conc, Response = Response)


# Make Modell
library(drc)

drm <- drm(Response ~ Conc, data = df, fct = LL2.4())

plot(drm, main = paste("ED(drm, 50):", ED(drm, 50)[[1]]))


# Calculate ED50 (IC50)
ED(drm, 50)
WitheShadow
  • 522
  • 3
  • 14
  • 1
    The ED is reported in log-space. Normally this would be in base 10 but for you it looks to be in base 2. 2**-1.43 ≈ 0.37. So your intuition is correct. – emilliman5 Jan 23 '17 at 14:11
  • Thanks for your help. I was trying base 10, but did not thought, it could be an other base. Why are you sure it ist not base e (euler)? – WitheShadow Jan 25 '17 at 10:31
  • In all honestly I expected it to be in base 10 as well, and am not sure why it is not (I dug into the drm code too)... Why base 2, because your doses are a 2-fold serial dilution, which means log2 produces integer doses... I really do not have a better answer than that. Sorry – emilliman5 Jan 25 '17 at 19:05

1 Answers1

0

I hope it is not too late for this response.

Because you called LL2.4() without fixing c and d parameters, the algorithm comes out with estimates for them; therefore, the mid-point between the estimated c and d needs to be calculated.

I do not know why, but when you add logBase=exp(1) to the call of ed, the results align with your expectation.

Conc = c(0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000, 0.03125, 0.06250, 0.12500, 0.25000, 0.50000, 1.00000)
Response = c(167.11246201, 53.96960486, 128.42857143, 43.67173252, 4.51975684, 0.34042553, 120.10334347, 101.14589666, 155.17629179, 35.31306991, 8.56534954, 1.71124620, 146.34954407, 108.50151976, 163.60182371, 64.70212766, 2.88145897, 0.50759878, 82.92401216, 109.80547112, 116.69300912, 26.85410334, 3.01519757, 0.37386018, 87.06990881, 84.82978723, 118.36474164, 27.52279635, 2.34650456, 0.10638298, 89.47720365, 109.47112462, 85.43161094, 17.69300912, 2.31306991, 0.07294833)

df = data.frame(Conc = Conc, Response = Response)

ll24fit = drm(Response ~ Conc, data = df, fct = LL2.4())

# Calculate ED50 (IC50)
coefValues = summary(ll24fit)[["coefficients"]]
midPoint = mean(coefValues[2:3])
ed50 = ED(ll24fit, 50,logBase=exp(1), interval="delta")
plot(ll24fit, log="", main = paste("ED(drm, 50):", ED(ll24fit, 50, logBase=exp(1))[[1]]))
abline(v=ed50[1], col="#666666", lty="longdash", lwd=0.75)
abline(h=midPoint, col="#666666", lty="longdash", lwd=0.75)

enter image description here

James Silva
  • 62
  • 2