Implementing a Class of Permutation Tests: The coin Package

This description of the R package coin is a (slightly) modiﬁed version of Hothorn, Hornik, van de Wiel, and Zeileis (2008a) published in the Journal of Statistical Software . The R package coin implements a uniﬁed approach to permutation tests providing a huge class of independence tests for nominal, ordered, numeric, and censored data as well as multivariate data at mixed scales. Based on a rich and ﬂexible conceptual framework that embeds diﬀerent permutation test procedures into a common theory, a computational framework is established in coin that likewise embeds the corresponding R functionality in a common S4 class structure with associated generic functions. As a consequence, the computational tools in coin inherit the ﬂexibility of the underlying theory and conditional inference functions for important special cases can be set up easily. Conditional versions of classical tests—such as tests for location and scale problems in two or more samples, in-dependence in two-or three-way contingency tables, or association problems for censored, ordered categorical or multivariate data—can easily be implemented as special cases using this computational toolbox by choosing appropriate transformations of the observations. The paper gives a detailed exposition of both the internal structure of the package and the provided user interfaces along with examples on how to extend the implemented functionality.


Introduction
Conditioning on all admissible permutations of the data for testing independence hypotheses is a very old, yet very powerful and popular, idea (Fisher 1935;Ernst 2004).Conditional inference procedures, or simply permutation or re-randomization tests, are implemented in many different statistical computing environments.These implementations, for example wilcox.test()for the Wilcoxon-Mann-Whitney test or mantelhaen.test()for the Cochran-Mantel-Haenszel χ 2 test in R (R Core Team 2019) or the tools implemented in StatXact (Cytel Inc. 2003), LogXact (Cytel Inc. 2006), or Stata (StataCorp. 2003)-see Oster (2002Oster ( , 2003) ) for an overview-all follow the classical classification scheme of inference procedures and offer procedures for location problems, scale problems, correlation, or nominal and coin: Implementing a Class of Permutation Tests group time control 300 300 300 300 300 300 300 300 300 300 300 300 treatment 18 22 75 163 271 300 300 300 300 300 300 300 Table 1: The rotarod data: length of time on rotating cylinder by group.
ordered categorical data.Thus, each test procedure is implemented separately, maybe with the exception of conditional versions of linear rank statistics (Hájek, Šidák, and Sen 1999) in PROC NPAR1WAY as available in SAS (SAS Institute Inc. 2003).
Theoretical insights by Strasser and Weber (1999) open up the way to a unified treatment of a huge class of permutation tests.The coin package for conditional inference is the computational counterpart to this theoretical framework, implemented in the R system for statistical computing (R Core Team 2019).Hothorn, Hornik, van de Wiel, and Zeileis (2006) introduce the package and illustrate the transition from theory to practice.Here, we focus on the design principles upon which the coin implementation is based as well as on the more technical issues that need to be addressed in the implementation of such conceptual tools.Within package coin, formal S4 classes describe the data model and the conditional test procedures, consisting of multivariate linear statistics, univariate test statistics and a reference distribution.Thus, one can work with objects representing the theoretical entities nearly in the same way as one would work with mathematical symbols.Generic functions for obtaining statistics, conditional expectation and covariance matrices as well as p-value, distribution, density and quantile functions for the reference distribution help to extract information from these objects.The infrastructure is conveniently interfaced in the function independence_test(), thus providing the computational counterpart of the theoretical framework of Strasser and Weber (1999).
Here, we start out with an illustrative application of independence_test() to a small twosample problem (see Table 1) and then continue to introduce the underlying computational building blocks using the same data set.The data was used previously by Bergmann, Ludbrook, and Spooren (2000) as a challenging example in a comparison of test statistics and p-values of the Wilcoxon-Mann-Whitney rank sum test as reported by eleven statistical packages.More specifically, n = 24 rats received a fixed oral dose of a centrally acting muscle relaxant as active treatment or a saline solvent as control.The animals were placed on a rotating cylinder and the length of time each rat remained on the cylinder was measured, up to a maximum of 300 seconds.The rats were randomly assigned to the control and treatment groups, thus a re-randomization test as implemented in independence_test() is the appropriate way to investigate if the response is independent of the group assignment.The data are particularly challenging because of the many ties in the (right-censored) response (19 observations take the maximal value 300) and the quasi-complete separation (smaller values of time are only observed in the treatment group).Conceptually, this makes computation of the two-sided exact p-value for any two-sample contrast very simple: p = 2 • 19 12 / 24 12 = 0.03727.However, in software, this often makes computations of the exact p-value more difficult because several simple algorithms fail.
Utilizing coin, the hypothesis of independence of length of time and group assignment can be specified by a formula which, together with a data frame rotarod, serve as arguments to independence_test(): R> library("coin") R> data("rotarod", package = "coin") R> independence_test(time ~group, data = rotarod, + ytrafo = rank_trafo, distribution = exact()) Exact General Independence Test data: time by group (control, treatment) Z = 2.4389, p-value = 0.03727 alternative hypothesis: two.sided Here, the conditional Wilcoxon-Mann-Whitney test was performed via a rank transformation of the response, employing the exact distribution for obtaining the p-value (yielding the correct result outlined above).Users of R can easily interpret this output since it is represented in the same format as classical tests in the basic stats package.Based on the p-value derived from the exact conditional distribution of the test statistic Z, the independence of group assignment and time on the cylinder can be rejected.
Although the above piece of code looks embarrassingly simple, the underlying computations are much more sophisticated than visible at first sight: The data are pre-processed along with their transformations, deviations from independence are captured by a (possibly multivariate) linear statistic, standardized by conditional expectation and variance, and aggregated to a final test statistic.Subsequently, the requested reference distribution is computed, from which a p-value is derived, and everything is wrapped into an object that can be conveniently printed or queried for more detailed information.After briefly reviewing the underlying theory from Strasser and Weber (1999) in Section 2, we introduce formal S4 classes and methods capturing all the outlined steps in Section 3. Section 4 provides further information about high-level user interfaces and extensibility, and Section 5 further illustrates how to employ the software in practice using a categorical data example.Section 6 concludes the paper with a short discussion, some more details about the underlying theory can be found in an appendix.

