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

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

# Access the data
females.N[1 ]
#> [1] 44
females.N[ 2]
#> [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[-1]
#There is no data for 1958 so we caon't get the data.

Population growth rate: vectorized

TASK

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

We are now trying to find the lambda of the two ratios. It is to find the differences between the first derivitives.

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

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(). The reason why you are using length is because if we change a sequence, this vector will still work.

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 calculates the differences between each two of the data in the sequence in one list. It is different because it does not set the variable before we do calculation, it just use the origional name for it in the calculation.


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

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

TASK

What does this do? Briefly describe what the [-1] is doing. This has dropped the first number in the list.

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

length(females.Ntemp[2:10]/females.Ntemp[1:9])
#> [1] 9

“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] 44 47 46 44 46 45 46 40 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)]

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?
What does females.N[-length(females.N)?

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.) This code has changed lambda.i into a vector instead of sequence. This code added an NA on the last of the sequence to help the sequence to have the same length with others.

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

lambda_log <- log(lambda.i)

Assemble the dataframe

bear_N <- data.frame(census,
                year.t,
                females.N,
                lambda.i, 
                lambda_log)
length(census) 
#> [1] 39
length(year.t)
#> [1] 39
length(females.N)
#> [1] 39
length(lambda.i) 
#> [1] 39
length(lambda_log)
#> [1] 39

TASK

List 3 functions that allow you to examine this dataframe.

1.summary() 2.dim() 3.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, 
     xlab = "Year", 
     ylab = "Population index (females + cubs)",  
     type = "b")

Bears love to eat trash. Yellowstone closed the last garbage dump in 1970 https://www.yellowstonepark.com/things-to-do/yellowstone-bears-no-longer-get-garbage-treats

Challenge questions

We will cover this all in the next lecture. Feel free to explore this code yourself.

CHALLENGE TASK

Plot a vertical line at 1970. Write a sentence or indicating if you think the population was impacted by this.

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

# I think the population was impacted by this.

Plot lambda

CHALLENGE TASK

Make a histogram of population growth rates
Write a brief sentence describing the shape of these data.

hist(bear_N$lambda.i)

Plot ln(lambda)

CHALLENGE TASK

Make a histogram of natural log population growth rates
Write a brief sentence describing the shape of these data.

hist(bear_N$lambda_log)

Mean of log(lambda)

CHALLENGE TASK

Briefly describe what happens when you delete na.rm = T It removes the na values from the list.

mean(bear_N$lambda_lo, na.rm = T)
#> [1] 0.02134027

Histogram with mean

In statistics the mean is often represented as the Greek letter “mu”. This can be represented as “u”.

CHALLENGE TASK Save the mean to an object called u

u <- mean(bear_N$lambda_log, na.rm = T)

CHALLENGE TASK Make a histogram with the mean plotted on it

hist(bear_N$lambda_log)
abline(v = u)

Density dependence

CHALLENGE TASK Make a graph that indicates if “density dependence” is occurring.

plot(lambda.i ~ females.N, data = bear_N, 
     type = "p",
     ylab = "Population growth Rate",
     xlab = "Females with cubs")
abline(v = 1970)

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

This will be the core topic of next lecture.

We can find the density and population growth based on the graphs, thus we can determine that the wether the population is going to extinct.

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?

We have found the population growth and watch the hat which shows the growth rate of this population, thus we can figure out wether the pupulation is likely to go extinct.

This is the step that we are eliminating the na values in the lists.

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 sampling of lambdas

This is the step we are choosing and saving a random sample from the hat.

# Choosing random samples from a hat.
sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#> [1] 1.051282
#Pull lambda from the heat and saving it.
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)

Grab the data

In this step we want to grab data from the dataset.

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

In this step we predict the number and save it

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

Simulation the hard way

Just like excel but worse. Simulate the population change year by year.

#1997 - 1998
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.1998 <- lambda_rand.t*N.1997
#1998 - 1999
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.1999 <- lambda_rand.t*N.1998
#1999 - 2000
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2000 <- lambda_rand.t*N.1999
#2000 - 2001
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2001 <- lambda_rand.t*N.2000
#2001 - 2002
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2002 <- lambda_rand.t*N.2001
#2002 - 2003
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2003 <- lambda_rand.t*N.2002
#2003 - 2004
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2004 <- lambda_rand.t*N.2003
#2004 - 2005
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2005 <- lambda_rand.t*N.2004

Plot population change

Plot the population change from 1997 to 2004.

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

Do simulation year by year

We are trying to simulate the population growth of the bear by doing calculation year by year. It is a worse for() loop

# Initial condition

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 in a hard way
N.current <- N.initial

# This is where the for() loop would be
t <- 1
  
  # Grab random lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # Determine population size
  N.t <- N.current*lambda_rand.t
  
  # Determine the year we are at
  year.t <- 1997+t
  
  # Plot the population
  # points update on existing graph
  points(N.t ~ year.t)

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

Now we are running the real for loop


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

# Start at 1997 pop size
N.current <- N.1997

# Start of for loop
for(t in 1:50){
  
  # grab random lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # Determine the population size
  N.t <- N.current*lambda_rand.t
  
  # Determine the year
  year.t <- 1997+t
  
  # Add new points to the existing graph
  points(N.t ~ year.t)
  
  # Update N.current
  N.current <- N.t
}

Goopy R plotting code

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

Still the for() loop simulation

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

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

Challenge questions

Plot lambda

Plot ln(lambda)

Mean of log(lambda)

Histogram with mean

Density dependence

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

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

Random sampling of lambdas

Grab the data

One round of population simulate

Simulation the hard way

Plot population change

Do simulation year by year