Fitting logistic growth curves to data

Question

I've been attempting to fit logistic growth equations to data sets I have, with mixed results. I typically use a setup like this:

# 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))

The goal here is to model Diversity over time logistically; this is a species diversity curve that asymptotes. However, for particular datasets, I cannot get the model to work and can't for the life of me figure out why. As an example, in one iteration, the Diversity (and therefore, new8) value that is being pulled is

[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")

This data is beginning to reach its limit, yet I cannot fit a logistic curve to it. If I try, I get an exceeded max iterations error, no matter how high I crank up the iterations. I've played with the starting parameters over and over with no luck. Currently, for example the code looks like this:

# 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))

Any help is more than appreciated. Spent all day sitting on my couch stuck on this. If someone has a better way to coerce a logistic growth curve out of data, I'd love to hear it! As a side note, I've used SSlogis for these datasets with no luck, either.

Answer 1

Numerical instability is often a problem with models involving exponential terms. Try evaluating your model at your starting parameters:

> 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

With predicted data having such small values, it's likely that any moderate change in the parameters, as required by nls() to estimate gradients, will produce changes in the data that are very small, barely above or even below minFactor().

It's better to normalize your data so that its numerical range is within a nice friendly range, like 0 to 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)

With this way of normalizing diversity, your Po parameter can always be assumed to be 0. That means we can eliminate it from the model. The model now only has two degrees of freedom instead of three, which also makes fitting easier.

That leads us to the following model:

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))

Your data doesn't look like that great a fit to the model to me, but I'm not an expert in your field:

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