Exercise 1: length-weight relationship

Goal

Build a user-defined function

\[W = aL^b\]

  • \(a\) is very small (in the order of \(\times 10 ^{-3}\))
  • \(b\) ~ 3

Task

  • Write a function called LengthWeight()
  • Takes three input parameters: L, a, and b
  • returns W using this equation.
# This function will not work
LengthWeight <- function(L, a, b){
  W = a * (L ^ b)
}

Test the function

LengthWeight(L = 30, a = 0.0002, b = 3.05)

Doesn’t return anyhing. Not even an error or warning message… What happened?

Things to keep in mind

return statement

LengthWeight <- function(L, a, b){
  W = a * (L ^ b)
  return(W)
}

LengthWeight(L = 30, a = 0.0002, b = 3.05)
## [1] 6.401029
LengthWeight <- function(L, a, b){
  a * (L ^ b)
}

LengthWeight(L = 30, a = 0.0002, b = 3.05)
## [1] 6.401029

Define default values of parameters

LengthWeight <- function(L = 30, a = 0.0002, b = 3.05){
  W = a * (L ^ b)
  return(W)
}

LengthWeight()
## [1] 6.401029

Override fedault values of parameters

LengthWeight(L = 45, a = 0.001, b = 3.1)
## [1] 133.3395

Apply function to a vector of lengths

Imagine I sampled 10,000 fish from port, and was only able to take their lengths. Or that I did underwater surveys and recorded fish length. On any case, I don’t want to use the function 10,000 times. How do I use my LengthWeight function 10,000 times?

First let me just randomly simulate 1,000 fish lengths between 15 cm and 75 cm.

set.seed(43) # Set a random seed to obtain consistent random numbers
lengths <- runif(n = 10000, min = 15, max = 75) # Get 1000 lengths
hist(lengths)

for loop?

# WRONG!

ptm <- proc.time()
weights <- numeric(length = 10000)
for (i in 1:10000){
  weights[i] <- LengthWeight(L = lengths[i])
}

# for loops are VERY slow

proc.time() - ptm
##    user  system elapsed 
##   0.012   0.000   0.013

R is vectorized

Avoid for loops as much as you can! ALthough there are some cases where you HAVE to use them.

ptm <- proc.time()
weights <- LengthWeight(L = lengths)
proc.time() - ptm
##    user  system elapsed 
##   0.004   0.000   0.002

Exercise 2: von Bertalanffy growth equation

Goal

Build a more complex user-defined function

\[L = L_{\infty} \left(1 - e^{-k(t - t_0)}\right)\]

Task

  • Write a function VonBert -T akes four parameters: age (yr), k, Linfinity (cm), and t0 (yr)
  • Returns length (cm) from the von Bertalanffy equation
VonBert <- function(age, k = 0.2, Linfinity = 80, t0 = -0.2){
  length <- Linfinity * (1 - exp(-k * (age - t0)))
  return(length)
}

Test the function

VonBert(age = 5)
## [1] 51.72363

Things to keep in mind

You can use functions sequentially.

How much does teh fish above weight?

Approach 1

Save length-at-age into variable, then use that as input for length-weight function

length <- VonBert(age = 5)
LengthWeight(L = length)
## [1] 33.71175
# I could have saved the above into a variable, but nah

Approach 2

Nested functions

LengthWeight(L = VonBert(age = 5))
## [1] 33.71175

Approach 3

tidyverse approach using the pipe %>% (Ctrl / Cmd + Shift + M) from the magrittr package

library(magrittr)

VonBert(age = 5) %>% 
  LengthWeight()
## [1] 33.71175

Or even cooler

age <- 5

age %>% 
  VonBert() %>% 
  LengthWeight()
## [1] 33.71175

Vectorized and plotting directly

ages <- 1:60

ages %>% 
  VonBert() %>% 
  LengthWeight() %>% 
  plot(x = ages, type = "p", col = "steelblue", pch = 16)

Exercise 3: looping over age, printing age squared

