Computational Ecology: Population Viability Analysis 9/15/2020

Data: Figure 5 (plus 10 years)

Make data vectors, calculate lambda, and put together dataframe with all necessary data.

census

The census period; an index from 1 to 39 of how many years of data have been collected.

census <- 1:39

year t

The year: 1959 to 1997 (Dennis et al use 1959-1987)

year.t   <- 1959:1997

Population size

Population size is recorded as the number of females with …

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)

Population growth rate: example

Population growth rate is…

Enter the population size for each year

females.N.1959 <- 44
females.N.1960 <- 47

Calculate the ratio of the 2 population sizes

lambda.59_60 <- females.N.1960/females.N.1959
lambda.59_60
#> [1] 1.068182

Access the population sizes by using bracket notation rather than hard coding

# Access the data
females.N[2]
#> [1] 47
females.N[1]
#> [1] 44

# store in objects
females.N.1959 <- females.N[2]
females.N.1960 <- females.N[1]

# confirm the output
females.N.1960/females.N.1959
#> [1] 0.9361702

Calculate lambda using bracket notation

lambda.59_60 <- females.N[2]/females.N[1]

The first year of data is 1959. What is lambda for 1958 to 1959?

females.N[1]
#> [1] 44
#lambda.58_59 <- females.N[1]/females.N[ ]
#no data for 1958

Population growth rate: vectorized

TASK

Briefly describe (1-2 sentence) what this code is doing.

#give the data for the second and third year
females.N[2:3]
#> [1] 47 46
#give the data for the first and second year
females.N[1:2]
#> [1] 44 47
females.N[2:3]/females.N[1:2]
#> [1] 1.0681818 0.9787234
#calculate two lambda in a single line
#calculate population growth ratio between two data points

This is similar t the previous code chunk, just using all of the data (no need to describe)

length(females.N)
#> [1] 39
females.N[2:39]/females.N[1:38]
#>  [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

TASK What does this do? Briefly describe in 1 to 2 sentences why I am using length().

#make things more flexible instead of using constant numbers
len <- length(females.N)
females.N[2:len]/females.N[1:len-1]
#>  [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

TASK What does this do? Briefly describe in 1 to 2 sentences what is different about this code chunk from the previous one.

#this code does the same thing above but calculates length inside one function
females.N[2:length(females.N)]/females.N[1:length(females.N)-1]
#>  [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

Negative indexing

Make a short vector to play with; first 10 years

females.N[1:10]
#>  [1] 44 47 46 44 46 45 46 40 39 39
females.Ntemp <- females.N[1:10]

Check - are there 10 numbers

females.Ntemp
#>  [1] 44 47 46 44 46 45 46 40 39 39

TASK

What does this do? Briefly describe what the [-1] is doing.

#remove the first elements
females.Ntemp[-1]
#> [1] 47 46 44 46 45 46 40 39 39

TASK How many lambdas can I calculate using the first 10 years of data? 9 lambdas

females.Ntemp[2:10]/females.Ntemp[1:9]
#> [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#> [8] 0.9750000 1.0000000

“Negative indexing” allows you to drop a specific element from a vector.

TASK Drop the the first element

females.Ntemp[-1]
#> [1] 47 46 44 46 45 46 40 39 39

TASK Drop the second element

females.Ntemp[-2]
#> [1] 44 46 44 46 45 46 40 39 39

TASK

How do you drop the 10th element? Type in the code below.

females.Ntemp[-10]
#> [1] 44 47 46 44 46 45 46 40 39

TASK How do you access the last element? Do this in a general way without hard-coding.

females.Ntemp[10]
#> [1] 39
females.Ntemp[length(females.Ntemp)]
#> [1] 39

TASK How do DROP the last element? Do this in a general way without hard-coding. By general, I mean in a way that if the length of the vector females.Ntemp changed the code would still drop the correct element.

females.Ntemp[-10]
#> [1] 44 47 46 44 46 45 46 40 39
females.Ntemp[-length(females.Ntemp)]
#> [1] 44 47 46 44 46 45 46 40 39

TASK Calculate the first 9 lambdas.

lambda.i <- females.Ntemp[-1]/females.Ntemp[-10]

Converting between these 2 code chunks would be a good test question : )

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

Calcualte lambdas for all data

TASK

Below each bulleted line describe what the parts of the code do. Run the code to test it.

