Pop_growth

1. Modeling the basic exponential/geometric population growth model

Let’s think of a hypothetical population that you have observed over two years. Say there were 100 individuals in the first year, and 150 individuals in the second year. Then, we can say that the growth rate of the population over those two years is 𝑁2𝑁1=1.5 . More generally, we can think of the annual growth rate, 𝜆 as the number of individuals in one year compared to the year before:

𝜆=𝑁𝑡+1𝑁𝑡

If this rate of growth is constant over 𝑡 number of years and the initial population size is 𝑁0 , then the population size at year 𝑡 is:

𝑁𝑡=𝑁0𝜆𝑡 This is the model of geometric population growth, i.e., a population growing in discrete time steps with constant growth rate.

Exponential population growth occurs when a population is growing in continuous time with a constant growth rate, 𝑟 :

𝑁𝑡=𝑁0𝑒𝑟𝑡

Here, 𝑟 is the instantaneous rate of growth which is the birth rate (births per individual) minus the death rate (deaths per individual) at any given time scale.

Geometric and exponential growth models are very similar, except that they differ in several mathematical properties. But given a certain time unit, we can go between the two models pretty easily. You just have to remember that 𝜆=𝑒𝑟 and ln𝜆=𝑟

Let’s now project population growth over 10 years with starting population of 100 and 𝜆=1.5 . What we need to do is assign the constants (𝑁0 and 𝜆 ) and then project the population sizes at 𝑁𝑡 when 𝑡 ranges from 0 (intial population) to 10

N0 = 100
lambda = 1.5
t = 0:10
N = N0 * lambda^t
round(N, 0)

##  [1]  100  150  225  338  506  759 1139 1709 2563 3844 5767

Let’s plot this on a linear scale and log scale.

plot(t, N, type = "o", pch = 19, las = 1)

plot(t, N, type = "o", log = "y", pch = 19, las = 1)

These two plots describe exponential (or geometric) population growth: i.e., population under a constant growth rate.

2. Varying the growth rate

Now, let’s try projecting the population with different growth rates, 𝜆 . Let’s think about how to do this. We will now have two parameters that vary (𝑡 and 𝜆 ). This means that we want an output that is a matrix, with one parameter (say 𝑡 ) on the rows and the other parameter (𝜆 ) on the columns. We will do this by employing the sapply() function. Now, we will define a sequence numbers for 𝜆 (here, it will be 0.6, 0.8, 1, 1.2 and 1.4) and then use sapply() to run multiple versions of the population growth rate equation. This function will return a vector or matrix, depending on what is the simplest form.

N0 = 100
lambda = seq(0.6, 1.4, 0.2)
t = 0:10
N = sapply(lambda, function(lambda) N0 * lambda^t)
N

##              [,1]      [,2] [,3]     [,4]      [,5]
##  [1,] 100.0000000 100.00000  100 100.0000  100.0000
##  [2,]  60.0000000  80.00000  100 120.0000  140.0000
##  [3,]  36.0000000  64.00000  100 144.0000  196.0000
##  [4,]  21.6000000  51.20000  100 172.8000  274.4000
##  [5,]  12.9600000  40.96000  100 207.3600  384.1600
##  [6,]   7.7760000  32.76800  100 248.8320  537.8240
##  [7,]   4.6656000  26.21440  100 298.5984  752.9536
##  [8,]   2.7993600  20.97152  100 358.3181 1054.1350
##  [9,]   1.6796160  16.77722  100 429.9817 1475.7891
## [10,]   1.0077696  13.42177  100 515.9780 2066.1047
## [11,]   0.6046618  10.73742  100 619.1736 2892.5465

colors = brewer.pal(5, "Set1")
matplot(t, N, type = "o", las = 1, pch = 1:5, col = colors, lty = 1:5)
legend("topleft", legend = c("0.6", "0.8", "1", "1.2", "1.4"),
    title = "lambda", pch = 1:5, lty = 1:5, col = colors, cex = 0.8)

Ok, we can see that there are multiple trajectories here, but it’s kind of hard to see what is goin on at the bottom of the figure. Let’s try plotting the population sizes on a log scale:

colors = brewer.pal(5, "Set1")
matplot(t, N, type = "o", las = 1, pch = 1:5, col = colors, lty = 1:5,
    log = "y", yaxt = "n")