Permutation tests in a nutshell
In the following we give a brief overview of the general theory for permutation tests as developed by Strasser and Weber (1999) and implemented by Hothorn et al. (2006).
The task is to test the independence of two variables Y and X from sample spaces Y and X which may be measured at arbitrary scales and may be multivariate as well.In addition, B ∈ {1, . . ., r}, a nominal factor with r levels, indicates a certain block structure of the observations: for example study centers in a multi-center randomized clinical trial where only a re-randomization of observations within blocks is admissible.We are interested in testing the null hypothesis of conditional independence of Y and X within blocks B against arbitrary alternatives, for example shift or scale alternatives, linear trends, association in contingency tables etc. Strasser and Weber (1999) where the linear statistic for each block is given by The function I(•) is the indicator function and vec denotes the vec operator (which stacks the columns of a matrix).Here, g : X → R p×1 is a transformation of the X measurements and  Strasser and Weber (1999) and are given in Appendix A. Having the conditional expectation and covariance at hand we are able to standardize an observed linear statistic t ∈ R pq (of the form given in Equation 1) and aggregate it to some univariate test statistic c = c(t, µ, Σ).Various choices for c(t, µ, Σ) are conceivable, e.g., a quadratic form or a maximum type statistic (see Section 3).In the latter case, a natural first step is to standardize each of the pq statistics in t by its expectation and standard deviation: where σ ∈ R pq is the conditional variance, i.e., the diagonal elements of the conditional covariance Σ.
In the following, we describe a class structure for representing these theoretical objects along with possible choices of test statistics c and computations or approximations of the associated reference distributions.

A class structure for permutation tests
In this section, the theory of permutation tests, as briefly outlined in the previous section, is captured in a set of S4 classes and methods: In Section 3.1, we suggest objects representing the data and the independence problem, for which a test statistic is computed subsequently.Objects for the associated reference distribution are constructed in Section 3.2.Finally, in Section 3.3, everything is combined in a single object for the whole testing procedure.

