Chapter 1 Fitting a Logistic Model to Data

Example 1.1 We are given population data concerning the growth of a species of yeast in a closed vessel. Use the data given below and the solution \(u(t)\) to the logistic differential equation to find estimates of the parameters \(u_0\), \(r\) and \(K\).

\[u(t) = \frac{K}{1+e^{-rt} \left (\frac{K}{u_0}-1 \right )}\] where \(K\) is the carrying capacity of the population, \(r\) is the growth rate and \(u_0\) is the initial population size at time \(t=0\). For more details, please review the lectures where we covered all details related to the solution of the logistic equation.

We use a graphical approach for this problem to develop some intuition about fitting the logistic model parameters to given data.

  • What is \(u_0\)?
  • What is a good estimate for \(K\)? What do you predict as the maximum sustainable population (carrying capacity), based on this model?
  • What is a good estimate for \(r\)? Adjust the value of \(r\) manually until the logistic solution fits the data well.
The yeast population p (in millions) as a function of time t (in hours), is given by the data below: 
t<-0:17 # time in hours
# population in millions
p<-c(9.6, 18.3, 29.0, 47.2, 71.1, 119.1, 174.6, 257.3, 350.7, 
 441.0, 513.3, 559.7, 594.8, 629.4, 640.8, 651.1, 655.9, 659.6)

Hints:

  • Start by creating a scatterplot of the data of yeast population against time.
  • Find the horizontal asymptote to \(u(t)\) by finding \(\lim_{t\to\infty}u(t)\).
  • Does the solution to the logistic equation level out? What is that level?
  • Use the last observation above to make a good guess for the value of \(K\).

Problem Solution:

# load the required packages
library(mosaic)
library(mosaicCalc)

First, we plot the population size p against time t, using the data given above. We use the plotting function gf_point(), which is available when you load the mosaic package to visualize the data. Note that the plotting function takes a formula p ~ t, where p and t are the given data vectors, and it creates the scatterplot of p against t:

gf_point(p ~ t)

Next, we define the logistic solution as an R function u (the dependent variable) with argument t (the independent variable) using makeFun() from the mosaic package. Note that the R function u also depends on the parameters u0, r and K and it is created using the formula u(t) ~ t, which makes u a function of t:

u<-makeFun(K/(1+exp(-r*t)*(K/u0 - 1)) ~ t)

The initial yeast population at time zero is p[1]=9.6, based on the data, so the parameter for the initial population is u0=9.6.

To find the horizontal asymptote, we need to find \(\lim_{t\to\infty}u(t)\). As \(t\to\infty\) the term \(e^{-rt}\) goes to 0:

\[\lim_{t\to\infty}u(t) = \lim_{t\to\infty}\frac{K}{1+0\cdot (K/u_0 -1)} = K\] Thus, the horizontal asymptote is \(u=K\) (on the \(tu\)-plane). The yeast data appear to converge to the horizontal level of around 670, so we can take this as a reasonable estimate for \(K=670\) to be the carrying capacity. We can plot the horizontal asymptote \(u=670\) along with the data to make sure that it fits the data well:

gf_point(p ~ t) %>% 
  gf_hline(yintercept=670, col="red")

In the code above, after we create the scatterplot of the data using gf_point(p ~ t), we then use the pipe operator %>%, which implements a composition of functions, and it composes the first plotting function with the second plotting function gf_hline, which plots the horizontal line \(u=670\), for which we only need to specify the yintercept=670 (col="red" is optional).

Finally, to find a good estimate for \(r\), we can manually adjust the value of \(r\) until we get a good fit to the data. Note that we first create the scatter plot of the data as before, using gf_point(p ~ t), and then we use the pipe operator %>% to compose the two plotting functions. The second plotting function, composed with the first one, is the slice_plot() function from the mosaicCalc package. Note that slice_plot() takes the formula u(t,r=r, u0=u0, K=K) ~ t, where we give the parameters r, u0 and K specific numerical values and leave u to be a function of time t only. We then specify the domain over which to plot the graph of the function u(t) by giving a range for t, t=range(0,20) (col="blue" is optional).

u0=9.6 # initial population
K=670  # carrying capacity
r=0.54 # growth rate
gf_point(p ~ t) %>% 
  slice_plot(u(t,r=r, u0=u0, K=K) ~ t, domain(t=range(0,20)), 
             col="blue")

We then start to manually adjust the parameter r until the model fits the data well.

Manually adjusting the parameter \(r\) seems to produce a reasonable fit when \(r\approx 0.54\).