9 Functions

History channel aliens conspiracy guy. Top text: I've got functions inside functions...; Bottom text: ... inside functions.

Intended Learning Outcomes

Write your own custom functions with function()
Be able to test your custom functions

Setup

Open your reprores project
Create a new quarto file called 09-func.qmd
Update the YAML header
Replace the setup chunk with the one below:

```{r}
#| label: setup
#| include: false
library(tidyverse)  # loads purrr for iteration
library(testthat)   # unit testing
library(broom)      # converts test output to tidy tables

set.seed(8675309) # makes sure random numbers are reproducible
```

9.1 Custom functions

In addition to the built-in functions and functions you can access from packages, you can also write your own functions (and eventually even packages!).

9.1.1 Structuring a function

The general structure of a function is as follows:

function_name <- function(my_args) {
  # process the arguments
  # return some value
}

Here is a very simple function. Can you guess what it does?

add1 <- function(my_number) {
  my_number + 1
}

add1(10)

[1] 11

Let’s make a function that reports p-values in APA format (with “p = [rounded value]” when p >= .001 and “p < .001” when p < .001).

First, we have to name the function. You can name it anything, but try not to duplicate existing functions or you will overwrite them. For example, if you call your function rep, then you will need to use base::rep() to access the normal rep function. Let’s call our p-value function report_p and set up the framework of the function.

report_p <- function() {
}

9.1.2 Arguments

We need to add one argument, the p-value you want to report. The names you choose for the arguments are private to that argument, so it is not a problem if they conflict with other variables in your script. You put the arguments in the parentheses of function() in the order you want them to default (just like the built-in functions you’ve used before).

report_p <- function(p) {
}

9.1.3 Argument defaults

You can add a default value to any argument. If that argument is skipped, then the function uses the default argument. It probably doesn’t make sense to run this function without specifying the p-value, but we can add a second argument called digits that defaults to 3, so we can round p-values to any number of digits.

report_p <- function(p, digits = 3) {
}

Now we need to write some code inside the function to process the input arguments and turn them into a returned output. Put the output as the last item in function.

report_p <- function(p, digits = 3) {
  if (p < .001) {
    reported = "p < .001"
  } else {
    roundp <- round(p, digits)
    reported = paste("p =", roundp)
  }
  
  reported
}

You might also see the returned output inside of the return() function. This does the same thing.

report_p <- function(p, digits = 3) {
  if (p < .001) {
    reported = "p < .001"
  } else {
    roundp <- round(p, digits)
    reported = paste("p =", roundp)
  }
  
  return(reported)
}

When you run the code defining your function, it doesn’t output anything, but makes a new object in the Environment tab under Functions. Now you can run the function.

report_p(0.04869)
report_p(0.0000023)

[1] "p = 0.049"
[1] "p < .001"

9.1.4 Scope

What happens in a function stays in a function. You can change the value of a variable passed to a function, but that won’t change the value of the variable outside of the function, even if that variable has the same name as the one in the function.

reported <- "not changed"

# inside this function, reported == "p = 0.002"
report_p(0.0023) 

reported # still "not changed"

[1] "p = 0.002"
[1] "not changed"

9.1.5 Warnings and errors

What happens when you omit the argument for p? Or if you set p to 1.5 or “a”?

You might want to add a more specific warning and stop running the function code if someone enters a value that isn’t a number. You can do this with the stop() function.

If someone enters a number that isn’t possible for a p-value (0-1), you might want to warn them that this is probably not what they intended, but still continue with the function. You can do this with warning().

report_p <- function(p, digits = 3) {
  if (!is.numeric(p)) stop("p must be a number")
  if (p <= 0) warning("p-values are normally greater than 0")
  if (p >= 1) warning("p-values are normally less than 1")
  
  if (p < .001) {
    reported = "p < .001"
  } else {
    roundp <- round(p, digits)
    reported = paste("p =", roundp)
  }
  
  reported
}

report_p()

Error in `report_p()`:
! argument "p" is missing, with no default

report_p("a")

Error in `report_p()`:
! p must be a number

report_p(-2)

