剂量反应曲线绘图疑问：ED50未对应Y轴50%值的原因及修正方法

问题描述 #

你使用以下R代码绘制Diuron的剂量-反应曲线并展示ED50值，但发现ED50对应的Y轴数值并非50，用ED()函数计算ED50后仍存在这个问题：

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

问题原因 #

• LL.4模型的参数逻辑：你用的LL.4()是四参数逻辑斯蒂模型，公式为：

  f(x) = c + (d - c)/(1 + exp(b(log(x) - log(e))))
  

  其中e参数对应的是反应达到**中间值（(c+d)/2）**的剂量，而不是绝对的50%。你的数据中最小反应c≈6.67%，最大反应d=100%，所以中间值是(6.67+100)/2≈53.33%，这就是你看到的y50数值。

• ED函数的默认行为：ED(D,50)默认计算的是相对效应的50%，即从基线c到最大值d的中间点，而不是绝对的50%死亡率。所以默认得到的ED50对应的Y值自然不是50。

正确绘制方法 #

1. 正确计算绝对50%死亡率对应的ED50 #

使用ED()函数的type="absolute"参数，指定计算绝对数值50对应的剂量：

# 计算绝对50%死亡率对应的ED50
ed50_abs <- ED(D, 50, type = "absolute", logBase = exp(1), interval = "delta")
# 提取ED50数值
ed50_val <- as.numeric(ed50_abs[1])
# 预测该剂量对应的Y值（应该接近50）
y_50_abs <- predict(D, newdata = data.frame(diuron = ed50_val))

2. 修正完整绘图代码 #

以下是调整后的代码，确保ED50线指向Y=50的位置，同时优化坐标轴设置：

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]

# 处理0剂量的对数转换问题
df$diuron0 <- df$diuron
df$diuron0[df$diuron0 == 0] <- 0.5

# 计算绝对50%死亡率对应的ED50
ed50_abs <- ED(D, 50, type = "absolute", logBase = exp(1), interval = "delta")
ed50_val <- as.numeric(ed50_abs[1])

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") +
  geom_ribbon(data = newdat, aes(x = diuron, ymin = pmin, ymax = pmax), 
              alpha = 0.2, fill = "#3C6255") +
  coord_trans(x = "log") +
  # 绘制ED50的垂直线和水平线，指向Y=50
  geom_segment(aes(x = ed50_val, y = 0, xend = ed50_val, yend = 50),
               lty = 2, colour = "gray50") +
  geom_segment(aes(x = 0.5, y = 50, xend = ed50_val, yend = 50),
               lty = 2, colour = "gray50") +
  # 优化x轴刻度，自动包含ED50值
  scale_x_log10(
    limits = c(100, 10000), 
    breaks = c(100, 500, 1000, round(ed50_val), 5000, 10000),
    labels = c(100, 500, 1000, round(ed50_val), 5000, 10000)
  ) +
  ggtitle("Diuron Dose-Response Curve (ED50 for 50% Mortality)") +
  xlab("Diuron (ug/L)") +
  ylab("Mortality (%)") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  theme_classic()

关键改动说明 #

• 用type="absolute"获取绝对50%死亡率对应的ED50剂量
• 调整线段绘制，让水平线直接指向Y=50，而不是模型中间值
• 自动将ED50值加入x轴刻度，无需手动硬编码
• 添加了置信区间的填充，让曲线更完整

内容的提问来源于stack exchange，提问作者Hayoung Kwon