标签: drc

我被这个问题困住了，我的 R 能力显然不足以解决它。希望有人可以帮助我。

我目前正在编写一个脚本，该脚本使用包“drc”的函数“drm”。我想从 drm 使用不同给定函数生成的最佳模型中获取 EC10、20 和 50 值。我想为不与所有或部分这些函数收敛的数据集实现一个解决方案（例如，如果不收敛，则打印出“数据集不收敛”）。使用 tryCatch 处理错误不起作用。它只捕获“drmOpt”产生的错误而不是“optim”产生的错误，从而停止脚本。
这是我在没有 try 或 tryCatch 的情况下得到的错误：

优化错误（startVec，opfct，hessian = TRUE，method = optMethod，control = list（maxit = maxIt，：非有限有限差分值[2]
drmOpt（opfct，opfct1，startVecSc，optMethod，约束，warnVal）中的错误, : 收敛失败

使用 try 或 tryCatch 我只会得到第一个错误。
可以使用以下代码重现该错误（LL.2 会产生错误，LL.3 不会）：

library(drc)
library(data.table)
data <- data.table(Dose = c(0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.5, 0.5,  0.5,  0.5,  3.0,  3.0, 3.0,  3.0, 9.0,  9.0,  9.0,  9.0, 27.0, 27.0, 27.0, 27.0, 81.0, 81.0, 81.0, 81.0), Value = c(1, 3, 2, 7, 5, 7, 6, …

在生物学中,我们经常想要绘制剂量反应曲线.R包'drc'非常有用,基本图形可以轻松处理'drm models'.但是,我想将我的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") …

我正在尝试重现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)  

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

编辑

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

我不确定为什么 ggplot2 中的 mean_sdl 函数（来自 Hmisc）生成的误差线比手动生成的误差线宽得多并绘制均值 + sd 和均值 - sd。我的代码：

library(drc)
library(tidyverse)

test_dataset <- 
  structure(
    list(
      X = c(1e-10, 1e-08, 3e-08, 1e-07, 3e-07, 1e-06, 3e-06, 1e-05, 3e-05, 1e-04, 3e-04),
      AY1 = c(0, 11, 125, 190, 258, 322, 354, 348, NA, 412, NA),
      AY2 = c(3, 33, 141, 218, 289, 353, 359, 298, NA, 378, NA),
      AY3 = c(2, 25, 160, 196, 345, 328, 369, 372, NA, 399, NA),
      BY1 = c(3, NA, 11, 52, 80, 171, 289, 272, 359, …

当我尝试使用 drc 包和 drm 函数将数据与非线性回归模型拟合时，我对这条警告消息感到困惑。

我有

 N_obs <- c(1, 80, 80, 80, 81, 82, 83, 84, 84, 95, 102, 102, 102, 103, 104, 105, 105, 109, 111, 117, 120, 123, 123, 124, 126, 127, 128, 128, 129, 130)

 times <- c(3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32)

该模型是

  model.drm <- drm(N_obs ~ times, data = data.frame(N_obs = …

我想绘制相同数据的线性模型(LM)和非线性(GLM)模型.

LM和GLM之间的范围应在16%-84%之间,引用:第3.5节

我已经包含了更完整的代码块,因为我不确定在哪一点上我应该尝试削减线性模型.或者我搞砸了 - 我想用线性模型.

下面的代码会产生以下图像:

我的目标(摘自以前的引文链接).

这是我的数据:

mydata3 <- structure(list(
               dose = c(0, 0, 0, 3, 3, 3, 7.5, 7.5, 7.5, 10,     10, 10, 25, 25, 25, 50, 50, 50), 
               total = c(25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L), 
               affected = c(1, 0, 1.2, 2.8, 4.8, 9, 2.8, 12.8, 8.6, 4.8, 4.4, 10.2, 6, 20, 14, 12.8, 23.4, 21.6), 
               probability = c(0.04, 0, 0.048, 0.112, …

我正在绘制具有渐近尾部的剂量 - 反应曲线.我真的想在显示0的图中包括载体(对照)剂量

0通常计算为.0000000001的剂量 - 这些图中的常见做法.

我真的很喜欢下面的图像显示这个,图像是从pdf中获取的,如何使用该程序绘图:GraphPad:PRISM

旁注:我已经找到了如何使用基本图形,但不使用ggplot2.

关于matlab提出了类似但不同的SE问题:这里

我的R代码如下:

library(ggplot2)
library(scales)

ggplot(df, aes(x=dose,y=probability, group=model))+

  geom_ribbon(aes(ymin=Lower,ymax=Upper,x=dose,
                  fill=model, col=model,alpha=2))+

  #Axis log transformation:
  annotation_logticks(scaled = TRUE,sides="b") +
  scale_x_log10(breaks = 10^(-1:10),
                labels = trans_format("log10", math_format(10^.x)))+
  #Axes labels:
  labs(x="dosage (log scale)", y="response",size=1)

数据:

df<-structure(list(dose = c(1.0000001, 1.04737100217022, 1.09698590648847,1.14895111335032, 1.20337795869652, 1.26038305255123, 1.32008852886009,1.38262230716338, 1.44811836666478, 1.51671703328309, 61.5098612858473,64.4236386159454, 67.4754441930906, 70.6718165392165, 74.0196039119089,77.525978976861, 81.1984541753771, 85.0448978198478, 89.073550951683,93.2930449978201, 97.7124202636365, 102.341145301888, 107.189137199173,112.266782823381, 117.584961077656, 123.155066208544, 128.98903221828,135.099358433491, 141.499136285126, 148.202077356965, 155.222542762819,162.575573915347, 6294.98902185499, 6593.18830115198, 6905.51354792318,7232.63392192496, 7575.25028161186, 7934.09668573241, 8309.9419660568,8703.59137460616, 9115.88830891252, 9547.71611900591, 10000),probability = c(0.000541224108467882, …