# EEB 201 homework
# Lloyd-Smith 9/17/21
# Joanna Wu
# In class we saw how to simulate discrete-time models (such as the discrete logistic growth model) in R, using a for loop.
# Now we will work with an alternative model for density-dependent population growth in discrete time, the Ricker model:
# nt+1 =ntexpr1−nt K
# where K > 0 and r can take positive or negative values.
# (a) Write a function that runs the Ricker model, plots the result, and returns the time series as an output.
# At minimum, your function should take all parameter values and initial conditions as input arguments.
# (Hint: this should involve minimal modification of a modeling function we wrote in class.)
# first define parameters
# n = initial population size
# t = max time years
# r = growth rate
# k = carrying capacity
fun.ricker <- function(n0, tmax, r, k){
# input blank parameters
output <- matrix(NA, nrow=1, ncol=tmax+1)
output[1] <- n0
# equation
for(t in 1:tmax){
output[t+1] <- output[t]*exp(r*(1-output[t]/k))
}
# plot
plot(0:tmax, output, xlab = 'Year', ylab = 'Population',
main = paste0('Growth rate = ', r, ', K = ', k))
lines(0:tmax, output)
# return value
return(output)
}
# run function
fun.ricker(n0=100, tmax=10, r=1.05, k=200)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 100 169.0459 198.8756 200.0531 199.9973 200.0001 200 200 200 200 200
# (b) Explore the dynamics of the model. Try to find combinations of parameter values that yield the following patterns:
# • Population decreases to n = 0.
# try low growth rates
fun.ricker(n0=100, tmax=10, r=-2, k=200) # with r=-2, population nears 0 in 3 years
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 100 36.78794 7.192578 1.046002 0.1430495 0.01938736 0.002624302 0.00035517 4.806721e-05 6.505192e-06 8.803821e-07
# slow decay rate over long time
fun.ricker(n0=50, tmax=10, r=-0.5, k=200) # population decays to near 0 in 5 years
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 34.36446 22.71293 14.58096 9.172124 5.692213 3.501984 2.142738 1.306617 0.7950964 0.4832099
# • Population approaches stable equilibrium at n∗ = K, without oscillations.
# lower n0 than K, long time frame, slow growth rate
fun.ricker(n0=50, tmax=40, r=0.5, k=200)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## [1,] 50 72.74957 99.99789 128.4005 153.5695 172.471 184.7589 191.9346 195.8439 197.8894 198.9363 199.466 199.7325 199.8661 199.933 199.9665
## [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## [1,] 199.9832 199.9916 199.9958 199.9979 199.999 199.9995 199.9997 199.9999 199.9999 200 200 200 200 200 200 200 200 200
## [,35] [,36] [,37] [,38] [,39] [,40] [,41]
## [1,] 200 200 200 200 200 200 200
# • Decaying oscillations around n∗ = K.
# high rate to increase oscillation rate
# lower n0 than k
fun.ricker(n0=100, tmax=50, r=1.8, k=150)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## [1,] 100 182.2119 123.7951 169.5396 134.1034 162.2876 140.0387 157.8199 143.6839 154.9975 145.9755 153.1982 147.4301 152.0475 148.3572
## [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
## [1,] 151.3109 148.9493 150.8392 149.3278 150.5372 149.5699 150.3439 149.7248 150.2201 149.8239 150.1409 149.8873 150.0902 149.9279 150.0577
## [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45]
## [1,] 149.9538 150.0369 149.9705 150.0236 149.9811 150.0151 149.9879 150.0097 149.9923 150.0062 149.995 150.004 149.9968 150.0025 149.998
## [,46] [,47] [,48] [,49] [,50] [,51]
## [1,] 150.0016 149.9987 150.001 149.9992 150.0007 149.9995
# • Persistent, regular oscillations.
fun.ricker(n0=150, tmax=50, r=2.1, k=200)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## [1,] 150 253.5688 144.4831 258.8089 139.5742 263.2403 135.5098 266.7177 132.3771 269.2647 130.1146 271.0254 128.5665 272.1894 127.5501
## [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
## [1,] 272.9348 126.9023 273.4021 126.4974 273.6908 126.2477 273.8677 126.0948 273.9754 126.0018 274.0409 125.9454 274.0805 125.9112 274.1045
## [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45]
## [1,] 125.8905 274.119 125.878 274.1278 125.8704 274.1331 125.8658 274.1363 125.8631 274.1382 125.8614 274.1394 125.8604 274.1401 125.8598
## [,46] [,47] [,48] [,49] [,50] [,51]
## [1,] 274.1405 125.8594 274.1408 125.8592 274.1409 125.8591
# • Crazy, random-looking fluctuations (chaos).
fun.ricker(n0=150, tmax=50, r=4, k=100) # increase r to beyond 3
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 150 20.30029 492.0713 7.604177e-05 0.004151727 0.226639 12.2624 409.9526 0.001691629 0.09235358 5.023742 224.354 1.551254 79.59996
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28]
## [1,] 180.0105 7.334547 298.6325 0.1058124 5.752763 249.5276 0.6303148 33.5572 478.647 0.0001265444 0.006909055 0.3771174 20.28165 491.9861
## [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42]
## [1,] 7.628799e-05 0.00416517 0.2273727 12.30174 410.6211 0.00164968 0.09006352 4.899619 219.8999 1.816978 92.2494 125.7782 44.85311 407.1853
## [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50] [,51]
## [1,] 0.001876878 0.1024664 5.571592 243.427 0.7848497 41.52696 430.6365 0.0007769019 0.04241609
# Which parameter is the key driver of these patterns?
# r, rate of growth, seems to be the key driver of these patterns.
# (c) Choose six interesting values of this parameter. Write a script that makes an array of six plots showing the
# model dynamics for each of these values.
# set up plot windows in rows and columns
par(mfrow=c(2,3)) # 2 rows x 3 cols
# create a model, changing r each time
fun.growth.rate <- function(rr){
# input blank parameters
output <- matrix(NA, nrow=1, ncol=20+1) # set tmax = 20
output[1] <- 50 # set n0 = 50
# equation
for(t in 1:20){
output[t+1] <- output[t]*exp(rr*(1-output[t]/75)) # set k = 75
}
# plot
plot(0:20, output, xlab = 'Year', ylab = 'Population',
main = paste0('r = ', rr.vec[ii]), col = 'darkorchid4')
# add labels pf growth rate on each plot
lines(0:20, output) # adds lines connecting outputs
}
# loop through the model, changing parameter values each time
# rr = 1
# fun.growth.rate(rr) # try once
# create a vector to loop
rr.vec <- c(-0.5, -0.004, 0.5, 1, 1.8, 4)
# loop and plot each output
for(ii in 1:length(rr.vec)){
fun.growth.rate(rr.vec[ii])
}
# (d) Imagine n0 = 20 and K = 1000 for a certain population of deer that is growing according to the Ricker model.
# You are a wildlife manager, and are concerned about how long it will take for the population to reach half of its
# carrying capacity. That is, you want to know tK/2, the first year that nt ≥ K/2. Suppose your output time series is called nVec.
# Write an R command that will determine the index of the first element of nVec that is ≥ K/2.
fun.k500 <- function(r){
# input parameters to start with 20 and t = 1
nVec <- 20
t <- 1
while(nVec <= 500){
# start with time 1
# add a year
# determine what the population is
nVec <- nVec*exp(r*(1-nVec/1000))
t <- t + 1
}
return(t) # return the value where population exceeds 500
}
# run function, defining your own growth rate
fun.k500(r=.04) # t = 98
## [1] 98
fun.k500(r=1.2) # t = 4 for a faster growth rate
## [1] 4
# fun.k500(r=0) # continuous loop, 0 growth rate
# (e) Write a script that runs the necessary simulations and collects the necessary
# data to plot how tK/2 depends on r, for the range of r from 0.1 to 0.9.
# set parameters
r.vec <- seq(0.1:0.9, by=0.1)
n0 <- 20
kk <- 1000
n.vec <- rep(NA, 1, length(r.vec)) # repeat NA r.vec times
# write a for loop
for(ii in 1:length(r.vec)){
# create individual outputs for log results
log.output[ii] <- fun.k500(r.vec[ii])
# store each log result in position of nvec
n.vec[ii] <- log.output[ii]
}
# plot
par(mfrow=c(1,1))
plot(r.vec, n.vec, xlab = 'Instantaneous growth rate, r_d',
ylab = 'Time it takes to reach K/2 (n = 500)',
main = 'Sensitivity to value of growth rate',
type = 'b', col = 'darkcyan')
# (f) Write pseudo-code to do a joint sensitivity analysis for two parameters. That is, choose a vector of values you’d like
# to consider for both r and K, and choose a simple output from your model (e.g. population size at t = 10). Run the model for all possible combinations of these values and collect the results in a matrix with appropriate dimensions. Plot the results in some way. (Hint: the straight-forward way to do this uses a ‘nested loop’.)
# parameters - vary r and K
# goal: plot population size at t = 10 for various r and K
# use fun.ricker written earlier: fun.ricker(n0=100, tmax=10, r=1.05, k=200)
# choose 4 K values
k.vec <- c(20, 50, 75, 500)
# choose 4 r values
r.vec <- c(-1, 0.5, 1, 2.2)
# for each K value, loop fun.ricker over each r value
# set up for 4x4 output plots
par(mfrow=c(4,4))
# edit the equation
fun.ricker <- function(n0, tmax, r, k){
# input blank parameters
output <- matrix(NA, nrow=1, ncol=tmax+1)
output[1] <- n0
# equation
for(t in 1:tmax){
output[t+1] <- output[t]*exp(r*(1-output[t]/k))
}
# plot
plot(0:tmax, output, xlab = 'Year', ylab = 'Population',
main = paste0('Growth rate = ', r, ', K = ', k))
lines(0:tmax, output)
# return value
print(paste0('Growth rate = ', r, ' and K = ', k))
print(output)
}
# nested for loop to create plot
# n0 = 50
# tmax = 10
for(rr in r.vec){
for(kk in k.vec){
fun.ricker(n0 = 50, tmax = 10, r = rr, k = kk)
}
}
## [1] "Growth rate = -1 and K = 20"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 224.0845 6054096 Inf Inf Inf Inf Inf Inf Inf Inf
## [1] "Growth rate = -1 and K = 50"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 50 50 50 50 50 50 50 50 50 50
## [1] "Growth rate = -1 and K = 75"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 35.82657 21.25043 10.37829 4.384579 1.710105 0.6436218 0.2388159 0.08813565 0.03246142 0.01194706
## [1] "Growth rate = -1 and K = 500"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 20.32848 7.788747 2.910304 1.076891 0.3970202 0.1461716 0.05378924 0.01979009 0.007280654 0.002678442
## [1] "Growth rate = 0.5 and K = 20"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 23.61833 21.57564 20.74227 20.36091 20.17803 20.08842 20.04406 20.022 20.01099 20.00549
## [1] "Growth rate = 0.5 and K = 50"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 50 50 50 50 50 50 50 50 50 50
## [1] "Growth rate = 0.5 and K = 75"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 59.06802 65.68712 69.89462 72.31449 73.62083 74.30085 74.64798 74.82337 74.91153 74.95572
## [1] "Growth rate = 0.5 and K = 500"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 78.41561 119.5348 174.8752 242.0629 313.2917 377.6028 426.7678 459.1938 478.3193 488.8028
## [1] "Growth rate = 1 and K = 20"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 11.15651 17.3605 19.80972 19.99909 20 20 20 20 20 20
## [1] "Growth rate = 1 and K = 50"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 50 50 50 50 50 50 50 50 50 50
## [1] "Growth rate = 1 and K = 75"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 69.78062 74.80974 74.99976 75 75 75 75 75 75 75
## [1] "Growth rate = 1 and K = 500"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 122.9802 261.4028 421.2629 493.1095 499.9521 500 500 500 500 500
## [1] "Growth rate = 2.2 and K = 20"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 1.844158 13.58771 27.50912 12.04337 28.89717 10.85966 29.68056 10.23313 29.96369 10.01397
## [1] "Growth rate = 2.2 and K = 50"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 50 50 50 50 50 50 50 50 50 50
## [1] "Growth rate = 2.2 and K = 75"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 104.1005 44.33361 108.9937 40.211 111.5662 38.16848 112.4379 37.49555 112.6576 37.32752
## [1] "Growth rate = 2.2 and K = 500"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 50 362.1371 664.2241 322.4776 704.2492 286.6998 732.8613 263.0559 746.1455 252.6185 750.2144
# (g) (BONUS, OPTIONAL, PRETTY STRAIGHTFORWARD) Convert your pseudo-code into an R script to do a two-dimensional sensitivity
# analysis, and create a visually appealing plot to summarize the results. (You will need to google for plotting commands, e.g. contour or surface plots).
# (h) (BONUS, OPTIONAL, PRETTY HARD) Write a script to make a bifurcation plot for this model, showing how its long-term
# dynamics change as the parameter r is varied. You will need to collect a set of values reflecting the long-term dynamics of N for each value of r, where r falls between 0 and 4. Plot these N-values as points on the y-axis, versus the corresponding value of r on the x-axis. Hint: you may need to look up the R command matplot.
# knit into r markdown file for sharing
#library(knitr)
#spin('~/Documents/EEB Orientation/joannawu_lloyd-smith_homework.R')