permutation hypothesis test in r

• Data Visualization
• Statistics in R
• Machine Learning in R
• Data Science in R
• Packages in R

Permutation Hypothesis Test in R Programming

In simple words, the permutation hypothesis test in R is a way of comparing a numerical value of 2 groups. The permutation Hypothesis test is an alternative to: 

• Independent two-sample t-test 
• Mann-Whitney U aka Wilcoxon Rank-Sum Test

Let’s implement this test in R programming .

Why use the Permutation Hypothesis Test?  

• Small Sample Size. 
• Assumptions(for parametric approach) not met. 
• Test something other than classic approaches comparing Means and Medians. 
• Difficult to estimate the SE for test-statistic.

Permutation Hypothesis Test Steps

• Specify a hypothesis 
• Choose test-stat(Eg: Mean, Median, etc. ) 
• Determine Distribution of test-stat 
• Convert test-stat to P-value 

Note: P-value = No. of permutations having a test-stat value greater than observed test-stat value/ No. of permutations.

Implementation in R

• Dataset: Chicken Diet Data. This dataset is a subset of the “chickwts” data in the “R dataset package”. Download the data set here .
• Hypothesis: The weight of the chicken is independent of the type of diet.

Test-Statistics

• Test-Statistics #1: The absolute value of the difference in mean weights for the two diets | Y 1 – Y 2 | . This is the same test statistics as the independent two-sided two-sample t-test.
• Test-Statistics #2: The absolute value of the difference in median weights for the two diets | Median 1 – Median 2 |

Please Login to comment...

Similar reads.

• R-Statistics

Improve your Coding Skills with Practice

What kind of Experience do you want to share?

Statistical inference with permutation tests

To learn about the basics of permutation tests and statistical resampling from an excellent textbook, see @resampling-book. For a primer on hypothesis testing with permutation tests in the context of topological data analysis, see @hyptest. Since the distribution of topological features has not been well characterized yet, statistical inference on persistent homology must be nonparametric. Given two sets of data, \(X\) and \(Y\), conventional statistical inference generally involves comparison of the parameters of each population with the following null and alternative hypotheses:

\[ \begin{aligned} H_0&: \mu_X=\mu_Y \ H_A&: \mu_X\neq\mu_Y \end{aligned} \] If we define a function \(T\) that returns the persistent homology of a point cloud, then given two point clouds, \(C\) and \(D\), we can use a permutation test to conduct analogous statistical inference with the following null and alternative hypotheses:

\[ \begin{aligned} H_0&:T©=T(D) \ H_A&:T©\neq T(D) \end{aligned} \] TDAstats uses the Wasserstein distance (aka Earth-mover's distance) as a similarity metric between persistent homologies of two point clouds [@wasserstein-calc]. Although visual analysis of plots (topological barcodes and persistence diagrams) is essential, a formal statistical procedure adds objectivity to the analysis. The case study below highlights the main features of TDAstats pertaining to statistical inference. For practice, perform the steps of the case study to the unif3d and sphere3d datasets.

Case study: unif2d versus circle2d

To ensure that all the code output in this section is reproducible, we set a seed for R's pseudorandom number generator. We are also going to need the unif2d and circle2d datasets provided with TDAstats, so we load them right after setting the seed.

The unif2d dataset is a numeric matrix with 100 rows and 2 columns containing the Cartesian x- and y-coordinates (columns 1 and 2, respectively) for 100 points (1 per row). The points are uniformly distributed within the unit square with corners \((0, 0)\), \((0, 1)\), \((1, 1)\), and \((1, 0)\). We confirm this with the following scatterplot.

The points do appear uniformly distributed as described above. Next, we take a look at the circle2d dataset, which is also a numeric matrix with 100 rows and 2 columns. However, circle2d contains the Cartesian x- and y-coordinates for 100 points uniformly distributed on the circumference of a unit circle centered at the origin. Like we did with unif2d , we confirm this with a scatterplot.

The points indeed appear to be uniformly distributed on a unit circle.

Before we use a permutation test to see if unif2d and circle2d exhibit distinct persistent homologies, we should take a look at the topological barcodes of each. Since we have 2-dimensional data, we are primarily concerned with the presence of 0-cycles and 1-cycles. If points were connected to each other by edges in a distance-dependent manner, then the resulting graphs (assuming a “good” distance-dependence) for unif2d and circle2d would have a single major component. Thus, we do not expect interesting behavior in the 0-cycles for either dataset. There also does not appear to be a prominent 1-cycle for the points in unif2d . However, the circle2d dataset was intentionally designed to have a single prominent 1-cycle containing all the points in the dataset. Thus, when we plot the topological barcodes for circle2d we should see a persistent 1-cycle that we do not see in the barcode for unif2d . We confirm our expectations with the following code.

We note two aspects of the topological barcodes above: (1) the limits of the horizontal axis are very different making direct comparison difficult; (2) it could be confusing to tell which barcode corresponds to which dataset. To fix these issues and demonstrate how the topological barcodes can be modified with ggplot2 functions ( plot_barcode returns a ggplot2 object), we run the following code.

We can safely ignore the warnings printed by ggplot2. Rescaling the horizontal axis had two major effects. First, we notice that the 0-cycles which appeared far more persistent for unif2d than for circle2d are now comparable. Second, the 1-cycles in unif2d are not persistent after the rescaling operation. Since the only prominent 1-cycle is now in circle2d , our expectations with respect to the topological barcodes were correct. We can now run a permutation test on the two datasets to confirm that the persistent homologies of the two are, in fact, distinct. To do this, all we have to do is use the permutation_test function in TDAstats, and specify the number of iterations. Increasing the number of iterations improves how well the permutation test approximates the distribution of all point permutations between the two groups, but also comes at the cost of speed. Thus, a number of iterations that is sufficiently large to properly approximate the permutation distribution but not too large to be computed is required. Almost certainly, the ideal number of iterations will change as the available computing power changes.

Note that the printed p-values for each set of cycles are unadjusted p-values. To see how p-values can be adjusted for permutation tests, see @resampling-book. You may also want to look at the null distributions generated by the permutation test for each dimension as follows.

Given that both vertical lines are far right of the plotted histograms (corresponding to the p-values of zero), we can conclude safely that the permutation test has given us sufficient evidence to reject the null hypothesis. Thus, the persistent homologies of unif2d and circle2d appear to be significantly different.

N.B.: persistence diagrams (using the plot_persist function) could replace the topological barcodes above. However, since the vertical and horizontal axes are important in persistence diagrams, the ylim ggplot2 function would also have to be used to rescale axes.

For practice, you can repeat the case study for the unif3d and sphere3d datasets. Keep in mind that the dim parameter in the calculate_homology function would likely have to be changed and that you will have a third permutation distribution generated that would need to be plotted.

Want to create or adapt books like this? Learn more about how Pressbooks supports open publishing practices.

Group Comparisons

16 Permutation Tests

Learning Objectives

To explore the theory behind permutation-based tests.

To illustrate how permutation tests can be conducted in R.

Key Packages

require(tidyverse)

Introduction

A statistical test involves the calculation of a test statistic followed by an assessment of how likely the calculated value of the test statistic would be if the data were randomly distributed.  In the case of ANOVA , the test statistic is the F -statistic, and it is compared to the theoretical distribution of F -values with the same degrees of freedom.

We will consider a number of tests using a range of test statistics.  However, there is no theoretical distribution for these test statistics.  Rather, the distribution of the test statistic will be derived from the data.   Generally, this reference distribution is generated by permuting (i.e., randomly reordering) the group identities, recalculating the test statistic, saving that value, and repeating this process many times (Legendre & Legendre 2012).

If the patterns observed in the data are unlikely to have arisen by chance, then the actual value of the test statistic should differ from the set of values obtained from the permutations.

Key Takeaways

The number of sample units directly affects the number of permutations.

Some permutations are functionally equivalent, so there are fewer combinations for a given sample size.

Number of Permutations (and Combinations)

The number of samples directly affects the number of possible permutations and combinations .

A permutation is a re-ordering of the sample units.  From a total of n sample units, the number of possible permutations ( P ) is:

[latex]P = n![/latex]

The number of permutations rises rapidly with sample size – see the table below.  To explore this, enumerate the permutations of the letters {a, b, C, D}.  There are four values, so the number of permutations is:

[latex]P = n! = 4 \cdot 3 \cdot 2 \cdot 1 = 24[/latex]

Assume that the letters {a, b, C, D} represent sample units, with lower and upper cases representing different groups.  In other words, the first two sample units are in one group and the last two are in the other group.  One permutation of these letters is {a, C, b, D}.  In this permutation, we assign the first and third sample units to one group and the second and fourth sample units to the other group.  Note that this ‘assignment’ is temporary and only for the purpose of this permutation.

When we think about group comparisons, it is helpful to recognize that some permutations are functionally equivalent.  For example, consider the permutations {a, C, b, D} and {C, a, D, b}.  In both permutations, the first and third sample units are assigned to one group and the second and fourth sample units are assigned to the other group.  These permutations represent the same combination.  A combination is a unique set of sample units, irrespective of sample order within each set.  The number of combinations ( C ) of size r is:

[latex]C_{r}^{n} = \frac{n!}{r! (n - r)!}[/latex]

(this equation is from Burt & Barber 1996).  This is for equally-sized groups; the calculations are more complicated for groups of different sizes.

For our simple example of four sample units, there are

[latex]C_{r}^{n} = \frac{n!}{r! (n - r)!} = \frac{4!}{2! (4 - 2)!} = \frac{4 \cdot 3 \cdot 2 \cdot 1}{2 \cdot 1 (2 \cdot 1)} = 6[/latex]

The following table shows how the number of permutations and combinations (into two equally-sized groups) rises rapidly with the number of sample units.

Probabilities

Permutation-based probabilities are calculated as the proportion of permutations in which the computed value of the test statistic is equal to or more extreme than the actual value.  This calculation can be made with any number of permutations, though it is easier to do so mentally if the denominator is a multiple of ten, such as 1,000.

The actual sequence of group identities is one of the possible permutations and therefore is included in the denominator of the probability calculation.   This is why it is common to do, for example, 999 permutations – once the actual sequence is included, the denominator of the probability calculation is 1,000.

The minimum possible P -value is partly a function of the number of permutations.  For example, consider a scenario in which the test statistic is larger when calculated with the real data than when calculated for any of the permutations.  If we had only done 9 permutations, we would calculate [latex]P = \frac{1}{9 + 1} = 0.1[/latex] However, if we had done 999 permutations, we would calculate [latex]P = \frac{1}{999 + 1} = 0.001[/latex]

Would it make sense to declare that the effect is significant in the second case but not in the first?

For studies with reasonably large sample sizes, there are many more permutations than we can reasonably consider.  For example, there are 8.2 × 10 47 permutations of 40 samples as reported in the above table.  If we considered one every millisecond (i.e., 1000 per second), it would still take us 2.6 × 10 37 years to consider them all!

Given the large number of possible permutations, we usually assess only a small fraction of them.  T his means that permutation-based probability estimates are subject to sampling error and vary from run to run.  The variation between runs declines as the number of permutations increases; more permutations will result in more consistent estimates of the probability associated with a test statistic.  Legendre & Legendre (2012) offer the following recommendations about how many permutations to compute:

• Use 500 to 1000 permutations during exploratory data analyses
• Rerun with an increased number of permutations if the computed probability is close to the preselected significance level (either above or below)
• Use more permutations (~10,000) for final, published results

The system.time() function can be used to measure how long a series of permutations requires.

The statistical significance of a permutation-based test is the proportion of permutations in which the computed value of the test statistic is equal to or more extreme than the actual value.

The actual value of the test statistic is unaffected by the number of permutations.

It’s ok to use small numbers of permutations during exploratory analyses, but use a large number (~10,000) for final analyses.

Exchangeable Units

I mentioned two assumptions of ANOVA / MANOVA in that chapter ; a third foundational assumption of both techniques is that the sample units are independent .  Analyses may be suspect when this is not the case or when the lack of independence is not properly accounted for.  Failure to account for lack of independence is a type of pseudoreplication (Hurlbert 1984).

The assumption of independence also applies for permutation tests – it is what justifies exchangeability in a permutation test.  This means that it is possible to analyze a permutation-based test incorrectly.   Permutations need to be restricted when sample units are not exchangeable.    The correct way of permuting data depends on the structure of the study and the hypotheses being tested.  The basic idea is that  the exchangeable units that would form the denominator when testing a term in a conventional ANOVA are those that should be permuted during a permutation test of that term .

Questions about independence and exchangeability are particularly pertinent for data obtained from complex designs that include multiple explanatory variables simultaneously.  See Anderson & ter Braak (2003), Anderson et al. (2008), and Legendre & Legendre (2012) for details on how to identify the correct exchangeable units for a permutation test.  

For example, in a split-plot design one factor can be applied to whole plots and another factor to split plots (i.e., within the whole plots).  Each of these factors would require a different error term.

• Analyses of the whole plot factor use the unexplained variation among whole plots as the error term. This is evident in the fact that the df for the whole plot error term is based on the number of whole plots, regardless of how many measurements were made within them.  In a permutation test, variation among whole plots is assessed by restricting permutations such that all observations from the same whole plot are permuted together.  Variation within whole plots is ignored when analyzing whole plot effects.
• Analyses of the split-plot factor use the residual as the error term and therefore do not require restricted permutations. However, they do require the inclusion of a term that uniquely identifies each whole plot so that the variation among whole plots is accounted for.  Doing so allows the analysis to focus on the variation within whole plots.  If a model included interactions with the split-plot factor, these would also be tested at this scale.

More information on this topic is provided in the chapter about complex models .

Implementation in R

The sample() function can be used to permute data.  If needed, you can use the size argument to create a subset, and the replace argument to specify whether to sample with replacement (by default, this is FALSE ).

The vegan package, drawing on the permute package, includes a number of options for conducting permutations.  This topic is explained in more detail in the chapters about controlling permutations and restricting permutations .

Simple Example, Graphically

Since our simple example only has two response variables, it is easily visualized: library(tidyverse) ggplot(data = perm.eg, aes(x = Resp1, y = Resp2)) +   geom_point(aes(colour = Group, shape = Group), size = 5) +   labs(title = "Real Data") +   theme_bw() ggsave("graphics/main.png", width = 3, height = 2.5, units = "in", dpi = 300) (did you notice that we saved the image, and where we saved it?)

To conduct a permutation, we can permute either the grouping factor or the data.  Can you see why these are equivalent?  We would not permute both at the same time … do you see why that is?

We’ll permute the grouping factor: perm.eg$perm1 <- sample(perm.eg$Group)

We’ve added the permutation as a new column within the perm.eg object.  View the object to compare the permutation with the original grouping factor.  Note that the number of occurrences of each group remains the same in permutations as in the original.

Let’s visualize this permutation of the data.  We can use the same code as above with a few changes:

• Name of column identifying the groups used for colour and shape  within geom_point() .
• Title of figure
• Name of file to which image is saved

ggplot(data = perm.eg, aes(x = Resp1, y = Resp2)) +   geom_point(aes(colour = perm1, shape = perm1), size = 5) +   labs(title = "Permutation 1") +   theme_bw() ggsave("graphics/perm1.png", width = 3, height = 2.5, units = "in", dpi = 300)

It is possible but somewhat unlikely that the group identities in your graph match this one.  Be sure you understand why!

Simple Example, Distance Matrix

Analyses will be based on the distance matrix so let’s consider this.  Here it is :

Resp.dist <- perm.eg |>   dplyr::select(Resp1, Resp2) |>   dist() round(Resp.dist, 3)

Group identity was not part of the distance matrix calculation.  This means that permuting the group identities doesn’t change the distance matrix itself.

Permuting group identities does change which distances connect sample units assigned to the same group.  For example, plots 1, 2, and 3 are all in group A in our real data but plots 2, 4, and 6 were assigned to group A in Permutation 1 above.

Note : To keep our example simple, we did not relativize the data and calculated Euclidean distances.  These decisions do not affect the permutations but, as discussed before, these decisions should be based on the nature of your data and your research questions.

Grazing Example

Our grouping factor for this example is current grazing status (Yes, No).  The sample() function can also be applied here: sample(grazing)

We’ll use this example in more detail in upcoming chapters.

Anderson, M.J., R.N. Gorley, and K.R. Clarke. 2008. PERMANOVA+ for PRIMER: guide to software and statistical methods . PRIMER-E Ltd, Plymouth Marine Laboratory, Plymouth, UK. 214 p.

Anderson, M.J., and C.J.F. ter Braak. 2003. Permutation tests for multi-factorial analysis of variance. Journal of Statistical Computation and Simulation 73:85-113.

Burt, J.E., and G.M. Barber. 1996. Elementary Statistics for Geographers . 2nd edition. Guilford Publications.

Hurlbert, S.H. 1984. Pseudoreplication and the design of ecological field experiments. Ecological Monographs 54:187-211.

Legendre, P., and L. Legendre. 2012. Numerical ecology . 3rd English edition. Elsevier, Amsterdam, The Netherlands.

Media Attributions

Applied Multivariate Statistics in R Copyright © 2024 by Jonathan D. Bakker is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License , except where otherwise noted.

Share This Book

Permutation Tests

An increasingly common statistical tool for constructing sampling distributions is the permutation test (or sometimes called a randomization test). Like bootstrapping, a permutation test builds - rather than assumes - sampling distribution (called the “permutation distribution”) by resampling the observed data. Specifically, we can “shuffle” or permute the observed data (e.g., by assigning different outcome values to each observation from among the set of actually observed outcomes). Unlike bootstrapping, we do this without replacement.

Permutation tests are particularly relevant in experimental studies, where we are often interested in the sharp null hypothesis of no difference between treatment groups. In these situations, the permutation test perfectly represents our process of inference because our null hypothesis is that the two treatment groups do not differ on the outcome (i.e., that the outcome is observed independently of treatment assignment). When we permute the outcome values during the test, we therefore see all of the possible alternative treatment assignments we could have had and where the mean-difference in our observed data falls relative to all of the differences we could have seen if the outcome was independent of treatment assignment. While a permutation test requires that we see all possible permutations of the data (which can become quite large), we can easily conduct “approximate permutation tests” by simply conducting a vary large number of resamples. That process should, in expectation, approximate the permutation distribution.

For example, if we have only n=20 units in our study, the number of permutations is:

That number exceeds what we can reasonably compute. But we can randomly sample from that permutation distribution to obtain the approximate permutation distribution, simply by running a large number of resamples. Let's look at this as an example using some made up data:

The difference in means is, as we would expect (given we made it up), about 1:

To obtain a single permutation of the data, we simply resample without replacement and calculate the difference again:

Here we use the permuted treatment vector s instead of tr to calculate the difference and find a very small difference. If we repeat this process a large number of times, we can build our approximate permutation distribution (i.e., the sampling distribution for the mean-difference). We'll use replicate do repeat our permutation process. The result will be a vector of the differences from each permutation (i.e., our distribution):

We can look at our distribution using hist and draw a vertical line for our observed difference:

At face value, it seems that our null hypothesis can probably be rejected. Our observed mean-difference appears to be quite extreme in terms of the distribution of possible mean-differences observable were the outcome independent of treatment assignment. But we can use the distribution to obtain a p-value for our mean-difference by counting how many permuted mean-differences are larger than the one we observed in our actual data. We can then divide this by the number of items in our permutation distribution (i.e., 2000 from our call to replicate , above):

Using either the one-tailed test or the two-tailed test, our difference is unlikely to be due to chance variation observable in a world where the outcome is independent of treatment assignment.

library(coin)

We don't always need to build our own permutation distributions (though it is good to know how to do it). R provides a package to conduct permutation tests called coin . We can compare our p-value (and associated inference) from above with the result from coin :

Clearly, our approximate permutation distribution provided the same inference and a nearly identical p-value. coin provides other permutation tests for different kinds of comparisons, as well. Almost anything that you can address in a parametric framework can also be done in a permutation framework (if substantively appropriate). and anything that coin doesn't provide, you can build by hand with the basic permutation logic of resampling.

Permutation t-Test

Description.

Performs one and two sample permutation t-tests on vectors of data.

The implemented test corresponds to the proposal of Chapter 15 of Efron and Tibshirani (1993) for equal variances as well as Janssen (1997) respectively Chung and Romano (2013) for unequal variances.

The function returns permutation p values and confidence intervals as well as the results ot the t-test without permutations.

The formula interface is only applicable for the 2-sample tests.

alternative = "greater" is the alternative that x has a larger mean than y .

If paired is TRUE then both x and y must be specified and they must be the same length. Missing values are silently removed (in pairs if paired is TRUE ). If var.equal is TRUE then the pooled estimate of the variance is used. By default, if var.equal is FALSE then the variance is estimated separately for both groups and the Welch modification to the degrees of freedom is used.

If the input data are effectively constant (compared to the larger of the two means) an error is generated.

A list with class "perm.htest" (derived from class htest ) containing the following components:

Code and documentation are for large parts identical to function t.test .

B. Efron, R.J. Tibshirani. An Introduction to the Bootstrap . Chapman and Hall/CRC 1993.

A. Janssen (1997). Studentized permutation tests for non-i.i.d, hypotheses and the generalized Behrens-Fisher problem. Statistics and Probability Letters , 36 , 9-21.

E. Chung, J.P. Romano (2013). Exact and asymptotically robust permutation tests. The Annals of Statistics , 41 (2), 484-507.

t.test , meanCI , meanDiffCI , boot.t.test