Goal

Practice writing for loops

Task

  • Write a function AgeSquare()
  • Takes a vector ages as an input
  • Print out the square of the ages

Inside the function, you will:

  • Use the length() function to calculate how many elements there are in the ages vector
  • Store this in a variable nages
  • Then loop through the values in the ages vector using a for loop
  • Use print() to output the square of each length
AgeSquare <- function(ages){
  nages <- length(ages)
  for (i in 1:nages){
    print(ages[i]^2)
  }
}

Test it

AgeSquare(ages=1:10)
## [1] 1
## [1] 4
## [1] 9
## [1] 16
## [1] 25
## [1] 36
## [1] 49
## [1] 64
## [1] 81
## [1] 100

Things to keep in mind

This is actually what we would call a command. functions should return a value, while commands execute an operation. write.csv, print, cat… are commands.

Can I assign function output to a variable?

age2 <- AgeSquare(ages=1:10)
## [1] 1
## [1] 4
## [1] 9
## [1] 16
## [1] 25
## [1] 36
## [1] 49
## [1] 64
## [1] 81
## [1] 100
age2
## NULL

This function is not very useful, but it is meant to make you think about how to set up for loops for the few times where you will need to use them.

Exercise 4: looping over age, printing length

Goal

Call a user-defined function within a user-defined function (simila but not quite the same as nesting functions)

Task

  • Copy your AgeSquare() function, and rename it AgeLengths()
  • It will now take as parameters a vector ages, and parameters k, Linfinity, and t0 for the von Bertalanffy growth equation

Inside the for loop:

  • Instead of calculating and printing the square of the ages, call the VonBert() function you created in Exercise 2 to calculate and print the length at each age in the vector
AgeLengths <- function(ages, k = 0.2, Linfinity = 80, t0 = -0.2){
  nages <- length(ages)
  for (i in 1:nages){
    print(VonBert(ages[i], k = k, Linfinity = Linfinity, t0 = t0))
  }
}
AgeLengths(ages=1:10)
## [1] 17.06977
## [1] 28.47709
## [1] 37.81661
## [1] 45.46316
## [1] 51.72363
## [1] 56.84926
## [1] 61.04578
## [1] 64.4816
## [1] 67.29461
## [1] 69.5977

Things to keep in mind

You can do this same operation in other ways, depending on what you want. The above will work just fine, but it i) uses a for loop that we don’t need, and ii) doesn’t let you save the lengths to a variable. Let’s address both problems.

Avoid for loop by vectorizing

AgeLengths <- function(ages, k = 0.2, Linfinity = 80, t0 = -0.2){
    print(VonBert(ages, k = k, Linfinity = Linfinity, t0 - t0))
}

Test it

AgeLengths(ages=1:10)
##  [1] 14.50154 26.37440 36.09507 44.05368 50.56964 55.90446 60.27224
##  [8] 63.84828 66.77609 69.17318

Note that this function prints a long thing, rather than a tall thing. Also, notice the braketed number on the side. This printed a vector, instead of just a nummber.

Can I save to variable?

lengths <- AgeLengths(ages=1:10)
##  [1] 14.50154 26.37440 36.09507 44.05368 50.56964 55.90446 60.27224
##  [8] 63.84828 66.77609 69.17318
lengths
##  [1] 14.50154 26.37440 36.09507 44.05368 50.56964 55.90446 60.27224
##  [8] 63.84828 66.77609 69.17318

Yes! But it’s still printing, and I don’t like that.

Save output without printing

AgeLengths <- function(ages, k = 0.2, Linfinity = 80, t0 = -0.2){
    return(VonBert(ages, k = k, Linfinity = Linfinity, t0 - t0))
}
AgeLengths(ages=1:10)
##  [1] 14.50154 26.37440 36.09507 44.05368 50.56964 55.90446 60.27224
##  [8] 63.84828 66.77609 69.17318

Exercise 5: looping over age, storing and returning a vector of lengths