Warning in report_p(-2): p-values are normally greater than 0

report_p(2)

Warning in report_p(2): p-values are normally less than 1

[1] "p < .001"
[1] "p = 2"

9.2 Unit Tests

Unit testing is the practice of checking that your functions provide expected output given specific input. You can do unit testing informally, by running your function in the console and deciding if the output makes sense, or more formally, by recording code that compared the actual and expected output.

Let’s take the example of out report_p function. Let’s check that it gives the expected output of “p = 0.05” when you set p to .05. We’ll use the term object to refer to the output of the function and expected to refer to what the function should produce, as this will map onto the terminology used in more formal unit testing. Here, we can use == to check that they are equivalent.

object <- report_p(.05)
expected <- "p = 0.05"
object == expected

[1] TRUE

Now let’s check that you get the expected output for a number that needs to be rounded.

object <- report_p(.123456789)
expected <- "p = 0.123"
object == expected

[1] TRUE

We should check that the digits argument works as expected, too.

object <- report_p(.123456789, digits = 2)
expected <- "p = 0.12"
object == expected

[1] TRUE

You should also check that your function handles incorrect input in the way you want it to. We set the function to give a warning when p is greater than 1 or less than 0, but still return a formatted string.

object <- report_p(1.5)

Warning in report_p(1.5): p-values are normally less than 1

expected <- "p = 1.5"
object == expected

[1] TRUE

9.2.1 Testthat

We can see the warning, but how do we test that it’s what we expect? You can set this up with base R, but at this point it’s easier to move to more formal unit testing with the testthat package. This package provides a lot of function that you can use to test your functions both informally and in package development.

The function expect_equal() is one of the more commonly used functions. The arguments are object and expected (you can probably fuigure out what to put here).

object <- report_p(.05)
expected <- "p = 0.05"

expect_equal(object, expected)

This function does nothing overt when the test passes, but throws an error if the test doesn’t pass, including some explanation of the differences. When you run this in a development environment for creating an R package, it can also provide a summary of the success of all your tests.

object <- report_p(.05)
expected <- "p = 0.5"

expect_equal(object, expected)

Error:
! Expected `object` to equal `expected`.
Differences:
1/1 mismatches
x[1]: "p = 0.05"
y[1]: "p = 0.5"

If you want to check for a specific warning or error, wrap the code in expect_warning() or expect_error().

expect_warning(object <- report_p(1.5))

If your code does produce a warning, now you won’t see anything. If the code doesn’t produce an expected warning, you will get an error report.

expect_warning(object <- report_p(.05))

Error:
! Expected `object <- report_p(0.05)` to produce warnings.

You can check that the warning is what you expect by setting the second argument to the text. This value is expected to be a regular expression, so if it’s just the literal exact text, set fixed = TRUE.

expect_warning(object <- report_p(1.5),
               "p-values are normally less than 1", 
               fixed = TRUE)

9.2.2 Test Driven Development

How should your function handle different kinds of input? What if p is a vector of numbers?

p <- c(0.05, 1/3, .5)
object <- report_p(p)

Error in `if (p <= 0) ...`:
! the condition has length > 1

This gives an error, so we need to fix the function. Test driven development is a coding philosophy where all development of a function is preceeded by writing a failing unit test. Here, we want the report_p() function to return a vector of formatted strings. So we should write a test that would pass if we achieved this goal.

p <- c(0.05, 1/3, .5)
object <- report_p(p)

Error in `if (p <= 0) ...`:
! the condition has length > 1

expected <- c("p = 0.05", "p = 0.333", "p = 0.5")
expect_equal(object, expected)

Error:
! Expected `object` to equal `expected`.
Differences:
Lengths differ: 1 is not 3

Now, we have a very clear goal and way to know when we’ve achieved it. (And if, like me, you tend to go down rabbit holes and get distracted when you code, the failing test will remind you to come back to this.)

Let’s vectorise this function to let it handle input with a length greater than 1. This uses some iteration.

