Computational Ecology: Simulating population dynamics

Introduction

The explores code needed to run population simulations. All populations vary in size from year to year. Even if they are consistently growing, consistently shrinking, or staying stable, there will be variation year to year. If we want to project population dynamics into the future we need to take into account this variation. We do this by randomly selecting values for key parameters. But the values aren’t totally random - they need to be grounded in reality.

For a basic population simulation we need 2 things

An initial population size
A distribution of population growth rates (lambda)

We can get get our set of population growth rates from the observed past dynamics of the population. We can then make the assumption that what happened in the past will resemble what will happen in the future - a strong assumption, but we’ll start here - and we’ll simulate future possibilities by randomly selecting population growth rates from those we’ve observed in the past.

Preliminaries

This code builds our dataframe.

# the vectors of data
census <- 1:39
year.i   <- 1959:1997
females.N <- c(44,47,46,44,46,
               45,46,40,39,39,
               42,39,41,40,33,
               36,34,39,35,34,
               38,36,37,41,39,
               51,47,57,48,60,
               65,74,69,65,57,
               70,81,99,99)

# lambda
lambda.i <- females.N[-1]/females.N[-length(females.N)]

# make lambda the same length as the rest of the data
lambda.i <- c(lambda.i,NA)

# log of lambda
lambda_log <- log(lambda.i)

# make the dataframe
bear_N <- data.frame(census,
                year.i,
                females.N,
                lambda.i, 
                lambda_log)

Always plot your data to make sure it looks right!

plot(females.N ~ year.i, data = bear_N, 
     type = "b",
     ylab = "Population index (females + cubs)",
     xlab = "Year")
abline(v = 1970)

How do we determine if a population is likely to go extinct?

TASK: ADD 1-2 sentences of notes here

First, we create a vector of the lambda values. This will allow us to “pull numbers out of a hat” to simulate population dynamics.

hat_of_lambdas <- bear_N$lambda.i

We can use is.na() to check if there are any NA values. We aren’t going to want to use this in our simulations. Note that the last value in the vector is TRUE

is.na(hat_of_lambdas)
#>  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [37] FALSE FALSE  TRUE

This is a more complicated bit of code for those who are curious.

any(is.na(hat_of_lambdas) == TRUE)
#> [1] TRUE

We can drop the NA.

Recall that the length of the vector is 39

length(hat_of_lambdas)
#> [1] 39

Confirm that the last value is NA

hat_of_lambdas[39]
#> [1] NA

Make a vector without the NA using negative indexing

hat_of_lambdas[-39]
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000

Drop the NA using generalized code instead of hard coding

hat_of_lambdas[-length(hat_of_lambdas)]
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000

Advanced: we can use the function na.omit() also.

na.omit(hat_of_lambdas)
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000
#> attr(,"na.action")
#> [1] 39
#> attr(,"class")
#> [1] "omit"

Create a vector of just lambdas, no NAs

hat_of_lambdas <- hat_of_lambdas[-length(hat_of_lambdas)]

Let’s look at a histogram of the lambdas

hist(hat_of_lambdas)

Randomly sampling lambdas

TASK: ADD 1-2 sentences of notes here

To simulate into the future we need to randomly draw population growth rates from those we have observed.

# TASK: add 1 comment to the code below

# Examples:

# sample() draws random subsets
## from a vector
sample(x = hat_of_lambdas,  #input: lambdas
       size = 1,            # how many lambdas: just 1
       replace = TRUE)      # sample with replacement: yes
#> [1] 0.9782609

#store the lambdsa
lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)

Initial population size

TASK: ADD 1-2 sentences of notes here

We’ll start out population simulations with the final population size of bears in our data set.