Goal

Practice creating vectors, storing values in them, and returning the vector

Task

  • Copy your AgeLengths() function, and rename it AgeLengthsVec()
  • Instead of using print() to report the calculated lengths, it will now return a vector lengths

Inside the function you will need to:

  • Create an empty vector lengths using the numeric() function
  • In the for loop store the calculated length in element i of the lengths vector, lengths[i]
  • At the end of the function, return the calculated vector using the return() command
AgeLengthsVec <- function(ages, k = 0.2, Linfinity = 80, t0 = -0.2){
  nages <- length(ages)
  lengths <- numeric(length = nages)
  for (i in 1:nages){
    lengths[i] <- VonBert(ages[i], k = k, Linfinity = Linfinity, t0 = t0)
  }
  return(lengths)
}
lengths <- AgeLengthsVec(ages=1:10)
lengths
##  [1] 17.06977 28.47709 37.81661 45.46316 51.72363 56.84926 61.04578
##  [8] 64.48160 67.29461 69.59770

Things to keep in mind

Again, this function can be made faster by vectorizing

AgeLengthsVec <- function(ages, k = 0.2, Linfinity = 80, t0 = -0.2){
  lengths<- VonBert(ages, k = k, Linfinity = Linfinity, t0 = t0)
  return(lengths)
}

If you think about it, this function is not doing anything special other than taking the exact same inputs and passing them to another function in the exact same way… but you get the jist.

No more loops from now on

Exercise 6: looping over age, length and weight

Goal

Learn how to compartmentalize your code

Task

Having built functions that convert age to length, and length to weight, now we can loop through vectors of age and predict the length and weight for fish at that age

  • Create a function AgeValues()
  • That takes as input a vector of ages age.vec, and the parameters in the above functions
  • And uses the functions you created in exercises 1 and 2 to return length at each age, and weight at each age.
AgeValues <- function(age.vec, k = 0.2, Linfinity = 80, t0 = -0.2, a = 0.0002, b = 3.05){
  lengths <- VonBert(age = age.vec, k = k, Linfinity = Linfinity, t0 = t0)
  weights <- LengthWeight(L = lengths, a = a, b = b)
  
  return(list(lengths, weights))
}
AgeValues(0:10)
## [[1]]
##  [1]  3.136845 17.069771 28.477086 37.816606 45.463158 51.723625 56.849263
##  [8] 61.045779 64.481597 67.294606 69.597703
## 
## [[2]]
##  [1]  0.006536327  1.146369107  5.460613369 12.970650946 22.745305722
##  [6] 33.711751572 44.971849015 55.882998976 66.040314381 75.226311999
## [11] 83.357458078

Things to keep in mind

You can have different classes of return statement, depending on what you want.

Unless you jump into the apply and purrr world, lists are difficult to work with. For example, we can’t just plot from the list:

AgeValues(0:10) %>% 
  plot()
## Error in xy.coords(x, y, xlabel, ylabel, log): 'x' is a list, but does not have components 'x' and 'y'

The above returns a list, because we told it to. You might want other things

Return a matrix

AgeValues <- function(age.vec, k = 0.2, Linfinity = 80, t0 = -0.2, a = 0.0002, b = 3.05){
  lengths <- VonBert(age = age.vec, k = k, Linfinity = Linfinity, t0 = t0)
  weights <- LengthWeight(L = lengths, a = a, b = b)
  
  return(matrix(c(lengths, weights), ncol = 2))
}
AgeValues(0:10)
##            [,1]         [,2]
##  [1,]  3.136845  0.006536327
##  [2,] 17.069771  1.146369107
##  [3,] 28.477086  5.460613369
##  [4,] 37.816606 12.970650946
##  [5,] 45.463158 22.745305722
##  [6,] 51.723625 33.711751572
##  [7,] 56.849263 44.971849015
##  [8,] 61.045779 55.882998976
##  [9,] 64.481597 66.040314381
## [10,] 67.294606 75.226311999
## [11,] 69.597703 83.357458078
AgeValues(0:10) %>% 
  plot(type = "p", col = "steelblue", pch = 16)

