drc :: drc plot with ggplot2

MYa*_*208 5 r ggplot2 drc

我正在尝试重现drc情节ggplot2.这是我的第一次尝试(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)  
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述

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()
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述

编辑

我遵循这里给出的方法(有一些变化).

dww*_*dww 8

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

查看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)
}
Run Code Online (Sandbox Code Playgroud)

retData包含拟合模型线的坐标,因此我们可以使用它来绘制plot.drc使用的相同模型

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()
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述

这与您使用predict(object = chickweed.m1)在ggplot中创建的版本相同.因此,差异不在于模型线,而在于绘制数据点的位置.我们可以通过改变功能的最后一行输出从刚果(金)::: plot.drc数据点invisible(retData)list(retData, plotPoints).为方便起见,我将drc ::: plot.drc的整个代码复制到一个新函数中.需要注意的是,如果你想复制这一步,还有由未在刚果(金)的命名空间导出drcplot叫了几个功能,所以drc:::需要预先考虑到的函数的调用parFct,addAxes,brokenAxis,和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)
}
Run Code Online (Sandbox Code Playgroud)

并使用您的数据运行它

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")
Run Code Online (Sandbox Code Playgroud)

现在,使用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()
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述

最后,将此plot(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()
Run Code Online (Sandbox Code Playgroud)