用ggplot将S形曲线拟合到点

问题描述

我有一个简单的数据框，用于测量不同剂量药物治疗的反应：

drug <- c("drug_1","drug_1","drug_2","drug_2")

conc <- c(100.00,33.33,11.11,3.70,1.23,0.41,0.14,0.05,100.00,0.05)

mean_response <- c(1156,1833,1744,1256,1244,1088,678,489,2322,1867,1333,944,567,356,200,177)

std_dev <- c(117,317,440,134,38,183,153,719,218,185,117,166,167,88,50)

df <- data.frame(drug,conc,mean_response,std_dev)

我可以使用以下代码绘制这些点，并获得我想要的可视化的基本基础：

p <- ggplot(data=df,aes(y=mean_response,x= conc,color = drug)) +
  geom_pointrange(aes(ymax = (mean_response + std_dev),ymin = (mean_response - std_dev))) +
  scale_x_log10()

p

plot

接下来我想对这些数据进行的操作是在曲线上添加一个S型曲线，以适合每种药物的曲线点。之后，我想计算该曲线的EC50。我意识到我的数据可能没有S形曲线的整个范围，但我希望对已有的数据能得到最好的估计。同样，drug_1的终点未遵循S形曲线的预期趋势，但这实际上并不出乎意料，因为该药物所在的溶液可以抑制高浓度时的反应（每种药物在不同的溶液中）。我想从数据中排除这一点。

在将S形曲线拟合到我的数据时，我陷入了困境。我看过其他一些解决方案，可将S型曲线拟合到数据中，但似乎无济于事。

与我的问题非常接近的一篇文章是： (sigmoid) curve fitting glm in r

基于此，我尝试了：

p + geom_smooth(method = "glm",family = binomial,se = FALSE)

这会产生以下错误，并且似乎默认为绘制直线：

`geom_smooth()` using formula 'y ~ x'
Warning message:
Ignoring unkNown parameters: family

我也尝试过此链接中的解决方案： Fitting a sigmoidal curve to this oxy-Hb data

在这种情况下，出现以下错误：

computation Failed in `stat_smooth()`:
Convergence failure: singular convergence (7)

并且没有线添加到绘图中。

我尝试查找这两个错误，但似乎找不到适合我的数据的原因。

任何帮助将不胜感激！

解决方法

正如我在评论中所说，对于一个非常简单的问题，我只会使用geom_smooth()；一旦遇到麻烦，我就会改用nls。

我的答案与@Duck的答案非常相似，但有以下区别：

我同时显示了未加权和（反方差）加权拟合。
为了使加权拟合工作，我不得不使用nls2软件包，该软件包提供了更健壮的算法
我使用SSlogis()进行自动（自动启动）初始参数选择
我在ggplot2之外进行所有预测，然后将其输入geom_line()

p1 <- nls(mean_response~SSlogis(conc,Asym,xmid,scal),data=df,subset=(drug=="drug_1" & conc<100)
        ##,weights=1/std_dev^2  ## error in qr.default: NA/NaN/Inf ...
          )

library(nls2)
p1B <- nls2(mean_response~SSlogis(conc,subset=(drug=="drug_1" & conc<100),weights=1/std_dev^2)

p2 <- update(p1,subset=(drug=="drug_2"))
p2B <- update(p1B,subset=(drug=="drug_2"))

pframe0 <- data.frame(conc=10^seq(log10(min(df$conc)),log10(max(df$conc)),length.out=100))
pp <- rbind(
    data.frame(pframe0,mean_response=predict(p1,pframe0),drug="drug_1",wts=FALSE),data.frame(pframe0,mean_response=predict(p2,drug="drug_2",mean_response=predict(p1B,wts=TRUE),mean_response=predict(p2B,wts=TRUE)
)

library(ggplot2); theme_set(theme_bw())
(ggplot(df,aes(conc,mean_response,colour=drug)) +
 geom_pointrange(aes(ymin=mean_response-std_dev,ymax=mean_response+std_dev)) +
 scale_x_log10() +
 geom_line(data=pp,aes(linetype=wts),size=2)
)

我相信EC50等效于xmid参数...请注意加权和未加权估计之间的巨大差异...

我建议采用与您想要的接近的下一种方法。我也尝试过使用binomial系列为您的数据进行设置，但是关于0到1之间的值存在一些问题。在这种情况下，您需要一个附加变量才能确定各个比例。以下各行中的代码使用非线性近似来绘制输出。

最初，数据：

library(ggplot2)
#Data
df <- structure(list(drug = c("drug_1","drug_1","drug_2","drug_2"),conc = c(100,33.33,11.11,3.7,1.23,0.41,0.14,0.05,100,0.05),mean_response = c(1156,1833,1744,1256,1244,1088,678,489,2322,1867,1333,944,567,356,200,177),std_dev = c(117,317,440,134,38,183,153,719,218,185,117,166,167,88,50)),class = "data.frame",row.names = c(NA,-16L))

在非线性最小二乘法中，您需要定义初始值以搜索理想参数。我们使用带有基本功能nls()的下一个代码来获取这些初始值：

#Drug 1
fm1 <- nls(log(mean_response) ~ log(a/(1+exp(-b*(conc-c)))),df[df$drug=='drug_1',],start = c(a = 1,b = 1,c = 1))
#Drug 2
fm2 <- nls(log(mean_response) ~ log(a/(1+exp(-b*(conc-c)))),df[df$drug=='drug_2',c = 1))

使用这种初始的参数方法，我们使用geom_smooth()绘制图。我们再次使用nls()查找正确的参数：

#Plot
ggplot(data=df,aes(y=mean_response,x= conc,color = drug)) +
  geom_pointrange(aes(ymax = (mean_response + std_dev),ymin = (mean_response - std_dev))) +
  geom_smooth(data = df[df$drug=='drug_1',method = "nls",se = FALSE,formula = y ~ a/(1+exp(-b*(x-c))),method.args = list(start = coef(fm1),algorithm='port'),color = "tomato")+
  geom_smooth(data = df[df$drug=='drug_2',method.args = list(start = coef(fm0),color = "cyan3")

输出：

curve-fitting ggplot2 model-fitting r r