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[]
# No data for 1958

Population growth rate: vectorized

TASK

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

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

The above code is finding lambda values for 59-60 and 60-61.

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

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

Gets all of the lambda values for each year in the data set. You are using length so that if more data is added to the set (more years, for example), you do not need to change the code as it already accounts for the whole data set.

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


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

It is the same, except this code chunk does not store the length in a variable before doing the calculation, it just uses the formula in place of len.

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.

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

It is getting rid of the first number in the array.

TASK How many lambdas can I calculate using the first 10 years of data?

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

9 lambdas.

“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[2:10]/females.Ntemp[1:9]

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

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

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

Not sure what it does, but I would guess that it creeates a blannk columnn?

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)

TASK

List 3 functions that allow you to examine this dataframe.

dim()
summary()
colnames()

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="Pop index (females+club)", xlab="Yr")

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="Pop index (females+club)", xlab="Yr")
abline(v=1970)

Plot lambda

CHALLENGE TASK

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

Plot ln(lambda)

CHALLENGE TASK

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

Mean of log(lambda)

CHALLENGE TASK

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

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

CHALLENGE TASK Make a histogram with the mean plotted on it

Density dependence

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

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

We are going to randomgy generate population growth rates.

hat_of_lambdas <- bear_N$lambda.i

Determine if something is 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

Get rid of the NA value

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)

Sample One Value

This code chunk samples a single value from a theoretical ‘hat’ and then saves it to an object.

#pulling from hat
sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#> [1] 1.25
#Saving value
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)

Looking at Data frames

These are a few functions you can use to evaluate the content of data frames. Head gives you the first x values and tail gives you the last x values.

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

We arer multiplying the lambda by the first sampled population size to estimate the population size of the next year

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

Simulation the ‘hard way’

Every 2 lines represents sampling the lambda value, multiplying it by the population of the previous year to estimate the next year’s population. This is not a very efficient way because we are just hard-coding the vaules.

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

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

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2000 <- lambda_rand.t*N.1999

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

Plotting population change

Plots population change per year

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

Plotting with for loops

This is plotting points of population depending on the year using a pseudo for() loop method.

# Initial conditions

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

# making a new plot
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

# at this point, this is where the for loop would be
t <- 1
  
  # sample lamda
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # get population size like we did above
  N.t <- N.current*lambda_rand.t
  
  # increment the year
  year.t <- 1997+t
  
  # plot population
  points(N.t ~ year.t)

  
  # update Ncurrent to 'loop'
  N.current <- N.t

This is a repeat of the above^ code except this time we are using a traditional for loop.


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

#initial conditions (1997s pop size)
N.current <- N.1997

# for loop header to loop t from 1 to 50
for(t in 1:50){
  
  # sample rrandom lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # estimate pop size
  N.t <- N.current*lambda_rand.t
  
  # increement year
  year.t <- 1997+t
  
  # plot the point
  points(N.t ~ year.t)
  
  # update Ncurrent
  N.current <- N.t
}

Format area for multi planel plot support (possibly), mostly just goofy code

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

END HERE BECAUSE THIS IS WHERE I MADE THE MULTI PANEL PLOT

————————————————–

ADD 1-2 sentences note here

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
}

ADD TITLE HERE

Add 1-2 sentence of notes here

par(mfrow = c(10,10), mar = c(0,0,0,0), xaxt = "n", yaxt = "n")

Add 1-2 sentence of notes here


plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550), cex = 0.5)
abline(h = N.1997)
abline(h = 0, col = "red")
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, cex = 0.5)
  
  N.current <- N.t
}

ADD TITLE HERE

Add 1-2 sentence of notes here

## ADD TITLE HERE
num.years <- 50

## ADD TITLE HERE
N.vector <- rep(NA, num.years)

## ADD TITLE HERE
N.vector[1] <- N.1997

t<- 1
  # print the year
  print(t)
#> [1] 1
  
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t

ADD TITLE HERE

Add 1-2 sentence of notes here

## ADD TITLE HERE
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD TITLE HERE
#don't uncomment print() in notebook mode!
for(t in 1:num.years){
  # print the year
  #print(t)
  
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t
  
  
}

ADD TITLE HERE

par(mfrow = c(1,1), mar = c(3,3,3,3), xaxt = "s", yaxt = "s")

df <- data.frame(N.vector, year = seq(1997, 1997+50))

plot(N.vector ~ year, data = df)


plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

ADD TITLE HERE

ADD 2-3 sentences of notes here

## ADD TITLE
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD TITLE
lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)


## ADD TITLE
t <- 1


  print(t)
#> [1] 1
  
  lambda_rand.t <- lambdas_vector[t]
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t
  
  t <- t+1
  N.vector[1:t]
#> [1]  99.00 107.25

ADD TITLE HERE

## ADD NOTES
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD NOTES
lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)


## ADD NOTES
for(t in 1:num.years){
  
  ## ADD NOTES
  lambda_rand.t <- lambdas_vector[t]
  
  ## ADD NOTES
  N.t <- N.vector[t]*lambda_rand.t
  
  ## ADD NOTES
  N.vector[t+1] <- N.t
  
}

## ADD NOTES
df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

How can I shorten the code?

lambdas_vector[t]
N.t

num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)

for(t in 1:num.years){
  lambda_rand.t <- lambdas_vector[t]
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t
  
}

df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

Shortened code

num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)

for(t in 1:num.years){
  N.vector[t+1] <- N.vector[t]*lambdas_vector[t]
  
}

df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

Portfolio 1: PVA Assignment

Jason Walker

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?

Sample One Value

Looking at Data frames

One round of population simulate

Simulation the ‘hard way’

Plotting population change

Plotting with for loops

END HERE BECAUSE THIS IS WHERE I MADE THE MULTI PANEL PLOT

————————————————–

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

Shortened code

Parameteric sampling