有没有办法在一个图上绘制具有不同衰减常数的指数衰减曲线？

问题描述

我一直试图在一张图上绘制不同的指数衰减曲线。最初我认为这会很容易，但结果却令人沮丧。

我想得到什么：

Picture of what I am trying to get.

nlsplot(k_data_nls,model = 6,start = c(a= 603.3,b= -0.03812),xlab = "hours",ylab = "copies")

nlsplot(r4,model=6,start=c(a=25.5487,b=-0.5723),ylab = "copies")

这里是一些额外的数据代码：

df4 <- data.frame(hours=c(0,1,3,5,12,24,48,96,168,336,504,720),copies=c(603.3,406,588,393.27,458.47,501.67,767.53,444.13,340.6,298.47,61.42,51.6))
nlsfit(df4,start=c(a=603.3,b=-0.009955831526))
d4plot <- nlsplot(df4,b=-0.009955831526))

r4 <- data.frame(hours=c(0,copies=c(26,13.44,4.57,3.12,6.89,0.71,0.47,0.24,0.48))
nlsLM(copies ~ a*exp(b*hours),data=r4,start=list(a=26,b=-0.65986))
r4plot <- nlsplot(r4,b=-0.5723))

从本质上讲，我希望能够在一张图上获得这两个图。我是 R 的新手，所以我不太确定我可以从哪里开始。谢谢！

解决方法

我不知道这是否真的有用，因为它太具体了，但我会这样做（使用 ggplot2）。首先，您需要要绘制的函数的数据。取 x 代表您要显示的所有值，并将您的函数和系数应用于数据。您需要有数据点，而不仅仅是一个函数来绘制数据。

df_simulated <- data.frame("x" = rep(1:100,2),"class"= rep(c("DNA","RNA"),each = 100))
df_simulated$y <- c(1683.7 * exp(-0.103 * 1:100),# DNA
                    578.7455 * exp(-0.156 * 1:100))  # RNA

但是，由于我从未使用过您使用的软件包，因此我不知道如何从模型中提取值，因此我采用了您的示例图中的值。重要的是，两组的“模拟”值都在一个数据框中，并且您有一列将点归因于相应的组（RNA 或 DNA）。如果你这样做，至少它会更容易。然后您需要一个数据框，其中包含您对点的实际观察结果。我又发明了数据：

df_observed <- data.frame("x" = c(12,13,25,26,50,51),"y" = c(500,250,5),"class" = rep(c("DNA",3))

然后您可以创建绘图。使用 color=class 指定数据点将按“类”分组并相应地着色。（“apple”和“banana”只是用来演示换行的虚拟词）

ggplot() +
  geom_line(data = df_simulated,aes(x = x,y = y,color = class),size = 1,linetype = "dashed") +
  geom_point(data = df_observed,size = 4,pch = 1) +
  annotate("text",x = 50,y = 1250,label = "DNA\napple",color = "tomato",hjust = 0) +
  annotate("text",y = 750,label  ="RNA\nbanana",color = "steelblue",hjust = 0) +
  ggtitle(expression(~italic("Styela clava")~"(isolated)")) +
  ylab("COI copies per 1ml") +
  xlab("Time since removal of organisms (hours)") +
  theme_classic() +
  theme(legend.position = "none") +
  scale_color_manual(values = c("DNA" = "tomato","RNA" = "steelblue"))

这是输出：

首先要注意，R平方通常用于线性模型而不是非线性模型，所以the use of this statistic is suspect here；然而，下面我们无论如何都会展示它，因为这似乎是我们所要求的。经常使用的另一种拟合优度测量是残差标准误差。如果 fm 是来自 nls 的拟合模型，则 sigma(fm) 是残差标准误差。较小的值更有利。 summary(fm) 也报告这个值。

对于 df4 和 r4 中的每一个，我们使用 lm 来获取起始值（对两边取对数，我们得到一个在 log(a) 和 b 中呈线性的模型），运行 nls 拟合并得到系数。

现在绘制点并添加拟合线和图例。（请注意，在设置图表时，我们使用 rbind，它假定 df4 和 r4 具有相同的列名，它们确实如此。）

请注意，问题中提供的数据与问题图片中显示的数据大不相同。

下面的代码不需要起始值，因为它使用 lm 来获取它们，运行 nls 并自动提取图形所需的任何信息。

1) 经典图形 在这个替代方案中，没有使用任何包。

