8

在生物学中，我们经常想要绘制剂量反应曲线。R 包“drc”非常有用，基础图形可以轻松处理“drm 模型”。但是，我想将我的 drm 曲线添加到 ggplot2。

我的数据集：

 library("drc")
 library("reshape2")
 library("ggplot2")
 demo=structure(list(X = c(0, 1e-08, 3e-08, 1e-07, 3e-07, 1e-06, 3e-06, 
 1e-05, 3e-05, 1e-04, 3e-04), Y1 = c(0, 1, 12, 19, 28, 32, 35, 
 39, NA, 39, NA), Y2 = c(0, 0, 10, 18, 30, 35, 41, 43, NA, 43, 
 NA), Y3 = c(0, 4, 15, 22, 28, 35, 38, 44, NA, 44, NA)), .Names = c("X", 
"Y1", "Y2", "Y3"), class = "data.frame", row.names = c(NA, -11L
))

使用基本图形：

plot(drm(data = reshape2::melt(demo,id.vars = "X"),value~X,fct=LL.4(),na.action = na.omit),type="bars")

产生一个很好的 4 参数剂量反应图。

试图在 ggplot2 中绘制相同的情节，我偶然发现了 2 个问题。

  1. 没有办法直接添加drm模型曲线。我需要将 4-PL 重写为一个函数，并以 stat_function 的形式添加，至少可以说很麻烦。

    ggplot(reshape2::melt(demo,id.vars = "X"),aes(X,value)) + 
  geom_point() + 
  stat_function(fun = function(x){
    drm_y=function(x, drm){
      coef(drm)[2]+((coef(drm)[3]-coef(drm)[2])/(1+exp((coef(drm)[1]*(log(x)-log(coef(drm)[4]))))))
    }
+ drm_y(x,drm = drm(data = reshape2::melt(demo,id.vars = "X"), value~X, fct=LL.4(), na.action = na.omit))
 })

  2. 如果这还不够，它仅在 scale_x 是连续的情况下才有效。如果我想添加scale_x_log10()，我得到： Warning message: In log(x): NaNs produced。

我意识到这一点，log10(0) = -Inf但有一些方法可以解决这个问题。无论是（与 plot.drc 的情况一样），x=0 值都绘制在 x 轴上，基本上是之前最低 x 值的 1/100。( demo$X[which.min(demo$X)+1]/100) 或在 GraphPad Prism 中，0 完全从剂量反应曲线中省略。

我的问题是：

  1. 有没有办法直接在 ggplot2 中绘制 drm 模型？

  2. 如何将数据集与其相应的 4-PL 曲线拟合链接起来，以便它们以相同的颜色绘制？

4

3 回答 3

17

该软件包的作者最近发表的一篇论文drc包括提取参数以供 ggplot2 使用的说明。它们不在 ggplot2 中工作，而是从模型中提取数据。这是他们应用于您的数据的解决方案。

demo1 <- reshape2::melt(demo,id.vars = "X") # get numbers ready for use.
demo.LL.4 <- drm(data = demo1,value~X,fct=LL.4(),na.action = na.omit) # run model.

该函数可以从模型predict中提取参数。drm它与使用curveid.

# predictions and confidence intervals.
demo.fits <- expand.grid(conc=exp(seq(log(1.00e-04), log(1.00e-09), length=100))) 
# new data with predictions
pm <- predict(demo.LL.4, newdata=demo.fits, interval="confidence") 
    demo.fits$p <- pm[,1]
    demo.fits$pmin <- pm[,2]
    demo.fits$pmax <- pm[,3]

他们建议改变零浓度以避免 coord_trans 出现问题。

demo1$XX <- demo1$X
demo1$XX[demo1$XX == 0] <- 1.00e-09

然后是绘制曲线，省略会geom_ribbon阻止绘制错误。

ggplot(demo1, aes(x = XX, y = value)) +
  geom_point() +
  geom_ribbon(data=demo.fits, aes(x=conc, y=p, ymin=pmin, ymax=pmax), alpha=0.2) +
  geom_line(data=demo.fits, aes(x=conc, y=p)) +
  coord_trans(x="log")

在此处输入图像描述

为了将多条曲线绘制在一起，可以重复该过程。为每个集合添加 ID。

demo.fits_1 <- data.frame(label = "curve1", demo.fits)

然后使用rbind来组合所有提取的参数。从那里 ggplot 可以处理颜色。

于 2016-05-05T06:21:28.437 回答
4

我将回答我自己的问题，希望这将帮助其他面临同样问题的人。

当然，可以使用 ggplot2 和 drc 包绘制剂量反应曲线，只需添加 geom_ 或stat_smooth (method=drm, fct=LL.4(),se=FALSE)在线性刻度或 geom_ 上绘制或stat_smooth (method=drm, fct=L.4(),se=FALSE)添加scale_x_log10()。

为了能够使用 log10 比例，我将数据转换为：

demo <- demo %>% 
      mutate(X = 
       ifelse(X == 0, 
              yes = (sort(demo$X[which.min(sort(demo$X)) + 1]/100)),
              no = X
              )
            )         #looks for the pre-lowest value in X and divides it by 100

在本例中，我已将 X = 0 值替换为 X = 上一个 X 值的 1/100（在本例中为 1e-10）。但是，您可以像 Prism 一样，通过从数据集中完全省略它来轻松删除弄乱对数绘图的 0 值。正如我发现的那样，需要注意的一件事是 ggplot 首先缩放轴然后添加数据，这就是代码在尝试 log10(0) 时中断的原因。

另一个微妙之处是 stat_smooth 函数完全能够使用处理 drm 模型，method = drm但它不知道如何拟合“SE”置信区间。因此选择se = FALSE可以进行绘图，并且以我的拙见，无论如何都会使绘图变得不那么混乱 - 只需添加误差线。

最后，将 更改fct = LL.4()为fct = L.4()允许在 log10 比例上绘图，因为再次首先选择比例，然后进行拟合。因此，即使轴值是非对数的，ggplot 实际上已将数据集转换为 log10，因此拟合函数现在只需 logit-4P（即 L.4()）而不是 log-logit-4P（LL .4())。

geom_smooth() 和 stat_smooth() 函数自然会采用与数据集相同的颜色，无需调整拟合函数的颜色以对应数据点的颜色。

总之：

demo <- demo %>% 
      mutate(X = 
       ifelse(X == 0, 
              yes = (sort(demo$X[which.min(sort(demo$X)) + 1]/100)),
              no = X
              )
            )
demo.long <- reshape2::melt(demo,id.vars = "X") #reshapes the demo dataset to long format
ggplot(data = demo.long,
       aes(x = X, y = value, col = variable)
      ) + 
   geom_point() + 
   geom_smooth(method = drm, fct = L.4(), se = FALSE) +
   scale_x_log10() #plots out the dataset with the corresponding 4-parameter log-logit dose response curves
于 2016-04-29T20:34:31.637 回答
0

更新的答案：geom_smooth(method = drm, method.args = list(fct = L.4()), se = FALSE)非常有帮助！

于 2022-01-07T19:27:07.457 回答