首页 文章

drc :: drc用ggplot2绘图

提问于
浏览
5

我正在尝试使用 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)

enter image description here

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

enter image description here

Edited

我按照here的方法(有一些更改) .

3 回答

  • 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)
    }
    

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

    enter image description here

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

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

    enter image description here

    最后,将此图( 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()
    
  • 2

    从@dww的答案获得直觉,我不得不在原始代码中进行两处小改动 . 只需将 start!=0 替换为 end!=Inf in

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

    给出正确的图表 .

  • 1

    我真的很喜欢dww提出的解决方案 . 我可能会建议对此解决方案进行推广 . 通过将以下行添加到自编写的 drc:::plotdrc() 版本中,您可以概括解决方案 . 该函数接受 drc:::plotdrc() 函数的输入,但输出具有相同规格的ggplot对象作为原始函数的默认基本图输出 .

    只需将 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)`
    

相关问题