用 ggplot 将 sigmoidal 曲线拟合到点

Posted 2023-03-12

【中文标题】用 ggplot 将 sigmoidal 曲线拟合到点

【英文标题】：Fitting a sigmoidal curve to points with ggplot

【发布时间】：2020-12-13 13:18:53

【问题描述】：

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

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

conc <- c(100.00, 33.33, 11.11, 3.70, 1.23, 0.41, 0.14, 
        0.05, 100.00, 33.33, 11.11, 3.70, 1.23, 0.41, 0.14, 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, 200, 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

我想对这些数据做的下一件事是在绘图中添加一条 sigmoidal 曲线，它适合每种药物的绘制点。之后，我想计算这条曲线的 EC50。 我意识到我的数据中可能没有 S 形曲线的整个范围，但我希望能用我所拥有的得到最好的估计。此外，drug_1 的最终点并未遵循 S 形曲线的预期趋势，但这实际上并不意外，因为药物所在的溶液可以抑制高浓度的反应（每种药物在不同的溶液中）。我想从数据中排除这一点。

我在将 sigmoidal 曲线拟合到我的数据的步骤中遇到了困难。我查看了其他一些将 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) 

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

我已尝试查找这两个错误，但似乎找不到对我的数据有意义的原因。

任何帮助将不胜感激！

【问题讨论】：

我不建议尝试在 inside ggplot 开始这样做：它在流程中添加了另一层，使得故障排除更加困难。

【参考方案1】：

我会建议下一个接近你想要的方法。我还尝试使用binomial family 设置您的数据，但在 0 和 1 之间的值存在一些问题。在这种情况下，您需要一个额外的变量来确定各自的比例。以下行中的代码使用非线性近似来绘制您的输出。

最初，数据：

library(ggplot2)
#Data
df <- structure(list(drug = c("drug_1", "drug_1", "drug_1", "drug_1", 
"drug_1", "drug_1", "drug_1", "drug_1", "drug_2", "drug_2", "drug_2", 
"drug_2", "drug_2", "drug_2", "drug_2", "drug_2"), conc = c(100, 
33.33, 11.11, 3.7, 1.23, 0.41, 0.14, 0.05, 100, 33.33, 11.11, 
3.7, 1.23, 0.41, 0.14, 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, 200, 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',], start = c(a = 1, b = 1, 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 = "nls", se = FALSE,
              formula = y ~ a/(1+exp(-b*(x-c))),
              method.args = list(start = coef(fm0),
                                 algorithm='port'),
              color = "cyan3")

输出：

【参考方案2】：

正如我在评论中所说，我只会使用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,Asym,xmid,scal),data=df,
            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,pframe0),
               drug="drug_2",wts=FALSE),
    data.frame(pframe0,mean_response=predict(p1B,pframe0),
               drug="drug_1",wts=TRUE),
    data.frame(pframe0,mean_response=predict(p2B,pframe0),
               drug="drug_2",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 参数...注意加权和未加权估计之间的巨大差异...

【讨论】：

嗯。这是一个很好的观点（我猜值得一个错误 [infelicity] 报告？似乎没有记录？）但是......根据经验，我从 nls-without-weights 和 nls2-with-weights 得到不同的答案.我想。我想我可以仔细检查代码和答案......（从来没有弄清楚为什么 nls 在包含权重时会遇到麻烦......我确实发现它不会因 1/sd 权重而失败（其中没有统计意义）但使用 1/sd^2 权重...

看起来 nls2 + weights 的补丁是 2013 年 12 月 16 日的 supplied here（但 nls2 的 CRAN 版本可以追溯到 2013-03-07！）

在 github 上有一个新版本 (0.3) 的 nls2 修复了大多数已知问题，包括带有子集= 和 weights= 参数的问题（并添加了一种新的采样方法）。 devtools::install_github("ggrothendieck/nls2")