Lecture 3: Functions

If we only had one data set to analyze, it would probably be faster to load the file into a spreadsheet and use that to plot some simple statistics. But we have twelve files to check, and may have more in the future. In this lesson, we’ll learn how to write a function so that we can repeat several operations with a single command.

Defining a Function

You can write your own functions in order to make repetitive operations using a single command. Let’s start by defining your function “my_function” and the input parameter(s) that the user will feed to the function. Afterwards you will define the operation that you desire to program in the body of the function within curly braces ({}). Finally, you need to assign the result (or output) of your function in the return statement.

Now let’s see this process with an example. We are going to define a function fahrenheit_to_celsius that converts temperatures from Fahrenheit to Celsius:

fahrenheit_to_celsius <- function(temp_F) {
  temp_C <- (temp_F - 32) * 5 / 9
  return(temp_C)
}

We define fahrenheit_to_celsius by assigning it to the output of function. The list of argument names are contained within parentheses. Next, the body of the function–the statements that are executed when it runs–is contained within curly braces ({}). The statements in the body are indented by two spaces, which makes the code easier to read but does not affect how the code operates.

When we call the function, the values we pass to it are assigned to those variables so that we can use them inside the function. Inside the function, we use a return statement to send a result back to whoever asked for it.

Automatic Returns

In R, it is not necessary to include the return statement. R automatically returns whichever variable is on the last line of the body of the function. While in the learning phase, we will explicitly define the return statement.

Let’s try running our function. Calling our own function is no different from calling any other function:

# freezing point of water
fahrenheit_to_celsius(32)

## [1] 0

# boiling point of water
fahrenheit_to_celsius(212)

## [1] 100

# pizza oven temperature
fahrenheit_to_celsius(700)

## [1] 371.1111

We’ve successfully called the function that we defined, and we have access to the value that we returned.

Exercise: Composing Functions

Now that we’ve seen how to turn Fahrenheit into Celsius, it’s easy to turn Celsius into Kelvin: Hint: temp_C + 273.15 Calculate freezing point of water in Kelvin

## [1] 273.15

What about converting Fahrenheit to Kelvin? We could write out the formula, but we don’t need to. Instead, we can compose the two functions we have already created:

fahrenheit_to_kelvin <- function(temp_F) {
  temp_C <- fahrenheit_to_celsius(temp_F)
  temp_K <- celsius_to_kelvin(temp_C)
  return(temp_K)
}

# freezing point of water in Kelvin
fahrenheit_to_kelvin(32.0)

## [1] 273.15

This is our first taste of how larger programs are built: we define basic operations, then combine them in ever-larger chunks to get the effect we want. Real-life functions will usually be larger than the ones shown here–typically half a dozen to a few dozen lines–but they shouldn’t ever be much longer than that, or the next person who reads it won’t be able to understand what’s going on.

Exercise: Named Variables and the Scope of Variables

Functions can accept arguments explicitly assigned to a variable name in the function call functionName(variable = value), as well as arguments by order:

input_1 <- 20
mySum <- function(input_1, input_2 = 10) {
  output <- input_1 + input_2
  return(output)
}

1. Given the above code was run, which value does mySum(input_1 = 1, 3) produce?
   1. 4
   2. 11
   3. 23
   4. 30
2. If mySum(3) returns 13, why does mySum(input_2 = 3) return an error?

Testing, Error Handling, and Documenting

Once we start putting things in functions so that we can re-use them, we need to start testing that those functions are working correctly. To see how to do this, let’s write a function to center a dataset around a particular midpoint:

center <- function(data, midpoint) {
  new_data <- (data - mean(data)) + midpoint
  return(new_data)
}

We could test this on our actual data, but since we don’t know what the values ought to be, it will be hard to tell if the result was correct. Instead, let’s create a vector of 0s and then center that around 3. This will make it simple to see if our function is working as expected:

z <- c(0, 0, 0, 0)
z

## [1] 0 0 0 0

center(z, 3)

## [1] 3 3 3 3

Dataset

We are studying inflammation in patients who have been given a new treatment for arthritis, and need to analyze the first dozen data sets. The data sets are stored in comma-separated values (CSV) format. Each row holds the observations for just one patient. Each column holds the inflammation measured in a day, so we have a set of values in successive days.

That looks right, so let’s try center on our real data. We’ll center the inflammation data from day 4 around 0: download this

dat <- read.csv(file = "data/inflammation-01.csv", header = FALSE)
centered <- center(dat[, 4], 0)
head(centered)

## [1]  1.25 -0.75  1.25 -1.75  1.25  0.25

