logistic增长曲线与数据的拟合

Question

我一直在尝试将logistic增长方程与我拥有的数据集进行拟合，结果喜忧参半。我通常使用这样的设置：

# Post PT
time <- 1:48

Diversity <- new8

plot(time, Diversity,log="y",las=1, pch=16, type="l")

logisticModel <- nls(Diversity~K/(1+exp(Po+r*time)), start=list(Po=25, r=-1.6, K=200),control=list(maxiter=1000,minFactor=.00000000001))

这里的目标是在逻辑上模拟随时间变化的Diversity；这是一条渐近的物种多样性曲线。然而，对于特定的数据集，我不能让模型工作，也不能找出原因。例如，在一次迭代中，拉取的分集(因此，new8)值为

[1]  25  22  68  72 126 141  82  61  97 126 101 110 173 164 160 137 122 113 104 104 109 102 107 122 149 127 137 146 185 188 114  91 102 132 147
[36] 148 151 154 165 215 216 206 205 207 207 220 200 204

# plot via this, and it is a nice species diversity curve beginning to level off

plot(Diversity,type="l")

这些数据开始达到它的极限，但我不能用逻辑曲线来拟合它。如果我尝试，我会得到一个超过最大迭代次数的错误，无论我把迭代次数调得多高。我一次又一次地尝试启动参数，但都没有成功。例如，当前的代码如下所示：

# Post PT
time <- 1:48

Diversity <- new8

plot(time, Diversity,log="y",las=1, pch=16, type="l")

logisticModel <- nls(Diversity~K/(1+exp(Po+r*time)), start=list(Po=25, r=-1.6, K=200),control=list(maxiter=1000,minFactor=.00000000001))

任何帮助都是非常感谢的。整天坐在我的沙发上。如果有人有更好的方法从数据中得出逻辑增长曲线，我很想听听！顺便说一句，我对这些数据集使用过SSlogis，但也没有成功。

Answer 1

对于包含指数项的模型，数值不稳定性通常是一个问题。尝试以您的起始参数评估您的模型：

> 200/(1+exp(25-1.6*df$norm_time))
 [1] 2.871735e-09 2.969073e-09 3.069710e-09 3.173759e-09 3.281333e-09 3.392555e-09 3.507546e-09 3.626434e-09 3.749353e-09
[10] 3.876437e-09 4.007830e-09 4.143676e-09 4.284126e-09 4.429337e-09 4.579470e-09 4.734691e-09 4.895174e-09 5.061097e-09
[19] 5.232643e-09 5.410004e-09 5.593377e-09 5.782965e-09 5.978979e-09 6.181637e-09 6.391165e-09 6.607794e-09 6.831766e-09
[28] 7.063329e-09 7.302742e-09 7.550269e-09 7.806186e-09 8.070778e-09 8.344338e-09 8.627170e-09 8.919589e-09 9.221919e-09
[37] 9.534497e-09 9.857670e-09 1.019180e-08 1.053725e-08 1.089441e-08 1.126368e-08 1.164546e-08 1.204019e-08 1.244829e-08
[46] 1.287023e-08 1.330646e-08 1.375749e-08

由于预测数据具有如此小的值，因此nls()估计梯度所需的参数中的任何适度更改都可能会在数据中产生非常小的更改，仅略高于甚至低于minFactor()。

最好对数据进行标准化，使其数值范围在一个友好的范围内，如0到1。

require(stringr)
require(ggplot2)
new8 <- '25  22  68  72 126 141  82  61  97 126 101 110 173 164 160 137 122 113 104 104 109 102 107 122 149 127 137 146 185 188 114  91 102 132 147 148 151 154 165 215 216 206 205 207 207 220 200 204'
Diversity = as.numeric(str_split(new8, '[ ]+')[[1]])
time <- 1:48
df = data.frame(time=time, div=Diversity)

# normalize time
df$norm_time <- df$time / max(df$time)

# normalize diversity
df$norm_div <- (df$div - min(df$div)) / max(df$div)

使用这种标准化多样性的方法，您的Po参数始终可以被假定为0。这意味着我们可以从模型中消除它。模型现在只有两个自由度，而不是三个，这也使得拟合变得更容易。

这将我们引向以下模型：

logisticModel <- nls(norm_div~K/(1+exp(r*norm_time)), data=df, 
                     start=list(K=1, r=-1.6), 
                     control=list(maxiter=1000, minFactor=.00000000001))

你的数据看起来并不适合我的模型，但我不是你这个领域的专家：

ggplot(data=df, aes(x=norm_time, y=norm_div)) + 
geom_point(log='y') + 
  geom_line(aes(x=norm_time, y=predict(logisticModel)), color='red') +
  theme_bw()

quartz.save('~/Desktop/SO_31236153.png', type='png')

summary(logisticModel)

Formula: norm_div ~ K/(1 + exp(r * norm_time))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
K   0.6940     0.1454   4.772 1.88e-05 ***
r  -2.6742     2.4222  -1.104    0.275    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1693 on 46 degrees of freedom

Number of iterations to convergence: 20 
Achieved convergence tolerance: 5.895e-06

​

​