Functions (Completed Notebook)

R for Data Science

Set-up

library(tidyverse)

Acknowledgements

This notebook is based on Chapter 19 of R for Data Science.

1 What is a function?

You’ve already used tons of functions.

Take for example sum(). It sums the elements of a numerical vector:

# vector with elements 1 and 2
x <- c(1,2) 
# add them up
sum(x)

[1] 3

1.1 Vending machines

In general a function is an object that converts an input into an output.

A vending machine is a function:

• input: money
• output: tasty snack

When you write functions think of yourself as a vending-machine designer. Why? Because vending machines are easy to use. The input is obvious. The output is obvious. Nobody has to read a manual to use a vending machine.

So when you design a function, ask yourself:

• Does it have an easy-to-understand name?
• Does it have easy-to-understand inputs (arguments)?
• Does it have an easy-to-understand output?

1.2 Who cares?

In some cases these desiderata must be satisfied so other people can use your code without fuss.

But in most cases you should think about your future self. They review code you wrote a few weeks ago. myfunc(x). What the hell is that?!

1.3 Our focus

We will focus on writing functions with different inputs:

• vectors of length 1 (e.g., a number)
• vectors of length \(n\) (e.g., many numbers)

Our goal is to build up to iteration: applying a function many times to different objects (e.g., to many columns of a tibble/data frame). Efficiently running the same calculation on different data is where programming really becomes important.

1.4 When should you write a function?

The book says it best:

You should consider writing a function whenever you’ve copied and pasted a block of code more than twice (i.e. you now have three copies of the same code).

Here are some other thoughts:

• Writing functions helps you think more about design (“How should this code be used?”).
• Writing functions is a great way to learn stats/machine learning. By coding up your own estimators you can really understand what goes on under the hood of canned estimators. This will make you a better scientist.
  • that said you should still use canned estimators in real work because they typically are optimized (e.g., written in C)

2 The function template

Here is a template for writing a function:

informative_name <- function(x) {
  output <- # do stuff and create an output
  return(output) # return the output
}

The function informative_name() takes an input x and returns output. If output is a specific thing (e.g., the sum of squares of a vector) you could give it a specific name (e.g., sum_of_squares).

As always in programming there are multiple ways of doing the same thing.

For instance, some people prefer a non-explicit return call:

informative_name2 <- function(x) {
  output <- # do stuff and create an output
  output # return the output
}

This function does the same thing: it returns output.

The book makes a few arguments (no pun intended!) as to why you should or shouldn’t explicitly call return() (see section 19.6.1). Mostly these are edge cases and you won’t run into them. At the end of the day you navigate it like any other convention. Do what everybody else does. That way your code can enjoy network effects.

In the case of return() it does not really matter. The advantage of return() is that it is explicit. Looking at informative_name() it is crystal clear what is returned. But at the end of the day this is a situation where you can pick one or the other. What matters is sticking with a choice. Consistency is key.

3 Arguments

3.1 Arguments of length 1

Everything in R is a vector.

For instance the number 2 is a vector:

class(2)

[1] "numeric"

A numeric vector of length 1.

Consider this function:

add_two <- function(x=2) { 
  output <- x + 2
  return(output)
}

add_two() does what it says on the tin: takes a number and adds two to it.

Press the “play” button to compile it. You will see it appear in the Environment tab. Once compiled a function can be used.

For instance:

# input: 3
# output: 3 + 2 = 5
add_two(x = 3)

[1] 5

Notice how we wrote function(x=2). That means the argument x is optional. If you don’t feed a number to add_two(), it will feed itself with its default value of 2:

add_two()

[1] 4

3.2 Checkpoint

Write a function called square_sum that takes a numerical input x (defaulting to 3) and

• squares it, then;
• adds the square to x to create an object output, then;
• returns output:

square_sum <- function(x=3){
  square_x = x^2
  output <- x + square_x
  return(output)
}

