Population_Viability

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

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 the change in growth rate from one year to the next.

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[ERROR]
# It is not possible to calculate a lambda for this since the vector of years does not include 1958.

Population growth rate: vectorized

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

The code is creating a small vector of values from the original vector females.N, and then dividing the two vectors.

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

Length() does the same thing as the above code but rather than hardcoding the value 39 the value is stored in a variable called 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

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

Instead of the length of the vector being stored as a variable, the function length() is being called directly into the equation. It will do the same thing as the above code.


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.

[-1] subtracts the first value of the vector. This is known as negative indexing. I don’t think this is a thing in java as far as I know.

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?

You can calculate 9 lambdas from 10 years, since each lambda needs 2 years at least.

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.

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)]
lambda.i
#> [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#> [8] 0.9750000 1.0000000

Calculating 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? Subtracts the first value of the vector females.N
What does females.N[-length(females.N)]? Subtracts the last value of the vector 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)]
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.)

This code creates a vector with lambda.i. I assume the NA is to add another list named NA or there is no other list.

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!) Calculates natural log.

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.

Summary()
Mean()
Median()

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.

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

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, col = "red")

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

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

We determine if a population is likely to go extinct if the species shows a general downward trend or negative growth rate. We randomly pull population growth rates out of a “hat” to build data about population growth.

hat_of_lambdas <- bear_N$lambda.i
hist(bear_N$lambda.i)

NA NAN -INF INF NULL

Above is basically if the data is null is.na() == true or false if it is NA/null

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)

Data Sampling Random Growth Rates

We are now sampling a random number from our “hat” hat_of_lambdas vector. We take a random value from the data frame and then save it into a variable.

#Sampling with hat_of_lambdas, a vector of growth rate values. Sampling with replacement; pulling a random value from the vector
#Replace = TRUE means that you can recycle values if size is greater than the vector length
sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#> [1] 1.027778

#Saves random value into a variable
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
lambda_rand.t
#> [1] 1.090909

Statistics of the Data Framem Bear_N

These allow us to look at statistics of our data frame. Head looks at the first 6 values, tail looks at the last 6 values, summary looks at stats such as mean, median, etc.

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
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
N.1997 <- 99

One Round of Population Simulation

This is a prediction of how many bears there will be the next year given the number of bears the previous year and a growth rate (lambda).

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

Simulating the Data Poorly

This is a much more tedious way to simulate data than the above code. We are essentially hard-coding a simulation.

#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

Plotting the Data Above

This is a plot of the theoretical simulation that we ran above.

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

Simulating the Data Using a For Loop

Rather than do it in a type of “spreadsheet” type of way, which is rather inefficient, an easier way to do it is by using for loops. This only plots a few points at a time because t is equal to only 1.

# Starting year

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

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

# Hard coded For loop
N.current <- N.initial

# For looping with 1 value
t <- 1
  
  # Grabbing a lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # Population size
  N.t <- N.current*lambda_rand.t
  
  # Plot population
  year.t <- 1997+t
  
  # Creating a coordinate
  points(N.t ~ year.t)

  
  # Updating the current N
  N.current <- N.t

The proper way to for loop a simulation, for looping 50 samples for the plot.


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

N.current <- N.1997

# For looping where t repeats itself incrementally
for(t in 1:50){
  
  # Grabbing a random lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # Determine population size
  N.t <- N.current*lambda_rand.t
  
  # Plotting the population
  year.t <- 1997+t
  
  # Plots coordinate
  points(N.t ~ year.t)
  
  # Update N
  N.current <- N.t
}

Goofy R-plottin - this code allows us to create multiple graphs to view simulations.

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

This runs 50 samples of lambda and then plots it; best example of for looping data.

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
}

Population_Viability

Thomas

9/20/2020

Data: Figure 5 (plus 10 years)

census

year t

Population size

Population growth rate: example

Population growth rate: vectorized

Negative indexing

Calculating 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?

Data Sampling Random Growth Rates

Statistics of the Data Framem Bear_N

One Round of Population Simulation

Simulating the Data Poorly

Plotting the Data Above

Simulating the Data Using a For Loop