head(bear_N)
#>   census year.i females.N  lambda.i  lambda_log
#> 1      1   1959        44 1.0681818  0.06595797
#> 2      2   1960        47 0.9787234 -0.02150621
#> 3      3   1961        46 0.9565217 -0.04445176
#> 4      4   1962        44 1.0454545  0.04445176
#> 5      5   1963        46 0.9782609 -0.02197891
#> 6      6   1964        45 1.0222222  0.02197891
tail(bear_N)
#>    census year.i females.N  lambda.i lambda_log
#> 34     34   1992        65 0.8769231 -0.1313360
#> 35     35   1993        57 1.2280702  0.2054440
#> 36     36   1994        70 1.1571429  0.1459539
#> 37     37   1995        81 1.2222222  0.2006707
#> 38     38   1996        99 1.0000000  0.0000000
#> 39     39   1997        99        NA         NA
N.1997 <- 99

Sample population growth from 1997 to 1998

TASK: ADD 1-2 sentences of notes here

Simulate the population 1 year into the future.

#draw a random lambda
sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#> [1] 1.25

#save a random lambda
lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)

#multiple the random lambda by
## the initial population size
## to determine simulated population size
## in 1998
lambda_rand.i*N.1997
#> [1] 109.7027

N.1998 <- lambda_rand.i*N.1997

Simulate the population many years into the future

TASK: ADD 1-2 sentences of notes here

In the code chunk below I carry out a short population simulate “by hand” from 1997 to 2005. This is not so different than how we would do this in a spreadsheet like Excel.

The basic workflow is

Pick a random lambda
Multiple lambda by last year’s population size to get this years population size
Move ahead in time 1 year
Go back to 1 / repeat

# 1st iteration
## Pick a random lambda
lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)

## Multiple lambda by last year's population size to get this years population size
N.1998 <- lambda_rand.i*N.1997

## Move ahead in time 1 year
### (this is implied by changing the names of the vectors)

## repeat the previous steps


lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.1999 <- lambda_rand.i*N.1998

lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2000 <- lambda_rand.i*N.1999

lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2001 <- lambda_rand.i*N.2000

lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2002 <- lambda_rand.i*N.2001

lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2003 <- lambda_rand.i*N.2002

lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2004 <- lambda_rand.i*N.2003

lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2005 <- lambda_rand.i*N.2004

Plot “by hand” simulation

TASK: ADD 1-2 sentences of notes here

Let’s gather together our seven or so years of data an plot it.

year <- seq(1997, 2004)
N.rand <- c(N.1998,N.1999,N.2000,N.2001,N.2002,N.2003,N.2004,N.2005)
df.rand <- data.frame(N.rand, year)
plot(N.rand ~ year, data = df.rand, type = "b")

Simplify the simulation

TASK: ADD 1-2 sentences of notes here

In this code we’ll skip the step of writing the code out many time. Instead, we’ll use a variable called “i” to keep track of how many years into the future we have gone, and initially assign it the value 1. We can then run the code that does the 1st year of simulated population change, then change i to 2. This is an outline of what will become a for() loop.

NOTE: It can be very tempting to use the letter t instead of i when you’re simulating a time series. However, t() is a function in R, and if you assign t a value, the t() function won’t work. Its good practice to never overwrite and R function with your own information.

# Simulatieon by hand

## Preliminaries

N.1997 <- 99
N.initial <- 99

# Plot the initial population dynamcs population

## this plot won't work
plot(N.1997 ~ c(1997))

## let's set the y  axis limits
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50))

## and the x limits
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))

If you run this you’ll get a graph with 2 points

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))

# The actual simulation
#
N.current <- N.initial

# Sei initila year
i <- 1
  
  # sample random lambda
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # calculate population size next year
  N.i <- N.current*lambda_rand.i
  
  # determine the current year
  year.i <- 1997+i
  
  # plot the current population size
  points(N.i ~ year.i)

  
  # Update the population size
  N.current <- N.i

Now change i to 2. You’ll get a new graph with 3 points.

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
points(N.i ~ year.i)

# Sei initila year
i <- 2
  
  # sample random lambda
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # calculate population size next year
  N.i <- N.current*lambda_rand.i
  
  # determine the current year
  year.i <- 1997+i
  
  # plot the current population size
  points(N.i ~ year.i)

  
  # Update the population size
  N.current <- N.i

