Functions

Packages for this section

library(tidyverse)
library(broom) # some regression stuff later

Don’t repeat yourself xxx

  • Suppose you have a whole bunch of numbers you want to multiply by 5 and add 3:
a <- 50
b <- 11
d <- 3
as <- 5 * a + 3
as
[1] 253
bs <- 5 * b + 3
bs
[1] 58
ds <- 5 * d + 3
ds
[1] 18

What’s the problem?

  • Same calculation done three different times, by copying, pasting and editing.

  • Dangerous: what if you forget to change something after you pasted?

  • Programming principle: “don’t repeat yourself”.

  • Hadley Wickham: don’t copy-paste more than twice.

  • Instead: write a function.

Anatomy of function

  • Header line with function name and input value(s).
  • Body with calculation of values to output/return.
  • Return value: the output from function. In our case:
linear <- function(x) {
  5 * x + 3
}
  • If last line of function calculates value without saving it, that value is returned.
  • R also has a return but that is not usually needed.

About the input; testing 1/2

  • The input to a function can be called anything. Here we called it x. This is the name used inside the function.
  • The function is a “machine” for calculating the new version of its input. It doesn’t do anything until you call it:
linear(50)
[1] 253
linear(11)
[1] 58
linear(3)
[1] 18

Testing 2/2

q <- 17
linear(q)
[1] 88
linear("text")
Error in `5 * x`:
! non-numeric argument to binary operator
  • It works! (At least, it works when it should and fails when it should.)

Vectorization 1/2

  • We conceived our function to work on numbers

  • but it actually works on vectors too, as a free bonus of R:

linear(c(50, 11, 3))
[1] 253  58  18

More than one input

  • Allow the values to be multiplied and added (so far 5 and 3) also to be input to the function.
  • In the spirit of the equation of a line, call them slope and intercept:
linear2 <- function(x, slope, intercept) {
  x * slope + intercept
}

Calling the function

  • Call the function with the inputs in the right order:
linear2(10, 5, 3)
[1] 53
  • or give the inputs names, in which case they can be in any order:
linear2(intercept = 3, slope = 5, x = 51)
[1] 258

Defaults 1/2

  • Many R functions have values that you can change if you want to, but usually you don’t want to, for example:
x <- c(3, 4, 5, NA, 6, 7)
mean(x)
[1] NA
mean(x, na.rm = TRUE)
[1] 5
  • Mean of data with missing value is missing, unless you specify na.rm=TRUE: then, missing values removed first.

  • Or, na.rm has a default value of FALSE: that’s what it will be unless you change it.

Defaults 2/2

  • In our function, set default values for intercept and slope like this:
linear3 <- function(x, slope = 5, intercept = 3) {
  x * slope + intercept
}

Using defaults (or not)

  • If you specify a value for either or both of slope and intercept, they will be used. If you don’t, values 5 for slope and/or 3 for intercept will be used instead:
linear3(10, 2, -1) # 10 * 2 - 1
[1] 19
linear3(10, intercept = 12) # 10 * 5 + 12
[1] 62
linear3(10) # 10 * 5 + 3
[1] 53

but…

the value of x is not optional:

linear3(slope = 1, intercept = 4)
Error in `linear3()`:
! argument "x" is missing, with no default

Using R’s built-ins

  • When you write a function, you can use anything built-in to R, or even any functions that you defined before.
  • For example: you have data x and y, and you want to obtain a predicted value of y for an input value x_pred based on the regression of y on x.
  • Steps:
    • run the regression
    • pull out the intercept and slope
    • use function linear3 to use x_pred, the slope, and the intercept to work out the prediction.

Without a function

  • In worksheet 1, you met a dataframe mtcars which contained data on some cars, including gas mileage mpg and weight wt. What would be the predicted gas mileage of a car that weighs 4 tons?

  • Regression, but not obvious how to get slope and intercept:

mtcars.1 <- lm(mpg ~ wt, data = mtcars)
summary(mtcars.1)

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Without a function, cont’d

  • Use tidy from broom:
m <- tidy(mtcars.1)
m
  • Slope is second number in estimate column and intercept is first:
slope <- m$estimate[2]
intercept <- m$estimate[1]

Without a function, cont’d

  • and now use linear3:
x_pred <- 4
linear3(x_pred, slope, intercept)
[1] 15.90724

Making a function of this

  • Use “generic” names for things:
make_prediction <- function(x, y, x_pred) {
  reg.1 <- lm(y ~ x)
  m <- tidy(reg.1)
  my_slope <- m$estimate[2]
  my_intercept <- m$estimate[1]
  linear3(x_pred, my_slope, my_intercept)
}

Test function on mtcars data:

make_prediction(mtcars$wt, mtcars$mpg, 4)
[1] 15.90724
  • The answer is correctly what we got before.

Passing things on

  • lm has a lot of options, with defaults, that we might want to change.
  • Two of these are formula, the y ~ x thing, and data, the dataframe.
  • Instead of supplying x and y, we can make lm do the heavy lifting for us, thus:
make_prediction2 <- function(x_pred, ...) {
  reg.1 <- lm(...)
  m <- tidy(reg.1)
  my_slope <- m$estimate[2]
  my_intercept <- m$estimate[1]
  linear3(x_pred, my_slope, my_intercept)
}
  • The ... in the header line means “accept any other input”, and the ... in the lm line means “pass anything other than x_pred straight on to lm”.

Using ...

To use make_prediction2, instead of supplying an x and y, we supply a model formula and a dataframe, so that using this function looks a lot more like using lm:

make_prediction2(x_pred = 4, formula = mpg ~ wt, data = mtcars)
[1] 15.90724

and the answer is the same as before, but more aesthetically obtained, and the function will work with any variables in any dataframe.

Running a function for each of several inputs

  • We might want to do several predictions of mpg, for different values of wt, maybe these:
wts <- c(2, 3, 4)
wts
[1] 2 3 4
  • you might think of using a for loop for this
  • but R has a “functional programming” method called map, in this case map_dbl because our function make_predictions2 returns a decimal number (dbl). It goes like this:
map_dbl(wts, \(x) make_prediction2(x_pred = x,
                                   formula = mpg ~ wt,
                                   data = mtcars))

To run the map_dbl:

map_dbl(wts, \(x) make_prediction2(x_pred = x,
                                   formula = mpg ~ wt,
                                   data = mtcars))
[1] 26.59618 21.25171 15.90724
  • “for each of the values in wts, run make_predictions2 on that value as shown, combining the results together”
  • The last result is the we got before (predicted gas mileage for a car weighing 4 tons).
  • The answers look sensible (as weight increases, gas mileage in miles per gallon decreases)
  • you can check the others by running make_prediction2 yourself.

What if function returns more than one thing?

  • For example, finding quartiles:
quartiles <- function(x) {
  quantile(x, c(0.25, 0.75))
}
quartiles(mtcars$mpg)
   25%    75% 
15.425 22.800 
  • When function returns more than one thing, run repeatedly with map (or map_df) instead of map_dbl.

To work out quartiles of the wt and mpg columns in mtcars

mtcars %>% 
  select(wt, mpg) -> mtcars2
  • and then, “for each thing (column) in mtcars2, work out the quartiles of it”.
  • Method 1, gives a list:
map(mtcars2, \(x) quartiles(x))
$wt
    25%     75% 
2.58125 3.61000 

$mpg
   25%    75% 
15.425 22.800 

or…

Method 2, gives a dataframe:

map_df(mtcars2, \(x) quartiles(x))
  • first row: quartiles of wt
  • second row: quartiles of mpg
  • Principle: pretend quartiles returns one-row dataframe.

Map in data frames with mutate

  • map can also be used within data frames to calculate new columns.
  • calculate the predicted mpg for each wt in dataframe mtcars2:
mtcars2 %>% 
  mutate(pred_mpg = map_dbl(wt, \(w) make_prediction2(
    x_pred = w,
    formula = mpg ~ wt, 
    data = mtcars2
  )))

