Computational Ecology: Population Viability Analysis 9/15/2020
Jacob Krasner
9/15/2020
Data: Figure 5 (plus 10 years)
Make data vectors, calculate lambda, and put together datafram with all necessary dataa
Census
The ccensus period; an index from 1 to 39 of how many years of data have been collected
census <- 1:39Year t
The year: 1959 to 1997(Dennis et al use 1959 - 1987)
year.t <- 1959:1997Population Size
Population size is reccorded as the number of females wih…
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…
females.N.1959 <- 44
females.N.1960 <- 47lambda.59_60 <- females.N.1960/females.N.1959
lambda.59_60
#> [1] 1.068182
(47-44)/44
#> [1] 0.06818182females.N[2] #equivalent to 1960 because second vaalue in vector
#> [1] 47
females.N[1] #equivalent to 1959 because first value in vector
#> [1] 44
females.N.1959 <- females.N[1] #we can store thhose into the yea
females.N.1960 <- females.N[2] #vectors
females.N.1960 / females.N.1959
#> [1] 1.068182The firs year of data is 1959. What is lambda for 1958 o 1959 Even though there are 39 years, we can only calculate 38 lambdas because we calculate based off of difference between two years
#the below three are equivalent, all print out the second and third values
females.N[2:3]
#> [1] 47 46
females.N[c(2,3)]
#> [1] 47 46
females.N[c(2:3)]
#> [1] 47 46
#Prints ratio of second and third lambda over the first and second lambdas
females.N[2:3] / females.N[1:2]
#> [1] 1.0681818 0.9787234Calculate all lambdas at once
length(females.N)
#> [1] 39
#compares higher and lower for each value
#Starts with 2 and divides with 1
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###What does this do
#This further simplifies because we dont even need to know how long something is
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.0000000females.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]
females.Ntemp <- females.N[seq(1,10)]females.Ntemp
#> [1] 44 47 46 44 46 45 46 40 39 39females.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.0000000Negative Indexing allows you to drop a specific element from a vector
The first element
females.Ntemp[-1]
#> [1] 47 46 44 46 45 46 40 39 39The second element
females.Ntemp[-2]
#> [1] 44 46 44 46 45 46 40 39 39Drops the 10th element
females.Ntemp[-10]
#> [1] 44 47 46 44 46 45 46 40 39females.Ntemp[-c(1:2)]
#> [1] 46 44 46 45 46 40 39 39
females.Ntemp[-c(1,2)]
#> [1] 46 44 46 45 46 40 39 39
females.Ntemp[-1:-2]
#> [1] 46 44 46 45 46 40 39 39####Access and remove the last element
females.Ntemp[length(females.Ntemp)]
#> [1] 39
females.Ntemp[-length(females.Ntemp)]
#> [1] 44 47 46 44 46 45 46 40 39Calculate the first 9 lambdas
lambda.i <- females.Ntemp[2:(2+8)] / females.Ntemp[1:(1+8)]
lambda.i <- females.Ntemp[-1] / females.Ntemp[-10]
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.0000000census <- 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)
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
#> 7 7 1965 46 0.8695652 -0.13976194
#> 8 8 1966 40 0.9750000 -0.02531781
#> 9 9 1967 39 1.0000000 0.00000000
#> 10 10 1968 39 1.0769231 0.07410797
#> 11 11 1969 42 0.9285714 -0.07410797
#> 12 12 1970 39 1.0512821 0.05001042
#> 13 13 1971 41 0.9756098 -0.02469261
#> 14 14 1972 40 0.8250000 -0.19237189
#> 15 15 1973 33 1.0909091 0.08701138
#> 16 16 1974 36 0.9444444 -0.05715841
#> 17 17 1975 34 1.1470588 0.13720112
#> 18 18 1976 39 0.8974359 -0.10821358
#> 19 19 1977 35 0.9714286 -0.02898754
#> 20 20 1978 34 1.1176471 0.11122564
#> 21 21 1979 38 0.9473684 -0.05406722
#> 22 22 1980 36 1.0277778 0.02739897
#> 23 23 1981 37 1.1081081 0.10265415
#> 24 24 1982 41 0.9512195 -0.05001042
#> 25 25 1983 39 1.3076923 0.26826399
#> 26 26 1984 51 0.9215686 -0.08167803
#> 27 27 1985 47 1.2127660 0.19290367
#> 28 28 1986 57 0.8421053 -0.17185026
#> 29 29 1987 48 1.2500000 0.22314355
#> 30 30 1988 60 1.0833333 0.08004271
#> 31 31 1989 65 1.1384615 0.12967782
#> 32 32 1990 74 0.9324324 -0.06995859
#> 33 33 1991 69 0.9420290 -0.05971923
#> 34 34 1992 65 0.8769231 -0.13133600
#> 35 35 1993 57 1.2280702 0.20544397
#> 36 36 1994 70 1.1571429 0.14595391
#> 37 37 1995 81 1.2222222 0.20067070
#> 38 38 1996 99 1.0000000 0.00000000
#> 39 39 1997 99 NA NAplot(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?
Below is an example of a for loop
How do we best characterize what will happen in the future
Throw numbers into a hat and draw them out randomly
Using this randoness created by a hat, we can momdel population dynamics
hist(bear_N$lambda.i)
hat_of_lambdas <- bear_N$lambda.iNote that NA is missing data or not applicable so the NA values are simply empty
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 NAThis checks if each value within hat_of_lambdas is na or not… False if not
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 TRUEChecks if any of the values that we have in the array NA
any(is.na(hat_of_lambdas) == TRUE)
#> [1] TRUEDrop the NA Convoluted way of doing the below
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.0000000Omits all NA’s from the array
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"This is a vector contining the last value of the array
hat_of_lambdas <- hat_of_lambdas[-length(hat_of_lambdas)]hist(hat_of_lambdas)
Sampling and Storing from hat of lambdas
Pulled lambda from hat usnig sample Saved it to a function Note that replace is TRUE so we put it back but we start over no memory
#pull one lambda from hat and return it
sample(x = hat_of_lambdas,
size = 1,
replace = TRUE)
#> [1] 0.975
#run again but this time store it... runs independent of above
lambda_rand.t <- sample(x = hat_of_lambdas,
size = 1,
replace = TRUE)##Initial Population Size
Set up for population simulation Getting top (6) and bottom (6) of the data frame
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##One Round of population simulaton
Note that generally you should not hard code Pull random value (random growth rate)
#random growth rate and population in a certain year
1.22807*99
#> [1] 121.5789
#random growth rate * # of bears in a certain year
lambda_rand.t*N.1997
#> [1] 93.78947
#store into next years value of bears
N.1998 <- lambda_rand.t*N.1997##Simulation the Hard way
Using this method has us rewrite for each simulation and is highly inneficienty
Shows trajectory based on random variation
#1997 - 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
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.2004Plot
Set the range of years you want to go through and put this into years
Then get the number of bears from 1998 - 2005 and store in a vector in N.rand
Then we make a data frame with comparison of the number of bears and the year (year on the x)
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")
##For Loop with above data
We are doing 1 year of simulation at a time … but in this scenario, we have to change the value of t, or years passed, each time
#Initial Conditions
N.1997 <- 99
N.initial <- 99
#Explore xlim = argument -- sets size of x axis
plot(N.1997 ~ c(1997))
plot(N.1997 ~ c(1997),
xlim = c(1997, 1997+50))
#sets x limits and y limits
plot(N.1997 ~ c(1997),
xlim = c(1997, 1997+50),
ylim = c(0, 550))
#difficult for loop
#
N.current <- N.initial
# shouldnt usually use t because its actually a function but this is the staryt of the loop
t <- 1
#get a laambda
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#determien population size
N.t <- N.current*lambda_rand.t
#add 1 to the year through t
year.t <- 1997+t
# plot the point
points(N.t ~ year.t)
# update the current value
N.current <- N.tADD 1-2 sentences note here
# Create a plot with the below parameters and set the current # of bears into the current value
plot(N.1997 ~ c(1997),
xlim = c(1997, 1997+50),
ylim = c(0, 550))
N.current <- N.1997
#This is the for loop, do it while t is ccontained int he sequence 1:50
for(t in 1:50){
# get a lmabda
lambda_rand.t <- sample(x = hat_of_lambdas,
size = 1,
replace = TRUE)
# determine population size
N.t <- N.current*lambda_rand.t
# add value of yeaars held in t to the initial year to get current year
year.t <- 1997+t
#plot the point on the graph
points(N.t ~ year.t)
# update the current value
N.current <- N.t
}
Declares paramenters for R plottting
par(mfrow = c(3,3),
mar = c(1,1,1,1))Using the above and below code, we can keep track of 9 simulations
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
}
##Further expanding for loops
Allows for simulation collections….still not a nested for loop
par(mfrow = c(10,10), mar = c(0,0,0,0), xaxt = "n", yaxt = "n")We now plot and repeat run as many times as allowed in parameter above
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.0 93.5ADD 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)))