Open the leslie_projection.Rmd file in docs/. All code shown below for Part 0 is in that file. You should be able to run it out-of-the-box. In EDS 212, we learned about Leslie matrices for projecting populations based on mortality and reproduction for different organism life stages. Let’s return to one example we checked out.
A population of fish we’re interested in has four life stages: eggs (E), fry (F), juvenile (J), breeding adult (A). You are told that annually:
# Make the Leslie Matrix:
fish_leslie <- matrix(c(0, 0, 0, 200, 0.05, 0, 0, 0, 0, 0.10, 0, 0, 0, 0, 0.16, 0.9), nrow = 4, ncol = 4, byrow = TRUE)
# Check it out:
fish_leslie [,1] [,2] [,3] [,4]
[1,] 0.00 0.0 0.00 200.0
[2,] 0.05 0.0 0.00 0.0
[3,] 0.00 0.1 0.00 0.0
[4,] 0.00 0.0 0.16 0.9Recall, we are able to use the dot product (%*%) to estimate the population at the next year (we’ll draw this on the board to remember what it looks like).
The initial population structure is given by:
# Initial population structure (Roe, Fry, Juvenile, Adult):
fish_year0 <- c(5000, 8000, 600, 400)
# ------ WHAT HAPPENS HERE? ------ #
time_yr <- seq(from = 0, to = 8, by = 1)
proj_year <- matrix(ncol = length(fish_year0), nrow = length(time_yr))
proj_year[1, ] <- fish_year0
for (i in 2:length(time_yr)) {
proj_year[i,] <- fish_leslie %*% proj_year[i-1,]
}
# -------------------------------- #
# The rest is wrangling & visualization:
colnames(proj_year) <- c("eggs", "fry", "juvenile", "adult")
proj_df <- data.frame(time_yr, proj_year) %>%
pivot_longer(cols = -time_yr, names_to = "lifestage", values_to = "stage_population")
ggplot(data = proj_df, aes(x = time_yr, y = stage_population)) +
geom_line(aes(color = lifestage)) +
scale_y_log10()
{testthat}fun_testing.Rmdtestthat package (install.packages("testthat"))devtools package (install.packages("devtools"))remotes packages (install.packages("remotes"))As we develop algorithms, we’ll change our code. We want a way to automatically check our work to make sure it’s behaving as expected. The testthat package “tries to make testing as fun as possible.”
Unit test: A unit test is an automated check of a piece (“unit”) of your code
Let’s consider an example. We’re writing a function to find the mean value of each column, then return the lowest and highest mean calculated (in that order). We would expect:
Let’s write the function first, then create some accompanying tests to make sure they’re working.
mean_range <- function(df) {
col_means <- apply(X = df, MARGIN = 2, FUN = mean, na.rm = TRUE) # Returns column means as a vector
col_mean_max <- max(col_means) # Looks for the maximum value in the vector
col_mean_min <- min(col_means) # Looks for the minimum value in the vector
return(c(col_mean_min, col_mean_max)) # Prints the vector with minimum & maximum
}
# Try it out:
mean_range(df = mtcars)[1] 0.40625 230.72188OK great. But we’re doing some work on this code, and we don’t want to have to try a bunch of different things manually each time we change it to see how it breaks. Instead, let’s write some automated tests for this function that help us avoid that.
testthatThe testthat package contains a whole bunch of helper functions built for testing. Take a look at options in the package documentation.
Here are some examples of the types of built-in functions for testing:
expect_length(): does the function return a vector of the expected length?expect_equal(): does the function return an expected value?expect_true(): does the code return TRUE or FALSE?For example, let’s use expect_length() to test our function above, which should return a vector of length 2.
expect_length(mean_range(mtcars), 2)Nothing is returned - that’s passing. You’ll get an error message if the test fails. What if we say the expected length of the vector returned is 3?
expect_length(mean_range(mtcars), 3)Error: mean_range(mtcars) has length 2, not length 3.Let’s do another one!
For this function to be working correctly, the output of mean_range(mtcars) should have two values and we expect the first element to be smaller than the second element.
Our test might look something like:
expect_true(mean_range(mtcars)[1] < mean_range(mtcars)[2])That passes! What happens if you change the alligator mouth direction?
expect_true(mean_range(mtcars)[1] > mean_range(mtcars)[2])Error: mean_range(mtcars)[1] > mean_range(mtcars)[2] is not TRUE
`actual`: FALSE
`expected`: TRUE That test fails. This is just a little flavor of what unit tests do.
When you’re writing functions (or working with someone else’s functions), often they’re not stored in the actual script or notebook you’re working in. If that’s the case, you can still use them by sourcing a function.
Let’s make a function that’s stored in a script, save it in our project working directory, and use it in an R Markdown document.
In our new R script, let’s write a function that, when given inputs of the user’s favorite food and animal, returns the name of their future food cart.
name_cart <- function(food, animal) {
print(paste("Mc", stringr::str_to_title(animal), "'s ", stringr::str_to_title(food), "Mart", sep = ""))
}
# Try it out:
name_cart(food = "burrito", animal = "pika")Save the script with a descriptive name, e.g. food_cart_functions.R (this could, in theory, contain multiple functions). Atop our .Rmd (in the setup chunk), we’ll source the script so we have access to its functions. Mine will look like this:
source("food_cart_functions.R")
You may have your scripts in a subfolder (e.g. R/ or something) - you can use here::here() within the source() function to point to the right place within your R Project as needed.
Now, add a code chunk in your R Markdown document that uses the name_cart() function you’ve created. It works over here, too!
name_cart(food = "strawberry", animal = "marmot")[1] "McMarmot's StrawberryMart"This gives us a valuable tool for organizing things within a project. Often, you’ll have a separate folder called R where you have your scripts / functions, and those are sourced as needed for use in your notebook or by another script you’re working on.
Mariola Sánchez-González, Margarida Tomé, Gregorio Montero. Modelling height and diameter growth of dominant cork oak trees in Spain. Annals of Forest Science, Springer Nature (since 2011)/EDP Science (until 2010), 2005, 62 (7), pp.633-643.
In the src/ folder there is a script cork_oak_growth.R that contains one function height_t2 that estimates the height of a cork oak tree at some age (time 2), given a prior height of the tree, the age when that height was recorded (years), and the age you’re trying to estimate the new height for (these ages are t1 and t2, respectively, in the model).
Create a new .Rmd. Source the script by running source(here::here("path", "to", "script.R")) in your setup chunk. Ensure that you can use the function in your .Rmd whenever you want to.
Documentation is a really important part of data science. Remember: “past you doesn’t respond to emails” (Wickham). I’ll also add that present you will forget things immediately, and future you is a grump. So make sure, for your sake and the sake of your collaborators, you document your functions (along with any useful annotation throughout your work) to make them easier for you to use later on.
It’s important to note here that documentation can mean different things.
Informally, this can mean adding relevant comments to your script so that you know:
Formally, in R, the {roxygen2} package can help guide our documentation.
Let’s try creating some documentation for a little function. Write a function to calculate the maximum sustained harvest (\(H\)), given a fishery carrying capacity (\(K\)) and intrinsic growth rate, \(r\):
\[H = \frac{Kr}{4}\]
# Write the functio
max_sustainable_harvest <- function (K, r) {
harvest <- (K * r) / 4
return(harvest)
}
# Try it out:
max_sustainable_harvest(K = 36000, r = 0.31)[1] 2790Now, how can we document this function?
We’ll use what’s called a Roxygen Skeleton - a skeleton set of documentation that prompts us to input important information about our function. Note: this exercise is just so you get a feel for what a roxygen skeleton looks like.
To add a Roxygen Skeleton:
#' Title
#'
#' @param
#' @param
#'
#' @return
#' @export
#'
#' @examplesUpdate the fields with descriptions and types
@param: information about the parameters (arguments) in your function@return: information about what the function returns to the user@export: we’ll return to this, but this lets an R package know that this is a function available for the user to use (not an internal function)@examples: add helpful examples for the user#' Maximum sustainable harvest
#'
#' A little function to return the maximum sustainable harvest given the carrying capacity and growth rate.
#'
#' @param K A number indicating the carrying capacity of the fishery (in individuals)
#' @param r A number indicating the intrinsic growth rate of the fishery (in 1 / yr)
#'
#' @return A number indicating the maximum annual sustainable harvest for the fishery
#' @export
#'
#' @examples
#' max_sustainable_yield(K = 3.8e5, r = 0.62)So overall, our function code chunk might look like this:
#' Maximum sustainable harvest
#'
#' A little function to return the maximum sustainable harvest given the carrying capacity and growth rate.
#'
#' @param K A number indicating the carrying capacity of the fishery (in individuals)
#' @param r A number indicating the intrinsic growth rate of the fishery (in 1 / yr)
#'
#' @return A number indicating the calculated maximum annual sustainable harvest for the fishery
#' @export
#'
#' @examples
#' max_sustainable_harvest(K = 3.8e5, r = 0.62)
max_sustainable_harvest <- function (K, r) {
harvest <- (K * r) / 4
return(harvest)
}
# Try it out:
max_sustainable_harvest(K = 36000, r = 0.31)[1] 2790Now, this might seem weird because it’s not obvious how that odd formatting of the Roxygen Skeleton actually helps. Why not just make them normal comments, and it words that a human can easily read?
The answer is: because this structure and syntax will make it possible for us to store our documentation as part of an R package containing our function. And that’s what we want to do next…
Let’s put a bunch of what we’ve been working on together into an R package that you can use, and share with others so they can use your amazing function!
Step-by-step:
In RStudio, new R Project > New directory > R package
Name your R package something awesome, like your favorite food and favorite animal (e.g. tacopika)
Save your new package project in your EDS 221 folder, and open your project in a new session
Hey, you made an R package!
hello() that prints the worlds “Hello, world!”R folderWait you made an R package? Yup try it out.
library(your_package_name) in the Console)hello(). Try it!OK but let’s really make it your own…
Create a new R script in your project
Create a function within the script of your choosing (you’ll share this with a classmate). Some requirements:
Save the script in the R subfolder with the same name as your function. For example, if your function is called jurassic_park(), then the script should be saved in the R folder as jurassic_park.R.
Click ‘Install and restart’ in the Build tab again
Try your new function!
So that’s cool that your new function exists within the package and you can use it. But now try checking the documentation on your function using ?function_name. What shows up? Zilch, right?
We want to add documentation for our function that we can check it out whenever we use the function, and so that other users will have critical information about what it does and how it works.
devtools::document() (note: if you see ‘Warning: Namespace not generated by Roxygen2’, you should delete the existing NAMESPACE file, then run devtools::document() again to resolve)man/ folder and see that an .Rd (R documentation) file has been added that is associated with your function name, and that contains your updated documentationThings to remember:
devtools::document() any time you update your function or documentationInstall and restart to see the changes when you testPutting it on GitHub
Your friends and family will definitely want to use the amazing R package you’ve just created. How can they install and use it?
Collaborators can install your package straight from a GitHub repo. Let’s put it up there.
usethis::use_git()usethis::use_github()This should take you to your new GitHub repo for the package.
Installing a package from GitHub
Once your package is in a public repo on GitHub, anyone can install its contents by running:
`remotes::install_github("username/reponame")`For example, if my package is in a repo called tacopika, then someone could install that package of mine using:
`remotes::install_github("ryoliver/tacopika")`You can check all functions in a package by running lsf.str("package:packagename") in the console. For example, after installing the package tacopika above, running lsf.str("package:tacopika") in the Console will return all functions in it.
Now that you’ve made your own package, let’s dive into packages that other developers have made so that we do more data science and less reinventing the wheel with code to do common useful things.
Let’s spend a bit of time learning how to search for, choose, investigate put to use some R & Python Packages.
CRAN Task Views (search by topic): https://cran.r-project.org/web/views/
Microsoft R Application Network (MRAN) (search by keywords): https://mran.microsoft.com/packages
Some trusted & awesome package sources:
Other great places to learn about package updates, releases, use:
Explore packages:
Dig a bit deeper:
Installation
We’ve seen the two main ways you’ll install R packages:
install.packages("packagename"))devtools::install_github("username/packagename") or remotes::install_github("username/packagename")There are other ways to get R packages. For example, Bioconductor provides “open source software for bioinformatics” and is home to a lot of great software & packages for genomics work. See the Bioconductor installation page for information about how to install their packages & tools.
What about Python packages? We haven’t dealt with this at all, since most of the packages we need are installed with Anaconda. What id we need something different?
Follow along to make sure you can install packages with pip:
Open the Terminal
Check for pip by running the appropriate commands (make sure you choose the correct one for your OS)
Install pip if necessary (follow along with commands here)
Install Python packages by running the following in the command line, replacing “project_name”:
python3 -m pip install "project_name"py -m pip install "project_name"For example, I have pip installed. To install the seaborn package, I run: python3 -m pip install "seaborn"
Go ahead and try it: Use the appropriate command for your OS to install the Python version of the palmerpenguins package (also named palmerpenguins).
Use and responsibility
Anyone can make and publish an R or Python package to GitHub. That means it is your responsibility to determine if it is safe to use for your science. Once you’ve sufficiently explored, vetted and tested a package, the next step is to use it.
Citing packages
Open source software is a major contribution to science and deserves to be credited. Often, developers are creating and maintaining packages voluntarily and without pay or support. Please remember to cite software and tools that helps you with your work.
Luckily, packages (at least those on CRAN) come with a built-in citation to make it easy for you to see the citation. Use citation(packagename) to return package citation.
For example:
citation("palmerpenguins")To cite palmerpenguins in publications use:
Horst AM, Hill AP, Gorman KB (2020). palmerpenguins: Palmer
Archipelago (Antarctica) penguin data. R package version 0.1.0.
https://allisonhorst.github.io/palmerpenguins/. doi:
10.5281/zenodo.3960218.
A BibTeX entry for LaTeX users is
@Manual{,
title = {palmerpenguins: Palmer Archipelago (Antarctica) penguin data},
author = {Allison Marie Horst and Alison Presmanes Hill and Kristen B Gorman},
year = {2020},
note = {R package version 0.1.0},
doi = {10.5281/zenodo.3960218},
url = {https://allisonhorst.github.io/palmerpenguins/},
}And an aside: R Markdown now makes it pretty nice to interface between your Zotero libraries and your R Markdown document. In the Visual Editor mode for R Markdown, you should see an `@` symbol. Clicking on it (if you have Zotero) should bring up your libraries.
Select the citation you want to insert. For example, data used in this example are from the palmerpenguins R package (Horst, Presmanes Hill, and Gorman 2020).
Click on the citation text for a preview of how it’ll show up in the final document!
Learn more about citations in R Markdown: https://blog.rstudio.com/2020/11/09/rstudio-1-4-preview-citations/
End Interactive Session 5
Extra unit test example (from lecture):
fish_mass <- function(alpha, beta, fish_length) {
if (fish_length <= 0) {
stop("Fish length must be greater than or equal to 0.")
}
mass <- alpha * (fish_length ^ beta)
return(mass)
}
# Unit tests
expect_length(fish_mass(alpha = 1, beta = 2, fish_length = 3), 2)Error: fish_mass(alpha = 1, beta = 2, fish_length = 3) has length 1, not length 2.Extra Roxygen example:
#' Calculate the distance traveled in miles
#'
#' @param rate a number in units of miles per hour
#' @param time a number in units of hours
#'
#' @return a number indicating the distance traveled in miles
#' @export
#'
#' @examples
distance <- function(rate, time) {
dist <- rate * time
return(dist)
}