Write a function first and then map it

  • If the “for each” part is simple, go ahead and use map_-whatever.
  • If not, write a function to do the complicated thing first.
  • Example: “half or triple plus one”: if the input is an even number, halve it; if it is an odd number, multiply it by three and add one.
  • This is hard to do as a one-liner: first we have to figure out whether the input is odd or even, and then we have to do the right thing with it.

Odd or even?

  • Odd or even? Work out the remainder when dividing by 2:
6 %% 2
[1] 0
5 %% 2
[1] 1
  • 5 has remainder 1 so it is odd.

Write the function

  • First test for integerness, then test for odd or even, and then do the appropriate calculation:
hotpo <- function(x) {
  stopifnot(round(x) == x) # passes if input an integer
  remainder <- x %% 2
  if (remainder == 1) { # odd number
    ans <- 3 * x + 1
  }
  else { # even number
    ans <- x %/% 2 # integer division
  }
  ans
}

Test it

hotpo(3)
[1] 10
hotpo(12)
[1] 6
hotpo(4.5)
Error in `hotpo()`:
! round(x) == x is not TRUE

One through ten

  • Use a data frame of numbers 1 through 10 again:
tibble(x = 1:10) %>% mutate(y = map_int(x, \(x) hotpo(x)))

Until I get to 1 (if I ever do)

  • If I start from a number, find hotpo of it, then find hotpo of that, and keep going, what happens?
  • If I get to 4, 2, 1, 4, 2, 1 I’ll repeat for ever, so let’s stop when we get to 1:
hotpo_seq <- function(x) {
  ans <- x
  while (x != 1) {
    x <- hotpo(x)
    ans <- c(ans, x)
  }
  ans
}
  • Strategy: keep looping “while x is not 1”.
  • Each new x: add to the end of ans. When I hit 1, I break out of the while and return the whole ans.

Trying it 1/2

  • Start at 6:
hotpo_seq(6)
[1]  6  3 10  5 16  8  4  2  1

Trying it 2/2

  • Start at 27:
hotpo_seq(27)
  [1]   27   82   41  124   62   31   94   47  142   71  214
 [12]  107  322  161  484  242  121  364  182   91  274  137
 [23]  412  206  103  310  155  466  233  700  350  175  526
 [34]  263  790  395 1186  593 1780  890  445 1336  668  334
 [45]  167  502  251  754  377 1132  566  283  850  425 1276
 [56]  638  319  958  479 1438  719 2158 1079 3238 1619 4858
 [67] 2429 7288 3644 1822  911 2734 1367 4102 2051 6154 3077
 [78] 9232 4616 2308 1154  577 1732  866  433 1300  650  325
 [89]  976  488  244  122   61  184   92   46   23   70   35
[100]  106   53  160   80   40   20   10    5   16    8    4
[111]    2    1

Which starting points have the longest sequences?

  • The length of the vector returned from hotpo_seq says how long it took to get to 1.
  • Out of the starting points 1 to 100, which one has the longest sequence?

Top 10 longest sequences

tibble(start = 1:100) %>%
  mutate(seq_length = map_int(
    start, \(start) length(hotpo_seq(start)))) %>%
  slice_max(seq_length, n = 10)
  • 27 is an unusually low starting point to have such a long sequence.

What happens if we save the entire sequence?

tibble(start = 1:7) %>%
  mutate(sequence = map(start, \(start) hotpo_seq(start)))
  • Each entry in sequence is itself a vector. sequence is a “list-column”.

Using the whole sequence to find its length and its max

tibble(start = 1:7) %>%
  mutate(sequence = map(start, \(start) hotpo_seq(start))) %>%
  mutate(
    seq_length = map_int(sequence, \(sequence) length(sequence)),
    seq_max = map_int(sequence, \(sequence) max(sequence))
  )

Does it work with rowwise?

tibble(start=1:7) %>% 
  rowwise() %>% 
  mutate(sequence = list(hotpo_seq(start))) %>% 
  mutate(seq_length = length(sequence)) %>% 
  mutate(seq_max = max(sequence))

It does.

Final thoughts on this

  • Called the Collatz conjecture.
  • Nobody knows whether the sequence always gets to 1.
  • Nobody has found an \(n\) for which it doesn’t.
  • A tree (link).