It’s hard to tell from the default output whether the result is correct, but there are a few simple tests that will reassure us:

dat <- read.csv(file = "data/inflammation-01.csv", header = FALSE)
# original mean
mean(dat[, 4])

## [1] 1.75

# centered mean
mean(centered)

## [1] 0

That seems right: the original mean was about 1.75 and the mean of the centered data is 0. We can even go further and check that the standard deviation hasn’t changed:

# original standard deviation
sd(dat[, 4])

## [1] 1.067628

# centered standard deviation
sd(centered)

## [1] 1.067628

Those values look the same, but we probably wouldn’t notice if they were different in the sixth decimal place. Let’s do this instead:

# difference in standard deviations before and after
sd(dat[, 4]) - sd(centered)

## [1] 0

Sometimes, a very small difference can be detected due to rounding at very low decimal places. R has a useful function for comparing two objects allowing for rounding errors, all.equal:

all.equal(sd(dat[, 4]), sd(centered))

## [1] TRUE

It’s still possible that our function is wrong, but it seems unlikely enough that we should probably get back to doing our analysis. However, there are two other important tasks to consider: 1) we should ensure our function can provide informative errors when needed, and 2) we should write some documentation for our function to remind ourselves later what it’s for and how to use it.

dat[, 4] <- centered
write.table(dat, file = "inflammation-centered.csv", col.names = FALSE, row.names = FALSE, quote = FALSE, sep = ",")

Documentation

A common way to put documentation in software is to add comments like this:

center <- function(data, midpoint) {
  # return a new vector containing the original data centered around the
  # midpoint.
  # Example: center(c(1, 2, 3), 0) => c(-1, 0, 1)
  new_data <- (data - mean(data)) + midpoint
  return(new_data)
}

Writing Documentation

Formal documentation for R functions is written in separate .Rd using a markup language similar to LaTeX. You see the result of this documentation when you look at the help file for a given function, e.g. ?read.csv. The roxygen2 package allows R coders to write documentation alongside the function code and then process it into the appropriate .Rd files. You will want to switch to this more formal method of writing documentation when you start writing more complicated R projects.

Exercise: Functions to Create Graphs

Write a function called analyze that takes a filename as an argument and displays the three graphs produced in the (average, min and max inflammation over time). analyze("data/inflammation-01.csv") should produce the graphs already shown, while analyze("data/inflammation-02.csv") should produce corresponding graphs for the second data set. Be sure to document your function with comments. download this

Saving Plots to a File

So far, we have built a function analyze to plot summary statistics of the inflammation data:

analyze <- function(filename) {
  # Plots the average, min, and max inflammation over time.
  # Input is character string of a csv file.
  dat <- read.csv(file = filename, header = FALSE)
  avg_day_inflammation <- apply(dat, 2, mean)
  plot(avg_day_inflammation)
  max_day_inflammation <- apply(dat, 2, max)
  plot(max_day_inflammation)
  min_day_inflammation <- apply(dat, 2, min)
  plot(min_day_inflammation)
}

And also built the function analyze_all to automate the processing of each data file:

analyze_all <- function(folder = "data", pattern) {
  # Runs the function analyze for each file in the given folder
  # that contains the given pattern.
  filenames <- list.files(path = folder, pattern = pattern, full.names = TRUE)
  for (f in filenames) {
    analyze(f)
  }
}

While these are useful in an interactive R session, what if we want to send our results to our collaborators? Since we currently have 12 data sets, running analyze_all creates 36 plots. Saving each of these individually would be tedious and error-prone. And in the likely situation that we want to change how the data is processed or the look of the plots, we would have to once again save all 36 before sharing the updated results with our collaborators.

Here’s how we can save all three plots of the first inflammation data set in a pdf file:

pdf("inflammation-01.pdf")
analyze("data/inflammation-01.csv")
dev.off()

The function pdf redirects all the plots generated by R into a pdf file, which in this case we have named “inflammation-01.pdf”. After we are done generating the plots to be saved in the pdf file, we stop R from redirecting plots with the function dev.off.

Overwriting Plots

If you run pdf multiple times without running dev.off, you will save plots to the most recently opened file. However, you won’t be able to open the previous pdf files because the connections were not closed. In order to get out of this situation, you’ll need to run dev.off until all the pdf connections are closed. You can check your current status using the function dev.cur. If it says “pdf”, all your plots are being saved in the last pdf specified. If it says “null device” or “RStudioGD”, the plots will be visualized normally. {: .callout}

We can update the analyze function so that it always saves the plots in a pdf. But that would make it more difficult to interactively test out new changes. It would be ideal if analyze would either save or not save the plots based on its input.