Although not a main focus of this course, population models are a strong motivation for what we do. Briefly population models:
Give us a good way to illustrate how parameters such as abundance, survival, and recruitment can potentially interact to determine population growth
Give us a framework for making predictions about population growth based on empirical estimates of these parameters
Potentially allow us to interpret in a conservation decision making framework the impacts of management actions on influencing key parameters, and thus influencing growth.
In addition, several of the statistical models we will deal with explicitly or implicitly incorporate assumptions from simple population models into the structure of statistical model. These include (but are not limited to)
Jolly-Seber and related models for open marked populations
Robust Design models
State -space models
There are many ways to construct population models in R (or other platforms including spreasheets) and there is no one right way. Fundamentally the approaches boil down to
Models based on an explicit equation solving where the population will be at some time t, given initial conditions and assumed parameter values
Models based on discrete iteration over a series of time steps from some starting point and assumed parameter values
We will start out with simple growth models that involve neither stochastic effects nor density feedback, and then create more complex and realistic models with random effects,density, or both. We will also briefly consider age structure and models based on direct modeling of stochastic processes.
Exponential growth model
The exponential growth model assumes that the population starts out at time t=0 at some initial N(0), and growth at constant exponential rate r. If r>0 this results in theoretically unlimited growth over time, if r<0 eventual asymptotic decline to 0, and if r=0 the population stays where it is, N0. The fundamental equation of growth is
N(t+1)= N(t) exp(r)
This can be solve to find a value for N at any future t, given N0 and r as
N(t)= N(0) exp(rt)
There are thus 2 basic ways to simulate growth: directly (by the 2nd equation) or iteratively (by the first). The 2 give equivalent results (this can also be proven mathematically). Both approaches are illustrated in the attached code, which you can run and compare. You can also change the values of N0, r, and T (number of years simulated) in any sensible (no negative N or time please!) fashion you wish.
Adding stochastic effects
As a (small) concession to realism, suppose that we assume that population growth rates are not constant every year, but instead allow rates to vary stochastically, as they would under environmental influences. A simple way to do this is to specify a mean (average) growth rate e.g. r.mean=0.1, a standard deviation e.g. r.sd= 0.2, and then draw annual values of r from a Normal distribution with this mean and sd. The attached code does this, using the iterative form of the growth model. This is a sensible way to proceed since the stochastic influence needs to propagate through the values of N over time, and I don't really see how to do that in a 1-step solution.
I have added a couple of other modifications to this model that might make it more useful in some situations, particularly involving small populations. First, I've rounded N to integer values, and second I've assigned an extinction threshold of N=2. In combination these features avoid nonsensical situations such as N is a small fractional value, or N=1, both equivalent to extinction in sexually reproducing organisms. Otherwise it could be possible for a population to 'recover' numerically from what should be extinction.
Adding density feedback
For additional realism, I have code 2 versions of the logistic model of density dependent growth. This model involves the parameters
N0 initial population size
rmax maximum growth rate
K an upper limit to growth (so-called carrying capacity)
and T number of years simulated
The first version is the exact, continuous model based on integration of differential equations, and the second is the discrete-time approximation for the same model (see Williams et al. 2002 Chapter 8 for more). The second will give results that are numerically very close but not the same as the first, as you can see by running each under the same conditions. The second approach tends to be used more in simulation, because it is easier to add features such as stochastic growth (as in the second example).
Modeling stochastic processes directly
The above models are fine for many applications, particularly where were dealing with large populations. They have some undesirable features however, such as allowance of inadmissible values such as noninteger values of N. In addition, they somewhat "black box" the actual processes that are occurring in nature, which are fundamentally stochastic birth and death processes. An alternative approach which is very flexible is to model these processes directly, for instance F(t) as a Poisson process based on current N(t) and some average recruitment rate, and D(t) as a binomial process with some average rate. Looked at this way we start with some current N(t) and get
F(t) ~ Poisson( f*N(t))
D(t) ~Binomial (N(t), d)
where f and d are respectively per-capita average recruitment and mortality probability. Then we get next year's N as
N(t+1) = N(t) +F(t) -D(t)
This approach works well at simulating small populations and totally avoids the non-integer and zeros problems. See the attached code for a simple example.
Incorporating age effects
It is well known that populations with age structure may grow at variable rates, even without stochastic effects or density feedbacks, due to instability in the age structure. Entire volumes have been written on models involving age (and size ) structure, eg. by Hal Caswell.; see the Wikipedia summary and references therein. We briefly cover age structured population growth in Chapter 3 of Conroy and Carroll (2009); a more in-depth coverage is provided by Williams et al. (2002).
Briefly, the idea is that the population is now stratified into a vector of abundance, with the number in each element representing number alive in an age class at time t. If for example we had 3 possible age classes (birth year, one year , and two year) we might have
N (t)' = [ 100, 50, 25]
If the population is experiencing constant (but age-specific) survival and recruitment, we describe growth by multiplying the current population vector by a projection matrix A containing age- specific recruitment and survival rates. Examples of this are provided in the accompanying R code. Then growth rate is very succinctly described by
N (t+1)= A N (t)
We then try to find a point a which both the proportions of N in each age class as well as overall growth become constant, known as the stable age distribution. As with simpler population growth models there is a iterative solution to finding this point: we can start with some initial N (0) and iteratively find the point where stability occurs. Alternatively, we can use eigenanalysis (see also) to find a direct, analytical solution. Both approaches are illustrated in the attached code.
Code from class discussion showing how to put population simulation inside a replicate loop and how to do this all as a user-defined function