# data read in
library(here)
toad_path <- here("data/cane-toad-wildnet-to-1999.parquet")
library(arrow)
toads_raw <- read_parquet(file = toad_path)
head(toads_raw)
#> scientificName decimalLatitude decimalLongitude eventDate
#> 1 Rhinella marina -26.75672 153.1261 1992-07-06
#> 2 Rhinella marina -23.15955 150.4811 1992-03-01
#> 3 Rhinella marina -25.28674 151.9127 1996-09-09
#> 4 Rhinella marina -27.14006 153.0594 1994-05-21
#> 5 Rhinella marina -19.36499 143.1678 1972-01-01
#> 6 Rhinella marina -27.59299 152.5157 1996-05-18
#> coordinateUncertaintyInMeters dataResourceName
#> 1 1300 WildNet - Queensland Wildlife Data
#> 2 450 WildNet - Queensland Wildlife Data
#> 3 500 WildNet - Queensland Wildlife Data
#> 4 1300 WildNet - Queensland Wildlife Data
#> 5 3600 WildNet - Queensland Wildlife Data
#> 6 100 WildNet - Queensland Wildlife Data
#> basisOfRecord license year
#> 1 OCCURRENCE CC-BY 4.0 (Int) 1992
#> 2 OCCURRENCE CC-BY 4.0 (Int) 1992
#> 3 OCCURRENCE CC-BY 4.0 (Int) 1996
#> 4 OCCURRENCE CC-BY 4.0 (Int) 1994
#> 5 OCCURRENCE CC-BY 4.0 (Int) 1972
#> 6 OCCURRENCE CC-BY 4.0 (Int) 19961 Why pipelines?
In this section we will discuss why pipelines are a useful idea. We will use an example analysis of cane toad data to help motivate this.
Many Australians, particularly those in Queensland, are familiar with the almost-folklore tale of the cane toad. For the unfamiliar, cane toads were introduced to Australia in 1935 as a form of pest control, to eat the cane beetle, who were eating sugar cane. The cane toad was not very good at eating the cane beetle, who lived near the top of the sugar cane, and the cane toads were poor climbers. The cane toad was a prolific breeder, and quickly spread across wider Queensland. It has severely impacted local fauna, as it has poison glands that will kill most animals. Even it’s tadpoles are highly toxic.
So, this analysis follows a story of using various releases of existing cane toad data. Our aim of the analysis is to plot where they are, and work out how fast they are moving.
Overview
Duration 90 minutes
Questions
- Why would I want a pipeline at all?
- Where were the cane toads, and how far did they move?
- What changes when my script becomes a document?
- What happens when more data arrives?
- How do I know if my outputs are current?
What you need this session
- A session of RStudio open
- The packages from Setup
- The script and data I send you at the start of the session
Photo by Brian Gratwicke Found on wikimedia, licensed under CC BY 2.0
NoteYour Turn: Discussion on pipelines and reproducibility
Let’s discuss:
- What do you know about reproducibility? What do you think of when you hear “reproducibility”?
- Have you heard of pipelines before?
- What sort of issues have you encountered with reproducibility?
- What are you hoping to get out of this course?
- What are your concerns with applying this material?
1.1 Why pipelines?
Data analysis is iterative. Cleaning, exploring and modelling rarely run start to finish without problems. Every step forward you take, sometimes takes you a couple back. E.g., finding a problem in a plot takes you to a problem with data cleaning.
I find a common pattern then is to re-run the code - from the top. It’s easier to just start “carte blanche” each time, because it’s too hard to work out hiwch parts of the analysis need re-running.
This can be fine. If you’re lucky! Often times your analysis can be so big - in terms of taking time, or managing complexity, that re-running it costs real time.
So then you start saving the intermediate steps. Now you have another problem…how do you know when to re-run one part of the analysis, which other parts depend upon changes to another?
This class of problems is what pipeline tools are designed to solve. They watch which files change, and which part of the code depends on those files. Then, they can run only the parts that depend upon them.
In R, targets is one of these pipeline toolkits, and it’s what we will be learning how to use.
The rest of this chapter is mostly hands-on programming. We will start with some toad data.
1.2 Live coding
- Reading in cane toad data
- Exploring it
- Save the script
1.3 The question
Toads were released at one place, on one date, and they spread west. So how far west did the edge get, decade by decade, and how far did it move each time?
Somebody has emailed you some files: A script, and some data.
1.4 Replicate it
Below you can view the entire analysis that I will live code:
NoteThe live code
Here’s the script, a piece at a time. The code below is read straight out of ch1/01-flat-script.R, so what you see here is what’s in the file.
One row is one report of one species, at one place, on one date. It isn’t a survey, and it isn’t a count. Somebody saw a toad and wrote it down.
The column names are long, and they are in camel case. {janitor} tidies them up.
toads_raw |> names()
#> [1] "scientificName" "decimalLatitude"
#> [3] "decimalLongitude" "eventDate"
#> [5] "coordinateUncertaintyInMeters" "dataResourceName"
#> [7] "basisOfRecord" "license"
#> [9] "year"
library(janitor)
toads_raw |>
clean_names() |>
names()
#> [1] "scientific_name" "decimal_latitude"
#> [3] "decimal_longitude" "event_date"
#> [5] "coordinate_uncertainty_in_meters" "data_resource_name"
#> [7] "basis_of_record" "license"
#> [9] "year"Then shorter names for the columns we use most, and a year and a decade to group by.
library(tidyverse)
toads <- toads_raw |>
clean_names() |>
rename(
lat = decimal_latitude,
lon = decimal_longitude,
date = event_date,
coord_var_m = coordinate_uncertainty_in_meters,
resource_name = data_resource_name
) |>
mutate(
year = year(date),
decade = floor(year / 10) * 10,
.after = date
)
head(toads)
#> # A tibble: 6 × 10
#> scientific_name lat lon date decade coord_var_m
#> <chr> <dbl> <dbl> <dttm> <dbl> <dbl>
#> 1 Rhinella marina -26.8 153. 1992-07-06 00:00:00 1990 1300
#> 2 Rhinella marina -23.2 150. 1992-03-01 00:00:00 1990 450
#> 3 Rhinella marina -25.3 152. 1996-09-09 00:00:00 1990 500
#> 4 Rhinella marina -27.1 153. 1994-05-21 00:00:00 1990 1300
#> 5 Rhinella marina -19.4 143. 1972-01-01 00:00:00 1970 3600
#> 6 Rhinella marina -27.6 153. 1996-05-18 00:00:00 1990 100
#> # ℹ 4 more variables: resource_name <chr>, basis_of_record <chr>,
#> # license <chr>, year <dbl>What have we got?
Before trusting anything that comes out of this, look at it.
# EDA
glimpse(toads)
#> Rows: 4,494
#> Columns: 10
#> $ scientific_name <chr> "Rhinella marina", "Rhinella marina", "Rhinella marina…
#> $ lat <dbl> -26.75672, -23.15955, -25.28674, -27.14006, -19.36499,…
#> $ lon <dbl> 153.1261, 150.4811, 151.9127, 153.0594, 143.1678, 152.…
#> $ date <dttm> 1992-07-06, 1992-03-01, 1996-09-09, 1994-05-21, 1972-…
#> $ decade <dbl> 1990, 1990, 1990, 1990, 1970, 1990, 1960, 1970, 1990, …
#> $ coord_var_m <dbl> 1300, 450, 500, 1300, 3600, 100, 1800, 1800, 1300, 130…
#> $ resource_name <chr> "WildNet - Queensland Wildlife Data", "WildNet - Queen…
#> $ basis_of_record <chr> "OCCURRENCE", "OCCURRENCE", "OCCURRENCE", "OCCURRENCE"…
#> $ license <chr> "CC-BY 4.0 (Int)", "CC-BY 4.0 (Int)", "CC-BY 4.0 (Int)…
#> $ year <dbl> 1992, 1992, 1996, 1994, 1972, 1996, 1964, 1974, 1992, …
library(visdat)
vis_dat(toads)
That’s how often somebody wrote a toad down, not how many toads there were. More records in the 1990s mostly means more people looking.
Where are they?
First, just the points. Then the points on a map of Australia.
# What's in ozmap_states
ozmap_states
#> Simple feature collection with 9 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 105.5507 ymin: -43.63203 xmax: 167.9969 ymax: -9.229287
#> Geodetic CRS: GDA94
#> # A tibble: 9 × 2
#> NAME geometry
#> * <chr> <MULTIPOLYGON [°]>
#> 1 New South Wales (((150.7016 -35.12286, 150.6611 -35.11782, 150.6…
#> 2 Victoria (((146.6196 -38.70196, 146.6721 -38.70259, 146.6…
#> 3 Queensland (((148.8473 -20.3457, 148.8722 -20.37575, 148.85…
#> 4 South Australia (((137.3481 -34.48242, 137.3749 -34.46885, 137.3…
#> 5 Western Australia (((126.3868 -14.01168, 126.3625 -13.98264, 126.3…
#> 6 Tasmania (((147.8397 -40.29844, 147.8902 -40.30258, 147.8…
#> 7 Northern Territory (((136.3669 -13.84237, 136.3339 -13.83922, 136.3…
#> 8 Australian Capital Territory (((149.2317 -35.222, 149.2346 -35.24047, 149.271…
#> 9 Other Territories (((167.9333 -29.05421, 167.9188 -29.0344, 167.93…They’re all in one state, so zoom in on Queensland.
You can watch them come down the coast and head inland.
Find the front
The front is the westernmost record in each decade. Longitude gets smaller as you go west, so that’s the smallest lon.
How far did it move?
{geodist} measures distances across the surface of the earth, in metres. It takes a matrix of coordinates.
sequential = TRUE measures from each row to the next, rather than every row to every other row. There’s no step into the first decade, so pad = TRUE puts an NA there.
# How far is it from one decade's edge to the next?
library(geodist)
# use {geodist} to calculate the distance
# it takes a matrix of inputs:
dist_mat <- cbind(lon = toad_west_front$lon, lat = toad_west_front$lat)
dist_mat
#> lon lat
#> [1,] 145.3678 -16.46512
#> [2,] 145.1178 -17.68165
#> [3,] 144.3178 -18.14832
#> [4,] 139.5511 -17.74832
#> [5,] 139.1511 -17.03165
#> [6,] 138.0679 -17.18183
#> [7,] 138.5428 -18.69089
# geodist() walks down the rows and measures each step across the surface of the
# earth, in metres. `sequential = TRUE` is what makes it row-to-row rather than
# every-pair. There is no step into the first decade, so we pad it out:
distances_m <- dist_mat |>
geodist(measure = "geodesic", sequential = TRUE, pad = TRUE)
distances_m
#> [1] NA 137238.08 99261.06 506881.80 89987.55 116477.00 174435.68
toad_speed <- toad_west_front |>
mutate(
distance_km = distances_m / 1000,
.after = scientific_name
)
toad_speed
#> # A tibble: 7 × 11
#> scientific_name distance_km lat lon date decade coord_var_m
#> <chr> <dbl> <dbl> <dbl> <dttm> <dbl> <dbl>
#> 1 Rhinella marina NA -16.5 145. 1937-01-01 00:00:00 1930 10000
#> 2 Rhinella marina 137. -17.7 145. 1943-01-01 00:00:00 1940 1800
#> 3 Rhinella marina 99.3 -18.1 144. 1950-01-01 00:00:00 1950 1800
#> 4 Rhinella marina 507. -17.7 140. 1968-01-01 00:00:00 1960 3600
#> 5 Rhinella marina 90.0 -17.0 139. 1979-01-01 00:00:00 1970 1800
#> 6 Rhinella marina 116. -17.2 138. 1983-10-01 00:00:00 1980 900
#> 7 Rhinella marina 174. -18.7 139. 1996-08-14 00:00:00 1990 1000
#> # ℹ 4 more variables: resource_name <chr>, basis_of_record <chr>,
#> # license <chr>, year <dbl>
ggplot(toad_speed, aes(x = decade, y = distance_km)) +
geom_col()
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_col()`).
The 1960s are the big one, at about 507 km.
Look closely at the 1990s, though. The westernmost record is further east than it was in the 1980s, and the distance is still positive. A distance doesn’t know which way it went.
NoteYour Turn
- Run the script top to bottom. Do you get the same distances?
- Which decade moved furthest? Find it on the map.
- The front is a single record each decade. What could go wrong with that?
1.5 From a script to a document
Someone wants to read this! Can you convert this into a quarto document?
So put it in a Quarto document. Each part of the script goes in a chunk, with a sentence or two around it saying what it does.
NoteYour Turn
Turn 01-flat-script.R into a Quarto document, and render it.
1.6 Tidy it up
Now make it read like a document.
- Libraries at the top. Put all
library()calls at the top. - Tidy up the iterations. The
View(), thenames()beforeclean_names(), the plot without a map, the plot with the whole of Australia. Those were you working it out. They were useful at the time.
Let’s make this a shorter, tidy document.
NoteYour Turn
- Move every
library()call into one chunk at the top. - Remove the steps that were you exploring.
- Render it again. Does it still give the same distances?
1.7 Someone wants the data
Your work is getting noticed!
A colleague asks you for the data with the distances added.
Sure, you say, you can add this to the quarto file?
Render, and the CSV is there.
NoteYour Turn
Add the write_csv() chunk and render.
- Change a single word in the text and render again. Did the CSV change?
- If your colleague opened the CSV tomorrow, how would they know which render wrote it?
1.8 Rendering a document to write a dataset?
It’s a lot!
But something about this feels a bit strange? You are rendering a document, and you get:
- plots, HTML file
- a written CSV file
Creating the CSV is a side effect of rendering.
It’s only as current as your last render, and renders happen for reasons that have nothing to do with the data:
Fix a typo in a sentence? CSV gets rewritten. If the render fails halfway, the file might get re-written? Depending on where it is?
And like, how do you know which version of the code produced it? How do you know it is up to date? Can we trust that this data came from the new data? Or is it the old one? Should we just run it again to be sure?
I don’t think that rendering a document to save data a good idea.
I’m not saying don’t save data! That’s a good thing.
But I am saying that writing a file should be a step that something keeps track of, rather than something that happens on the way past.
1.9 How else could you lay this project out?
Is there another way you could lay the project out?
One long script. That’s what you were sent. It works, and it’s where everybody starts.
Numbered scripts. Split it into 01-read.R, 02-clean.R, 03-front.R and 04-distance.R.
- Better! The numbers give the project an entry point. Someone new knows where to start.
- But they don’t tell you what needs running again. Change
02-clean.R, and you have to remember that 03 and 04 depend on it.
One Quarto document.
- Everything in one place, the report comes for free.
Functions, with a driver script.
- R functions go in
R/ - A short
analysis.Rcalls them in order
toad-analysis/
├── data/ what arrived, never edited by hand
├── R/ functions: code that gets used
├── analysis.R code that gets run, in order
└── output/ anything you could delete and make againThis is my favourite. We’ll talk about it more later.
How to pick a good one
Your system is good, if you can answer this:
Can you delete output/ and make all of it again by running your code?
For every layout above, the answer is yes. With a caveat: only you know how.
Your project might not be perfectly organised, but it should be possible for someone to understand how it is structured, and what to run first.
NoteRead more
This section is adapted from the project organisation chapter of R Best Practices, which goes through each of these layouts, and a few more, in detail.
1.10 More data arrives
Your boss (me) emails you the records up to 2010, cane-toad-wildnet-to-2010.parquet.
Can you run it again?
1.11 Stop
Stop here, and look at what’s on your disk.
You have a script, a document, a rendered report, and a CSV. Some of them were made from the 1999 data and some from the 2010 data.
NoteYour Turn
In pairs, and without running anything.
- Does
output/toad-speed.csvhave seven decades in it, or nine? - Is the rendered report from the same run as the CSV?
- How would you find out, without rendering everything again?
Nothing errored at any point. Every file looks fine. You just can’t tell which of them are current.
1.12 What does a pipeline look like?
Here’s one I prepared earlier.
ch1/ch1-targets/, written as a targets pipeline.
# To run:
# Sys.setenv(TAR_PROJECT = "ch1")
# tar_make()
source(here::here("ch1/ch1-targets/packages.R"))
source(here::here("ch1/ch1-targets/functions.R"))
tar_assign({
toad_path <- here("data/cane-toad-wildnet-to-1999.parquet") |> tar_file()
toads_raw <- read_parquet(file = toad_path) |> tar_target()
toads <- tidy_toads(toads_raw) |> tar_target()
# Let's just look at qld
qld <- ozmap_states |>
filter(NAME == "Queensland") |>
tar_target()
# ---- find the front ---------------------------------------------------------
# Find the most western toads per decade
toad_west_front <- toads |>
group_by(decade) |>
# smallest lon == most westerly
slice_min(lon, n = 1) |>
ungroup() |>
tar_target()
toad_speed <- toad_west_front |>
add_distance(lon = "lon", lat = "lat") |>
tar_target()
report <- tar_quarto(path = "ch1/ch1-targets/tar-ch-1.qmd")
})Each step is a target, with a name.
The targets package works out which ones depend on which, and when something changes it reruns only those.
You never have to work out what to run first. The answer is always tar_make().
Over the rest oif the course we will focus on these key ideas:
Functions. Each step of the analysis does one job.
Results worth keeping. Each target gets saved, and is only built again when something it depends on changes.
The dependency graph. Which steps depended on the file you changed? That question has an answer, it’s a real object, and
targetswill draw it for you.
TipWhere we are going
By the end of this course you’ll change one line in a cleaning function, run tar_make(), and watch the things that depend on it rebuild while everything else stays put.
1.13 But first: your code
We won’t jump straight into targets, not yet!
There are mechanics to it, and tools of the trade, and introducing all of it now would be too much at once.
So let’s start somewhere smaller. Next session we take the script apart and rebuild it as functions.
Summary
- You replicated an analysis from a script and a data file.
- A document is for reading. Libraries at the top, and the exploring taken out.
- A document that writes files has side effects you cannot see from the files.
- Functions, results worth keeping, and the dependency graph are the three things this course is about.
NoteWhere this comes from
The section on laying out an analysis draws on R Best Practices, which covers project organisation, file paths and naming in more depth.
The function-writing and debugging sessions draw on Introduction to Functions and Debugging, which covers both in more depth.