使用ggplot将S型曲线拟合到数据点

Question

我有一个简单的数据框，记录了不同剂量下药物治疗的反应测量值：

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

我接下来想要做的事情是在这些数据上添加一个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 

我也尝试了这个链接提供的解决方案：将sigmoid曲线拟合到oxy-Hb数据上。

在这种情况下，我收到以下错误：

Computation failed in `stat_smooth()`:
Convergence failure: singular convergence (7) 

同时也没有任何线条被添加到图表中。

我已经试图查找这些错误的原因，但似乎无法找到与我的数据相符合的原因。

非常感谢任何提供帮助的人！

- FriendlyFred

---

我不建议一开始就在 ggplot 中尝试这样做：这会增加另一个步骤，使故障排除更加困难。 - Ben Bolker

Answer 1

正如我在评论中所说的，我只会在遇到非常简单的问题时使用 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参数...请注意加权和未加权估计之间的巨大差异...

- Ben Bolker

---

哦，这是一个非常好的观点（值得一份错误报告吧？似乎没有记录？），但是……根据经验，我从nls-without-weights和nls2-with-weights得到了不同的答案。我想。我可以再次检查代码和答案……（从未弄清楚为什么nls在包含权重时会遇到麻烦……我确实发现它不会因为1/sd权重（这些权重没有统计意义）而失败，但是会因为1/sd^2权重而失败…… - Ben Bolker

似乎在2013年12月16日提供了nls2 + weights的补丁（请参阅此处），但是CRAN版本的nls2可以追溯至2013-03-07！ - Ben Bolker

在 Github 上有一个新版本（0.3）的 nls2，修复了大多数已知问题，包括 subset= 和 weights= 参数的问题（并添加了一种新的抽样方法）。devtools::install_github("ggrothendieck/nls2") - G. Grothendieck

Answer 2

我建议采用接近您想要的方法。我也尝试使用二项式系列来设置您的数据，但是在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")

输出：

- Duck