drc:: ggplot2 的 drc 绘图答案

【问题标题】：drc:: drc plot with ggplot2drc:: ggplot2 的 drc 绘图
【发布时间】：2016-11-05 07:16:49
【问题描述】：

我正在尝试用ggplot2 重现drc 的情节。这是我的第一次尝试（下面给出了 MWE）。但是，我的 ggplot2 与基本 R 图有点不同。我想知道我是否在这里遗漏了什么？

library(drc)
chickweed.m1 <- drm(count~start+end, data = chickweed, fct = LL.3(), type = "event")

plot(chickweed.m1, xlab = "Time (hours)", ylab = "Proportion germinated", 
xlim=c(0, 340), ylim=c(0, 0.25), log="", lwd=2, cex=1.2)

library(data.table)
dt1 <- data.table(chickweed)

dt1Means1 <- dt1[, .(Germinated=mean(count)/200), by=.(start)]
dt1Means2 <- dt1Means1[, .(start=start, Germinated=cumsum(Germinated))]
dt1Means  <- data.table(dt1Means2[start!=0], Pred=predict(object=chickweed.m1))

library(ggplot2)
ggplot(data= dt1Means, mapping=aes(x=start, y=Germinated)) + 
    geom_point() +
    geom_line(aes(y = Pred)) +
    lims(y=c(0, 0.25)) +
    theme_bw()

已编辑

我遵循here 给出的方法（做了一些更改）。

【问题讨论】：

您可以随时查看getAnywhere(plot.drc) 了解如何在基础图中计算数据。
你想在 ggplot 中改变什么？两个图表中的数据似乎相同。
@MYaseen208：所以给他赏金吧。保持它打开没有任何意义。
@AndrewBrēza：请看dww的回答。
当然@42-，赏金是给dww的。

标签： r ggplot2 drc

【解决方案1】：

注意，您可以跳到最后一段以获得简单的答案。该答案的其余部分记录了我是如何得出该解决方案的

查看drc:::plot.drc的代码，我们可以看到最后一行无形地返回了一个data.frame retData

function (x, ..., add = FALSE, level = NULL, type = c("average", 
                                                      "all", "bars", "none", "obs", "confidence"), broken = FALSE, 
          bp, bcontrol = NULL, conName = NULL, axes = TRUE, gridsize = 100, 
          log = "x", xtsty, xttrim = TRUE, xt = NULL, xtlab = NULL, 
          xlab, xlim, yt = NULL, ytlab = NULL, ylab, ylim, cex, cex.axis = 1, 
          col = FALSE, lty, pch, legend, legendText, legendPos, cex.legend = 1, 
          normal = FALSE, normRef = 1, confidence.level = 0.95) 
{
  # ...lot of lines omitted...
  invisible(retData)
}

retData 包含拟合模型线的坐标，因此我们可以使用它来对 plot.drc 使用的相同模型进行 ggplot

pl <- plot(chickweed.m1, xlab = "Time (hours)", ylab = "Proportion germinated", 
        xlim=c(0, 340), ylim=c(0, 0.25), log="", lwd=2, cex=1.2)
names(pl) <- c("x", "y")
ggplot(data= dt1Means, mapping=aes(x=start, y=Germinated)) + 
  geom_point() +
  geom_line(data=pl, aes(x=x, y = y)) +
  lims(y=c(0, 0.25)) +
  theme_bw()

这与您在 ggplot 中使用 predict(object=chickweed.m1) 创建的版本相同。因此，差异不在于模型线，而在于绘制数据点的位置。我们可以通过将函数的最后一行从invisible(retData) 更改为list(retData, plotPoints)，从drc:::plot.drc 导出数据点。为方便起见，我将 drc:::plot.drc 的整个代码复制到了一个新函数中。请注意，如果您希望复制此步骤，drcplot 调用的一些函数未在 drc 命名空间中导出，因此需要在所有对函数 parFct、addAxes、@ 的调用之前添加 drc::: 987654334@、makeLegend。

drcplot <- function (x, ..., add = FALSE, level = NULL, type = c("average", 
                                                      "all", "bars", "none", "obs", "confidence"), broken = FALSE, 
          bp, bcontrol = NULL, conName = NULL, axes = TRUE, gridsize = 100, 
          log = "x", xtsty, xttrim = TRUE, xt = NULL, xtlab = NULL, 
          xlab, xlim, yt = NULL, ytlab = NULL, ylab, ylim, cex, cex.axis = 1, 
          col = FALSE, lty, pch, legend, legendText, legendPos, cex.legend = 1, 
          normal = FALSE, normRef = 1, confidence.level = 0.95) 
{
  # ...lot of lines omitted...
  list(retData, plotPoints)
}