Data structure
We are provided with n realizations of (Y i , X i , B i , w i ).In addition to the variables X, Y, and B, it is convenient (for example to efficiently represent large contingency tables) to include case weights w i , defaulting to w i ≡ 1.This data structure is represented by class "IndependenceProblem": Class "IndependenceProblem"

Slot Class
x "data.frame"y "data.frame"block "factor" weights "numeric" Note that objects of this class implicitly define the null distribution H 0 and all admissible permutations of observations within blocks.
For our illustrating rotating rats example, we specify the hypothesis of independence of variables time and group by initializing a new independence problem with the corresponding observations:

Independence problems and linear statistics
The transformation functions g and h as well as the transformed observations g(X i ) and h(Y i ), i = 1, . . ., n, are added to the data structure by extending class "IndependenceProblem": Class "IndependenceTestProblem" Contains "IndependenceProblem"

Slot Class
xtrans "matrix" ytrans "matrix" xtrafo "function" ytrafo "function" The ytrafo and xtrafo slots correspond to the transformations h and g, respectively.The ith row of the n×q matrix ytrans corresponds to h(Y i ).Similarly, the rows of xtrans (n×p) correspond to g(X i ).Note that, in addition to the data, hypothesis and permutation scheme, the test statistic T is defined by objects of class "IndependenceTestProblem" as well.
In the case of both X and Y being univariate factors with p and q levels, g and h are the corresponding dummy codings and the linear statistic T is the (vectorized) p × q contingency table of X and Y.In the rats example, the default dummy coding for factor group is employed and a rank transformation (via rank_trafo()) is applied to time.

coin: Implementing a Class of Permutation Tests
R> itp <-new("IndependenceTestProblem", ip, ytrafo = rank_trafo) The linear statistic and its conditional expectation and covariance are stored in objects of class "IndependenceLinearStatistic": Class "IndependenceLinearStatistic" Contains "IndependenceTestProblem" Slot Class linearstatistic "matrix" expectation "matrix" covariance "matrix" The jth column of the pq × r matrices linearstatistic and expectation corresponds to T j and µ j , respectively.The covariance slot is a pq(pq + 1)/2 × r matrix where each column corresponds to the lower triangular elements of Σ j , but the full covariance matrix can be extracted as needed (see below).In the rotating rats example, such an object is easily created via

R> ils <-new("IndependenceLinearStatistic", itp)
Using methods for suitable generic functions (see Table 2), the linear statistic for the rotating rats can be extracted via

R> statistic(ils, type = "linear") control 180
This is simply the sum of the average ranks in the control group, i.e., 180 = 12 • 15 because all 12 observations have the maximal average rank 15.Additionally, the associated conditional mean and variance under H 0 can be computed using

R> variance(ils)
control 151.3043 based upon which we can now set up a test statistic.

Test statistics
The specification of the inference procedure is completed by the definition of a univariate test statistic c which is represented by a virtual class "IndependenceTestStatistic": Function Description statistic(object, type, partial) Extraction of the aggregated (partial = FALSE) or partial (partial = TRUE) linear statistic (type = "linear") t or t j , centered statistic (type = "centered") t − µ or t j − µ j , and standardized statistic (type = "standardized") z or z j , respectively.For objects inheriting from "IndependenceTestStatistic", the univariate test statistic c (type = "test").expectation(object, partial) Extraction of the aggregated (partial = FALSE) or partial (partial = TRUE) conditional expectation µ or µ j , respectively.variance(object, partial) Extraction of the aggregated (partial = FALSE) or partial (partial = TRUE) conditional variance σ or σ j , respectively.covariance(object, invert, partial) Extraction of the aggregated (partial = FALSE) or partial (partial = TRUE) conditional covariance (invert = FALSE) Σ or Σ j , and its Moore-Penrose inverse (invert = TRUE) Σ + or Σ + j , respectively.
Table 2: List of generic functions with methods for classes inheriting from "IndependenceLinearStatistic".
Class "IndependenceTestStatistic" Contains "IndependenceLinearStatistic" Slot Class teststatistic "numeric" standardizedlinearstatistic "numeric" The slot standardizedlinearstatistic contains z, the (possibly multivariate) linear statistic standardized by its conditional expectation and variance (Equation 2).Slot teststatistic is for storing univariate test statistics c.Methods for all generics listed in Table 2 are also available for objects of this class.coin implements three "IndependenceTestStatistic" subclasses with associated univariate test statistics: scalar test statistics c scalar , maximum-type statistics c max and quadratic forms c quad .In case of univariate linear statistics t, i.e., for pq = 1, a natural test statistic c is simply the standardized linear statistic (Equation 2) A special class is available for this Class "ScalarIndependenceTestStatistic" Contains "IndependenceTestStatistic"

