Skip to content

Functions and conditional statements

Learning objectives

  • Identify parts of a function
  • Write custom functions
  • Use conditional statements to control outputs
  • Refine functions to improve efficiency
  • Access function source code from other packages

Anatomy of functions

It's not a secret that for as long you've used R, you've been using functions. Here, we're going to make them! Functions are, like all things in R, a type of object (implications discussed in a later episode). They are the "verbs" of the R programming language where, instead of storing data, functions store code that defines how it behaves and operates on other objects.

In an R script, functions are created and assigned using the following syntax:

f <- function(<arguments>) {

  <do_something>

}

Lets break down the anatomy of the pseudocode:

  • f is the assigned name of the function.
  • function() initialises a function call, telling R that
    • We are writing a function
    • Anything inside round brackets are arguments for the function
    • Anything after the closing round bracket and on the same line is the code that defines the function
  • { <do_something> } is the expression that defines how the function works

Let's write a function to give some context as to what each part does.

\[ f(x) = \sqrt{2(x + 1)} \]

The above is an equation for a function called \(f\) that takes a value \(x\) as input and outputs a value based on \(x\). We can translate that into R code:

code

1
2
3
f1 <- function(x) {
  sqrt(2 * (x + 1))
}

Here, we've "wrapped" an arithmetic operation in a function.

Alternative syntax

We can also use a shorthand for the function() call by replacing function with a \:

code

f1 <- \(x) sqrt(2 * (x + 1))

This syntax is useful for on-the-fly functions that can be applied using map() functions.

Notice also that the function expression is on the same line as the function call. This syntax is acceptable for short expressions but not encouraged for longer ones.

Infix functions

Most objects we think of as functions are prefix functions, where arguments are placed after the function call. However, there are also infix functions, where arguments are placed on either side of the function call. For example:

code

x <- 2 + 3

Both <- and + are infix functions for assignment and addition, respectively. All basic arithmetic (e.g., +, -, *, and /), assignment (including =), logical (e.g., |, &, &&, etc.), and relational (e.g., >, !=, ==, etc.) operators are primitive (i.e., foundational part of the language) infix functions.

The pipe operator %>% introduced by the magrittr package is a user-generated infix function. As this is a custom infix function, the % that surrounds the operator is necessary. It does not have to be a symbol or operator and can be represented by anything really, (e.g., %in%) even emojis!

We can define our own infix function:

code

1
2
3
4
5
6
`%<I>%` <- function(lhs, rhs) {
  # Get the minima from the left hand side and the maxima from the right hand side
  c(min(lhs), max(rhs))
}

rnorm(10) %<I>% rnorm(200)

Control flow 1: Conditional execution

We will take a little detour from functions to learn about control flow. These are constructs that code for decision making (this episode) and repetition (next episode). Decision making/conditional execution helps to direct the "flow" of operations using Boolean logic (TRUE/FALSE). In R, conditional execution is coded using if-else statments:

1
2
3
4
5
6
7
8
9
if (<predicate>) {

  <expression_1>

} else {

  <expression_2>

}

The pseudocode above means if the <predicate> is (or returns a) TRUE, execute <expression_1>, otherwise execute <expression_2>.

For example:

code

1
2
3
4
5
if (x < 10) {
  print("x is less than 10")
} else {
  print("x is more than 10")
}

Predicate and predicate functions

A predicate is a statement or function that evaluation that returns one TRUE or FALSE value. For example, this is a predicate that evaluates if x is larger than 5:

x > 5

If we set x <- 10, the code will return a TRUE.

R has many built-in predicate functions (i.e., functions that evaluate a statement or object and returns a TRUE/FALSE). They are often prefixed using is.<something>(). For example:

code

x <- 9
tim <- "Tim is eating"
y <- 1:6

# Predicate function
is.character(x)
is.character(tim)
is.character(y)
is.numeric(y)

Notice that regardless of the length of the vector, the predicate functions return only one Boolean variable. This is important for conditional statements as it can only evaluate predicates and not a vector of Boolean variables.

