STM1001: Computer Lab 9B

Data Science Module

Topic 9B: Writing R Functions

Welcome to the ninth Data Science module computer lab.

By now, you will have gained some experience using a variety of functions in R, ranging from inbuilt ‘base R’ functions, to package-specific functions, such as those contained in the plotly package.

In this computer lab we will extend our R coding skills and focus on writing functions in R.

By the end of this lab, you should understand the R code syntax required to create functions, and feel comfortable writing simple and intermediate functions. Let’s get started!

1 R Functions Overview

Usually, if you need to carry out a specific statistical process in R, there will already be a relevant function for you to use! However, when you are conducting your own coding or research, sometimes it can be extremely beneficial to write your own custom functions. These can range from simple functions of a few lines of code, to extremely detailed, multi-layered functions that call other functions and packages. The creative possibilities are endless.

Recall that we briefly discussed writing R functions several weeks ago, in Section 4.2 of the Data Visualisation in R supplement.

Before you continue with this lab, please read over the R function content in this supplement, to refresh your memory. It may also be helpful to keep this content open in a separate tab while you work through the lab material.

2 Writing Simple Functions

Once you have reread the content in Section 4.2 of the Data Visualisation in R supplement, let’s try to write some simple functions in R. These functions actually already exist in base R, but this will be good practice for subsequent work.

2.1 Mean Function

Suppose that you would like to write an R function that computes the mean. Recall that the mean \(\bar{x}\) of a set of \(n\) values \(x_1, x_2, \cdots, x_n\) is computed using the formula \[\bar{x} = \displaystyle \frac{\sum_{i=1}^n x_i}{n} .\]

There are several ways that we could write this function. We have made a head-start - take a look at the partially completed code below, and see if you can complete it by filling in the missing ... parts.

mean_func <- function(values){
  # Argument:
  # values: This is our list of values
  
  n <- length(values)
  sum(...) / ...

}                      

Note that, based on this format, our function can compute the mean of a string of numbers of any non-zero length.

2.1.1

Verify that your mean_func function works by computing the mean of the set of values 2, 3, 4, 5, 6, 7, 8. If your function is correct, you should receive an answer of 5.

2.2 Sample standard deviation Function

Now that you have had some practice writing R functions, let’s tackle a harder problem.

Recall the equation for the sample standard deviation \(s\) (introduced in Topic 2):

\[s = \displaystyle \sqrt{ \frac{1}{n-1} \sum_{i=1}^n (x_i - \bar{x})^2} .\] Using your mean_func function, and your code from 2.1 as a guide, write a function called sd_func that computes the sample standard deviation of a set of \(n\) values \(x_1, x_2, \cdots, x_n\).

Hints:

• You can use the same input as for your mean_func function.
• When you are multiplying parts of a function together, you need to use the \(\ast\) symbol.
• You can use the built-in R function sqrt() to compute the square root.
• You can use your mean_func function within your new function, when computing the \(\bar{x}\) value.
• Be careful to include brackets to separate different parts of your calculations.
• Don’t forget the square term in the equation.

2.2.1

Verify that your sample standard deviation function works by computing the sample standard deviation of the set of values 2, 3, 4, 5, 6, 7, 8. If your function is correct, you should receive an answer of 2.160247.

2.2.2

Check how your function compares to the inbuilt R function sd, which takes the one argument - the set of values. If your function is written correctly, you should obtain an identical result!

3 Writing a \(t\)-test Function

Now that we have had some practice writing R functions, let’s move on to a slightly harder scenario, and try to write an R function for computing a \(t\)-test.

Recall that in Topic 5 we introduced the \(t\)-test as a common type of hypothesis test.

The test statistic \(t\) used for the one-sample \(t\)-test has the following formula:

\[t = \displaystyle \frac{\bar{x} - \mu_0}{s / \sqrt{n}},\] where \(\bar{x}\), \(s\) and \(n\) are the sample mean, sample standard deviation, and sample size respectively, as previously defined, and \(\mu_0\) denotes the mean under the null hypothesis.

Note that \(s / \sqrt{n}\) is the standard error of \(t\), and is often denoted by \(SE\).

3.1 \(t\)-test test statistic

Using your functions mean_func and sd_func, write a new R function called t_test_func that computes the test statistic \(t\) for a one-sample \(t\)-test with an arbitrary null hypothesis mean \(\mu_0\).

Hints:

• Most of the hard work is already done, since you can use your two functions.
• You will need to provide one additional input argument.

3.2

To test your test statistic function, evaluate it for the set of values 2, 3, 4, 5, 6, 7, 8, with the null hypothesis mean set to 4.

Compare your result to the result of the inbuilt R function t.test, by running the code t.test(c(2:8), mu = 4). If you have written your function correctly, the t.test \(t=\) output should be identical to the value of \(t\) computed by your t_test_func function.

Hint: Don’t worry if your t_test_func function gives a slightly different value - as long as it is correct to 3 or 4 decimal places that is fine. If it is very different, just go back over your function, and check for errors such as insufficient or misplaced brackets. It’s very easy to make mistakes like this.

3.3 Degrees of freedom

