如何优化R中drc与ggplot2绘制的剂量-反应曲线？

问题描述 #

作为剂量-反应建模新手，使用drc和ggplot2绘制曲线时遇到两个问题：

• 数据点与坐标轴之间存在间隙
• x轴数值聚集在末端，无法均匀间隔分布

原始数据 #

structure(list(nickel = c(0.5, 0.5, 0.5, 0.5, 50, 50, 75, 75, 
                          113, 113, 169, 169, 253, 253, 380, 380, 570, 570, 854, 854, 1281, 
                          1281, 1922, 1922, 2883, 2883, 4325, 4325), growth = c(100L, 104L, 
                          101L, 93L, 101L, 100L, 94L, 75L, 74L, 71L, 60L, 85L, 72L, 52L, 
                          53L, 71L, 65L, 43L, 55L, 49L, 38L, 40L, 50L, 35L, 20L, 0L, 0L, 
                          17L)), row.names = c(NA, -28L), class = "data.frame")

原始代码 #

Ni.m1 <- drm(dose_response.data$growth~dose_response.data$nickel, data = dose_response.data, fct = LL.3())

newdata <- expand.grid(conc=exp(seq(log(0.1), log(6000), length=200)))
pm <- predict(Ni.m1, newdata=newdata, interval="prediction")

newdata$p <- pm[,1]
newdata$pmin <- pm[,2]
newdata$pmax <- pm[,3]

dose_response.data$nickel <- dose_response.data$nickel
dose_response.data$nickel[dose_response.data$nickel == 0] <- 0.5

ggplot(dose_response.data, aes(x = nickel, y = growth))+
  geom_point(pch=1,cex=2) +
  geom_ribbon(data=newdata, aes(x=conc, y=p, ymin=pmin, ymax=pmax), alpha=0.2) +
  geom_line(data=newdata, aes(x=conc, y=p)) +
  coord_trans(x="log") +
  xlab("Nickel"~(µg.L^{-1})) + 
  ylab("Growth (% control)") +
  ylim (-20,130) +
#manually add lines representing the EC10 concentration
  geom_segment(aes(x = 42.5, y = 0, xend = 42.5, yend = 90), linetype = 2, lwd=1) +
  geom_segment(aes(x = 0.1, y = 90, xend = 42.5, yend = 90), linetype = 2, lwd=1) 

原始绘图结果 #

解决方案 #

1. 消除数据与坐标轴的间隙 #

coord_trans(x="log")是先绘制图形再对坐标轴做对数转换，容易导致坐标轴范围和实际数据范围错位，产生间隙。改用scale_x_log10()，它会先对x轴数据做对数转换再绘制，同时设置expand = c(0, 0)可以彻底消除间隙。

2. 让x轴刻度均匀分布 #

在scale_x_log10()中使用breaks参数指定对数刻度，比如用scales::log_breaks()自动生成均匀分布的对数刻度，或者手动指定刻度值，确保刻度在对数轴上均匀分布。

修改后的完整代码 #

library(drc)
library(ggplot2)
library(scales)

# 数据处理
dose_response.data <- structure(list(nickel = c(0.5, 0.5, 0.5, 0.5, 50, 50, 75, 75, 
                          113, 113, 169, 169, 253, 253, 380, 380, 570, 570, 854, 854, 1281, 
                          1281, 1922, 1922, 2883, 2883, 4325, 4325), growth = c(100L, 104L, 
                          101L, 93L, 101L, 100L, 94L, 75L, 74L, 71L, 60L, 85L, 72L, 52L, 
                          53L, 71L, 65L, 43L, 55L, 49L, 38L, 40L, 50L, 35L, 20L, 0L, 0L, 
                          17L)), row.names = c(NA, -28L), class = "data.frame")

Ni.m1 <- drm(growth ~ nickel, data = dose_response.data, fct = LL.3())

newdata <- expand.grid(conc=exp(seq(log(0.1), log(6000), length=200)))
pm <- predict(Ni.m1, newdata=newdata, interval="prediction")

newdata$p <- pm[,1]
newdata$pmin <- pm[,2]
newdata$pmax <- pm[,3]

dose_response.data$nickel[dose_response.data$nickel == 0] <- 0.5

# ggplot绘图
ggplot(dose_response.data, aes(x = nickel, y = growth))+
  geom_point(pch=1, cex=2) +
  geom_ribbon(data=newdata, aes(x=conc, y=p, ymin=pmin, ymax=pmax), alpha=0.2) +
  geom_line(data=newdata, aes(x=conc, y=p)) +
  # 替换coord_trans为scale_x_log10，设置刻度和消除间隙
  scale_x_log10(
    breaks = log_breaks(n = 8),  # 自动生成8个均匀分布的对数刻度
    expand = c(0, 0),            # 消除坐标轴间隙
    labels = label_number()      # 保持刻度为常规数字格式
  ) +
  xlab(expression(Nickel ~ (µg.L^{-1}))) + 
  ylab("Growth (% control)") +
  ylim (-20,130) +
  # EC10标注线
  geom_segment(aes(x = 42.5, y = 0, xend = 42.5, yend = 90), linetype = 2, lwd=1) +
  geom_segment(aes(x = 0.1, y = 90, xend = 42.5, yend = 90), linetype = 2, lwd=1) +
  theme_bw()

关键修改说明 #

• 替换coord_trans(x="log")为scale_x_log10()，从根本上解决对数转换导致的坐标轴错位问题
• expand = c(0, 0)：移除x轴两侧的空白间隙，让数据贴合坐标轴
• breaks = log_breaks(n = 8)：自动生成均匀分布的对数刻度，避免刻度聚集在末端
• 优化了drm函数的公式写法，更符合R语法规范

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