Return a data.frame

AgeValues <- function(age.vec, k = 0.2, Linfinity = 80, t0 = -0.2, a = 0.0002, b = 3.05){
  lengths <- VonBert(age = age.vec, k = k, Linfinity = Linfinity, t0 = t0)
  weights <- LengthWeight(L = lengths, a = a, b = b)
  
  return(data.frame(lengths, weights))
}
AgeValues(0:10)
##      lengths      weights
## 1   3.136845  0.006536327
## 2  17.069771  1.146369107
## 3  28.477086  5.460613369
## 4  37.816606 12.970650946
## 5  45.463158 22.745305722
## 6  51.723625 33.711751572
## 7  56.849263 44.971849015
## 8  61.045779 55.882998976
## 9  64.481597 66.040314381
## 10 67.294606 75.226311999
## 11 69.597703 83.357458078
AgeValues(0:10) %>% 
  plot(col = "steelblue", pch = 16)

Exercise 7: calculating Beverton-Holt steepness curves

The Beverton-Holt stock-recruitment curve can be parameterized with a steepness parameter \(h\), together with \(R_0\) (unfished recruitment) and \(S_0\) (unfished spawning biomass) using the following set of equations:

\[ S_0 = SPR_0 \times R_0 \]

\[ \alpha = \frac{1 - h}{4hR_0}S_0 \]

\[ \beta = \frac{5h - 1}{4hR_0} \]

\[ R_{t+1} = \frac{S_t}{\alpha + \beta S_t} \] ## Goal

Run the same function many times to compare parameter values

Task

  • Assume that the spawner-per-recruit in the absence of fishing (\(SPR_0\)) has been calculated already
  • The parameters of the Beverton-Holt equation will be \(h\), \(R_0\) and \(S_0\)
  • Create a function BevHolt()
  • Takes as an input these three parameters, and loops over (remember, not looping any more)
  • 100 equally-spaced values of spawning biomass \(S_t\) from 0 to \(S_0\), using the length.out option of seq()
  • For each value of spawning biomass, calculates and stores the corresponding recruitment \(R_{t+1}\)
  • Store the values of \(S_t\) and \(R_{t+1}\) in vectors Svec and Rvec
  • Create a plot to illustrate the relationship between the two
  • Vary h to examine how the curve looks at different steepness values

Note: if you set \(S_0 = 1\) and \(R_0 = 1\) then the plot will be relative to unfished levels

BevHolt <- function(h = 0.5, R0 = 1, S0 = 1){
  alpha <- ((1 - h) / (4 * h * R0)) * S0
  beta <- ((5 * h) - 1) / (4 * h * R0)
  
  Svec <- seq(from = 0, to = S0, length.out = 100)
  Rvec <- Svec / (alpha + (beta * Svec))

  plot(x = Svec, y = Rvec, type = "l", col = "steelblue")
}
BevHolt()

Things to keep in mind

Functions should only do one thing (i.e. Generate data, transform data, calculate a variable, plot something…)

Why? well, I can’t get the sample plot by just using the function above. i might want to generate lots of data for different values of \(h\) and then look at it together.

BevHolt returns values as data.frame