report_p <- function(p, digits = 3) {
  if (length(p) > 1) {
    reported_vec <- map_chr(p, report_p)
    return(reported_vec)
  }
  
  if (!is.numeric(p)) stop("p must be a number")
  if (p <= 0) warning("p-values are normally greater than 0")
  if (p >= 1) warning("p-values are normally less than 1")
  
  if (p < .001) {
    reported = "p < .001"
  } else {
    roundp <- round(p, digits)
    reported = paste("p =", roundp)
  }
  
  reported
}

Now, let’s rerun the test code that failed before.

p <- c(0.05, 1/3, .5)
object <- report_p(p)
expected <- c("p = 0.05", "p = 0.333", "p = 0.5")
expect_equal(object, expected)

9.2.3 Deciding what to test

How do you decide what to test? Do you need to test that you get the correct output for p = 0.05 and p = 0.1 and p = 0.15 and p = 0.20 and p = 0.25, ad infinitum? Of course not! But it’s difficult to give you a specific rule; this requires judgement, experience, and knowledge about what you want the function to do.

In general, I try to cover the following categories in my tests:

Typical, simple input
Input that should raise a warning
Input that should raise an error
NULL, NA or vector of length 0 input is handled as expected
More complicated input
Change argument defaults (do each argument separately)
Combinations of arguments (if they might interact)

Make some more unit tests for report_p(). What happens if p is NULL, NA or c()?

9.3 Iterating your own functions

9.3.1 Build code

First, let’s build up the code that we want to iterate.

9.3.1.1 Simulate and structure data

Create a vector of 20 random numbers drawn from a normal distribution with a mean of 5 and standard deviation of 1 using the rnorm() function and store them in the variable A.

A <- rnorm(20, mean = 5, sd = 1)

A tibble is a type of table or data.frame. The function tibble::tibble() creates a tibble with a column for each argument. Each argument takes the form column_name = data_vector.

Create a table called dat including two vectors: A that is a vector of 20 random normally distributed numbers with a mean of 5 and SD of 1, and B that is a vector of 20 random normally distributed numbers with a mean of 5.5 and SD of 1.

dat <- tibble(
  A = rnorm(20, 5, 1),
  B = rnorm(20, 5.5, 1)
)

9.3.1.2 Statistical test

You can run a Welch two-sample t-test by including the two samples you made as the first two arguments to the function t.test. You can reference one column of a table by its names using the format table_name$column_name

t.test(dat$A, dat$B)


    Welch Two Sample t-test

data:  dat$A and dat$B
t = -1.4715, df = 37.977, p-value = 0.1494
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.1408729  0.1804332
sample estimates:
mean of x mean of y 
 4.829928  5.310148

You can also convert the table to long format using the gather function and specify the t-test using the format dv_column~grouping_column.

longdat <- gather(dat, group, score, A:B)

t.test(score~group, data = longdat)


    Welch Two Sample t-test

data:  score by group
t = -1.4715, df = 37.977, p-value = 0.1494
alternative hypothesis: true difference in means between group A and group B is not equal to 0
95 percent confidence interval:
 -1.1408729  0.1804332
sample estimates:
mean in group A mean in group B 
       4.829928        5.310148

9.3.1.3 Tidy output

You can use the function broom::tidy() to extract the data from a statistical test in a table format. The example below pipes everything together.

tibble(
  A = rnorm(20, 5, 1),
  B = rnorm(20, 5.5, 1)
) |>
  gather(group, score, A:B) |>
  t.test(score~group, data = _) |>
  broom::tidy()

estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
-0.7629601	4.725315	5.488276	-2.625878	0.0130202	32.86047	-1.354193	-0.1717276	Welch Two Sample t-test	two.sided

Note

In the pipeline above, t.test(score~group, data = _) uses the _ notation to change the location of the piped-in data table from it’s default position as the first argument to a different position.

9.3.1.4 Extract important values

Finally, we can extract a single value from this results table using pull().

tibble(
  A = rnorm(20, 5, 1),
  B = rnorm(20, 5.5, 1)
) |>
  gather(group, score, A:B) |>
  t.test(score~group, data = _) |>
  broom::tidy() |>
  pull(p.value)

[1] 0.07260504