用你的数据运行它

pl <- drcplot(chickweed.m1, xlab = "Time (hours)", ylab = "Proportion germinated", 
          xlim=c(0, 340), ylim=c(0, 0.25), log="", lwd=2, cex=1.2)

germ.points <- as.data.frame(pl[[2]])
drc.fit <- as.data.frame(pl[[1]])
names(germ.points) <- c("x", "y")
names(drc.fit) <- c("x", "y")

现在，用 ggplot2 绘制这些得到你想要的东西

ggplot(data= dt1Means, mapping=aes(x=start, y=Germinated)) + 
  geom_point(data=germ.points, aes(x=x, y = y)) +
  geom_line(data=drc.fit, aes(x=x, y = y)) +
  lims(y=c(0, 0.25)) +
  theme_bw()

最后，将此图 (germ.points) 的数据点值与原始 ggplot (dt1Means) 中的数据点值进行比较，显示了差异的原因。 dt1Means 中的计算点相对于 plot.drc 中的点提前了一个时间段。换句话说，plot.drc 将事件分配给它们发生的时间段的结束时间，而您将发芽事件分配给它们发生的时间间隔的开始时间。您可以简单地调整它，例如，使用

dt1 <- data.table(chickweed)
dt1[, Germinated := mean(count)/200, by=start]
dt1[, cum_Germinated := cumsum(Germinated)]
dt1[, Pred := c(predict(object=chickweed.m1), NA)]  # Note that the final time period which ends at `Inf` can not be predicted by the model, therefore added `NA` in the final row

ggplot(data= dt1, mapping=aes(x=end, y=cum_Germinated)) + 
  geom_point() +
  geom_line(aes(y = Pred)) +
  lims(y=c(0, 0.25)) +
  theme_bw()

【讨论】：

感谢@dww 非常有帮助的回答。对于您的回答，据我所知，dt1Means2 <- dt1Means1[, .(start=start, Germinated=c(0,cumsum(Germinated)))] 就足够了。但是，它没有给我所需的输出。请有任何想法。
只使用dt1Means2 <- dt1Means1[, .(start=start, Germinated=c(0,cumsum(Germinated)))] 不是解决方案。它给出以下消息In as.data.table.list(jval) : Item 1 is of size 35 but maximum size is 36 (recycled leaving a remainder of 1 items)。
@MYaseen208，该消息是警告，而不是错误。您可以放心地忽略它。也就是说，这只是我可以对您的代码进行的最小更改以获得您想要的正确绘图 - 但不是最整洁的方式。更好的是使用结束时间而不是开始时间正确构建您的数据（这是我在答案中已经建议的）。我已经编辑以展示如何做到这一点。
当curveid = Temp 时，我正在努力如何获得ggplot2 版本。如果您查看此question，将不胜感激。谢谢

【解决方案2】：

从@dww 的回答中得到直觉，我不得不对我的原始代码做两个小改动。只需在

中将start!=0 替换为end!=Inf

dt1Means1 <- dt1[, .(Germinated=mean(count)/200), by=.(start, end)]
dt1Means  <- data.table(dt1Means2[start!=0], Pred=predict(object=chickweed.m1))

给出正确的图表。

【讨论】：

【解决方案3】：

我真的很喜欢 dww 提供的解决方案。我是否可以建议对这个解决方案进行概括。通过将以下几行添加到drc:::plotdrc() 的自写版本中，您可以概括解决方案。该函数接受drc:::plotdrc() 函数的输入，但输出一个与原始函数的默认基本图输出具有相同规格的ggplot-object。

只需将invisible(retData, plotPoints) 替换为

result <- list(retData, plotPoints) 
points <- as.data.frame(result[[2]])
drc.fit <- as.data.frame(result[[1]]) 
names(points) <- c("x", "y")
names(drc.fit) <- c("x", "y")` 

gg_plot <- ggplot2::ggplot(data=points, aes(x = x, y = y)) + 
geom_point() +
geom_line(data=drc.fit, aes(x = x, y = y)) +
scale_x_continuous(trans='log10', limits = xlim) +
ylab(ylab) +
xlab(xlab) +
lims(y = ylim) +
theme_bw()
return(gg_plot)`

【讨论】：