BevHolt <- function(h = 0.5, R0 = 1, S0 = 1){
  alpha <- ((1 - h) / (4 * h * R0)) * S0
  beta <- ((5 * h) - 1) / (4 * h * R0)
  
  Svec <- seq(from = 0, to = S0, length.out = 100)
  Rvec <- Svec / (alpha + (beta * Svec))

  results <- data.frame(h = h, Svec = Svec, Rvec = Rvec)
  
  return(results)
}
BevHolt()
##       h       Svec       Rvec
## 1   0.5 0.00000000 0.00000000
## 2   0.5 0.01010101 0.03921569
## 3   0.5 0.02020202 0.07619048
## 4   0.5 0.03030303 0.11111111
## 5   0.5 0.04040404 0.14414414
## 6   0.5 0.05050505 0.17543860
## 7   0.5 0.06060606 0.20512821
## 8   0.5 0.07070707 0.23333333
## 9   0.5 0.08080808 0.26016260
## 10  0.5 0.09090909 0.28571429
## 11  0.5 0.10101010 0.31007752
## 12  0.5 0.11111111 0.33333333
## 13  0.5 0.12121212 0.35555556
## 14  0.5 0.13131313 0.37681159
## 15  0.5 0.14141414 0.39716312
## 16  0.5 0.15151515 0.41666667
## 17  0.5 0.16161616 0.43537415
## 18  0.5 0.17171717 0.45333333
## 19  0.5 0.18181818 0.47058824
## 20  0.5 0.19191919 0.48717949
## 21  0.5 0.20202020 0.50314465
## 22  0.5 0.21212121 0.51851852
## 23  0.5 0.22222222 0.53333333
## 24  0.5 0.23232323 0.54761905
## 25  0.5 0.24242424 0.56140351
## 26  0.5 0.25252525 0.57471264
## 27  0.5 0.26262626 0.58757062
## 28  0.5 0.27272727 0.60000000
## 29  0.5 0.28282828 0.61202186
## 30  0.5 0.29292929 0.62365591
## 31  0.5 0.30303030 0.63492063
## 32  0.5 0.31313131 0.64583333
## 33  0.5 0.32323232 0.65641026
## 34  0.5 0.33333333 0.66666667
## 35  0.5 0.34343434 0.67661692
## 36  0.5 0.35353535 0.68627451
## 37  0.5 0.36363636 0.69565217
## 38  0.5 0.37373737 0.70476190
## 39  0.5 0.38383838 0.71361502
## 40  0.5 0.39393939 0.72222222
## 41  0.5 0.40404040 0.73059361
## 42  0.5 0.41414141 0.73873874
## 43  0.5 0.42424242 0.74666667
## 44  0.5 0.43434343 0.75438596
## 45  0.5 0.44444444 0.76190476
## 46  0.5 0.45454545 0.76923077
## 47  0.5 0.46464646 0.77637131
## 48  0.5 0.47474747 0.78333333
## 49  0.5 0.48484848 0.79012346
## 50  0.5 0.49494949 0.79674797
## 51  0.5 0.50505051 0.80321285
## 52  0.5 0.51515152 0.80952381
## 53  0.5 0.52525253 0.81568627
## 54  0.5 0.53535354 0.82170543
## 55  0.5 0.54545455 0.82758621
## 56  0.5 0.55555556 0.83333333
## 57  0.5 0.56565657 0.83895131
## 58  0.5 0.57575758 0.84444444
## 59  0.5 0.58585859 0.84981685
## 60  0.5 0.59595960 0.85507246
## 61  0.5 0.60606061 0.86021505
## 62  0.5 0.61616162 0.86524823
## 63  0.5 0.62626263 0.87017544
## 64  0.5 0.63636364 0.87500000
## 65  0.5 0.64646465 0.87972509
## 66  0.5 0.65656566 0.88435374
## 67  0.5 0.66666667 0.88888889
## 68  0.5 0.67676768 0.89333333
## 69  0.5 0.68686869 0.89768977
## 70  0.5 0.69696970 0.90196078
## 71  0.5 0.70707071 0.90614887
## 72  0.5 0.71717172 0.91025641
## 73  0.5 0.72727273 0.91428571
## 74  0.5 0.73737374 0.91823899
## 75  0.5 0.74747475 0.92211838
## 76  0.5 0.75757576 0.92592593
## 77  0.5 0.76767677 0.92966361
## 78  0.5 0.77777778 0.93333333
## 79  0.5 0.78787879 0.93693694
## 80  0.5 0.79797980 0.94047619
## 81  0.5 0.80808081 0.94395280
## 82  0.5 0.81818182 0.94736842
## 83  0.5 0.82828283 0.95072464
## 84  0.5 0.83838384 0.95402299
## 85  0.5 0.84848485 0.95726496
## 86  0.5 0.85858586 0.96045198
## 87  0.5 0.86868687 0.96358543
## 88  0.5 0.87878788 0.96666667
## 89  0.5 0.88888889 0.96969697
## 90  0.5 0.89898990 0.97267760
## 91  0.5 0.90909091 0.97560976
## 92  0.5 0.91919192 0.97849462
## 93  0.5 0.92929293 0.98133333
## 94  0.5 0.93939394 0.98412698
## 95  0.5 0.94949495 0.98687664
## 96  0.5 0.95959596 0.98958333
## 97  0.5 0.96969697 0.99224806
## 98  0.5 0.97979798 0.99487179
## 99  0.5 0.98989899 0.99745547
## 100 0.5 1.00000000 1.00000000