The else statement is optional (and often unnecessary). The expression will not be evaluated if the predicate returns FALSE.

code

1
2
3
4
5
6
7
8
9
# Predicate returns TRUE, expression is evaluated
if (x < 10) {
  print("x is less than 10")
}

# Predicate returns FALSE, expression is not evaluated
if (x > 10) {
  print("x is greater than 10")
}

For complex predicate evaluations (e.g., intermediate thresholds) or evaluations that lead to different paths (e.g., hierarchical decision making), we can chain and/or nest multiple statments together:

if (<predicate_1>) {

  <expression_1>

} else if (<predicate_2>) { # Chained statements: an if-else after an if-else

  if (<predicate_2a>) { # Nested statements: an if-else inside and if-else

    <expression_2a>

  } else {

    <expression_2b>

  }

} else {

  <expression_3>

}

Vectorised conditional execution

The above conditional statements only operate when there is only one Boolean value. What if we wanted to generate different outcomes depending on the elements in a vector? In this case, we can use these two functions depending on how many outcomes we need:

If we needed one outcome each for TRUE or FALSE, we can use if_else(), a vectorised if-else statement that evaluates a Boolean vector (or a condition that returns one) and produces outcomes depending on whether each element is a TRUE or FALSE. For example:

code

1
2
3
4
5
6
7
fridge <- c("banana", "apple", "kale", "durian", "pear", "spinach")

if_else(
  fridge %in% fruit,
  str_c(str_to_upper(fridge), " is a fruit"),
  str_c(fridge, " is not a fruit")
)

Base R equivalent: ifelse()

Base R also has a similar implementation: ifelse(). However, dplyr::if_else() strictly preserves the data type (e.g., the output will not be coerced to a character if the input is a factor) and can evaluate missing values; base::ifelse() does not enforce data type preservation and returns missing values when it encounters them.

If we required different outcomes based on multiple, hierarchical conditions, we can use case_when(), a general vectorised if-else analogous to chained if-else statements. For example:

code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
fridge <- c("cheese", fridge, "milk")
dairy <- c("cheese", "yoghurt", "butter", "milk")
vegetable <- c("beetroot", "beans", "kale")

case_when(
  fridge %in% fruit ~ str_c(fridge, " is a fruit"),
  fridge %in% vegetable ~ str_c(fridge, " is a vegetable"),
  fridge %in% dairy ~ str_c(fridge, " is a dairy product"),
  .default = str_glue("You don't even like {fridge}!")
)

Alpha diversity

Time to apply all of the above to write a function! For this section, we will be writing a function that calculates the \(\alpha\)-diversity (a measure local diversity) of a sample. The simplest measure is richness (often denoted \(q\) or \(S\)), which represents the number of species/taxonomic units. While richness is simple and intuitive to understand, it is heavily biased by rare members of the sample community. This is always the case with microbiome data. Other \(\alpha\)-diversity indices were constructed to account for the uneven distribution of taxonomic units by incorporating relative abundance information. Popular options are:

Shannon's index, \(H\)

\[ H = -\sum_{i=1}^{q} p_{i} \log p_{i} \]

Where \(p_{i}\) is the proportion (or relative frequency) of the \(i^{th}\) organism relative to the sample total.

Simpson's diversity, \(D\)

\[ D = \frac{1}{\lambda} \]

Where the \(\lambda\) is Simpson's concentration index, defined as

\[ \lambda = \sum_{i=1}^{q} p_{i}^2 \]

For the mathematically inclined, you might observe that these indices:

  • Gives more weight to abundant taxa (\(D\) moreso than \(H\))
  • Reach their maxima when all taxa are exactly evenly distributed in the sample
  • Reach their minima when one taxa that dominate
  • Are not on the same scale (\(H\) is on a log scale, \(q\) and \(D\) are on the "number of taxonomic units" scale)

For now, it is enough to know that they are mathematical formulas that we are going to translate into code.