marks = c(1, 5, 10, 50, 100, 500, 1000)
axis(side = 2, at = marks, labels = format(marks, drop0trailing = TRUE),
    las = 1)
legend("bottomleft", legend = c("0.6", "0.8", "1", "1.2", "1.4"),
    title = "lambda", pch = 1:5, lty = 1:5, col = colors, cex = 0.8)

You can see that two lines go up (𝜆>1 ), two lines go down (𝜆<1 ) and one line stays flat (𝜆=1 ). That is, populations increase when 𝜆>1 and vice versa.

Since the instantaneous rate of growth, 𝑟=ln𝜆 :

Population does not change when 𝜆=1 or 𝑟=0 Population increases when 𝜆>1 or 𝑟>0 Population decreases when 𝜆<1 or 𝑟<0

3. Simulating stochasticity

Thus far, we have assumed that population growth rate is constant through time. This is clearly not the case in most natural populations. Even disregarding density-dependent mechanisms, populationg growth rates fluctuate due to different sources of stochasticity (i.e., randomness). One major source of randomness in population growth rate is environmental stochasticity: i.e., random fluctuations in population growth rate due to good and bad years.

Let’s say that a population has a stable population growth rate–i.e., average instantaneous growth rate, 𝑟¯=0 . But the actual growth rate fluctuates around this average with some variance, 𝜎2𝑟=0.01 . Hence, for each time step (say, 1,000), we could randomly draw a number from a normal distribution with mean of 0 and variance of 0.01 (or 𝑠𝑑=𝑣𝑎𝑟‾‾‾‾√=0.1 ), like so:

set.seed(2)
rs = rnorm(1000, mean = 0, sd = 0.1)
hist(rs, xlab = "r")

For simplicity, we can also just take the distribution of lambdas from a log-normal distribution directly (output not shown):

set.seed(2)
lambdas = rlnorm(100, meanlog = 0, sdlog = 0.1)
round(lambdas, 2)

##   [1] 0.91 1.02 1.17 0.89 0.99 1.01 1.07 0.98 1.22 0.99 1.04 1.10 0.96 0.90 1.20
##  [16] 0.79 1.09 1.00 1.11 1.04 1.23 0.89 1.17 1.22 1.00 0.78 1.05 0.94 1.08 1.03
##  [31] 1.08 1.03 1.11 0.97 0.93 0.94 0.84 0.91 0.95 0.98 0.96 0.82 0.92 1.21 1.06
##  [46] 1.22 0.97 0.99 0.98 0.89 0.92 1.23 0.95 1.14 0.90 0.82 0.97 1.10 1.12 1.18
##  [61] 0.84 1.23 0.93 1.02 1.05 0.92 0.82 0.95 1.01 0.91 0.91 1.03 0.99 1.04 0.99
##  [76] 0.91 1.14 1.08 1.11 0.87 1.10 0.84 0.95 0.87 0.80 1.20 0.94 0.97 0.96 1.04
##  [91] 1.17 1.18 0.89 0.87 0.86 0.88 1.22 1.00 0.92 0.94

Ok, now let’s use this method of generating stochasticity in population growth rates to simulate population growth. Let’s remind ourselves what the population projection over 20 years looks like under geometric population growth with 𝜆=1.2 and 𝑁0=100 .

Here, let’s try generating this plot slightly differently–we will use a for-loop to calculate 𝑁𝑡 for each time step 𝑡 and store it as a vector.

N0 = 100  #initial population size
times = 20  #number of years into the future
N = vector(length = times)  #empty vector to store pop. sizes
N[1] = N0  #initial population size should be the first N
lambda = 1.2  #growth rate
# start loop: Take previous year's N and multiply by lambda
for (t in 2:times) {
    N[t] = N[t - 1] * lambda
}
plot(1:times, N, type = "o", las = 1)

This is a bit more round-about than before, but this lays the groundwork for our stochasticity simulation. Now, let’s simulate a population with the same mean growth rate1 but with standard deviation of 0.1.

set.seed(2)
N0 = 100  #initial population size
times = 20  #number of years into the future
N = vector(length = times)  #empty vector to store pop. sizes
N[1] = N0  #initial population size should be the first N

# lambdas--we only need 19 numbers because growth only
# happens between 2 years.
lambda = rlnorm(times - 1, meanlog = 0, sdlog = 0.1)