Test it out by running:

• square_sum(x=3)
• square_sum(x=100)

square_sum(x = 3)

[1] 12

3.3 Local variables

Inside square_sum you defined an object output. But if you call it you’ll get an error:

output

Error: object 'output' not found

This is because output is a local variable. It only “exists” (to R) inside square_sum. The big idea here is that functions are environments. What happens in Vegas stays in Vegas – unless you return() it!

3.4 Multiple “single” numerical arguments

Functions of course can take many arguments. Just specify them inside function():

informative_name <- function(x, y, z) {
  output <- # do stuff with x, y and z to create an output
  return(output) # return the output
}

3.5 Checkpoint

Write a function sum_squares that takes two arguments x and y and:

• squares x and y, then;
• adds the squares to create output, then;
• returns output

sum_squares <- function(x, y){
  square_x = x^2
  square_y = y^2
  output = square_x + square_y
  return(output)
}

Test it out for x=10 and y=100:

sum_squares(2, 3)

[1] 13

class(sum_squares)

[1] "function"

3.6 Arguments of length \(n\)

sum() can add a single number:

sum(2)

[1] 2

which of course is not useful at all.

What is useful is sum() adding up many numbers:

set.seed(123) # so we get the same random numbers
# create a vector of 1000 random standard normal numbers (mean 0, sd 1)
random_normal_numbers <- rnorm(n = 1000, mean = 0, sd = 1)
# add them up
sum(random_normal_numbers)

[1] 16.12787

To sum() the argument x is always a numerical vector of any length.

3.7 Checkpoint

Re-writesquare_sum so that it takes in a vector, squares each value, adds them up, and returns:

square_sum <- function(x){
  square_x = x^2
  output = sum(square_x)
  return(output)
}

Test it out on random_normal_numbers:

square_sum(x = random_normal_numbers)

[1] 982.7356

4 Functions to do math

Writing functions is useful when you want to do math.

Most math is usually just a bunch of small, easy steps (addition, subtraction, etc.) chained together.

You might find yourself in this position if you need to code up a stats estimator.

4.1 Checkpoint

Write a function called my_sd() that calculates the standard deviation of a vector. The formula is:

\[ s = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (x_i - \bar{x})^2} \]

Think about the steps when coding it up:

1. calculate the deviations \(x_i - \bar{x}\) where \(\bar{x}\) is the sample mean
2. square the deviations
3. add up the squared deviations
4. multiply by \(\frac{1}{n-1}\) (use length()) to get the length of x, i.e. \(n\)
5. take the square root

my_sd = function(x){
  # step 1
  deviations = x - mean(x)
  # step 2
  squared_deviations = deviations^2
  # step 3
  sum_squared_deviations = sum(squared_deviations)
  # step 4
  n = length(x)
  variance = (1/(n-1))*sum_squared_deviations
  # step 5
  sd = sqrt(variance)
  # return!!!!!
  return(sd)
}

Test it out on random_normal_numbers:

my_sd(random_normal_numbers)

[1] 0.991695

Compare it to the canned routine sd():

sd(random_normal_numbers)

[1] 0.991695

Side note: Why is the result close to 1?

Because we drew those numbers from a distribution with a standard deviation of 1. So my_sd() and sd() recover the parameter.

This idea of “recovering” paramaters is the bedrock of stats. The goal is to “recover” the population parameter from a sample.

4.2 Checkpoint

Most often you will use functions on observational data rather than simulated data.

Load the diamonds data set:

data("diamonds")

The standard deviation of price is just sd() on the price vector inside the diamonds tibble:

diamonds %>% 
  summarise(sd(price))

Replicate the results from sd() with my_sd().

diamonds %>% 
  summarise(my_sd(price))

4.3 Checkpoint

Do carats influence the price of diamonds? We can build a model:

\[ \begin{aligned} \text{prices} &= f(\text{carats}) + \epsilon \\ &= \beta_0 + \beta_1(\text{carats}) + \epsilon \end{aligned} \] where \(\beta_1\) captures the average effect of carats on prices. If \(\hat{\beta_1}\) (the estimate) is positive, then more carats = higher prices, on average. The formula for \(\beta_1\) is:

\[ \hat{\beta_1} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2} \] We can estimate the parameters with lm():

lm(formula = price ~ carat, data = diamonds) %>% 
  summary()

Call:
lm(formula = price ~ carat, data = diamonds)

Residuals:
     Min       1Q   Median       3Q      Max 
-18585.3   -804.8    -18.9    537.4  12731.7 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2256.36      13.06  -172.8   <2e-16 ***
carat        7756.43      14.07   551.4   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1549 on 53938 degrees of freedom
Multiple R-squared:  0.8493,    Adjusted R-squared:  0.8493 
F-statistic: 3.041e+05 on 1 and 53938 DF,  p-value: < 2.2e-16

Write a function called my_lm() that takes an x and y argument and calculates the regression coefficient of x on y. Then replicate the results from lm(formula = price ~ carat, data = diamonds). Think about what code you can re-use from my_sd().

my_lm = function(x, y){
  # 1. NUMERATOR
  ## step 1
  x_deviations = x - mean(x)
  y_deviations = y - mean(y)
  ## step 2
  product_deviations = x_deviations * y_deviations
  ## step 3
  numerator = sum(product_deviations)
  # 2. DENOMINATOR
  ## step 1
  x_deviations_squared = x_deviations^2
  ## step 2
  denominator = sum(x_deviations_squared)
  # 3. put it all together
  beta = numerator / denominator
  # RETURN!!
  return(beta)
}

my_lm(x = diamonds$carat, y = diamonds$price)

[1] 7756.426

5 Conditional execution

What if you want your function to do different things depending on the input?

Use an if statement!

This function checks if the argument is equal to two. If it is, add two and return.

add_two_if <- function(x) { 
  if(x == 2){
    output <- x + 2
    return(output)
  }
}

Try it out:

add_two_if(x = 1)

The if(x==2) is the logical test, a test that only returns TRUE or FALSE:

2 == 2

[1] TRUE

2 == 3

[1] FALSE

So when if(x==2) evaluates to TRUE the function proceeds.

5.0.1 Logical tests

Like in dplyr::filter() we can use the following conditionals inside if():

• ==: “is equal to”
• !=: “is not equal to”
• &: “and”
• |: “or”
• >, >=: “greater than”, “less than or equal to”
• <, <=: “less than”, “less than or equal to”

5.0.2 else

But what if we pass an argument other than two?

add_two_if(x = 3)

Nothing happens.

What we need is a way to handle the exception – a way to deal with arguments that are not equal to 2.

This is where else comes in.

The general structure of an if-else statement is:

if(TRUE) { # if the test is true
  # do this thing
} else { # otherwise
  # do this other thing
}

Let’s modify add_two_if to add three if an argument is not two:

add_two_if <- function(x) { 
  if(x > 2 && x < 4){
    output <- x + 2
    # do some function
  } else {
    output <- x - 2
    # do some other function
  }
  return(output)
}

add_two_if(x = 2)

[1] 0

Notice how the return() statement is outside the if...else scope. (Anything in R wrapped in {} is a “scope”).

In theory you can chain together as many if...else statements as you want by adding more logical tests:

add_two_if <- function(x) { 
  if(x == 2){
    output <- x + 2
  } else if(x == 3) {
    output <- x + 3
  } else if(x == 4){
    output <- x + 3
  } else { 
    print("this is a bad function")
  }
  return(output)
}

But you want to avoid this. It’s hard to read and probably not necessary. The chapter has a few suggestions for alternatives. For now we will stick to one logical test.

5.1 Checkpoint