Starting small: richness

We will start by writing a function to calculate richness (the number of taxonomic units in the sample). Applied in R, this means the number of non-zero elements in a numeric vector.

code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# Test vector
test_vector <- asv$AS1A3

# Richness
calc_alpha_diversity <- function(x) {
  sum(x > 0)
}

# Test function
calc_alpha_diversity(test_vector)

Adding arguments: Shannon's and Simpsons's diversity

Next, lets add the other indices to be part of the function. We will use conditional execution to control output depending on the index is requested.

code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
calc_alpha_diversity <- function(x, index) {

  # Richness
  q <- sum(x > 0)

  # Shannon's index
  p <- x / sum(x)
  H <- - sum(p * log(p))

  # Simpson's diversity
  p <- x / sum(x)
  D <- 1 / sum(p ^ 2)

  if (index == "richness") {
    return(q)
  } else if (index == "shannon") {
    return(H)
  } else if (index == "simpson") {
    return(D)
  }

}

# Test that it works
calc_alpha_diversity(test_vector, index = "richness")
calc_alpha_diversity(test_vector, index = "shannon")
calc_alpha_diversity(test_vector, index = "simpson")

Here, we're calculating all possible outcomes, then returning values based on the index requested.

return()

The return() function returns/outputs the value specified in it and stops evaluation of any code below it. In calc_alpha_diversity(), if the index requested is richness, and subsequent code to calculate the other indices are not evaluated. In the case where a function code has no return(), the function will output the last evaluation.

The function we wrote works for \(q\) and \(D\), but not for \(H\). Lets revise that.

Troubleshooting and refining code 🔧

Lets troubleshoot our function first. Why did the code not work for Shannon's index?

A little math

log(0) is undefined. However, in R, it is given the value -Inf for bases > 0 and Inf for bases \< 0. If any p is 0, the operation 0 * log(0) will return NaN (Not a Number). Therefore, the summation of that will return NA.

We will also consider other areas for refinement:

Relative frequency p is not required for richness calculations and is shared for calculations of Shannon's and Simpson's indices.

Assigning of \(H\) and \(D\) helps us, the human, to differentiate the indices in mathematical notation. However, R does not care which name is used. The decision to return either index is already managed by conditional statements.

Chained statements are not necessary here. Conditional execution (and return()) should occur earlier to prevent unnecessary code evaluation.

Some people may have a reason to use another base for the log() function used in calculating Shannon's index. While the natural logarithm (log in base \(e\)) is often the default choice, there should be an option to change that if desired. To set a default value in a custom function, add the desired value to the argument within the function call (see below).

With that, lets revise the function:

code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
calc_alpha_diversity <- function(
  x, index, 
  log.base = exp(1) # Provide natural logarithm as default
) {

  # Richness
  if (index == "richness") {
    return(sum(x > 0))
  }

  # Relative frequency
  p <- x / sum(x)

  # Shannon's index
  if (index == "shannon") {
    H <- - sum(p * log(p, base = log.base), # Provide option for changing log base
               na.rm = TRUE)                # Removes NA prior to summation
  }

  if (index == "simpson") {
    H <- 1 / sum(p ^ 2)
  }

  H

}

# Test the function
calc_alpha_diversity(test_vector, index = "richness")
calc_alpha_diversity(test_vector, index = "shannon")
calc_alpha_diversity(test_vector, index = "simpson")

It works! Note that without a return, the function simply outputs the last evaluted code (in this case, this is the value H).

Defending our function 🛡

What would happen if we used our functions with non-ideal inputs? For this, we will take hints from the concepts in defensive programming. Defensive programming is more about design choices than cybersecurity, where we approach designing code that:

  • Has predictable outputs, regardless of inputs
  • "Fails fast", thus saving time and computational resources
  • Reads well, so others can audit the code and understand what each line is doing

Can you identify some scenarios where our function might not do so well?

Pitfalls and food for thought

What does NA mean?

