How do you program a custom hypothesis test in R?

Question

There are a number of interesting hypothesis tests programmed in R, where the output of the test appears in a nice user-friendly format. For example, here is the output of Welch's T-test using the t.test function in the stats package.

#Run a T-test on some example data
X <- c(30, 32, 40, 28, 29, 35, 30, 34, 31, 39);
Y <- c(19, 20, 44, 45, 8, 29, 26, 59, 35, 50);
TEST <- stats::t.test(X,Y);

#Print the TEST object
TEST;

        Welch Two Sample t-test

data:  X and Y
t = -0.13444, df = 10.204, p-value = 0.8957
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -12.27046  10.87046
sample estimates:
mean of x mean of y 
     32.8      33.5

As you can see, the output for this test is in a user-friendly format that gives all the required information relating to the output of the test. It gives this output in a format that is different to the standard output when you call an R object.

---

Question: Suppose you encounter a type of hypothesis test that is not in any of the existing R packages, and you want to program this test as a function so that you can run it easily on new data, and get a nice user-friendly output like the one above. How do you program this?

Answer 1

Here's what you need to do in general

Hypothesis testing functions in R create and output a list object of class h.test. This type of object has a specific set of required components set out in its documentation, and it also has a special method of printing under the print.htest setting in the global environment. That printing method draws out information from the list, but prints it in the user-friendly way you see in the output in the question. The list should contain the components set out below, including naming several of the objects with a names attribute. (You are some other optional components shown in the linked documentation.)

---

Textual description of test

• method: A character string giving the name of the hypothesis test. This will appear as the first sentence of the print output.

• data.name: A character string giving a description of the data, which usually includes reference to the names of the data vectors used in the test. For this part it is useful to use the substitute and deparse functions to extract the names of the user inputs to the function as the appropriate names (example shown below).

---

Specification of hypotheses

• null.value: A numeric variable giving the value of the parameter under the null hypothesis (with a names attribute).

• alternative: A character string set to greater, less or two-sided, to specify the direction of the alternative hypothesis relative to the null value.

---

Test statistic and p-value

• estimate: The estimated value of the parameter (with a names attribute). This value will be a function of the data inputs for the testing function.

• statistic: The value of the test statistic (with a names attribute). This value will either be a direct function of the data inputs for the testing function, or a function of the parameter estimate.

• p.value: The p-value for the test (which should be a number between zero and one). This value will be a function of the test statistic.

---

Confidence interval (optional)

• conf.int: A confidence interval represented by a vector with two elements, where the first is the lower bound and the second is the upper bound (with a conf.level attribute giving the confidence level). If you are using this component, it is desirable to require the function to take a significance level as an input, in order to specify the desired confidence level.

---

In order to create a custom hypothesis-testing function, you will need to create a function that produces a list containing the required components shown above, customised to your particular test. For the substantive parts of the test (i.e., the estimate, test statistic, p-value, and confidence interval), you will need to use the appropriate formulae for your particular test. Note that you can put these elements in any order in your list, so long as all the required elements are there.

You can also add other components to the list if you wish. It is good practice to add an initial part of your function to check the inputs to the function, to ensure that they are of the correct form, and to stop the function and give error messages if the input is defective in some way. Once your list is created, you set the class of the object to h.test and output the object at the end of the function.

---

Here's an example of implementation for a particular test

In a related question I gave an example of code for a hypothesis test taken from Tarone (1979). Below is a slightly modified version of that code that serves as an example for how you can program a function for a custom hypothesis test.

Observe that the code first checks the inputs, and then builds up each of the required components of the test, using the appropriate names and formulae for that specific test. Once these components have been computed, we create a list object called TEST, composed of these elements, and we set its class to h.test. We output this object at the end of the function. (It is also worth observing the code for data.name, which extracts the variable names that are input by the user.)

Tarone.test <- function(N, M) {
    
    #Check validity of inputs
    if(!(all(N == as.integer(N)))) { stop("Error: Number of trials should be integers"); }
    if(min(N) < 1) { stop("Error: Number of trials should be positive"); }
    if(!(all(M == as.integer(M)))) { stop("Error: Count values should be integers"); }
    if(min(M) < 0) { stop("Error: Count values cannot be negative"); }
    if(any(M > N)) { stop("Error: Observed count value exceeds number of trials"); }
    
    #Set description of test and data
    method      <- "Tarone's Z test";
    data.name   <- paste0(deparse(substitute(M)), " successes from ", 
                          deparse(substitute(N)), " trials");
    
    #Set null and alternative hypotheses
    null.value  <- 0;
    attr(null.value, "names") <- "dispersion parameter";
    alternative <- "greater";
    
    #Calculate test statistics
    estimate    <- sum(M)/sum(N);
    attr(estimate, "names") <- "proportion parameter";
    S           <- ifelse(estimate == 1, sum(N),
                          sum((M - N*estimate)^2/(estimate*(1 - estimate))));
    statistic   <- (S - sum(N))/sqrt(2*sum(N*(N-1))); 
    attr(statistic, "names") <- "z";
    
    #Calculate p-value
    p.value     <- 2*pnorm(-abs(statistic), 0, 1);
    attr(p.value, "names") <- NULL;
    
    #Create htest object
    TEST        <- list(method = method, data.name = data.name,
                        null.value = null.value, alternative = alternative,
                        estimate = estimate, statistic = statistic, p.value = p.value);
    class(TEST) <- "htest";
    TEST; }

Below we create some count data to implement this test and see what the output looks like. As you can see, the output is the same user-friendly output you get for other hypothesis tests in R, where the components of the test have been pulled out of the list and presented in a nice simple manner. The output shows the name of the test and describes the data, and then it gives the test statistic and p-value for the test. It also describes the alternative hypothesis and gives the sample estimate of the parameter.

#Generate example data
TRIALS <- c(30, 32, 40, 28, 29, 35, 30, 34, 31, 39);
COUNTS <- c( 9, 10, 22, 15,  8, 19, 16, 19, 15, 10);

#Apply Tarone's test to the example data
TEST <- Tarone.test(TRIALS, COUNTS);
TEST;

        Tarone's Z test

data:  COUNTS successes from TRIALS trials
z = 2.5988, p-value = 0.009355
alternative hypothesis: true dispersion parameter is greater than 0
sample estimates:
proportion parameter 
           0.4359756 

Answer 2

The hypothesis test functions in the stats package use classic S3 object-orientated programming. You write a function that creates a "htest" object, which is a list with a standard set of components, and R has a built-in print method for objects of that class. The user-level function is traditionally called something like yourname.test but can have any name. It can have any appropriate arguments.

• Type ?t.test to see the definition of a "htest" object.
• See stats:::t.test.default to see an example of a function that creates a "htest" object.
• See stats:::print.htest to see how the user-friendly output is created.

Here is a toy example that performs a very simple chisquare test:

demo.test <- function(s2, df=1)
{
  pval <- pchisq(s2, df, lower.tail=FALSE)
  out <- list(
    statistic=s2,
    parameter=NULL,
    p.value=pval,
    null.value=NULL,
    alternative="greater",
    method="demo",
    data.name="s2")
  class(out) <- "htest"
  out
}

Then

> TEST <- demo.test(30, df=10)
> TEST

        demo

data:  s2
= 30, p-value = 0.0008566
alternative hypothesis: greater

If you want to be fancier, you can make your function S3 generic (like the stats package functions) in order to handle different types of input (e.g., a formula instead of data vectors). But an ordinary function like the above example might satisfy your needs.