9.3.2 Custom function

Next, we can group the code above inside a function.

First, name your function t_sim and wrap the code above in a function with no arguments.

t_sim <- function() {
  tibble(
    A = rnorm(20, 5, 1),
    B = rnorm(20, 5.5, 1)
  ) |>
    gather(group, score, A:B) |>
    t.test(score~group, data = _) |>
    broom::tidy() |>
    pull(p.value) 
}

Run it a few times to see what happens.

t_sim()

[1] 0.004556227

9.3.2.1 Iterate

Let’s run the t_sim function 1000 times, assign the resulting p-values to a vector called reps, and check what proportion of p-values are lower than alpha (e.g., .05). This number is the power for this analysis.

reps <- replicate(1000, t_sim())
alpha <- .05
power <- mean(reps < alpha)
power

[1] 0.346

9.3.2.2 Set seed

You can use the set.seed function before you run a function that uses random numbers to make sure that you get the same random data back each time. You can use any integer you like as the seed.

set.seed(90201)

Warning

Make sure you don’t ever use set.seed() inside of a simulation function, or you will just simulate the exact same data over and over again.

An alignment chart meme: Lawful Good = set.seed(1234), Neutral Good = *picks random number, sets it as seed*, Chaotic Good = *pop culture reference*, Lawful Neutral = set.seed(420), True Neutral = set.seed(1), Chaotic Neutral = str <-'I am awesome.'; set.seed(str_length(str)), Lawful Evil =set.seed(666), Neutral Evil: *does not set a seed*, Chaotic Evil = set.seed(.Random.seed[1]) — Figure 9.1: Meme from @KellyBodwin: if you can explain why this is funny, you understand seeds

9.3.2.3 Add arguments

You can just edit your function each time you want to calculate power for a different sample n, but it is more efficient to build this into your function as an arguments. Redefine t_sim, setting arguments for the mean and SD of group A, the mean and SD of group B, and the number of subjects per group. Give them all default values.

t_sim <- function(n = 10, m1=0, sd1=1, m2=0, sd2=1) {
  tibble(
    A = rnorm(n, m1, sd1),
    B = rnorm(n, m2, sd2)
  ) |>
    gather(group, score, A:B) |>
    t.test(score~group, data = _) |>
    broom::tidy() |>
    pull(p.value) 
}

9.3.3 Test your function

Test your function with some different values to see if the results make sense.

t_sim(100)
t_sim(100, 0, 1, 0.5, 1)

[1] 0.8460064
[1] 0.0002446404

Use replicate to calculate power for 100 subjects/group with an effect size of 0.2 (e.g., A: m = 0, SD = 1; B: m = 0.2, SD = 1). Use 1000 replications.

reps <- replicate(1000, t_sim(100, 0, 1, 0.2, 1))
power <- mean(reps < .05)
power

[1] 0.284

Compare this to power calculated from the power.t.test function.

power.t.test(n = 100, delta = 0.2, sd = 1, type="two.sample")


     Two-sample t test power calculation 

              n = 100
          delta = 0.2
             sd = 1
      sig.level = 0.05
          power = 0.2902664
    alternative = two.sided

NOTE: n is number in *each* group

Calculate power via simulation and power.t.test for the following tests:

20 subjects/group, A: m = 0, SD = 1; B: m = 0.2, SD = 1
40 subjects/group, A: m = 0, SD = 1; B: m = 0.2, SD = 1
20 subjects/group, A: m = 10, SD = 1; B: m = 12, SD = 1.5

9.4 Exercises

Create a custom function that takes a string of text as input and returns that text in title case (with the first letter of each word uppercase)
Does this function work if you give it a vector or 3 strings? How about if you give it a list of 3 strings? What happens if you give it a non-character data type, like numbers?
Add error checking and appropriate messages for when the user gives this function inappropriate input.
Create some tests for this function to make sure the output is as expected. What kinds of difficult edge cases can you think of?

Glossary

term	definition
argument	A variable that provides input to a function.

Further Resources

Chapters 19 of R for Data Science
Unit tests in Intro to R Packages.