Lets consider how these tables were constructed. They were originally sequences that were processed by a software and then reads were counted for each taxonomic unit. If a taxonomic unit was not found in a particular sample read file, the software would have returned a 0. Therefore, NA could not have arisen from processing this particular batch.

However, we can imagine a scenario where NA could be introduced from a joining step during data cleaning or pre-processing. This can happen when processing samples by batches.

Critically, regardless of the reasons for missingness, is a missing count functionally equivalent to 0?

What is the effect of NA?

We used sum() quite often in the code. Furthermore, we even added na.rm = TRUE in the expression used to calculate Shannon's index to account for the NA introduced when evaluating \(0 \log 0\). Would we handle NA in all parts where sum() is used in this way?

Do negative numbers make sense?

Remember that this is a table of counts. What do negative numbers mean in that context? Does it make sense?

How does the function handle negative numbers?

\(q\) would ignore negative numbers. If the negative counts were an oversight or accidental, it would result in an undercount of richness.

\(H\) would not produce a numeric output (the log of a negative number is undefined).

\(D\) would be an underestimate (the sum would be an undercount of the whole).

Does non-numeric data make sense for this function?

These indices are meant for count data (i.e., numeric data). Suppose that human error led to the mis-typing of inputs, would we want to coerce the data into numeric type using as.numeric()? Are we guessing the intentions and errors of the user? Coercing characters that are not numbers will cause as.numeric() to return NA.

Is it fit for purpose?

We established that the function is meant to provide information only on numeric data. Albeit being a small function that runs quickly, we can imagine a scenario where we might either (1) write a function that has many complex steps or (2) run a substantially larger data set through the diversity function. In either scenario, would we want to spend time and computational resources on non-target inputs?

Let's test our function through these non-ideal scenarios. Create vectors with non-ideal properties and try them on the current version of the function.

code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# Randomly introduce NA
test_missing <- test_vector
test_missing[sample.int(length(test_missing), 5)] <- NA

# Randomly introduce negative numbers
test_negative <- ifelse(
  test_vector > 0 & runif(length(test_vector)) > 0.95,
  -test_vector, test_vector
)

# Coerce to character
test_character <- as.character(test_vector)

Random replacements

The code above uses two ways to generate replacements at random. In line 3, a random sampler sample.int() was used to indicate 5 random indices to which we wanted to inject NA. In lines 6-9, a conditional statement was constructed to replace positive, non-zero values, picked based on a numeric (double) vector picked from a uniform distribution, with negative equivalents.

With those pitfalls in mind, we will need to make some design decisions and implement them in code. Here are my choices:

Functionally, they are the same as not observing a taxa in a sample (equivalent to 0). Warn the user when this happens but continue with evaluating input.

These data types are not intended for use with diversity indices here. Lets not spend computing time and resources on them. To achieve this, we will "error out" of the function where if negative numbers and non-numeric data are encountered, the function will stop further computation and produce an error message.

Feel free to implement your own design choices by modifying the code below!

code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
calc_alpha_diversity <- function(
  x, index, 
  log.base = exp(1) # Natural logarithm as default
) {

  # Warn for missing values
  if (anyNA(x)) {
    warning("x contains NAs which are removed prior to calculations.")
    x <- x[!is.na(x)]
  }

  # Error on non-numeric and negative data
  if (!is.numeric(x) | any(x < 0)) {
    stop("x must be positive numeric values!")
  }

  # Richness
  if (index == "richness") {
    return(sum(x > 0))
  }

  # Relative frequency
  p <- x / sum(x)

  # Shannon's index
  if (index == "shannon") {
    H <- - sum(p * log(p, base = log.base), # Provide option for changing log base
               na.rm = TRUE)                # Removes NA from 0 * log(0)
  }

  if (index == "simpson") {
    H <- 1 / sum(p ^ 2)
  }

  H

}

Try the new function on the test vectors:

code