Test different values of \(h\) and plot with ggplot2

library(ggplot2)
## Error in library(ggplot2): there is no package called 'ggplot2'
multiple_h <- BevHolt(h = 0.2) %>% 
  rbind(BevHolt(h = 0.3)) %>%
  rbind(BevHolt(h = 0.4)) %>%
  rbind(BevHolt(h = 0.5)) %>%
  rbind(BevHolt(h = 0.6)) %>%
  rbind(BevHolt(h = 0.7)) %>%
  rbind(BevHolt(h = 0.8)) %>%
  rbind(BevHolt(h = 0.9))


ggplot(data = multiple_h, mapping = aes(x = Svec, y = Rvec, color = h, group = h)) +
  geom_line() +
  theme_bw() +
  scale_color_viridis_c()
## Error in ggplot(data = multiple_h, mapping = aes(x = Svec, y = Rvec, color = h, : could not find function "ggplot"

A sort of “sensitivity analysis”

ggplot(data = multiple_h, mapping = aes(x = Svec, y = h, fill = Rvec, z = Rvec)) +
  geom_raster(interpolate = T) +
  geom_contour(color = "black", linetype = "dashed") +
  theme_bw() +
  scale_fill_viridis_c() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  geom_point(x = 0.4, y = 0.6, size = 3, fill = "red", color = "black", shape = 21)
## Error in ggplot(data = multiple_h, mapping = aes(x = Svec, y = h, fill = Rvec, : could not find function "ggplot"

Plotting inside functions is not adviced, unless it makes a lot of sense. But still, make a “data” function and then a “plotting” function

PlotBevHolt <- function(data){
  data %>% 
    ggplot(aes(x = Svec, y = Rvec)) +
    geom_line(linetype = "dashed", size = 1, color = "steelblue") +
    theme_bw()
}
BevHolt() %>% 
  PlotBevHolt()
## Error in ggplot(., aes(x = Svec, y = Rvec)): could not find function "ggplot"

Exercise 8: logistic selectivity

In many fisheries, selectivity of the fishing gear increases with length, from zero at small lengths to one at high lengths. It is often convenient to use a logistic curve with parameters \(L_50\) (length at 50% selectivity) and \(L_95\) (length at 95% selectivity) as popularized by Punt and Kennedy (1997):

\[ S = \frac{1}{1 + e^{\left(-ln(19)\frac{L - L_{50}}{L_{95} - L_{50}} \right)}} \]

Task

  • Write a function LogisticSel()
  • Takes three parameters: a vector of lengths Lvec (cm), L50 (cm), and L95 (cm)
  • Returns the selectivity of the gear at that length (ranging from 0 to 1)
LogisticSel <- function(Lvec = 55, L50 = 40, L95 = 60){
  S <- 1 / (1 + exp(-1* log(19) * ((Lvec - L50) / (L95 - L50))))
  return(S)
}
LogisticSel()
## [1] 0.900995
LogisticSel(Lvec = 0:100) %>% 
  plot(col = "steelblue", pch = 16)