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 females with lil’ cubs

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 expressed as a ratio of population sizes

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

We get 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[1]
females.N.1960 <- females.N[2]

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

1.068182, confirmed.

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

Trick question, we don’t have the data for the 58th year. There will only be 28 lambdas calculated

Population growth rate: vectorized

TASK

Briefly describe (1-2 sentence) what this code is doing. Using vectorization, two of the lambdas will be calculated at once in line 117.

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

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

This is similar to 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(). We are working on trying to compute all the 38 lambdas at once in a general fashion (so we must subtract one from the second term’s ratio)

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
len 
#> [1] 39

TASK What does this do? Briefly describe in 1 to 2 sentences what is different about this code chunk from the previous one. Here, we are generalizing even further with the length command.

2:39 
#>  [1]  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
#> [26] 27 28 29 30 31 32 33 34 35 36 37 38 39
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

2:39 = seq(2,39)

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

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. The negative removes a value and the 1 describes that it’s removing the first element

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 = (n-1) = (10-1)

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[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. Cutting out the certain chunk of that vector that will allow us to align and calculate lambda.i

lambda.i <- females.Ntemp[-1 ]/females.Ntemp[-10]
lambda.i
#> [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#> [8] 0.9750000 1.0000000

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? Drops the first element in the females.N vector
What does females.N[-length(females.N)? Drops the last element in the females.N vector

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) ]
lambda.i
#>  [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

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.) Does it vectorize all the lambdas? IDK.

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!) By default, the natural log.

lambda_log <- log(lambda.i)
?log()

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.summary 2.view 3.length

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)
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
#> 7       7   1965        46 0.8695652 -0.13976194
#> 8       8   1966        40 0.9750000 -0.02531781
#> 9       9   1967        39 1.0000000  0.00000000
#> 10     10   1968        39 1.0769231  0.07410797
#> 11     11   1969        42 0.9285714 -0.07410797
#> 12     12   1970        39 1.0512821  0.05001042
#> 13     13   1971        41 0.9756098 -0.02469261
#> 14     14   1972        40 0.8250000 -0.19237189
#> 15     15   1973        33 1.0909091  0.08701138
#> 16     16   1974        36 0.9444444 -0.05715841
#> 17     17   1975        34 1.1470588  0.13720112
#> 18     18   1976        39 0.8974359 -0.10821358
#> 19     19   1977        35 0.9714286 -0.02898754
#> 20     20   1978        34 1.1176471  0.11122564
#> 21     21   1979        38 0.9473684 -0.05406722
#> 22     22   1980        36 1.0277778  0.02739897
#> 23     23   1981        37 1.1081081  0.10265415
#> 24     24   1982        41 0.9512195 -0.05001042
#> 25     25   1983        39 1.3076923  0.26826399
#> 26     26   1984        51 0.9215686 -0.08167803
#> 27     27   1985        47 1.2127660  0.19290367
#> 28     28   1986        57 0.8421053 -0.17185026
#> 29     29   1987        48 1.2500000  0.22314355
#> 30     30   1988        60 1.0833333  0.08004271
#> 31     31   1989        65 1.1384615  0.12967782
#> 32     32   1990        74 0.9324324 -0.06995859
#> 33     33   1991        69 0.9420290 -0.05971923
#> 34     34   1992        65 0.8769231 -0.13133600
#> 35     35   1993        57 1.2280702  0.20544397
#> 36     36   1994        70 1.1571429  0.14595391
#> 37     37   1995        81 1.2222222  0.20067070
#> 38     38   1996        99 1.0000000  0.00000000
#> 39     39   1997        99        NA          NA

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

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

We randomly pull growth rates to make predictions on the future. Averaging the growth rates is one way to gain insight on the likelihood of extinction

hist(bear_N$lambda.i)

hat_of_lambdas <- bear_N$lambda.i
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

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 lambda sampling

The function sample will allow us to randomly generate a value from our hat of lambdas. The argument claiming truth just means that the value in the hat can be recycled if size exceeds the size of our hat vector.

# pulled a random lambda from the hat
sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#> [1] 1.076923
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)

Getting inital population size

We want to look at some of the properties our bear vector has

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
summary(bear_N)
#>      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

One round of population simulated

Multiplying random growth rate by number of bears from previous year to predict new population size next year

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

Complex coding

We’re basically simulating what we did in the previous code chunk, only now we’re doing it for more years/making more predictions

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

#1998-99
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-01
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2001 <- lambda_rand.t*N.2000

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

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

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

#2004-05
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2005 <- lambda_rand.t*N.2004
N.2005
#> [1] 105.7996

Plot change in population

Saving our data each run and plotting it to visually see any population size patterns

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 loops

A few steps closer to using for loops to make predictions

# Initial conditions

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

# Explore xlim = argument 
#xlim and ylim sets the axis

plot(N.1997 ~ c(1997))

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

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

# For loop the hard way
#
N.current <- N.initial

# for loop would start here
t <- 1
  
  # get a lambda from the hat
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # determine pop size
  N.t <- N.current*lambda_rand.t
  
  # for the next years
  year.t <- 1997+t
  
  # plot the pop

  points(N.t ~ year.t)

  
  # update N.current
  N.current <- N.t
  # we've done one simulation from 1997-98 and to keep the process going we should increment t by one and re run this code chunk

Now we want to start making a for loop to mass produce data.


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

N.current <- N.1997

# doing what we did in the code chunk above, except now all at once
for(t in 1:50){
  
  # all at once running this code below
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # each new growth rate
  N.t <- N.current*lambda_rand.t
  
  # each consecutive year 
  year.t <- 1997+t
  
  # plotting the points to save 
  points(N.t ~ year.t)
  
  # new population 
  N.current <- N.t
}

Funky r plotting code

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

Creating 9 simulations of what we just did, we see that there are very distinct predictions.

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

Kaylee Wanamaker

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 lambda sampling

Getting inital population size

One round of population simulated

Complex coding

Plot change in population

For loops