# seascape models

## A ‘tidyverse’ approach to simulation models

I am a big fan of Rstudio’s (mainly led by Hadley Wickham) series of R packages known as the ‘tidyverse’. The packages modernise much of R’s functionality to better deal with moderately large databases and programming.

There is even a tidyverse package now, so you can hit up all their packages with one short `library(tidyverse)` call.

If you are into programming (ie doing more than just running a couple of linear models) then I recommend you check out the purrr package. It provides many functions that basically make looping easier and faster to code.

A stochastic simulation of population size generated using a loop (black line) and the accumulate function from purrr package (red points)

Here I explore how to use `purrr` to create a simple simulation model where a state variable at a time depends on its past state.

Skipping to my main point, it turns out that `purrr` may be more convenient in that it saves typing over writing loops. However, in this case the `purrr` function required more computational time.

The model we will use to test `purrr` will create stochastic simulations of population size that look something like is pictured.

First load the tidyverse package (or just purrr) and specify some parameters.

``````library(tidyverse)
#timesteps
tmax <- 100
#pop growth
r <- 1.2
#capacity
K <- 20
#abundance sd
sdev <- 0.1
``````

Now we generate a series of random numbers. We will multiply these by population size at each time-step to create stochastic variation:

``````set.seed(42) #so we get the same numbers
pmult <- exp(rnorm(100, sd = 0.1))
``````

Now, write a function that specifies how abundance changes at a single time point:

``````popmod <- function(N, popmult, r, K){
N * (r  - r*(N / K))*popmult
}
``````

The standard for loop was faster than accumulate

Now, let’s implement a function that uses a normal for loop to loop over abundances at time.

``````f <- function(pmult, tmax){
x2 <- rep(NA, tmax)
x2[1] <- 0.1
for (t in 2:tmax){
x2[t] <-
popmod(x2[t-1],pmult[t-1],  r, K)
}
x2
}
``````

We can use our new function to simulate abundance at time like this (as pictured):

``````loopabund <- f(pmult, tmax)
``````

Now for the `purrr` equivalent. We will use the accumulate function, which loops over a vector sequentially. Basically we would want to do this because it saves typing. Arguably it can also make your code more readable (assuming your reader knows how accumulate works).

``````accumulate(pmult, ~popmod(.x, .y, r = r, K = K), .init = 0.1)
``````

Accumulate works by taking a vector then iteratively applying a function. In this case we have applied it to our `popmod` function.

We have written our function as a formula, proceeded by `~`, for convenience. This means we can use the `.x` and `.y` arguments. The `.x` is the argument that will be accumulated (ie population size). The `.y` is the argument that we iteratively apply the function to - in this case the vector of multiples.

We can wrap our use of `accumulate` in a function, for convenience later.

``````f2 <- function(pmult) {
x3 <- accumulate(pmult, ~popmod(.x, .y, r = r, K = K), .init = 0.1)
x3
}
``````

Now we have two ways to simulate the same population. Let’s check they are the same by plotting them over each other:

``````accabund <- f2(pmult)[1:tmax]

#put our two simulations in a data.frame
datp <- data.frame(time = 1:tmax,
loopabund = loopabund,
accabund = accabund)

ggplot(datp, aes(x = time, y = loopabund)) +
geom_line() +
geom_point(aes(y = accabund), color = 'red') +
ylab("Abundance")
``````

So they are equivalent. We can then test for computational efficiency using the `microbenchmark` package:

``````microbenchmark::microbenchmark(f(pmult, tmax), f2(pmult)) %>% boxplot()
``````

It turns out that the for loop was about 40% faster than the purrr function aggregate. Such a time difference is trivial for this problem, but for more complex models could make a big difference. So the choice between purrr and base R comes down to whether you want to emphasise ease of coding and readability or speed of computation.

#### Contact: Chris Brown

Designed by Chris Brown. Source on Github