# start loop: Take previous year's N and multiply by lambda
for (t in 2:times) {
    N[t] = N[t - 1] * lambda[t - 1]
}
plot(1:times, N, type = "o", las = 1)

Now, let’s try running this simulation multiple times (10x) and plotting the results in one figure. To do this, we will employ the sapply() function and run the for-loop sequence inside it. The nice thing about this is that the sapply() function will automatically organize the output as a matrix for us, so that we can directly use the matplot() function to plot the results.

# multiple simulations
set.seed(2)
sims = 10
outmat = sapply(1:sims, function(x) {
    times = 20
    N0 = 100
    N = vector(length = times)
    N[1] = N0
    lambda = rlnorm(times - 1, meanlog = 0, sdlog = 0.1)
    for (t in 2:times) {
        N[t] = N[t - 1] * lambda[t - 1]
    }
    N
})
matplot(1:times, outmat, type = "l", las = 1, ylab = "Population Size",
    xlab = "Years")
abline(h = 100, lty = 2)

##4. Population Viability Analysis: An over-simplified example Population Viability Analysis (PVA) is a population modeling method to assess extinction risk of a population or species. It is the mathematical foundation of conservation biology, and it is used worldwide to guide management decisions for threatened species. Here, I will walk through a simplified example of PVA using an example dataset.

4.1 The Song Sparrow Dataset Here, we will be using the sparrows dataset from the R package primer. You should have already installed and loaded the package at the beginning of the module. If you haven’t yet, run library(primer) to load it. The primer package is a package that is associated with a textbook called “A Primer of Ecology with R”, written by Henry Stevens3. The sparrow dataset is based on counts of Song Sparrows (Melospiza melodia) in Darrtown, OH from 1966 to 2003.

data(sparrows)
head(sparrows)

##   Year Count ObserverNumber
## 1 1966    34              1
## 2 1967    40              1
## 3 1968    42              1
## 4 1969    54              1
## 5 1970    49              1
## 6 1971    71              1

plot(Count ~ Year, data = sparrows, type = "o", las = 1)

sparrows %>%
    ggplot(aes(x = Year, y = Count)) + geom_line(color = "grey") +
    geom_point(shape = 21, color = "black", fill = "#69b3a2",
        size = 6) + theme_bw() + ggtitle("Sparrow population")

##4.2 Calculating lambda from counts The first thing we can do with this data is to extract the annual population growth rates 𝜆 for each pair of years. That is, we need to calculate 𝑁2𝑁1 , 𝑁3𝑁2 , 𝑁4𝑁3 , etc.

counts = sparrows$Count
l = counts[-1]/counts[-length(counts)]
round(l, 2)

##  [1] 1.18 1.05 1.29 0.91 1.45 0.93 0.85 1.29 1.07 0.48 0.86 0.81 0.85 1.41 0.90
## [16] 1.46 0.88 1.03 1.19 0.77 1.38 0.96 0.91 0.78 1.25 3.00 0.80 0.97 0.72 1.25
## [31] 1.07 0.41 0.92 1.00 1.26

hist(l, breaks = 20, main = "Histogram of lambdas", col = "gray")

mean(log(l))

## [1] 0.006709703

sd(log(l))

## [1] 0.3366946

You’ll notice that the average ln𝜆=𝑟 is really close to 0, but slightly above it. This means that, the population is pretty stable at the moment. However, the variance is fairly big.

Let’s just assume that we can estimate the distribution of annual population growth rates 𝜆 as a log-normal distribution with those mean-logs and sd-logs. With this, we can generate 50 projected 𝜆 values:

set.seed(2)
sim.l = rlnorm(50, meanlog = mean(log(l)), sdlog = sd(log(l)))
round(sim.l, 2)

##  [1] 0.74 1.07 1.72 0.69 0.98 1.05 1.28 0.93 1.96 0.96 1.16 1.40 0.88 0.71 1.83
## [16] 0.46 1.35 1.02 1.42 1.16 2.04 0.67 1.72 1.94 1.01 0.44 1.18 0.82 1.31 1.11
## [31] 1.29 1.12 1.45 0.91 0.78 0.82 0.56 0.74 0.83 0.93 0.88 0.52 0.76 1.91 1.24
## [46] 1.97 0.91 0.98 0.95 0.67

##4.3 Projecting the population over 20 years Let’ use this to simulate this Song Sparrow population over the next 20 years (or rather, from 2003 to 2023).

