使用 R 拟合 S 形曲线答案

【问题标题】：Using R to fit a Sigmoidal Curve使用 R 拟合 S 形曲线
【发布时间】：2016-01-07 02:26:45
【问题描述】：

我读过一篇文章（R 中的 Sigmoidal Curve Fit）。它被标记为重复，但我看不到与帖子相关的任何内容。为帖子给出的答案还不够。

我读到了webpage

与其他人类似，他使用这种格式来适应线条：

fitmodel <- nls(y~a/(1 + exp(-b * (x-c))), start=list(a=1,b=.5,c=25))

问题是，在大多数情况下都给出了 a,b,c，我不知道应该将哪一组 a,b,c 用于我的数据集。有人可以给我一些关于如何获取参数的建议吗？

这是我的一组数字：

x <- c(3.9637878,3.486667,3.0095444,2.5324231,2.0553019,1.5781806,1.1010594,0.6242821)
y <- c(6491.314,6190.092,2664.021,2686.414,724.707,791.243,1809.586,541.243)

【问题讨论】：

你必须猜到a, b, c。如果您知道曲线应该如何看起来总是有助于用随机系数绘制曲线（例如a = 20, b = 0.1, c = 0.2, curve(a/(1 + exp(-b * (x-c))), 0, 100) 并查看您的猜测如何检查x-c 是否正确。不应该是@ 987654327@
这是一种可能的方式，但还有其他可能的方式让我以更“统计上令人信服”的方式做到这一点吗？例如，是否可以创建一个循环来找到 a、b、c 的最佳组合集。或者，如果有的话，我可以使用一些函数或命令让程序为我计算吗？
曾经我也在寻找这个“gral”。我还没有找到。我不是统计学家，但我认为猜测起始值是在nlm 中估计方程参数的常用方法

标签： r curve-fitting

【解决方案1】：

幸运的是，R 为逻辑模型提供了一个自启动模型。它使用了轻微的重新参数化，但实际上与您的模型相同：Asym/(1+exp((xmid-input)/scal))

自启动模型可以为您估算好的起始值，因此您不必指定它们。

plot(y ~ x)
fit <- nls(y ~ SSlogis(x, Asym, xmid, scal), data = data.frame(x, y))

summary(fit)
#Formula: y ~ SSlogis(x, Asym, xmid, scal)
#
#Parameters:
#      Estimate Std. Error t value Pr(>|t|)
#Asym 1.473e+04  2.309e+04   0.638    0.551
#xmid 4.094e+00  2.739e+00   1.495    0.195
#scal 9.487e-01  5.851e-01   1.622    0.166
#
#Residual standard error: 941.9 on 5 degrees of freedom
#
#Number of iterations to convergence: 0 
#Achieved convergence tolerance: 4.928e-06

lines(seq(0.5, 4, length.out = 100), 
      predict(fit, newdata = data.frame(x = seq(0.5, 4, length.out = 100))))

当然，您的数据并不真正支持该模型。估计的中点正好在您的数据范围的正确极限，因此参数估计（特别是渐近线）非常不确定。

【讨论】：

谢谢，您的回答对我很有用。还有一个问题：为什么曲线的形状不是“S”形？对我来说，它更像是一条指数曲线。我看到你的方程和我的一样，但问题出在哪里？
模型是S形曲线。您只是看不到 S 的上部，因为它超出了绘图限制。

【解决方案2】：

我用来拟合你的数据的代码：

 df <- data.frame(x=c(3.9637878,3.486667,3.0095444,2.5324231,2.0553019,1.5781806,1.1010594,0.6242821),                     
                  y=c(6491.314,6190.092,2664.021,2686.414,724.707,791.243,1809.586,541.243))

library(drc)
fm <- drm(y ~ x, data = df, fct = G.3())

plot(fm)
summary(fm)

拟合后的样子：

【讨论】：

【解决方案3】：

我看到了两个问题。

nls 的默认算法对起始参数非常敏感。在您的示例数据中，我发现使用algorithm='port' 很有用。或者切换到“健壮”的实现也可能会有所帮助。
它有助于理解参数在模型中的作用。

您的模型的简单解释是： sigmoid 在 y 中从 0 到 a。它在 x=c 处到达“中途”点。 b 具有斜率的作用，如果为负，则模型将从 a 变为 0。

具体到您发布的测试数据，我将估计起始值如下：

我注意到的第一件事 - 您的数据并非完全“接近”为零，因此添加大约 1000 的偏移量 d 可能会很有用。
a 等于或大于 5000
c 大于 2 - 可能是 3
b 需要猜测 - 从 x 2 到 3.5，您的信号从 1000 跳到 6000 会产生 5000 的差异 - 除以 a - 1/1.5 = 0.66 或更大的斜率...让四舍五入。

所以最终使用公式

fitmodel <- nls(y ~a/(1 + exp(-b * (x-c)) ) + d, start=list(a=5000,b=1,c=3, d=1000))

适合（没有 d 也可以）。尝试了一下，我发现设置 algorithm='port' 使命令对起始值更加不敏感。

【讨论】：