You could continue this an let the population grow.

Simulating with a for() loop

TASK: ADD 1-2 sentences of notes here

Now let’s let R handle the process of running the code, updating the year, and re-running the code. This is done with a for() loop.

All for loops look like this

for(i in 1:ending.number){
  ...
  # do something cool
  #save the output
  x <- y
  
}

When the for() loop is done the output is something you want, like a time series of simulated population dynamics.

i is a variable that R will change from 1, to 2, to 3, to … the ending.number. ending.number is how many iterations or cycles you want the for() loop to run for. That is, how many time do you want the code to be repeated.

1:ending.number makes a vector of integers from 1 to ending.number. If ending numbers was 10, then it would be

1 2 3 4 5 … 10

To run a for() loop for 50 cycles you’d do this

for(i in 1:50){
  ...
  # do something cool
  #save the output
  x <- y
  
}

Let’s run our simulation as a for() loop for 50 years.


# Ploi initila conditions in 1997.
## x and y axis limits are set to not cut the graph off
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))

# Set  starting population size.
N.current <- N.1997

# The for loop
#  
for(i in 1:50){
  
  # Draw a random lambda
  lambda_rand.i <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # Calcualte new population size
  N.i <- N.current*lambda_rand.i
  
  # increase the year by 1 year
  ## the for loop says "for(i in 1:50){...}
  ## So i starts out as 1
  year.i <- 1997+i
  
  # Plot the new population size
  points(N.i ~ year.i)
  
  # Update the current population size
  N.current <- N.i
}

TASK: ADD 1-2 sentences of notes here

Set the plotting window to have a 3 by 3 grid.

par(mfrow = c(3,3), mar = c(1,1,1,1))

TASK: ADD 1-2 sentences of notes here

Run the code several times to simulate different populations. The will produce 9 replicated runs of the simulate.

# Run it once
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
N.current <- N.1997
for(i in 1:50){
  
  lambda_rand.i <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  N.i <- N.current*lambda_rand.i
  
  year.i <- 1997+i
  
  points(N.i ~ year.i)
  
  N.current <- N.i
}


# Run it a 2nd time
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
N.current <- N.1997
for(i in 1:50){
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
   N.i <- N.current*lambda_rand.i
  year.i <- 1997+i
  points(N.i ~ year.i)
  N.current <- N.i
}


# Run it a 3rd time
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
N.current <- N.1997
for(i in 1:50){
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
   N.i <- N.current*lambda_rand.i; year.i <- 1997+i
  points(N.i ~ year.i); N.current <- N.i
}


# Run it a 4th time
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550)); N.current <- N.1997
for(i in 1:50){
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
   N.i <- N.current*lambda_rand.i; year.i <- 1997+i ;points(N.i ~ year.i)
   N.current <- N.i}


# Run it a 5th time
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550)); N.current <- N.1997
for(i in 1:50){lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
   N.i <- N.current*lambda_rand.i; year.i <- 1997+i ;points(N.i ~ year.i)
   N.current <- N.i}


# Run it a 6th time
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550)); N.current <- N.1997
for(i in 1:50){lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
   N.i <- N.current*lambda_rand.i; year.i <- 1997+i ;points(N.i ~ year.i)
   N.current <- N.i}


# Run it a 7th time
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550)); N.current <- N.1997
for(i in 1:50){lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
   N.i <- N.current*lambda_rand.i; year.i <- 1997+i ;points(N.i ~ year.i)
   N.current <- N.i}


# Run it a 8th time
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550)); N.current <- N.1997
for(i in 1:50){lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
   N.i <- N.current*lambda_rand.i; year.i <- 1997+i ;points(N.i ~ year.i)
   N.current <- N.i}


# Run it a 9th time
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550)); N.current <- N.1997
for(i in 1:50){lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
   N.i <- N.current*lambda_rand.i; year.i <- 1997+i ;points(N.i ~ year.i)
   N.current <- N.i}

Plot a 10 x 10 grid of population

TASK: ADD 1-2 sentences of notes here

