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
census <- c(1:39)
census <- seq(1, 39)
census <- seq(1, 39, by = 1)

year t

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

year.t   <- 1959:1997
year.t   <- seq(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 available 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
#the above two lines of code are accessing the population sizes starting from and stopping at certain points in the vector  

females.N[2:3]/females.N[1:2]
#> [1] 1.0681818 0.9787234
#this is dividing the first value in each vector and the second. so 47/44 and 46/47

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
#since females.N is not a matrix you can get the dimensions of the vector by using the length function
#the length is stored in len and that is used as the stopping point(max value) in the sequence. The values in the  full length of the vector starting at position 2 are divided by the values of the vector starting at position 1 going until 1 less the total length

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
#Same operation performed as previous one. Here the length is used directly instead of being assigned to a variable

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[seq(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
#the output showed the vector without the first value. [-1] dropped that first value

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

“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)]

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

lambda.x <- c(lambda.i, NA)
length(lambda.x)
#> [1] 39
lambda.x
#>  [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
# NA- not available?, creates a vector with lambda.i and NA, to lengthen lambda so it is 39 elements.

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.x)
lambda_lug <- log(lambda.x)
length(lambda_lug)
#> [1] 39
#computes natural logarithms - log base 10

Assemble the dataframe

bear_N <- data.frame(census,year.t,
                females.N,
                lambda.x, 
                lambda_lug)

TASK

List 3 functions that allow you to examine this dataframe.

dim()
summary()
str()

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

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 = "Population index (females + cubs)",
     xlab = "Year")
abline(v = 1970)

Plot lambda

CHALLENGE TASK

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

hist(bear_N$lambda.x)

Plot ln(lambda)

CHALLENGE TASK

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

hist(bear_N$lambda_lug)

Mean of log(lambda)

CHALLENGE TASK

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

mean(bear_N$lambda_lu, na.rm = T)
#> [1] 0.02134027
#outputs 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_lug, na.rm = T)

CHALLENGE TASK Make a histogram with the mean plotted on it

hist(bear_N$lambda_lug)
abline(v = u)

Density dependence

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

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

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

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

#where bear population is
plot(females.N ~ year.t, data = bear_N, 
     type = "b",
     ylab = "Population index (females + cubs)",
     xlab = "Year")
# y = bx + a
abline(v = 1970)
abline( v = 1987, col = "red")

#population in yellowstone didn't have a chance to persist due to high variability, time stopped at red line
#many bears before black line, were eating poorly and getting close to people (trash, campgrounds)
#after red line, population continued to expand, doe pva and see if hypothesis is true

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

ADD 1-2 sentences note here simulate population of bears and preduct in future what will happen with that population, variations occur, don’t know how conditions will play out. run multiple simulations in the end - don’t know how random vairaitions will play out. take key parameters and characterize how variable they are and draw them out randomly to model population dynamics we randomly pull population growth rates out of a hat. simulation -

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

#what is the risk that the population will die out, how do you best characterize what will happen in the future
#key parameter is lambda - population growth rate , 37 ex. of how population grows, large distribution, mean and median - around 1, popuation variation- from .8 to 1.3
#draw numbers bear_N$lambda.i randomly

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
#Na is present in dataset- empty slot in vector, NA, NaN, Inf- math operation that doesn't result in real number, -Inf, Null
#is.na- looks at every value in lambda and returns T/F if NA is present

any(is.na(hat_of_lambdas) == TRUE) 
#> [1] TRUE
#check are there any values that are equal to 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)] #another way to remove NA, better way, more general form
#>  [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)# drop NA in a single function - na.omit
#>  [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)]
#make a new vector without that NA value

hist(hat_of_lambdas)

ADD TITLE HERE Perform a random sample, choose a random value from the parameter, random sampling of lambdas

ADD 1-2 sentences note here