1
2
3
calc_alpha_diversity(test_missing, index = "shannon")
calc_alpha_diversity(test_negative, index = "simpson")
calc_alpha_diversity(test_character, index = "richness")

Improve function design

Read other functions

One of the best ways to improve how we write functions is to read other people's functions! Not only do we get to see how others design their functions, we can also learn new and useful functions along the way!

We can access the source code of a function using getAnywhere():

code

getAnywhere(lm)

Some functions may not have directly accessible source code. For example:

code

getAnywhere(anova)

Notice that the body of the function consists only of UseMethod("anova"). This is because anova() can take many different kinds of inputs and run different code depending on the input (i.e., it has different methods() for the same function). We can see what those methods are and inspect the source code for any of them:

code

# What methods are available for anova?
methods(anova)

# Inspect anova.lm
getAnywhere(anova.lm)

In some cases, a function could have the same name but are implemented differently because they are part of a different package. For example:

code

getAnywhere(filter)

To access the filter() source code from the stats package, we need pick one by its index:

code

# Get source code for filter in stats package
getAnywhere(filter)[2]

return() early

Use return() early in the source code especially if other arguments will lead to more time consuming evaluations. This prevents wasteful use of time and resources to compute variables that will not even be used as outputs!

switch() evaluations depending on argument input

swtich() is a useful function that evaluates different expressions depending on the argument inputs. It can also be given a default output when no arguments match those specified which can be coded to output a warning or message.

code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
calc_alpha_diversity <- function(x, index, log.base = exp(1)) {

  # Warn for missing values
  if (anyNA(x)) {
    warning("x contains NAs which are removed prior to calculations.")
    x <- x[!is.na(x)]
  }

  # Error on non-numeric and negative data
  if (!is.numeric(x) | any(x < 0)) {
    stop("x must be positive numeric values!")
  }

  # Richness
  if (index == "richness") {
    return(sum(x > 0))
  }

  # Relative frequency
  p <- x / sum(x)

  # Indices
  switch(
    index, # Value to be matched against in the list below
    shannon = - sum(p * log(p, base = log.base), na.rm = TRUE),
    simpson = 1 / sum(p ^ 2),
    warning("Unknown index. Should be one of 'richness', 'shannon', or 'simpson'.") # Default output if unmatched
  )

}

Think of the users

Who is the target user for your functions?

If it's just you, do whatever you want! However, the fact that you are assigning code to functions implies that you probably want to reuse it for other things/projects or even share them with your colleagues/the world!

If you intend to share your functions with other people/computers/robots, it is important to think about dependencies. Other people will have different set-ups. Thus, it is advisable to write functions using other functions in base R. If you must use other packages, make sure to advise and install early in the script.

While we write functions to be as predictable, elegant, and concise as possible, there are situations where we want to add flexibility. Keep in mind that flexible functions often involve more code. To achieve this, we can write functions that call other custom functions. This is called modularisation. Modular code also helps with troubleshooting, readability, and benchmarking as pieces of it can be tested. This is especially important for production-level code.

What stage are you at?

Assuming you are writing functions for your analyses, the syntax of your functions are likely geared towards your specific needs at a particular given stage. At the exploratory analysis stage, functions are written with flexbility in mind, perhaps even a little messy or less efficient that it could be.

As you progress to production level code and/or need to use complex processes (permutation tests, sub-sampling, Bayesian methods, etc.), you're likely going to need to think about:

✅ Code safety Predictable and non-redundant inputs and outpus, and preventing ad inifinitum runs

✅ Efficiency Use benchmarked algorithms, reduce/remove side effects, parallel processing, reduce overheads

✅ Troubleshooting Informative errors and warnings, modularised functions

✅ Readability Commented code, verbose syntax, appropriate assignments (as opposed to long piped processes)

Break your own code!

The best way to learn from your own functions? Break them! How would it handle nonsensical inputs? What are its limits? How much data is too much data? Is it taking seconds, minutes, days to run? We can learn so much about how the function (and R in general) works when we stress test them.