Computational Ecology: The Bears of Yellowstone National Park

Introduction

This document introduces a classic dataset used in computational population ecology - the population dynamics of Grizzly Bears in Yellowstone National Park. That data were first published in Figure 5 of Dennis et al (1991), who extracted the data from other sources. Morris and Doak (2002) then added additional years of data. I will build up the dataset piece by piece and explain the R code used to create all the parts.

Preliminaries

First, we’ll build up individual vectors that will make up the dataframe. Then we’ll calculate the key data we are considering: the population growth rate (lambda) of the population.

Building up the dataframe

Census period

The census period; an index from 1 to 39 of how many years of data have been collected.

This can be build several different ways. Each of these is equivalent:

census <- 1:39
census <- c(1:39)
census <- seq(1, 39)
census <- seq(1, 39, 1)
census <- seq(1, 39, by = 1)

year t

The actual year: 1959 to 1997 (Dennis et al use 1959-1987). This again can be build several different ways.

year.t   <- 1959:1997
census   <- c(1959:1997)
year.t   <- seq(1959, 1997)
year.t   <- seq(1959, 1997, 1)
year.t   <- seq(1959, 1997, by = 1)

Population size (N)

Population size is recorded as the number of females with cubs. This does not represent the actual full population size, but is correlated with it.

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: hard-coding example

Population growth rate (lambda) can be calculated as a ratio between the population size in two subsequent years.

Enter the population size for two subsequent years

## Population size in 1959
### The first year of the study
females.N.1959 <- 44

## Population size in 1960
### The second year of the study
females.N.1960 <- 47

Calculate the ratio of the 2 subsequent population sizes. This is lambda.

lambda.59_60 <- females.N.1960/females.N.1959
lambda.59_60
#> [1] 1.068182

What we just did gives us lambda, but it requires hard coding the information.

Population growth rate: a slightly better way

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

# Check:
## Access the data
females.N[1]
#> [1] 44
females.N[2]
#> [1] 47

# store data 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

females.N.1960/females.N.1959 == lambda.59_60
#> [1] TRUE

Population growth rate: the best way

Calculate lambda using bracket notation and without creating intermediary object.

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
# can't be done; there is no date for 1958
# lambda.58_59 <- females.N[1]/females.N[ ]

Population growth rate: vectorized

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

This code is calculating two lambdas, from the first to the second year, and from the second to the third year.

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

# the denominator
females.N[1:2]
#> [1] 44 47

# two lambdas
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 from all of the years. Note that I am hard-coding the length of the vector.

# how many data points do we have?
length(females.N)
#> [1] 39

# the numerators:
females.N[2:39]
#>  [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

# the denominators
females.N[1:38]
#>  [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

# the lambdas
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 next line of code do? Briefly describe in 1 to 2 sentences why I am using length().

Instead of hard coding the number of years of data (39), this code calculates the length the vector females.N and stores it as the object len. Then instead of using the number 39 in the code, len is used. Instead of 38, len-1 is used.

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.

In this code, instead of creating an object to hold the length, everything is done on a single line of 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: the first 10 years of data.

