In order to understand functions better, let’s take a look at what they consist of.

Typing the name of a function shows you the code that runs when you call it. This is the rt function, which generates random numbers from a t-distribution:^[20]

rt

## function (n, df, ncp)
## {
##     if (missing(ncp))
##         .External(C_rt, n, df)
##     else rnorm(n, ncp)/sqrt(rchisq(n, df)/df)
## }
## <bytecode: 0x0000000019738e10>
## <environment: namespace:stats>

As you can see, rt takes up to three input arguments: n is the number of random numbers to generate, df is the number of degrees of freedom, and ncp is an optional noncentrality parameter. To be technical, the three arguments n, df, and ncp are the formal arguments of rt. When you are calling the function and passing values to it, those values are just called arguments.

In between the curly braces, you can see the lines of code that constitute the body of the function. This is the code that is executed each time you call rt.

Notice that there is no explicit “return” keyword to state which value should be returned from the function. In R, the last value that is calculated in the function is automatically returned. In the case of rt, if the ncp argument is omitted, some C code is called to generate the random numbers, and those are returned. Otherwise, the function calls the rnorm, rchisq, and sqrt functions to generate the numbers, and those are returned.

To create our own functions, we just assign them as we would any other variable. As an example, let’s create a function to calculate the length of the hypotenuse of a right-angled triangle (for simplicity, we’ll use the obvious algorithm; for real-world code, this doesn’t work well with very big and very small numbers, so you shouldn’t calculate hypotenuses this way):

hypotenuse <- function(x, y)
{
  sqrt(x ^ 2 + y ^ 2)
}

Here, hypotenuse is the name of the function we are creating, x and y are its (formal) arguments, and the contents of the braces are the function body.

Actually, since our function body is only one line of code, we can omit the braces:

hypotenuse <- function(x, y) sqrt(x ^ 2 + y ^ 2)   #same as before

We can now call this function as we would any other:

hypotenuse(3, 4)

## [1] 5

hypotenuse(y = 24, x = 7)

## [1] 25

When we call a function, if we don’t name the arguments, then R will match them based on position. In the case of hypotenuse(3, 4), 3 comes first so it is mapped to x, and 4 comes second so it is mapped to y.

If we want to change the order that we pass the arguments, or omit some of them, then we can pass named arguments. In the case of hypotenuse(y = 24, x = 7), although we pass the variables in the “wrong” order, R still correctly determines which variable should be mapped to x, and which to y.

It doesn’t make much sense for a hypotenuse-calculating function, but if we wanted, we could provide default values for x and y. In this new version, if we don’t pass anything to the function, x takes the value 5 and y takes the value 12:

hypotenuse <- function(x = 5, y = 12)
{
  sqrt(x ^ 2 + y ^ 2)
}
hypotenuse() #equivalent to hypotenuse(5, 12)

## [1] 13

We’ve already seen the formals function for retrieving the arguments of a function as a (pair)list. The args function does the same thing in a more human-readable, but less programming-friendly, way. formalArgs returns a character vector of the names of the arguments:

formals(hypotenuse)

## $x
## [1] 5
##
## $y
## [1] 12

args(hypotenuse)

## function (x = 5, y = 12)
## NULL

formalArgs(hypotenuse)

## [1] "x" "y"

The body of a function is retrieved using the body function. This isn’t often very useful on its own, but we may sometimes want to examine it as text—to find functions that call another function, for example. We can use deparse to achieve this:

(body_of_hypotenuse <- body(hypotenuse))

## {
##     sqrt(x^2 + y^2)
## }

deparse(body_of_hypotenuse)

## [1] "{"                   "    sqrt(x^2 + y^2)" "}"

The default values given to formal arguments of functions can be more than just constant values—we can pass any R code into them, and even use other formal arguments. The following function, normalize, scales a vector. The arguments m and s are, by default, the mean and standard deviation of the first argument, so that the returned vector will have mean 0 and standard deviation 1:

normalize <- function(x, m = mean(x), s = sd(x))
{
  (x - m) / s
}
normalized <- normalize(c(1, 3, 6, 10, 15))
mean(normalized)        #almost 0!

## [1] -5.573e-18

sd(normalized)

## [1] 1

There is a little problem with our normalize function, though, which we can see if some of the elements of x are missing:

normalize(c(1, 3, 6, 10, NA))

## [1] NA NA NA NA NA

If any elements of a vector are missing, then by default, mean and sd will both return NA. Consequently, our normalize function returns NA values everywhere. It might be preferable to have the option of only returning NA values where the input was NA. Both mean and sd have an argument, na.rm, that lets us remove missing values before any calculations occur. To avoid all the NA values, we could include such an argument in normalize:

normalize <- function(x, m = mean(x, na.rm = na.rm),
  s = sd(x, na.rm = na.rm), na.rm = FALSE)
{
  (x - m) / s
}
normalize(c(1, 3, 6, 10, NA))

## [1] NA NA NA NA NA

normalize(c(1, 3, 6, 10, NA), na.rm = TRUE)

## [1] -1.0215 -0.5108  0.2554  1.2769      NA

This works, but the syntax is a little clunky. To save us having to explicitly type the names of arguments that aren’t actually used by the function (na.rm is only being passed to mean and sd), R has a special argument, ..., that contains all the arguments that aren’t matched by position or name:

normalize <- function(x, m = mean(x, ...), s = sd(x, ...), ...)
{
  (x - m) / s
}
normalize(c(1, 3, 6, 10, NA))

## [1] NA NA NA NA NA

normalize(c(1, 3, 6, 10, NA), na.rm = TRUE)

## [1] -1.0215 -0.5108  0.2554  1.2769      NA

Now in the call normalize(c(1, 3, 6, 10, NA), na.rm = TRUE), the argument na.rm does not match any of the formal arguments of normalize, since it isn’t x or m or s. That means that it gets stored in the ... argument of normalize. When we evaluate m, the expression mean(x, ...) is now mean(x, na.rm = TRUE).

If this isn’t clear right now, don’t worry. How this works is an advanced topic, and most of the time we don’t need to worry about it. For now, you just need to know that ... can be used to pass arguments to subfunctions.

Functions can be used just like other variable types, so we can pass them as arguments to other functions, and return them from functions. One common example of a function that takes another function as an argument is do.call. This function provides an alternative syntax for calling other functions, letting us pass the arguments as a list, rather than one at a time:

do.call(hypotenuse, list(x = 3, y = 4)) #same as hypotenuse(3, 4)

## [1] 5

Perhaps the most common use case for do.call is with rbind. You can use these two functions together to concatenate several data frames or matrices together at once:

dfr1 <- data.frame(x = 1:5, y = rt(5, 1))
dfr2 <- data.frame(x = 6:10, y = rf(5, 1, 1))
dfr3 <- data.frame(x = 11:15, y = rbeta(5, 1, 1))
do.call(rbind, list(dfr1, dfr2, dfr3)) #same as rbind(dfr1, dfr2, dfr3)

##     x        y
## 1   1  1.10440
## 2   2  0.87931
## 3   3 -1.18288
## 4   4 -1.04847
## 5   5  0.90335
## 6   6  0.27186
## 7   7  2.49953
## 8   8  0.89534
## 9   9  4.21537
## 10 10  0.07751
## 11 11  0.31153
## 12 12  0.29114
## 13 13  0.01079
## 14 14  0.97188
## 15 15  0.53498

It is worth spending some time getting comfortable with this idea. In Chapter 9, we’re going to make a lot of use of passing functions to other functions with apply and its derivatives.

When using functions as arguments, it isn’t necessary to assign them first. In the same way that we could simplify this:

menage <- c(1, 0, 0, 1, 2, 13, 80) #See http://oeis.org/A000179
mean(menage)

## [1] 13.86

