# 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]
## [1,] 100 36.78794 7.192578 1.046002 0.1430495 0.01938736 0.002624302 0.00035517 4.806721e-05 6.505192e-06
## [,11]
## [1,] 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]
## [1,] 50 72.74957 99.99789 128.4005 153.5695 172.471 184.7589 191.9346 195.8439 197.8894 198.9363 199.466
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 199.7325 199.8661 199.933 199.9665 199.9832 199.9916 199.9958 199.9979 199.999 199.9995 199.9997 199.9999
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
## [1,] 199.9999 200 200 200 200 200 200 200 200 200 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]
## [1,] 100 182.2119 123.7951 169.5396 134.1034 162.2876 140.0387 157.8199 143.6839 154.9975 145.9755 153.1982
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 147.4301 152.0475 148.3572 151.3109 148.9493 150.8392 149.3278 150.5372 149.5699 150.3439 149.7248 150.2201
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## [1,] 149.8239 150.1409 149.8873 150.0902 149.9279 150.0577 149.9538 150.0369 149.9705 150.0236 149.9811 150.0151
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## [1,] 149.9879 150.0097 149.9923 150.0062 149.995 150.004 149.9968 150.0025 149.998 150.0016 149.9987 150.001
## [,49] [,50] [,51]
## [1,] 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]
## [1,] 150 253.5688 144.4831 258.8089 139.5742 263.2403 135.5098 266.7177 132.3771 269.2647 130.1146 271.0254
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 128.5665 272.1894 127.5501 272.9348 126.9023 273.4021 126.4974 273.6908 126.2477 273.8677 126.0948 273.9754
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## [1,] 126.0018 274.0409 125.9454 274.0805 125.9112 274.1045 125.8905 274.119 125.878 274.1278 125.8704 274.1331
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## [1,] 125.8658 274.1363 125.8631 274.1382 125.8614 274.1394 125.8604 274.1401 125.8598 274.1405 125.8594 274.1408
## [,49] [,50] [,51]
## [1,] 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]
## [1,] 150 20.30029 492.0713 7.604177e-05 0.004151727 0.226639 12.2624 409.9526 0.001691629 0.09235358 5.023742
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## [1,] 224.354 1.551254 79.59996 180.0105 7.334547 298.6325 0.1058124 5.752763 249.5276 0.6303148 33.5572 478.647
## [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## [1,] 0.0001265444 0.006909055 0.3771174 20.28165 491.9861 7.628799e-05 0.00416517 0.2273727 12.30174 410.6211
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [1,] 0.00164968 0.09006352 4.899619 219.8999 1.816978 92.2494 125.7782 44.85311 407.1853 0.001876878 0.1024664
## [,45] [,46] [,47] [,48] [,49] [,50] [,51]
## [1,] 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
# store each log result in position of nvec
n.vec[ii] <- fun.k500(r.vec[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))
tmax = 10
# input blank parameters
output <- matrix(NA, nrow=1, ncol=tmax+1)
output[1] <- n0
# edit the equation
fun.ricker <- function(n0, tmax, rr, kk){
# equation
for(t in 1:tmax){
output[t+1] <- output[t]*exp(rr*(1-output[t]/kk))
}
# plot
plot(0:tmax, output, xlab = 'Year', ylab = 'Population',
main = paste0('Growth rate = ', rr, ', K = ', kk))
lines(0:tmax, output)
# return value
print(paste0('Growth rate = ', rr, ' and K = ', kk))
print(output)
}
# initialize matrix
out.mat <- matrix(rep(NA), nrow = length(r.vec), ncol = length(k.vec))
# use a nested for loop to create plot
# n0 = 50
# tmax = 10
for(rr in 1:length(r.vec)){
for(kk in 1:length(k.vec)){
# store output in a vector
output <- fun.ricker(n0 = 50, tmax = 10, rr = r.vec[rr], kk = k.vec[kk])
pop.t10 <- output[1,11]
# store each result in position of the matrix
out.mat[rr,kk] <- pop.t10
}
}
## [1] "Growth rate = -1 and K = 20"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 20 20 20 20 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,] 20 10.97623 5.02918 2.045906 0.7840822 0.2930067 0.1084247 0.0399738 0.0147173 0.005415786 0.001992572
## [1] "Growth rate = -1 and K = 75"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 9.606106 4.016779 1.558989 0.5855662 0.2171062 0.08010045 0.0294988 0.01085627 0.003994377 0.001469527
## [1] "Growth rate = -1 and K = 500"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 7.657858 2.860647 1.058412 0.390193 0.143656 0.05286329 0.01944937 0.007155303 0.002632326 0.0009683839
## [1] "Growth rate = 0.5 and K = 20"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 20 20 20 20 20 20 20 20 20 20
## [1] "Growth rate = 0.5 and K = 50"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 26.99718 33.97961 39.88357 44.12951 46.79768 48.32055 49.13892 49.56387 49.78051 49.88989
## [1] "Growth rate = 0.5 and K = 75"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 28.85834 39.2523 49.81559 58.92256 65.58893 69.83585 72.282 73.60369 74.29204 74.64351
## [1] "Growth rate = 0.5 and K = 500"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 32.32149 51.59428 80.78704 122.858 179.1407 246.9117 318.0215 381.4949 429.4918 460.8676
## [1] "Growth rate = 1 and K = 20"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 20 20 20 20 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,] 20 36.44238 47.79327 49.94985 49.99997 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,] 20 41.64018 64.96594 74.2658 74.99638 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,] 20 52.23393 127.902 269.2016 427.1167 494.1425 499.9654 500 500 500 500
## [1] "Growth rate = 2.2 and K = 20"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 20 20 20 20 20 20 20 20 20 20
## [1] "Growth rate = 2.2 and K = 50"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 74.86843 25.06623 75.08399 24.90109 75.13328 24.86346 75.14406 24.85524 75.14639 24.85346
## [1] "Growth rate = 2.2 and K = 75"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 100.3903 47.6691 106.2707 42.46649 110.281 39.17828 112.0441 37.79832 112.5631 37.39971
## [1] "Growth rate = 2.2 and K = 500"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 20 165.2951 720.852 272.7853 741.318 256.3722 748.8904 250.504 750.8886 248.9738 751.3435
# output matrix of 4x4 r and k values
out.mat
## [,1] [,2] [,3] [,4]
## [1,] 20 0.001992572 0.001469527 9.683839e-04
## [2,] 20 49.889891282 74.643509869 4.608676e+02
## [3,] 20 50.000000000 75.000000000 5.000000e+02
## [4,] 20 24.853460740 37.399713942 7.513435e+02
# (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')