r2 <- function(fm,digits = 3) {
  y <- fitted(fm) + resid(fm) 
  r2 <- 1 - deviance(fm) / sum((y - mean(y))^2)
  if (is.numeric(digits)) r2 <- round(r2,digits)
  r2
}

fo <- copies ~ a * exp(b * hours)  # formula used in nls

# get nls fitted model and coefficients for df4
co_d0 <- coef(lm(log(copies) ~ hours,df4,subset = copies > 0))
fmd <- nls(fo,start = list(a = exp(co_d0[[1]]),b = co_d0[[2]]))
co_d <- round(coef(fmd),4)

# get nls fitted model and coefficients for r4
co_r0 <- coef(lm(log(copies) ~ hours,r4,subset = copies > 0))
fmr <- nls(fo,start = list(a = exp(co_r0[[1]]),b = co_r0[[2]]))
co_r <- round(coef(fmr),4)

both <- rbind(cbind(df4,col = "red"),cbind(r4,col = "blue"))
plot(both[1:2],col = both$col,xlab = "Time since removal of organisms",ylab = "COI copies per 1ml",main = "C)" ~ italic("Styela clava") ~ "(isolated)",adj = 0)
lines(fitted(fmd) ~ hours,col = "red",lty = 2)
lines(fitted(fmr) ~ hours,col = "blue",lty = 2)

legend <- c(bquote(DNA),bquote(y == .(co_d[[1]]) * e ^ {.(co_d[[2]])*x}),bquote(R^2 == .(r2(fmd))),bquote(),bquote(RNA),bquote(y == .(co_r[[1]]) * e ^ {.(co_r[[2]])*x}),bquote(R^2 == .(r2(fmr))))
legend("right",legend = as.expression(legend),bty = "n",text.col = c("red","red",NA,"blue","blue"))

2) ggplot2 这使用了 ggplot2 和 gridtext。 r2、fmd、fmr、co_d 和 co_r 均取自 (1)。我们使用 gridtext 中的 richtest_grob 为图例创建自定义 grob，并使用 annotate_custom 传递它。

library(gridtext)
library(ggplot2)

txt <- sprintf(
  "<span style='color:red'>DNA
    <br>y = %.3f*e<sup>%.3fx</sup> 
    <br>R<sup>2</sup> = %.3f</span>
    <br><br><span style='color:blue'>RNA
    <br>y = %.3f*e<sup>%.3fx</sup>
    <br>R<sup>2</sup> = %.3f</span>",co_d[[1]],co_d[[2]],r2(fmd),co_r[[1]],co_r[[2]],r2(fmr))

both2 <- rbind(cbind(df4,fitted = fitted(fmd)),fitted = fitted(fmr)))
ggplot(both2,aes(hours,copies,col = I(col))) +
  geom_point() +
  geom_line(aes(y = fitted),linetype = 2) +
  annotation_custom(richtext_grob(txt,hjust = 0)) +
  theme(legend.position = "none") +
  labs(x = "Time since removal of organisms",y = "COI copies per 1ml") +
  ggtitle(("C)" ~ italic("Styela clava") ~ "(isolated)"))

3) 点阵

这使用了 (1) 中的图例和 (2) 中的 both2。首先为数据点创建一个图。它还将包含图例、轴和标签。然后为拟合线添加一个图层。 main.settings 指定主标题应左对齐和粗体，并改编自 this page。

library(latticeExtra)

main.settings <- list(par.main.text = list(font = 2,just = "left",x = grid::unit(25,"mm")))

xyplot(copies ~ hours,both2,col = both2$col,adj = 0,key = list(text = list(as.expression(legend),col = c("red","blue")),x = 0.65,y = 0.65,columns = 1),par.settings = main.settings) +
  as.layer(xyplot(fitted ~ hours,groups = col,type = "l",lty = 2))

graphics graphics nls r r