Grizzly Bear Submission

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 <- c(1:39)

year t

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

year.t <- c(1959:1997)

Population size

Population size is recorded as the number of females with offspring

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 rate at which the individuals in a population increases in a given time period, expressed as a fraction of the initial population

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?

#there is no data for 1948, so it cannot be done
females.N[1]
#> [1] 44
lambda.58_59 <- females.N[1]/females.N[ ]

Population growth rate: vectorized

TASK

This code is assigning two years to one vector (1959 and 1960) and two years to another (1960 and 1961), and then dividing the population growth rates of the respective years by each other (1960/1959 and 1961/1960).

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 to the previous code chunk, just using all of the data (no need to describe) # [2:39] same as seq(2,39)

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 This takes the length of females.N, which is the total number of years for which there are data points, and assigns it to len. Then, it divides numbers ranging from the 2nd year to the last year (len) by their respective 1-year-difference points from year 1 to 38 (len-1).

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 This does the same but with one less line of code by using the length() command inside the statement instead of assigning it to a new variable.


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

This removes the first index (44) from females.Ntemp

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

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]

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?

This removes the first element from the vector

What does females.N[-length(females.N)?

This removes the last element from the 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)]

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 makes lambda.i a mutable list with no data in it? I don’t know

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

log() calculates the 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.

head()
tail()
summary()

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)
abline(v = 1987, col = "red")

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

We determine if a population is likely to go extinct by simulating population dynamics with computational tools. We use a key parameter and characterize how variable it is, then randomly select population growth rates.

hat_of_lambdas <- bear_N$lambda.i

is.na tels you true/false is something 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 drops null values

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)

Pulling First Lambda From Hat

No seed is set, so a different lambda is selected every time this chunk is run. Because replacement is true, if you pull more than your sample size, lambdas will be recycled.

# chooses a random lambda w/ replacement
sample(x = hat_of_lambdas, 
       size = 1,           #size can be set to pull a number of lambdas at a time
       replace = TRUE)
#> [1] 1.090909

# chooses a random lambda and saves it to an object
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1, replace = TRUE)

Get Initial Population Size

Gets the population sizes at the beginning (head) and end (tail) of the data.

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
dim(bear_N)
#> [1] 39  5

N.1997 <- 99

Simulation of One Round of Population Growth

Random population growth rate * population in 1997 to extrapolate predicted population in 1998.

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

Simulation… the Hard Way

Multiply a given year’s population by a random lambda for each year to get next year’s simulated population size.

#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

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 of Population Change over 7 years

A plot of one instance of a simulated change in population size over 7 years based on randomly selected lambdas.

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

One year of simulation (a slightly better way)

This is not the most efficient way to simulate because it only does 1 year at a time, but it is better than the previous, hard-coded way.

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

# for() loop the hard way

N.current <- N.initial

# Spot where the for() loop would be
t <- 1
  
  # Grab a random lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # Determine current population size
  N.t <- N.current*lambda_rand.t
  
  # Increment year by 1
  year.t <- 1997+t
  
  # Add new point to plot for population size
  #points() updates an existing graph
  points(N.t ~ year.t)

  
  # Update value for current population size
  N.current <- N.t

Every time you run the simulation, you get a new result. Therefore, it needs to be run many many times to get useful results.

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

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

# increment t from 1 to 50, adding 1 after each round
for(t in 1:50){
  
  # selects a random labmda
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # calculates new population size with current N and the random lambda
  N.t <- N.current*lambda_rand.t
  
  # increments the year
  year.t <- 1997+t
  
  # plots point
  points(N.t ~ year.t)
  
  # assigns newly calculated N to current N
  N.current <- N.t
}

Goofy R plotting code / magic

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

Run this chunk 9 times to fill each of the 9 spots of the grid.

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
}

Grizzly Bear Submission

Brandon Duncan

9/22/2020

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?

Pulling First Lambda From Hat

Get Initial Population Size

Simulation of One Round of Population Growth

Simulation… the Hard Way

Plot of Population Change over 7 years

One year of simulation (a slightly better way)