Write a function truth_or_dare that takes a vector and checks whether the third entry is positive. If yes, print “truth”. If no, print “dare”.

Update We changed this to also demo stop().

truth_or_dare = function(x){
  if(length(x) == 1) stop("vector is too short")
  if(x[3] %% 2 == 0) {
    warning("The challenge is to do ", round(rnorm(1, mean = 10, sd =5),0), " dares")
    #print(paste())
  } else {
    print("truth")
  }
}

vec = c(2,3,4,5)
truth_or_dare(2)

Error in truth_or_dare(2) : vector is too short

5.2 ifelse()

The if...else scope only works for vectors of length 1.

Run this chunk:

another_vector <- c(1, NA, 3, NA)
if(is.na(another_vector)) {
  "NA :("
} else{
  "number"
}

the condition has length > 1 and only the first element will be used

[1] "number"

if stopped after the first entry of another_vector.

Use ifelse() to run an if..else scope over an entire vector:

ifelse(test = is.na(another_vector), yes = 0, no = "a")

[1] "a" "0" "a" "0"

6 Appendix

6.1 Dot-dot-dot

You may have noticed that many functions in R have the argument ....

For instance, the arguments to mean() are:

mean(x, ...)

The argument x is mandatory. The argument ... is optional. What is cool about ... is that you can pass any number of arbitrary arguments to it.

For instance, what if you have NA in a vector?

v <- c(1,2,NA)

mean(x) will return NA:

mean(v)

[1] NA

but you can pass na.rm to ... so that mean() only looks at numerical values:

mean(v, na.rm = TRUE)

[1] 1.5

Ditto sd(). Without na.rm:

sd(v)

[1] NA

With na.rm:

sd(v, na.rm = TRUE)

[1] 0.7071068

This underscores the real value of ...: passing optional arguments to other functions that have optional arguments – like mean() and sd().

6.1.1 Checkpoint

Say you have this vector with NAs:

# create a vector of random numbers
na_vector <- rnorm(n = 1000, mean = 0, sd = 1)
# randomly replace some of the values as NA
na_vector[sample(1:1000, 100)] <- NA

Write a function called z_score that calculates the z-score of each element of na_vector. The formula is:

\[ z_i = \frac{x_i - \bar{x}}{s} \]

where \(\bar{x}\) is the sample mean and \(s\) is the sample standard deviation. Make sure your function can account for NAs by including ... and passing that to mean(x, ...) and sd(x, ...).

z_score = function(x, ...){
  z = (x - mean(x, ...)) / sd(x, ...)
  return(z)
}

Test it out on na_vector:

z_score(x = na_vector, na.rm = TRUE) %>% 
  mean(x = ., na.rm = TRUE)

[1] -7.998159e-18

Section 19.5.3 links to other ways to use .... They’re handy when you write functions that take entire data sets as arguments. For more see Chapters 18-21 of the Tidyverse Design Guide.

6.1.2 Aside: randomly sampling a vector

Let’s randomly replace some of the values in random_normal_numbers using indexing and sample().

Recall that [] indexes vectors by position. For instance, the first three elements of random_normal_numbers:

random_normal_numbers[1:3]

[1] -0.5604756 -0.2301775  1.5587083

and sample() will randomly sample any vector. For instance, randomly sample two integers between 1 and 10:

sample(x = 1:10, size = 3)

[1] 1 9 3

For instance I could randomly sample a vector:

# vector of 1, 2, .., 10
x <- 1:10
# the length of the vector
n <- length(x)
# randomly draw three elements from x
x[sample(x = 1:n, size = 3)]

[1] 4 1 9

so we can combine indexing and sample to randomly switch values in a vector to NA:

# get the length of random_normal_numbers (1000)
n <- length(random_normal_numbers)
# from the positional index of 1:n, randomly choose 100 of them, 
## and reassign their corresponding values in random_normal_numbers to NA
random_normal_numbers[sample(1:n, 100)] <- NA