Slot Class
alternative "character" paired "logical" that also defines a character vector specifying the alternative to test against ("two.sided","greater" and "less").Thus, the construction of a scalar test statistic corresponds to the construction of a suitable object via R> sits <-new("ScalarIndependenceTestStatistic", ils, + alternative = "two.sided")R> statistic(sits, type = "standardized") control 2.438909 which yields the standardized Wilcoxon-Mann-Whitney statistic reported in the introductory example.In the multivariate case (pq > 1), a natural extension is to employ a maximum-type statistic of the form where the definition reflects the associated alternative with name given in quotes.Again a special class "MaxTypeIndependenceTestStatistic" is available for this: Class "MaxTypeIndependenceTestStatistic" Contains "IndependenceTestStatistic"
It is computationally more expensive because the Moore-Penrose inverse Σ + of Σ is involved.Such statistics are represented by objects of class "QuadTypeIndependenceTestStatistic" defining slots for the lower triangular part of Σ + and its rank (degrees of freedom): Class "QuadTypeIndependenceTestStatistic" Contains "IndependenceTestStatistic" Slot Class covarianceplus "numeric" df "numeric" paired "logical" A slot alternative is not needed because, by construction, quadratic forms cannot be applied to one-sided hypotheses.

Representation of conditional null distributions
The conditional distribution (or an approximation thereof) and thus the p-value corresponding to the statistic c(t, µ, Σ) can be computed in several different ways.For some special forms of the linear statistic, the exact distribution of the test statistic is tractable.For twosample problems, the shift algorithm by Streitberg andRöhmel (1986, 1987) and the split-up algorithm by van de Wiel ( 2001) are implemented as part of the package.
Conditional Monte Carlo procedures can always be used to approximate the exact distribution.In this case, within each block, a sufficiently large number of random samples from all admissible permutations of the observations is drawn.The test statistic is computed for the permuted Y values and the distribution of these test statistics is an approximation to the conditional reference distribution.When p-values are computed, confidence intervals are available from the binomial distribution.
Strasser and Weber (1999, Theorem 2.3) showed that the conditional distribution of linear statistics T with conditional expectation µ and covariance Σ tends to a multivariate normal distribution with parameters µ and Σ as n i=1 I(B i = j)w i → ∞ for all j = 1, . . ., r.Thus, the asymptotic conditional distribution of the standardized linear statistic z is normal and therefore, p-values for scalar or maximum-type univariate statistics can be computed directly in the univariate case (pq = 1) or approximated by numerical algorithms (Genz 1992) as implemented in package mvtnorm (Genz, Bretz, and Hothorn 2008) in the multivariate setting.For quadratic forms c quad which follow a χ 2 distribution with degrees of freedom given by the rank of Σ (see Johnson and Kotz 1970, Chapter 29), asymptotic probabilities can be computed straightforwardly.
A null distribution is represented by either a distribution (and p-value) function only Class "PValue" Slot Class pvalue "function" p "function" name "character" or, where possible, augmented by, e.g., its density and quantile functions: Class "NullDistribution" Contains "PValue" Slot Class size "function" pvalueinterval "function" midpvalue "function" q "function" d "function" support "function" parameters "list" Currently, there are three classes extending "NullDistribution": "ExactNullDistribution", "ApproxNullDistribution" (adding slots seed and nresample containing, respectively, the current state of the random number generator and the number of Monte Carlo replicates) and "AsymptNullDistribution" (adding a slot seed).All of them can be queried for probabilities, quantiles, etc., using suitable methods (see Table 3 for an overview).New methods for computing or approximating the conditional distribution can be integrated into the framework via suitable inheritance from "PValue" (an example is given in Section 4).
The support() function returns the support of a discrete distribution.Computation of the p-value interval based on an observed test statistic c and its conditional null distribution.size(object, alpha, type) Computation of the test size at the nominal significance level α using either the rejection region corresponding to the p-value (type = "p-value") or the mid-p-value (type = "mid-p-value").dperm(object, x) Evaluation of the probability density function at x. pperm(object, q) Evaluation of the cumulative distribution function for quantile q. qperm(object, p) Evaluation of the quantile function for probability p. rperm(object, n) Generation of n random numbers from the null distribution.

