Search code examples
rmodelnls

How to work out the time of when the asymptote occurs in a logistic model in R?

I have fitted logistic growth models in R using nls function.

pop.ss <- nls(MRDRSLT ~ SSlogis(TIME, phi1, phi2, phi3), data = testing, na.action = na.omit)

Formula: MRDRSLT ~ SSlogis(TIME, phi1, phi2, phi3)

Parameters:
      Estimate Std. Error t value Pr(>|t|)   
phi1   0.23179    0.03317   6.988  0.00602 **
phi2 431.16641   36.68846  11.752  0.00132 **
phi3  79.58386   29.09809   2.735  0.07164 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01824 on 3 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 7.55e-06

PH1 is the growth rate where the model starts to plateau. Is there away of find the time this occurred easily ?

Or do I need to store the coefficients and rearrange to make t the subject ?

Solution

  • Let's try to clarify our thinking here by creating our own logistic curve using our own parameters.

    We will look at a time interval of 1:100, and make time = 50 the midpoint of our curve. On the y axis we will have our curve go from 0 on the left and plateau at 20 on the right. Notice that because of the shape of the logistic curve, the plateau is only ever "reached" at time = infinity.

    We will also add some random noise so we don't get singularities in our nls fit.

    phi1  <- 20
phi2  <- 50
phi3  <- 0.1
TIME  <- 1:100

testing <- data.frame(
  TIME,
  MRDRSLT = phi1/(1 + exp(phi3 * (phi2 - TIME))) + rnorm(length(TIME))
)

plot(testing$TIME, testing$MRDRSLT)

    enter image description here

    This clearly has a logistic shape. Now let's see if we can reverse engineer those parameters:

    pop.ss <- nls(MRDRSLT ~ SSlogis(TIME, phi1, phi2, phi3), 
              data = testing, na.action = na.omit)

summary(pop.ss)
#> 
#> Formula: MRDRSLT ~ SSlogis(TIME, phi1, phi2, phi3)
#> 
#> Parameters:
#>      Estimate Std. Error t value Pr(>|t|)    
#> phi1  20.3840     0.2751   74.09   <2e-16 ***
#> phi2  50.0940     0.5863   85.45   <2e-16 ***
#> phi3  10.4841     0.4722   22.20   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.057 on 97 degrees of freedom
#> 
#> Number of iterations to convergence: 0 
#> Achieved convergence tolerance: 4.528e-06

    Now let us plot what those parameters actually represent on the data:

    plot(testing$TIME, testing$MRDRSLT)
abline(h = summary(pop.ss)$coef["phi1", 1], col = "red")
abline(v = summary(pop.ss)$coef["phi2", 1], col = "red")

    We can see that phi1 is the value of y at the upper plateau, and that phi2 is the midpoint where the inflection point occurs.

    Now, you are looking for the point in time at which the curve "starts to plateau". This is a slippery concept, but technically speaking the point at which the curve starts to plateau is the point at which the second derivative of the curve becomes negative. That occurs at the inflection point, which is equal to phi2. The curve may not "look like" it is plateauing at that point, but any other point that you choose on the curve to call the "start" of the plateau is entirely arbitrary.

    Created on 2020-09-21 by the reprex package (v0.3.0)