Writing functions in R

library(tidyverse)

Today we will be learning to write our own functions in R. Consider the discrete Logistic growth function:

\[x_{t+1} = x_t + r x_t \left(1 - \frac{x_t}{K} \right)\]

Here we exress this function in R:

growth <- function(xt, r, K){
  r * xt * (1 - xt / K)
}

We can simulate the function as follows:

simulate <- function(max_time, x0, r, K){

  x <- numeric(max_time)
  x[1] <- x0
  
  for (t in 1:(max_time - 1)) {
    x[t + 1] <- x[t] + growth(x[t], r = r, K)
  }
  
  pop_data <- data_frame(time = 1:max_time, pop_size = x)
  pop_data
}

We can perform a few simulations with different values of r:

simulate(50, 50, 1, 100) %>%
  ggplot(aes(time, pop_size)) + geom_line()

simulate(50, 50, 2.2, 100) %>%
  ggplot(aes(time, pop_size)) + geom_line()

simulate(50, 50, 2.7, 100) %>%
  ggplot(aes(time, pop_size)) + geom_line()

Wow, things get pretty crazy as r gets large. Let’s try and get a more detailed picture of the long-term behavior of the system as a function of the r value:

## Try all these values of r
r_range <- seq(1,3, length.out = 400)


## An initial simulation. 
df <- simulate(300, 10, r_range[1], 100) %>% 
    filter(time > 200)  %>% # use only the long-term behavior
    mutate(r = r_range[1])  # store the value of r that we used.

# Okay, now repeat this for all values of r: 
for (r in r_range) {
  
  df1 <- simulate(300, 10, r, 100) %>% 
    filter(time > 200) %>% 
    mutate(r = r)
  
  ## technical note: this could be slow for really big data, technically we should "initialize" this first.
  ## Instead, we'll later learn how to automate that with functional programming. 
  df <- bind_rows(df, df1)
  
}

Let’s plot the results showing the equilibrium population size(s) as a function of the value of r we chose:

df %>% 
  ggplot(aes(r, pop_size)) + 
  geom_point(size = 0.1)