First, let’s look at the population size in 2003, the last year in this dataset:

tail(sparrows)

##    Year Count ObserverNumber
## 31 1997    84              2
## 32 1998    90              2
## 33 1999    37              3
## 34 2000    34              3
## 35 2001    34              3
## 36 2003    43              3

So we will use N0 = 43 as our starting point, and simulate the population over the next 20 years.

# Single simulation
set.seed(2)
time = 21
N0 = 43
N = vector(length = time)
N[1] = N0
sim.l = rlnorm(time, meanlog = mean(log(l)), sdlog = sd(log(l)))
for (t in 2:time) {
    N[t] = N[t - 1] * sim.l[t - 1]
}
par(mar = c(4, 4, 1, 4))
plot(1:(time), N, type = "o", las = 1, xaxt = "n")
axis(side = 1, at = c(1, 6, 11, 16, 20), labels = c(2003, 2008,
    2013, 2018, 2023))

projected <- data.frame(Year = seq(2004, 2024, by = 1), Count = N,
    ObserverNumber = "Projected")

sparrows2 <- rbind(sparrows, projected)

sparrows2$ObserverNumber <- as.factor(sparrows2$ObserverNumber)

sparrows2 %>%
    ggplot(aes(x = Year, y = Count, fill = ObserverNumber)) +
    geom_line(color = "grey") + geom_point(shape = 21, color = "black",
    size = 6) + theme_bw() + ggtitle("Sparrow population")

Using the same principles as the stochastic simulations in Section 3, we can also conduct multiple simulations and plot the results

set.seed(2)
sims = 10
outmat = sapply(1:sims, function(x) {
    time = 21
    N0 = 43
    N = vector(length = time)
    N[1] = N0
    sim.l = rlnorm(time, meanlog = mean(log(l)), sdlog = sd(log(l)))
    for (t in 2:time) {
        N[t] = N[t - 1] * sim.l[t - 1]
    }
    N
})
par(mar = c(4, 4, 1, 4))
matplot(1:time, outmat, type = "l", las = 1, lty = 5, ylim = c(0,
    200), ylab = "N", xaxt = "n", xlab = "Year")
axis(side = 1, at = c(1, 6, 11, 16, 20), labels = c(2003, 2008,
    2013, 2018, 2023))

You can see that there are scenarios in which the population increases to very large numbers, and other scenarios in which the opulation declines.

##4.4 The PVA: What is the extinction risk? Now, let’s say we want to ask a question about extinction risk of this population. For example, we might ask what are the chances that this population goes extinct within 100 years, given the observed pattern in population sizes? To figure this out, we need to conduct a large number of simulations (in this case, let’s do 1,000x) projecting the population for 100 years (or 101, if starting from year 0). Then, we want to ask: in how many iterations of this model does the population size dip below 2 (the minimum number of individuals needed to reproduce)?

Let’s start by doing the simulations in the same way as above, but changing sims=1000 and time=101.

set.seed(2)
sims = 1000
outmat = sapply(1:sims, function(x) {
    time = 101
    N0 = 43
    N = vector(length = time)
    N[1] = N0
    sim.l = rlnorm(time, meanlog = mean(log(l)), sdlog = sd(log(l)))
    for (t in 2:time) {
        N[t] = N[t - 1] * sim.l[t - 1]
    }
    N
})

dim(outmat)

## [1]  101 1000

time = 101
par(mar = c(4, 4, 1, 4))
matplot(1:time, outmat, type = "l", las = 1, lty = 5, ylab = "N",
    xaxt = "n", xlab = "Year")
axis(side = 1, at = c(1, 6, 11, 16, 20), labels = c(2003, 2008,
    2013, 2018, 2023))

So, what we have now is a matrix with 101 rows (year 0 to year 100), and 1,000 columns, each representing a single iteration of the simulation. Thus, we could figure out how many columns has at least one time point in which N < 2 by using the apply() function to ask whether the minimum of each column is less than 2.

# which columns have at least one value less than 2?
minpop = apply(outmat, 2, function(x) min(x) < 2)
sum(minpop + 0)/sims  #proportion of columns with TRUE

## [1] 0.279

Thus, there is approximately a 28% chance that the population will go extinct within 100 years due purely to stochasticity.