females.N[1:10]
#>  [1] 44 47 46 44 46 45 46 40 39 39
females.N[seq(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) # length = 10 : )
#> [1] 10

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

females.Ntemp[-1] prints 9 numbers. It is dropping off the first element of the vector.

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?

With 10 years of data, you can calculate 9 lambdas.

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. This is a useful way to avoid hard coding.

TASK: Drop the the first element of this vector

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

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

TASK: Drop the second element of this vector

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

#answer
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.

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

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

TASK: How do you access the last element of the vector? Do this in a general way without hard-coding.

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

#answer
## length gives us the total length
## which is the same as the index of the last
## element of the vector
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.

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

# answer - note the "-" in front of length
females.Ntemp[-length(females.Ntemp)]
#> [1] 44 47 46 44 46 45 46 40 39

TASK: Calculate the first 9 lambdas in the vector females.Ntemp

# All lambdas: 44 47 46 44 46 45 46 40 39 39

# Drop 1st    [-1]47 46 44 46 45 46 40 39 39 

# Drop last    44 47 46 44 46 45 46 40 39 [-10]

# Numerator        47 46 44 46 45 46 40 39 39   
#                 /  /  / ...
# Denominator   44 47 46 44 46 45 46 40 39 

# 47/44 is 1st lambda
# 46/47 is 2nd
# 44/46 is 3rd

# Answer
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?

Drops the first element of the vector

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

Drops the last element of the vector

TASK: Calculate lambdas for all of the data

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

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

# all lambdas
lambda.i <- females.N[-1]/females.N[-length(females.N)]

Finish putting together dataframe

We have made 3 things

an index of the survey, 1 to 39
the raw time series of counts of the bears
a vector of lambda

Now we’ll put things together.

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

There are 39 years of counts of bears. Lambdas is ratio of years and represents a transition between years. For 39 years, there are 38 transitions. For example, if we had 3 years of data with these counts

10, 12, 34

We’d have 2 transitions

10 to 12
12 to 34

If we had 10 years of data, we’d have 9 transitions.

So, the vector of the raw data (counts) is longer than the vector of lambdas

length(lambda.i) == length(females.N)
#> [1] FALSE

I want to put all my data - counts (N) and lambdas into a single dataframe. So I need to add an extra value to my lambdas to make them the same length. Otherwise, R throws an error.

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

TASK: Check the help file; what type of log does log() calculate

# the natural log, not log base 10
lambda_log <- log(lambda.i)

#log base 10 would be
log10(lambda.i)
#>  [1]  0.028645181 -0.009340026 -0.019305155  0.019305155 -0.009545318
#>  [6]  0.009545318 -0.060697840 -0.010995384  0.000000000  0.032184683
#> [11] -0.032184683  0.021719250 -0.010723865 -0.083546051  0.037788561
#> [16] -0.024823584  0.059585690 -0.046996563 -0.012589127  0.048304680
#> [21] -0.023481096  0.011899223  0.044582133 -0.021719250  0.116505569
#> [26] -0.035472318  0.083776998 -0.074633618  0.096910013  0.034762106
#> [31]  0.056318363 -0.030382629 -0.025935734 -0.057038501  0.089223184
#> [36]  0.063386979  0.087150176  0.000000000           NA

Assemble the dataframe

Put all the pieces together using data.frame()

bear_N <- data.frame(census,
                year.t,
                females.N,
                lambda.i, 
                lambda_log)

TASK: List at least 3 functions that allow you to examine this dataframe.

dim()
nrow()
ncol()
summary()
head()
tail()
is()
class()

Examine 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,                         # y ~ x
     data = bear_N,                              #data = the dataframe
     type = "b",                                 # b = plots both points and lines
     ylab = "Population index (females + cubs)", #ylab = y axis label
     xlab = "Year")                              #xlab = x axis label

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

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.

This is done with abline()

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

Plot lambda

CHALLENGE TASK

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

Histograms are made with hist(). These data are normal-ish, but trail off to the right (are skewed).

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.

Sometimes taking the log of things makes them more normal. Not so much here.

hist(bear_N$lambda_log)

hist(log(bear_N$lambda.i))

Mean of log(lambda)

CHALLENGE TASK: Briefly describe what happens when you delete na.rm = T

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

#this just gives you an NA
mean(bear_N$lambda_lo)
#> [1] NA

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.

This was not covered in lecture and is just here for reference. Density dependence occurs when a vital rate of an organism (survival, birth rate, lambda) changes with the density of an organism. Often, when density gets high, vital rates go down.

In these data, there doesn’t seem to be an density dependence.

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

Computational Ecology: The Bears of Yellowstone National Park

Nathan Brouwer

9/15/2020

Introduction

Preliminaries

Building up the dataframe

Census period

year t

Population size (N)

Population growth rate: hard-coding example

Population growth rate: a slightly better way

Population growth rate: the best way

Population growth rate: vectorized

Negative indexing

Calculate lambdas for all data

Finish putting together dataframe

Create special columns

Assemble the dataframe

Examine the population growth rates

Plotting the raw data

Challenge questions

Plot lambda

Plot ln(lambda)

Mean of log(lambda)

Histogram with mean

Density dependence