support(object)
Extraction of the support of the null distribution.
Table 3: Classes and methods for conditional null distributions.
Using these tools, the exact p-value for the independence test on the rotating rats along with the complete exact reference distribution (allowing for the computation of quantiles, for example) can be derived via

Objects for conditional tests
A conditional test is represented by a test statistic of class "IndependenceTestStatistic" and its conditional null distribution inheriting from class "PValue".In addition, a character string giving the name of the test procedure is defined in class "IndependenceTest": Class "IndependenceTest" Slot Class distribution "PValue" statistic "IndependenceTestStatistic" estimates "list" method "character" call "call" Remember that objects of class "IndependenceTestStatistic" represent the data, hypothesis, linear statistic and test statistic along with conditional expectation and covariance matrix.
The estimates slot may contain parameter estimates where available, for example an estimate and corresponding confidence interval for a shift parameter derived from a conditional Wilcoxon-Mann-Whitney test.
A complete description of the conditional independence test for the rotating rats data is given by an object of class "IndependenceTest" which is conveniently represented by the corresponding show() method: R> new("IndependenceTest", statistic = sits, distribution = end) Exact General Independence Test data: time by group (control, treatment) c = 2.4389, p-value = 0.03727 Of course, the methods previously defined in this section (see Tables 2 and 3) are defined for objects of class "IndependenceTest" as well.Thus, all theoretical entities introduced in Section 2 are now captured in a single object of class "IndependenceTest" and all methods for extracting information from it are readily available.

Interfaces to permutation inference
In Section 3, all the necessary computational building blocks are introduced for implementing the general class of permutation tests outlined in Section 2. However, one rarely needs to exploit the full flexibility of each component of the framework.More often, one wants to employ sensible defaults for most (if not all) steps in the analysis but preserving the possibility to extend a few steps based on user-supplied methods.For this purpose, coin provides the function independence_test() as the main user interface for performing independence tests.Many of its arguments have flexible defaults or can be specified by a simple character string while still allowing to plug in much more complex user-defined objects, e.g., for data preparation, computation of the null distribution, or transformation functions.

A convenient user interface
Via the generic function independence_test(), all steps described in Section 3 can be carried out using a single command.The main workhorse behind it is the method for objects of class "IndependenceProblem": independence_test(object, teststat = c("maximum", "quadratic", "scalar"), distribution = c("asymptotic", "approximate", "exact"), alternative = c("two.sided","less", "greater"), xtrafo = trafo, ytrafo = trafo, scores = NULL, check = NULL, ...) Thus, object describes the data and the null hypothesis.Arguments xtrafo and ytrafo refer to the transformations g and h: Both are by default set to function trafo() which chooses suitable transformations based on the scale of the considered variables (see below).The three types of univariate test statistics discussed above are hard-coded and can be specified by a simple string.Similarly, the reference distribution and the alternative hypothesis can be supplied as strings.In addition to these simple specifications, more flexible specifications for the data (object), the transformations (xtrafo, ytrafo), and the distribution are available, as discussed in the following.The scores argument takes a named list of numeric vectors to be used as scores for ordered factors.Validity checks for objects of class "IndependenceProblem" specified to argument check can be used to test for certain aspects of the data, e.g., when one has to make sure that two independent samples are present.

