In IBAHCM/RPiR: Reproducible Programming in R

library(learnr)
tutorial_options(exercise.reveal_solution = FALSE)
gradethis::gradethis_setup()

knitr::opts_chunk$set(error = TRUE)
set.seed(123)

Overview

In this session we generalise the code from a growth model to a model with births and deaths.

Background

There are many refinements to our code that we could make, and indeed we will work through several of them during the course, but our current framework is already sufficiently sophisticated to handle some much more complex population dynamics. We are therefore ready to start to thinking about our models rather than the mechanics of what we are doing.

Tasks

In the modelling lecture we discussed a simple population growth model of the form

$$N(t+1) = N(t) + \lambda \times N(t)$$

which we can state roughly in English as

"the next population size is the current size plus the new additions (lambda times the current size)"

and this can be rewritten in turn in R as:

next.count <- latest$count + lambda * latest$count

This use of greek letters is fine for mathematicians, but we use growth.rate for lambda in our code as we want the code to read more like English. We will also write this on multiple lines as we calculate the new additions separately and then add them in. Note that, for consistency, throughout these practicals our code will refer to the data frame with the current latest population counts in as latest as in Practical 1-4, and these are referred to in the equations such as those on this page as $N(t)$. The updated values, referred to above as $N(t+1)$, appear as next.something in the functions if they have a name at all, and updated.something in the main script just to be distinct.

We now want to update the model to include both births and deaths, i.e.

$$N(t+1) = N(t) + b \times N(t) - d \times N(t)$$ which can be similarly translated as

"the next population size is the current size plus the new births (b times the current size) minus the new deaths (d times the current size)"

where $b$ represents the birth rate and $d$ represents the death rate.

In this practical, we will reuse the RPiR package while adapting 0104-run-growth.R and 0104-step-growth.R appropriately. For convenience, we've generated 0105-step-birth-death.R (from 0104-step-growth.R) for you, with blanks for you to fill in:

#' ### Function: step_deterministic_birth_death() 
#' Run one step of a simple deterministic exponential birth and death model. 
#'
#' Arguments: 
#' - latest -- the population count now
#' - birth.rate -- the birth rate 
#' - death.rate -- the death rate 
#'
#' Returns:
#' - the updated population count
#'
step_deterministic_birth_death <- function(latest, ___, ___) {
  # Calculate changes to population
  new.births <- ___ * latest$count 
  ___ <- ___ * latest$count 
  # Calculate population at next timestep
  next.count <- latest$count + new.births - ___
  # Create a data frame with the updated counts and return it
  data.frame(count = next.count) 
}

**Hint:** Note that we give the names of the arguments to the function in the comments at the top.

#' ### Function: step_deterministic_birth_death() 
#' Run one step of a simple deterministic exponential birth and death model. 
#'
#' Arguments: 
#' - latest -- the population count now
#' - birth.rate -- the birth rate 
#' - death.rate -- the death rate 
#'
#' Returns:
#' - the updated population count
#'
step_deterministic_birth_death <- function(latest, birth.rate, death.rate) {
  # Calculate changes to population
  new.births <- birth.rate * latest$count 
  new.deaths <- death.rate * latest$count 
  # Calculate population at next timestep
  next.count <- latest$count + new.births - new.deaths
  # Create a data frame with the updated counts and return it
  data.frame(count = next.count) 
}

grade_this_code()

Likewise, we've generated 0105-run-birth-death.R (from 0104-run-growth.R), with added blanks for you to fill in. Again,source() has been commented out as it won't work in the learnr environment. Instead, the necessary files have been preloaded for you. Try to fill in the blanks and compare the outputs as you vary the birth and death rates:

library(RPiR)
step_deterministic_birth_death <- function(latest, birth.rate, death.rate) {
  # Calculate changes to population
  new.births <- birth.rate * latest$count 
  new.deaths <- death.rate * latest$count 
  # Calculate population at next timestep
  next.count <- latest$count + new.births - new.deaths
  # Create a data frame with the updated counts and return it
  data.frame(count = next.count) 
}

library(RPiR)
# Load the step_deterministic_birth_death() function into the global environment (my workspace)
# Note this next line would be needed if you were running this in RStudio
# source("0105-step-birth-death.R")

#' Set up the simulation parameters
#' --------------------------------
#' First we set up the parameters for the simulation.

# Set the birth rate
human.annual.birth <- 0.015

# Set the death rate
___ <- ___

# Starting population size
initial.count <- 7000000000

# And setting times
start.time <- 0
end.time <- 100

#' Run the simplest possible simulation
#' ------------------------------------
#' Then run it so that we can get the output we need

# Set up the population starting size (at the first timestep)
population.df <- data.frame(count = initial.count)

# The timesteps that the simulation will run through
timesteps <- seq(from = start.time + 1, to = end.time)

# Now we loop through the time itself (starting at the second timestep)
for (new.time in timesteps) {
  # Calculate population at next timestep
  updated.population <- ___(latest = tail(population.df, 1),
                            ___ = human.annual.birth,
                            ___ = human.annual.death) 
  # Add new element onto end of population vector
  population.df <- rbind(population.df, updated.population)
}

#' Plot the results
#' ----------------
#' And finally we output the results.

# Now we can plot the timesteps against the population vector
population.df$time <- c(start.time, timesteps)
plot_populations(population.df)

#' Set up a second population
#' --------------------------
___

#' Run the simulation
#' ------------------
___

#' Plot the results on the same plot
#' ---------------------------------
___

You can plot a second (or third, etc.) result on top of a first with commands like:

plot_populations(population.b, new.graph = FALSE, col = "red")

However, if you want to do that, you should put the plot commands together, because the #' comment closes the current plot and you can't then draw on top of it. Remember that col sets the colour of the line in the plot -- use names "red", "blue", etc.

quiz(caption = "Quiz",
     question("When do you get the same answer in two different examples?",
              answer("When the death rate is minus the birth rate"),
              answer("When the death rate is zero"),
              answer("When the difference between the death rate and the birth rate is the same", 
                     correct = TRUE),
              answer("When the death rate plus the birth rate is zero"),
              allow_retry = TRUE,
     correct = "Correct! You should find that you get the same output if the difference between the birth and death rate is equal to the net growth rate in the previous model."),
     question("When do you get the same answer with the birth death model as the growth model?",
              answer("When the birth rate is zero"),
              answer("When the birth rate minus the death rate is equal to the growth rate", 
                     correct = TRUE),
              answer("When the death rate is double the birth rate"),
              answer("When the death rate is equal to the birth rate"),
              allow_retry = TRUE,
     correct = "Correct! You should find that you get the same output if the difference between the birth and death rate is equal to the net growth rate in the previous model.")
     )