CODEWITHSUNDEEP
rggplot2drmdrc
Hayoung Kwon
Hayoung Kwon

Reputation: 1

Dose-response curve, ED50

I used the following code to plot the dose-response curve and display the ED50 value. However, the ED50 value was displayed at a value other than 50 on the y-axis. I used the ED function to find the ED50 value, but why am I getting these results? Is there a way to plot the graph correctly?

library(drc)
library(tidyverse)

response_diuron <- as.numeric(c("6.666667", "10", "20", "23.33333", "100", "100"))
diuron <- as.numeric(c("0", "100", "500", "1000", "5000", "10000"))
df <- data.frame(diuron, response_diuron)
D <- drm(response_diuron ~ diuron, data = df, fct = LL.4())
newdat <- expand.grid(diuron = exp(seq(log(0.5), log(10000), length = 500)))
pm <- predict(D, newdata = newdat, interval = "confidence")
newdat$p <- pm[, 1]
newdat$pmin <- pm[, 2]
newdat$pmax <- pm[, 3]

df$diuron0 <- df$diuron
df$diuron0[df$diuron0 == 0] <- 0.5

ed50 <- ED(D, 50, logBase = exp(1), interval = "delta")
coefs <- setNames(coef(D), c("b", "c", "d", "e"))
y50 <- predict(D, newdata = data.frame(diuron = coefs["e"]))

ggplot(df, aes(x = diuron0, y = response_diuron)) +
  geom_point(shape = 21, size = 3, stroke = 1, colour = "#3C6255") +
  geom_line(data = newdat, aes(x = diuron, y = p), lwd = 1, color = "#3C6255") +
  coord_trans(x = "log") +
  geom_segment(aes(x = coefs["e"], y = 0, xend = coefs["e"], yend = y50),
    lty = 2, colour = "gray50"
  ) +
  geom_segment(aes(x = coefs["e"], y = y50, xend = 0, yend = y50),
    lty = 2, colour = "gray50"
  ) +
  scale_x_log10(
    limits = c(100, 10000), breaks = c(100, 500, 1000, 1780, 5000, 10000),
    labels = c(100, 500, 1000, 1780, 5000, 10000)
  ) +
  ggtitle("Diuron") +
  xlab("Diuron (ug/L)") +
  ylab("Mortality (%)") +
  theme(axis.text.x = element_text(angle = 90)) +
  theme_classic()

Upvotes: 0

Views: 204

Answers (1)

TheN
TheN

Reputation: 513

I have checked your code, the parameter given in geom_segment() (both the horizontal and vertical lines) is y50, which has a value of 56.6. That's why the value is not displayed as 50 on the y-axis.

Edit 15 June 2023:

With reference to this post (Wrong ED50 with Dose Response Model using DRM), I have edited your code to get the intended results. As I am not familiar with pharmacology, I am not entirely sure of what the numbers represent though, i.e., whether they really refer to the number that you are looking for that makes scientific sense.

ed50 <- ED(D, 50, type = "absolute", interval = "delta")
coefs <- setNames(ed50[1], "e")
y50 <- predict(D, newdata = data.frame(diuron = ed50))

ggplot(df, aes(x = diuron0, y = response_diuron)) +
  geom_point(shape = 21, size = 3, stroke = 1, colour = "#3C6255") +
  geom_line(data = newdat, aes(x = diuron, y = p), lwd = 1, color = "#3C6255") +
  coord_trans(x = "log") +
  geom_segment(aes(x = coefs["e"], y = 0, xend = coefs["e"], yend = y50),
               lty = 2, colour = "gray50"
  ) +
  geom_segment(aes(x = coefs["e"], y = y50, xend = 0, yend = y50),
               lty = 2, colour = "gray50"
  ) +
  scale_x_log10(
    limits = c(100, 10000), breaks = c(100, 500, 1000, round(coefs["e"]), 5000, 10000),
    labels = c(100, 500, 1000, round(coefs["e"]), 5000, 10000)
  ) +
  ggtitle("Diuron") +
  xlab("Diuron (ug/L)") +
  ylab("Mortality (%)") +
  theme(axis.text.x = element_text(angle = 90)) +
  theme_classic()

Upvotes: 1

Related Questions