Data specification
The standard way of specifying relationships between variables in R are formulas in combination with a data frame.Hence, a "formula" method for independence_test() is provided that interprets the left hand side variables of a formula as Y variables (univariate or possibly multivariate), the right hand side as X variables (univariate or multivariate as well).An optional blocking factor can specified after a vertical bar, e.g.,

y1 + y2 ~x1 + x2 | block
This specifies an independence problem between two Y variables and two X variables (in case all variables are numeric the linear statistic is four-dimensional with p = 2 and q = 2) for each level in block.As usual, data, weights and subset arguments can be specified as well.
Based on all these arguments the "IndependenceProblem" is built and simply passed on to the independence_test() method described above.
For (simple) categorical data, there is yet another way of specifying the independence problem, namely via the "table" method of independence_test().Its first argument is allowed to be a two-or three-dimensional table : The first two margins are interpreted as univariate categorical X and Y variables and the optional third margin is taken to be the blocking factor.Again, an "IndependenceProblem" object is derived from this and passed on to the associated method.

Null distributions
In the simplest case, the distribution argument of the independence_test() methods can be specified by a simple character string.However, this is not flexible enough in some situations and hence it is also possible to supply a function that can compute a "PValue" object from an "IndependenceTestStatistic" object.
For the most important special cases, suitable function generators are provided in coin.
For example, the function approximate(nresample = 1000) returns a Monte Carlo function that draws nresample (default: 10000) random permutations.Similarly, exact() and approximate() return functions computing the exact or asymptotic null distributions, respectively.Again, computational details in the computation of the null distribution can be controlled via arguments of the function generators.
Additionally, it is also possible to set distribution to a user-supplied algorithm for computing the conditional null distribution.It just has to be provided in the form of a function taking an object inheriting from "IndependenceTestStatistic" and returning an object inheriting from class "PValue".
As an example, consider the computation of the exact p-value for testing independence of two continuous random variables.The identity transformation is used for both g and h, thus a conditional version of the test for zero Pearson correlation is constructed.We use a tiny artificial example where enumeration and evaluation of all possible permutations is still feasible.The function sexact() extracts both variables, computes all permutations using the function permutations() from e1071 (Dimitriadou, Hornik, Leisch, Meyer, and Weingessel 2008), computes the linear test statistic for each permutation, standardizes it by expectation and variance, and then sets up the distribution function.
Due to this flexibility, almost all special-purpose functionality implemented in packages ex-actRankTests (Hothorn 2001;Hothorn andHornik 2002, 2006) and maxstat (Hothorn and Lausen 2002;Hothorn 2007) can be conveniently provided within the coin framework, so that both packages will become deprecated in the future.

