R the Software Engineering Way: Chapter 1

Chapter Zero

Chapter Two

Chapter 1: Writing programs in R

In the last chapter, we've worked on setting up a basic R project with a development environment, a development container, version control and all the tools that we need to write reproducible, consistent R code. We also talked a bit about command line tools and general principles for how to learn tech stuff. What we didn't do, however, is actually write any R, so in this chapter we're going to learn how to use R to write some basic programs.

In my experience, statisticians without programming experience tend to write R as long scripts of, basically, one instruction after another, almost all built-in. This works pretty well for doing stats: you have a nice long list of everything you've done in one spot, and without any real need for control flow in most statistical applications, it's not too bad from the perspective of understanding the code. However, it's an absolute disaster if you want to reuse the code you've written for something else, you need code that does different things based on different conditions (not too important if you're just using stats packages, very important for most software engineering applications) or you want to let other people use your code for their own purposes. In all of these cases, being able to factor your code into functions and structure it effectively becomes quite important.

This chapter then aims to cover two important programming concept that software engineers use all the time but that statisticians tend not to use so much: functions and control flow. We'll cover what a function in R actually is and how it works, what good functions look like, and discuss conditionals, for loops and while loops. In the process, we'll also discuss why we organise our code as packages, learn how to write unit tests for code and why we should do that, and finally we'll show you how to actually use the development environment that we spent the last chapter building.

The commit for this chapter can, as previously, be found at https://gitlab.com/irishenceaway/r.the.software.way: the commit for this chapter is tagged chapter1.

Organising your code: packages

You'll remember that in the last chapter I had you install the devtools package for R and scaffold your project directory as an R package. I didn't specify why we did that at the time, but we should probably look to explain it now.

First and foremost, structuring your work as a package makes you write the bulk of your code as functions. This is best practice in almost all languages, but in R, especially if you've trained as a statistician, it's easy to get into the habit of writing ad-hoc scripts that aren't reusable or reproducible. Structuring your code as a package means that you have to write better code, and it's also much more like other programming languages than the usual script-based approach. It also has the advantage of enforcing a canonical project structure: all of your functions go into project files in the /R directory, all of your data lives in /data and so on and so forth.

Structuring your work as a package also means that you can build it and distribute it to other people if you find that the package has general use. This might be on CRAN, or it might be via a public repository or any number of other options: packages, in general are quite a bit more portable than scripts.

Finally, the development tools for packages are much more developed than the tools available for ad-hoc scripting (where RStudio is basically your only option). We've already used the devtools package, which provides scaffolding tools for packages and also allows us to load, document and test packages from the interactive R shell. Unit testing in particular is vital for maintaining reasonable levels of code quality, and it's basically only possible when working with a package framework.

The general rule for packages is that all the R code in a package should be in the /R directory, and that it should consist entirely and only of function definitions. That's probably about all you need to know for now, though: you'll pick up more as we work through the rest of the chapter.

Functions

You will probably have found that in the process of writing R code, you often reuse very similar code that encapsulates the same logic, applied to slightly different sets of data each time. A common principle in software engineering is DRY (don't repeat yourself), and this is because, if you're copy-pasting boilerplate code with different changes each time, it has a whole lot of undesirable effects. In the first instance, it clutters your files: I remember one nightmare scenario where a codebase had 20,000 lines of R in one file, most of which was one code block with a few variables changed each time. Having code disorganised in this way also makes it much easier for bugs to find their way into the code: after all, if you're editing each of the repeated blocks manually, all it takes is for you to forget to make one change, and your data is all of a sudden in danger of being corrupted. Ideally, you'd want to reuse the same code each time, changing the arguments as needed, and the way you do that is with functions.

A function in R is rather broader in scope than the mathematical definition you're probably familiar with, which is a bijective map between a domain and a codomain. While you can write functions in R that behave like that (and wherever possible, you should), a function in R is simply a block of code that's assigned to a name. In its most general form, a function definition in R might look like this:

heaviside_step <- function (x) {
  if (x < 0){
    return(0)
  }else {
    return(1)
  }
}

This is a simple implementation of the Heaviside step function. Looking at the implementation here, you can see that the name of the function is simply a variable name being assigned a value in the same way that you might read in a CSV as a dataframe and assign it to a variable. In software engineering, we'd describe this by saying that functions in R are first-class objects, which means that we can treat them more or less as any other data type and even use them as inputs to other functions. After the name and the assignment, we have the phrase function (x). This tells us that a) what we're assigning is a function and b) that this function has one argument, which is x. Functions in R can have an arbitrary number of arguments, and they can be any type of object: vectors, dataframes, matrices, strings, scalars and even other functions.

After that we have the body of the function, which is a simple block of code. You can put almost anything in here. In this case, we have an if/else statement (on which we'll expand more in our section on control flow) that returns zero if the argument x is less than zero and the value one if the argument is greater than zero. A function in R can have one or zero return values (in general, if a function has zero return values, R will return a value in any case, but there's nothing that says you have to use it or do anything with it). If a function has zero return values, it usually means that it's doing something that has side effects (more on this later). If you want to return multiple values from a function, by contrast, the usual way to do that in R is to bundle your values into a list or similar object and return that list as your return value.

R style guides tend to encourage you to use an implicit return statement rather than explicitly specifying the function's return value. A fair number of clever people believe in this, and while I'm loath to contradict the likes of Hadley Wickham, I personally think this is a mistake. It makes it difficult to reason about what any given function will do or return, it makes debugging a living hell and it is profoundly not how any other languages behave. An explicit return statement is, in general, a good idea.

Having written a function and assigned a name to it, you can then load the function into memory by executing the code. There are a couple of ways of doing this: in an interactive shell, you might just press enter. If you've written the function in a file, you can run source('filename.R'), or if your function is defined in an installed package, you can run library(package_name) to attach the package. Finally, if, as we're doing in this exercise, you're writing your code as a package which is currently not built, you can run devtools::load_all() to load your current package functions into memory.

Finally, having your function in memory and ready to use, you can call the function as follows:

z <- heaviside_step(5)

This will call the heaviside_step function with the argument x taking the value 5, calculate the return value using the code in the body of heaviside_step and assign the value to the variable z. And that's the basics of how functions in R work!

Pure functions and side effects

Unlike mathematical functions, R functions can have what we in software call "side effects". For example, you can write a function in R that returns no values, but that generates a plot in ggplot2 that then gets written to your filesystem. We say that this has side effects in opposition to what we'd call a "pure function", which is more or less the mathematical definition of a function. If you call a pure function with the same arguments, you'll get the same result each time. With our hypothetical plotting function, by contrast, what happens when you run the function depends on the existing state of the filesystem when you run it: if you already have a graph in the target directory, you'll have something different happen to what happens when the directory's empty. Functions with side effects are thus a bit of a pain to work with, and in general we should avoid them whenever possible: that said, it's basically impossible to do anything useful without introducing side effects, as basically any I/O operation or any interaction with an API or a filesystem introduces side effects. Best practice is generally to decouple code with side effects from any logic.

Variable scope in functions

As I'm writing this book for people who are generally already capable R users, I've not felt the need to do the "variables 101" thing in this book: you all know what a variable is, even if you might not be familiar with functions. However, now that we've introduced functions in a bit more depth than statisticians using R usually go into, it's worth going into a bit of depth as to what variables in R actually work and how they do what they do. To do that, let's examine the following piece of code:

phi <- 2

zeta <- function (x, y, z) {
  phi <- 3
  return((x + y + z) / phi)
}

print(zeta(2, 4, 8))

We've defined here a number of variables. The first definition of phi lives in the global scope: this is where variables live if you define them in the top-level environment. In this case, if you evaluate phi in the global environment, it will return 2. The second definition of phi, however, lives in the function scope. When a function executes, it creates an environment for itself as a child of the global environment that can have its own variable definitions independent of the global environment. That means that when the function executes, the value of phi will be set to 3, and the return value of the call below will be (2 + 4 + 8)/3, or 4.66... . If you were to delete the internal definition though, the function would find that phi isn't defined in function scope and look in the parent scope for a definition. It'd then find the value of phi in the global scope, and the final print statement would return 14/2, or 7.

A lot of statisticians use this kind of pattern to pass data between functions and global scope, but for a software engineer this causes real problems, simply because you don't know what the parent scope of the function is: if someone imports the function zeta into a script somewhere, they could set phi in the global scope to anything. Then you're fucked. In general, if you're using a variable in a function, it should always, always be a variable defined in the function scope. The only way you should ever pass external data to a function is through an argument, and the only way anything should ever leave a function is through a return statement (or, well, there are some exceptions for side effects, but the general principle stands).

Writing good functions

As a rule, you want to write small, self-contained functions, and when more complex behaviour is called for, build that up by composing multiple smaller functions rather than trying to do it all in one big one.

In general, you shouldn't write functions that don't return anything: the primary reason for this is that it makes it very difficult to write unit tests for them (more on this later). Of course, sometimes you need to write code that interacts with the machine or otherwise causes side effects, and that's inevitable, but that code should, wherever possible, be isolated from the rest of the codebase.

As stated before, functions should only ever communicate with other code through arguments and return values.

Functions should be testable: given an input, you should be able to look at a given function and work out what it's going to return so that you can write a unit test for them (more on this later). If you write code that's difficult to test, read and think about, it's easy to get into some real trouble.

Exercise

In this exercise you'll write a simple function for evaluating the values of an arbitrary quadratic polynomial of form f(x) = ax^2 + bx + c.

• Create a file in your project's /R directory and name it something like eval-quadratic.R. Use the command line for this.
• Define a function eval_quadratic that takes as arguments values x, a, b, c and returns the value of f(x). Set sensible default values for a, b, c.
• Start an interactive R shell, and use the load_all() function from devtools to make the function quadratic_val available in your interactive prompt and try it out with a few different values for x, a, b, c. Confirm that it behaves in the way you'd want it to.
• Using git and the command line, stage and commit your changes with an appropriate commit message.

Container interlude: using your development container

You'll probably have noticed a fairly glaring issue towards the end of the last exercise: we put all of that effort last week into building and running our container, but when we write code in our local repository, there's no obvious way of running the code in the container and you probably had to resort to running it locally. Given that I had you go to all the trouble of writing and building the Dockerfile, you would be reasonable in expecting that there's a solution to this. And there is!

The first option here might be to copy all the files into the container in the Dockerfile at build time. Now, there are circumstances where this might make sense: you might want to distribute a container image with everything in it so that another researcher can just pull the image from a container registry, run one command and have full access to all of your analysis. This is, in fact, a way in which I encourage researchers to distribute their findings. However, for development purposes this leaves a bit to be desired. After all, it means that you have to rebuild your image every time you change your source code, and as I'm sure you're aware by now, images take a moment to build. This goes double for R images, because R package installs are often very slow. We thus need a different solution to work with while we're developing the thing.

The way we actually solve this problem is with bind mounts. With a bind mount, you can bind a directory on your local machine to a directory in the container, temporarily overwriting the contents of the directory. You'll note that I had you create a directory in the Dockerfile called r-software-engineering: this is what it's for.

While it is possible to mount the directory by adding a flag to the docker run command you use to run the image (the exact command would be sudo docker run --rm -ti --volume .:/r-software-engineering iris-meredith/r-the-software-way), you've probably noticed that the command is already getting pretty long, and if you aren't careful about it, you can end up cluttering your system with unused containers, which is a pain. What we'll use, then, is the docker compose tool. We begin with this by creating a compose.yaml file in our project root with the following contents:

services:
  r-development-container:
    image: iris-meredith/r-the-software-way
    stdin_open: true
    tty: true
    volumes:
    - type: bind
      source: .
      target: /r-software-engineering

In here we've defined a service called r-development-container based on our local development image (it's worth noting that if you change the image, you'll need to rebuild it in the usual way before running this: it is possible to configure compose to build the image by default, but who has time to wait for the image to build before every time you want to use R). We specify that we want to keep standard input open in the image and allocate a pseudo-TTY to the container, attaching our local terminal to the standard I/O of the container (this is done by the -it flag in the docker run command). We then define a bind volume, binding the project root on our machine to the directory /r-software-engineering on our container. We can then simply run sudo docker compose run --rm r-development-container to spin up our development environment, without having to specify any of the more finicky options (we do still have to specify --rm to tell our system to get rid of the container once we disconnect from it, as unfortunately I've not found a way to make compose do it). Compose is a very powerful tool, and we'll use it for more things down the line (including its intended uses, even: this is a slight abuse of the tool), but for now, this lets you make use of the development container in the intended manner.