Recall that the degrees of freedom for a one-sample \(t\)-test is equal to \(n-1\). It would be good to include this information in your t_test_func output.

As noted in Section 4.2.2 of the Data Visualisation in R supplement, if you want multiple results to be provided as the output of a function, you will need to combine these results into one output object.

There are several ways in which we can do this. The simplest option is to use the combine function, c(). For example, if your test statistic \(t\) result is stored in the object t.val, and your degrees of freedom result is stored in the object df, you could use c(t.val, df) to output your results.

• However, this won’t provide any extra information about your output - that might be fine if you are the only person using your function, but if you intend your function to be used by others, then adding more detail to the output might be desirable.

• As a simple extra step, continuing the example above you could use c("test statistic value" = t.val, "degrees of freedom" = df) to include some helpful additional details about your function output.

Using the information above, update your t_test_func function so that it outputs the test statistic \(t\) and the degrees of freedom, with clear labels.

Test your updated function, using the inputs discussed in 3.2.

3.4 The concatenate Function

As a further improvement, rather than using e.g. c("test statistic value" = t.val, "degrees of freedom" = df), we could instead use the concatenate function cat to provide more detailed output that is easier to read and interpret.

An example application of the cat function is provided below:

cat("The test statistic is", round(t.val, 4))

This code will output a single character vector, e.g. The test statistic is 1.2247, which is nicer to read than the c function method we used in 3.3.

Here, anything that is not an R code result is surrounded by quotation marks "...". A comma is used to separate sections of R code results and text (or text related code).

If we would like to include several statements in our output, we can use "\n" to specify a new line. For example, the layout for two lines could be:

cat("Some text here", some R code here, "\n", 
    "some more text", some more R code)

Using the information above, incorporate the cat function into your t_test_func so that the output is in a clear, easy to read format.

Test your updated function, using the inputs discussed in 3.2.

Note: This brief introduction does not cover all aspects of the concatenate function.

3.5 \(p\)-value Function

Suppose that we would also like to compute the \(p\)-value that corresponds to the \(t\) test statistic. Recall from Topic 5 that for a two-sided hypothesis test, the \(p\)-value is equal to the probability that we observe a test statistic value as extreme, or more extreme, than the observed test statistic we have computed using our sample.

In other words, we have that \[\text{$p$-value} = 1 - P(-t \leq T \leq t) ,\] i.e. \[\text{$p$-value} = P(T \leq -t) + P(T \geq t) .\]

Recall that in section 2 of Computer Lab 4 we discussed the use of functions like dnorm, pnorm and qnorm for computing the density, probability distribution function and quantile value respectively of normally distributed random variables.

Similarly, we can use the inbuilt R function pt to compute the probability that the test statistic \(T\) (which follows a Student’s \(t\) distribution) takes a value less than or equal to an arbitrary value. We will however also need to include the degrees of freedom as an argument.

For example, if we would like to compute \(P(T \leq 5)\), and our degrees of freedom is 4, we can use the code pt(5, 4). If we would like to compute \(P(T \geq 5)\), we can be clever and use the code 1 - pt(5, 4)¹.

Based on this information, add extra code to your t_test_func function so that it also provides (as a third output) the corresponding \(p\)-value for the test statistic.

Hint: Try to use the symmetry property of the Student’s \(t\)-distribution when including code for the \(p\)-value calculations in your function.

3.5.1

Run both your updated t_test_func function and the inbuilt t.test R function, using the inputs discussed in 3.2. Verify that your t_test_func function \(p\)-value is correct, by comparing it to the output of the t.test R function.

3.5.1.1 We now have written functions that, combined, can be used to conduct a one-sample \(t\)-test for a two-sided hypothesis. Well done!

4 Extension: Adding Additional Details to our \(t\)-test Function

To finish up this computer lab, we could add some additional details to our function.

4.1

Add comments with the hash symbol \(\#\) within your t_test_func function to explain what each line of code is doing - e.g.:

t_test_func <- function(...){
  
  ...
  
  # this code computes the degrees of freedom
  ...
  
  ...

}   

4.2

Use the list function to further modify your t_test_func function so that your output is a list, with each output named appropriately. Also provide the sample mean as a fourth output.

You can provide names for the objects you include in your list very easily, e.g. mean = mean.val will assign the name mean to your mean.val result.

4.2.1

Verify your updated code is working as intended using the inputs discussed in 3.2, and use the names function to check what results are available from your new t_test_func output.

Try using the $ operator to check specific values from your output, e.g. t.test.output$mean.

Great job, that concludes our R function writing.

Hopefully this lab has demystified R function creation, and shown you that writing functions in R is perhaps easier than you might have thought.

These notes have been prepared by Rupert Kuveke and Amanda Shaker. The copyright for the material in these notes resides with the authors named above, with the Department of Mathematical and Physical Sciences and with La Trobe University. Copyright in this work is vested in La Trobe University including all La Trobe University branding and naming. Unless otherwise stated, material within this work is licensed under a Creative Commons Attribution-Non Commercial-Non Derivatives License BY-NC-ND.

---

1. since we are dealing with a continuous random variable↩︎