Bear Population Growth
Theresa Hutchins
9/25/2020
#the census period
census <- 1:39
census <- c(1:39)
census <- seq(1, 39)
census <- seq(1, 39, by = 1)
#1959-1997
year.t <- 1959:1997
#population size recorded as number of females with cubs
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)
#females.N.1959 <- 44
#females.N.1960 <- 47
#access the data
#females.N[1] #44
#females.N[2] #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.068182
#pop growth rate: vectorized
#females.N[2:3] #47 46
#females.N[1:2] # 44 47
#females.N[2:3]/females.N[1:2] #1.07 0.98
#length(females.N) #39
#females.N[2:39]/females.N[1:38]
#len <- length(females.N) #length is used so the proper amount of indexes are calculated
#females.N[2:len]/females.N[1:len-1]
#calculates lambda for the endtire year span receorded
#females.N[2:length(females.N)]/females.N[1:length(females.N)-1]
##
##Negative indexing
#females.N[1:10]
#females.N[seq(1,10)]
#females.Ntemp <- females.N[seq(1,10)]
#length(females.Ntemp) #10
#females.Ntemp[-1] #returns all indexs besides the first one
#females.Ntemp[2:10]/females.Ntemp[1:9] #9 lambdas returned
#females.Ntemp[-1] #drop first elem
#females.Ntemp[-2] #drops second element
#females.Ntemp[-10] #drops 10th elem
#females.Ntemp[length(females.Ntemp)] #accessing the last element
#females.Ntemp[-length(females.Ntemp)] #drops last element
#calculate first 9 lambdas
#calc lambda for the ratio of two pop sizes
#lambda.59_60 <- females.N.1960/females.N.1959 #bracket notation
#lambda.59_60 <- females.N[2]/females.N[1]
#correct numbered lambdas
lambda.i <- females.N[-1]/females.N[-length(females.N)]
lambda.i <- c(lambda.i,NA)
lambda_log <- log(lambda.i)
#assembling the dataframw
#can get min, max and the summary
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)
How do we determine if a population is likely to go extinct?
hat_of_lambdas- want to predict the future of when and what will happen to the bear population. We want to simulate a bear population, and predict what will happen in the future. Take key parameters, apply how variable the rates are and spit out random predictions.In this case we are interested in pop growth rate.
hist(bear_N$lambda.i)
#right skewed histogram of log
hist(bear_N$lambda_log)
#mean of log(lambda)
u <- mean(bear_N$lambda_lo, na.rm = T) #.02134
hist(bear_N$lambda_log)
abline(v = u) #mean plotted on graph
#density dependance
plot(lambda.i ~ females.N, data = bear_N, type = "p",
ylab = "Population growth Rate", xlab = "Females with Cubs")
abline(v = 1970)
#pop growth rates
hat_of_lambdas <- bear_N$lambda.i#gives the set of logical statements determining if an index 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 returns if there are any values equal to true
any(is.na(hat_of_lambdas) == TRUE)
#> [1] TRUEDropping 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)] #drop last index
#> [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#allows to drop the NAs
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)
Random Sampling of Lambdas
Every time this code is run a new number is chosen. Starts form same hat every time. If the size increases for the limited years of data numbers are recycled (if replace = true).
# chooses a random lambda from the "hat"
sample(x = hat_of_lambdas, size = 1,replace = TRUE)
#> [1] 1.222222
#pull and save to object, dif value every time
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
lambda_rand.t
#> [1] 0.9756098Initializing the populations.
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 <- 99One round of population simulation
Predicts the amount of bears for the given year.
1.22807*99
#> [1] 121.5789
lambda_rand.t*N.1997
#> [1] 96.58537
N.1998 <- lambda_rand.t*N.1997Simulation the Hard Way
Want to simulate a hypothetical future using lambda and multiplying them.
#pulls random lamdba, multiplied by pop size that year
#1997-1998
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.1998 <- lambda_rand.t*N.1997
#1998-1999
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.1999 <- lambda_rand.t*N.1998
#1999-2000
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2000 <- lambda_rand.t*N.1999
#200-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.2004PLotting future random populations
Creating a graph of years 1997-2004 of randimly generated poulations.
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 individual points
Will be using a for loop to take out some laborous coding.
# Initial conditions
N.1997 <- 99
N.initial <- 99
# explore xlim = argument , sets size of xaxis
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))
# the not a for() loop
#for loop the hard way
N.current <- N.initial
# this is where the foor loop would be
t <- 1
# grab a lambda
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace =TRUE)
# determine pop size
N.t <- N.current*lambda_rand.t
# 1997 + t = 1998 poop growth
year.t <- 1997+t
#plot the population
#points() update an existing graph
points(N.t ~ year.t)
# update N.current
N.current <- N.tUsing a for loop to generate future lambdas to predict future growth. A pattern in data doesnt mean a process in biology.
# make a new plot
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
#start at 1997 pop size
N.current <- N.1997
# for loop, avoids individually incrementing the years
for(t in 1:50){
# choses random lambda
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
}
goofy r plotting code 3 by 3 panel of graphs to see variability
par(mfrow = c(3,3), mar = c(1,1,1,1))Running more 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
}
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.000 81.675ADD 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)))