Permutation tests in practice: A categorical data example
The job satisfaction of African-American males in the USA (Agresti 2002,  The positive diagonal and (mostly) negative off-diagonal elements convey that higher income categories seem to be associated with higher job satisfaction.Thus, to direct power against ordered alternatives, a linear-by-linear association statistic (Agresti 2002) should be used instead of the omnibus χ 2 statistic.This can be conveniently performed within coin, e.g., using simple equidistant scores and a Monte Carlo approximation of the null distribution: R> it <-independence_test(js, distribution = approximate(nresample = 10000), + scores = list(Job.Satisfaction = 1:4, Income = 1:4)) R> pvalue(it) [1] 0.0116 99 percent confidence interval: 0.009024885 0.014651019 Using this strategy, the null hypothesis of independence of job satisfaction and income can be rejected in favor of a positive association of both variables.Other choices of scores are also conceivable.Especially when there is an underlying numeric scale, interval midpoints are often used (see Hothorn et al. 2006, for an example).
For other patterns of dependence-e.g., when only a few cells in a large table deviate from independence-a maximum-type statistic is also useful for contingency tables.To complete our tour of coin tools for categorical data, we briefly illustrate this approach using the job satisfaction data again (even though a maximum-type statistic is clearly not very powerful for the dependence pattern in this data set).The maximum-type test is set up easily:

R> independence_test(js, teststat = "maximum")
Asymptotic General Independence Test data: Job.Satisfaction by Income (<5000, 5000-15000, 15000-25000, >25000) stratified by Gender maxT = 1.6258, p-value = 0.7215 alternative hypothesis: two.sided with its conditional asymptotic null distribution being available immediately (due to the joint multivariate normal distribution for the contingency table T).Single-step adjusted p-values for each cell of the contingency table corresponding to this maximum test can be computed via was used for the computations, Figure 1 was created using the vcd package (Meyer, Zeileis, and Hornik 2006).Streitberg B, Röhmel J (1987)

A. Expectation and covariance
The conditional expectation and covariance matrix of linear statistics T as given in Equation 1in Section 2 are computed as follows.Let w •j = n i=1 I(B i = j)w i denote the sum of the weights in block j and S j the set of all permutations of the observations in block j.The conditional expectation of the transformation h in block j is with corresponding conditional covariance matrix This is the basis for computing the conditional expectation and covariance of the linear statistic T j in block j: where ⊗ is the Kronecker product.The conditional expectation and covariance of T, aggregated over all r blocks, are then given by E(T|S j ) = µ = This function can then be passed to independence_test() for computing the distribution function and p-value: R> independence_test(y ~x, data = correxample, alternative = "less"independence_test() chooses the transformations g and h via the wrapper function trafo().This, in turn, chooses the actual transformation based on the scale of each variable (individually) in a data frame data: trafo(data, numeric_trafo = id_trafo, factor_trafo = f_trafo,  ordered_trafo = of_trafo, surv_trafo = logrank_trafo,  var_trafo = NULL, block = NULL) i = j)w i g(X i )  E(h|S j ) ⊤  and j |S j ), COV(T|S j ) = Σ = r j=1 Σ j = r j=1 COV(T j |S j ).
suggest deriving scalar test statistics for testing H 0 from multivariate linear is also called influence function and may depend on the full vector of responses (Y 1 , . .., Y n ), however only in a permutation symmetric way, i.e., the value of the function must not depend on the order in which Y 1 , . .., Y n appear.The case weights w i are assumed to be integervalued, indicating that w i observations of Y i , X i and B i , i = 1, . .., n, are available, with default w i ≡ 1.The distribution of T depends on the joint distribution of Y and X, which is unknown under almost all practical circumstances.At least under the null hypothesis one can dispose of this dependency by fixing X 1 , . . ., X n and conditioning on all possible permutations of the responses Y 1 , . . ., Y n within block j, where j = 1, . . ., r.The conditional expectation µ ∈ R pq and covariance Σ ∈ R pq×pq of T under H 0 given all permutations S of the responses are derived by Similarly to the Mood test, the conditional counterpart of many other classical tests (some of them implemented in package stats) are easily available through independence_test() by specifying the appropriate xtrafo, ytrafo and teststat.This includes the Wilcoxon-Mann-Whitney or Cochran-Mantel-Haenszel tests, but also many other well-known tests as shown in Table

Table 7
Thus, the test does not indicate significant departure from independence.However, ordering of the factors is not exploited by the Cochran-Mantel-Haenszel test.As some positive correlation of the two factors would seem natural, it is worth having a closer look at the data and the test result.The underlying linear statistic is This is exactly the original contingency table aggregated over the block factor Gender: .8) is described by measures of income and reported job satisfaction in four (ordinal) classifications.It seems natural to surmise that job satisfaction increases with income.The data (see Figure1) is given as a three-dimensional "table" with variables Income and Job.Satisfaction according to Gender (labels slightly modified for convenience): . "Exakte Verteilungen für Rang-und Randomisierungstests im allgemeinen c-Stichprobenfall."EDV in Medizin und Biologie, 18(1), 12-19.van de Wiel MA (2001)."The Split-Up Algorithm: A Fast Symbolic Method for Computing p-Values of Distribution-Free Statistics."Computational Statistics, 16, 519-538.Westfall PH, Young SS (1993).Resampling-Based Multiple Testing.John Wiley & Sons, New York.