Robust Maximum Association Between Data Sets: The R Package ccaPP

This package vignette is an up-to-date version of Alfons, Croux, and Filzmoser (2016), published in the Austrian Journal of Statistics. An intuitive measure of association between two multivariate data sets can be defined as the maximal value that a bivariate association measure between any one-dimensional projections of each data set can attain. Rank correlation measures thereby have the advantage that they combine good robustness properties with good efficiency. The software package ccaPP provides fast implementations of such maximum association measures for the statistical computing environment R. We demonstrate how to use ccaPP to compute the maximum association measures, as well as how to assess their significance via permutation tests.


Introduction
Projection pursuit allows to introduce intuitive and therefore appealing association measures between two multivariate data sets.Suppose that the data sets X and Y consist of p and q variables, respectively.A measure of multivariate association between X and Y can be defined by looking for linear combinations Xα and Y β having maximal association.Expressed in mathematical terms, we define an estimator ρR (X, Y ) = max where R is an estimator of a bivariate association measure R such as the Pearson correlation, or the Spearman or Kendall rank correlation.Using the projection pursuit terminology, R is the projection index to maximize.The projection directions corresponding to the maximum association are called weighting vectors and are estimated by R(Xα, Y β).
(2) Alfons, Croux, and Filzmoser (2016) developed the alternate grid algorithm for the computation of such maximum association estimators and studied their theoretical properties for various association measures.It turns out that the Spearman and Kendall rank correlation yield maximum association estimators with good robustness properties and good efficiency.
This paper is a companion paper to Alfons et al. (2016) that demonstrates how to apply the maximum association estimators in the statistical environment R (R Core Team 2015) using the add-on package ccaPP (Alfons 2015).The package is freely available on CRAN (Comprehensive R Archive Network, http://CRAN.R-project.org).
Note that using the Pearson correlation as the projection index of the maximum association estimator corresponds to the first step of canonical correlation analysis (CCA; see, e.g., Johnson and Wichern 2002), hence the package name ccaPP.Since CCA is a widely applied statistical technique, various algorithms and extensions are implemented in R packages on CRAN.Two important examples are briefly discussed in the following.The package CCA (González, Déjean, Martin, and Baccini 2008;González and Déjean 2012) extends the built-in R function cancor() with additional numerical and graphical output.Moreover, it provides a regularized version of CCA for data sets containing a large number of variables.Bayesian models and inference methods for CCA are implemented in the package CCAGFA (Klami, Virtanen, and Kaski 2013;Virtanen, Leppaaho, and Klami 2015).
The remainder of the paper is organized as follows.In Section 2, the design and implementation of the package are briefly discussed.Section 3 demonstrates how to compute the maximum association estimators, and Section 4 illustrates how to test for their significance.
A comparison of computation times is given in Section 5.The final Section 6 concludes the paper.

Design and implementation
Various bivariate association measures and the alternate grid algorithm for the maximum association estimators are implemented in C++, and integrated into R via the package RcppArmadillo (Eddelbuettel and Sanderson 2014;Eddelbuettel, François, and Bates 2015).
The following bivariate association measures are available in the package ccaPP: It should be noted that these are barebones implementations without proper handling of missing values.Hence the first three functions come with a substantial speed gain compared to R's built-in function cor().Moreover, the fast O(n log(n)) algorithm for the Kendall correlation (Knight 1966) is implemented in corKendall(), whereas cor() uses the naive O(n 2 ) algorithm.
The alternate grid algorithm for the maximum association estimators is implemented in the function maxCorGrid().Any of the bivariate association measures above can be used as projection index, with the Spearman rank correlation being the default.We do not recommend to use the quadrant correlation since its influence function is not smooth, which may result in unstable estimates of the weighting vectors.For more details on the theoretical properties of the maximum association estimators, the reader is referred to Alfons et al. (2016).
To assess the significance of a maximum association estimate, a permutation test is provided via the function permTest().Parallel computing to increase computational performance is implemented via the package parallel, which is part of R since version 2.14.0.

Maximum association measures
In this section, we show how to apply the function maxCorGrid() from the package ccaPP to compute the maximum association estimators.We thereby use the classic diabetes data (Andrews and Herzberg 1985, page 215), which are included as example data in the package.
First we load the package and the data.All measurements are taken for a group of n = 76 persons.
library("ccaPP") data("diabetes") x <-diabetes$x y <-diabetes$y Component x consists of p = 2 variables measuring relative weight and fasting plasma glucose, while component y consists of q = 3 variables measuring glucose intolerance, insulin response to oral glucose and insulin resistance.It is of medical interest to establish a relation between the two data sets.
The function maxCorGrid() by default uses the Spearman rank correlation as projection index.
The estimated weighting vectors can be accessed through components a and b of the returned object, respectively.
With the argument method, another bivariate association measure can be set as projection index, e.g., the Kendall rank correlation, the M-association or the Pearson correlation.
Note that the Spearman and Kendall rank correlation estimate different population quantities than the Pearson correlation.Thus the above values of the different maximum association measures are not directly comparable.The argument consistent can be used for the former two methods to get consistent estimates of the maximum correlation under normal distributions.
The M-association measure is consistent at the normal model and estimates the same population quantity as the Pearson correlation.

Permutation tests
To assess the significance of maximum association estimates, permutation tests can be performed with the function permTest().The number of random permutations to be used can be set with the argument R, which defaults to 1000.On machines with multiple processor cores, only the argument nCores needs to be set to take advantage of parallel computing in order to reduce computation time.If nCores is set to NA, all available processor cores are used.
In the examples in this section, we use 2 processor cores.Furthermore, we set the seed of the random number generator via the argument seed for reproducibility of the results.Since we employ parallel computing, ccaPP uses random number streams (L'Ecuyer, Simard, Chen, and Kelton 2002) from the package parallel rather than the default R random number generator.
permTest(x, y, nCores = 2, seed = 2014) ## ## Permutation test for no association ## ## r = 0.534699, p-value = 0.001000 ## R = 1000 random permuations ## Alternative hypothesis: true maximum correlation is not equal to 0 Again, the Spearman rank correlation is used as projection index by default.A different bivariate association measure can be specified via the argument method, which is passed down to the function maxCorGrid().
Since the focus of ccaPP is on robustness, we introduce an outlier into the diabetes data as in Taskinen, Kankainen, and Oja (2003).More precisely, we replace the value 0.81 of the first observation of variable glucose intolerance by 8.1, i.e., by a simple shift of the comma.

Computation times
This section analyzes the computation times of the methods implemented in ccaPP.All computations are performed in R version 3.2.2 on a machine with an Intel Xeon X5670 CPU.The computation times are recorded with the R package microbenchmark (Mersmann 2014).
First, we compare the barebones implementations of the Pearson, Spearman and Kendall correlations (functions corPearson(), corSpearman() and corKendall() in ccaPP) with their counterparts from the base R function cor().We also include the M-association measure from the function corM() in the comparison.The bivariate association measures are computed for 10 random draws from a bivariate normal distribution with true correlation ρ = 0.5 and sample size n = 100, 1 000, 10 000, 100 000.For each random sample, computation times from 10 independent runs are recorded.
Table 1 contains the average computation times of the bivariate association measures.Clearly, the fast O(n log(n)) algorithm for the Kendall correlation (Knight 1966) in ccaPP is a huge improvement over the naive O(n 2 ) algorithm in base R. Time savings for the Spearman and Pearson correlation are also substantial, considering that they are only due to a lack of missing data handling.For the M-association, the computation time is somewhat higher than that of the Spearman and Kendall correlation.
Since the projection pursuit algorithm for the maximum association measures involves computing a large number of bivariate associations (see Alfons et al. 2016), the faster barebones implementations are crucial to keep the computation of the maximum association feasible.We employ the same procedure as above to record the computation time of the maximum association measures, except that each of the random samples is drawn from a multivariate normal distribution such that the true maximum correlation is ρ = 0.5 and the corresponding weighting vectors are α = (1, 0, . . ., 0) and β = (1, 0, . . ., 0) .The sample size is set to n = 100, 1 000, 10 000, the dimension of X is p = 5, 10, 50, and the dimension of Y is q = 1, 5, 10, 50.
Inspired by canonical correlation analysis (CCA), we also compute other association measures for comparison.In CCA, the first canonical correlation is given by the square root of the largest eigenvalue of the matrix where Σ XX = Cov(X), Σ Y Y = Cov(Y ), Σ XY = Cov(X, Y ) and Σ Y X = Σ XY (see, e.g., Johnson and Wichern 2002).This is of course identical to the maximum association measure with the Pearson correlation as projection index.Other association measures are obtained by plugging different scatter matrices into (3).However, such a measure is in general different from the maximum association measure based on the corresponding bivariate association, with the maximum association being much easier to interpret.Here we plug in scatter matrices ccaPP: Robust Maximum Association corresponding to the Pearson, Spearman and Kendall correlation.For the Pearson correlation, the corresponding scatter matrix is the sample covariance matrix.For the Spearman and Kendall correlation, the scatter matrices are given by the respective pairwise associations multiplied with scale estimates of the corresponding variables.Furthermore, since a multivariate M-estimator of the covariance matrix is not robust, we instead use the minimum covariance determinant estimator (MCD; see Rousseeuw and Van Driessen 1999).
Table 2 lists average computation times for various values of n, p and q.The function maxCorGrid() is thereby used with the default values for all control parameters of the algorithm (see the corresponding R help file).For the maximum association measures, the number of bivariate associations that have to be computed clearly takes a toll on computation time compared to the association measures based on the full scatter matrices.Note that the Kendall correlation is the exception, as the computation of the full scatter matrix uses R's built-in cor() function, and therefore the naive O(n 2 ) algorithm.Also note that computing the full MCD scatter matrix requires more observations than variables, i.e., n > p + q, hence it cannot be computed for n = 100 and p = q = 50.
For the Pearson correlation, the projection pursuit algorithm to find the maximum association cannot be recommended since the first canonical correlation is much faster to compute.However, the focus of ccaPP is on the Spearman and Kendall rank correlation, for which the maximum association measures are much more intuitive than the association measures based on the full scatter matrix.In our opinion, the gain of easy interpretability outweighs the increased computational cost.In any case, the maximum association measures are still reasonably fast to compute for many problem sizes due to our C++ implementation.
It is also worth noting that the association measures based on a full scatter matrix require the number of observations to be larger than the number of variables in each of the two data sets, i.e., n > max(p, q).The maximum association measures do not have this limitation, although computation time increases considerably in high dimensions.

Conclusions
The package ccaPP provides functionality for the statistical computing environment R to compute intuitive measures of association between two data sets.These maximum association measures seek the maximal value of a bivariate association measure between one-dimensional projections of each data set.We recommend the maximum Spearman and Kendall rank correlation measures because of their good robustness properties and efficiency.For details on the theoretical properties of the estimators, as well as the alternate grid algorithm and extensive numerical results, the reader is referred to Alfons et al. (2016).
Due to our C++ implementation, the maximum association measures are reasonably fast to compute.The significance of maximum association estimates can be assessed via permutation tests, which allow for parallel computing to decrease computation time.In addition, the corresponding functions in ccaPP are easy to use.

Table 1 :
Average computation time (in milliseconds) of the bivariate association measures in base R and the package ccaPP.

Table 2 :
Average computation time (in seconds) of the maximum association measures in package ccaPP, as well as association measures based on corresponding full correlation matrix.