VSCodium allows you to automate a fair bit of this work by setting up a Development Container. That way, whenever you open the project, you'll automatically load into the container and execute everything in there. This has advantages for larger teams and in terms of convenience, but I've elected not to use the technology here (mostly because I like to understand how all of my stack works, and Dev Containers hide a lot of complexity). If you read the documentation for Dev Containers and like what you see, you might wish to set them up for yourself.

Exercise

In your project directory:

• Create a compose.yaml file with a service analogous to the one above, with your image names and suchlike replacing the ones I've used.
• Run sudo docker compose run -rm r-development container and confirm that it drops you into an R interactive shell as expected.
• Run devtools::load_all() in the interactive shell and confirm that you can use the quadratic function you just wrote and that it runs correctly.
• Write the mtcars dataset to a file named mtcars.csv or something similar, then quit the interactive shell. Confirm that mtcars.csv was successfully created in your root project directory. Delete the file afterwards (you don't want it cluttering your repository).
• Commit your changes using git, being sure to add a sensible commit message.

Unit testing

Having written your first function and being able to run it in your development container, the next step is to immediately write a unit test for the function. A unit test in this context is an automated test script that executes the function with given inputs and checks that the outputs are what they should be. I cannot stress enough the importance of having comprehensive unit tests for your code: it is quite possibly the single biggest thing that you can do to maintain the quality of your codebase and prevent it from disintegrating into the kinds of messes that we so often see.

It might not immediately be clear why you should automate tests rather than just test things ad-hoc: after all, loading the functions into your session and trying them out like that has been good enough so far. And yes, it does work, but it has a considerable number of ways in which it fails that makes it somewhat inadequate.

The first thing to note is that ad-hoc testing of this kind is something that you have to repeat every time you make a change to the function. I'm sure you've had the experience of looking back at some code you'd written in the past and having absolutely no idea what it does or how it works any more, and if you then need to modify that function, you can end up doing a lot of work to get back to your old understanding and possibly get it wrong and break things. The problem becomes significantly bigger if you've written other functions that use that function, because then if you accidentally change the expected outputs for your function, you might well break other parts of your code in the process.

Automated unit tests, in this case, enforce a given interface for your code: a given function with a given set of inputs will have to consistently produce the same set of outputs for tests to pass. You don't have to remember interfaces or what function calls what other function or anything: you can simply run the tests and identify breaking changes, and if the tests fail, either change the code so that they pass, or change the tests so that they accurately reflect the new interface. This makes maintaining or changing your code much faster down the line, saving you the time that it took to write the tests many times over.

Consistent unit testing also helps significantly with writing better code. It is, after all, in the nature of unit tests that they push towards the functions you've written being simple: you need to reason carefully about what the correct output for a given function input would actually be. If you're writing a function and find it hard to write tests for it, for example, it's probably a sign that the function is poorly-designed and you need to break it down into simpler pieces. It also encourages you to avoid code with side-effects where possible: while it is possible to test code with side effects, and I'll show you some techniques for doing so down the line, it's much more of a pain than testing pure functions.

This encourages you to break down your code into small, self-contained functions that are easy for people to understand, which massively improves the quality of the code you write. It also encourages you to isolate side effects, which as we've already discussed, makes it much much easier to identify where bugs actually are. Finally, unit tests make it much easier for multiple people to collaborate on the same codebase: you can make changes with confidence, knowing that unit tests will catch a lot of the worst disasters. Between all of these things, automated testing becomes an extremely useful barrier that prevents sloppiness and helps keep the codebase in much better shape than might otherwise be expected.

Every language tends to have its own frameworks and best practices for unit testing: in R, the relevant framework is called testthat, and it's included in the devtools package. To begin using the framework, start an R shell in your project directory, attach devtools and run usethis::use_test(). This will create a tests directory in your project directory, with one subdirectory called testthat (this is where your tests will go) and an R file called testthat.R that you shouldn't touch (it contains some scaffolding for automated testing). To write tests for your eval-quadratic function, create a file in the testthat subdirectory called test-eval-quadratic.R (the pattern for test names is that tests for a given file should be in a file called test-<filename>.R). The contents of a test file looks like this:

test_that("the quadratic evaluation function produces expected results", {
  expect_equal(eval_quadratic(0, 1, 0, 0), 0)
}
)

A test in R is a call to the function test_that, which takes as arguments a note about what the test is for, and then a set of assertions about what the function should output given a set of inputs. In this case, we're testing that f(x) = 1*x^2 should return a value of 0 for f(0). We can then run this test by running devtools::test() in an interactive R shell and get a set of results: all going well, the test will pass, and if the test fails, it will tell you where it failed and how, so that you can locate the bug and fix it. In general, you should write a test whenever you write a function or wherever you find a bug: this prevents the pattern of "fixing the same bug multiple times" that consistently happens when you write code without automated testing.

Exercise

In your project directory:

• Start an interactive R session using your development container and run usethis::use_test() to scaffold your test directory.
• Write a unit test for the quadratic function you built in the last section: use the example code as a model, and add some more expectations to get better coverage. Think carefully about potential ways in which the results the function gives could be wrong.
• Run the unit test using devtools::test and confirm that all of your tests pass.
• Commit your changes, being sure to add a coherent commit message.

Control flow

There are two core things that really distinguish the practice of writing a statistical script from the practice of writing programs as such. The first is the use of functions in the way we're using them now: as packets of functionality that go into a larger program rather than as tools that you call on their own in an interactive session. The second is control flow. In a statistical script, you can often get away without much control flow: you're often executing a pretty linear set of steps on some data where you don't need to really make decisions or repeat things. And that's all well and good for statistics. When writing software, however, control flow suddenly becomes quite important: programs often have to behave in different ways depending on external conditions, or repeat a process for a given number of times or until a given condition is met. Hence, while it's possible to get by with R without having an amazing grasp of control flow, if you want to write software in R you need to have a strong grasp of it. There are two major control flow constructs in R: conditionals and loops.

Conditionals

The conditional statement of relevance here is the if else block. In R, the syntax for that looks like this:

if (condition1) {
  statement1
} else if (condition2) {
  statement2
} else {
  statement3
}

Conditions in the above code are simply R code that evaluate to a TRUE or FALSE value. Then, if condition1 evaluates as being true, we execute statement1 and leave the block. If condition1 is false but condition2 is true, we execute statement2 and leave the block. If neither condition is true, we execute statement3 and exit.

The statements can do more or less anything you like: one obvious use for this might be switching between data sources based on some value. For example, if we've written a function that takes an arbitrary data source as an argument, we can write some code to read a .csv file if that's what the user passes, a .spss file if the user passes an SPSS dataset, or write some code to pull a table from a database if the user passes a database connection string. Using a conditional, in this case, allows you to write significantly more powerful and flexible code than would be possible with the usual scripting approach.

Loops

Loops in R aren't used frequently: this is in part because R is heavily vectorised so it's often easier to use vector operations. Thanks to the way in which R handles objects, it's also easy to write quite inefficient and slow R code using loops, so most authorities tend to discourage the use of loops. And again, for statistics, that makes sense. Unfortunately, even in stats, there are situations where nothing but a loop will do, especially in Bayesian statistics where such things as Markov Chain Monte Carlo rely on you sampling from a random process repeatedly. More prosaically, vectorised operations are in fact loops under the hood: just heavily optimised ones. And in software engineering more generally, you will find yourself using loops a lot, especially in less vectorised languages. It's thus a good idea to have a strong understanding of how loops work.

Loops, in essence, repeat a certain set of instructions until a given condition is met. R has two flavours (well, three, but the last one is weird and overcomplicates things) of loop: the for loop and the while loop.

for loops

The for loop in R runs a certain set of code over every element in a vector or list. The syntax for it is as follows:

for (x in collection) {
    statement(x)
}

In the simplest form, this is very similar to the apply family of functions in R: for each element in the collection, you do something to the element and then output the result. If you set collection <- 1:100, for example, you can run the same code over every integer between 1 and 100. Obviously, if that's all you're doing, you should probably just use apply. The for loop, however, allows you to do more powerful and sophisticated things. Take, for example, this code:

x_old <- NULL
for (x in collection) {
  if (!is.null(x_old)) {
    print(x - x_old)
  }
  x_old <- x
}

This is a short loop that calculates forward differences for a numerical vector. It can strictly speaking be vectorised, but it's a real pain in the behind to do so and I don't think it's a particularly natural way to express it compared to the loop approach. The loop is easier for a human to read, easier to test and generally less of a headache to work with.

The main pitfall with for loops in R is object allocation. When you create a vector in R, you have to allocate memory for that, which is a pretty expensive operation, but not in itself an issue. However, in other languages people will often construct a vector by doing something like this:

# Bad code: do not use
i_squared <- vector()
for (i in 1:1000){
  i_squared <- c(i_squared, i^2)
}

Because in R, when we create a new object from an old one, we make a copy, this code reallocates a whole new vector each time we run it, meaning that we repeat the slow process a thousand times, which makes for slow code. A lot of the perceived slowness of R as a language is a direct result of code like this. A more efficient way of writing the same loop might look like this:

i_squared <- rep(0, times=1000)
for (i in 1:1000){
  i_squared[i] <- i^2
}

By preallocating the vector before the loop, we only allocate memory once and then assign by element, drastically speeding the loop execution up. While you'll probably still use apply for most things that you could do with a loop, this is a useful tool to know.

while loops

The second kind of loop that R supports is the while loop, which repeats a certain set of instructions until a given condition is met. The syntax for this is as follows:

while (condition){
  statement
}

The condition here, as with our earlier discussion of conditional statements, is simply an R expression that evaluates to TRUE or FALSE: if the condition is true, the statement executes and then the condition is retested. This is useful for a wide range of processes that don't necessarily have a clear "do this for everything in the vector" kind of phrasing: for statistics in particular, they're especially useful for processes that involve some kind of convergence, or where you stop once some metric hits a certain level. An example of this might be approximating the exponential function via Taylor series:

# Code snippet for demonstration purposes: probably won't run without modification
while (estimate_new - estimate_old > 1e-6){
 estimate_new <- estimate_old + (x^i)/factorial(i)
 i <- i + 1
}

In this kind of situation, we don't know how many times the loop is going to run: we don't break out of the loop until the difference between consecutive estimates is less than 1e-6, or some other small number of our choice. Of course, if the process fails to converge or you write a condition that never evaluates to false, the loop will never terminate and your process will block forever: this is important to note. While loops are a lot more flexible and powerful than for loops, and they're a technique well worth having in your arsenal.

You can combine conditionals and loops to create some very complex behaviours: however, keeping the complexity of your functions to a minimum is generally a good idea. In general, you probably shouldn't nest loops and conditionals more than two or three deep: if you're hitting that point, consider factoring out some of the logic into a separate function.

Exercise

The first few parts of this exercise involve implementing the FizzBuzz problem, which is a popular question in software engineering interviews and tests. It's also very useful for teaching control flow. So, in your project, and using both version control and your development container:

• Write a function that takes an integer as an argument and returns the string "Fizz" if the integer is a multiple of 3, "Buzz" if it's a multiple of 5, "FizzBuzz" if it's a multiple of both, or returns the original integer if none of these cases hold. Write a unit test for it, taking care to consider all possible behaviours, and run the test to confirm that it behaves as you'd expect.
• Write a second function that calls the first function. This function should accept as an argument an integer and return a vector of strings, where the value of any given element corresponds to the return value of the first function for the index of that element. Write a unit test for that function, and run the tests, confirming that they pass.
• Using the quadratic function you wrote earlier, write an implementation of Newton's Method to find the roots of a quadratic equation. Write a unit test for the function and confirm that it passes. A few hints: functions are first-class objects, so you can pass them to other functions, and (because I'm not trying to teach you calculus) the derivative of a quadratic function is a linear function. You'll also need to choose how accurate you want your approximation to be. Don't worry too much about complex roots and suchlike for now: we'll deal with them later.
• Think about the Newton's method implementation that we wrote above. How could it fail? What cases does it not take into account, and how might you modify the implementation to do so (we'll have a go at doing that in future exercises).
• Commit your changes and push them to your remote repository.

Between chapters zero and one, you now have all the basic skills that you need to write code and do statistical analysis in a reproducible, robust way that a software engineer would be proud of. In the next chapter, we'll show you how all of this knowledge comes together in a non-trivial statistical analysis.