to:

mean(c(1, 0, 0, 1, 2, 13, 80))

## [1] 13.86

we can also pass functions anonymously:

x_plus_y <- function(x, y) x + y
do.call(x_plus_y, list(1:5, 5:1))

## [1] 6 6 6 6 6

#is the same as
do.call(function(x, y) x + y, list(1:5, 5:1))

## [1] 6 6 6 6 6

Functions that return functions are rarer, but no less valid for it. The ecdf function returns the empirical cumulative distribution function of a vector, as seen in Figure 6-1:

(emp_cum_dist_fn <- ecdf(rnorm(50)))

## Empirical CDF
## Call: ecdf(rnorm(50))
##  x[1:50] = -2.2, -2.1,  -2,  ..., 1.9, 2.6

is.function(emp_cum_dist_fn)

## [1] TRUE

plot(emp_cum_dist_fn)

Figure 6-1. An empirical cumulative distribution function

A variable’s scope is the set of places from which you can see the variable. For example, when you define a variable inside a function, the rest of the statements in that function will have access to that variable. In R (but not S), subfunctions will also have access to that variable. In this next example, the function f takes a variable x and passes it to the function g. f also defines a variable y, which is within the scope of g, since g is a subfunction of f. So, even though y isn’t defined inside g, the example works:

f <- function(x)
{
  y <- 1
  g <- function(x)
  {
    (x + y) / 2 #y is used, but is not a formal argument of g
  }
  g(x)
}
f(sqrt(5))      #It works! y is magically found in the environment of f

## [1] 1.618

If we modify the example to define g outside of f, so it is not a subfunction of f, the example will throw an error, since R cannot find y:

f <- function(x)
{
  y <- 1
  g(x)
}
g <- function(x)
{
  (x + y) / 2
}
f(sqrt(5))

##  January February    March    April      May
##   0.6494   1.4838   0.9665   0.4527   0.7752

In the section Environments, we saw that the get and exists functions look for variables in parent environments as well as the current one. Variable scope works in exactly the same way: R will try to find variables in the current environment, and if it doesn’t find them it will look in the parent environment, and then that environment’s parent, and so on until it reaches the global environment. Variables defined in the global environment can be seen from anywhere, which is why they are called global variables.

In our first example, the environment belonging to f is the parent environment of the environment belonging to g, which is why y can be found. In the second example, the parent environment of g is the global environment, which doesn’t contain a variable y, which is why an error is thrown.

This system of scoping where variables can be found in parent environments is often useful, but also brings the potential for mischief and awful, unmaintainable code. Consider the following function, h:

h <- function(x)
{
  x * y
}

It looks like it shouldn’t work, since it accepts a single argument, x, but uses two arguments, x and y, in its body. Let’s try it, with a clean user workspace:

h(9)

##  January February    March    April      May
##   -8.436    6.583   -2.727  -11.976   -6.171

So far, our intuition holds. y is not defined, so the function throws an error. Now look at what happens if we define y in the user workspace:

y <- 16
h(9)

## [1] 144

When R fails to find a variable named y in the environment belonging to h, it looks in h’s parent—the user workspace (a.k.a. global environment), where y is defined—and the product is correctly calculated.

Global variables should be used sparingly, since they make it very easy to write appalling code. In this modified function, h2, y is randomly locally defined half the time. With y defined in the user workspace, when we evaluate it y will be randomly local or global!

h2 <- function(x)
{
  if(runif(1) > 0.5) y <- 12
  x * y
}

Let’s use replicate to run the code several times to see the result:

replicate(10, h2(9))

##  [1] 144 144 144 108 144 108 108 144 108 108

When the uniform random number (between 0 and 1) generated by runif is greater than 0.5, a local variable y is assigned the value 12. Otherwise, the global value of 16 is used.

As I’m sure you’ve noticed, it is very easy to create obscure bugs in code by doing things like this. Usually it is better to explicitly pass all the variables that we need into a function.