使用 R 中的逻辑函数拟合 sigmoid 曲线

Question

我有遵循 sigmoid 曲线的数据，我想拟合一个逻辑函数来提取每个参与者的三个（或两个）参数。我在网上找到了一些方法，但我不确定哪个是正确的选择。

This tutorial 解释说你应该像这样使用 nls() 函数：

fitmodel <- nls(y~a/(1 + exp(-b * (x-c))), start=list(a=1,b=.5,c=25))
 ## get the coefficients using the coef function 
params=coef(fitmodel)

...您显然需要起始值来找到最佳拟合值 (?)。

然后this post 解释说，要获取起始值，您可以使用“自启动模型可以为您估算好的起始值，因此您不必指定它们”：

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

但是，我还在其他地方读到您应该使用 SSlogis 函数来拟合逻辑函数。请有人确认这两个步骤是否是最好的方法吗？还是应该使用从以前的类似数据中提取的值作为起始值？

另外，如果我根本不希望逻辑函数由渐近线定义，该怎么办？

谢谢！

Answer 1

没有最好的方法，但 SSlogis 确实消除了设置起始值的麻烦，而如果您指定公式，您可以更好地控制参数化。

如果问题真的是如何将a 固定在预定水平，这里的值1，无需重写公式，然后在运行nls 之前设置a 并从起始值中省略它。

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

或者，这会将 a=1 代入公式 fo，得到 fo2，而无需自己重写公式。

fo2 <- do.call("substitute", list(fo, list(a = 1)))
nls(fo2, start = list(b = 0.5, c = 25))

Answer 2

作为@G。格洛腾迪克写道，没有通用的“最佳方式”，它总是取决于你的特定目标。使用SSLogis 是个好主意，因为您不需要指定起始值，但定义自己的函数更灵活。请参阅下面的示例，我们使用启发式方法自己派生起始值，而不是手动指定它们。然后，我们拟合了一个逻辑模型，并作为一个小小的奖励，使用具有显式滞后阶段的 Baranyi 增长模型。

# time (t)
x <- c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20)

# Algae cell counts (Mio cells per ml)
y <- c(0.88, 1.02, 1.43, 2.79, 4.61, 7.12,
       6.47, 8.16, 7.28, 5.67, 6.91)

## we now plot the data linearly and logarithmically
## the layout function is another way to subdivide the plotting area
nf <- layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE), respect = TRUE)
layout.show(nf) # this shows how the plotting area is subdivided

plot(x, y)
plot(x, log(y))

## we see that the first points show the steepest increase,
## so we can estimate a start value of the growth rate
r <- (log(y[5]) - log(y[1])) / (x[5] - x[1])
abline(a=log(y[1]), b=r)

## this way, we have a heuristics for all start parameters:
## r:  steepest increase of y in log scale
## K:  maximum value
## N0: first value

## we can check this by plotting the function with the start values
f <- function(x, r, K, N0) {K /(1 + (K/N0 - 1) * exp(-r *x))}
plot(x, y, pch=16, xlab="time (days)", ylab="algae (Mio cells)")
lines(x, f(x, r=r, K=max(y), N0=y[1]), col="blue")

pstart <- c(r=r, K=max(y), N0=y[1])
aFit   <- nls(y ~ f(x, r, K,N0), start = pstart, trace=TRUE)

x1 <- seq(0, 25, length = 100)
lines(x1, predict(aFit, data.frame(x = x1)), col = "red")
legend("topleft",
       legend = c("data", "start parameters", "fitted parameters"),
       col = c("black", "blue", "red"),
       lty = c(0, 1, 1),
       pch = c(16, NA, NA))

summary(aFit)
(Rsquared <- 1 - var(residuals(aFit))/var(y))

## =============================================================================
## Approach with Baranyi-Roberts model
## =============================================================================

## sometimes, a logistic is not good enough. In this case, use another growth
## model
baranyi <- function(x, r, K, N0, h0) {
  A <- x + 1/r * log(exp(-r * x) + exp(-h0) - exp(-r * x - h0))
  y <- exp(log(N0) + r * A - log(1 + (exp(r * A) - 1)/exp(log(K) - log(N0))))
  y
}

pstart <- c(r=0.5, K=7, N0=1, h0=2)
fit2   <- nls(y ~ baranyi(x, r, K, N0, h0), start = pstart, trace=TRUE)

lines(x1, predict(fit2, data.frame(x = x1)), col = "forestgreen", lwd=2)

legend("topleft",
       legend = c("data", "logistic model", "Baranyi-Roberts model"),
       col = c("black", "red", "forestgreen"),
       lty = c(0, 1, 1),
       pch = c(16, NA, NA))