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 reproductive mature females with cubs in the park

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 (or can be expressed as) a ration of population sizes.

Enter the population size for each year

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

#(1960-1959)/1959
#3bears/44 we started with
(47-44)/44
#> [1] 0.06818182

Calculate the ratio of the 2 population sizes

lambda.59_60 <- females.N.1960/females.N.1959

44*1.068182
#> [1] 47.00001

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

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

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

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

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[ ]
# can't be done; there is no data for 1958

39 years of data but only 38 lambdas being calculated (lambda represents transition)

Population growth rate: vectorized

TASK

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

females.N[2:3]
#> [1] 47 46
females.N[1:2]
#> [1] 44 47

females.N[2:3]/females.N[1:2]
#> [1] 1.0681818 0.9787234

The first lines of code will print the second and third year of data, and the first and second year of data, respectively. The second line of code gives value for lambda, and calculates population growth rate.

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

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

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

By using length, you are allowing R to do the work for you, rendering the code without any human miscalculations. Generalization of the code above.

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


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

Generalizing even further, just calculating the length inside the function.

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[seq(1,10)]

Check - are there 10 numbers

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

TASK

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

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

Drops the first element from the vector.

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

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

You can calculate 9 lambdas

“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[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[-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)]

Calculate 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? Removes the first element from the vector.

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

What does females.N[-length(females.N)? Removes the number of elements in the vector.

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

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

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

NA means not applicable or missing data.

TASK

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

lambda_log <- log(lambda.i)

log10

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.

dim()
summary()
str()

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")
abline(v = 1970)

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

Essentially characterizing what could happen in the future. To model population dynamics we randomly pull population growth rates (lambda) out of a hat.

hist(bear_N$lambda.i)

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() is a handy function

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 Sampling of Lambda

You can pull a lambda from the hat using the sample function. It is relatively inefficient because it’s one at a time.

# pulled a lambda from the hat
sample(x = hat_of_lambdas, 
       size = 1,               #size = ...
       replace = TRUE)         #replace = TRUE aka recycle = TRUE
#> [1] 0.825

# pull lambda from hat and save to object
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)

Getting initial population size

Different functions can give you different ways to view and analyze the data frame.

head(bear_N)      # head() shows you top 6 rows in data frame
#>   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)      # tail() shows you bottom 6 rows of data
#>    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
dim(bear_N)       # dim() shows you size of the data frame
#> [1] 39  5
summary(bear_N)   # summary() shows you statistical summary like mean, median, etc.
#>      census         year.t       females.N        lambda.i     
#>  Min.   : 1.0   Min.   :1959   Min.   :33.00   Min.   :0.8250  
#>  1st Qu.:10.5   1st Qu.:1968   1st Qu.:39.00   1st Qu.:0.9452  
#>  Median :20.0   Median :1978   Median :44.00   Median :1.0000  
#>  Mean   :20.0   Mean   :1978   Mean   :49.79   Mean   :1.0281  
#>  3rd Qu.:29.5   3rd Qu.:1988   3rd Qu.:57.00   3rd Qu.:1.1038  
#>  Max.   :39.0   Max.   :1997   Max.   :99.00   Max.   :1.3077  
#>                                                NA's   :1       
#>    lambda_log      
#>  Min.   :-0.19237  
#>  1st Qu.:-0.05639  
#>  Median : 0.00000  
#>  Mean   : 0.02134  
#>  3rd Qu.: 0.09874  
#>  Max.   : 0.26826  
#>  NA's   :1

N.1997 <- 99

One Round of Population simulation

Pulling a random value out of a hat, and simulating how many bears there are in the year

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

Simulation the hard way

How you would do it in Excel… but worse.

Modeling a trajectory with random variation showing how a population declines while ignoring carrying capacity. Class is performing a crowd sourced simulation of bears.

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

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

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

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

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

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

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

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

Plot population change

Plotting the same kind of population change simulation performed in the past code chunk on a graph.

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

For() loop

Step by step code outlining how to make a for() loop. First outlines the long and hard way of making one.

# Initial conditions

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

# Explore xlim = argument
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

# Where the for() loop would be
t <- 1
  
  # grab a lambda...
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # determine population size
  N.t <- N.current*lambda_rand.t
  
  # establish the year
  year.t <- 1997+t
  
  # plot the population
  # points() adds a point to an existing plot / updates an existing graph
  points(N.t ~ year.t)

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

  # run and repeat for each year you want to project population

Easier way of writing the loop, without having to repeat for every year


# make a new plot with starting population size
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550), xlab = "year.t", ylab = "N.current")

N.current <- N.1997

# start at 1997s population size
for(t in 1:50){
  
  # grab a lambda...
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # determine population size
  N.t <- N.current*lambda_rand.t
  
  # establish the year
  year.t <- 1997+t
  
  # plot the population
  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))

Can make multi-panel plot of random simulations by running code chunk over and over as many times as you need

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
}

Computational Ecology: Population Viability Analysis 9/15/2020 - 9/17/2020

Nathan Brouwer … edited by Dara Czernikowski

Data: Figure 5 (plus 10 years)

census

year t

Population size

Population growth rate: example

Population growth rate: vectorized

Negative indexing

Calculate 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 Sampling of Lambda

Getting initial population size

One Round of Population simulation

Simulation the hard way

Plot population change

For() loop