What does females.N[-1] do? remove the first elements
What does females.N[-length(females.N)? remove the last elements

TASK Calculate lambdas for all of the data

females.N[-1]
#>  [1] 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
#> [26] 47 57 48 60 65 74 69 65 57 70 81 99 99
females.N[-length(females.N)]
#>  [1] 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
#> [26] 51 47 57 48 60 65 74 69 65 57 70 81 99

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

Finish putting together dataframe

Create special columns

TASK

What does this code do? Why do I include NA in the code? (I didn’t cover this in lecture, so just type 1 line - your best guess. “I don’t know” is fine.)

It means the value is not available.

lambda.i <- c(lambda.i,NA)

TASK

Check the help file; what type of log does log() calculate (I forgot to put this question on the test!)

It computes logarithms with base base. log(x, base)

lambda_log <- log(lambda.i)

Assemble the dataframe

bear_N <- data.frame(census,
                year.t,
                females.N,
                lambda.i, 
                lambda_log)

TASK

List 3 functions that allow you to examine this dataframe.

1.dim() 2.summary() 3.head()

Examing the population growth rates

Plotting the raw data

TASK

Plot a time series graph of the number of bears (y) versus time (x)
Label the y axis “Population index (females + cubs)”
Label the x axis “Year”
Change the plot to type = “b” so that both points and dots are shown.

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

census <- 1:39
year.t   <- 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.i <- females.N[-1]/females.N[-length(females.N)]
lambda.i <- c(lambda.i,NA)
lambda_log <- log(lambda.i)
bear_N <- data.frame(census,
                year.t,
                females.N,
                lambda.i, 
                lambda_log)

plot(females.N ~ year.t, 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?

To model population dynamics we randomly pull population growth rates out of a hat

hat_of_lambdas <- bear_N$lambda.i

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

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

Drop the NA

length(hat_of_lambdas)
#> [1] 39
hat_of_lambdas[39]
#> [1] NA
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
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

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"

hat_of_lambdas <- hat_of_lambdas[-length(hat_of_lambdas)]

hist(hat_of_lambdas)

Random Samplings of Lambda

Pulled a lambda from the hat, and the pulled lambda from hat and save it to object

# pulled a lambda from the hat
sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#> [1] 0.8421053
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE) #can't be false

Initialize population simulation

head(bear_N)
#>   census year.t 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.t 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

One round of population simulate

Pulled out a random number, and simulate how many bears in 1998

1.22807*99
#> [1] 121.5789
lambda_rand.t*N.1997
#> [1] 108
N.1998 <- lambda_rand.t*N.1997

Simulation the hard way

Pulled out a random lambda to get population of bears. Repeat to get population for the next year

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

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

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

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

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

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

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

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

Plot population change

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")

Population simulation

hard way

# initial conditions

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

# explore xlim = arguments
plot(N.1997 ~ c(1997))

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

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

# the not a for() loop for loop
#for loop the hard way
N.current <- N.initial

# this is where the for() loop would be
t <- 1
  
  # grab a lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # determine pop size
  N.t <- N.current*lambda_rand.t
  
  # get years
  year.t <- 1997+t
  
  # plot pop 
  #points() updates an existing graph
  points(N.t ~ year.t)

  
  # update N.current
  N.current <- N.t

for() loop


# make a new plot with starting pop size
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))

#start at 1997s pop size
N.current <- N.1997

# use for loop to generate pop simulation
for(t in 1:50){
  
  # grab a lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # determine pop size
  N.t <- N.current*lambda_rand.t
  
  # get years
  year.t <- 1997+t
  
  # points() updates an existing graph
  points(N.t ~ year.t)
  
  # update N.current
  N.current <- N.t
}

Goofy R plotting code/magic

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

keep on track of simulation to see the variability

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

Plot pop change

**Goofy R plotting code/magic with 10*10**

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

generate 100 simulation by using loop


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(t in 1:50){
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  N.t <- N.current*lambda_rand.t
  year.t <- 1997+t
  points(N.t ~ year.t, cex = 0.5)
  
  N.current <- N.t
}

print year

## initial years
num.years <- 50

## repeat numbers
N.vector <- rep(NA, num.years)

## assign it to object
N.vector[1] <- N.1997

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

Print year

use for loop to print year

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

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

Goofy R plotting code/magic with 1*1 and plot pop

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)))

Print year

hard way print year

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

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


## print year
t <- 1


  print(t)
#> [1] 1
  
  lambda_rand.t <- lambdas_vector[t]
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t
  
  t <- t+1
  N.vector[1:t]
#> [1]  99.0 101.2

Plot pop change

## 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)


## use for loop
for(t in 1:num.years){
  
  ## grab lambda
  lambda_rand.t <- lambdas_vector[t]
  
  ## get value
  N.t <- N.vector[t]*lambda_rand.t
  
  ## update value
  N.vector[t+1] <- N.t
  
}

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

How can I shorten the code?

lambdas_vector[t]
N.t

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(t in 1:num.years){
  lambda_rand.t <- lambdas_vector[t]
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t
  
}

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(t in 1:num.years){
  N.vector[t+1] <- N.vector[t]*lambdas_vector[t]
  
}

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: Population Viability Analysis 9/15/2020

Nathan Brouwer

9/15/2020

Data: Figure 5 (plus 10 years)

census

year t

Population size

Population growth rate: example

Population growth rate: vectorized

Negative indexing

Calcualte lambdas for all data

Finish putting together dataframe

Create special columns

Assemble the dataframe

Examing the population growth rates

Plotting the raw data

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

Random Samplings of Lambda

Initialize population simulation

One round of population simulate

Simulation the hard way

Plot population change

Population simulation

Plot pop change

print year

Print year

Goofy R plotting code/magic with 1*1 and plot pop

Print year

Plot pop change

Shortened code