英文:
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()
英文:
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()
答案1
得分: 1
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 (https://stackoverflow.com/questions/57568882/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()
英文:
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 (https://stackoverflow.com/questions/57568882/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()
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论