The ED50 value differs in drc when calculated vs. plotted

Question

I have a dose-response curve and I need to calculate and plot a 4-parameter logistic regression using the R package drc and ggplot. It needs to have 1 as the upper limit and 0 as the lower limit. It works fine exept that the calculated ED50 value is different from the ED50 value when read from the graph created by ggplot.

Here is a minimal example:

library(ggplot2)
library(drc)

df <- data.frame(dose = c("2000", "666.67", "222.22", 
                          "74.04", "24.69", "8.23", 
                          "2.74", "0.91", "0.30", 
                          "0.10", "0.03", "0.01"),
                 response = c("0.569767442", "0.709302326", "0.767441860",
                 "0.712209302", "0.747093023", "0.723837209", 
                 "0.71802326", "0.7558140", "0.7906977",
                 "0.7616279", "0.8197674", "1"))

df$dose <- as.numeric(df$dose)
df$response <- as.numeric(df$response)

# calculating the ED 50
logistic <- drm(response~dose, 
                data = df, 
                fct = LL.4(fixed = c(NA, 0, 1, NA)))      
ED50 <- ED(logistic, 50)[1]

# plotting
ggplot(data = df, aes(x = dose, y = response))+
  scale_x_continuous(trans = "log10", limits = c(1E-2, 2E6))+
  scale_y_continuous()+
  geom_smooth(method = drm, 
              method.args = list(fct = L.4(fixed = c(NA, 0, 1, NA)), 
              # "b", "c", "d", "e" 
              # L.4 and not LL.4 because the x scale is on a log10 
              se=FALSE,
              fullrange = TRUE))+ 
  geom_point(alpha = 0.5)+
  geom_point(aes(x = ED50, y = 0.5), color = "red") +
  coord_cartesian(xlim = c(1E-2,2E6), ylim = c(0,1))

It gives the following graph (the calculated ED50 value is marked in red).

Obviously, the calculated ED50 does not match the ED50 value suggested by the graph even though I calulated both using the LL.4 method. The plot uses a L.4 on a log axis, which should be identical to LL.4 on a non-log scale as suggested here.

I have no clue what the problem is or what I am missing. Thank you for your help and time in advance!

score 0 · Answer 1 · answered Oct 17 '22 at 11:54

This doesn't really give an answer to the question why the methods give you different ED50's, that might be more suited for CrossValidated. However, you can circumvent the issue of reconstituting your model on a log axis, by just using stat_function() + predict() instead of using geom_smooth(). That way, you can be sure that the plotted data comes from the model that also gave you the ED50 estimate.

library(ggplot2)
library(drc)

df <- data.frame(dose = c("2000", "666.67", "222.22", 
                          "74.04", "24.69", "8.23", 
                          "2.74", "0.91", "0.30", 
                          "0.10", "0.03", "0.01"),
                 response = c("0.569767442", "0.709302326", "0.767441860",
                              "0.712209302", "0.747093023", "0.723837209", 
                              "0.71802326", "0.7558140", "0.7906977",
                              "0.7616279", "0.8197674", "1"))

df$dose <- as.numeric(df$dose)
df$response <- as.numeric(df$response)

# calculating the ED 50
logistic <- drm(response~dose, 
                data = df, 
                fct = LL.4(fixed = c(NA, 0, 1, NA)))      
ED50 <- ED(logistic, 50)[1]


ggplot(df, aes(dose, response)) +
  stat_function(
    fun = function(x) predict(logistic, newdata = data.frame(dose = x))
  ) +
  geom_point(alpha = 0.5) +
  geom_point(
    aes(x = ED50, y = 0.5),
    colour = "red"
  ) +
  scale_x_continuous(trans = "log10")

^{Created on 2022-10-17 by the reprex package (v2.0.0)}

The ED50 value differs in drc when calculated vs. plotted

1 Answers1