Shows how to simulate ecological time series using the purrr accumulate functions and dplyr.

dplyr_simulations.R

# This example shows how to use purrr to simulate ecological values. The 
# accumulate function in purrr can be quite flexible as it can take lists as 
# entries, and pass extra parameters on to the function. The key feature
# to allow for time-varying parameters here is that the input for the function
# has to include both the present state of the population(s), but also a time
# index, so the function knows what values of the parameters to refer to. 

library(dplyr)
library(purrr)
library(ggplot2)

# Creating the population growth function. Note: the first two arguments
# have to be .x (the prior vector of population(s)) and .y, which
# would be the current state vector. There are extra parameters added
# after them. These have to be vectors, as the function indexes them. 
# This example function is a Ricker population growth model:
logistic_growth = function(.x, .y, growth, comp) {
  pop = .x[["pop"]] #gets the population parameter from the last time step (.x)
  time = .y[["time"]] #need to get the time of the current step (.y)
  #Note that the new_pop calculation uses the parameters indexed at the current time
  new_pop = pop*exp(growth[time] - pop*comp[time])
  return(c(time=time,pop=new_pop))
}


# Starting parameters: the number of time steps to simulate, intial population size,
# growth rate and intraspecific competition rate
n_steps  = 100
pop_init = 1
growth = 0.5
comp = 0.05
#First test: fixed growth rates
test1 = data_frame(time = 1:n_steps,pop_init = pop_init, 
                   growth=growth,comp =comp)


# note where sim_data is specified here: this takes the columns of your raw data
# and converts them into a list column where each item in the list is a named 
# list with time and initial population added. Then the accumulate function
# takes the arguments sim_data (the set of state vectors) and whatever other
# parameters you specified
out1 = test1 %>%
  mutate(
    sim_data = simplify_all(transpose(list(time=time,pop=pop_init))), 
    sims = accumulate(sim_data,logistic_growth, 
                      growth=growth,comp=comp))%>%
  bind_cols(simplify_all(transpose(.$sims)))



# This is the same example, except I drew the growth rates from a normal distribution
# with a mean equal to the mean growth rate and a std. dev. of 0.1
test2 = data_frame(time = 1:n_steps,pop_init = pop_init, 
                   growth=rnorm(n_steps, growth,0.1),comp=comp)

out2 = test2 %>%
  mutate(
    sim_data = simplify_all(transpose(list(time=time,pop=pop_init))),
    sims = accumulate(sim_data,logistic_growth, 
                      growth=growth,comp=comp))%>%
  bind_cols(simplify_all(transpose(.$sims)))


# This demostrates how to use this approach to simulate replicates using dplyr
# Note the crossing function creates all combinations of its input values
test3 = crossing(rep = 1:10, time = 1:n_steps,pop_init = pop_init, comp=comp) %>%
  mutate(growth=rnorm(n_steps*10, growth,0.1))

out3 = test3 %>%
  group_by(rep)%>%
  mutate(
    sim_data = simplify_all(transpose(list(time=time,pop=pop_init))),
    sims = accumulate(sim_data,logistic_growth, 
                      growth=growth,comp=comp))%>%
  bind_cols(simplify_all(transpose(.$sims)))

# Plots the output simulations
print(qplot(time, pop, data=out1)+
        geom_line() +
        geom_point(data= out2, col="red")+
        geom_line(data=out2, col="red")+
        geom_point(data=out3, col="red", alpha=0.1)+
        geom_line(data=out3, col="red", alpha=0.1,aes(group=rep)))