# choose random value from hat via sample function
sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#> [1] 1
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE) # pull lambda from hat and save
#population will grow at x, then pull another lambda
#start with same hat
#change size to 50 - didn't put 50 values in hat, only put 37,values will start to recycle
#replace = TRUE = recycle = TRUE
#if replace is set to FALSE- in sample is larger than population - error, since values cannot be recycled

get initial population size

ADD 1-2 sentences note here initialize simulation for population head, tail - functions that allow you to view dataframes

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 #initial population size,

one round of population simulation

ADD 1-2 sentences note here pull 1.22807 randomly from hat- poplation growth rate
in 1997, want to predict number of bears in 1998

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

Simulation the hard way

ADD 1-2 sentences note here like excel, but worse…
99 bears in 1997 - simulate hypothetical scenario- how many bears in 50 years taek population of bears, random lambda- multiply both and see how many bears you get. take that number, pull random lambda, multiply them and get the next year’s population

modeling how a population could decline ignoring carrying capacity- keeps population from getting to high

#1997  to 1998
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.1998 <- lambda_rand.t*N.1997 #multiply random lambda by population size to get next year's population
#start with 113

#1989 to 1999
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 #184

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 #197 - peak

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2005 <- lambda_rand.t*N.2004 #get 166 as my population, overall growth, slight declines in between years

ADD TITLE HERE Plotting of simulation data

ADD 1-2 sentences note here In my particular plot: see a sharp increase until 2000, decreases until 2002, then increases until 2003, decrease in 2004

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

ADD TITLE HERE Prelimary for loop stes

ADD 1-2 sentences note here

# Initial conditions

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

# Explore xlim = argument
# xlim - sets size of x axis
plot(N.1997 ~ c(1997)) #plot single point

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50)) #orient point to the left

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

# for loop, the hard way
#
N.current <- N.initial #have current population start at initial

#  where for loop would be
t <- 1 #start at year 1 of 50 year simulation, incrementer
  
  # Sample random lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # ADD COMMENT HERE determine population size - 99 multiplied by lambda to get new population size
  N.t <- N.current*lambda_rand.t
   
  # ADD COMMENT HERE go to the next year
  year.t <- 1997+t
  
  # ADD COMMENT HERE plot that point (1998)
  #points() updates existing graph
  points(N.t ~ year.t)

  
  # ADD COMMENT HERE update N.current, have gone through 1 year of simulation
  N.current <- N.t
  
  #to go the next year, change t to 2 and then run code again to get 3rd data point

ADD 1-2 sentences note here For loop


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

#start at 1997 population size
N.current <- N.1997

# For loop
for(t in 1:50) {#instead of having to induvidually increment, starting at 1 and going to 50- each time, add to t and run inside code
  
  # ADD COMMENT HERE
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # ADD COMMENT HERE
  N.t <- N.current*lambda_rand.t
  
  # ADD COMMENT HERE
  year.t <- 1997+t
  
  # ADD COMMENT HERE
  points(N.t ~ year.t)
  
  # ADD COMMENT HERE
  N.current <- N.t
}

#graph in screen shot- logistic graph, have generated density dependence and logistic population growth
#pattern in data doesn't mean process in biology

ADD 1-2 sentences note here plotting code/magic declares parameters for plotting 3 by 3 panel of graphs

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

ADD 1-2 sentences note here keep track of 9 simulations to see variability

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 expanding parameters

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

Add 1-2 sentence of notes here add red line to keep track of each simulation still need to click for each simulation


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 nested for loop

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.000 81.675

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

Computational Ecology: Population Viability Analysis 9/15/2020

Nathan Brouwer

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?

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

ADD TITLE HERE Perform a random sample, choose a random value from the parameter, random sampling of lambdas

get initial population size

one round of population simulation

Simulation the hard way

ADD TITLE HERE Plotting of simulation data

ADD TITLE HERE Prelimary for loop stes

ADD TITLE HERE

ADD TITLE HERE nested for loop

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

Shortened code

Parameteric sampling