We can expand this if we want. I’ll set up a 10 x 10 grid, but only plot one simulation. I could add more if I kept re-running the for() loop.

First, set up the plotting panel.

par(mfrow = c(10,10), mar = c(0,0,0,0), xaxt = "n", yaxt = "n")

TASK: ADD 1-2 sentences of notes here

Now run the simulation.


plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550), cex = 0.5)
abline(h = N.1997)
abline(h = 0, col = "red")
N.current <- N.1997
for(i in 1:50){
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  N.i <- N.current*lambda_rand.i
  year.i <- 1997+i
  points(N.i ~ year.i, cex = 0.5)
  
  N.current <- N.i
}

NOTE: Lecture ended here

The lecture ended here. Below are some challenge questions for you to figure out if you want what I’m doing.

CHALLENGE ADD TITLE HERE

CHALLENGE TASK: ADD 1-2 sentences of notes here

## ADD TITLE HERE
num.years <- 50

## ADD TITLE HERE
N.vector <- rep(NA, num.years)

## ADD TITLE HERE
N.vector[1] <- N.1997

t<- 1
  # print the year
  print(t)
#> [1] 1
  
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  N.i <- N.vector[i]*lambda_rand.i
  
  N.vector[i + 1] <- N.i

CHALLENGE: ADD TITLE HERE

TASK: ADD 1-2 sentences of notes here

## ADD TITLE HERE
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD TITLE HERE
#don't uncomment print() in notebook mode!
for(i in 1:num.years){
  # print the year
  #print(t)
  
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  N.i <- N.vector[i]*lambda_rand.i
  
  N.vector[i + 1] <- N.i
  
  
}

CHALLENGE ADD TITLE HERE

par(mfrow = c(1,1), mar = c(3,3,3,3), xaxt = "s", yaxt = "s")

df <- data.frame(N.vector, year = seq(1997, 1997+50))

plot(N.vector ~ year, data = df)


plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

CHALLENGE ADD TITLE HERE

TASK: ADD 1-2 sentences of notes here

## ADD TITLE
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD TITLE
lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)


## ADD TITLE
t <- 1


  print(t)
#> [1] 1
  
  lambda_rand.i <- lambdas_vector[i]
  
  N.i <- N.vector[i]*lambda_rand.i
  
  N.vector[i + 1] <- N.i
  
  t <- i + 1
  N.vector[1:t]
#>  [1] 99 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [51] NA

CHALLENGE ADD TITLE HERE

## ADD NOTES
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD NOTES
lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)


## ADD NOTES
for(i in 1:num.years){
  
  ## ADD NOTES
  lambda_rand.i <- lambdas_vector[i]
  
  ## ADD NOTES
  N.i <- N.vector[i]*lambda_rand.i
  
  ## ADD NOTES
  N.vector[i + 1] <- N.i
  
}

## ADD NOTES
df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

CHALLENGE: How can I shorten the code?

lambdas_vector[i]
N.i

num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)

for(i in 1:num.years){
  lambda_rand.i <- lambdas_vector[i]
  
  N.i <- N.vector[i]*lambda_rand.i
  
  N.vector[i + 1] <- N.i
  
}

df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

Shortened code

num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)

for(i in 1:num.years){
  N.vector[i + 1] <- N.vector[i]*lambdas_vector[i]
  
}

df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

Computational Ecology: Simulating population dynamics

Nathan Brouwer

9/15/2020

Introduction

Preliminaries

How do we determine if a population is likely to go extinct?

Randomly sampling lambdas

Initial population size

Sample population growth from 1997 to 1998

Simulate the population many years into the future

Plot “by hand” simulation

Simplify the simulation

Simulating with a for() loop

Plot a 10 x 10 grid of population

NOTE: Lecture ended here

CHALLENGE ADD TITLE HERE

CHALLENGE: ADD TITLE HERE

CHALLENGE ADD TITLE HERE

CHALLENGE ADD TITLE HERE

CHALLENGE ADD TITLE HERE

CHALLENGE: How can I shorten the code?

Shortened code