Now, this is a really simplified version of PVA–it does not account for uncertainty in the estimates of mean and sd of ln𝜆, nor does it take into account any other realistic factors such as temporal autocorrelation of growth rates, density-dependent factors, habitat change/loss, etc. etc. But I hope that you got a sense for how to use R for simulating complex dynamics.

Stevens, M.H.H. (2009) A Primer of Ecology with R, Springer, 2nd Printing.↩︎

Ok now, so far we projected populations, but we assumed

Exponential growth model: ‘r’ is constant Stochastic growth model: ‘r’ is a random variable

Now we say that ‘r’ is regulated by the population density (e.g., because ressources are limited,..) (https://rpubs.com/vrchka26/465217)

# Paramters: Logistic growth models assume linear
# density-dependence of 'r' Every population has a 'K', and
# when N = K, 'r' ~ 0

# K = carrying capacity N= Population size
#'r' = per capita growth rate
# r max= maximum growth rate

# fiExample: Australian fox

### Create some fake population data, of let's say newly
### introduced fox population
t = (seq(1900, 1920, 1))  #####time vector
N = (sort(round((runif(21, 200, 2489)), 0), decreasing = F))  ###population vector, rounded and sorted in ascending order
pop <- as.data.frame(cbind(t, N))  ###bind columns
(pop)

##       t    N
## 1  1900  264
## 2  1901  266
## 3  1902  373
## 4  1903  411
## 5  1904  461
## 6  1905  631
## 7  1906  764
## 8  1907  767
## 9  1908  945
## 10 1909  957
## 11 1910 1379
## 12 1911 1429
## 13 1912 1468
## 14 1913 1621
## 15 1914 1657
## 16 1915 1820
## 17 1916 1971
## 18 1917 2056
## 19 1918 2085
## 20 1919 2253
## 21 1920 2261

plot(pop$t, pop$N)
model = nls(N ~ SSlogis(t, a, b, c), data = pop)
plot(pop$t, pop$N)
lines(pop$t, predict(model))

summary(model)

## 
## Formula: N ~ SSlogis(t, a, b, c)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## a 2540.6276    97.0999   26.16 8.92e-16 ***
## b 1910.4013     0.4713 4053.62  < 2e-16 ***
## c    4.6306     0.2823   16.40 2.86e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 63.57 on 18 degrees of freedom
## 
## Number of iterations to convergence: 0 
## Achieved convergence tolerance: 3.678e-07

See if the population is growing as ‘density-dependent’

############## use ratio method to find 'r'
len = length(pop$t)
## calculate r as r = log(Nt+1/Nt)
r1 = pop$N[-1]/pop$N[-len]
r = log(r1)
r

##  [1] 0.007547206 0.338082111 0.097014795 0.114804829 0.313907820 0.191261927
##  [7] 0.003919012 0.208698126 0.012618464 0.365310486 0.035616300 0.026926031
## [13] 0.099142313 0.021965496 0.093827763 0.079704527 0.042221319 0.014006508
## [19] 0.077493806 0.003544532

Fit linear model

y = r  ##where y=r, x= pop$N[1len-1]
x = pop$N[1:len - 1]
fit.lm = lm(y ~ x)
plot(x, y, xlab = "Population", ylab = "Population growth rate")
abline(fit.lm)

##### looking at the plot it looks like this population has
##### a logistic growth trend

####### find rmax and K
rmax = fit.lm$coefficients[1]
rmax

## (Intercept) 
##   0.2059058

slope = fit.lm$coefficients[2]
############## K #carrying capacity
K = -rmax/slope
K

## (Intercept) 
##    2463.761

Now we will see how well the model fits our data We will project our population using logistic growth model

N0 = pop$N[1]  #starting population
time = 1:20
N_t = rep(NA, length(t))  ##create a empty vector to store results
N_t[1] = N0


##### Write your own function to project population You can
##### give different values of N0, R, K and time

c.logist <- function(N0, r, K, time) {
    N_t = rep(NA, length(t))  #empty vector
    N_t[1] = N0  #initial pop size
    for (i in 2:length(time)) {
        N_t[i] = K/(1 + ((K - N0)/N0) * exp(-r * time[i]))
    }

    return(N_t)
}

N_t = c.logist(N0 = N0, r = 0.32, time = time + 1, K = K)


pop$Nt <- N_t

ggplot() + geom_line(data = pop, aes(x = t, y = N)) + geom_point(data = pop,
    aes(x = t, y = Nt, color = "Nt"))

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).